I found somewhere the following formula for the quantum Fisher information of a multimode Gaussian state $(\mathbf d_\theta, \sigma_\theta)$: $$\tag{1} H(\theta) = \frac{1}{2} \text{Tr}\left[ \left( \sigma_\theta^{-1} \partial_\theta \sigma_\theta \right)^2 \right] + \left( \partial_\theta \mathbf{d}_\theta \right)^\top \sigma_\theta^{-1} \left( \partial_\theta \mathbf{d}_\theta \right). $$ I'm far from an expert on this topic, but it seems that other sources list seemingly much more complicated expressions. For instance, this paper recommends the formula for the QFIM $$ H_{ij}(\boldsymbol\theta) = \sum_{k,l=1}^{N} \left[ \frac{(\lambda_k - \lambda_l)^2}{\lambda_k \lambda_l - 1} \text{Re}(R^i_{kl} R^j_{kl}) + \frac{(\lambda_k + \lambda_l)^2}{\lambda_k \lambda_l + 1} \text{Re}(Q^i_{kl} Q^j_{kl}) \right] + \sum_{k=1}^{N} \frac{\partial_i \lambda_k \partial_j \lambda_k}{\lambda_k^2 - 1} + 2 \partial_i d^\dagger \sigma^{-1} \partial_j d. $$ in terms of the symplectic eigenvalues of the covariance matrix and certain matrices $R_i, Q_i$ computed from the diagonalizing transformation. Alternatively, the formula $$ H_{ij}(\boldsymbol\theta) = \frac{1}{2} \text{vec}[\partial_i \sigma]^\dagger M^{-1} \text{vec}[\partial_j \sigma] + 2 \partial_i d^\dagger \sigma^{-1} \partial_j d, $$ for $M = \bar\sigma\otimes\sigma - \Omega \otimes \Omega$ given the symplectic form $\Omega$. This formula is only valid as stated when none of the modes is pure; in order to generalized to any Gaussian state (which is the case I am interested in) one needs to first compute the QFI for the covariance matrix $\nu\sigma_\theta$, where $\nu>0$ is a regularizing parameter, and then take the limit for $\nu\to 1$.
Even in the case of a single parameter, which is the one I am interested in, it is not really obvious to me how my original formula related to these general ones. The second term looks close enough, assuming the displacement vector is real and ignoring for a moment the factor of two, but I don't quite believe the first term is equivalent. When is $(1)$ valid, if it is ever valid at all?
I should add that my goal is to write some code to compute the QFI given any pair $(d,\sigma)$ and their derivatives, computed for a specific parametric Gaussian channel. Numerical efficiency is a concern for me, which is why both the Williamson's decomposition approach and the limiting procedure look quite unappealing.
