Variational polaron transform
This note follows a detailed derivation of the partial / variational polaron transformation for a Frenkel exciton model with local harmonic baths, and the minimization of the Gibbs–Bogoliubov–Feynman (GBF) free energy to obtain optimal displacement parameters.
1. Partial polaron transform
We start from a general Frenkel exciton Hamiltonian with site-local baths:
\[ \begin{aligned}{} \hat{H} &= \sum_n \left( \varepsilon_n + \sum_\nu \frac{g_{n,\nu}^2}{2\omega_{n,\nu}^2} \right) \lvert n\rangle\langle n\rvert + \sum_n \sum_{m>n} \Delta_{nm} \big( \lvert n\rangle\langle m\rvert + \lvert m\rangle\langle n\rvert \big) \\ &\quad + \sum_{n,\nu} \frac{1}{2} \left( \hat{P}_{n,\nu}^2 + \omega_{n,\nu}^2 \hat{R}_{n,\nu}^2 \right) - \sum_{n,\nu} g_{n,\nu}\,\hat{R}_{n,\nu} \lvert n\rangle\langle n\rvert. \end{aligned} \tag{1} \]
The diagonal part (system + reorganization) plus bath and system–bath couplings are therefore already of the usual displaced-oscillator form. We define a generalized (partial) polaron transformation
\[ \hat{S} = \sum_{k,\mu} i f_{k,\mu} \hat{P}_{k,\mu} \lvert k\rangle\langle k\rvert. \tag{2} \]
The diagonal projectors commute with this generator:
\[ e^{\hat{S}} \lvert n\rangle\langle n\rvert e^{-\hat{S}} = \lvert n\rangle\langle n\rvert. \tag{3} \]
For off-diagonal operators we use the Baker–Campbell–Hausdorff series:
\[ e^{\hat{S}} \lvert n\rangle\langle m\rvert e^{-\hat{S}} = \lvert n\rangle\langle m\rvert + [\hat{S},\lvert n\rangle\langle m\rvert] + \frac{1}{2!} [\hat{S},\lvert n\rangle\langle m\rvert]^2 + \cdots. \tag{4} \]
The first commutator is
\[ \begin{aligned}{} [\hat{S},\lvert n\rangle\langle m\rvert] &= \sum_{k,\mu} i f_{k,\mu} \big[ \hat{P}_{k,\mu}\lvert k\rangle\langle k\rvert, \lvert n\rangle\langle m\rvert \big] = \sum_{k,\mu} i f_{k,\mu} \hat{P}_{k,\mu} \big[ \lvert k\rangle\langle k\rvert, \lvert n\rangle\langle m\rvert \big] \\ &= \sum_{k,\mu} i f_{k,\mu} \hat{P}_{k,\mu} \big( \lvert k\rangle\langle k\rvert n\rangle\langle m\rvert - \lvert n\rangle\langle m\rvert k\rangle\langle k\rvert \big) \\ &= \sum_{k,\mu} i f_{k,\mu} \hat{P}_{k,\mu} \big( \lvert k\rangle \delta_{kn} \langle m\rvert - \lvert n\rangle \delta_{mk} \langle k\rvert \big) \\ &= \sum_\mu i f_{n,\mu} \hat{P}_{n,\mu} \lvert n\rangle\langle m\rvert - \sum_\mu i f_{m,\mu} \hat{P}_{m,\mu} \lvert n\rangle\langle m\rvert \\ &= i \left( \sum_\mu f_{n,\mu} \hat{P}_{n,\mu} - \sum_\mu f_{m,\mu} \hat{P}_{m,\mu} \right) \lvert n\rangle\langle m\rvert. \end{aligned} \tag{5} \]
Define
\[ \hat{X}_n = \sum_\mu f_{n,\mu}\hat{P}_{n,\mu}, \qquad \hat{X}_{nm} = \hat{X}_n - \hat{X}_m, \]
then
\[ [\hat{S},\lvert n\rangle\langle m\rvert] = i \hat{X}_{nm}\,\lvert n\rangle\langle m\rvert. \]
The second nested commutator becomes
\[ \begin{aligned}{} [\hat{S},\lvert n\rangle\langle m\rvert]^2 &= \big[ \hat{S}, i\hat{X}_{nm}\lvert n\rangle\langle m\rvert \big] = i \sum_{k,\mu} f_{k,\mu} \big[ \hat{P}_{k,\mu}\lvert k\rangle\langle k\rvert, \hat{X}_{nm}\lvert n\rangle\langle m\rvert \big] \\ &= i \sum_{k,\mu} f_{k,\mu} \hat{X}_{nm}\hat{P}_{k,\mu} \big[ \lvert k\rangle\langle k\rvert, \lvert n\rangle\langle m\rvert \big] \\ &= i \sum_{k,\mu} f_{k,\mu}\hat{X}_{nm}\hat{P}_{k,\mu} \big( \lvert k\rangle\delta_{kn}\langle m\rvert - \lvert n\rangle\delta_{mk}\langle k\rvert \big) \\ &= i \sum_\mu f_{n,\mu}\hat{X}_{nm}\hat{P}_{n,\mu}\lvert n\rangle\langle m\rvert + i \sum_\mu f_{m,\mu}\hat{X}_{nm}\hat{P}_{m,\mu}\lvert n\rangle\langle m\rvert \\ &= i \hat{X}_{nm}^2 \lvert n\rangle\langle m\rvert. \end{aligned} \tag{6} \]
Continuing this pattern, the series resums to
\[ \begin{aligned}{} e^{\hat{S}} \lvert n\rangle\langle m\rvert e^{-\hat{S}} &= \sum_{a=0}^\infty \frac{1}{a!} (i\hat{X}_{nm})^a \lvert n\rangle\langle m\rvert = e^{i\hat{X}_{nm}} \lvert n\rangle\langle m\rvert. \end{aligned} \tag{7} \]
2. Bath operator transforms
Momentum operators commute with each other and with the projectors, so they are invariant:
\[ e^{\hat{S}} \hat{P}_{n,\nu} e^{-\hat{S}} = \hat{P}_{n,\nu}. \tag{8} \]
For the bath position operator,
\[ e^{\hat{S}} \hat{R}_{n,\nu} e^{-\hat{S}} = \hat{R}_{n,\nu} + [\hat{S},\hat{R}_{n,\nu}] + \frac{1}{2!} [\hat{S},\hat{R}_{n,\nu}]^2 + \cdots. \tag{9} \]
The first commutator is
\[ \begin{aligned}{} [\hat{S},\hat{R}_{n,\nu}] &= \sum_{k,\mu} i f_{k,\mu} \big[ \hat{P}_{k,\mu}\lvert k\rangle\langle k\rvert, \hat{R}_{n,\nu} \big] = \sum_{k,\mu} i f_{k,\mu} [\hat{P}_{k,\mu},\hat{R}_{n,\nu}] \lvert k\rangle\langle k\rvert \\ &= \sum_{k,\mu} i f_{k,\mu}(-i) \delta_{kn}\delta_{\mu\nu} \lvert n\rangle\langle n\rvert = f_{n,\nu} \lvert n\rangle\langle n\rvert. \end{aligned} \tag{10} \]
Since this contains no further bath operators and commutes with $\hat{S}$, higher commutators vanish and
\[ e^{\hat{S}} \hat{R}_{n,\nu} e^{-\hat{S}} = \hat{R}_{n,\nu} + f_{n,\nu} \lvert n\rangle\langle n\rvert. \tag{11} \]
Transforming the bath Hamiltonian:
\[ \begin{aligned}{} e^{\hat{S}}\hat{H}_b e^{-\hat{S}} &= \frac{1}{2} \sum_{n,\nu} e^{\hat{S}} \big( \hat{P}_{n,\nu}^2 + \omega_{n,\nu}^2 \hat{R}_{n,\nu}^2 \big) e^{-\hat{S}} \\ &= \frac{1}{2} \sum_{n,\nu} \left[ \hat{P}_{n,\nu}^2 + \omega_{n,\nu}^2 \big( \hat{R}_{n,\nu} + f_{n,\nu}\lvert n\rangle\langle n\rvert \big)^2 \right] \\ &= \frac{1}{2} \sum_{n,\nu} \left[ \hat{P}_{n,\nu}^2 + \omega_{n,\nu}^2 \big( \hat{R}_{n,\nu}^2 + 2 f_{n,\nu}\hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert + f_{n,\nu}^2 \lvert n\rangle\langle n\rvert \big) \right] \\ &= \frac{1}{2} \sum_{n,\nu} \big( \hat{P}_{n,\nu}^2 + \omega_{n,\nu}^2\hat{R}_{n,\nu}^2 \big) + \sum_{n,\nu} f_{n,\nu}\omega_{n,\nu}^2 \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert \\ &\qquad + \frac{1}{2} \sum_{n,\nu} \omega_{n,\nu}^2 f_{n,\nu}^2 \lvert n\rangle\langle n\rvert \\ &= \hat{H}_b + \sum_{n,\nu} f_{n,\nu}\omega_{n,\nu}^2 \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert + \frac{1}{2} \sum_{n,\nu} \omega_{n,\nu}^2 f_{n,\nu}^2 \lvert n\rangle\langle n\rvert. \end{aligned} \tag{12} \]
The system–bath Hamiltonian transforms as
\[ \begin{aligned}{} e^{\hat{S}}\hat{H}_{sb}e^{-\hat{S}} &= -\sum_{n,\nu} g_{n,\nu} \big( e^{\hat{S}} \hat{R}_{n,\nu} e^{-\hat{S}} \big) \lvert n\rangle\langle n\rvert \\ &= -\sum_{n,\nu} g_{n,\nu} \big( \hat{R}_{n,\nu} + f_{n,\nu}\lvert n\rangle\langle n\rvert \big) \lvert n\rangle\langle n\rvert \\ &= -\sum_{n,\nu} g_{n,\nu} \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert - \sum_{n,\nu} g_{n,\nu}f_{n,\nu} \lvert n\rangle\langle n\rvert. \end{aligned} \tag{13} \]
Combining, we obtain
\[ \begin{aligned}{} e^{\hat{S}} (\hat{H}_b + \hat{H}_{sb}) e^{-\hat{S}} &= \hat{H}_b + \sum_{n,\nu} (f_{n,\nu}\omega_{n,\nu}^2 - g_{n,\nu}) \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert \\ &\qquad + \sum_{n,\nu} \left( \frac{1}{2} \omega_{n,\nu}^2 f_{n,\nu}^2 - g_{n,\nu}f_{n,\nu} \right) \lvert n\rangle\langle n\rvert. \end{aligned} \tag{14} \]
Completing the square in the last term:
\[ \begin{aligned}{} \frac{1}{2}\omega_{n,\nu}^2 f_{n,\nu}^2 - g_{n,\nu}f_{n,\nu} &= \frac{\omega_{n,\nu}^2}{2} \left( f_{n,\nu}^2 -2\frac{g_{n,\nu}}{\omega_{n,\nu}^2} f_{n,\nu} + \frac{g_{n,\nu}^2}{\omega_{n,\nu}^4} \right) - \frac{g_{n,\nu}^2}{2\omega_{n,\nu}^2} \\ &= \frac{\omega_{n,\nu}^2}{2} \left( f_{n,\nu} - \frac{g_{n,\nu}}{\omega_{n,\nu}^2} \right)^2 - \frac{g_{n,\nu}^2}{2\omega_{n,\nu}^2}. \end{aligned} \]
Hence
\[ \begin{aligned}{} e^{\hat{S}} (\hat{H}_b + \hat{H}_{sb}) e^{-\hat{S}} &= \hat{H}_b + \sum_{n,\nu} (f_{n,\nu}\omega_{n,\nu}^2 - g_{n,\nu}) \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert \\ &\quad + \sum_{n,\nu} \frac{\omega_{n,\nu}^2}{2} \left( f_{n,\nu} - \frac{g_{n,\nu}}{\omega_{n,\nu}^2} \right)^2 \lvert n\rangle\langle n\rvert \\ &\quad - \sum_{n,\nu} \frac{g_{n,\nu}^2}{2\omega_{n,\nu}^2} \lvert n\rangle\langle n\rvert. \end{aligned} \tag{15} \]
Adding this to the diagonal system part cancels the original reorganization energy and yields the transformed diagonal Hamiltonian:
\[ \begin{aligned}{} e^{\hat{S}}\hat{H}_d e^{-\hat{S}} &= \sum_n \left( \varepsilon_n + \sum_\nu \frac{g_{n,\nu}^2}{2\omega_{n,\nu}^2} \right) \lvert n\rangle\langle n\rvert \\ &\quad + \sum_{n,\nu} \left[ \frac{\omega_{n,\nu}^2}{2} \left( f_{n,\nu} - \frac{g_{n,\nu}}{\omega_{n,\nu}^2} \right)^2 - \frac{g_{n,\nu}^2}{2\omega_{n,\nu}^2} \right] \lvert n\rangle\langle n\rvert \\ &\quad + \hat{H}_b + \sum_{n,\nu} (f_{n,\nu}\omega_{n,\nu}^2 - g_{n,\nu}) \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert \\ &= \sum_n \left[ \varepsilon_n + \sum_\nu \frac{\omega_{n,\nu}^2}{2} \left( f_{n,\nu} - \frac{g_{n,\nu}}{\omega_{n,\nu}^2} \right)^2 \right] \lvert n\rangle\langle n\rvert \\ &\quad + \sum_{n,\nu} (f_{n,\nu}\omega_{n,\nu}^2 - g_{n,\nu}) \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert + \hat{H}_b. \end{aligned} \tag{16} \]
The off-diagonal part transforms as
\[ \begin{aligned}{} e^{\hat{S}}\hat{H}_\mathrm{off}e^{-\hat{S}} &= \sum_n\sum_{m>n} \Delta_{n,m} \left( e^{i\hat{X}_{nm}}\lvert n\rangle\langle m\rvert + e^{-i\hat{X}_{nm}}\lvert m\rangle\langle n\rvert \right). \end{aligned} \tag{17} \]
Collecting, the partially polaron-transformed Frenkel exciton Hamiltonian is
\[ \begin{aligned}{} \hat{H}_\mathrm{PT} &= e^{\hat{S}}\hat{H}e^{-\hat{S}} \\ &= \sum_n \left[ \varepsilon_n + \sum_\nu \frac{\omega_{n,\nu}^2}{2} \left( f_{n,\nu} - \frac{g_{n,\nu}}{\omega_{n,\nu}^2} \right)^2 \right] \lvert n\rangle\langle n\rvert \\ &\quad + \sum_{n,\nu} (f_{n,\nu}\omega_{n,\nu}^2 - g_{n,\nu}) \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert \\ &\quad + \sum_n\sum_{m>n} \Delta_{n,m} \left( e^{i\hat{X}_{nm}}\lvert n\rangle\langle m\rvert + e^{-i\hat{X}_{nm}}\lvert m\rangle\langle n\rvert \right) \\ &\quad + \sum_{n,\nu} \frac{1}{2} \left( \hat{P}_{n,\nu}^2 + \omega_{n,\nu}^2 \hat{R}_{n,\nu}^2 \right). \end{aligned} \tag{18} \]
3. Gibbs–Bogoliubov–Feynman free energy
We split the transformed Hamiltonian into non-interacting and interacting parts, $\hat{H}_\mathrm{PT} = \hat{H}_0 + \hat{H}_I$, choosing
\[ \begin{aligned}{} \hat{H}_0 &= \sum_n \left[ \varepsilon_n + \sum_\nu \frac{\omega_{n,\nu}^2}{2} \left( f_{n,\nu} - \frac{g_{n,\nu}}{\omega_{n,\nu}^2} \right)^2 \right] \lvert n\rangle\langle n\rvert \\ &\quad + \sum_n\sum_{m>n} \kappa_{n,m}\Delta_{n,m} \big( \lvert n\rangle\langle m\rvert + \lvert m\rangle\langle n\rvert \big) + \hat{H}_b, \end{aligned} \tag{19} \]
and
\[ \begin{aligned}{} \hat{H}_I &= \sum_{n,\nu} (f_{n,\nu}\omega_{n,\nu}^2 - g_{n,\nu}) \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert \\ &\quad + \sum_n\sum_{m>n} \Delta_{n,m} \left( e^{i\hat{X}_{nm}}\lvert n\rangle\langle m\rvert + e^{-i\hat{X}_{nm}}\lvert m\rangle\langle n\rvert \right) \\ &\quad - \sum_n\sum_{m>n} \kappa_{n,m}\Delta_{n,m} \big( \lvert n\rangle\langle m\rvert + \lvert m\rangle\langle n\rvert \big), \end{aligned} \tag{20} \]
where the thermal-averaged renormalization factor is
\[ \kappa_{n,m} = \mathrm{Tr} \big[ e^{i\hat{X}_{nm}} e^{-\beta\hat{H}_b} \big] = \langle e^{\pm i\hat{X}_{nm}}\rangle_{\hat{H}_b}. \tag{21} \]
By construction,
\[ \langle \hat{H}_I \rangle_{\hat{H}_b} = \mathrm{Tr} \big[ \hat{H}_I e^{-\beta\hat{H}_b} \big] = 0. \tag{22} \]
The GBF upper bound on the free energy is
\[ A_B = -\frac{1}{\beta} \ln \mathrm{Tr} \big[ e^{-\beta\hat{H}_0} \big] + \langle\hat{H}_I\rangle_{\hat{H}_b}. \tag{23} \]
Using Eq. (22),
\[ A_B = -\frac{1}{\beta} \ln \mathrm{Tr} \big[ e^{-\beta\hat{H}_0} \big]. \tag{24} \]
For a two-site (Frenkel biexciton) model, where $\hat{H}_b$ is separable, $\hat{H}_b = \hat{H}_{1,b}\otimes \mathbf{1}_2 + \mathbf{1}_1\otimes\hat{H}_{2,b}$, one finds
\[ \begin{aligned}{} \kappa_{1,2} &= \mathrm{Tr} \big[ e^{i(\hat{X}_1\otimes\mathbf{1}_2 - \mathbf{1}_1\otimes\hat{X}_2)} e^{-\beta(\hat{H}_{1,b}\otimes\mathbf{1}_2 + \mathbf{1}_1\otimes\hat{H}_{2,b})} \big] \\ &= \mathrm{Tr} \big[ e^{i\hat{X}_1 - \beta\hat{H}_{1,b}} \big] \mathrm{Tr} \big[ e^{-i\hat{X}_2 - \beta\hat{H}_{2,b}} \big] \\ &= \kappa_1\kappa_2, \end{aligned} \tag{25–29} \]
where $\kappa_n = \langle e^{\pm i\hat{X}_n}\rangle_{\hat{H}_{n,b}}$ is the single-site renormalization.
4. Renormalization of inter-site coupling
For site 1,
\[ e^{i\hat{X}_1} = \exp\!\left( i\sum_\nu f_{1,\nu}\hat{P}_{1,\nu} \right) = \prod_\nu \exp\big(i f_{1,\nu}\hat{P}_{1,\nu}\big) = \prod_\nu \hat{D}_\nu(i f_{1,\nu}), \tag{30} \]
where $\hat{D}_\nu$ is the displacement operator for mode $\nu$. To evaluate $\kappa_n = \langle e^{i\hat{X}_n}\rangle$, we first compute the position-space matrix elements of a thermal state for a single harmonic mode:
\[ \begin{aligned}{} \langle x \vert \hat{\rho}_\mathrm{th} \vert x' \rangle &= \frac{1}{Z} \langle x \vert e^{-\beta \hat{H}_B} \vert x'\rangle = \frac{1}{Z} \sum_{n,m} \langle x\vert n\rangle \langle n\vert e^{-\beta\hat{H}_B}\vert m\rangle \langle m\vert x'\rangle \\ &= \frac{1}{Z} \sum_{n} e^{-\beta\omega(n+\tfrac{1}{2})} \psi_n(x)\psi_n^\ast(x'). \end{aligned} \tag{31} \]
Using the harmonic oscillator wavefunctions
\[ \psi_n(x) = \left(\frac{\omega}{\pi}\right)^{1/4} \left(\frac{1}{2^n n!}\right)^{1/2} e^{-\omega x^2/2} H_n\big(\sqrt{\omega}\,x\big), \tag{32} \]
we obtain
\[ \begin{aligned}{} \langle x \vert \hat{\rho}_\mathrm{th} \vert x' \rangle &= \frac{e^{-\beta\omega/2}}{Z} \left(\frac{\omega}{\pi}\right)^{1/2} e^{-\tfrac{\omega}{2}(x^2 + x'^2)} \sum_{n=0}^\infty \left( e^{-\beta\omega} \right)^n \frac{1}{2^n n!} H_n(\sqrt{\omega}\,x) H_n(\sqrt{\omega}\,x'). \end{aligned} \tag{33} \]
Using Mehler's identity
\[ \sum_{n=0}^\infty \frac{(\rho/2)^{2n}}{n!} H_n(a)H_n(b) = \frac{1}{\sqrt{1-\rho^2}} \exp\left[ -\frac{\rho^2}{1-\rho^2}(a^2+b^2) + \frac{2\rho}{1-\rho^2}ab \right], \tag{34} \]
and identifying $\rho = e^{-\beta\omega}$, $a = \sqrt{\omega}x$, $b=\sqrt{\omega}x'$, the density matrix becomes
\[ \langle x \vert \hat{\rho}_\mathrm{th} \vert x' \rangle = \left( \frac{\omega}{\pi} \tanh \frac{\beta\omega}{2} \right)^{1/2} \exp\left[ -\frac{\omega}{4} \left( \tanh \frac{\beta\omega}{2} (x+x')^2 + \coth \frac{\beta\omega}{2} (x-x')^2 \right) \right]. \tag{38} \]
The thermal average of a momentum displacement is
\[ \begin{aligned}{} \langle e^{i f \hat{P}}\rangle &= \mathrm{Tr} \big[ e^{i f\hat{P}}\hat{\rho}_\mathrm{th} \big] = \int dx\, \langle x\vert e^{i f\hat{P}} \hat{\rho}_\mathrm{th}\vert x\rangle = \int dx\, \langle x+f \vert \hat{\rho}_\mathrm{th}\vert x\rangle. \end{aligned} \tag{39} \]
Using the form in Eq. (38) and performing the Gaussian integral (Eqs. (40)–(41)), we arrive at
\[ \langle e^{i f \hat{P}}\rangle = \exp\left[ -\frac{1}{4} f^2 \omega \coth\left(\frac{\beta\omega}{2}\right) \right]. \tag{42} \]
For a given site $n$ with modes $\{(n,\nu)\}$,
\[ \begin{aligned}{} \kappa_n &= \left\langle e^{i\hat{X}_n}\right\rangle = \left\langle \prod_\nu e^{i f_{n,\nu}\hat{P}_{n,\nu}} \right\rangle = \prod_\nu \left\langle e^{i f_{n,\nu}\hat{P}_{n,\nu}}\right\rangle \\ &= \prod_\nu \exp\left[ -\frac{1}{4} f_{n,\nu}^2 \omega_{n,\nu} \coth\left(\frac{\beta\omega_{n,\nu}}{2}\right) \right] \\ &= \exp\left[ -\frac{1}{4} \sum_\nu f_{n,\nu}^2 \omega_{n,\nu} \coth\left(\frac{\beta\omega_{n,\nu}}{2}\right) \right]. \end{aligned} \tag{43} \]
For the two-site case,
\[ \kappa_{nm} = \kappa_n \kappa_m = \exp\left[ -\frac{1}{4} \sum_{a\in\{n,m\}} \sum_\nu f_{a,\nu}^2 \omega_{a,\nu} \coth\left(\frac{\beta\omega_{a,\nu}}{2}\right) \right]. \tag{44} \]
5. Minimization of the GBF free energy
For the non-interacting Hamiltonian $\hat{H}_0$ (system + bath) we have
\[ \begin{aligned}{} A_B &= -\frac{1}{\beta} \ln \mathrm{Tr} \big[ e^{-\beta(\hat{H}_S + \hat{H}_B)} \big] = -\frac{1}{\beta} \ln \Big( \mathrm{Tr}[e^{-\beta\hat{H}_S}] \, \mathrm{Tr}[e^{-\beta\hat{H}_B}] \Big) \\ &= -\frac{1}{\beta} \ln \mathrm{Tr}[e^{-\beta\hat{H}_S}] - \frac{1}{\beta}\ln Z, \end{aligned} \tag{45} \]
where $Z = \mathrm{Tr}[e^{-\beta\hat{H}_B}]$ is the bath partition function. For the Frenkel biexciton, we diagonalize the effective system Hamiltonian
\[ \hat{H}_S = \begin{pmatrix} \varepsilon_1 + \delta_1 & \Delta_{12}\kappa_{12} \\ \Delta_{12}\kappa_{12} & \varepsilon_2 + \delta_2 \end{pmatrix} \equiv \begin{pmatrix} \tilde{\varepsilon}_1 & \Delta_\kappa \\ \Delta_\kappa & \tilde{\varepsilon}_2 \end{pmatrix}, \tag{47} \]
with
\[ \Delta_\kappa = \Delta_{12}\kappa_{12}, \qquad \tilde{\varepsilon}_m = \varepsilon_m + \delta_m, \]
and
\[ \delta_n = \sum_\nu \frac{\omega_{n,\nu}^2}{2} \left( f_{n,\nu} - \frac{g_{n,\nu}}{\omega_{n,\nu}^2} \right)^2. \tag{48} \]
The eigenvalues of $\hat{H}_S$ are obtained from
\[ \det(\hat{H}_S - \lambda\mathbf{1}) = 0 \;\Rightarrow\; (\tilde{\varepsilon}_1 - \lambda) (\tilde{\varepsilon}_2 - \lambda) - \Delta_\kappa^2 = 0. \tag{49} \]
The roots are
\[ \lambda_\pm = \frac{\tilde{\varepsilon}_1 + \tilde{\varepsilon}_2}{2} \pm \frac{1}{2} \sqrt{ (\tilde{\varepsilon}_1 - \tilde{\varepsilon}_2)^2 + 4\Delta_\kappa^2 } \equiv \bar{\varepsilon} \pm \frac{\eta}{2}, \tag{50} \]
where
\[ \bar{\varepsilon} = \frac{\tilde{\varepsilon}_1 + \tilde{\varepsilon}_2}{2}, \qquad \eta = \sqrt{ (\tilde{\varepsilon}_1 - \tilde{\varepsilon}_2)^2 + 4\Delta_\kappa^2 }. \]
Hence
\[ \begin{aligned}{} \mathrm{Tr}[e^{-\beta\hat{H}_S}] &= e^{-\beta\lambda_+} + e^{-\beta\lambda_-} = e^{-\beta\bar{\varepsilon}} \big( e^{-\beta\eta/2} + e^{\beta\eta/2} \big) \\ &= 2 e^{-\beta\bar{\varepsilon}} \cosh\left(\frac{\beta\eta}{2}\right). \end{aligned} \tag{51} \]
The GBF free energy becomes
\[ \begin{aligned}{} A_B &= -\frac{1}{\beta} \ln\left[ 2 e^{-\beta\bar{\varepsilon}} \cosh\left(\frac{\beta\eta}{2}\right) \right] - \frac{1}{\beta}\ln Z \\ &= \bar{\varepsilon} - \frac{1}{\beta} \ln\left[ \cosh\left(\frac{\beta\eta}{2}\right) \right] - \frac{1}{\beta}\ln(2Z). \end{aligned} \tag{52} \]
It is useful to introduce the dimensionless parameters
\[ \phi_{n,\nu} = \frac{f_{n,\nu}}{g_{n,\nu}/\omega_{n,\nu}^2}. \tag{53} \]
Then the onsite energy shifts are
\[ \delta_n = \sum_\nu \frac{g_{n,\nu}^2}{2\omega_{n,\nu}^2} (\phi_{n,\nu} - 1)^2. \tag{54} \]
Assuming both excitons have the same spectral density, Eq. (44) becomes
\[ \begin{aligned}{} \kappa_{12} &= \exp\left[ -\frac{1}{4} \sum_{a\in\{1,2\}} \sum_\nu f_{a,\nu}^2 \omega_{a,\nu} \coth\left(\frac{\beta\omega_{a,\nu}}{2}\right) \right] \\ &= \exp\left[ -\sum_\nu \big( \phi_{1,\nu}^2 + \phi_{2,\nu}^2 \big) \left( \frac{g_{1,\nu}}{2\omega_{1,\nu}} \right)^2 \omega_{1,\nu} \coth\left(\frac{\beta\omega_{1,\nu}}{2}\right) \right]. \end{aligned} \tag{55} \]
The diagonal part of the variationally polaron-transformed system Hamiltonian can then be written as
\[ \hat{H}_S^\mathrm{vPT} = \sum_n \left[ \varepsilon_n + \sum_\nu \frac{g_{n,\nu}^2}{2\omega_{n,\nu}^2} (\phi_{n,\nu}-1)^2 \right] \lvert n\rangle\langle n\rvert, \tag{56} \]
the residual system–bath coupling is
\[ \hat{H}_{SB}^\mathrm{vPT} = \sum_{n,\nu} g_{n,\nu}(\phi_{n,\nu}-1) \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert, \tag{57} \]
and the collective momentum operators become
\[ \hat{X}_n = \sum_\nu f_{n,\nu}\hat{P}_{n,\nu} = \sum_\nu \phi_{n,\nu} \left( \frac{g_{n,\nu}}{\omega_{n,\nu}^2} \right) \hat{P}_{n,\nu}, \tag{58} \]
so that
\[ \hat{X}_{12} = \hat{X}_1 - \hat{X}_2 = \sum_\nu \left(\frac{g_{1,\nu}}{\omega_{1,\nu}^2}\right) \big( \phi_{1,\nu}\hat{P}_{1,\nu} - \phi_{2,\nu}\hat{P}_{2,\nu} \big). \tag{59} \]
The full variationally transformed Hamiltonian for the biexciton becomes
\[ \begin{aligned}{} \hat{H}_\mathrm{vPT} &= \sum_n \left[ \varepsilon_n + \sum_\nu \frac{g_{n,\nu}^2}{2\omega_{n,\nu}^2} (\phi_{n,\nu}-1)^2 \right] \lvert n\rangle\langle n\rvert \\ &\quad + \sum_{n,\nu} g_{n,\nu}(\phi_{n,\nu}-1) \hat{R}_{n,\nu}\lvert n\rangle\langle n\rvert \\ &\quad + \Delta_{12} \left( e^{i\hat{X}_{12}}\lvert 1\rangle\langle 2\rvert + e^{-i\hat{X}_{12}}\lvert 2\rangle\langle 1\rvert \right) \\ &\quad + \sum_{n,\nu} \frac{1}{2} \left( \hat{P}_{n,\nu}^2 + \omega_{n,\nu}^2 \hat{R}_{n,\nu}^2 \right). \end{aligned} \tag{60} \]
Instead of directly minimizing Eq. (52), one can set $\partial A_B / \partial f_{n,\nu} = 0$ for each $(n,\nu)$ and solve the resulting coupled equations for the optimal $f_{n,\nu}$.
5.1 Minimization with respect to exciton 1
Taking the derivative of $A_B$ with respect to $f_{1,\nu}$,
\[ \begin{aligned}{} \frac{\partial A_B}{\partial f_{1,\nu}} &= \frac{\partial \bar{\varepsilon}}{\partial f_{1,\nu}} - \frac{1}{\beta} \frac{\partial}{\partial f_{1,\nu}} \ln\left[ \cosh\left(\frac{\beta\eta}{2}\right) \right] \\ &= \frac{1}{2} \frac{\partial \tilde{\varepsilon}_1}{\partial f_{1,\nu}} - \frac{1}{2} \tanh\left(\frac{\beta\eta}{2}\right) \frac{\partial \eta}{\partial f_{1,\nu}}. \end{aligned} \tag{62} \]
The first term is
\[ \begin{aligned}{} \frac{1}{2} \frac{\partial \delta_1}{\partial f_{1,\nu}} &= \frac{1}{2} \frac{\partial}{\partial f_{1,\nu}} \sum_\mu \frac{\omega_{1,\mu}^2}{2} \left( f_{1,\mu} - \frac{g_{1,\mu}}{\omega_{1,\mu}^2} \right)^2 \\ &= \sum_\mu \frac{\omega_{1,\mu}^2}{2} \left( f_{1,\mu} - \frac{g_{1,\mu}}{\omega_{1,\mu}^2} \right) \delta_{\mu,\nu} \\ &= \frac{\omega_{1,\nu}^2}{2} \left( f_{1,\nu} - \frac{g_{1,\nu}}{\omega_{1,\nu}^2} \right). \end{aligned} \tag{63} \]
For $\eta$,
\[ \begin{aligned}{} \frac{\partial \eta}{\partial f_{1,\nu}} &= \frac{\partial}{\partial f_{1,\nu}} \sqrt{ (\tilde{\varepsilon}_1 - \tilde{\varepsilon}_2)^2 + 4\Delta_\kappa^2 } \\ &= \frac{1}{2\eta} \frac{\partial}{\partial f_{1,\nu}} \left[ (\tilde{\varepsilon}_1 - \tilde{\varepsilon}_2)^2 + 4\Delta_\kappa^2 \right] \\ &= \frac{1}{\eta} (\tilde{\varepsilon}_1 - \tilde{\varepsilon}_2) \frac{\partial\tilde{\varepsilon}_1}{\partial f_{1,\nu}} + \frac{4\Delta_\kappa}{\eta} \frac{\partial\Delta_\kappa}{\partial f_{1,\nu}}. \end{aligned} \tag{64} \]
Using $\Delta_\kappa = \Delta_{12}\kappa_{12}$ with $\kappa_{12}=\kappa_1\kappa_2$,
\[ \frac{\partial \kappa_1}{\partial f_{1,\nu}} = -\frac{\kappa_1}{2} f_{1,\nu}\omega_{1,\nu} \coth\left(\frac{\beta\omega_{1,\nu}}{2}\right), \tag{65} \]
so
\[ \begin{aligned}{} \frac{\partial \eta}{\partial f_{1,\nu}} &= \frac{\tilde{\varepsilon}_1 - \tilde{\varepsilon}_2}{\eta} \frac{\partial \delta_1}{\partial f_{1,\nu}} - \frac{2\Delta_\kappa^2}{\eta} f_{1,\nu}\omega_{1,\nu} \coth\left(\frac{\beta\omega_{1,\nu}}{2}\right). \end{aligned} \tag{66} \]
Substituting Eqs. (63) and (66) into Eq. (62),
\[ \begin{aligned}{} \frac{\partial A_B}{\partial f_{1,\nu}} &= \frac{\omega_{1,\nu}^2}{2} \left( f_{1,\nu} - \frac{g_{1,\nu}}{\omega_{1,\nu}^2} \right) \\ &\quad - \frac{1}{2} \tanh\left(\frac{\beta\eta}{2}\right) \left[ \frac{\tilde{\varepsilon}_1 - \tilde{\varepsilon}_2}{\eta} \frac{\partial \delta_1}{\partial f_{1,\nu}} - \frac{2\Delta_\kappa^2}{\eta} f_{1,\nu}\omega_{1,\nu} \coth\left(\frac{\beta\omega_{1,\nu}}{2}\right) \right] \\ &= \frac{\omega_{1,\nu}^2}{2} \left( f_{1,\nu} - \frac{g_{1,\nu}}{\omega_{1,\nu}^2} \right) \\ &\quad - \frac{1}{2} \tanh\left(\frac{\beta\eta}{2}\right) \left[ \frac{\tilde{\varepsilon}_1 - \tilde{\varepsilon}_2}{\eta} \omega_{1,\nu}^2 \left( f_{1,\nu} - \frac{g_{1,\nu}}{\omega_{1,\nu}^2} \right) - \frac{2\Delta_\kappa^2}{\eta} f_{1,\nu}\omega_{1,\nu} \coth\left(\frac{\beta\omega_{1,\nu}}{2}\right) \right]. \end{aligned} \tag{67} \]
Setting $\partial A_B / \partial f_{1,\nu} = 0$ and rearranging gives
\[ \begin{aligned}{} f_{1,\nu} &= \frac{g_{1,\nu}}{\omega_{1,\nu}^2} \, \frac{ 1 - \displaystyle \frac{\Delta\tilde{\varepsilon}}{\eta} \tanh\left(\frac{\beta\eta}{2}\right) }{ 1 - \displaystyle \frac{1}{\eta} \tanh\left(\frac{\beta\eta}{2}\right) \left[ \Delta\tilde{\varepsilon} - 2 \frac{\Delta_\kappa^2}{\omega_{1,\nu}} \coth\left(\frac{\beta\omega_{1,\nu}}{2}\right) \right] }, \end{aligned} \tag{72} \]
where $\Delta\tilde{\varepsilon} = \tilde{\varepsilon}_1 - \tilde{\varepsilon}_2$.
5.2 Minimization with respect to exciton 2
Similarly, for $f_{2,\nu}$ one finds
\[ \frac{\partial A_B}{\partial f_{2,\nu}} = \frac{1}{2} \frac{\partial\delta_2}{\partial f_{2,\nu}} - \frac{1}{2} \tanh\left(\frac{\beta\eta}{2}\right) \frac{\partial\eta}{\partial f_{2,\nu}}, \tag{73} \]
with
\[ \frac{1}{2}\frac{\partial\delta_2}{\partial f_{2,\nu}} = \frac{\omega_{2,\nu}^2}{2} \left( f_{2,\nu} - \frac{g_{2,\nu}}{\omega_{2,\nu}^2} \right), \tag{74} \]
and
\[ \begin{aligned}{} \frac{\partial\eta}{\partial f_{2,\nu}} &= \frac{\tilde{\varepsilon}_2 - \tilde{\varepsilon}_1}{\eta} \frac{\partial\delta_2}{\partial f_{2,\nu}} - \frac{2\Delta_\kappa^2}{\eta} f_{2,\nu}\omega_{2,\nu} \coth\left(\frac{\beta\omega_{2,\nu}}{2}\right). \end{aligned} \tag{77} \]
Substituting these into Eq. (73), setting $\partial A_B / \partial f_{2,\nu} = 0$, and rearranging yields
\[ \begin{aligned}{} f_{2,\nu} &= \frac{g_{2,\nu}}{\omega_{2,\nu}^2} \, \frac{ 1 + \displaystyle \frac{\Delta\tilde{\varepsilon}}{\eta} \tanh\left(\frac{\beta\eta}{2}\right) }{ 1 + \displaystyle \frac{1}{\eta} \tanh\left(\frac{\beta\eta}{2}\right) \left[ \Delta\tilde{\varepsilon} + 2 \frac{\Delta_\kappa^2}{\omega_{2,\nu}} \coth\left(\frac{\beta\omega_{2,\nu}}{2}\right) \right] }. \end{aligned} \tag{81} \]
Eqs. (72) and (81) together define the optimal variational displacements $\{f_{n,\nu}\}$ that minimize the GBF free energy $A_B$ for the variational polaron frame.