Let $G = (V, E)$ be a finite graph and let $\Delta$ denote its Laplacian (here we take the convention that the Laplacian is positive semidefinite). For simplicity name the vertices $1, 2, ... n$. Then $\Delta$ is the matrix associated to the quadratic form
$$q(x_1, ... x_n) = \sum_{(i, j) \in E} (x_i - x_j)^2$$
which one might call the Dirichlet energy. Recall that, for an appropriate change of variables, such a quadratic form can be written as
$$q(z_1, ... z_n) = \sum_i \lambda_i z_i^2$$
where $\lambda_i$ are the eigenvalues of $\Delta$. Thus to determine the multiplicity of the zero eigenvalue it suffices to determine when $q$ can be zero. But from the first expression it should be clear that $q = 0$ if and only if $x_i = x_j$ whenever $(i, j) \in E$; that is, whenever $(x_1, ... x_n)$ determines a function
$$x : V \to \mathbb{R}$$
which is constant on each connected component of $G$. The dimension of this space of functions is clearly the number of connected components of $G$, and the conclusion follows.
What does it mean "orthogonal eigenvectors of eigenvalue zero"? Are those eigenvalues really equal 0?
Eigenvector for eigenvalue zero is a vector $\mathbf{a}$ such that $L\mathbf{a}=\mathbf{0}$. In this case
$$\begin{bmatrix}
L_{G_{1}} & 0\\
0 & L_{G_{2}}
\end{bmatrix}
\begin{bmatrix}
0\\
1
\end{bmatrix}
=
\begin{bmatrix}
0\\
0
\end{bmatrix}
\text{ and }
\begin{bmatrix}
L_{G_{1}} & 0\\
0 & L_{G_{2}}
\end{bmatrix}
\begin{bmatrix}
1\\
0
\end{bmatrix}=
\begin{bmatrix}
0\\
0
\end{bmatrix}$$
hence the vectors $\mathbf{a}_1=[1,\dots,1,0,\dots,0]^T$ and $\mathbf{a}_2=[0,\dots,0,1,\dots,1]^T$ are eigenvectors corresponding to eigenvalue 0.
They are orthogonal since $\mathbf{a}_1^T\mathbf{a}_2=0$. Since we know that these vectors are orthogonal, they must be linearly independent.
How do we know that these eigenvalues correspond to $\lambda_2$?
Since the eigenvalues are ordered increasingly and we have two linearly independent eigenvectors, i.e. at least two eigenvalues are zero, we know that $\lambda_1=\lambda_2=0$.
I suspect that all $\lambda$'s are ordered by multiplicity. But how can I prove that every n by n symmetric matrix has n eigenvalues, counted with multiplicity.
Eigenvalues are roots of characteristic polynomial, which has degree $n$. So, counting multiplicity, there is $n$ of them. Here we are speaking about algebraic multiplicity, see wikipedia.
Why actually it's true if $u$ and $v$ are connected by an edge, we have $x(u)=x(v)$, it's not obvious for me, even if $x$ is eigenvector, and what happens in the general case?
Since $\sum\limits_{(u,v)\in E}^{} (x(u)-x(v))^{2} = 0$ is a sum of non-negative real numbers, each of them must be zero. So we have $(x(u)-x(v))^{2}=0$ and, consequently, $(x(u)-x(v))=0$, for each $(u,v)\in E$.
OK, by now we prove that $\lambda_{1}=0$, but what about $\lambda_{2}$?
It was shown in the preceding part of the proof, that (if $G$ is connected) every eigenvector for zero is multiple of $[1,1,\dots,1]^T$. This shows that eigenspace corresponding to zero eigenvalue has dimension 1, i.e., the geometric multiplicity of this eigenvalue is 1.
EDIT: What I wrote here before was unnecessary complicated. Laplacian matrix is symmetric, hence it is diagonalizable. Thus algebraic and geometric multiplicity are the same.
However, we need to show that algebraic multiplicity of $\lambda_1=0$ is 1. Suppose this would be not true. It would imply that there exists a vector $\mathbf a$ such that $L \mathbf a = [1,1,\dots,1]^T$. This would imply that vector $[1,1,\dots,1]^T$ belongs to the column space of the matrix $L$. But directly from the definition of of Laplacian matrix we can see that this vector is orthogonal to each column of $L$. (At least for non-weighted graphs - I did not think about the weighted case.)
I agree with you that this last argument seems to be missing in the pdf-file you linked. (And probably different explanations are possible.)
Just a few words of explanation. (But it is better to read it with more details at wiki or in some linear algebra textbook.)
If $\mathbf{a}$ is an eigenvector of a matrix $A$, it means that $(A-\lambda I)\mathbf a=\mathbf 0$. If algebraic multiplicity (=the multiplicity of $\lambda$ as a root of characteristic polynomial) is greater than the geometric multiplicity (=dimension of the corresponding eigenspace), it implies that there is an eigenvector $\mathbf a$ and a generalized eigenvector $\mathbf b$ such that $(A-\lambda I)\mathbf b=\mathbf a$ .
This was the property which we used in the last paragraph.
Geometric multiplicity is precisely the number of Jordan blocks in the Jordan normal form.
It might be worth mentioning, that in the book Bapat: Graphs and Matrices, Universitext, Springer 2010, the Laplacian matrix is defined as follows.
The author first defines
the (vertex-edge) incidence matrix of G denoted by $Q(G)$, which is the $n\times m$ matrix with entries $0$ if the given vertex and edge are incident, $1$ if the edge starts in this vertex and $-1$ if it ends there. It is relatively easy to show that rank of $Q(G)$ is $n-k$, where $k$ is the number of connected components of this graph, see p.11-12 of this book.
Then the Laplacian matrix can be obtained as $L(G)=Q(G)Q(G)^T$, see p.45
Using this fact we get $\operatorname{rank}(L)=\operatorname{rank}(QQ^T)=\operatorname{rank}(Q)$, see Lemma 4.3, p.46.
The fact that matrix $Q$ and its Gram matrix $QQ^T$ have the same rank is mentioned in the same book on p.10
as an exercise:
Let $A$ be an $m\times n$ matrix. Show that $A$ and $A^TA$ have the same null space. Hence
conclude that $\operatorname{rank} A = \operatorname{rank}A^TA$.
I.e., it was shown that $n-\operatorname{rank}(G)$ is the number of connected components.
Since $L(G)$ is symmetric, it is diagonalizable, and thus this is the same as the number of non-zero eigenvalues.
Again, I do not know whether something similar can be done in the weighted case.
Best Answer
In general, the idea with upper bounds on Laplacian eigenvalues is to use the following characterization: if $\lambda_1(L) \le \dots \le \lambda_n(L)$ are the eigenvalues of a symmetric matrix $L$, then $$ \lambda_1(L) = \inf\left\{\frac{\mathbf x^{\mathsf T}L\mathbf x}{\mathbf x^{\mathsf T} \mathbf x}: \mathbf x \ne \mathbf 0\right\} $$ and if $\mathbf u^{(1)}$ is the eigenvector corresponding to $\lambda_1(L)$, then $$ \lambda_2(L) = \inf\left\{\frac{\mathbf x^{\mathsf T}L\mathbf x}{\mathbf x^{\mathsf T} \mathbf x}: \mathbf x \ne \mathbf 0, \mathbf x^{\mathsf T}\mathbf u^{(1)} = 0\right\}. $$ The minimum in both cases is achieved by taking $\mathbf x$ to be the corresponding eigenvector.
If we pick a specific value of $\mathbf x$ (satisfying $\mathbf x \ne \mathbf 0, \mathbf x^{\mathsf T}\mathbf u^{(1)} = 0$) that we think is going to be a decent approximation of the eigenvector to $\lambda_2(L)$, and compute $\frac{\mathbf x^{\mathsf T}L\mathbf x}{\mathbf x^{\mathsf T} \mathbf x}$, that gives us an upper bound on the actual value of $\lambda_2(L)$.
In the case where $L$ is the Laplacian matrix of a graph, $\lambda_1 = 0$ and the corresponding $\mathbf u^{(1)}$ is the all-$1$ vector, which makes the optimization problem for $\lambda_2$ easier to state. Being orthogonal to $\mathbf u^{(1)}$ just means all entries of the eigenvector have to add up to $0$.
The bound that $\lambda_2(G) \le \lambda_2(G-S) + |S|$ is obtained in exactly this way. Let $G$ be an $n$-vertex graph, and let $\mathbf x \in \mathbb R^{n-|S|}$ be the second eigenvector of the Laplacian matrix of $G-S$. Then the vector $\mathbf y \in \mathbb R^n$ obtained by padding $\mathbf x$ with $|S|$ zeroes is a good guess at the eigenvector of $\lambda_2$ for the Laplacian matrix of $G$. (And it satisfies the orthogonality condition because the sum of $\mathbf y$'s entries is the same as the sum of $\mathbf x$'s entries, which is $0$.)
Knowing that $$ \frac{\mathbf x^{\mathsf T}L_{G-S}\mathbf x}{\mathbf x^{\mathsf T}\mathbf x} = \lambda_2(G-S) $$ we first want to estimate $$ \frac{\mathbf y^{\mathsf T}L_{G}\mathbf y}{\mathbf y^{\mathsf T}\mathbf y}. $$ In the denominator, $\mathbf y^{\mathsf T}\mathbf y = \mathbf x^{\mathsf T}\mathbf x$ because $\mathbf y$ is just $\mathbf x$ with some zero entries. In the numerator, things are a bit more complicated.
The zero entries of $\mathbf y$ means that the extra rows and columns of $L_G$ that weren't in $L_{G-S}$ (corresponding to vertices in $S$) go away, and $L_G$ mostly agrees with $L_{G-S}$ on the rows and columns they do share. The exception is the diagonal: the diagonal entries of $L_{G-S}$ measure degree in $G-S$, and the diagonal entries of $L_G$ measure degree in $G$. So when we go from $\mathbf x^{\mathsf T}L_{G-S}\mathbf x$ to $\mathbf y^{\mathsf T}L_{G}\mathbf y$, we gain an extra $x_i^2$ (or $y_i^2$) for every extra edge that vertex $i$ had to $S$: increasing the degree of vertex $i$ in $G$ beyond what it had in $G-S$.
But there are at most $|S|$ such edges, and so the extra contribution is at most $|S|x_i^2$ for vertex $i$, and $|S|\mathbf x^{\mathsf T}\mathbf x$ for all vertices of $G-S$ combined. Therefore $$ \frac{\mathbf y^{\mathsf T}L_{G}\mathbf y}{\mathbf y^{\mathsf T}\mathbf y} \le \frac{\mathbf x^{\mathsf T}L_{G-S}\mathbf x + |S|\mathbf x^{\mathsf T}\mathbf x}{\mathbf x^{\mathsf T}\mathbf x} = \frac{\mathbf x^{\mathsf T}L_{G-S}\mathbf x}{\mathbf x^{\mathsf T}\mathbf x} + |S|. $$ Since $\frac{\mathbf y^{\mathsf T}L_{G}\mathbf y}{\mathbf y^{\mathsf T}\mathbf y}$ is an upper bound on $\lambda_2(G)$, and $\frac{\mathbf x^{\mathsf T}L_{G-S}\mathbf x}{\mathbf x^{\mathsf T}\mathbf x}$ is exactly $\lambda_2(G-S)$, we get $$ \lambda_2(G) \le \lambda_2(G-S) + |S|. $$ In particular, if $S$ is a vertex cut for $G$, then $\lambda_2(G-S)=0$, so $\lambda_2(G) \le |S|$.