Ground state degeneracy of a variation of Toric Code model We know that the ground state degeneracy of Toric Code model is 4. An easy way of seeing this is the following: 
Consider a 2D spin model where all the spins live on the links. The Hamiltonian is $H=-\sum_v{A_v}-\sum_p{B_p}$, where $A_v=\sigma_{v1}^z\sigma_{v2}^z\sigma_{v3}^z\sigma_{v4}^z$ is the product of 4 spins around a vertex $v$, and $B_p=\sigma_{p1}^x\sigma_{p2}^x\sigma_{p3}^x\sigma_{p4}^x$ is the product of 4 spins around a plaquette. All the $A_v$ and $B_p$'s commute with each other and they all square to one; thus each $A_v$ and $B_p$'s can take values $\pm1$. 
To see the ground state degeneracy, we notice that $\prod_vA_v=1$ and $\prod_pB_p=1$, being two constraints on the $A_v$ and $B_p$'s. Thus different $A_v$ and $B_p$'s can now only label $2^{2N}/4$ different configurations ($N$ being the total vertex number) because of the 2 constraints. But the total Hilbert space is still $2^{2N}$-dimensional, thus on average each configuration has a 4-fold degeneracy.
Now suppose we add a third term in Toric Code model, making it $H=-\sum_v{A_v}-\sum_p{B_p}-\sum_p{C_p}$, where $C_p=\sigma_{p1}^y\sigma_{p2}^y\sigma_{p3}^y\sigma_{p4}^y$ being the product of 4 $\sigma_y$'s  around a plaquette. What happens to the ground state degeneracy?
Note that after adding the third term, the Hamiltonian is still a sum of commuting projectors, thus the counting method should still work; I'm just having trouble counting the number of constraints. 
My guess is that after adding the $C$ term, the ground state degeneracy reduces to 2. Anyone can prove/disprove me?
 A: That is an interesting problem. I think your conjecture is partially correct. After adding the $H_C=-\sum_pC_p$ term to the Hamiltonian, the ground state degeneracy will be reduced for sure. The remaining degeneracy can be 2-fold in some cases. But for the most general case, the degeneracy can be completely lifted, leaving unique ground state.
I think it would be easier to prove it using the loop algebra method developed in this work [Y.-Z. You and X.-G. Wen, arXiv:1204.0113] (which you may have already known). Admittedly, the constraint-counting method you mentioned is more rigorous, and is successful for the original toric code model. However for your case with $H_C$ term, after thinking for a while, I could not figure out how to count the constraints either.
Here is how the loop algebra works.

First, let us define 4 loop operators according to the above figure (red - $\sigma^x$, green - $\sigma^y$, blue - $\sigma^z$)
$$Q_1=\prod_{l\in\text{x-line}}\sigma_l^x, Q_2=\prod_{l\in\text{y-line}}\sigma_l^x, 
P_1=\prod_{l\in\text{x-line}}\sigma_l^z, P_2=\prod_{l\in\text{y-line}}\sigma_l^z.$$
Assuming periodic boundary condition, the loops wind around the torus. Operators $Q_1$ and $Q_2$ measure the $\mathbb{Z}_2$ flux through the two torus holes respectively. While $P_1$ and $P_2$ are responsible to generate or to remove such fluxes (by moving a pair of $m$-excitations around the torus). In the original toric code model, the 4 degenerated ground states on the torus can be labeled by the eigen values of $Q_1$ and $Q_2$ (i.e. the fluxes through the holes), and different ground states are connected by the action of $P_1$ and $P_2$. Therefore the ground state Hilbert space is just the representation space of these loop operators.
It can be easily verified that these loop operators follows the algebra: $Q_1P_2=-P_2Q_1$, $Q_2P_1=-P_1Q_2$, $[Q_1,Q_2]=[P_1,P_2]=[Q_1,P_1]=[Q_2,P_2]=0$ (by checking the bound overlaps). The 4 loop operators can be divided into two anti-commuting pairs. To construct the representation for these operators, we start from $\{Q_1,P_2\}=0$. The anti-commutation relation requires $Q_1$ and $P_2$ to be represented by two Pauli matrices, say $Q_1\bumpeq\sigma_3$ and $P_2\bumpeq\sigma_1$. The same applies for $\{Q_2,P_1\}=0$. However the representation space must be enlarged such that the commutation relations between the anti-commuting pairs are also satisfied:
$$Q_1\bumpeq\sigma_3\otimes\sigma_0, Q_2\bumpeq\sigma_0\otimes\sigma_3, P_1\bumpeq\sigma_0\otimes\sigma_1, P_1\bumpeq\sigma_1\otimes\sigma_0.$$
For the original toric code model, it is easy to show that the Hamiltonian commutes with all the 4 loop operators, therefore it must act trivially in the representation space of loop operators, $H\bumpeq\sigma_0\otimes\sigma_0$. That is the algebraic reason that the ground states (and in fact each energy level) must be (at least) 4-fold degenerated.
Now let us consider the effect of adding $H_C$ term. All the $C_p$ operators can be divided into two classes: those away from $Q_{1,2}$ loops as $C_{p1}$ in the figure, and those adjacent to $Q_{1,2}$ loops as $C_{p2}$ in the figure. The Hamiltonian is the sum of these two classes of $C_p$ operators: $H_C=H_{1}+H_{2}$, with $H_{1,2}=-\sum C_{p1,p2}$. Note that $H_1$ commutes with all the loop operators still. However $H_2$ anti-commutes with $Q_1$ or $Q_2$ while commutes with $P_1$ and $P_2$. Therefore $H_1$ and $H_2$ are represented differently: $H_1\bumpeq\sigma_0\otimes\sigma_0$, $H_2\bumpeq a\sigma_1\otimes\sigma_0+ b\sigma_0\otimes\sigma_1+ c\sigma_1\otimes\sigma_1$ (combined with some coefficients $a,b,c$ to be determined case by case). $H_2$ induces mixing between the degenerated ground states. According to the idea of degenerate perturbation, the degeneracy will be lifted, and in the worst case, completely lifted.  
