$$ |f(A+XX^T) - f(A)| \le \sigma^2\eta(A) + t\sigma^2\eta(1) = (1+t)\sigma^2\eta(A). \tag{+} $$

$$ |f(A+XX^T) - f(A)| \le \sigma^2\eta(A) + t\sigma^2\eta(A) = (1+t)\sigma^2\eta(A). \tag{+} $$

• $\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 \sum_{i=1}^n (\lambda_i(A))^{-1} \le \frac{n}{\lambda_n(A)}. \end{split} $$

• $\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} $$

Improved non-asymptotic tail bound for $|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$.

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$.