This answer is based on rchilton1980 answer and the comments that have followed.
Let $K = L L^T$ be the Cholesky factorization of $K$.
Let $Y = L^{-1} Z$ , computed by forward substitution.
Let $\mathcal{I} (A)$ denote the number of negative eigenvalues of a given matrix $A$ (negative index of inertia).
We seek to compute $\mathcal{I}_\sigma = \mathcal{I} (K - \sigma M) = \mathcal{I} (L L^T - \sigma Z Z^T)$ .
Since the inertia is invariant under a change of basis (see https://en.wikipedia.org/wiki/Matrix_congruence and https://en.wikipedia.org/wiki/Sylvester%27s_law_of_inertia), inertia count $\mathcal{I}_\sigma$ can be written as
$\mathcal{I}_\sigma = \mathcal{I} \left( L^{-1} \left( L L^T - \sigma Z Z^T \right) L^{-T} \right) = \mathcal{I} \left( I_N - \sigma Y Y^T \right) = \mathcal{I} \left( \frac{1}{\sigma} I_N - Y Y^T \right)$ .
Now, it should be noted that $Y Y^T$ and $Y^T Y$ have the same nonzero eigenvalues (see for instance https://en.wikipedia.org/wiki/Singular-value_decomposition). The eigendecomposition
$\left( Y Y^T \right) \Phi_Y = \Phi_Y \begin{bmatrix} \Lambda_Y & 0\\ 0 & 0 \end{bmatrix}$ such that $Y Y^T = \Phi_Y \begin{bmatrix} \Lambda_Y & 0\\ 0 & 0 \end{bmatrix} \Phi_Y^T$ and $\Phi_Y \Phi_Y^T = I_N$ yields
$\mathcal{I}_\sigma = \mathcal{I} \left( \Phi_Y \left( \frac{1}{\sigma} I_N - \begin{bmatrix} \Lambda_Y & 0\\ 0 & 0 \end{bmatrix} \right) \Phi_Y^T \right)$ and consequently,
$\mathcal{I}_\sigma = \mathcal{I} \left( \frac{1}{\sigma} I_N - \begin{bmatrix} \Lambda_Y & 0\\ 0 & 0 \end{bmatrix} \right)$ , which is given by $\mathcal{I}_\sigma = \mathcal{I} \left( \frac{1}{\sigma} I_r - \Lambda_Y \right)$ .
Matrix $\Lambda_Y$ is the diagonal matrix of the eigenvalues of $Y^T Y$. Once $\Lambda_Y$ is calculated (independent of $\sigma$), then $\mathcal{I}_\sigma$ can be computed at no cost for any $\sigma > 0$.
Alternatively, we now show how it is possible to compute all the $r$ finite eigenvalues and associated eigenvectors of the generalized eigenvalue problem (GEP) $K \Phi = M \Phi \Lambda$ , with the same cost.
Recalling that $M = Z Z^T$, it can be deduced that the eigenvectors belong to the image (column space) of $(N \times r)$ real matrix $S = K^{-1} Z$ (see for instance Large-scale generalized eigenvalue problem with low rank LHS matrix).
Applying Rayleigh-Ritz (see https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Ritz_method) to above GEP (Galerkin projection onto $S$), a $(r \times r)$ GEP is obtained,
$(S^T K S) X = (S^T M S) X \Lambda$.
That is to say, an eigenvector $\boldsymbol{\varphi}$ of the $(N \times N)$ GEP is written as $\boldsymbol{\varphi} = S \textbf{x}$ with $\textbf{x}$ an eigenvector of the $(r \times r)$ GEP (that is, $\Phi = S X$). This procedure yields exact eigenpairs since $S$ spans the sought eigenspace.
Now, it can be observed that $S^T K S = Y^T Y$ and $S^T M S = {(Y^T Y)}^2$. Introducing the notation $A = Y^T Y$, the $(r \times r)$ GEP writes $A X = A^2 X \Lambda$, which is equivalent to the standard eigenvalue problem (SEP) $A X = X \Lambda^{-1}$ . It can be recognized that $\Lambda = \Lambda_Y^{-1}$.
As a conclusion, the proposed approach for inertia count is as costly as directly computing all the eigenvalues of the associated GEP. The cost (dense eigenvalue problem with dimension $r$) rapidly increases with $r$.
It should be noted that the cost for computing $\mathcal{I}_\sigma$ for one given $\sigma > 0$ can be less than that of computing all the eigenvalues, as one also has $\mathcal{I}_\sigma = \mathcal{I} \left( \frac{1}{\sigma} I_r - Y^T Y \right)$, which can be obtained by LDL factorization of $\frac{1}{\sigma} I_r - Y^T Y$.
My next question would be how to determine the inertia count of $K_\sigma - \sigma Z Z^T$ with $K_\sigma = K - \sigma M_0$ , in which $M_0$ is a sparse $(N \times N)$ symmetric positive-semidefinite matrix. The difference with the above is that $K_\sigma$ is not positive definite. I am presently using the above developments to answer this question.
Edit: See below for the answer to this question.
We seek to compute $\mathcal{I}_\sigma = \mathcal{I} \left( K - \sigma M \right)$ with $M = M_0 + Z Z^T$.
We have $\mathcal{I}_\sigma = \mathcal{I} \left( K_\sigma - \sigma Z Z^T \right)$ with $K_\sigma = K - \sigma M_0$ a full-rank sparse symmetric matrix.
The LDL factorization of $K_\sigma$ is performed, which yields $K_\sigma = L_\sigma D_\sigma L_\sigma^T$.
Thus, we obtain $\mathcal{I}_\sigma = \mathcal{I} \left( L_\sigma^{-1} \left( K_\sigma - \sigma Z Z^T \right) L_\sigma^{-T} \right) = \mathcal{I} \left( D_\sigma - \sigma Y_\sigma Y_\sigma^T \right)$ with $Y_\sigma = L_\sigma^{-1} Z$. Which gives $\mathcal{I}_\sigma = \left( \frac{D_\sigma}{\sigma} - Y_\sigma Y_\sigma^T \right)$.
Introducing the $(N \times N)$ GEP $\left( \frac{D_\sigma}{\sigma} \right) \Phi = \left( Y_\sigma Y_\sigma^T \right) \Phi \Lambda$ such that $\Phi^T \left( \frac{D_\sigma}{\sigma} \right) \Phi = \Lambda$ and $\Phi^T \left( Y_\sigma Y_\sigma^T \right) \Phi = I_N $ (unit mass eigenvector normalization), we obtain $\mathcal{I}_\sigma = \mathcal{I} \left( \Phi^T \left( \frac{D_\sigma}{\sigma} - Y_\sigma Y_\sigma^T \right) \Phi \right) = \mathcal{I} \left( \Lambda - I_N \right)$. This GEP cannot be solved.
Similarly to before, the $r$ eigenvectors associated with the $r$ finite eigenvalues belong to the column space of $S_\sigma = {\left( \frac{D\sigma}{\sigma} \right)}^{-1} Y_\sigma$. The $(r \times r)$ GEP $\left( S_\sigma^T \left( \frac{D_\sigma}{\sigma} \right) S_\sigma \right) X = \left( S_\sigma^T \left( Y_\sigma Y_\sigma^T \right) S_\sigma \right) X \Lambda_f$ allows to obtain the finite eigenvalues $\Lambda_f$ of the $(N \times N)$ GEP. Let $A_\sigma = S_\sigma^T \left( \frac{D_\sigma}{\sigma} \right) S_\sigma$. One has $A_\sigma = Y_\sigma^T {\left( \frac{D_\sigma}{\sigma} \right)}^{-1} Y_\sigma$ and $S_\sigma^T \left( Y_\sigma Y_\sigma^T \right) S_\sigma = A_\sigma^2$. Thus, the $r$ finite eigenvalues $\Lambda_f$ can be deduced from the $(r \times r)$ SEP $A_\sigma X = X \Lambda_f^{-1}$.
Let $\Lambda_i$ denote the matrix of the $N - r$ infinite eigenvalues.
We have $\mathcal{I}_\sigma = \mathcal{I} \left( \Lambda - I_N \right) = \mathcal{I} \left( \Lambda_f - I_r \right) + \mathcal{I} \left( \Lambda_i - I_{N - r} \right)$. Let $\mathcal{I}_f = \mathcal{I} \left( \Lambda_f - I_r \right)$ and $\mathcal{I}_i = \mathcal{I} \left( \Lambda_i - I_{N - r} \right)$. We have $\mathcal{I}_f = \mathcal{I} \left( X^T A_\sigma^{-1} X - X^T X \right)$ because $X^{-1} = X^T$, and consequently, $\mathcal{I}_f = \mathcal{I} \left( A_\sigma^{-1} - I_r \right)$. Since for a real scalar $a$ we have $a^{-1} \in ] - \infty , 1 [ \,\,\, \Rightarrow \,\, a \in ] - \infty , 0 [ \, \cup \,] 1 , + \infty [\,$, it can be seen that $\mathcal{I}_f = \mathcal{I} \left( A_\sigma^{-1} - I_r \right) = \mathcal{I} \left( A_\sigma \right) + \mathcal{I} \left( I_r - A_\sigma \right)$.
On the other hand, we have $\mathcal{I}_i = \mathcal{I} \left( \Lambda_i - I_{N - r} \right) = \mathcal{I} \left( \Lambda_i \right)$ ($\mathcal{I}_i$ is given by the number of negative infinite eigenvalues). Since $\mathcal{I} \left( \Lambda \right) = \mathcal{I} \left( \Lambda_f \right) + \mathcal{I} \left( \Lambda_i \right)$ and $\mathcal{I} \left( \Lambda \right) = \mathcal{I} \left( \Phi^T \left( \frac{D_\sigma}{\sigma} \right) \Phi \right) = \mathcal{I} \left( \frac{D_\sigma}{\sigma} \right) = \mathcal{I} \left( D_\sigma \right)$, we obtain $\mathcal{I}_i = \mathcal{I} \left( D_\sigma \right) - \mathcal{I} \left( \Lambda_f \right)$. In addition, we have $\mathcal{I} \left( \Lambda_f \right) = \mathcal{I} \left( \Lambda_f^{-1} \right) = \mathcal{I} \left( A_\sigma \right)$. Therefore, we obtain $\mathcal{I}_i = \mathcal{I} \left( D_\sigma \right) - \mathcal{I} \left( A_\sigma \right)$.
As a conclusion, it can be deduced that $\mathcal{I}_\sigma = \mathcal{I} \left( D_\sigma \right) + \mathcal{I} \left( I_r - A_\sigma \right)$.
Summary of the operations:
LDL factorization of $(N \times N)$ sparse matrix $K_\sigma = L_\sigma D_\sigma L_\sigma^T$
Solve $Y_\sigma = L_\sigma^{-1} Z$ by forward substitutions
Assemble $A_\sigma = Y_\sigma^T {\left( \frac{D_\sigma}{\sigma} \right)}^{-1} Y_\sigma$
LDL factorization of $(r \times r)$ dense matrix $ I_r - A_\sigma$
