Below, I provide a "high-probability" non-asymptotic bound (see (+) below) based on non-linear Berry-Esseen theory developed by Iosif Pinelis. I'd be grateful if someone would kindly check that I didn't screw up anything. Thanks in advance!
Main tool: non-linear Berry-Esseen theory
Let $ \mathcal H$ be a (possibly infinite-dimensional) Hilbert space with topological dual $\mathcal H^*$, and let $g:\mathcal H \rightarrow \mathbb R$ be a function such that
- $g(0)=0$,
- $g$ has linear approximant $L \in \mathcal H^*$ valid on a nonzero neighborhood of the origin, i.e there exist $\epsilon > 0,M>0$ such that
$$
|g(z) - L(z)| \le \frac{M}{2}\|z\|^2,\;\forall \|z\| \le \epsilon. \tag{1}
$$
Note that under such conditions, we automatically have that $g$ is (Fréchet) differentable at $0$ and $L=\nabla g(0)$.
Let $Z$ be a random vector on $\mathcal H$ such that $E Z = 0$ (i.e $Z$ is centered) and suppose $\tilde{\sigma} := \|L(Z)\|_{Z,2} := (E_Z|L(Z)|^2)^{1/2} < \infty$. The quantity $\tilde{\sigma}^2/k$ will play the role of a "proxy variance" for the random variable $g(Z)$. Finally, let $Z_1,\ldots,Z_k$ be iid copies of $Z$ and set $\overline{V} := (1/k)\sum_{i=1}^kZ_i$.
Theorem (Corollary 3.7 of reference paper). Let $p \in (2,3]$ such that $\|Z\|_{Z,p} := (E_Z\|Z\|_p)^{1/p} < \infty$. Then
$$
\sup_{t \in \mathbb R}\left|\mathbb P\left(\frac{g(\overline{V})}{\tilde{\sigma}/\sqrt{k}}\ge t\right)-\Phi(t)\right| \le Ck^{1-p/2},
$$
where $C$ is a constant which only depends on the distribution of $Z$ (in particular, $C$ is independent of the sample size $k$ and $p$).
Application: our problem
For our own business, we let
$\mathcal H$ be the euclidean space of $n$ by $n$ matrices equipped with the Frobenius trace inner product $\langle B,\tilde{B}\rangle_{Fro} := \text{trace}(B^T\tilde{B})$. We denote induced norm by $\|B\|_{Fro} := \langle B,B\rangle_{Fro}^{1/2}$. Note that this space is nothing but euclidean $\mathbb R^{n^2}$ in disguise.
Define
$$g(B) := \begin{cases}f(B + A + \sigma^2 I_n) - f(A+\sigma^2 I_n),&\mbox{ if }B + A + \sigma^2 I_n \succ 0,\\0,&\mbox{ else.}\end{cases}
$$
where $f(B) := \sum_{j=1}^n\log\lambda_j(B)$.
It's not hard to see that $g$ is differentiable at $0$ with derivative $L:=g'(0) := (A+\sigma^2 I_n)^{-1}$. Moreover, $g$ has Lipschitz continuous gradient at $0$ and and so (1) is satisfied.
For our random vector, we take $Z = RR^T-\sigma^2 I_n$, where $R \in \mathbb R^n$ is a random vector with iid entries from $N(0,\sigma^2 I_n)$. Thus $XX^T-\sigma^2 I_n = (1/k)\sum_{i=1}^k Z_j =: \overline{V}$, where $Z_j = X_jX_j^T$ are iid copies of $Z$ (where $X_i$ is the $i$th column of the random $n$ by $k$ matrix $X$ in our original problem). Of course, $\|Z\|_{Z,p} < \infty$ for all $p$ (Gaussians have finite moments).
By Cauchy-Schwarz inequality, one computes
$$
\begin{split}
\tilde{\sigma}^2 &= \|L(Z)\|_Z = \|L\|_{Fro}^2\|Z\|_{Z,2}^2 = \text{trace}((A+\sigma^2 I_n)^{-2})E_Z\|Z\|_2^2\\
& \le (\sum_{j=1}^n(\lambda_j(A) + \sigma^2)^{-1})^2\cdot n \sigma^4 \le n\sigma^4\eta(A)^2 = (\sqrt{n}\sigma^2\eta(A))^2,
\end{split}
$$
where
$$
\begin{split}
\eta(A) &= \eta(A; \sigma^2) := \text{trace}(A+\sigma^2 I_n)^{-1} = \sum_{i=1}^n(\lambda_i(A)+\sigma^2)^{-1}
\le n\min(\sigma^{-2},\lambda_n(A)^{-1}).
\end{split}
$$
Putting everything together, and applying the above theorem, we obtain that for every $t \ge 0$, w.p at least $\Phi(t\sqrt{k})+\mathcal O(k^{-1/2})$,
$$
|f(A+XX^T) - f(A + \sigma^2 I_n)| = |f(\overline{V}+A+\sigma^2 I_n) - f(A + \sigma^2 I_n)| \le t\sqrt{n}\sigma^2\eta(A).
$$
Combining with (*) then gives: for all $t \ge 0$, it holds w.p at least $\Phi(t\sqrt{k/n})+\mathcal O(k^{-1/2})$ that
$$
|f(A+XX^T) - f(A)| \le \sigma^2\eta(A) + t\sigma^2\eta(A) = (1+t)\sigma^2\eta(A). \tag{+}
$$
Edit: Improved non-asymptotic bound for the term $|f(A+XX^T)-f(A+\sigma^2 I_n)|$
It has been noted in the comments section of the question that Berry-Esseen sounds like an over-kill (and might produce sub-optimal bounds for my particular problem in which everything is sub-Gaussian...). To remedy this, I'll use well-known concentration results for estimation of covariance matrices, and then Weyl's inequality. Interestingly, just as in the Berry-Esseen method above, the key requirement is that the entries of $X$ be iid with finite moments of order $2 + \epsilon$.
By Proposition 2.1 of HOW CLOSE IS THE SAMPLE COVARIANCE MATRIX TO
THE ACTUAL COVARIANCE MATRIX?, we know that for any $t \ge 0$, it holds with probability $1-2e^{-ct^2}$ that
$$
\|XX^T-\sigma^2I\| = \left\|\frac{1}{k}\sum_{j=1}^kR_iR_i^T-\sigma^2 I_n\right\| \le t\sigma^2\sqrt{\frac{n}{k}}.
$$
Thus by Weyl's inequality, it holds with probability $1-2e^{-ckt^2/n}$ that
$$
\max_{1 \le i \le n} |\lambda_i(A+XX^T)-\lambda_i(A+\sigma^2)| \le \|XX^T-\sigma^2I_n\| = \left\|\frac{1}{k}\sum_{j=1}^kR_iR_i^T-\sigma^2 I_n\right\| \le t\sigma^2,
$$
and so by the triangle inequality, it holds with probability $1-2e^{-ckt^2/n}$ that
$$
\begin{split}
|f(A+XX^T)-f(A+\sigma^2 I_n)| &= \left| \sum_{i=1}^n \log\left(\frac{\lambda_i(A+XX^T)}{\lambda_i(A + \sigma^2 I_n)}\right)\right|\\
& \le \sum_{i=1}^n \left|\log\left(\frac{\lambda_i(A+XX^T)}{\lambda_i(A + \sigma^2 I_n)}\right)\right| \\
& \le \sum_{i=1}^n\log\left(1 + \frac{t\sigma^2}{\lambda_i(A)+\sigma^2}\right)
\\
&\le t\sigma^2\sum_{i=1}^n (\lambda_i(A)+\sigma^2)^{-1}\\
&= t\sigma^2\eta(A).
\end{split}
$$
