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} \]