# MBD expressions

I present three different yet equivalent expressions for the MBD energy (MBD1–MBD3 below). By MBD energy, I mean the interaction energy of a system of Drude oscillators under dipole approximation. (In our Chemical Reviews paper and my thesis, I started from the fluctuation–dissipation theorem for the electrons, but I start directly with oscillators here to simplify the discussion.)

First two general formulas that I'll use. For any Hamiltonian of the form $\hat H=\hat T+\hat U+\hat V$ (kinetic energy, one-particle potential, particle interaction), the interaction energy can be obtained via adiabatic connection by generalizing the Hamiltonian as $\hat H(\lambda)=\hat T+\hat U+\lambda\hat V$ and using the Hellmann–Feynman theorem, $E_\text{int}=E(1)-E(0)=\int_0^1\mathrm d\lambda\frac{\partial E}{\partial\lambda}=\int_0^1\mathrm d\lambda\langle\Psi(\lambda)|\hat V|\Psi(\lambda)\rangle \label{eq:adiabatic-connection}$ The zero-temperature fluctuation–dissipartion theorem for multiple quantities, $\hat A_i$, connects the correlations in their spontaneous fluctuations to the corresponding generalized susceptibilities, $\partial\langle\hat A_i\rangle/\partial f_j(u)$, which desribe response to the generalized external forces $f_j$ at freqeuncy $u$, $\langle\Psi|(\hat A_i-\langle\hat A_i\rangle)(\hat A_j-\langle\hat A_j\rangle)|\Psi\rangle=-\frac1\pi\int_0^\infty\mathrm du\operatorname{Im}\frac{\partial\langle\hat A_i\rangle}{\partial f_j(u)}\equiv-\frac1\pi\int_0^\infty\mathrm du\frac{\partial\langle\hat A_i\rangle}{\partial f_j(\mathrm iu)} \label{eq:fluctuation-dissipation}$ Now going to MBD, the Hamiltonian is $\hat H_\text{MBD}=\hat T+\hat U+\hat V=\sum_i\frac{\mathbf{\hat p}_i^2}2+\sum_i\frac12m_i\omega_i^2|\mathbf{\hat r}_i-\mathbf R_i|^2+\sum_{i<j}\boldsymbol{\hat\mu}_i\cdot\mathbf T_{ij}\boldsymbol{\hat\mu}_j$ where $\mathbf p_i$ is momentum, $m_i$ is mass, $\omega_i$ is frequency, $\mathbf r_i$ is position of the Drude particle, $\mathbf R_i$ is the center of the harmonic potential and (fixed) position of the compensating charge, $\boldsymbol\mu_i=q_i(\mathbf r_i-\mathbf R_i)$ is the dipole moment, and $\mathbf T_{ij}$ is the Gaussian-screened dipole–dipole interaction matrix, $\mathbf T_{ij}=\boldsymbol\nabla_i\otimes\boldsymbol\nabla_jv_\text{GG}(|\mathbf R_i-\mathbf R_j|)$ Introducing generalized coordinates, $\boldsymbol\xi_i=\sqrt{m_i}(\mathbf r_i-\mathbf R_i)$, and using the expression for the static dipole polarizability of a Drude oscillator, $\alpha(u)=q^2/m(\omega^2-u^2)$, the Hamiltonian can be rewritten as \begin{aligned} \hat H_\text{MBD}&=\sum_i\frac{\mathbf{\hat p}_i^2}2+\sum_i\frac12\omega_i^2\hat\xi_i^2+\sum_{i<j}\sqrt{\alpha_i(0)\alpha_j(0)}\omega_i\omega_j\boldsymbol{\hat\xi}_i\cdot\mathbf T_{ij}\boldsymbol{\hat\xi}_j \\ &=\sum_i\frac{\mathbf{\hat p}_i^2}2+\frac12\sum_{ij}\boldsymbol{\hat\xi}_i\cdot\mathbf Q_{ij}\boldsymbol{\hat\xi}_j \end{aligned} where we have introduced $\mathbf Q_{ij}=\delta_{ij}\omega_i^2\mathbf I+\sqrt{\alpha_i(0)\alpha_j(0)}\omega_i\omega_j\mathbf T_{ij}$ and $\sum_{i\neq j}$ was translated into $\mathbf T_{ij}=0$. The Hamiltonian can now be solved exactly by diagonalizing the $3N$-by-$3N$ real symmetric matrix $\mathbf Q$, $\mathbf Q\equiv\mathbf C\operatorname{diag}(\{\tilde\omega_k^2\})\mathbf C^\text T$, which leads to new coupled momenta and generalized coordinates, $\mathbf{\tilde p}=\mathbf C\mathbf p$, $\boldsymbol{\tilde\xi}=\mathbf C^\text T\boldsymbol\xi$, in which the Hamiltonian is uncoupled and the ground-state energy is simply $\sum_{k=1}^{3N}\tilde\omega_k/2$. The MBD energy is therefore $E_\text{MBD}=\sum_{k=1}^{3N}\frac{\tilde\omega_k}2-3\sum_{i=1}^N\frac{\omega_i}2\equiv\frac12\operatorname{Tr}\Big(\sqrt{\mathbf Q(\{\alpha_i(0),\omega_i\},\mathbf T)}\Big)-3\sum_{i=1}^N\frac{\omega_i}2 \tag{MBD1}\label{eq:MBD1}$ Obviously, there is a problem when $\mathbf Q$ is not positive-definite and its square root does not exist. This signifies a failure of the dipole approximation.

The fluctuation–dissipation theorem in $\eqref{eq:fluctuation-dissipation}$ can be applied to the MBD problem by identifying $A_i$ with the dipole moments $\boldsymbol\mu_i$ and the generalized susceptibility with the dipole–dipole polarizability $\boldsymbol\alpha_{ij}(u)\equiv-\partial\langle\boldsymbol\mu_i\rangle/\partial\mathbf E_{\text{ext},j}(u)$, where $\mathbf E_{\text{ext},j}$ is the external electric field on the $j$-th atom. Since $\langle\boldsymbol\mu_i\rangle=0$ in the unperturbed case, we have $\langle\Psi|\boldsymbol{\hat\mu}_i\otimes\boldsymbol{\hat\mu}_j|\Psi\rangle=\frac1\pi\int_0^\infty\mathrm du\,\boldsymbol\alpha_{ij}(\mathrm iu)$ From the adiabatic connection in $\eqref{eq:adiabatic-connection}$, we then have ($\operatorname{Tr}_{\mathbb R^3}$ denotes trace over the three Cartesian coordinates, $x$, $y$, $z$.) \begin{aligned} E_\text{MBD}&=\frac12\int_0^1\mathrm d\lambda\langle\Psi(\lambda)\big|\sum_{i\neq j}\boldsymbol{\hat\mu}_i\cdot\mathbf T_{ij}\boldsymbol{\hat\mu}_j\big|\Psi(\lambda)\rangle \\ &=\frac12\int_0^1\mathrm d\lambda\sum_{i\neq j}\operatorname{Tr}_{\mathbb R^3\!}\big(\langle\Psi(\lambda)|\boldsymbol{\hat\mu}_i\otimes\boldsymbol{\hat\mu}_j|\Psi(\lambda)\rangle\mathbf T_{ij}\big) \\ &=\frac1{2\pi}\int_0^1\mathrm d\lambda\sum_{i\neq j}\operatorname{Tr}_{\mathbb R^3\!}\big(\textstyle\int_0^\infty\mathrm du\,\boldsymbol\alpha_{ij}(\mathrm iu;\lambda)\mathbf T_{ij}\big) \\ &=\frac1{2\pi}\int_0^1\mathrm d\lambda\int_0^\infty\mathrm du\operatorname{Tr}_{i,\mathbb R^3\!}\big(\boldsymbol\alpha(\mathrm iu;\lambda)\mathbf T\big) \end{aligned} \tag{MBD2}\label{eq:MBD2} The interacting nonlocal polarizability can be expressed with a self-consistent equation ($\boldsymbol\alpha_{ij}(u;\lambda=0)\equiv\boldsymbol\alpha_{ij}(u)=\alpha_i(u)\delta_{ij}\mathbf I$), \begin{aligned} \boldsymbol\alpha(u;\lambda)&=\boldsymbol\alpha(u)-\boldsymbol\alpha(u)\lambda\mathbf T\boldsymbol\alpha(u;\lambda) \\ &=\big(\boldsymbol\alpha(u)^{-1}+\lambda\mathbf T\big)^{-1} \\ \end{aligned} The $\lambda$-integration can now be performed analytically, $\int_0^1\mathrm d\lambda\,\boldsymbol\alpha(u;\lambda)=\ln\big(\mathbf I+\boldsymbol\alpha(u)\mathbf T\big)\mathbf T^{-1}$ Plugging into $\eqref{eq:MBD2}$, we have \begin{aligned} E_\text{MBD}&=\frac1{2\pi}\int_0^\infty\mathrm du\operatorname{Tr}_{i,\mathbb R^3}\ln\big(\mathbf I+\boldsymbol\alpha(\mathrm iu)\mathbf T\big) \\ &=\frac1{2\pi}\int_0^\infty\mathrm du\ln\det_{i,\mathbb R^3}\big(\mathbf I+\boldsymbol\alpha(\mathrm iu)\mathbf T\big) \end{aligned} \tag{MBD3}\label{eq:MBD3} Finally, I show that $\eqref{eq:MBD3}$ is equivalent to $\eqref{eq:MBD1}$, as it should, by performing the $u$-integration analytically ($\boldsymbol\omega_{ij}=\omega_i\delta_{ij}\mathbf I$), \begin{aligned} E_\text{MBD}&=\frac1{2\pi}\int_0^\infty\mathrm du\operatorname{Tr}_{i,\mathbb R^3}\ln\big(1+\boldsymbol\alpha(\mathrm iu)\mathbf T\big) \\ &=\frac1{2\pi}\int_0^\infty\mathrm du\operatorname{Tr}_{i,\mathbb R^3}\ln\Big( \frac {\boldsymbol\omega^2 +\boldsymbol\alpha(0)^{\frac12}\boldsymbol\omega \mathbf T\boldsymbol\alpha(0)^{\frac12}\boldsymbol\omega +u^2\mathbf I} {\boldsymbol\omega^2+u^2\mathbf I}\Big) \\ &=\frac1{2\pi}\int_0^\infty\mathrm du\operatorname{Tr}_{i,\mathbb R^3}\Big(\ln\big(\mathbf Q(\{\alpha_j(0),\omega_j\},\mathbf T)+u^2\mathbf I\big)-\ln(\boldsymbol\omega^2+u^2\mathbf I)\big) \\ &=\frac1{2\pi}\int_0^\infty\mathrm du\operatorname{Tr}_{i,\mathbb R^3}\big(\ln(\boldsymbol{\tilde\omega}^2+\mathbf Iu^2)-\ln(\boldsymbol\omega^2+\mathbf Iu^2)\big)\\ &=\frac1{2\pi}\sum_{k=1}^{3N}\int_0^\infty\mathrm du\ln\left(\frac{\tilde\omega_k^2+u^2}{\omega_k^2+u^2}\right) \\ &=\sum_{k=1}^{3N}\frac{\tilde\omega_k-\omega_k}2\equiv\sum_{k=1}^{3N}\frac{\tilde\omega_k}2-3\sum_{i=1}^N\frac{\omega_i}2 \end{aligned}