  We propose a novel prescription for calculating the entanglement entropy of
the $SU(N)$ Yang-Mills gauge theories on the lattice under the strong coupling
expansion in powers of $\beta=2N/g^{2}$, where $g$ is the coupling constant.
Using the replica method, our Lagrangian formalism maintains gauge invariance
on the lattice. At $O(\beta^{2})$ and $O(\beta^{3})$, the entanglement entropy
is solely contributed by the central plaquettes enclosing the conical
singularity of the $n$-sheeted Riemann surface. The area law emerges naturally
to the highest order $O(\beta^{3})$ of our calculation. The leading $O(\beta)$
term is negative, which could in principle be canceled by taking into account
the "cosmological constant" living in interface of the two entangled
subregions. This unknown cosmological constant resembles the ambiguity of edge
modes in the Hamiltonian formalism. We further speculate this unknown
cosmological constant can show up in the entanglement entropy of scalar and
spinor field theories as well. Furthermore, it could play the role of a
counterterm to absorb the ultraviolet divergence of entanglement entropy and
make entanglement entropy a finite physical quantity.
