Why is the Chern number an integer? I have relatively limited knowledge on the Chern number, and I know that there exists high-level math proofs that the Chern number is an integer, but let me try to focus on the case I have in mind. $\DeclareMathOperator{\Tr}{Tr}$
Let $\Omega(\phi)$ denote the normalized states in a Hilbert space where $\phi\mapsto \Omega(\phi)$ is smooth on some closed orientable 2D manifold $M$. Let $P(\phi)$ denote the corresponding projection operator, i.e., $P(\phi) = |\Omega\rangle \langle \Omega|$. Let me define the connection $A = -i\langle \Omega|d\Omega\rangle$. Then by Stoke's theorem, it seems that
$$
\oint_\gamma A=\int_\text{interior}dA=-\int_\text{exterior}dA
$$
where the interior and exterior of $\gamma$ in $M$ are determined by the orientation of $M$ and $\gamma$. Notice that $idA = \Tr (PdP\wedge dP)$ is gauge-invariant.
Now, the usual argument is that $\oint_\gamma A$ is the Berry phase and only defined module $2\pi$ and thus $\int_M dA \in 2\pi \mathbb{Z}$. However, I don't quite understand why. Indeed, $A$ is a well-defined differential form, and thus is well-defined without any knowledge of the Berry phase. Indeed, we seem to have
$$
\int_M dA= \int_\text{interior} dA +\int_\text{exterior} dA=0
$$
And not $2\pi \mathbb{Z}$. Of course, I know this can't be quite true since there are simple examples such as the Bloch sphere. However, I'm having trouble seeing where I went wrong.
Follow-up. @Heidar' pointed out that $A$ is only locally defined. This can also be seen in the example of the Bloch sphere (my incomplete answer).
Now it seems that the only missing puzzle is why is the Chern number quantized. Indeed, in the Bloch sphere example, I can move the "singularity" of $A$ freely via a gauge transform. In this case, we can use properties of the Berry phase to rigorously show that the Chern number is quantized. However, in general, how would one "move such singularities"?
 A: Here are some comments that might be useful. It does not address the quantization of the integrals, but a related potential misunderstanding that might clear things up. It was too long to post as comments.
You say that $A$ is a well-defined differential form, I think this is wrong. In general if $$dF=0,$$ that does not imply that $$F=dA.$$ In general you also have $$F=dA + \omega$$ where $\omega$ is a so-called Harmonic form (look up Hodge decomposition in Nakahara for example) and captures the non-trivial topology of the underlying manifold (de Rahm cohomology).
In this formulation, both $A$ and $\omega$ are well-defined differential forms and you have (with appropriate normalization) $$\oint F = \oint (dA+\omega) = \oint \omega \in 2\pi\mathbb Z.$$
Your argument that $\oint dA = 0$ when $A$ is a globally well-defined differential form is correct. However, the harmonic form plays a crucial role.
But usually we don't of think of this in terms of harmonic forms, instead we consider $F = dA$ locally on each local patch and glue things together everywhere on the Manifold (transition functions on the intersection of local patches, look up principal bundles). This means that in this formulation, the field strength $F$ is a globally defined differential form, but the gauge field $A$ is only a differential form locally. It is NOT a well-defined differential form globally. You can consider it as a well-defined object globally, but this object is not a differential form, it's instead a Deligne-Beilinson cocycle (or, equivalently, a Cheeger-Simons differential). See https://arxiv.org/abs/1401.0740 for some more information. An ordinary differential form is the special case where the differential form is defined globally, and all transition forms are trivial.
Let's take an elementary example. Consider a vector field $\mathbf V(x,y)$ on a 2-torus, ie $$\mathbf V(x+2\pi,y) = \mathbf V(x,y+2\pi) = \mathbf V(x,y).$$
Assume that (similar to $dF=0$)
$$ \nabla\times \mathbf V = \mathbf 0.$$
A general vector field on the torus, satisfying the above, can be written as
$$ \mathbf V(x,y) = \nabla\phi(x,y) + \frac 1{2\pi}\mathbf\alpha,$$
where $\phi(x,y)$ is a well-defined $\mathbb R$-valued function on the torus, ie
$$ \phi(x+2\pi,y) = \phi(x,y+2\pi) = \phi(x,y), $$
and
$$ \mathbf \alpha = (\alpha_x, \alpha_y)\in\mathbb R^2,$$
is a constant vector. Here $\mathbf \alpha$ is the "Harmonic form". Integrating along, say, the x-cycle of the torus we get
$$ \oint_x \mathbf V\cdot d\mathbf x = \frac 1{2\pi}\int_0^{2\pi}\alpha_x = \alpha_x.$$
This is non-zero, despite $\nabla\times\mathbf V = \mathbf 0$, since the first De Rahm cohomology class of $T^2$ is non-trivial. You could ask why not include $\mathbf\alpha$ inside the definition of $\phi$ and write
$$ \mathbf V(x,y) = \nabla\tilde\phi(x,y) = \nabla\left[\phi(x,y) + x\alpha_x + y\alpha_y\right]. $$
This would be perfectly fine on $\mathbb R^2$, but not on $T^2$ since the function is $\tilde\phi(x,y)$ is not periodic and thus not globally defined on the full manifold. Essentially, it's not single-valued. But on local patches of the torus, this is well-defined. So you could define the vector field as $\mathbf V(x,y) = \nabla\tilde\phi(x,y)$ locally, and have a function $\tilde\phi(x,y)$ on each local patch of the torus such that they glue consistently together to define $\mathbf V$ globally. In this case you should be able to recover vector fields that satisfy $\nabla\times\mathbf V=\mathbf 0$ but $\oint \mathbf V\cdot d\mathbf r \not=0$ (when integrated long non-contractable homology cycles).
A: I think I understand why my reasoning is incorrect, though I would be happy if someone confirmed this (and possibly answer the follow-up question)
Indeed, I asked a related question here in math stackexchange on the example of the Bloch sphere, and it seems that in general $\Omega$ cannot be smooth on the entire closed 2-manifold, and thus I cannot apply Stoke's theorem on both the interior and exterior.
Indeed, consider the example of the Bloch sphere where
$$
\Omega(\theta,\phi) = \left(\cos{\frac{\theta}{2}}, \sin {\frac{\theta}{2}} e^{i\phi} \right)
$$
so that the local connection is given
$$
A = \frac{1}{2}(1-\cos\theta)d\phi
$$
It's clear that $\Omega$ and correspondingly $A$ is not quite smooth at the south pole $\theta=\pi$ due to the indeterminant $\phi$ value. Therefore, if $\gamma$ is a simple closed contour in the 2-sphere such as $\gamma(\phi) = (\theta_0,\phi)$ where $0<\theta_0<\pi$, we see that Stoke's theorem can only be applied to the upper hemisphere, i.e. $0<\theta<\theta_0$, i.e.,
$$
\oint_\gamma A=\int_\text{upper} dA
$$
Notice that by a gauge transform, I can move the singularity at the south pole to the north pole, e.g., consider
$$
\Omega'(\theta,\phi)=\left(\cos{\frac{\theta}{2}}  e^{-i\phi}, \sin {\frac{\theta}{2}} \right)
$$
Then the corresponding $A'$ is smooth on $\mathbb{S}^2$ except at the north pole, and thus we can apply Stoke's theorem on the lower hemisphere, i.e., $\theta_0 <\theta<\pi$, so that
$$
\oint_\gamma A'=-\int_\text{lower}dA'
$$
Notice that $dA=dA'$ since it is gauge-invariant. Since the Berry phase is also gauge-invariant, we have
$$
\exp \left(-i\oint_\gamma A'\right) =\exp\left(-i \oint_\gamma A\right)
$$
Hence,
\begin{align}
\oint_\gamma A -\oint_\gamma A' &\in 2\pi \mathbb{Z} \\
\int_\text{upper} dA-\int_\text{lower} dA &\in 2\pi \mathbb{Z}\\
\int_{\mathbb{S}^2} dA &\in 2\pi \mathbb{Z}
\end{align}
A: The reason why the Chern number is not always zero has been addressed in comments and other answers - namely, the Berry curvature $F$ is generally not exact over all of $M$.  If it is globally exact, then one has that
$$\mathrm{Ch} := \frac{1}{2\pi}\int_M F = \frac{1}{2\pi} \int_M \mathrm dA \overset{\text{Stokes}}{=} \frac{1}{2\pi}\int_{\partial M} A = 0$$
where we've used that $\partial M = \emptyset$.  In the more general case, we need to utilize a rather beautiful generalization to the Stokes theorem which utilizes the holonomy of the Berry connection.

Consider a space $N$ with boundary $\partial N$ as shown below,which is covered by two charts $U_1$ and $U_2$.  The cross-hatched region denotes the overlap $U_1\cap U_2$.

The holonomy of the Berry connection around a closed loop $\gamma$, denoted $\mathrm{hol}(\gamma)$, is the $U(1)$-valued object which yields the phase factor obtained by parallel transporting around $\gamma$.  If the connection 1-form is $\mathcal A = iA$, then
$$\mathrm{hol}(\gamma) := \exp\left[\int_\gamma \mathcal A\right]$$
In physics we identify this as the Berry phase obtained via transport around $\gamma$.
Let's say that the Berry curvature $\mathcal F$ is not exact on $N$ - however, it can be written $\mathcal F = \mathrm d\mathcal A_1$ on $U_1$, and $\mathcal F = \mathrm d\mathcal A_2$ on $U_2$.  Because $F$ is not globally exact, we cannot say that
$$\int_N \mathcal F = \int_{\partial N} \mathcal A$$
However, consider the introduction of a new loop $\gamma$ shown below, which starts and ends at a point $p$.  Note that $\gamma$ separates $N$ into an interior region $\mathrm{int}$ and an exterior region $\mathrm{ext}$.

We then have that
$$\mathrm{hol}(\partial N)\cdot \mathrm{hol}(\gamma) \cdot \mathrm{hol}(\gamma)^{-1} =\mathrm{hol}(\partial N)\tag{$\star$}$$ where $\mathrm{hol}(\gamma)^{-1}$ is the holonomy obtained by traversing $\gamma$ in the opposite direction.  Note that
$$\mathrm{hol}(\partial N)\cdot \mathrm{hol}(\gamma) = \exp\left[\int_{\partial N}\mathcal A_1\right] \cdot \exp\left[\int_{\gamma} \mathcal A_1\right] = \exp\left[\int_{\partial \mathrm{ext}} A_1\right]\overset{\text{Stokes}}{=} \exp\left[\int_\mathrm{ext} \mathcal F\right]$$
$$\mathrm{hol}(\gamma)^{-1} = \exp\left[\int_{\partial\mathrm{int}} \mathcal A \right] \overset{\text{Stokes}}{=} \exp\left[\int_{\mathrm{int}} \mathcal F\right]$$
and so therefore
$$\mathrm{hol}(\partial N)\cdot \mathrm{hol}(\gamma) \cdot \mathrm{hol}(\gamma)^{-1}= \exp\left[\int_\mathrm{ext} \mathcal F\right]\cdot \exp\left[\int_{\mathrm{int}} \mathcal F\right] = \exp\left[\int_N \mathcal F\right]$$
But from the above, this is simply equal to $\mathrm{hol}(\gamma) = \exp\left[\int_{\partial N} \mathcal A\right]$, and so we see that while Stokes' theorem itself does not hold, its "exponentiated" version does:
$$\exp\left[\int_{N} \mathcal F\right] = \exp\left[\int_{\partial N} \mathcal A\right]$$
This argument can be extended to any number of patches you'd like.  For each additional patch $U'$, you'd simply add a curve $\gamma'$ which encircles it through the overlap region and then add $\mathrm{hol}(\gamma') \cdot \mathrm{hol}(\gamma')^{-1}$ to $(\star)$.

With this result in hand, we can now address the question.  Let $M$ be a closed manifold and choose $\gamma$ to be a contractible loop in $M$. Assume that in a neighborhood of $\gamma$, we can write $\mathcal F= \mathrm d\mathcal A$.  Note that $\gamma$ divides $\mathcal M$ into two regions, which we arbitrarily call the interior and exterior.  On the interior (whose boundary is $\gamma$), we have
$$\exp\left[\int_\gamma \mathcal A\right] = \exp\left[\int_{\mathrm{int}} \mathcal F\right]$$
On the exterior (whose boundary is $-\gamma$), we have
$$\exp\left[\int_{-\gamma} \mathcal A\right] = \exp\left[\int_{\mathrm{ext}}\mathcal F\right]$$
but since $\exp\left[\int_{-\gamma}\mathcal A\right]=\exp\left[-\int_\gamma\mathcal A\right]$, it follows immediately that
$$\exp\left[\int_\mathrm{int} \mathcal F\right] \cdot \exp\left[\int_\mathrm{ext}\mathcal F\right]= \exp\left[\int_M \mathcal F \right] = 1$$
$$\implies \int_M \mathcal F = i\int_{M} F = 2\pi i \cdot n$$
$$\implies \mathrm{Ch} := \frac{1}{2\pi} \int_M F = n\in \mathbb Z$$
A: There are two different types  characteristic classes named after Chern. In index theorems one meets the numbers ${\rm ch}_n$ associated to the Chern Characters
$$
{\rm ch}[F]\equiv  {\rm tr}\left\{ \frac 1{n!} \left(\frac{F}{2\pi i}\right)^n\right\},
$$
where  $F$ is the curvature tensor of a complex vector bundle.
If we integrate over a manifold $M$ we get the numbers
$$
{\rm ch}_n[M]= \int_M {\rm tr}\left\{ \frac 1{n!} \left(\frac{F}{2\pi i}\right)^n\right\}.
$$
These are what physicist often call "Chern numbers".  They are not necessarily integers however. They can be rational numbers. The ${\rm ch}_n[M]$ are  are only guaranteed to be integers if  manifold $M$ is a spin manifold i.e one  on which one can consistently define a Dirac operator.
The other set of "Chern" objects are the Chern Classes
$$
c_1[F] = {\rm ch}_1[F]\\
c_2[F] =  \frac 12  {\rm ch}_1[F] \wedge {\rm ch}_1[F]- {\rm ch}_2[F]\\
c_3[F]= \frac 16 \left( ({\rm ch}_1[F])^3 - 6 {\rm ch}_1[F]{\rm ch}_2[F]+ 12 {\rm ch}_3[F]\right).
$$
The somewaht mysterious coefficients 6, 12 etc are derived by expressing the elementary symmetric functions of the Chern roots $z_i$ of $F$ in terms of the power symmetric functions $p_n=\sum z^n_i$. The integrals of the Chern-classes associated with a complex vector bundle over any manifold $M$ are guaranteed to be integers because they give the winding  numbers of the map from the manifold $M$ into  a subset of the  Schubert cycles in the complex Grassmanian ${\rm Gr}(n,N)$, where $n$ is the dimension of the complex vector space,  and $N$ is taken to be large enough that the result is indepenedent of $N$.
The set Schubert cycles can be labelled by Young diagrams, and so can the possible symmetric polynomials in  $F$.
