Hamiltonian Dynamics and Fundamental Phenomena in Biophysics: A Review

1. Introduction

A fundamental challenge in biophysics is to understand how biochemical reactions in living cells proceed with the speed and selectivity observed experimentally, given the crowded, thermally noisy environment of the cytoplasm. Standard models invoke short-range forces, electrostatic, van der Waals, hydrophobic, and diffusion-driven random encounters. These mechanisms are undoubtedly operative, but they do not obviously account for the remarkably high efficiency and rapidity with which cognate molecular partners find one another among thousands of competing species in a sub-micron space [1,2,3].

A minority but persistent view in theoretical biophysics, originating with H. Fröhlich in the late 1960s [4,5,6,7], proposes that biological macromolecules driven out of thermodynamic equilibrium by metabolic energy could undergo a collective ordering transition, a phonon condensation, in which the input energy concentrates into the lowest-frequency vibrational mode of a macromolecule, producing a coherent giant oscillating dipole moment. This coherent oscillation was predicted to activate resonant, selective, long-range electrodynamic forces between macromolecules vibrating at the same frequency [6], forces that could accelerate and direct molecular recognition processes.

Three key obstacles have historically limited the impact of this programme. First, the original quantum Wu–Austin microscopic Hamiltonian [8,9] suffers from an instability that requires stabilisation. Second, the connection between the rate-equation level and direct Hamiltonian dynamics had not been established in a classical framework accessible to modern molecular dynamics. Third, and most critically, no experimental evidence for out-of-equilibrium collective oscillations of macromolecules existed until recently.

All these obstacles have now been overcome. The present work reviews the following pivotal topics: i) A classical Hamiltonian in action-angle variables is derived from the quantum Wu–Austin model via the TDVP applied to product coherent states (Sec. 2–4). Nonlinear rate equations follow from the Koopman–von Neumann (KvN) reformulation of the Liouville equation. Their stationary solutions exhibit classical phonon condensation [10]. ii) classical electrodynamic Hamiltonian for two coupled oscillating dipoles is analysed via normal-mode theory (Sec. 5). The result is a long-range ($\sim 1/r^3$) interaction potential active only out of thermal equilibrium [11]. iii) A position-space classical Hamiltonian dynamics, closer to standard molecular dynamics, is simulated directly (Sec. 6) recovering the outcomes of the Fröhlich rate equations, and showing the formation of strong condensates at room temperature are found under resonant bath-mediated coupling [12]. iv) The TDVP is applied to a Davydov–Holstein–Fröhlich Hamiltonian for electron–phonon motion along a specific DNA sequence and its cognate restriction enzyme EcoRI (Sec. 7), and the associate dynamics yields a sharp cross-spectral co-resonance between the electron currents of the two molecules appears for the canonical recognition sequence and disappears upon randomisation [13].

Experimental evidence is reviewed in Sec. 8. Sec. 9 gives an outlook of future developments.

2. The Wu–Austin Quantum Model and the Classical Limit

2.1. Physical Picture

Fröhlich [4,5] modelled a macromolecule as an open system: an ensemble of coupled vibrational modes in contact simultaneously with a thermal bath at temperature $T_B$ and with an energy source at effective temperature $T_S\gg T_B$. He argued that when the energy-input rate exceeds a threshold, the system undergoes a nonequilibrium phase transition: instead of distributing the supplied energy among all modes (equipartition), it channels nearly all of it into the lowest-frequency mode, producing a macroscopic coherent oscillation.

2.1.1. Fr öhlich Rate Equations and Condensation

At the phenomenological level, Fröhlich described the dynamics in terms of the occupation numbers $n_k\left(t\right)=\langle {\widehat{a}}_{{\omega }_k}^{†}{\widehat{a}}_{{\omega }_k}\rangle$ of the vibrational modes [7]. Consider a set of z polar modes with frequencies ${\omega }_1\le {\omega }_k\le {\omega }_z$, $k=1,\dots ,z$, coupled to a heat bath at inverse temperature ${\beta }_B=1/k_B{T}_B$ and driven by an external source that injects quanta into mode k at rate $s_k$. Two classes of bath-mediated processes are retained.

First-order processes exchange a single quantum between mode k and the bath, tending to restore thermal equilibrium. Second-order processes transfer excitations between two different system modes k and j with the bath absorbing or supplying the energy difference; crucially, they conserve the total number of system quanta, ${\sum }_k{\dot{n}}_k^{\left(2\right)}=0$. Denoting by $\phi$ the linear relaxation coefficient and by $\chi$ the nonlinear intermode coupling, the complete rate equations are (for $k=1,\dots ,z$)

$${\dot{n}}_k=s_k-\phi \left[n_k{e}^{{\beta }_B\hslash {\omega }_k}-(n_k+1)\right]-\chi \sum _{j\ne k}\left[n_k(1+n_j)e^{{\beta }_B\hslash {\omega }_k}-n_j(1+n_k)e^{{\beta }_B\hslash {\omega }_j}\right].$$

(1)

The first term injects energy from the source; the second drives each mode toward thermal equilibrium with the bath; the third redistributes quanta across the spectrum without altering their total number.

Effective chemical potential.

This nonlinear redistribution is the key to the Fröhlich effect. Because downward transitions in frequency are statistically favoured by the Bose factors $(1+n_j)$, pumped energy is not simply equipartitioned. Once the total input rate $S={\sum }_k{s}_k$ exceeds a threshold set by the competition between pumping, dissipation, and nonlinear scattering, the excess population can no longer be absorbed by the higher-frequency modes and accumulates at the bottom of the spectrum. In the stationary regime one introduces an effective nonequilibrium chemical potential $\mu$, defined self-consistently by

$$e^{-{\beta }_B\mu }={\displaystyle \frac{{\displaystyle \phi +\chi \sum _j(1+n_j)}}{{\displaystyle \phi +\chi \sum _j{n}_j{e}^{{\beta }_B\hslash {\omega }_j}}}},$$

(2)

which satisfies $0\le \mu <\hslash {\omega }_1$ whenever $S>0$. The formal stationary solution then takes the Bose-like form

$$n_k=\left(1+{\displaystyle \frac{\phi }{\chi }}{\displaystyle \frac{s_k}{S}}\left(1-e^{-{\beta }_B\mu }\right)\right){\displaystyle \frac{1}{e^{{\beta }_B(\hslash {\omega }_k-\mu )}-1}},$$

(3)

so that as S increases, $\mu$ approaches $\hslash {\omega }_1$ from below and $n_1$ grows without bound while all other occupations remain comparatively limited.

Fröhlich condensation.

This is the Fröhlich condensation phenomenon: a nonequilibrium ordered state in which a single low-frequency polar mode becomes macroscopically populated and behaves as a collective coherent oscillation of the macromolecule. The analogy with Bose–Einstein condensation is instructive: in both cases the discrete sum over modes cannot be replaced by an integral once $\mu$ saturates at the band edge, forcing a macroscopic occupation of the ground mode. The mechanism is, however, intrinsically a nonequilibrium one. In an equilibrium Bose gas the condensed state is reached by lowering the temperature; here it is reached by maintaining a sufficiently strong energy throughput in an open, dissipative system at fixed $T_B$. When the nonlinear coupling vanishes ($\chi =0$), equation (2) gives $\mu =0$ and (3) reduces to a modified Planck distribution with no condensation; it is therefore the nonlinear bath-mediated intermode scattering that is the essential ingredient.

The predicted coherent oscillation lies in the sub-THz range, $\nu \sim {10}^{11}$ Hz, set by acoustic-like breathing vibrations of proteins or membrane sections of linear size $\ell \sim {10}^{-6}$ cm with effective sound velocity $v\sim {10}^5$ cm/s [7]. Once the lowest polar mode is macroscopically occupied it acts as an oscillating giant dipole; cell water with small ions in solution, with typical Debye length of 10Å, screens static charges efficiently, but does not screen oscillating fields at these frequencies. Two such coherently oscillating dipoles at separation R interact via a frequency-selective, attractive potential that in the near-resonance limit falls off as [7]

$$I\simeq -U{\displaystyle \frac{e^2{\left(z_e{z}_s\right)}^{1/2}}{M{\omega }_e^2{R}^3}}\propto -{\displaystyle \frac{1}{R^3}},$$

(4)

where $U=\hslash n_{-}{\omega }_e$ is the excitation energy, $z_e,z_s$ the numbers of coherently oscillating charges, M their mass, and $n_{-}$ the occupation of the lower normal mode of the coupled pair. This long-range force is proposed as the mechanism for frequency-selective recognition between enzymes and their substrates and for the attraction of homologous chromosomes during meiosis [7].

2.1.2. Microscopic Derivation: Wu–Austin Hamiltonian

A microscopic quantum Hamiltonian from which Fröhlich’s rate equations can be derived was given by Wu and Austin [8,9]. The system consists of three sets of quantum harmonic oscillators: normal modes of the macromolecule (frequencies ${\omega }_i\in I_{\mathrm{sys}}$, operators ${\widehat{a}}_{{\omega }_i}$, ${\widehat{a}}_{{\omega }_i}^{†}$); a thermal bath at $T_B$ (frequencies ${\Omega }_j\in I_{\mathrm{bth}}$, operators ${\widehat{b}}_{{\Omega }_j}$, ${\widehat{b}}_{{\Omega }_j}^{†}$); and an energy source modelled as a second bath at $T_S\gg T_B$ (frequencies ${\Omega }_k^{\prime }\in I_{\mathrm{src}}$, operators ${\widehat{c}}_{{\Omega }_k^{\prime }}$, ${\widehat{c}}_{{\Omega }_k^{\prime }}^{†}$). The total Hamiltonian is ${\widehat{H}}_{\mathrm{Tot}}={\widehat{H}}_0+{\widehat{H}}_{\mathrm{int}}$, where

$${\widehat{H}}_0=\sum _{{\omega }_i}\hslash {\omega }_i{\widehat{a}}_{{\omega }_i}^{†}{\widehat{a}}_{{\omega }_i}+\sum _{{\Omega }_j}\hslash {\Omega }_j{\widehat{b}}_{{\Omega }_j}^{†}{\widehat{b}}_{{\Omega }_j}+\sum _{{\Omega }_k^{\prime }}\hslash {\Omega }_k^{\prime }{\widehat{c}}_{{\Omega }_k^{\prime }}^{†}{\widehat{c}}_{{\Omega }_k^{\prime }},$$

(5)

the interaction contains the Wu–Austin trilinear couplings

$${\widehat{H}}_{\mathrm{int}}^{\mathrm{WA}}=\sum _{{\omega }_i,{\Omega }_j}{\eta }_{{\omega }_i{\Omega }_j}{\widehat{a}}_{{\omega }_i}^{†}{\widehat{b}}_{{\Omega }_j}+\sum _{{\omega }_i,{\Omega }_k^{\prime }}{\xi }_{{\omega }_i{\Omega }_k^{\prime }}{\widehat{a}}_{{\omega }_i}^{†}{\widehat{c}}_{{\Omega }_k^{\prime }}+\sum _{{\omega }_i,{\omega }_j,{\Omega }_k}{\chi }_{{\omega }_i{\omega }_j{\Omega }_k}{\widehat{a}}_{{\omega }_i}^{†}{\widehat{a}}_{{\omega }_j}{\widehat{b}}_{{\Omega }_k}^{†}+\mathrm{H}.\mathrm{c}.,$$

(6)

supplemented by a stabilising quartic term [14]

$${\widehat{H}}_{\mathrm{int}}^Q=\sum _{{\omega }_i,{\omega }_j,{\omega }_k,{\omega }_l}\left[{\kappa }_{{\omega }_i{\omega }_j{\omega }_k{\omega }_l}^{\left(1\right)}{\widehat{a}}_{{\omega }_i}^{†}{\widehat{a}}_{{\omega }_j}^{†}{\widehat{a}}_{{\omega }_k}{\widehat{a}}_{{\omega }_l}+{\kappa }_{{\omega }_i{\omega }_j{\omega }_k{\omega }_l}^{\left(2\right)}{\widehat{a}}_{{\omega }_i}^{†}{\widehat{a}}_{{\omega }_j}^{†}{\widehat{a}}_{{\omega }_k}^{†}{\widehat{a}}_{{\omega }_l}+\mathrm{H}.\mathrm{c}.\right].$$

(7)

The quartic term removes the pathology of the pure Wu–Austin model and corresponds physically to anharmonic interactions among the macromolecular modes. Fröhlich’s phenomenological rate equations (1) are recovered from this Hamiltonian by applying the quantum Liouville–von Neumann equation $i\hslash {\partial }_t\widehat{\rho }=[{\widehat{H}}_{\mathrm{Tot}},\widehat{\rho }]$ to the full density matrix $\widehat{\rho }$, tracing over the bath and source degrees of freedom in the Born–Markov approximation, and identifying $n_k=\mathrm{Tr}\left({\widehat{a}}_{{\omega }_k}^{†}{\widehat{a}}_{{\omega }_k}{\widehat{\rho }}_{\mathrm{sys}}\right)$. The trilinear terms in ${\widehat{H}}_{\mathrm{int}}^{\mathrm{WA}}$ generate the linear relaxation ($\phi$) and the pumping ($s_k$), while the trilinear system–bath coupling ${\chi }_{{\omega }_i{\omega }_j{\Omega }_k}$ produces the nonlinear intermode redistribution term ($\chi$) in (1); the quartic term ${\widehat{H}}_{\mathrm{int}}^Q$ ensures the thermodynamic stability of the condensed state.

2.1.3. Why a Classical Description Suffices

At room temperature ($T\approx 300$ K) and at the sub-THz frequencies relevant to collective protein vibrations ($\omega /2\pi \sim 0.1$–1 THz), $k_{B}T/\hslash \omega \gg 1$, so the field is in the classical regime and the Bose–Einstein mean occupation number ${\langle n\rangle }_{\mathrm{BE}}$ is well approximated by the Boltzmann estimate $k_{B}T/\hslash \omega$, ranging from $\approx 6$ at 1 THz to $\approx 62$ at $0.1$ THz. When driven out of equilibrium the effective temperature of each mode increases further, making the classical approximation even better. A classical framework is therefore not only legitimate but natural, and it connects directly to the language of molecular dynamics simulation.

3. Dequantization via the Time-Dependent Variational Principle

The time-dependent variational principle (TDVP) [15,16,17] extracts a classical Hamiltonian from a quantum one by restricting dynamics to a finite-dimensional manifold of trial states $\left|\psi \right(x_1,\dots ,x_N)\rangle$ parametrised by real variables $\left\{x_i\left(t\right)\right\}$. The equations of motion follow from

$$\delta S=0,S={\int }_0^{t}L(\psi ,\overline{\psi })d{t}^{\prime },$$

(8)

with the Lagrangian

$$L=i\hslash {\displaystyle \frac{\langle \psi |\dot{\psi }\rangle -\langle \dot{\psi }|\psi \rangle }{2\langle \psi |\psi \rangle }}-{\displaystyle \frac{\langle \psi \left|{\widehat{H}}_{\mathrm{Tot}}\right|\psi \rangle }{\langle \psi |\psi \rangle }}.$$

(9)

The resulting classical Hamiltonian is simply the expectation value $H_{\mathrm{Tot}}=\langle \psi |{\widehat{H}}_{\mathrm{Tot}}|\psi \rangle$. This procedure has been advocated as a “dequantization” method, a kind of inverse of geometric quantization [17]. Since ${\widehat{H}}_{\mathrm{Tot}}$ is built from bosonic creation and annihilation operators, the natural trial states are products of coherent states over all modes in $I_S=I_{\mathrm{sys}}\cup I_{\mathrm{bth}}\cup I_{\mathrm{src}}$:

$$|\Psi \left(t\right)\rangle =\underset{{\omega }_i\in I_{\mathrm{sys}}}{⨂}|z_{{\omega }_i}\left(t\right)\rangle \otimes \underset{{\Omega }_j\in I_{\mathrm{bth}}}{⨂}|z_{{\Omega }_j}\left(t\right)\rangle \otimes \underset{{\Omega }_k^{\prime }\in I_{\mathrm{src}}}{⨂}|z_{{\Omega }_k^{\prime }}\left(t\right)\rangle ,$$

(10)

where $|z\rangle =e^{-{\left|z\right|}^2/2}{\sum }_{k=0}^{\infty }{z}^k/\sqrt{k!}|k\rangle$ satisfies $\langle \widehat{n}\rangle ={\left|z\right|}^2$. Writing $z_i=n_i^{1/2}{e}^{-i{\theta }_i}$ introduces real pairs $(n_i,{\theta }_i)$ for every mode. A direct computation of the symplectic tensor $W_{ij}={\partial }_i{w}_j-{\partial }_j{w}_i$, with $w_i=i\hslash \langle \psi |{\partial }_{x_i}\psi \rangle$, yields

$$W_{n_i{\theta }_k}=-W_{{\theta }_k{n}_i}=\hslash {\delta }_{ik},W_{{\theta }_i{\theta }_k}=W_{n_i{n}_k}=0.$$

(11)

Hence $J_{{\omega }_i}\equiv \hslash n_{{\omega }_i}$ and ${\theta }_{{\omega }_i}$ are canonically conjugated variables, $\{J_{{\omega }_i},{\theta }_{{\omega }_k}\}={\delta }_{ik}$. The action $J_{\omega }=\hslash \langle {\widehat{n}}_{\omega }\rangle$ equals ℏ times the mean occupation number.

3.1. The Classical Hamiltonian in Action-Angle Variables

Evaluating $H_{\mathrm{Tot}}=\langle \Psi |{\widehat{H}}_{\mathrm{Tot}}|\Psi \rangle$ in $(J_{\omega },{\theta }_{\omega })$ variables one finds [10]:

$$\begin{array}{cc}\hfill H_{\mathrm{Tot}}& =\sum _{{\omega }_i}{\omega }_i{J}_{{\omega }_i}+\sum _{{\Omega }_j}{\Omega }_j{J}_{{\Omega }_j}+\sum _{{\Omega }_k^{\prime }}{\Omega }_k^{\prime }{J}_{{\Omega }_k^{\prime }}+\sum _{{\omega }_i,{\Omega }_j}{\eta }_{{\omega }_i{\Omega }_j}{J}_{{\omega }_i}^{1/2}{J}_{{\Omega }_j}^{1/2}cos({\theta }_{{\omega }_i}-{\theta }_{{\Omega }_j})\hfill \\ \hfill & +\sum _{{\omega }_i,{\Omega }_k^{\prime }}{\xi }_{{\omega }_i{\Omega }_k^{\prime }}{J}_{{\omega }_i}^{1/2}{J}_{{\Omega }_k^{\prime }}^{1/2}cos({\theta }_{{\omega }_i}-{\theta }_{{\Omega }_k^{\prime }})+\sum _{{\omega }_i,{\omega }_j,{\Omega }_k}{\chi }_{{\omega }_i{\omega }_j{\Omega }_k}{J}_{{\omega }_i}^{1/2}{J}_{{\omega }_j}^{1/2}{J}_{{\Omega }_k}^{1/2}cos({\theta }_{{\omega }_i}-{\theta }_{{\omega }_j}+{\theta }_{{\Omega }_k})\hfill \\ \hfill & +\sum _{{\omega }_i,{\omega }_j,{\omega }_k,{\omega }_l}{J}_{{\omega }_i}^{1/2}{J}_{{\omega }_j}^{1/2}{J}_{{\omega }_k}^{1/2}{J}_{{\omega }_l}^{1/2}\times \left[{\kappa }_{\dots }^{\left(1\right)}cos({\theta }_{{\omega }_i}+{\theta }_{{\omega }_j}-{\theta }_{{\omega }_k}-{\theta }_{{\omega }_l})+{\kappa }_{\dots }^{\left(2\right)}cos({\theta }_{{\omega }_i}+{\theta }_{{\omega }_j}+{\theta }_{{\omega }_k}-{\theta }_{{\omega }_l})\right].\hfill \end{array}$$

(12)

This is a fully classical, canonical Hamiltonian retaining the complete structure of the quantum model. Hamilton’s equations ${\dot{J}}_{\omega }=-\partial H/\partial {\theta }_{\omega }$, ${\dot{\theta }}_{\omega }=\partial H/\partial J_{\omega }$ govern all degrees of freedom.

3.2. Classical Rate Equations via the Koopman–Von Neumann Formalism

Rather than integrating Hamilton’s equations with an explicit bath, one works with the phase-space probability density $\rho (\left\{(J_{\omega },{\theta }_{\omega })\right\};t)$ satisfying the Liouville equation

$${\displaystyle \frac{\partial \rho }{\partial t}}=\{H_{\mathrm{Tot}},\rho \}\equiv -i{\widehat{\mathcal{L}}}_H\rho .$$

(13)

The Koopman–von Neumann (KvN) formalism [18,19] rewrites this as a Schrödinger-like equation $i{\partial }_t\psi ={\widehat{\mathcal{L}}}_H\psi$, $\rho ={\left|\psi \right|}^2$, on the Hilbert space $L^2\left(\Lambda \right)$, allowing interaction-picture perturbation theory to be applied directly to the classical problem. In this framework the action variables $J_{{\omega }_i}=\hslash n_{{\omega }_i}$ play the role of classical counterparts of the quantum occupation numbers, and the conjugate angles ${\theta }_{{\omega }_i}$ are the corresponding phases. Decomposing ${\widehat{\mathcal{L}}}_H={\widehat{\mathcal{L}}}_{H_0}+{\widehat{\mathcal{L}}}_{H_{\mathrm{int}}}$ and treating ${\widehat{\mathcal{L}}}_{H_{\mathrm{int}}}$ as a time-dependent perturbation adiabatically switched on from $t_0=0$, one works in the interaction picture and expands to second order. Averaging over thermal-bath initial conditions compatible with temperatures $T_B$ (bath) and $T_S\gg T_B$ (external source) then yields closed rate equations for $\langle J_{{\omega }_i}\rangle$, in full analogy with the quantum derivation of Fröhlich’s equations via time-dependent perturbation theory [8,9].

3.3. The Classical Fröhlich-like Rate Equations

The result of the second-order KvN expansion [10] is

$${\displaystyle \frac{d\langle J_{{\omega }_i}\rangle }{dt}}=s_{{\omega }_i}+b_{{\omega }_i}\left({\displaystyle \frac{k_B{T}_B}{{\omega }_i}}-\langle J_{{\omega }_i}\rangle \right)+\sum _{{\omega }_j\ne {\omega }_i}{c}_{{\omega }_i{\omega }_j}\left[\langle J_{{\omega }_j}\rangle -\langle J_{{\omega }_i}\rangle +{\displaystyle \frac{{\omega }_j-{\omega }_i}{k_B{T}_B}}\langle J_{{\omega }_i}\rangle \langle J_{{\omega }_j}\rangle \right]+\left(\mathrm{quartic}\right),$$

(14)

where $s_{{\omega }_i}>0$ is the energy-input rate from the external source, $b_{{\omega }_i}$ is the linear system-bath coupling constant (with characteristic thermalization timescale ${\tau }_{\mathrm{therm}}\approx b_{{\omega }_i}^{-1}$), and $c_{{\omega }_i{\omega }_j}$ encodes the nonlinear mode-mode coupling mediated by the bath. The quartic correction terms, absent in the original Wu–Austin model but necessary to bound the energy from below [14], involve the coupling constants ${\kappa }_{{\omega }_i{\omega }_j{\omega }_k{\omega }_l}^{(1,2,3)}$ of the anharmonic interaction among normal modes [10].

This equation has the following key properties. At equilibrium ($s_{{\omega }_i}=0$) the unique stationary solution is $\langle J_{{\omega }_i}\rangle =k_B{T}_B/{\omega }_i$ (energy equipartition), since $y_{{\omega }_i}={\omega }_i\langle J_{{\omega }_i}\rangle /\left(k_B{T}_B\right)=1$ for all modes. Equation (14) is functionally identical to Fröhlich’s original quantum rate equations, with $\langle J_{{\omega }_i}\rangle$ replacing the occupation number $\langle {\widehat{n}}_{{\omega }_i}\rangle$; this functional identity is the classical counterpart of the Bose-like condensation mechanism. As $s_{{\omega }_i}$ is increased above a threshold value $s_c$, the stationary solution bifurcates away from equipartition: an increasingly large fraction of the total energy is channelled into the lowest-frequency mode ${\omega }_0$, giving rise to classical phonon condensation. Numerical integration of the dimensionless form of Equation (14) confirms that this threshold sharpens with increasing number of modes N, approaching a sharp bifurcation in the thermodynamic limit [10].

Introducing $\tau =t{\omega }_0$, $y_{{\omega }_i}={\omega }_i\langle J_{{\omega }_i}\rangle /\left(k_B{T}_B\right)$, ${\alpha }_{{\omega }_i}={\omega }_i/{\omega }_0$ where ${\omega }_0={min}_{\omega }\omega$, the rate equations take the dimensionless form

$${\dot{y}}_{{\omega }_i}=S_{{\omega }_i}+B_{{\omega }_i}(1-y_{{\omega }_i})+\sum _{{\omega }_j\ne {\omega }_i}{C}_{{\omega }_i{\omega }_j}\left[{\displaystyle \frac{{\alpha }_{{\omega }_i}}{{\alpha }_{{\omega }_j}}}{y}_{{\omega }_j}-y_{{\omega }_i}+{\displaystyle \frac{{\alpha }_{{\omega }_j}-{\alpha }_{{\omega }_i}}{{\alpha }_{{\omega }_j}}}{y}_{{\omega }_i}{y}_{{\omega }_j}\right]+\left(\mathrm{quartic}\right),$$

(15)

with dimensionless control parameter $S=S_{{\omega }_i}$ (assumed mode-independent). The equilibrium solution is uniformly $y_{{\omega }_i}=1$; condensation corresponds to a stationary solution with $y_{{\omega }_0}\gg 1$ and $y_{{\omega }_i}\ll y_{{\omega }_0}$ for ${\omega }_i>{\omega }_0$.

4. Classical Phonon Condensation

4.1. Stationary States and Order Parameters

The stationary solutions (${\dot{y}}_{{\omega }_i}=0$) at small S are the equipartition fixed point $y_{{\omega }_i}=1$, invariant under permutations of mode labels. Above a threshold $S_c$ the nonlinear coupling breaks this permutation symmetry: energy flows preferentially toward ${\omega }_0$.

Two order parameters quantify this. The condensation index [20]

$$E_y={\displaystyle \frac{y_{{\omega }_0}/{\tilde{y}}_{{\omega }_0}}{{\sum }_{i=0}^N{y}_{{\omega }_i}}},{\tilde{y}}_{{\omega }_0}=1+S/B,$$

(16)

takes $E_y=0$ at equipartition and $E_y=1$ at full condensation. The ratio $p_1/p_0$ (where $p_i=y_{{\omega }_i}/{\sum }_j{y}_{{\omega }_j}$) goes from 1 to 0 at the same transition.

4.2. Numerical Evidence

The rate equations are numerically integrated for $N+1$ equally spaced modes ${\omega }_n={\omega }_0(1+n/N)$, starting from $y_{{\omega }_i}\left(0\right)=1$. Parameters: $B=1$, $C=0.1$, quartic coupling ${{\rm Y}}^{\left(1\right)}={{\rm Y}}^{\left(2\right)}={10}^{-4}$ [10].

Figure 1. Classical Fröhlich-like condensation. Normalized energy fractions $p_i$ in the normal modes vs. mode frequencies ${\omega }_i$, for $N=20$ modes. As the energy input rate S increases, deviations from equipartition grow progressively. Panels (a)–(d) correspond to $S=0.1$ (blue), 1 (green), 10 (purple), and 100 (pink). Equipartition gives equal bar heights; note the increasing concentration of energy in the lowest-frequency mode at large S. From Ref.[10]. 

Figure 1. Classical Fröhlich-like condensation. Normalized energy fractions $p_i$ in the normal modes vs. mode frequencies ${\omega }_i$, for $N=20$ modes. As the energy input rate S increases, deviations from equipartition grow progressively. Panels (a)–(d) correspond to $S=0.1$ (blue), 1 (green), 10 (purple), and 100 (pink). Equipartition gives equal bar heights; note the increasing concentration of energy in the lowest-frequency mode at large S. From Ref.[10]. 

As S increases, $p_0$ grows while all $p_i$ ($i>0$) decrease. The condensation index $E_y\left(S\right)$ and the ratio $(p_1/p_0)\left(S\right)$ display a crossover near $S\approx 10$ that sharpens progressively with N, with curves for different N crossing at a common point whose slope grows with N, suggestive of a sharp bifurcation in the large-N limit. This finite-size sharpening mirrors the behaviour of order-parameter curves near second-order equilibrium phase transitions.

Figure 2. Classical Fröhlich-like condensation. Condensation index ${\mathcal{E}}_y$ (left panel) and ratio $p_1/p_0$ (right panel) vs. energy input rate S, for increasing mode number $N_{\mathrm{sys}}$: from right to left curves $N=11,21,41,101,301$.At equipartition ${\mathcal{E}}_y=0$, $p_1/p_0=1$; at full condensation ${\mathcal{E}}_y=1$, $p_1/p_0=0$. The dashed oblique line marks the inflection tangent as a guide to a possible asymptotic bifurcation. From Ref.[10]. 

Figure 2. Classical Fröhlich-like condensation. Condensation index ${\mathcal{E}}_y$ (left panel) and ratio $p_1/p_0$ (right panel) vs. energy input rate S, for increasing mode number $N_{\mathrm{sys}}$: from right to left curves $N=11,21,41,101,301$.At equipartition ${\mathcal{E}}_y=0$, $p_1/p_0=1$; at full condensation ${\mathcal{E}}_y=1$, $p_1/p_0=0$. The dashed oblique line marks the inflection tangent as a guide to a possible asymptotic bifurcation. From Ref.[10]. 

The dynamics reveals a hierarchy: the highest-frequency modes lose energy first at increasing S; lower modes follow; eventually all modes deplete except ${\omega }_0$, which absorbs the surplus and saturates at $p_0\to 1$. The saturation is robust and mode-number independent at high drive.

4.3. Physical Interpretation

Classical phonon condensation is an instance of a nonequilibrium phase transition [21]: an open system with nonlinear internal couplings, dissipation, and external drive self-organises into a macroscopically ordered state when the energy-input rate exceeds a critical value. The genericity of this mechanism, depending only on nonlinear mode coupling, a thermal bath, and a source, not on molecular details, suggests the phenomenon should occur in real macromolecules. In the condensed state the lowest-frequency mode carries a coherent oscillation of the entire molecule, activating a giant oscillating dipole moment that is a necessary prerequisite for long-range electrodynamic forces [11].

5. Long-Range Resonant Electrodynamic Interactions

5.1. Classical Electrodynamic Hamiltonian for Two Dipoles

The electrodynamic interaction between two coherently oscillating macromolecules is modelled at the classical level as a pair of driven harmonic oscillators coupled through the retarded electromagnetic field of the intervening medium. Each macromolecule is represented by its dominant polar mode, an oscillating electric dipole ${\mu }_{A,B}\left(t\right)$ at natural frequency ${\omega }_{A,B}$ with effective charge-to-mass ratio ${\zeta }_{A,B}=Q_{A,B}^2/m_{A,B}$, where $Q_{A,B}$ and $m_{A,B}$ are the effective charge and mass of each dipole. This representation is appropriate because the low-frequency collective vibrations of macromolecules ($\omega /2\pi \sim 0.1$–1 THz) fall well within the classical regime $k_{B}T/\hslash \omega \gg 1$ at room temperature, and the dominant electrodynamic coupling involves the dipole moment of each molecule rather than higher multipoles.

The coupling is mediated by the electric field that each dipole radiates into the surrounding medium, here described by the frequency-dependent permittivity $\epsilon \left(\omega \right)$. The Debye screening that suppresses static electrostatic forces in the ionic cellular environment becomes ineffective at frequencies above $\sim 250$ MHz [11], so that the oscillating dipole fields can propagate without ionic screening at the sub-THz frequencies of interest. In the dipole approximation the equations of motion are [11]

$$\begin{array}{cc}\hfill {\ddot{\mu }}_A+{\omega }_A^2{\mu }_A& ={\zeta }_A{\mathbf{E}}_B({\mathbf{r}}_A,t),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\ddot{\mu }}_B+{\omega }_B^2{\mu }_B& ={\zeta }_B{\mathbf{E}}_A({\mathbf{r}}_B,t),\hfill \end{array}$$

(18)

where ${\mathbf{E}}_{A\left(B\right)}({\mathbf{r}}_{B\left(A\right)},t)$ is the electric field created by dipole A (resp. B) at the location of its partner. Damping terms ${\gamma }_{A,B}{\dot{\mu }}_{A,B}$ and anharmonic contributions are present in the full equations of motion [11] but are omitted here as the focus is on the conservative normal-mode structure from which the interaction potential is derived.

The field ${\mathbf{E}}_{A\left(B\right)}$ at the position of the partner is linear in the source dipole and is fully characterised, for a given separation r and frequency, by the electric susceptibility tensor ${\chi }_{ij}(r,\omega )$ of the medium, whose explicit form is obtained from the solution of the d’Alembert equation for the vector potential in Lorenz gauge [11]. Taking the z-axis along the intermolecular separation $\mathbf{r}$, the tensor is diagonal with elements

$${\chi }_{11}^{\prime }(r,\omega )={\chi }_{22}^{\prime }(r,\omega )=-{\displaystyle \frac{cos(\omega \sqrt{{\epsilon }^{\prime }}r/c)}{{\epsilon }^{\prime }\left(\omega \right)r^3}}\left(1-{\displaystyle \frac{{\omega }^2{\epsilon }^{\prime }\left(\omega \right)r^2}{c^2}}\right)+\dots ,{\chi }_{33}^{\prime }(r,\omega )=-2{\chi }_{11}^{\prime }(r,\omega )/\left(1-\dots \right),$$

(19)

where ${\chi }_{ii}^{\prime }$ denotes the real (in-phase) part of ${\chi }_{ii}$, which is the component that contributes to the conservative interaction potential. In the near zone ($r\ll c/\omega$) and for a non-dispersive background (${\epsilon }^{\prime }=\mathrm{const}$), Equation (19) reduces to the static dipole propagator ${\chi }_{ii}^{\prime }\sim {\sigma }_i/\left({\epsilon }^{\prime }{r}^3\right)$ (with ${\sigma }_{1,2}=-1$, ${\sigma }_3=2$), while at finite frequency the factor $\epsilon \left(\omega \right)$ renormalises this scaling and retardation corrections introduce oscillatory $r^{-2}$ and $r^{-1}$ terms whose importance grows with r.

The normal frequencies ${\omega }_N$ of the dissipation-free coupled system are obtained by substituting the Fourier ansatz ${\mu }_{A,B}\left(t\right)={\mu }_{A,B}^{\left(0\right)}{e}^{-i{\omega }_{N}t}$ into Eqs. (17)–() and requiring non-trivial solutions. The resulting secular condition is

$$\left({\omega }_A^2-{\omega }_N^2\right)\left({\omega }_B^2-{\omega }_N^2\right)={\zeta }_A{\zeta }_B{\left[{\chi }_{ii}^{\prime }(r,{\omega }_N)\right]}^2,$$

(20)

an implicit equation in ${\omega }_N$ because the right-hand side itself depends on ${\omega }_N$ through the medium response. This implicitness is not a minor technicality: as shown in Ref. [11], a naive explicit approximation of the right-hand side at zeroth order in the coupling leads to the erroneous conclusion, originally reached by Fröhlich, that long-range resonant interactions of the form $U\sim 1/r^3$ survive at thermal equilibrium. The exact treatment via suitable integration shows that these contributions cancel identically at thermal equilibrium, reducing U to a short-range $1/r^6$ potential.

The solutions of Equation (20) are obtained analytically by a contour-integral method in the complex plane: the Lagrange inversion theorem, combined with Rouché’s theorem, allows one to resolve the implicit ${\omega }_N$-dependence of the right-hand side and yields the normal-mode frequencies as controlled perturbative expansions in the coupling ${\zeta }_A{\zeta }_B{\left[{\chi }_{ii}^{\prime }\right]}^2$, uniformly valid across the resonant (${\omega }_A\simeq {\omega }_B$) and off-resonant (${\omega }_A\gg {\omega }_B$) regimes.. The key result is that the two regimes differ qualitatively in the range of the resulting interaction: off-resonance the potential scales as $1/r^6$ (short-range, van der Waals-like), whereas at resonance it scales as $1/r^3$ in the near zone and as $1/r$ in the far zone, long-range at all distances. It is this selectivity of the resonant interaction that makes it a candidate mechanism for specific biomolecular recognition.

5.2. Off-Resonance and Resonant Interaction Potentials

The normal-mode frequencies obtained from Equation (20) determine the interaction energy through the shift they impart to the bare frequencies ${\omega }_{A,B}$. The qualitative character of this shift depends critically on whether the two dipoles are in resonance or not.

Off-resonance regime.

When ${\omega }_A\gg {\omega }_B$, the right-hand side of Equation (20) is small compared with the unperturbed product ${({\omega }_A^2-{\omega }_B^2)}^2$, and the secular equation can be solved perturbatively. To leading order the frequency shift of each mode is proportional to ${\zeta }_A{\zeta }_B{\left[{\chi }_{ii}^{\prime }\right]}^2/({\omega }_A^2-{\omega }_B^2)$, i.e., second-order in the coupling ${\chi }_{ii}^{\prime }$. Since ${\chi }_{ii}^{\prime }\sim r^{-3}$ in the near zone, the resulting interaction potential scales as $U\left(r\right)\propto r^{-6}$ a short-range, van der Waals–like potential entirely analogous to the London dispersion interaction between non-resonant fluctuating dipoles.

Resonant regime.

At resonance (${\omega }_A\approx {\omega }_B={\omega }_0$) the structure of Equation (20) changes qualitatively. Setting ${\omega }_N={\omega }_0+\delta \omega$ and expanding to first order in $\delta \omega$, the left-hand side becomes ${({\omega }_0^2-{\omega }_N^2)}^2\approx 4{\omega }_0^2{\left(\delta \omega \right)}^2$, so the secular equation reduces to a perfect square and yields

$${\omega }_{i,\pm }\left(r\right)\simeq {\omega }_0\pm \sqrt{{\zeta }_A{\zeta }_B}{\displaystyle \frac{{\chi }_{ii}^{\prime }(r,{\omega }_0)}{2{\omega }_0}}.$$

(21)

The frequency shift is now first-order in ${\chi }_{ii}^{\prime }$, not quadratic. This is the key distinction from the off-resonance case: the square-root structure of the resonant secular equation lifts the second-order suppression and produces a qualitatively stronger coupling.

The total interaction energy is expressed in terms of the adiabatic action variables $J_{i,\pm }=E_{i,\pm }/{\omega }_{i,\pm }$ of the two normal modes, where $E_{i,\pm }$ is the energy stored in the ± branch of polarisation component i. Summing over the three spatial components,

$$U\left(r\right)={\displaystyle \sum _{i=1}^3}\Delta {\omega }_{0,i}\left(r\right)(J_{i,+}-J_{i,-}).$$

(22)

In the near zone ($r\ll c/{\omega }_0$) the susceptibility scales as ${\chi }_{ii}^{\prime }\sim r^{-3}$, so $\Delta {\omega }_{0,i}\sim r^{-3}$ and

$$U\left(r\right)\sim \pm {\displaystyle \frac{1}{r^3}},$$

(23)

a genuinely long-range interaction: the exponent drops from $-6$ in the off-resonance case to $-3$ at resonance, reflecting the fact that U is now linear rather than quadratic in the susceptibility ${\chi }_{ii}^{\prime }$. The sign and magnitude are controlled by the action difference $J_{i,+}-J_{i,-}$: at thermal equilibrium this difference vanishes at first order, and the long-range term disappears, recovering a $1/r^6$ free energy. Out of equilibrium, when phonon condensation concentrates energy in one normal mode, $J_{i,+}\gg J_{i,-}$ and the $1/r^3$ interaction is activated. The explicit behaviour of ${\chi }_{ii}^{\prime }(r,\omega )$ as a function of r for representative frequencies in the sub-THz and optical regimes is shown in Figure 1 of Ref. [11].

A crucial result is that resonant long-range interactions vanish at thermal equilibrium. At equilibrium the Boltzmann distribution weights the two normal-mode energies equally, causing the action difference $J_{i,+}-J_{i,-}$ in Equation (22) to vanish at first order, giving a short-range $1/r^6$ free energy [11]. Long-range resonant electrodynamic forces therefore require an out-of-equilibrium system, where energy concentration in one normal mode, precisely the condition delivered by phonon condensation, makes $J_{i,+}\gg J_{i,-}$.

6. Position-Space Classical Hamiltonian and Molecular Dynamics

The action-angle Hamiltonian of Sec. 3 is natural for perturbative rate equations but remote from $(q,p)$ variables used in molecular dynamics force fields. A complementary classical Hamiltonian was formulated directly in position space [12], opening a path toward connecting Fröhlich condensation to atomistic MD. This has motivated the study of the following model.

6.1. The Position-Space Hamiltonian

The system Hamiltonian is $H=H_0+H_{\mathrm{int}}$, with free part

$$H_0={\displaystyle \sum _{i=1}^N}\left({\displaystyle \frac{p_i^2}{2m_i}}+{\displaystyle \frac{1}{2}}{m}_i{\omega }_i^2{q}_i^2\right)+{\displaystyle \sum _{k=1}^{N_B}}\left({\displaystyle \frac{p_k^{\left(B\right)2}}{2m_k^{\left(B\right)}}}+{\displaystyle \frac{1}{2}}{m}_k^{\left(B\right)}{\omega }_k^{\left(B\right)2}{q}_k^{\left(B\right)2}\right)+{\displaystyle \sum _{l=1}^{N_S}}\left({\displaystyle \frac{p_l^{\left(S\right)2}}{2m_l^{\left(S\right)}}}+{\displaystyle \frac{1}{2}}{m}_l^{\left(S\right)}{\omega }_l^{\left(S\right)2}{q}_l^{\left(S\right)2}\right),$$

(24)

and interaction

$$H_{\mathrm{int}}=\sum _{ik}{\varphi }_{ik}{q}_i{q}_k^{\left(B\right)}+\sum _{il}{\xi }_{il}{q}_i{q}_l^{\left(S\right)}+\sum _{ijk}{\lambda }_{ijk}{q}_i{q}_j{q}_k^{\left(B\right)},$$

(25)

where ${\varphi }_{ik}$, ${\xi }_{il}$, ${\lambda }_{ijk}$ are linear protein-bath, linear protein-source, and nonlinear protein-protein-bath couplings. Second-order perturbation theory recovers the Fröhlich rate constants ${\Phi }_i$, ${\Xi }_i$, ${\Lambda }_{ij}$, all depending on resonance conditions [12]. Hamiltonian equations are integrated numerically with bath oscillators held at $T=300$ K by a Langevin thermostat [12]. The key findings are: i) Pure Fröhlich resonances (${\omega }_i-{\omega }_j\pm {\omega }_k^{\left(B\right)}=0$) produce robust condensation; Lifshits sum resonances (${\omega }_i+{\omega }_j-{\omega }_k^{\left(B\right)}=0$) destroy it. ii) Replacing the bath-mediated trilinear coupling ${\lambda }_{ijk}{q}_i{q}_j{q}_k^{\left(B\right)}$ by a protein-only coupling ${\lambda }_{ijk}{q}_i{q}_j{q}_k$ completely suppresses condensation: bath mediation is physically essential, not a technical device. iii) Strong condensates are found at $T=300$ K, contradicting an earlier study [20] that used constant coupling coefficients. The discrepancy traces to the coupling parametrisation: coefficients proportional to $\sqrt{m_i{m}_k^{\left(B\right)}{\omega }_i{\omega }_k^{\left(B\right)}}$ ensure all modes contribute equally to the interaction energy at equipartition.

Figure 3. Condensation index ρ of a Fröhlich system from Hamiltonian dynamics, for coupling coefficients ${\lambda }_{ijk}$ under three resonance conditions: ${\omega }_i-{\omega }_j\pm {\omega }_k^{\left(B\right)}=0$ (Fröhlich), ${\omega }_i+{\omega }_j-{\omega }_k^{\left(B\right)}=0$ (Lifshits), and ${\omega }_i\pm {\omega }_j\pm {\omega }_k^{\left(B\right)}=0$ (combined). From Ref.[12]. 

Figure 3. Condensation index ρ of a Fröhlich system from Hamiltonian dynamics, for coupling coefficients ${\lambda }_{ijk}$ under three resonance conditions: ${\omega }_i-{\omega }_j\pm {\omega }_k^{\left(B\right)}=0$ (Fröhlich), ${\omega }_i+{\omega }_j-{\omega }_k^{\left(B\right)}=0$ (Lifshits), and ${\omega }_i\pm {\omega }_j\pm {\omega }_k^{\left(B\right)}=0$ (combined). From Ref.[12]. 

Figure 4. Time evolution of total energies (kinetic + potential) in a Fröhlich system of nine protein modes (first five shown), computed as 300 ns moving averages from Hamiltonian dynamics. Parameters: bath $T=300$ K, source $T_S=3000$ K, protein frequencies ${\omega }_1=0.2$ THz to ${\omega }_9=1$ THz in 0.1 THz steps, ${\varphi }_{ik}$, ${\xi }_{il}$, ${\lambda }_{ijk}$ from Equation (5) of Ref.[12] with $\varphi =1.0$, $\xi =0.4$ (from 100 ns), unit masses, thermostat friction $0.1$ ps${-}^1$. Black curves: Fröhlich rate equations predictions from Equation (14) with $\alpha =0.02$ ps (solid: mode 1; dashed: higher modes). Left panel: $\lambda =0$; right panel: $\lambda =0.95$ K${-}^{1/2}$. From Ref.[12]. 

Figure 4. Time evolution of total energies (kinetic + potential) in a Fröhlich system of nine protein modes (first five shown), computed as 300 ns moving averages from Hamiltonian dynamics. Parameters: bath $T=300$ K, source $T_S=3000$ K, protein frequencies ${\omega }_1=0.2$ THz to ${\omega }_9=1$ THz in 0.1 THz steps, ${\varphi }_{ik}$, ${\xi }_{il}$, ${\lambda }_{ijk}$ from Equation (5) of Ref.[12] with $\varphi =1.0$, $\xi =0.4$ (from 100 ns), unit masses, thermostat friction $0.1$ ps${-}^1$. Black curves: Fröhlich rate equations predictions from Equation (14) with $\alpha =0.02$ ps (solid: mode 1; dashed: higher modes). Left panel: $\lambda =0$; right panel: $\lambda =0.95$ K${-}^{1/2}$. From Ref.[12]. 

7. Electron-Phonon Hamiltonian and DNA–Protein Co-Resonance

7.1. Motivation and Biological Context

Beyond collective vibrational (phonon) degrees of freedom, the electronic degrees of freedom of biomolecules provide a further and conceptually distinct channel through which selective long-range electrodynamic interactions can be generated [13]. This possibility is motivated by three converging lines of evidence. First, the Resonant Recognition Model (RRM) [22,23] showed empirically that a Fourier analysis of the electron-ion interaction potentials (EIIP) along a protein or DNA sequence produces spectral peaks that correlate with biological activity. Second, electron transport along DNA and along protein backbones is an experimentally confirmed phenomenon [24,25]. Third, a previous study [26] established that the electron current flowing along a DNA fragment under external energy supply can exhibit either a broad or a sharply peaked frequency spectrum depending on the excitation site and energy, suggesting that co-resonance between DNA and enzyme currents could provide a sequence-specific electrodynamic signal.

7.2. The Davydov–Holstein–Fröhlich Hamiltonian

To model electron–phonon motion along a biomolecular chain of N sites (nucleotides or amino acids), the following Hamiltonian is adopted [13]:

$$\widehat{H}={\widehat{H}}_{\mathrm{el}}+{\widehat{H}}_{\mathrm{ph}}+{\widehat{H}}_{\mathrm{int}},$$

(26)

with electronic part

$${\widehat{H}}_{\mathrm{el}}={\displaystyle \sum _{n=1}^N}\left[E_0{\widehat{B}}_n^{†}{\widehat{B}}_n+ϵ\langle {\widehat{B}}_n^{†}{\widehat{B}}_n\rangle {\widehat{B}}_n^{†}{\widehat{B}}_n+J_n\left({\widehat{B}}_n^{†}{\widehat{B}}_{n+1}+{\widehat{B}}_n^{†}{\widehat{B}}_{n-1}\right)\right],$$

(27)

phononic part (with anharmonic term)

$${\widehat{H}}_{\mathrm{ph}}={\displaystyle \frac{1}{2}}\sum _n\left[{\displaystyle \frac{{\widehat{p}}_n^2}{M_n}}+{\Omega }_n{({\widehat{u}}_{n+1}-{\widehat{u}}_n)}^2+{\displaystyle \frac{1}{2}}\mu {({\widehat{u}}_{n+1}-{\widehat{u}}_n)}^4\right],$$

(28)

and electron-phonon interaction

$${\widehat{H}}_{\mathrm{int}}=\sum _n{\chi }_n({\widehat{u}}_{n+1}-{\widehat{u}}_n){\widehat{B}}_n^{†}{\widehat{B}}_n.$$

(29)

Here ${\widehat{B}}_n$, ${\widehat{B}}_n^{†}$ are lowering and raising operators at lattice site n; $E_0$ is the initial electron excitation energy; $ϵ$ is the nonlinear electron self-trapping constant; $J_n$ is the site-dependent electron tunnelling amplitude; ${\widehat{u}}_n$, ${\widehat{p}}_n$ are the longitudinal displacement and momentum of the n-th nucleotide or amino acid, according to the kind of molecule, DNA or protein, respectively; ${\Omega }_n$ is the site-dependent spring constant; $M_n$ is the site mass; $\mu$ is the paqrticle-particle (nucleotide or amino acid) anharmonic coupling; and ${\chi }_n$ is the site-dependent electron-phonon coupling. The coupling parameters $J_n$ and ${\chi }_n$ are site-dependent, encoding the specific sequence of nucleotides (for DNA) or amino acids (for protein).

7.3. TDVP Applied to the Davydov Ansatz

The TDVP is again the instrument of choice. The trial wave function is written in the Davydov factorised form

$$\left|\psi \right(t)\rangle =|\Psi \left(t\right)\rangle \left|\Phi \right(t)\rangle ,$$

(30)

where the electronic part $\left|\Psi \right(t)\rangle$ is a single-excitation state

$$|\Psi \left(t\right)\rangle =\sum _n{C}_n\left(t\right){\widehat{B}}_n^{†}{|0\rangle }_{\mathrm{el}},$$

(31)

with probability amplitudes $C_n\left(t\right)$, and the phononic part $\left|\Phi \right(t)\rangle$ is a displaced coherent state

$$|\Phi \left(t\right)\rangle =e^{-{\displaystyle \frac{i}{\hslash }}{\sum }_n[{\beta }_n\left(t\right){\widehat{p}}_n-{\pi }_n\left(t\right){\widehat{u}}_n]}{|0\rangle }_{\mathrm{ph}},$$

(32)

so that $\langle \Phi |{\widehat{u}}_n|\Phi \rangle ={\beta }_n\left(t\right)$ and $\langle \Phi |{\widehat{p}}_n|\Phi \rangle ={\pi }_n\left(t\right)$.

Applying the stationary-action condition $\delta S=0$ with the Lagrangian $L\left(t\right)=i\hslash \langle \psi |{\partial }_t|\psi \rangle -\langle \psi |\widehat{H}|\psi \rangle$ yields equations of motion for the variational parameters $(C_n,C_n^{*},{\beta }_n,{\pi }_n)$:

$$i\hslash {\dot{C}}_n={\partial }_{C_n^{*}}H,{\dot{\beta }}_n={\partial }_{{\pi }_n}H,{\dot{\pi }}_n=-{\partial }_{{\beta }_n}H,$$

(33)

where $H=\langle \psi \left|\widehat{H}\right|\psi \rangle$ is the classical Hamiltonian obtained as the expectation value.

Evaluating H explicitly [13]:

$$\begin{array}{cc}\hfill H& =\sum _n[E_0|C_n{|}^2+ϵ{\left|C_n\right|}^4+J_n(C_n^{*}{C}_{n+1}+C_{n+1}^{*}{C}_n)\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{{\pi }_n^2}{M_n}}+{\Omega }_n{({\beta }_{n+1}-{\beta }_n)}^2+{\displaystyle \frac{\mu }{2}}{({\beta }_{n+1}-{\beta }_n)}^4\right)\hfill \\ \hfill & +{\chi }_n({\beta }_{n+1}-{\beta }_n){\left|C_n\right|}^2].\hfill \end{array}$$

(34)

The equations of motion that follow from Eqs. (33) and (34) are formally classical but describe the time evolution of quantum expectation values:

$$\begin{array}{cc}\hfill i\hslash {\dot{C}}_n& =\left[E_0+2ϵ{\left|C_n\right|}^2+{\chi }_n({\beta }_{n+1}-{\beta }_n)\right]C_n\hfill \\ \hfill & +J_n{C}_{n+1}+J_{n-1}{C}_{n-1},\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill M_n{\ddot{\beta }}_n& ={\Omega }_n{\beta }_{n+1}+{\Omega }_{n-1}{\beta }_{n-1}-{\Omega }_{n-1}{\beta }_n-{\Omega }_n{\beta }_n\hfill \\ \hfill & +{\chi }_n|C_n{|}^2-{\chi }_{n-1}{\left|C_{n-1}\right|}^2\hfill \\ \hfill & +\mu \left[{({\beta }_{n+1}-{\beta }_n)}^3-{({\beta }_n-{\beta }_{n-1})}^3\right].\hfill \end{array}$$

(36)

Equations (35)–(36) couple the quantum probability amplitudes for electron propagation to the classical phonon displacements of the chain, encoding the full electron–phonon dynamics of the specific biomolecular sequence.

7.4. Physical Parameters and Numerical Strategy

In Ref. [13] the DNA–protein system studied is a 66 bp DNA oligonucleotide interacting with the EcoRI restriction enzyme (276 amino acids). The site-dependent parameters $J_n$ and ${\chi }_n$ are determined from the electron-ion interaction potentials (EIIP) of each nucleotide (A, T, G, C) and each amino acid, tabulated in Refs. [22,23,27]. Tunnelling amplitudes are estimated from the transmission probability $P(n\to n\pm 1)$ through the potential barrier between adjacent sites; the site-dependent electron-phonon coupling is estimated as ${\chi }_n=dE/dx=(E_{n+1}-E_n)/a$, where a is the nearest-neighbour distance between particles ($a=3.4$ Å for nucleotides of DNA, $4.5$ Å for the amino acids of enzyme).

The particle displacements are initialised at thermal equilibrium at $T=310$ K:

$$\langle |b_n{\left(0\right)|\rangle }_n=\sqrt{{\displaystyle \frac{k_{B}T}{\hslash \omega {\Omega }^{\prime }}}},\langle |{\pi }_n\left(0\right){|\rangle }_n=\sqrt{{\displaystyle \frac{k_{B}T}{\hslash \omega }}}.$$

(37)

The electron wavefunction is initialised as a localised sech pulse centred at site $n_0$. Integration uses a symplectic leapfrog scheme for the phonon variables combined with an Euler predictor-corrector for the electron amplitudes, with relative energy conservation $\Delta E/E\sim {10}^{-6}$.

7.5. Cross-spectral Co-Resonance and Sequence Specificity

The primary observable is the cross Fourier spectrum of the electron current densities along the two chains. The average electron current along chain $\alpha$ ($\alpha =1$ for DNA, $\alpha =2$ for EcoRI) is

$$i_{\alpha }\left(t\right)={\displaystyle \frac{e\hslash }{2N_{\alpha }{a}_{\alpha }{m}_{e}i}}{\displaystyle \sum _{j=1}^{N_{\alpha }}}\left({\Psi }_{\alpha }^{*}(x_j,t){\displaystyle \frac{{\Psi }_{\alpha }(x_{j+1},t)-{\Psi }_{\alpha }(x_{j-1},t)}{2}}-\mathrm{c}.\mathrm{c}.\right),$$

(38)

and the cross-spectrum ${\tilde{i}}_1^{*}\left(\nu \right){\tilde{i}}_2\left(\nu \right)$ is computed from their Fourier transforms.

The central finding [13] is illustrated schematically as follows. When the EcoRI recognition site on the DNA is the canonical palindromic sequence $3^{\prime }$-CTTAA|G-$5^{\prime }$ (where | marks the cleavage site), the cross-spectrum displays a sharp co-resonance peak at $\nu \approx 20$–29 THz (depending on initial conditions), in quantitative agreement with the RRM prediction [22].

Figure 5. Cross frequency spectra between a DNA strand with $N_1=66$ nucleotides and theEcoRI enzyme with $N_2=276$ amino acids. Panels show DNA strand containing a) the canonical target CTTAAG recognition site, b) cyclic permutation of restriction site in a) to AGCTTA, c) one-nucleotide change from a) with its base-pair complement to CATAAG, and d) two-nucleotide change from a) with their base-pair complements to GTTAAC. From Ref.[13]. 

Figure 5. Cross frequency spectra between a DNA strand with $N_1=66$ nucleotides and theEcoRI enzyme with $N_2=276$ amino acids. Panels show DNA strand containing a) the canonical target CTTAAG recognition site, b) cyclic permutation of restriction site in a) to AGCTTA, c) one-nucleotide change from a) with its base-pair complement to CATAAG, and d) two-nucleotide change from a) with their base-pair complements to GTTAAC. From Ref.[13]. 

When the recognition sequence is randomised (e.g., AGCTTA, TCATGA, AGATCT), the sharp peak disappears and is replaced by a broad, noisy spectrum. When one nucleotide of the recognition site is exchanged with its complement (e.g., CATAAG), the peak undergoes a modest broadening; two exchanges (e.g., GTTAAC) broaden it more. The co-resonance is therefore robustly and specifically associated with the canonical recognition sequence.

Figure 6. Cross frequency spectra between a DNA strand and theEcoRI enzyme, under the same substrate and initial conditions of Figure 5. Panels show DNA strand containing: a) randomized restriction site TCATGA; b) one-nucleotide change from the canonical target CTTAAG to CTTAAC; c) two-nucleotide change from the canonical target CTTAAG to CATATG. From Ref.[13]. 

Figure 6. Cross frequency spectra between a DNA strand and theEcoRI enzyme, under the same substrate and initial conditions of Figure 5. Panels show DNA strand containing: a) randomized restriction site TCATGA; b) one-nucleotide change from the canonical target CTTAAG to CTTAAC; c) two-nucleotide change from the canonical target CTTAAG to CATATG. From Ref.[13]. 

Figure 7. Cross frequency spectra between a DNA strand and theEcoRI enzyme, under the same substrate and initial conditions of Figure 5. Panels show DNA strand containing a e) randomized restriction site AGATCT, f) one-nucleotide change from the canonical target CTTAAG to CATAAG, g) two-nucleotide change from the canonical target CTTAAG to GTTATG. From Ref.[13]. 

Figure 7. Cross frequency spectra between a DNA strand and theEcoRI enzyme, under the same substrate and initial conditions of Figure 5. Panels show DNA strand containing a e) randomized restriction site AGATCT, f) one-nucleotide change from the canonical target CTTAAG to CATAAG, g) two-nucleotide change from the canonical target CTTAAG to GTTATG. From Ref.[13]. 

The model also discriminates between point-mutation variants with known biological activity [13]: promiscuous mutations (Ala¹³⁸ →Thr, Glu¹⁹² →Lys, His¹¹⁴ →Tyr) that allow binding to near-cognate sequences still yield a sharp co-resonance peak, while the catalytically inactivating mutation Asp⁹¹ →Asn produces a sizeable decrease of the peak. These findings are coherent with the experimental biochemical data.

Figure 8. Cross frequency spectra between the original substrate and theEcoRI mutants with the single point mutations; a) $Al{a}^{138}\to Thr$; b) $Gl{u}^{192}\to Lys$; c) $Hi{s}^{114}\to Tyr$; and d) $As{p}^{91}\to Asn$. The initial conditions are the same of Figure 5. From Ref.[13]. 

Figure 8. Cross frequency spectra between the original substrate and theEcoRI mutants with the single point mutations; a) $Al{a}^{138}\to Thr$; b) $Gl{u}^{192}\to Lys$; c) $Hi{s}^{114}\to Tyr$; and d) $As{p}^{91}\to Asn$. The initial conditions are the same of Figure 5. From Ref.[13]. 

7.6. Physical Interpretation and Comparison with the Vibrational Channel

The TDVP-derived equations of motion (35)–(36) are formally classical but arise directly from a quantum Hamiltonian via the same variational dequantization procedure applied in Sec. 3 to the Wu–Austin model. The two applications of the TDVP operate on different physical degrees of freedom (collective vibrational modes vs. electron-phonon propagation) and at different energy scales (sub-THz vs. tens of THz), but share the same mathematical structure: a product ansatz over coherent (or single-excitation) states, a classical Lagrangian obtained as the expectation value of the Hamiltonian, and equations of motion for the classical parameters governing the trial state.

The co-resonance reported here constitutes a distinct and complementary mechanism for selective electrodynamic interaction between biomolecules. In the vibrational channel (Secs. Sec. 4Sec. 5), the relevant quantity is the oscillating dipole moment carried by a collective phonon mode driven into a condensed state by external energy input. In the electronic channel (this section), the relevant quantity is the oscillating electron current density whose spectral overlap with that of the partner molecule determines the strength of the co-resonant coupling.

Both channels require out-of-equilibrium conditions for their activation: the phonon channel requires energy input exceeding the condensation threshold; the electronic channel requires electron excitation, provided in the cell, for example, by redox reactions, ATP-driven ion currents, or biophoton emission [13]. It is plausible that both channels are simultaneously active in vivo, and that they reinforce each other in determining the selectivity and range of biomolecular recognition.

The mediating role of the aqueous environment is highlighted by the Jaynes–Cummings-like interaction Hamiltonian [13] $H_{\mathrm{int}}=\hslash \sqrt{N}\gamma (a^{†}{S}^{-}+a{S}^{+})$, where $a^{†}$, a are the creation and annihilation operators for the DNA radiative electric field, $S^{\pm }$ are raising/lowering operators for the collective orientational state of the N surrounding water dipoles, and the coupling scales as $\sqrt{N}$, suggesting that the large number of water molecules in a physiological domain provides a protective gap against thermalization for the long-range correlations [13,28].

8. Experimental Evidence

8.1. Experimental Strategy

The unprecedented experimental programme outlined in the present section has been focused on two observables predicted by theory. First, an out-of-equilibrium macromolecule undergoing phonon condensation should display a sharp sub-THz absorption feature absent at thermal equilibrium, with a threshold in the energy-input rate [10]. Second, if two such macromolecules interact via resonant electrodynamic forces, the collective oscillation frequency of each should shift as $\Delta \nu \propto 1/{\langle r\rangle }^3$ [11,30], and above a critical concentration a clustering phase transition should occur [30].

Figure 9. Long-range electrodynamic interactions - Principle and experimental approaches At thermal equilibrium, macromolecules show a Brownian diffusive motion in solution (left panel). By switching-on an external energy source, molecules are in an out-of-thermal equilibrium collective vibrational state that can generate ED forces through associated large dipolar resonant oscillations (right panel). From Ref.[30].

Figure 9. Long-range electrodynamic interactions - Principle and experimental approaches At thermal equilibrium, macromolecules show a Brownian diffusive motion in solution (left panel). By switching-on an external energy source, molecules are in an out-of-thermal equilibrium collective vibrational state that can generate ED forces through associated large dipolar resonant oscillations (right panel). From Ref.[30].

Out-of-equilibrium conditions are achieved by optical pumping of fluorochrome-labelled proteins with a continuous-wave Argon laser (488 nm), inducing “proteinquake” energy transfer to the vibrational modes. Two complementary THz near-field spectroscopy setups, a microwire probe and a plasmonic rectenna, operated in two independent laboratories, together with fluorescence correlation spectroscopy (FCS) serve as detectors.

8.2. Phonon Condensation in BSA

The first experimental test of the classical phonon condensation scenario was performed with bovine serum albumin (BSA), a globular protein of 67 kDa chosen as a well-characterised model system [10]. Because BSA has no endogenous chromophore suitable for efficient optical pumping, five to six Alexa 488 fluorochrome molecules were covalently attached to its lysine residues and excited by an argon-ion laser at 488 nm. The energy difference between absorbed and re-emitted photons ($\approx 0.19$ eV per fluorochrome per incident photon) provides the energy input that drives the protein out of thermal equilibrium via the proteinquake mechanism [29].

THz near-field absorption spectroscopy was performed simultaneously in two independent laboratories using complementary probe technologies: a $12\mu$m-diameter microcoaxial wire (Montpellier) and a plasmonic bow-tie rectenna coupled to a plasma-wave FET (Rome). Both setups operate in the $0.22$–$0.33$ THz range, with a spectral resolution below 300 Hz from continuous-wave Virginia Diodes sources. A principal resonance at $\nu =0.314$ THz was reproducibly observed in both laboratories only when (i) Alexa 488 was bound to BSA and (ii) the laser was switched on; control measurements on labelled BSA without laser illumination, on free dye without protein, and on unlabelled BSA with illumination all showed no spectral feature.

The resonance frequency is in quantitative agreement with the lowest spheroidal deformation mode of an elastic sphere,

$${\nu }_0={\displaystyle \frac{1}{2\pi }}{\left[{\displaystyle \frac{2(2l+1)(l-1)E}{\rho R_H^2}}\right]}^{1/2},l=2,$$

(39)

giving ${\nu }_0\approx 0.308$ THz using the independently measured Young modulus $E=6.75$ GPa and hydrodynamic radius $R_H=35$ Å [10].

A complementary, rough estimate treats the BSA as a sphere splitted into two parts of masses $m=33$ kDa joined by a spring of stiffness $k=E{A}_0/{\ell }_0$, giving $\nu \approx 0.300$ THz [10], within 5% of the observed resonance, providing independent confirmation that the feature is a global deformation mode of the whole molecule. The two weaker resonances observed at $0.278$ and $0.285$ THz can be tentatively assigned to torsional modes at frequencies ${\nu }_t={\nu }_0{[(2l+3)/\left(2(2l+1)\right)]}^{1/2}$ for $l=2$ and $l=3$, which yield $0.257$ and $0.246$ THz respectively; the remaining discrepancy is consistent with the non-spherical shape of BSA introducing different moments of inertia along different axes. The absorption line profile is well described by the Lorentzian expected for a damped harmonic oscillator (Equation 52 of Ref. [10]), with a quality factor $Q=\Delta \nu /\nu \approx 50$.

The threshold and saturation behaviour of the resonance intensity as a function of laser power matches the theoretical prediction from the nonlinear rate equations (15).

A consistency estimate of the onset timescale (several minutes of illumination before the resonance reaches its final amplitude) can be obtained from the energy balance

$${\displaystyle \frac{dE}{dt}}=-{\displaystyle \frac{2}{3}}{\displaystyle \frac{{\left(Ze\right)}^2{\left|\ddot{x}\right|}^2}{c^3}}-\Gamma +W,$$

(40)

where W is the optical energy input rate and the first term on the right is the bremsstrahlung (radiative) loss. With five to six Alexa 488 molecules per BSA each receiving roughly 120–300 photons per second at $500\mu$W of laser power, the upper bound for W is $3.8$–$9.5\times {10}^{-11}$ erg s${-}^1$, almost coinciding with the estimated radiative loss of $3.25\times {10}^{-12}$–$1.7\times {10}^{-11}$ erg s${-}^1$ [10]. Because the two rates are so nearly balanced, energy accumulates slowly in each molecule, and several minutes elapse before the condensation threshold is crossed. Once the condensed phase is established, the collective low-frequency oscillation carries a very large effective dipole moment: requiring the bremsstrahlung losses of the condensed oscillation to be balanced by W implies an oscillating dipole in the range $14,500$–$23,000$ Debye, corresponding to an effective charge $Z\approx 290$–460 elementary charges [10]. This is fully consistent with the strong THz absorption feature emerging against the water background ($\approx 2000$ dB/cm extinction).

Finally, we note the biological relevance of these estimates. In living cells, a protein molecule undergoes roughly ${10}^6$ collisions per second with ATP molecules at a standard intracellular concentration of $\sim 1$ mM. If even 1–2% of those collisions transfer the $\approx 50$ kJ/mol of ATP hydrolysis free energy to the protein, the available power is $\sim 8\times {10}^{-9}$–$1.6\times {10}^{-8}$ erg s${-}^1$ per molecule, at least two orders of magnitude larger than what was delivered by the laser in the in vitro experiment [10]. Phonon condensation in vivo could therefore be activated on a far shorter timescale than the minutes required in the laboratory.

Figure 10. Differential transmission and absorption spectra as functions of the frequency. Comparison of the two normalized spectra for the longest illumination durations. From Ref.[10]. 

Figure 10. Differential transmission and absorption spectra as functions of the frequency. Comparison of the two normalized spectra for the longest illumination durations. From Ref.[10]. 

Figure 11. Threshold-like behaviour of giant dipolar oscillations. (a) Intensity of the resonant peak measured at $0.314$ THz as a function of the optical laser power. (b) Normalized energy of the fundamental mode calculated as a function of the normalized source power. The different curves correspond to the reported numbers of normal modes of the BSA protein. Theory and experiment are inqualitativeagreement. From Ref.[10]. 

Figure 11. Threshold-like behaviour of giant dipolar oscillations. (a) Intensity of the resonant peak measured at $0.314$ THz as a function of the optical laser power. (b) Normalized energy of the fundamental mode calculated as a function of the normalized source power. The different curves correspond to the reported numbers of normal modes of the BSA protein. Theory and experiment are inqualitativeagreement. From Ref.[10]. 

8.3. Phonon Condensation in R-PE

R-phycoerythrin (R-PE), a hexameric light-harvesting protein of molecular weight 240 kDa derived from red algae, provides a particularly clean experimental system because it absorbs strongly at 488 nm through its 38 endogenous phycoerythrin and phycourobilin fluorochromes, without requiring external labelling [30]. Under laser illumination with a 488 nm source, a collective oscillation mode appears at ${\nu }_0=71$ GHz ($2.4$ cm${-}^1$), displaying both a threshold in the energy input rate and a saturation of the oscillation amplitude at high power, in full analogy with the BSA observations [10]. The collective mode frequency is consistent with the lowest extension mode of a torus of major radius $R=37.5$ Å and minor radius $r=30$ Å, yielding a Young modulus $E\approx 5.3$ GPa by inversion of the Blevins formula [30], a value that lies squarely between those of myoglobin ($3.5$ GPa) and BSA ($6.75$ GPa), both predominantly $\alpha$-helical proteins like R-PE.

Two collective extension modes appear at 71 GHz and 96 GHz, with the frequency ratio $96/71\approx 1.35$ matching the torus-mode theoretical ratio $\sqrt{2}\approx 1.41$ within $4\%$.

The two-mode structure carries an important mechanistic signature of phonon condensation. The second mode at 96 GHz is observable only after the first mode at 71 GHz has reached saturation [30]. This is precisely what the nonlinear rate equations predict: the condensation mechanism channels all injected energy into the fundamental mode first; only once that mode is saturated does energy become available to excite the next collective mode. In the in vitro experiment, the onset time for the 71 GHz peak is about 1 min (at 50 mW), while approximately 15 min of pumping is required to excite the 96 GHz peak. This hierarchy of timescales is fully consistent with the condensation picture and constitutes an additional internal consistency check on the interpretation.

Figure 12. R-PE coherent vibrational states. Comparison of the R-PE in saline solution (purple) and saline solution without R-PE (black). Two collective extension modes of R-PE appear at 71 GHz and 96 GHz. Experimental data (full circles); Lorentz fit (solid line). From Ref.[30].

Figure 12. R-PE coherent vibrational states. Comparison of the R-PE in saline solution (purple) and saline solution without R-PE (black). Two collective extension modes of R-PE appear at 71 GHz and 96 GHz. Experimental data (full circles); Lorentz fit (solid line). From Ref.[30].

8.4. Clustering Phase Transition

The most striking macroscopic consequence of the activated electrodynamic forces is a first-order clustering phase transition in R-PE solutions, detected by fluorescence correlation spectroscopy (FCS) [30]. The normalised diffusion coefficient $D/D_0$ measured by FCS is constant across all concentrations at $50\mu$W (below threshold): purely Brownian diffusion. At $100\mu$W and $150\mu$W, $D/D_0$ drops by several orders of magnitude at average intermolecular distances of $\approx 725$ Å and $\approx 925$ Å respectively. The transition point shifts to larger distances with increasing power (stronger oscillations, stronger forces), and is fully reversible upon switching off the laser.

Figure 13. Effect of protein concentration and laser power illumination on R-PE diffusion: Clustering phase transition. (a) Diffusion coefficients normalized to the Brownian $D_0$ values measured for each data series at 0.223 $\mu M$ ($\langle r\rangle \simeq 1950\mathring{A}$) and recorded at 50 $\mu W$ (blue circles), 100 $\mu W$ (green triangles), and 150 $\mu W$ (red triangles). Each point corresponds to the average of 5 independent experiments. (b) Clustering transition driven by long-range electrodynamic interparticle forces, obtained from Monte Carlo simulations by varying the initial mean separation $\langle r\rangle$. For ${\langle r\rangle }_{\mathrm{init}}\simeq 1000$ Å the system remains in the dispersed phase ($D/D_0=1$, Brownian diffusion). For ${\langle r\rangle }_{\mathrm{init}}\simeq 950$ Å it switches to the clustered phase ($D/D_0\ll 1$). Lower right panel: Molecular Dynamics computation of self-diffusion coefficient $D/D_0$ (normalized to the Brownian value $D_0$) versus intermolecular average distance $\langle d\rangle$ for a system of particles in a cubic box interacting through long-range electrodynamic forces. From Ref.[30].

Figure 13. Effect of protein concentration and laser power illumination on R-PE diffusion: Clustering phase transition. (a) Diffusion coefficients normalized to the Brownian $D_0$ values measured for each data series at 0.223 $\mu M$ ($\langle r\rangle \simeq 1950\mathring{A}$) and recorded at 50 $\mu W$ (blue circles), 100 $\mu W$ (green triangles), and 150 $\mu W$ (red triangles). Each point corresponds to the average of 5 independent experiments. (b) Clustering transition driven by long-range electrodynamic interparticle forces, obtained from Monte Carlo simulations by varying the initial mean separation $\langle r\rangle$. For ${\langle r\rangle }_{\mathrm{init}}\simeq 1000$ Å the system remains in the dispersed phase ($D/D_0=1$, Brownian diffusion). For ${\langle r\rangle }_{\mathrm{init}}\simeq 950$ Å it switches to the clustered phase ($D/D_0\ll 1$). Lower right panel: Molecular Dynamics computation of self-diffusion coefficient $D/D_0$ (normalized to the Brownian value $D_0$) versus intermolecular average distance $\langle d\rangle$ for a system of particles in a cubic box interacting through long-range electrodynamic forces. From Ref.[30].

The quantitative extent of the transition is remarkable: at an average intermolecular distance of 600 Å, the variance of fluorescence fluctuations increases by more than five orders of magnitude compared with the same quantity at 1950 Å, and the measured diffusion time ${\tau }_D$ jumps from $4.8\times {10}^{-5}$ s to $4.9$ s, a factor of $\sim {10}^5$ [30]. These values cannot be explained by any trivial heating effect: a thermal increase of diffusivity would shift $D_0$ uniformly and monotonically with concentration, producing no step-like feature. The molecular dynamics simulations with a $1/r^3$ interaction potential reproduce the measured critical distance at $150\mu$W quantitatively, providing independent theoretical support for the experimental interpretation. Monte Carlo computations confirm the formation of dense clusters below the critical intermolecular distance [30]. Direct visual evidence was obtained by confocal fluorescence microscopy: clusters of R-PE molecules appear and grow when the laser power exceeds the threshold, then rapidly dissolve when the laser is switched back to a power below the threshold for activating collective oscillations [30]. This reversibility on a timescale of seconds rules out irreversible photochemical aggregation and is fully consistent with the dynamical, energy-dependent nature of the underlying ED forces.

The fact that BSA does not show the diffusion-based clustering signal is itself informative and consistent with the theory. As discussed above, BSA requires approximately 10 min of optical pumping to activate a sharp collective oscillation. In the FCS geometry, the transit time of each protein through the confocal volume is far shorter than this activation time, so the collective oscillations cannot be established during the measurement. By contrast, R-PE, with its 38 endogenous fluorochromes and natural light-harvesting efficiency, reaches the condensed state in less than 3 s under comparable power density, making FCS measurement of its clustering feasible [30]. BSA is instead uniquely suited to the THz frequency-shift measurement, which uses long-exposure illumination of a static drop of solution rather than molecules in transit. Thus the two proteins play complementary roles and together constitute a self-consistent body of evidence.

These experimental outcomes are in quantitative agreement with the theoretical prediction of a clustering transition driven by the resonant $-1/r^3$ electrodynamic interaction potential, activated out of thermal equilibrium [30].

8.5. Frequency Shifts: Further Evidence for Long-Range Electrodynamic Forces

If resonant electrodynamic forces are active between collective dipoles oscillating at frequency ${\omega }_0$, the vibration frequency of each molecule shifts by $\Delta \nu \propto 1/{\langle r\rangle }^3\propto \mathcal{C}$ (concentration), where $\langle r\rangle ={\mathcal{C}}^{-1/3}$ is the average intermolecular distance [11,30]. This prediction is confirmed for both R-PE (at 71 GHz) and BSA (at $0.314$ THz): the measured frequency shift varies linearly with concentration, with a slope that increases with laser power (stronger oscillation amplitude implies larger dipole moment and hence stronger interaction). Experiments were performed in 200 mM NaCl solution, ensuring Debye screening of all electrostatic interactions.

For R-PE, frequency-shift measurements were performed at three laser powers: $31.5$, $39.5$, and 50 mW. In all cases the shift is linear in concentration, confirming the $1/{\langle r\rangle }^3$ dependence over the range of average intermolecular distances from roughly 700 to 1100 Å explored in the THz experiments [30]. The steepening of the slope with laser power is also quantitatively consistent with theory: a larger oscillation amplitude generates a larger oscillating dipole moment and hence a stronger electrodynamic coupling. The BSA frequency-shift measurement (at 40 mW) confirms the same linear law at $0.314$ THz, extending the evidence for long-range electrodynamic forces to a structurally unrelated protein. Crucially, the THz spectroscopy setup for these frequency-shift measurements used a lower spatial density of laser light than the FCS diffusion experiments, specifically to prevent a clustering transition during the measurement: the ED forces were deliberately weakened to remain below the threshold for cluster formation, so that the proteins kept diffusing freely while still experiencing the frequency shift [30]. This experimental design decouples the two observable signatures of ED forces and rules out cross-contamination of the two measurements.

Figure 14. Frequency shifts of the intramolecular collective vibrations of R-PE and BSA at different concentrations. Measurements were performed at room temperature in aqueous solution with 200 mM of NaCl. Panel (a) refers to R-PE. The shift is relative to the reference frequency measured at the lowest protein concentration. Measurements have been performed at different powers of the laser: 31.5 mW (blue circles), 39.5 mW (green triangles), 50 mW (red triangles). Purple squares and orange stars refer to theoretical outcomes worked out with different values of molecular dipole moments (see Supplementary Materials). Panel (b) refers to BSA at a laser power of 40 mW. Insets: The frequency shifts are measured through a Lorentz fitting of the experimental resonances. The different colors in the insets are just a visual help. From Ref.[30].

Figure 14. Frequency shifts of the intramolecular collective vibrations of R-PE and BSA at different concentrations. Measurements were performed at room temperature in aqueous solution with 200 mM of NaCl. Panel (a) refers to R-PE. The shift is relative to the reference frequency measured at the lowest protein concentration. Measurements have been performed at different powers of the laser: 31.5 mW (blue circles), 39.5 mW (green triangles), 50 mW (red triangles). Purple squares and orange stars refer to theoretical outcomes worked out with different values of molecular dipole moments (see Supplementary Materials). Panel (b) refers to BSA at a laser power of 40 mW. Insets: The frequency shifts are measured through a Lorentz fitting of the experimental resonances. The different colors in the insets are just a visual help. From Ref.[30].

9. Discussion and Outlook

A distinctive feature of the theoretical framework reviewed here is that the TDVP serves as a single unifying instrument across different physical scenarios. In Sec. 3, applied to product coherent states of bosonic operators, it produces a classical Hamiltonian in action-angle variables directly suitable for statistical (Liouville/KvN) analysis. In Sec. 7, applied to the Davydov factorised ansatz of electronic and phononic coherent states, it produces classical equations of motion for electron probability amplitudes coupled to classical phonon degrees of freedom, encoding the specific biomolecular sequence at the level of site-dependent coupling constants. In both cases the resulting dynamics is formally classical but directly inherits the quantum structure of the original Hamiltonian, without invoking any semiclassical approximation beyond the mean-field restriction to the chosen ansatz manifold.

9.1. Two Complementary Channels of Selective Long-Range Forces

The review presents two physically distinct but mathematically parallel channels through which selective long-range electrodynamic interactions can be activated between biomolecules:

The vibrational (phonon) channel: driven by metabolic or photonic energy input, a macromolecule undergoes classical Fröhlich condensation, concentrating energy into a coherent collective oscillation. The resulting giant oscillating dipole moment generates $1/r^3$ resonant forces, confirmed experimentally by THz spectroscopy and FCS.

The electronic channel: electron–phonon excitation along a specific DNA sequence or protein produces oscillating electron currents whose Fourier spectrum carries the “fingerprint” of the biomolecular sequence. The cross-spectrum of two interacting molecules displays a sharp co-resonance peak when and only when their sequences are cognate, providing a sequence-specific electrodynamic interaction mechanism.

The two channels operate at different energy and frequency scales (sub-THz vs. tens of THz) and are not mutually exclusive; in a living cell they may well operate simultaneously and synergistically.

9.2. Open Questions

From light pumping to ATP.

The experiments use laser excitation, which is a controlled but artificial energy supply. ATP hydrolysis provides an estimated power of $\sim {10}^{-8}$–${10}^{-7}$ erg s${-}^1$ per protein [10] orders of magnitude larger than the laser power available per protein in the experiments. The position-space Hamiltonian framework [12], via normal-mode analysis of force fields, offers a path toward modelling ATP-coupled condensation in atomistic detail.

Role of water.

The hydration shell plays an active role in both channels: it contributes to the effective dipole moment in the phonon channel and provides a Tavis–Cummings-like protective environment for electronic correlations in the DNA–protein channel [13,28]. A quantitative theory of the water contribution remains to be developed.

In vivo evidence.

Whether these phenomena operate in the complex, crowded, and dynamically heterogeneous cellular environment is the central open question. Selectivity of electrodynamic forces in a mixture of many molecular species, the role of molecular crowding, and the competition between ED forces and Brownian diffusion all require further investigation.

Quantum corrections.

The classical framework is quantitatively justified at sub-THz frequencies and room temperature, but quantum fluctuations in the condensed state [32] and quantum coherence in the electronic channel may be non-negligible and merit further study.

Connection to biomolecular condensates.

The first-order clustering transition driven by electrodynamic forces provides a physically distinct mechanism for the formation of reversible biomolecular assemblies [1], complementing models based on intrinsically disordered regions and multivalency. Whether this mechanism is relevant to phase-separated cellular compartments is an open and intriguing question.

9.3. Concluding Remarks

The programme reviewed here demonstrates that classical Hamiltonian mechanics, derived systematically from quantum models via the TDVP and implemented at the level of both rate equations and direct Hamiltonian simulation, produces rich and experimentally verifiable physics relevant to the organisation of biochemical reactions in living matter. Fröhlich’s half-century-old intuition has acquired rigorous Hamiltonian underpinning, direct experimental support, and a new electronic dimension through the DNA–protein co-resonance phenomenology. The next challenge is to establish whether and how these phenomena operate in the cellular environment, and what role they play in the extraordinary efficiency of biological molecular recognition.