I am looking for a computationally cheap way to compute $x$ such that $$(L L^T + \mu^2 I)x = y$$ where $L \in \mathbb{R}^{n \times n}$ is a lower triangular definite positive matrix (with some very small eigenvalues), $y \in \mathbb{R}^n$ and $\mu \in \mathbb{R}$ are known. If it is necessary, I can assume that $$\mu \ll 1$$ But $\mu^2$ is larger than the smallest eigenvalue of $LL^T$.
Basically, I would like to make the most of my knowledge of the Cholesky decomposition $L L^T$. Eventually, I hope to be able to compute $x$ in $\mathcal{O}(n^2)$. Approximate approaches are also welcomed.
I have seen here that this does not seem to be doable in a more general situation, but I hope the smallness of $\mu$ may help...
Any idea, reference or warning? Thanks for your help.