What exactly happens at the second-order phase transition of the 2D Toric code? For a 2D Toric code specified by
$$H = -J_s\sum_{s} \prod_{j\in s} \sigma^x_j - J_p\sum_{p} \prod_{j\in p} \sigma^z_p - h_x\sum_{l} \sigma^x_l - h_z\sum_{l} \sigma^z_l$$
where $s$ denotes stars, $p$ denotes plackets, and $l$ denotes links, there are known to be critical values of $h_x$ and $h_z$ that mark second-order phase transitions. In particular, the regime with the weaker external field is characterized by "confined excitations" while the regime with the stronger external field is characterized by "condensation" of those excitations. What do the wave functions and energies of states in these phases look like mathematically, and what causes this phase transition to occur?
When $h_x=0$, and $h_z=\infty$, I can see that the ground state is given by the "all-up" state in the $\sigma^z$ basis which is an equal superposition of all states in the $\sigma^x$ basis:
$$|\uparrow_1\dots \uparrow_N\rangle = (|\leftarrow\rangle + |\rightarrow\rangle)^{\otimes N}$$
so this state looks like a condensate of $x$-type bosons (electric charges). However, I am not sure what happens in the other regimes of each phase when $h_z$ is finite.
 A: The ground state of the toric code can be understand as a superposition of all loop configurations in the $z$ basis. The fact that these loops fluctuate at all length scales (and thus around the torus) leads to the topological order in the system. 
The $\sigma_z$ terms lead to a "tension" in the loops, penalizing long loops. Eventually, this tension will make very long loops impossible, giving rise to a "typical loop length" which is independent of the system size (similar to a correlation length) and lead to a phase transition to a trivial phase.
The $\sigma_x$ term does similar things in a dual basis.  (In some sense, it leads to a "breaking of loops", though one has to be careful with this picture.)
Note that this is only a qualitative, not an exact, picture of the ground state. In fact, already the Toric Code with only a $z$ field maps to the 3D classical Ising model for which no exact solution is known.
This phase transition has been studied by various means, see e.g. http://arxiv.org/abs/1012.1740 and the references therein.
