# Out-of-Equilibrium Collective Oscillation as Phonon Condensation in a Model Protein

Ilaria Nardecchia, Jeremie Torres, Mathias Lechelon, Valeria Giliberti, Michele Ortolani, Philippe Nouvel, Matteo Gori, Yoann Meriguet, Irene Donato, Jordane Preto, et al.

► To cite this version:

Ilaria Nardecchia, Jeremie Torres, Mathias Lechelon, Valeria Giliberti, Michele Ortolani, et al. Out-of-Equilibrium Collective Oscillation as Phonon Condensation in a Model Protein. Physical Review X, 2018, 8 (3), pp.031061. 10.1103/physrevx.8.031061. hal-01994697

HAL Id: hal-01994697

https://hal.science/hal-01994697

Submitted on 25 Jan 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Out-of-Equilibrium Collective Oscillation as Phonon Condensation in a Model Protein

Ilaria Nardecchia,<sup>1,2,3,*</sup> Jeremie Torres,<sup>4,†</sup> Mathias Lechelon,<sup>1,2,3,‡</sup> Valeria Giliberti,<sup>5,§</sup> Michele Ortolani,<sup>5,¶</sup>

Philippe Nouvel,<sup>4,**</sup> Matteo Gori,<sup>4,††</sup> Yoann Meriguet,<sup>4,‡‡</sup> Irene Donato,<sup>1,2,**</sup> Jordane Preto,<sup>6,¶¶</sup> Luca Varani,<sup>4,***</sup>
James Sturgis,<sup>1,7,†††</sup> and Marco Pettini<sup>1,2,‡‡‡</sup>

<sup>1</sup>Aix-Marseille Univ, Université de Toulon, CNRS, CPT, 13288 Marseille, France

<sup>2</sup>CNRS Centre de Physique Théorique UMR7332, 13288 Marseille, France

<sup>3</sup>Centre d'Immunologie de Marseille-Luminy, 13288 Marseille, France

<sup>4</sup>IES, University of Montpellier-CNRS, UMR 5214, 34000 Montpellier, France

<sup>5</sup>Department of Physics, Sapienza University of Rome, 00185 Rome, Italy

<sup>6</sup>Department of Physics and Department of Oncology, 3-336, Cross Cancer Institute, Edmonton, AB, T6G 1Z2, Canada

<sup>7</sup>Laboratoire d'Ingénierie des Systèmes Macromoléculaires UMR7255, 13402 Marseille, France

<figure><img src="image_2.png" /></figure>

(Received 1 November 2017; revised manuscript received 25 June 2018; published 10 September 2018)

We describe the activation of out-of-equilibrium collective oscillations of a macromolecule as a classical phonon condensation phenomenon. If a macromolecule is modeled as an open system—that is, it is subjected to an external energy supply and is in contact with a thermal bath to dissipate the excess energy— the internal nonlinear couplings among the normal modes make the system undergo a nonequilibrium phase transition when the energy input rate exceeds a threshold value. This transition takes place between a state where the energy is incoherently distributed among the normal modes and a state where the input energy is channeled into the lowest-frequency mode entailing a coherent oscillation of the entire molecule. The model put forward in the present work is derived as the classical counterpart of a quantum model proposed a long time ago by Fröhlich in an attempt to explain the huge speed of enzymatic reactions. We show that such a phenomenon is actually possible. Two different and complementary THz near-field spectroscopic techniques—a plasmonic rectenna and a microwire near-field probe—have been used in two different labs to eliminate artifacts. By considering an aqueous solution of a model protein, the bovine serum albumin, we find that this protein displays a remarkable absorption feature around 0.314 THz, when driven in a stationary out-of-thermal equilibrium state by means of optical pumping. The experimental outcomes are in very good qualitative agreement with the theory developed in the first part of the paper and in excellent quantitative agreement with the theoretical result, allowing us to identify the observed spectral feature with a collective oscillation of the entire molecule.

DOI: 10.1103/PhysRevX.8.031061

Subject Areas: Biological Physics, Nonlinear Dynamics, Soft Matter

## I. INTRODUCTION

Recent progress in terahertz technology has enabled us to look at biological systems with terahertz radiation, that is, in an energy domain (a few meV) of the order of the activation energy of many biological processes. These include the excitation of collective modes of vibration of biomolecules for which experimental evidence has originally been provided at thermal equilibrium by means of Raman spectroscopy [1] and, recently, by several terahertz spectroscopic studies mainly carried out using dry or low-hydrated powders because of the very strong absorption of water [2–5]. More recent studies also addressed solvated proteins [6,7]. Collective oscillations of biomolecules, possibly driven by metabolic activity in living matter by bringing about large oscillating dipole moments, could activate resonant (thus selective), intermolecular, attractive,

<footnote>* i.nardecchia@gmail.com</footnote>

<footnote>† jeremie.torres@umontpellier.fr</footnote>

<footnote>‡ mathias.lechelon@gmail.com</footnote>

<footnote>§ valeria.giliberti@roma1.infn.it</footnote>

<footnote>¶ michele.ortolani@roma1.infn.it</footnote>

<footnote>** philippe.nouvel@umontpellier.fr</footnote>

<footnote>†† gor16matteo@gmail.com</footnote>

<footnote>‡‡ yoann.meriguet@umontpellier.fr</footnote>

<footnote>§§ irene.irened@gmail.com</footnote>

<footnote>¶¶ jordane.preto@gmail.com</footnote>

<footnote>††† luca.varani@umontpellier.fr</footnote>

<footnote>†††† sturgis@imm.cnrs.fr</footnote>

<footnote>‡‡‡ pettini@cpt.univ-mrs.fr; corresponding author</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
electrodynamic forces acting at a long distance. However, we have shown [8] that this is not possible at thermal equilibrium, where all the mentioned experiments have been hitherto performed. Thus, understanding whether long-range intermolecular electrodynamic forces can really be activated requires that we figure out if out-of-thermal-equilibrium collective oscillations of biomolecules can be realized. The possibility of activating these forces could help explain the fast encounters of the cognate partners of biochemical reactions in living matter, even in the context of a low molecular concentration of one of the reactants, as might be the case of some ligand-receptor reaction or of the transcription machinery of DNA. From a theoretical viewpoint, collective oscillations of macromolecules have been described by Fröhlich [9–12] as a Bose-like condensation of their normal vibrational modes.

However, Fröhlich’s theory has been criticized and marginalized for several reasons. One reason is that the original theory is oversimplified and lacks many important explanations to qualify as relevant to understand basic mechanisms of life. Recently, remarkable progress has been made in the direction of assessing the possible relevance of Fröhlich condensates in a revisitation of the theory [13], which has led to the identification of three different kinds of Fröhlich condensates: weak, strong, and coherent. The authors of Ref. [13] explain that weak condensates may have profound effects on chemical and enzyme kinetics, while coherent condensates involve very large amounts of energy, are very fragile, and are not produced by the Wu-Austin Hamiltonian from which Fröhlich’s rate equations can be derived. The work of Ref. [13] is still theoretical, and it is developed in a quantum framework.

From the experimental side, very interesting observations of phonon condensates of quantum nature at thermal equilibrium have recently been reported in the literature. Remarkably, optical phonon condensates have been observed in hetero-structures at room temperature [14], arising from phase and frequency dynamical synchronization of phonons due to suitable scattering processes around defects. Another intriguing result concerns polariton Bose condensation at thermal equilibrium in a GaAs high-Q microcavity in a temperature range of 10–25 K [15], still a quite-high temperature with respect to Bose-Einstein condensation of alkali atoms. These results concern optical phonons or polaritons weakly interacting with phonons, respectively, in condensed matter systems. But what about low-frequency vibrations of biomolecules at room temperature? We may wonder whether a complementary step forward can be made by going beyond the quantum treatment of Fröhlich condensates, also after the revisitation of Ref. [13].

In fact, the frequency of collective oscillations of biomacromolecules is expected in the sub-THz domain, around $10^{11}$ Hz. Hence, at room temperature, it is $k_B T/\hbar\omega \gg 1$ (where $k_B$ and $\hbar$ are the Boltzmann and Planck constants, respectively), and the average number of phonons estimated at this frequency with the Bose-Einstein

statistics ($\langle n \rangle \approx 62.01$) is the same as that given by the Boltzmann statistics ($\langle n \rangle \approx 62.51$) within a 1% accuracy; as a consequence, a classical description seems reasonable. Moreover, when the system is brought out of thermal equilibrium by external energy injection, the effective temperature of each normal mode is higher, thus increasing the occupation number of each mode and making the classical approximation even more reasonable. Furthermore, for example, the vibrational properties of proteins are very well described by molecular dynamics simulations performed in a classical context. Therefore, the first important question addressed in the present work is whether a phonon condensation phenomenology can also be retrieved in a classical framework and out of thermal equilibrium. To this aim, by resorting to a dequantization method, we have worked out a classical version of the original Fröhlich model, finding that—remarkably—in a classical context too, Bose-like phonon condensation is possible [16]. This possibility requires us to consider a biomolecule as an open system—that is, far from thermal equilibrium with its environment—through which energy flows under the simultaneous actions of an external energy supply and of dissipation due to radiative, dielectric, and viscous energy losses. We find that the classical Bose-like condensation in the lowest vibrational mode occurs when the energy input rate exceeds some threshold value. Then, this *a priori* nonobvious result motivates an experimental effort to find if the theoretically predicted phenomenon is actually possible in the physical world. Two independent experiments, in geographically distinct laboratories, have confirmed the existence of out-of-equilibrium collective oscillations for a model protein. This is a proof of concept of which the most significant implication is that, in compliance with a theoretical prediction [8], a crucial prerequisite is fulfilled to excite intermolecular long-range electrodynamic interactions. In turn, as already mentioned above, these interactions could affect the biomolecular dynamics by contributing to drive the high efficiency and rapidity of mutual encounters of the partners of biochemical reactions in living matter, encounters that do not always appear to be the result of Brownian diffusion only.

The derivation of the mentioned classical model— together with the numerical evidence of phonon condensation—is reported in Sec. II. The results of two complementary experiments are reported and discussed in Sec. III. Materials and methods are described in Sec. IV. In Sec. V, we draw some conclusions.

## II. CLASSICAL OUT-OF-EQUILIBRIUM PHONON CONDENSATION

Some decades ago, the study of open systems far from thermodynamic equilibrium showed, under suitable conditions, the emergence of self-organization. Striking similarities were observed among very different physical systems, which have in common the fact that they are
all composed of many nonlinearly interacting subsystems. When a control parameter, typically the energy input rate, exceeds a critical value (that is, when the energy gain exceeds the energy losses), then the subsystems act cooperatively to self-organize in what is commonly referred to as a nonequilibrium phase transition [17,18]. This phenomenon is a generic consequence of the presence of the mentioned basic ingredients, that is, nonlinearly interacting subsystems, dissipation to a thermal bath, and external energy supply. Therefore, even idealized models of real systems are capable of correctly catching the occurrence of these collective behaviors, at least qualitatively. This fascinating topic was pioneered in the late 1960s by Fröhlich with the model mentioned in the Introduction. However, Fröhlich rate equations for the vibrational mode amplitudes, as occupation numbers in Fock's representation, of a generic biomolecule were originally put forward heuristically: The original formulation was lacking a microscopic model. A microscopic Hamiltonian—from which the Fröhlich rate equations can be derived—was later given by Wu and Austin [19–22] to model the dynamics of the normal modes of a macromolecule, of the thermal bath surrounding it, and of the external energy pump. The Wu-Austin Hamiltonian is formulated in the second quantization formalism.

In what follows, we indicate with $\hat{a}_{\omega_i}^\dagger$, $\hat{a}_{\omega_i}$ the quantum creation or annihilation operators for the vibrational normal modes of the main system (i.e., a biomolecule) with frequency $\omega_i \in \mathcal{I}_{\text{sys}}$.

Such a system is put in contact with a thermal bath that represents the degrees of freedom of the environment surrounding the protein and, possibly, other normal modes of the protein that can be considered at thermal equilibrium with the surrounding environment. The thermal bath is characterized by a temperature $T_B$, and it is represented by a collection of harmonic oscillators with characteristic frequencies $\Omega_j \in \mathcal{I}_{\text{bth}}$ whose annihilation or creation operators are $\hat{b}_{\Omega_j}$ and $\hat{b}_{\Omega_k}^\dagger$, respectively.

In order to put the system representing normal modes of a biomolecule out of thermal equilibrium, an external source of energy is necessary: Such an external source is represented as another thermal bath at a temperature $T_S \gg T_B$. Also, in this case, the corresponding thermal bath is described by a collection of harmonic oscillators with frequencies $\Omega_k' \in \mathcal{I}_{\text{src}}$, the quantum annihilation and creation operators of which are $\hat{c}_{\Omega_k'}$ and $\hat{c}_{\Omega_k}^\dagger$. These three sets of harmonic oscillators can be regarded as three subsystems of a larger isolated system $S$ (we coherently indicate with $\mathcal{I}_S$ the set of all the normal modes of the system) whose quantum dynamics is described by the Hamiltonian

$$ \hat{H}_{\text{Tot}} = \hat{H}_0 + \hat{H}_{\text{int}} \quad (1) $$

where $\hat{H}_0$ is the free Hamiltonian of the three sets of harmonic oscillators representing the molecular normal modes and the two heat baths

$$ \hat{H}_0 = \sum_{\omega_i \in \mathcal{I}_{\text{sys}}} \hbar \omega_i \hat{a}_{\omega_i}^{\dagger} \hat{a}_{\omega_i} + \sum_{\Omega_j \in \mathcal{I}_{\text{bth}}} \hbar \Omega_j \hat{b}_{\Omega_j}^{\dagger} \hat{b}_{\Omega_j} \\ + \sum_{\Omega_k' \in \mathcal{I}_{\text{src}}} \hbar \Omega_k' \hat{c}_{\Omega_k'}^{\dagger} \hat{c}_{\Omega_k'} \quad (2) $$

The interactions among normal modes are described by $\hat{H}_{\text{int}}$; in the original formulation by Wu and Austin, the interaction term has the form

$$ \begin{aligned} \hat{H}_{\text{intWA}} &= \hat{H}_{\text{sys-bth}} + \hat{H}_{\text{src-sys}} + \hat{H}_{\text{sys-bth-sys}} \\ &= \sum_{\omega_i \in \mathcal{I}_{\text{sys}}, \Omega_j \in \mathcal{I}_{\text{bth}}} \eta_{\omega_i, \Omega_j} \hat{a}_{\omega_i}^{\dagger} \hat{b}_{\Omega_j} \\ &\quad + \sum_{\omega_i \in \mathcal{I}_{\text{sys}}, \Omega_k' \in \mathcal{I}_{\text{src}}} \xi_{\omega_i, \Omega_k'} \hat{a}_{\omega_i}^{\dagger} \hat{c}_{\Omega_k'} \\ &\quad + \sum_{\omega_{A_i}, \omega_{A_j} \in \mathcal{I}_{\text{sys}}, \Omega_k \in \mathcal{I}_{\text{bth}}} \chi_{\omega_i, \omega_j, \Omega_k} \hat{a}_{\omega_i}^{\dagger} \hat{a}_{\omega_j} \hat{b}_{\Omega_k}^{\dagger} + \text{H.c.}, \end{aligned} \quad (3) $$

where $\eta_{\omega_i, \Omega_j}$, $\xi_{\omega_i, \Omega_k'}$, $\chi_{\omega_i, \omega_j, \Omega_k} \in \mathbb{C}$ are the coupling constants describing the linear interactions among the thermal bath modes and the biomolecule modes, the linear interactions between the external source and the biomolecule, and the mode-mode interactions among the biomolecule normal modes mediated by the thermal bath, respectively.

From these terms, it is possible to derive the Fröhlich rate equations by resorting to time-dependent perturbation theory; details are given in Refs. [19,20] and in the reference book in Ref. [23].

However, the mode coupling term corresponds to a potential energy unbounded from below, and consequently, this would give rise to dynamical instability of the system and, in the quantum context, to the absence of a finite energy ground state. This problem led to strong criticism against the Wu-Austin Hamiltonian model and also against the ensemble of Fröhlich condensation theory [24]. This problem can be easily fixed by adding a term with a quartic dependence on the creation and annihilation operators of the form

$$ \begin{aligned} \hat{H}_{\text{intQ}} = & \sum_{\omega_i, \omega_j, \omega_k, \omega_l \in \mathcal{I}_{\text{sys}}} [\kappa_{(1)} \omega_i \omega_j \omega_k \omega_l \hat{a}_{\omega_i}^{\dagger} \hat{a}_{\omega_j}^{\dagger} \hat{a}_{\omega_k} \hat{a}_{\omega_l} \\ & + \kappa_{(2)} \omega_i \omega_j \omega_k \omega_l \hat{a}_{\omega_i}^{\dagger} \hat{a}_{\omega_j}^{\dagger} \hat{a}_{\omega_k}^{\dagger} \hat{a}_{\omega_l} \\ & + \kappa_{(3)} \omega_i \omega_j \omega_k \omega_l \hat{a}_{\omega_i}^{\dagger} \hat{a}_{\omega_j}^{\dagger} \hat{a}_{\omega_k}^{\dagger} \hat{a}_{\omega_l}^{\dagger}] + \text{H.c.}, \end{aligned} \quad (4) $$

so the lower bound to the ground energy level does not go to $-\infty$ for large values of $N_{\omega_i}$. This quartic interaction stands for an anharmonic interaction among the normal modes of the biomolecule, a broadly studied topic of relevance to energy transport in biomolecules like proteins.
A. Dequantization of Wu-Austin Hamiltonian by time-dependent variational principle

In this section, a classical Hamiltonian system with its canonical equations is associated with the quantum system described by the quantum Hamiltonian of Eqs. (1), (3), and (4). This result is obtained by applying the time-dependent variational principle (TDVP) technique [25,26]. The same technique, which is deeper from both the mathematical and conceptual points of view, has been proposed as a “dequantization” technique in Ref. [27] as a kind of inverse procedure with respect to the geometrical quantization. This consists in the evaluation of the time-dependent operator action on coherent states of quantum harmonic oscillators. The scalar parameters describing the coherent states become generalized coordinates of a classical dynamical system whose equations of motion can be derived from a variational principle.

In more detail, one begins with the ansatz that the wave function depends on N parameters {$x_i$}$_{i=1,...,N}$,

$$|\psi\rangle = |\psi(x_1, \dots, x_N)\rangle, \quad (5)$$

where the parameters $x_i = x_i(t)$ are, in general, functions of time. For a quantum system with Hamiltonian $\hat{H}_{\text{Tot}}$, the equations of motion of $x_i$ can be derived using the following variational principle (equivalent to the least action principle):

$$\delta S = 0 \quad \text{with} \quad S = \int_0^t L(\psi, \bar{\psi}) dt', \quad (6)$$

where $L(\psi, \bar{\psi})$ is the Lagrangian associated with the system

$$L(\psi, \bar{\psi}) = \frac{i\hbar}{2} \left\langle \psi | \dot{\psi} \rangle - \left\langle \dot{\psi} | \psi \right\rangle - \frac{\langle \psi | \hat{H}_{\text{Tot}} | \psi \rangle}{\langle \psi | \psi \rangle} . \quad (7)$$

The equations of motions derived from Eq. (6) can be worked out in the framework of classical Hamiltonian dynamics.

The classical Hamiltonian is associated with the quantum one by simply taking the expectation value of the Hamiltonian operator $\hat{H}_{\text{Tot}}$ over the state $|\psi(x_1, \dots, x_N)\rangle$, that is,

$$H_{\text{Tot}} = \langle \psi(x_1, \dots, x_N) | \hat{H}_{\text{Tot}} | \psi(x_1, \dots, x_N) \rangle. \quad (8)$$

The Poisson brackets {$\cdot, \cdot$} depend only on the chosen parametrization for the wave function. Starting from the variables

$$w_i = i\hbar \left\langle \psi \left| \frac{\partial \psi}{\partial x_i} \right| \psi \right\rangle = -i\hbar \left\langle \frac{\partial \psi}{\partial x_i} \right| \psi \right\rangle, \quad (9)$$

one defines the antisymmetric tensor

$$W_{ij} = \frac{\partial w_j}{\partial x_i} - \frac{\partial w_i}{\partial x_j}, \quad (10)$$

so the equations of motion are implicitly given by

$$\sum_{j=1}^{N} W_{ij} \dot{x}_j = \frac{\partial H_{\text{Tot}}}{\partial x_i}. \quad (11)$$

If the condition det $W_{ij} \neq 0$ holds, then the matrix $W_{ij} = (W^{-1})_{ij}$ defines the Poisson brackets for the classical Hamiltonian system

$$\{f, g\} = \sum_{i,j}^{N} \frac{\partial f}{\partial x_i} \frac{\partial g}{\partial x_j}. \quad (12)$$

This formalism can be applied to the quantum system described by the quantum Hamiltonian of Eq. (1) to associate it with a classical Hamiltonian system. The choice of the parametrization for the wave function is quite arbitrary, and the TDVP, as any other variational principle, restricts the dynamics to a certain region of the Hilbert space. Since the Hamiltonian is expressed in terms of creation and annihilation operators of the quantum harmonic oscillators describing the system, the wave function is chosen as a product of the corresponding coherent states. In particular,

$$|\Psi(t)\rangle = \prod_{\omega_i \in \mathcal{I}_{\text{sys}}, \Omega_j \in \mathcal{I}_{\text{bth}}, \Omega'_k \in \mathcal{I}_{\text{src}}} |z_{\omega_i}(t)\rangle_{\text{sys}} \otimes |z_{\Omega_j}(t)\rangle_{\text{bth}} \otimes |z_{\Omega'_k}(t)\rangle_{\text{src}}, \quad (13)$$

where $|z_{\Omega_j}(t)\rangle_{\text{sys}}$, $|z_{\Omega_j}(t)\rangle_{\text{bth}}$, $|z_{\Omega'_k}(t)\rangle_{\text{src}}$ are normalized coherent states for the normal modes of the main system, the thermal bath, and the external source, respectively; their general form is given by

$$\begin{aligned} |z\rangle &= \exp \left[ -\frac{|z|^2}{2} \right] \sum_{k=0}^{+\infty} \frac{z^k}{\sqrt{k!}} |k\rangle \\ &= \exp \left[ -\frac{|z|^2}{2} \right] \sum_{k=0}^{+\infty} \frac{(z\hat{a}^\dagger)^k}{k!} |0\rangle, \end{aligned} \quad \text{where } z = z(t) \in \mathbb{C}. \quad (14)$$

From the definition of coherent states in Eq. (14), it follows that the expectation value for the occupation number n is given by the squared norm of z,

$$\langle n(t) \rangle = \langle z(t) | \hat{a}^\dagger \hat{a} | (t) \rangle = |z(t)|^2, \quad (15)$$

so, as we are interested in writing rate equations for these quantities, we parametrize the wave function $\Psi(t)$ with the set of real parameters {$n_i, \theta_i$}$_{i \in \mathcal{I}_S}$ such that
$$z_i = n_i^{1/2} \exp[-i\theta_i] \Rightarrow z_i^* = n_i^{1/2} \exp[i\theta_i], \quad n_i = |z_i|^2, \quad i \in \mathcal{I}_S. \qquad (16)$$

Using Eq. (9), it is possible to derive the Poisson brackets associated with the variables $\{(n_i, \theta_i)\}_{i \in \mathcal{I}_S}$:

$$ \begin{aligned} w_{n_i} &= i\hbar \left\langle \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} \left| \frac{\partial}{\partial n_i} \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} \right\rangle \right. \\ &= \frac{i\hbar}{\partial n_i / \partial |z_i|} \left\langle \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} \left| \frac{\partial}{\partial |z_i|} \otimes_{j \in \mathcal{I}_S} \left( \exp \left[ -\frac{|z_j|^2}{2} \right] \sum_{k=0}^{+\infty} \frac{|z_j|^k \exp[-ik\theta_j] (\hat{a}_j^\dagger)^k}{k!} |0\rangle \right) \right. \\ &= \frac{i\hbar}{2|z_i|} (-|z_i| \langle \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} | \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} \rangle + \exp[-i\theta_i] \langle \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} | \hat{a}_j^\dagger | \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} \rangle) \\ &= \frac{i\hbar}{2|z_i|} (-|z_i| + |z_i| \exp[i\theta_i] \exp[-i\theta_i]) = 0, \end{aligned} \qquad (17) $$

$$ \begin{aligned} w_{\theta_i} &= i\hbar \left\langle \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} \left| \frac{\partial}{\partial \theta_i} \right\rangle \left\langle \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} \right| \\ &= i\hbar \left\langle \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} \left| \frac{\partial}{\partial \theta_i} \otimes_{j \in \mathcal{I}_S} \left( \exp \left[ -\frac{|z_j|^2}{2} \right] \sum_{k=0}^{+\infty} \frac{|z_j|^k \exp[-ik\theta_j] (\hat{a}_j^\dagger)^k}{k!} |0\rangle \right) \right\rangle \\ &= i(-i)\hbar |z_i| \exp[-i\theta_i] \langle \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} | \hat{a}_i^\dagger | \{(n_j, \theta_j)\}_{j \in \mathcal{I}_S} \rangle = \hbar |z_i|^2 = \hbar n_i. \end{aligned} \qquad (18) $$

Consequently, using the definition (10), the entries of the W matrix are

$$W_{\theta_i \theta_k} = W_{n_i n_k} = 0, \qquad (19)$$

$$W_{n_i \theta_k} = -W_{\theta_k n_i} = \frac{\partial w_{n_k}}{\partial \theta_i} - \frac{\partial w_{\theta_i}}{\partial n_k} = \hbar \delta_{i,k}, \qquad (20)$$

and its inverse has the form

$$W_{\theta_i \theta_k} = W_{n_i n_k} = 0, \qquad (21)$$

$$W_{\theta_i n_k} = -W_{n_k \theta_i} = \frac{\delta_{i,k}}{\hbar}. \qquad (22)$$

Thus, it follows that the variables $J_\omega = \hbar n_\omega$ and $\theta_\omega$ are canonically conjugated variables. The classical Hamiltonian $H = H_0 + H_{\text{intWA}} + H_{\text{intQ}}$ for the variables $\{(\theta_\omega, J_\omega)\}_{\omega \in \mathcal{I}_S}$ is given by a free classical part

$$H_0 = \langle \Psi(\theta_\omega, J_\omega) | \hat{H}_0 | \Psi(\theta_\omega, J_\omega) \rangle = \sum_{\omega_i \in \mathcal{I}_{\text{sys}}} \omega_i J_{\omega_i} + \sum_{\Omega_k \in \mathcal{I}_{\text{bth}}} \Omega_k J_{\Omega_k} + \sum_{\Omega'_p \in \mathcal{I}_{\text{src}}} \Omega'_p J_{\Omega'_p}, \qquad (23)$$

by a semiclassical Wu and Austin interaction part

$$ \begin{aligned} H_{\text{intWA}} &= \langle \Psi(\theta_\omega, J_\omega) | \hat{H}_{\text{intWA}} | \Psi(\theta_\omega, J_\omega) \rangle \\ &= \sum_{\omega_i \in \mathcal{I}_{\text{sys}}} \sum_{\Omega_k \in \mathcal{I}_{\text{bth}}} \frac{|\eta_{\omega_i \Omega_k}|}{\hbar} J_{\omega_i}^{1/2} J_{\Omega_k}^{1/2} \cos(\theta_{\omega_i} - \theta_{\Omega_k} + \theta_{\eta_{\omega_i \Omega_k}}) \\ &\quad + \sum_{\omega_i \in \mathcal{I}_{\text{sys}}} \sum_{\Omega'_p \in \mathcal{I}_{\text{src}}} \frac{|\xi_{\omega_i \Omega'_p}|}{\hbar} J_{\omega_i}^{1/2} J_{\Omega'_p}^{1/2} \cos(\theta_{\omega_i} - \theta_{\Omega'_p} + \theta_{\xi_{\omega_i \Omega'_p}}) \\ &\quad + \sum_{\omega_i, \omega_j \in \mathcal{I}_{\text{sys}}} \sum_{\Omega_k \in \mathcal{I}_{\text{bth}}} \frac{|\chi_{\omega_i \omega_j \Omega_k}|}{\hbar^{3/2}} 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} + \theta_{\chi_{\omega_i \omega_j \Omega_k}}), \end{aligned} \qquad (24) $$
and by the quartic term

$$ H_{\text{intQ}} = \langle \Psi(\theta_{\omega}, J_{\omega}) | \hat{H}_{\text{intQ}} | \Psi(\theta_{\omega}, J_{\omega}) \rangle \\ = \sum_{\omega_i, \omega_j, \omega_k, \omega_l \in \mathcal{I}_{\text{sys}}} J_{\omega_i}^{1/2} J_{\omega_j}^{1/2} J_{\omega_k}^{1/2} J_{\omega_l}^{1/2} \left[ \frac{|\mathcal{K}(1)_{\omega_i \omega_j \omega_k \omega_l}|}{\hbar^2} \cos (\theta_{\omega_i} + \theta_{\omega_j} - \theta_{\omega_k} - \theta_{\omega_l} + \theta_{\mathcal{K}(1)\omega_i \omega_j \omega_k \omega_l}) \right. \\ \left. + \frac{|\mathcal{K}(2)_{\omega_i \omega_j \omega_k \omega_l}|}{\hbar^2} \cos (\theta_{\omega_i} + \theta_{\omega_j} + \theta_{\omega_k} - \theta_{\omega_l} + \theta_{\mathcal{K}(2)\omega_i \omega_j \omega_k \omega_l}) \right. \\ \left. + \frac{|\mathcal{K}(3)_{\omega_i \omega_j \omega_k \omega_l}|}{\hbar^2} \cos (\theta_{\omega_i} + \theta_{\omega_j} + \theta_{\omega_k} + \theta_{\omega_l} + \theta_{\mathcal{K}(3)\omega_i \omega_j \omega_k \omega_l}) \right], \quad (25) $$

where each complex coupling constant is given in the polar representation. In what follows, the coupling constants are considered real and rescaled such that

$$ \frac{|\eta_{\omega_i \Omega_k}|}{\hbar} \rightarrow \eta_{\omega_i \Omega_k}, \quad \frac{|\xi_{\omega_i \Omega'_p}|}{\hbar} \rightarrow \xi_{\omega_i \Omega'_p}, \quad \frac{|\chi_{\omega_i \Omega_j \Omega_k}|}{\hbar^{3/2}} \rightarrow \chi_{\omega_i \Omega_j \Omega_k}, \quad \frac{|\mathcal{K}(1,2,3)_{\omega_i \omega_j \omega_k \omega_l}|}{\hbar^2} \rightarrow \mathcal{K}(1,2,3)_{\omega_i \omega_j \omega_k \omega_l}. \quad (26) $$

With these choices, the total Hamiltonian of the system reads

$$ H_{\text{Tot}}(\{(J_{\omega_i}, \theta_{\omega_i})\}_{\omega_i \in \mathcal{I}_S}) = \sum_{\omega_i \in \mathcal{I}_{\text{sys}}} \omega_i J_i + \sum_{\Omega_j \in \mathcal{I}_{\text{bth}}} \Omega_j J_{\Omega_j} + \sum_{\Omega'_k \in \mathcal{I}_{\text{src}}} \Omega'_k J_{\Omega'_k} \\ + \sum_{\omega_i \in \mathcal{I}_{\text{sys}}} \sum_{\Omega_j \in \mathcal{I}_{\text{bth}}} \eta_{\omega_i \Omega_j} J_{\omega_i}^{1/2} J_{\Omega_j}^{1/2} \cos (\theta_{\omega_i} - \theta_{\Omega_j}) + \sum_{\omega_i \in \mathcal{I}_{\text{sys}}} \sum_{\Omega'_k \in \mathcal{I}_{\text{src}}} \xi_{\omega_i \Omega'_k} J_{\omega_i}^{1/2} J_{\Omega'_k}^{1/2} \cos (\theta_{\omega_i} - \theta_{\Omega'_k}) \\ + \sum_{\omega_i, \omega_j \in \mathcal{I}_{\text{sys}}} \sum_{\Omega_k \in \mathcal{I}_{\text{bth}}} \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}) \\ + \sum_{\omega_i, \omega_j, \omega_k, \omega_l \in \mathcal{I}_{\text{sys}}} J_{\omega_i}^{1/2} J_{\omega_j}^{1/2} J_{\omega_k}^{1/2} J_{\omega_l}^{1/2} [\mathcal{K}(1)_{\omega_i \omega_j \omega_k \omega_l} \cos (\theta_{\omega_i} + \theta_{\omega_j} - \theta_{\omega_k} - \theta_{\omega_l}) \\ + \mathcal{K}(2)_{\omega_i \omega_j \omega_k \omega_l} \cos (\theta_{\omega_i} + \theta_{\omega_j} + \theta_{\omega_k} - \theta_{\omega_l}) + \mathcal{K}(3)_{\omega_i \omega_j \omega_k \omega_l} \cos (\theta_{\omega_i} + \theta_{\omega_j} + \theta_{\omega_k} + \theta_{\omega_l})]. \quad (27) $$

In order to derive Fröhlich-like rate equations, the dynamics of the action variables $J_{\omega_i}$ of the system has to be studied. We could choose to investigate the dynamics of the system by letting observable quantities evolve in time (according to Hamilton's equations of motion) and then performing time averaging, and averaging on different initial conditions compatible with the assumption that the two subsystems $\mathcal{I}_{\text{bth}}$ and $\mathcal{I}_{\text{src}}$ are two thermal baths with different temperatures. However, this method has some disadvantages: The integration should be done numerically because of the presence of nonlinear interaction terms in the Hamiltonian, and a very large number of degrees of freedom would be necessary to adequately simulate the dynamics of the thermal baths. Moreover, long integration times would be necessary to attain the convergence of time averages of the observables, and, finally, this computational effort would not be worth it, as it would provide redundant information on the dynamics. Another way to derive rate equations for the $J_{\omega_i}$ consists in

a statistical approach: The relevant dynamical variable to consider is the phase-space distribution function $\rho(\{(J_{\omega_i}, \theta_{\omega_i})\}_{\omega_i \in \mathcal{I}_S}, t)$ so that the rate equations are written for the statistical averages of actions variables

$$ \langle J_{\omega_i}(t) \rangle = \int J_{\omega_i} \rho(\{(J_{\omega_i}, \theta_{\omega_i})\}_{\omega_i \in \mathcal{I}_S}, t) \prod_{\omega \in \mathcal{I}_S} dJ_{\omega} d\theta_{\omega}. \quad (28) $$

In the following section, classical Fröhlich-like rate equations are derived by resorting to the time evolution of the distribution function $\rho(\{(J_{\omega_i}, \theta_{\omega_i})\}_{\omega_i \in \mathcal{I}_S}, t)$ satisfying the Liouville equation.

## B. Derivation of classical rate equations using the Koopman–von Neumann formalism

Let $\rho(\{(J_{\omega_i}, \theta_{\omega_i})\}_{\omega_i \in \mathcal{I}_S}, t)$ be a probability density function for the whole system described by the Hamiltonian in Eq. (27); according to the Liouville theorem, the evolution of $\rho$ associated with this Hamiltonian is given by
$$ \frac{\partial \rho}{\partial t} = \{H, \rho\} = -i\mathcal{L}_H(\rho), \qquad (29) $$

where $t^2 = -1$ and $\{\cdot, \cdot\}$ are the canonical Poisson brackets

$$ \begin{aligned} \{f,g\} = & \sum_{\omega_i \in \mathcal{I}_{\text{sys}}} \left( \frac{\partial f}{\partial J_{\omega_i}} \frac{\partial g}{\partial \theta_{\omega_i}} - \frac{\partial g}{\partial J_{\omega_i}} \frac{\partial f}{\partial \theta_{\omega_i}} \right) \\ & + \sum_{\Omega_j \in \mathcal{I}_{\text{bth}}} \left( \frac{\partial f}{\partial J_{\Omega_j}} \frac{\partial g}{\partial \theta_{\Omega_j}} - \frac{\partial g}{\partial J_{\Omega_j}} \frac{\partial f}{\partial \theta_{\Omega_j}} \right) \\ & + \sum_{\Omega'_k \in \mathcal{I}_{\text{src}}} \left( \frac{\partial f}{\partial J_{\Omega'_k}} \frac{\partial g}{\partial \theta_{\Omega'_k}} - \frac{\partial g}{\partial J_{\Omega'_k}} \frac{\partial f}{\partial \theta_{\Omega'_k}} \right), \end{aligned} \qquad (30) $$

and $\mathcal{L}_H(\cdot) = i\{H, \cdot\}$ is the Liouville operator acting on functions defined on the phase space of the system. An interesting method to study and solve Liouville equations relies on the Koopman–von Neumann (KvN) formalism developed in the 1930s: A formal analogy among Liouville and Schrödinger equations is established so that classical mechanics can also be formulated in the framework of a Hilbert space of square integrable functions. In our case, the Hilbert space of complex square integrable functions in phase space is $L^2(\Lambda_{\{(J_\omega, \theta_\omega)\}_{\omega \in \mathcal{I}_S}})$, with the inner product defined by

$$ \begin{aligned} \langle f | g \rangle &= \int_{\Lambda_{\{(J_\omega, \theta_\omega)\}_{\omega \in \mathcal{I}_S}}} f^* g \prod_{\omega \in \mathcal{I}_S} dJ_\omega d\theta_\omega \\ &= \prod_{\omega \in \mathcal{I}_S} \int_0^{2\pi} d\theta_\omega \int_0^{+\infty} dJ_\omega f^* g, \end{aligned} \qquad (31) $$

with $f, g \in L^2(\Lambda_{\{(J_\omega, \theta_\omega)\}_{\omega \in \mathcal{I}_S}})$. On this space, we can define the action of the Liouville operator

$$ \hat{\mathcal{L}}_H |f\rangle = \mathcal{L}_H(f) \qquad (32) $$

and consider the domain $\mathcal{D}_{\hat{\mathcal{L}}_H} \subseteq L^2(\Lambda_{\{(J_\omega, \theta_\omega)\}_{\omega \in \mathcal{I}_S}})$ where the Liouville operator is self-adjoint, namely, $\hat{\mathcal{L}}_H^\dagger = \hat{\mathcal{L}}_H$ and $\mathcal{D}_{\hat{\mathcal{L}}_H} = \mathcal{D}_{\hat{\mathcal{L}}_H^\dagger}$. Let $\psi(\{(J_\omega, \theta_\omega)\}_{\omega \in \mathcal{I}_S}; t) \in \mathcal{D}_{\hat{\mathcal{L}}_H}$ be a normalized time-dependent function [28] such that

$$ i \frac{\partial \psi}{\partial t}(\mathbf{J}, \boldsymbol{\theta}; t) = \hat{\mathcal{L}}_H \psi(\mathbf{J}, \boldsymbol{\theta}; t); \qquad (33) $$

then, it can be proven that $\rho = \|\psi\|_{L^2(\Lambda_{\{(J_\omega, \theta_\omega)\}_{\omega \in \mathcal{I}_S}})} = \psi^*\psi$ is a normalized function for which Eq. (29) holds. Moreover, as $\hat{\mathcal{L}}_H$ is a self-adjoint operator, it represents the unitary time evolution of the “wave function” as

$$ \psi(\{(J_\omega, \theta_\omega)\}_{\omega \in \mathcal{I}_S}; t) = \exp[-it\hat{\mathcal{L}}_H]\psi(\{(J_\omega, \theta_\omega)\}_{\omega \in \mathcal{I}_S}; 0) \qquad (34) $$

in analogy with quantum mechanics. With this formalism, the rate equations for the average values of the actions

(which are the counterparts of quantum occupation numbers) associated with the normal modes $\omega_i$ of the main system are given by

$$ \begin{aligned} \frac{d}{dt} \langle J_{\omega_i} \rangle_t &= \frac{d}{dt} \langle \psi(t) | \hat{\mathcal{M}}_{J_{\omega_i}} | \psi(t) \rangle \\ &= i \langle \psi(t) | [\hat{\mathcal{L}}_H, \hat{\mathcal{M}}_{J_{\omega_i}}] | \psi(t) \rangle, \end{aligned} \qquad (35) $$

where $\hat{\mathcal{M}}_{J_{\omega_i}}$ is a multiplicative operator acting on $L^2(\Lambda_{\{(J_\omega, \theta_\omega)\}_{\omega \in \mathcal{I}_S}})$ as follows:

$$ \hat{\mathcal{M}}_{J_{\omega_i}} |\psi\rangle = |J_{\omega_i} \psi\rangle. \qquad (36) $$

The Liouville operator can be decomposed as

$$ \hat{\mathcal{L}}_H = \hat{\mathcal{L}}_{H_0} + \hat{\mathcal{L}}_{H_{\text{int}}}(t), \qquad (37) $$

and since the eigenfunctions of the operator $\hat{\mathcal{L}}_{H_0}$ are known, the action of the operator $\hat{\mathcal{L}}_{H_{\text{int}}}$ can be treated as a time-dependent perturbation, which is adiabatically turned on and off from $t_0 = 0$ to $t = +\infty$. Then, a classical analog of the interaction representation formalism is used. Thus, if $|\psi(t)\rangle_S$ is the wave function in the Schrödinger representation, then, in the interaction representation, $|\psi(t)\rangle_I$ reads

$$ |\psi(t)\rangle_I = \exp[\mathrm{it}\hat{\mathcal{L}}_{H_0}] |\psi(t)\rangle_S, \qquad (38) $$

and, given a generic operator $\hat{A}_S$ in the Schrödinger picture, its expression in the interaction picture $\hat{A}_I$ is

$$ \hat{A}_I(t) = \exp[\mathrm{it}\hat{\mathcal{L}}_{H_0}] \hat{A}_S \exp[-\mathrm{it}\hat{\mathcal{L}}_{H_0}]. \qquad (39) $$

With this formalism, the time evolution of $|\psi(t)\rangle_I$ can be written with the unitary evolution operator $\hat{U}(t; t_0)$ satisfying

$$ |\psi(t)\rangle_I = \hat{U}(t; t_0) |\psi(t_0)\rangle_I, \qquad (40) $$

$$ i \frac{\partial U(t; t_0)}{\partial t} = \hat{\mathcal{L}}_{H_{\text{int}}} (t) \hat{U}(t; t_0), \qquad (41) $$

and the formal solution of Eq. (40) is given by

$$ U(t; t_0) = \mathbb{I} - i \int_{t_0}^{t} \hat{\mathcal{L}}_{H_{\text{int}}} (t') \hat{U}(t'; t_0) dt'. \qquad (42) $$

At first order in $\hat{\mathcal{L}}_{H_{\text{int}}}(t')$, the unitary evolution operator $\hat{U}(t; t_0)$ in the right-hand side of Eq. (42) is substituted by the identity operator, meaning that the state $|\psi(t_0)\rangle_I$, if the perturbation is turned on at $t_0$, at the lowest order can be approximated by $|\psi(t_0)\rangle_I \approx |\psi(t_0)\rangle_S$ and is assumed to be coincident with the Schrödinger picture [i.e., $|\psi(0)\rangle = |\psi_0\rangle$]. Thus,
$$|\psi(t)\rangle_I \approx \left(\hat{\mathbb{I}} - i \int_0^t \hat{\mathcal{L}}_{H_{\text{int}}} (t') \mathrm{d}t'\right) |\psi_0\rangle_I \quad (43)$$

and, as $\hat{\mathcal{L}}_{H_{\text{int}}} (t')$ is self-adjoint, the time evolution for the “bra” has the form

$$\langle\psi(t)|_I \approx \langle\psi_0|_I \left(\hat{\mathbb{I}} + i \int_0^t \hat{\mathcal{L}}_{H_{\text{int}}} (t') \mathrm{d}t'\right). \quad (44)$$

The time derivative of the multiplicative operator $\hat{\mathcal{M}}_{J_{\omega_i}}$ in the interaction picture is derived according to Eq. (39),

$$\begin{aligned} (\hat{\mathcal{M}}_{J_{\omega_i}})_I &= i ([\hat{\mathcal{L}}_H, \hat{\mathcal{M}}_{J_{\omega_i}}])_I \\ &= i \exp [i t \hat{\mathcal{L}}_{H_0}] [\hat{\mathcal{L}}_H, \hat{\mathcal{M}}_{J_{\omega_i}}] \exp [-i t \hat{\mathcal{L}}_{H_0}]. \end{aligned} \quad (45)$$

This means that the average (35) can be entirely rewritten using the interaction picture as

$$\begin{aligned} \frac{d}{dt} \langle J_{\omega_i} \rangle_t &= I \langle \psi(t) | (\hat{\mathcal{M}}_{J_{\omega_i}})_I | \psi(t) \rangle_I \\ &= I \langle \psi(t) | ( [\hat{\mathcal{L}}_H, \hat{\mathcal{M}}_{J_{\omega_i}}] (t) )_I | \psi(t) \rangle_I \\ &\approx \langle \psi_0 | \left( \hat{\mathbb{I}} + \int_{-\infty}^t \hat{\mathcal{L}}_{H_{\text{int}}} (t') \mathrm{d}t' \right) ( [\hat{\mathcal{L}}_H, \hat{\mathcal{M}}_{J_{\omega_i}}] (t) )_I \left( \hat{\mathbb{I}} - \int_0^t \hat{\mathcal{L}}_{H_{\text{int}}} (t') \mathrm{d}t' \right) | \psi_0 \rangle \\ &= \langle \psi_0 | ( [\hat{\mathcal{L}}_H, \hat{\mathcal{M}}_{J_{\omega_i}}] (t) )_I | \psi_0 \rangle + \int_0^t \langle \psi_0 | [ ( [\hat{\mathcal{L}}_{H_{\text{int}}}, \hat{\mathcal{M}}_{J_{\omega_i}}] (t) )_I, ( \hat{\mathcal{L}}_{H_{\text{int}}} (t') )_I ] | \psi_0 \rangle \mathrm{d}t'. \end{aligned} \quad (46)$$

After lengthy computations (see Supplemental Material [29] for details), one finally arrives at the following rate equations for the expectation values of the actions (that is, the amplitudes of vibrational modes), which are worked out as follows:

$$\begin{aligned} \frac{d\langle J_{\omega_i} \rangle}{dt} = & s_{\omega_i} + b_{\omega_i} \left[ \frac{k_B T_B}{\omega_i} - \langle J_{\omega_i} \rangle \right] + \sum_{\substack{\omega_j \in \mathcal{I}_{\text{sys}} \\ \omega_j \neq \omega_i}} c_{\omega_i \omega_j} \left[ \langle (J_{\omega_j}) \rangle - \langle J_{\omega_i} \rangle \right] \\ & + \sum_{\substack{\omega_j, \omega_k, \omega_l \in \mathcal{I}_{\text{sys}} \\ \omega_i + \omega_j - \omega_k - \omega_l = 0}} \frac{16\pi\kappa_{(1)\omega_i\omega_j\omega_k\omega_l}^2}{\delta\omega_{\text{sys}}} \langle J_{\omega_l} \rangle (\langle J_{\omega_j} \rangle \langle J_{\omega_k} \rangle + \langle J_{\omega_i} \rangle \langle J_{\omega_k} \rangle - 2\langle J_{\omega_i} \rangle \langle J_{\omega_j} \rangle) \\ & + \frac{3\pi}{\delta\omega_{\text{sys}}} \left[ \sum_{\substack{\omega_j, \omega_k, \omega_l \in \mathcal{I}_{\text{sys}} \\ \omega_i + \omega_j + \omega_k - \omega_l = 0}} 3\kappa_{(2)\omega_i\omega_j\omega_k\omega_l}^2 \langle J_{\omega_l} \rangle (\langle J_{\omega_j} \rangle \langle J_{\omega_k} \rangle + 2\langle J_{\omega_i} \rangle \langle J_{\omega_k} \rangle - \langle J_{\omega_i} \rangle \langle J_{\omega_j} \rangle) \right. \\ & \left. + \sum_{\substack{\omega_j, \omega_k, \omega_l \in \mathcal{I}_{\text{sys}} \\ \omega_i + \omega_j + \omega_k - \omega_l = 0}} \kappa_{(2)\omega_i\omega_j\omega_k\omega_l}^2 \langle J_{\omega_l} \rangle \langle J_{\omega_k} \rangle (\langle J_{\omega_j} \rangle - 3\langle J_{\omega_i} \rangle) \right]. \end{aligned} \quad (47)$$

### C. Properties of the classical rate equations

The rate equations (47) that have been derived by solving the Liouville equation with the Hamiltonian (27) display a functionally similar structure to the original Fröhlich’s equations [9], apart from the additional quartic terms. As in the original formulation given by Fröhlich, the condensation phenomenon is found by considering the stationary solutions of the rate equations. It is convenient to rewrite Eq. (47) in a non-dimensional form after introducing the following variables:

$$\tau = t\omega_0 \text{ with } \omega_0 = \min_{\omega \in \mathcal{I}_{\text{sys}}} \omega, \quad y_{\omega_i} = \frac{\omega_i \langle J_{\omega_i} \rangle}{k_B T_B},$$

$$\alpha_{\omega_i} = \frac{\omega_i}{\omega_0}, \quad S_{\omega_i} = \alpha_{\omega_i} \frac{s_{\omega_i}}{k_B T_B}, \quad B_{\omega_i} = \frac{b_{\omega_i}}{\omega_0},$$

$$C_{\omega_i \omega} = \frac{c_{\omega_i \omega}}{\omega_0}, \quad \Upsilon_{(1)\omega_i \omega_j \omega_k \omega_l} = \frac{16\pi\kappa_{(1)\omega_i \omega_j \omega_k \omega_l} (k_B T_B)^2}{\delta\omega_{\text{sys}} \omega_0^3},$$

$$\Upsilon_{(2)\omega_i\omega_j\omega_k\omega_l} = \frac{3\pi\kappa_{(2)\omega_i\omega_j\omega_k\omega_l} (k_B T_B)^2}{\delta\omega_{\text{sys}}\omega_0^3}. \quad (48)$$

Thus, Eq. (47) reads
$$
\begin{equation}
\begin{aligned}
\dot{y}_{\omega_i} = \frac{dy_{\omega_i}}{dt} = & S_{\omega_i} + B_{\omega_i}(1 - y_{\omega_i}) + \sum_{\substack{\omega_j \in \mathcal{I}_{\text{sys}} \\ \omega_j \neq \omega_i}} C_{\omega_i \omega_j} \left[ \left( \frac{\alpha_{\omega_i}}{\alpha_{\omega_j}} y_{\omega_j} - y_{\omega_i} \right) + \frac{(\alpha_{\omega_j} - \alpha_{\omega_i})}{\alpha_{\omega_j}} y_{\omega_i} y_{\omega_j} \right] \\
& + \sum_{\substack{\omega_i, \omega_j, \omega_k, \omega_l \in \mathcal{I}_{\text{sys}} \\ \omega_i + \omega_j + \omega_k - \omega_l = 0}} \Upsilon_{(1)} \frac{y_{\omega_l}}{\alpha_{\omega_l}} \left( \frac{\alpha_{\omega_i}}{\alpha_{\omega_j} \alpha_{\omega_k}} y_{\omega_j} y_{\omega_k} + \frac{\alpha_i}{\alpha_k} y_{\omega_i} y_{\omega_k} - 2 \frac{\alpha_{\omega_i}}{\alpha_{\omega_j}} y_{\omega_i} y_{\omega_j} \right) \\
& + \left[ \sum_{\substack{\omega_i, \omega_j, \omega_k, \omega_l \in \mathcal{I}_{\text{sys}} \\ \omega_i + \omega_j + \omega_k - \omega_l = 0}} 3 \Upsilon_{(2)} \frac{y_{\omega_l}}{\alpha_{\omega_l}} \left( \frac{\alpha_{\omega_i}}{\alpha_{\omega_j} \alpha_{\omega_k}} y_{\omega_j} y_{\omega_k} + 2 \frac{1}{\alpha_{\omega_k}} y_{\omega_i} y_{\omega_k} - \frac{1}{\alpha_{\omega_j}} y_{\omega_i} y_{\omega_j} \right) \right. \\
& \left. + \sum_{\substack{\omega_i, \omega_j, \omega_k, \omega_l \in \mathcal{I}_{\text{sys}} \\ \omega_i - \omega_j - \omega_k = 0}} \Upsilon_{(2)} \frac{y_{\omega_k} y_{\omega_l}}{\alpha_{\omega_k} \alpha_{\omega_l}} \left( \frac{\alpha_{\omega_i}}{\alpha_{\omega_j}} y_{\omega_j} - 3 y_{\omega_i} \right) \right], \quad \omega_i \in \mathcal{I}_{\text{sys}}.
\end{aligned}
\tag{49}
\end{equation}
$$

By inspection, one can notice the following properties of the equations above:

(i) If the normal modes thermalize at the heat bath temperature $T_B$, so that $\langle J_{\omega_i} \rangle = k_B T_B / \omega_i$, it follows that the variables $y_{\omega_i}$ are equal and take the value $y_{\omega_i} = 1$.

(ii) By switching off the external source of energy, i.e., putting $S_{\omega_i} = 0$, the thermal solution above, that is, $y_{\omega_i} = 1$ for all the $\omega_i$, is a stationary solution of the system, namely, $\dot{y}_{\omega_i} = 0$.

In order to understand whether the phonon condensation phenomenon also exists in the above-defined classical framework, one has to work out nontrivial out-of-equilibrium stationary states of the model equations. Unfortunately, doing this analytically for a system with a large number of nonlinear equations is hopeless; therefore, one has to resort to numerical integration of the dynamical equations (49). A difficulty of the numerical approach is to provide a priori estimates of the coupling constants for a real biomolecule. To overcome this problem, we borrow from Ref. [30] the values of physically realistic nonlinear coupling parameters. Moreover, the nondimensional form of the equations partially simplifies the problem, as only the ratios among the coefficients have to be known. On the other hand, since we are only interested in a purely theoretical kind of qualitative study, some simplifying assumptions are in order. Thus, the coupling constants appearing in Eq. (49) should be independent of the frequency; in particular, the condition $S_{\omega_i} = S$ is also assumed (the energy injection rate is the same for all the normal modes). Then, we perform the numerical integration of the dynamical equations by changing the total number $N + 1$ of oscillators (normal modes) in a fixed range ($1 \le \omega_i \le 2$) of equally spaced frequencies given (in arbitrary units) by $\omega_n = 1 + n/N$, with $n = 0, 1, 2, \dots, N$. This choice models an increasing density of modes in a given fixed frequency interval. The classical energy

condensation phenomenon, numerically found and described below, is characterized by a strong deviation from energy equipartition among the normal modes in a stationary state of the dynamical system (47) [and, equivalently, of the system (49)]: As the energy injection rate increases, the system undergoes a major change in the energy distribution among its normal modes, resulting in a more “organized” phase (the energy is mainly channeled into the lowest-frequency mode). Such a condensation phenomenon can be viewed as associated with the breaking of the permutation symmetry for the set of variables $\{y_{\omega_0}, \dots, y_{\omega_N}\}$ in the stationary state, where $y_{\omega_i}$ is the energy content in the mode of frequency $\omega_i$ in $k_B T_B$ units. It is clear that such a symmetry holds at energy equipartition. In analogy with the characterization of symmetry breaking in the classical theory of equilibrium phase transitions, we have to define some (possibly scalar) indicators to detect the permutation symmetry breaking. In order to detect energy equipartition, and its violation in the presence of the condensation phenomenon, we borrow and adapt from Ref. [13] a condensation index $\mathcal{E}_y$, which reads

$$ \mathcal{E}_y = \frac{y_{\omega_0} - \tilde{y}_{\omega_0}}{\sum_{i=0}^{N} y_{\omega_i}}, \qquad (50) $$

where $y_{\omega_0}$ is the energy content of the fundamental mode and $\tilde{y}_{\omega_0} = 1 + S/B$ is the energy content of the fundamental mode given by the exact solution of Eq. (49) in the absence of all the nonlinear couplings (mode-mode and system-bath). Here, we replace the occupation numbers of the modes entering the condensation index of Ref. [13] with the mode energy contents. Note that $\mathcal{E}_y = 0$ corresponds to energy equipartition among the normal modes, and $\mathcal{E}_y = 1$ corresponds to full concentration of the energy into the lowest-frequency mode. In order to easily compare the results obtained at different
numbers *N* of modes, the energy fractions $y_{\omega_i}$ in each normal mode are normalized to give the variables $p_i = y_{\omega_i} / \sum_{\omega_j \in I_{sys}} y_{\omega_j}$.

Another possible choice for the indicator detecting the permutation symmetry breaking is any of the quantities $y_{\omega_i}/y_{\omega_0} = p_i/p_0$, with $i \neq 0$: These are expected to attain their maximum value equal to 1 at equipartition and to take values between zero and 1 for a Fröhlich-like condensation phenomenon. Numerical simulations have been performed in order to characterize this major change of the stationary states of the dynamical system (49), that is, to find the asymptotic stationary solutions of the variables $y_{\omega_i}(t)$ which correspond to the fixed points $\dot{y}_{\omega_i} = 0$ of the system. As a reasonable initial condition, at $t = 0$, the system is taken at thermal equilibrium with the heat bath; that is, all the $y_{\omega_i}(0)$ are set equal to 1.

D. Results of numerical simulations

In Figs. 1 and 2, the results of numerical integration of the rate equations (49) are reported for the set of parameters $B_{\omega_i} = B = 1$, $C_{\omega_i \omega_j} = C = 0.1$, with the initial conditions specified in the previous subsection. The coefficient *B* has

been chosen as a reference parameter because it is directly related to the characteristic thermalization timescale $\tau_{\text{therm}} \approx B^{-1}$. The quartic coupling constant has been chosen as $\Upsilon_{(1)} = \Upsilon_{(2)} = 10^{-4}$, which is sufficient to guarantee a finite bound of the energy from below. The choice $C/B = 0.1$ is coherent, with respect to what has been reported in the literature [23] concerning the already-investigated qualitative aspects of quantum Fröhlich condensation. The integration has been performed using a fourth-order Runge-Kutta algorithm, in each case for a time interval $\tau_{\text{int}} > 0$ sufficiently long so as to guarantee that $|\epsilon_{\omega_i}(\tau) - \epsilon_{\omega_i}(\tau_{\text{stat}})| < 10^{-5}$ for $\tau_{\text{stat}}$, $\tau \in [0, \tau_{\text{int}}]$. In Fig. 1, the reported results display the fraction of energy content in each mode for a system with $N = 20$, that is, 21 normal modes and different values of the energy injection rate. It is quite clear that the larger the energy injection rate *S*, the stronger the deviation from energy equipartition, with increasing energy concentration in the lowest-frequency mode $\omega_0$. This energy concentration phenomenon, which we can call phonon condensation, stems from the classical Liouville equation and displays several analogies with the phenomenology of the pseudo-Bose condensation predicted by Fröhlich in a quantum framework and

<figure><img src="image_3.png" /><figcaption>FIG. 1. Classical Fröhlich-like condensation. Normalized fractions of energy $p_i$ in normal modes vs the mode frequencies of the system $I_{sys}$ for a fixed number of modes: $N = 20$. The histograms display increasing deviations from energy equipartition among the normal modes as the energy input rate *S* increases. The spectra correspond to *S* = 0.1 (blue), *S* = 1 (green), *S* = 10 (purple), and *S* = 100 (pink). Equipartition corresponds to equal heights of the bars. Note the great difference among the energy fractions in the lowest normal mode and in the higher frequency modes.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Classical Fröhlich-like condensation. For an increasing number of modes entering the numerical simulations of Eq. (49), in the upper left panel, the condensation index $\mathcal{E}_y$, defined in Eq. (50), is reported versus the energy input rate $S$, and in the upper right panel, the ratio of the energy content $p_1/p_0$ is reported versus $S$. At equipartition, $\mathcal{E}_0 = 0$ and $p_1/p_0 = 1$; for the condensed state, $\mathcal{E}_y = 1$ and $p_1/p_0 = 0$. The colors of the lines correspond to $N_{sys} = 11$ (blue), $N_{sys} = 21$ (green), $N_{sys} = 41$ (brown), $N_{sys} = 101$ (orange), $N_{sys} = 201$ (purple), and $N_{sys} = 301$ (red). The dashed oblique line (tangent to the inflection point of the two highest $N$ curves) is a guide to the eye to mark a possible asymptotic bifurcation point. In the lower panel, the relative energy content of two groups of modes $-p_i$ for $i = 0-10$ (from blue to purple) and $i = 99, 100, 101$ (blue, green, magenta superposed)—are displayed as a function of $S$, having set $N_{sys} = 301$. The dashed line corresponds to energy equipartition.</figcaption></figure>

subsequently largely studied and debated in the same quantum framework.

A relevant feature of the Fröhlich condensation is the existence of a threshold in the energy injection rate, which corresponds to the breaking of the above-discussed per-mutation symmetry. In order to find out if a similar threshold effect is also displayed by our classical model, we study the behavior of the above-defined indicators of symmetry breaking (the condensation index $\mathcal{E}_y$ and the ratio $y_{\omega_i}/y_{\omega_0}$) as a function of $S$ and for different values of $N$. Figure 1 shows that, for the chosen values of $S$, it is $p_i > p_{i+1}$; hence, the parameter $p_1/p_0$ appears as a natural choice to detect a possible threshold effect. The results of numerical simulations of the order parameters have been reported in Fig. 2. Both indicators pass from their respective equipartition values ($\mathcal{E}_y \approx 0, p_1/p_0 \approx 1$) to their condensation values ($\mathcal{E}_y \approx 1, p_1/p_0 \approx 0$) by displaying a clear steepening at increasing $N$. In particular, at increasing $N$, the curves of the condensation index $\mathcal{E}_y(S, N)$ display a tendency to superpose at large $S$ values, and a sharpening of the transition knee at lower $S$ values, suggesting an

asymptotic convergence to a sharp bifurcation pattern. In fact, this sharpening of the transition curves with $N$ is reminiscent of the finite-size patterns of the order parameter as a function of the control parameter of equilibrium second-order phase transitions. An analogous phenomenon is displayed by the $S$ pattern ($p_1/p_0)(S)$ of the indicator $p_1/p_0$ at increasing $N$. Actually, the right panel of Fig. 2 shows that, around the value $S \approx 10$, the different curves ($p_1/p_0)(S, N)$ cross each other and get steeper at increasing $N$, again suggesting that the larger the number of modes the sharper the transition between almost equipartition among the modes and phonon condensation in the lowest-frequency mode. It is worth noting that this phenomenology is qualitatively confirmed by the outcomes reported in Fig. 5 (aimed at highlighting a saturation effect), where again a sharpening—at increasing $N$—of a threshold effect can be observed. The lower panel of Fig. 2, for the case of $N_{sys} = 201$, provides some interesting and more detailed information about the dynamics underlying the condensation phenomenon reported in the upper panels of the same figure. In fact, the lower panel shows how the
relative energy content $p_i$ of a randomly chosen subset of modes ($p_0, p_1, p_5, p_{10}, p_{100}$) deviates from energy equipartition at increasing $S$. We can observe that the highest frequency mode selected, $p_{100}$, deviates from equipartition at smaller values of $S$ than the lower frequency modes chosen; then, somewhat below $S = 10$, the energy content of all the selected modes starts decreasing, except the energy content $p_0$ of the lowest-frequency mode, which gathers more and more energy, approaching the value $p_0 = 1$, that is, full condensation at high $S$.

The main result of this theoretical part of the paper is an original derivation of a classical counterpart of the well-known Wu-Austin quantum Hamiltonian from which Fröhlich's rate equations can be derived. This classical counterpart consists of a system of nonlinear dynamical equations (49) for the (adimensionalized) energy content $y_{\omega_i}$ of each normal mode $\omega_i$ of an ensemble of oscillators coupled with two thermal baths at different temperatures. The ensemble of oscillators—ideally modeling the vibrations of the atoms or of groups of atoms in a macromolecule—is both linearly and nonlinearly coupled to the set of oscillators representing a thermal bath at temperature $T_B$; this thermal bath accounts for the environment (mainly water molecules) and for other fast internal oscillations of a macromolecule. The second thermal bath at a temperature $T_S \gg T_B$ models an external energy supply, in close analogy with the Wu-Austin model. The numerical study of Eq. (49) has shown that the condensation phenomenon is not peculiar of a quantum model, which *a priori* was not an obvious fact; hence, the occurrence of a phonon condensation phenomenon in a macromolecule at room temperature is more plausible in our classical model than in Fröhlich's quantum model. Therefore, even though Eq. (49) still represents a biomolecule in a very idealized way, the genericity of qualitative features of out-of-equilibrium phase transitions—mentioned at the beginning of this section—suggests that the condensation phenomenon should actually occur in real biomolecules, even though our model is not expected to predict the quantitative details of the phenomenon. This is enough to make the predicted phenomenology worth an effort for its experimental detection.

### III. EXPERIMENTAL DETECTION OF NONEQUILIBRIUM COLLECTIVE MODE

Even though experimental evidence of the existence of collective modes of vibration of biomolecules was provided at thermal equilibrium by means of Raman spectroscopy [1] many years ago and is still the object of many investigations [2–7,31], no experimental evidence is available for the possibility of exciting out-of-thermal-equilibrium collective oscillations of a biomolecule. Unveiling whether these can be activated amounts to understanding whether a necessary condition to activate long-range intermolecular electrostatic forces [8] can be fulfilled. This strongly motivated the joint theoretical and experimental work reported in the

present paper. The experimental counterpart of the theoretical work relies on two complementary experiments.

For both experiments, a model protein has been chosen: the BSA protein. This protein is mainly made out of α helices and is a “model” since it is largely studied in the biophysical chemistry literature. Our strategy to create a stationary out-of-thermal-equilibrium state of this molecule is to induce it by means of optical pumping—without involving any optical transition of the protein—through the excitation of some fluorochromes bound to each protein molecule. The optical excitation of these fluorochromes creates, on each protein, some “hot points” acting as the epicenters of a so-called “proteinquake” [32,33]—discussed in more detail later—and resulting in an energy transfer to the vibrational part of the protein. We used the Alexa488 fluorochrome, which is covalently bonded at the lysine residues of the BSA and which is excited by means of an argon laser (wavelength 488 nm). Some 0.19 eV per fluorochrome and per incident photon (the average energy difference between the absorbed and reemitted photons) is thus available for an energy transfer to the protein and, partly, to its environment. By attaching an average number of 5 fluorochromes per protein, a considerable amount of energy ($\gg k_B T$) can be continuously pumped into each protein. Two THz-near-field absorption spectroscopy setups of aqueous solutions of the protein (at 1 mg/mL concentration), operating in two distinct laboratories, have been used at room temperature [Figs. 3(a) and 3(c)]. In both experiments, THz radiation is produced by tunable, highly spectrally resolved ($<$ 300 Hz) and continuous-wave sources, with an average power of 1 mW, allowing an accurate detection of possible resonances. A typical experiment consists in three phases; during each of them, the aqueous solution of proteins is illuminated with THz radiation, performing a sweep in frequency and thus allowing us to measure the frequency dependence of the absorbed electromagnetic power (detected by the near-field probes) in the solution. During the first phase, no extra illumination is done with the argon laser; during the second phase, the argon laser is switched on to excite the fluorochromes bound to the proteins; finally, during the third phase, the laser is switched off to check whether some memory and irreversible photochemistry effect (photo-bleaching) or sample heating is present. The use of near-field coupling of metal probes to the sample eliminates the Fabry-Perot interferences often seen in optical spectra in fluidic cells [34].

The first setup [Figs. 3(a) and 3(b)] includes a micro-coaxial near-field probe inside a metallic rectangular waveguide, enabling a modal transition from $TM_{01}$ Sommerfeld's to $TE_{01}$ waveguide mode. The subwavelength diameter of the wire (12 µm) allows an extremely localized detection of the longitudinal component of the electric field at its end and in a volume of about 4 pL. The spectra were subtracted from the spectrum of pure water in
<figure><img src="image_5.png" /><figcaption>FIG. 3. Experimental setup of THz absorption near-field spectroscopy. (a) A drop of the biological sample placed under the near-field probe, which is directly immersed in the solution. (b) Picture of the near-field probe and its microwire. (c) A drop of the biological sample placed above the near-field rectenna. (d) Electron-beam microscopy picture of the bow-tie antenna with its integrated FET.</figcaption></figure>

order to remove artifacts coming either from the water absorption or from the geometry of the experimental setup. The second setup [35] [Figs. 3(c) and 3(d)] involves a near-field probe rectenna composed of a planar metal bow-tie antenna, with dimensions close to half a wavelength (at 0.3 THz), that enhanced the THz field in the feed-gap region (volume of about 0.2 pL) and a plasma-wave field-effect transistor (FET) integrated in the feed gap of the antenna. When illuminated by THz radiation, the antenna-coupled FET device provides a dc voltage—between source and drain contacts—proportional to the THz-field intensity [36–38]. A hemispherical silicon lens pressed on the back of the semiconductor substrate focuses the THz radiation on the antenna, simultaneously eliminating Fabry-Perot interference in the substrate [39]. The spectrum of the protein solution obtained in the absence of blue light illumination is subtracted from the spectrum obtained with blue light illumination.

Figure 4(a) presents the spectra obtained using the microcoaxial probe in the absence of blue light illumination (black squares) and in the presence of blue light illumination of different durations (from 3 to 9 min). In the former case (no illumination), there was no specific spectral feature in the studied frequency range, while in the latter case (with illumination), we observed spectral resonances that become more evident for an increasing duration of illumination. In particular, the strongest resonance appeared at 0.314 THz,

accompanied by two other minor resonances situated at 0.278 and 0.285 THz; these values did not depend on the time of illumination, and the strength of the resonances saturated after 9 min of illumination. These results have full reproducibility. Figure 4(b) presents the spectra obtained using the rectenna probe for two durations of blue light illumination. The spectra with illumination on versus illumination off were taken several times (two of them showed reproducibility). In this case, we also observed the appearance of evident resonances whose strength saturated at increasing durations of illumination (10 and 15 min in the two experimental runs reported). Long excitation times are needed because, under our experimental conditions, the energy dissipation rate and the energy supply rate are almost equal, so a long time is needed to accumulate enough energy into each protein in order to make intra-molecular nonlinear interaction terms sufficiently strong to activate the condensation phenomenon.

The spectra obtained using the two previously described methods and for the longer durations of illumination are compared in Fig. 4(c). The principal resonance at 0.314 THz is perfectly reproduced using two completely different and complementary setups. Since THz extinction in water is huge (2000 dB/cm), the emergence of this spectral feature of the protein out of the water absorption background must be associated with the activation of a giant dipole moment. This can happen only as a
<figure><img src="image_6.png" /><figcaption>FIG. 4. Differential transmission and absorption spectra as functions of the frequency. (a) Spectra obtained using the microwire probe after subtraction of the water spectrum in the absence of illumination and in the presence of illumination for the reported durations. (b) Spectra obtained using the rectenna, after subtraction of the protein solution without illumination, for the reported durations. (c) Comparison of the two normalized spectra for the longest illumination durations.</figcaption></figure>

consequence of the activation of a coherent oscillation of the whole molecule, possibly together with a dipole-moment strengthening due to an electretlike structuring of water dipoles of the first hydration layers of the proteins, as the hydration shell might contribute to the magnitude of the protein dipole moment [40].

Since the BSA is a heart-shaped globular protein, the possible lowest frequency of a global oscillatory mode is roughly estimated by schematizing the molecule as composed of two masses $m$, equal to half the total protein mass, joined by a spring of elastic constant $k$ given by $k = EA_0/l_0$, where $E$ is the Young modulus of the protein,

and $A_0$ and $l_0$ are its transverse section and length at rest, respectively. Using $m = 33$ kD, $A_0 \approx 1.2 \times 10^{-13}$ cm$^2$, $l_0 \approx 1.2 \times 10^{-7}$ cm, and $E = 6.75$ GPa, we find $\nu = (1/2\pi)\sqrt{k/m} \approx 0.300$ THz, which is close to the main resonance observed at 0.314 THz. Since the BSA molecule can be modeled to first order as a three-dimensional elastic nanoparticle [6], a more refined approximation is obtained by modeling the protein with an elastic sphere and then considering its vibrational frequencies.

The fundamental frequency of a spheroidal deformation mode of an elastic sphere is given by the formula [41]

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

which holds for $l \ge 2$. For the lowest mode ($l = 2$), we find the frequency $\nu_0 = 0.308$ THz, which agrees within an error of about 2% with the observed peak value at $\nu = 0.314$ THz, using the following data for the BSA protein: Young modulus $E = 6.75$ GPa obtained at room temperature using Brillouin light scattering of hydrated BSA proteins [42], hydrodynamic (Stokes) radius $R_H = 35$ Å, and specific volume $1/\rho = 0.74$ derived from small-angle x-ray scattering (SAXS) experiments [43]. Though such a modeling is unrealistic since it does not take into account the details of the protein structure and the associated normal modes [44], it nonetheless catches a relevant aspect of the global deformation dynamics of the BSA molecule, namely, the activation of a collective oscillation, also suggesting that the physical parameters adopted correspond quite well to the situation investigated [45]. Secondary resonances are also present in both spectra. We could tentatively give a possible explanation by considering torsional modes. These could be activated at the frequencies given by the relation [41]

$$ \nu_t = \nu_0 \left( \frac{(2l + 3)}{2(2l + 1)} \right)^{1/2}, \quad l \ge 2, $$

where $\nu_0$ is given by Eq. (51), whence, for $l = 2$ and $l = 3$, one finds $\nu = 0.257$ THz and $\nu = 0.246$ THz, respectively. These could be associated with the two weaker absorption lines observed at $\nu = 0.278$ THz and $\nu = 0.285$ THz. Here, the larger discrepancy can be attributed to the nonspherical shape of the BSA, which entails the existence of different moments of inertia according to the rotation axis, whereas the breathing mode is insensitive to this. Minor peaks are observed at higher frequencies with the rectenna (falling outside the accessible frequency range of the near-field probe) when the protein solution is illuminated with blue laser light. However, the blue light illumination could produce spurious signals from the 2D electron gas of the FET junction as a consequence of electron-hole pair excitation causing a change of the transistor channel conductivity. This effect is well known
and studied in the literature [46], so minor peaks could be
instrumental artifacts due to this electron-hole-pair creation
effect.

Let us stress an important point: Computational normal mode analysis for proteins has shown nearly continuous vibrational density of states [44], which have also been proven to be nearly uniformly optically active. Moreover, the coupling of these vibrational modes with water results in broad absorption features [47,48]. But this is true at thermal equilibrium, whereas in out-of-equilibrium phonon condensation, the energy content of all the normal modes is strongly depleted, with the exception of the collective mode, hence a narrow absorption feature. In fact, the computation of the function L(ω) in Sec. IV shows how a dipole actively oscillating at a given frequency entails an absorption feature of shape similar to the experimentally observed one.

According to our classical version of the Fröhlich model, it is also expected that the appearance of a collective oscillation should exhibit a thresholdlike behavior when increasing the energy flowing through the protein. Actually, Fig. 5(a) presents a clear threshold in the intensity of the resonance peak at 0.314 THz when the optical input power exceeds 10 µW. By using a classical formalism for the analysis of the out-of-equilibrium phonon condensation, we have calculated the intensity of the normal vibrational modes of the BSA protein as a function of the source power injected through the protein. Figure 5(b) highlights the thresholdlike behavior of the intensity of the fundamental mode that accumulates the energy at the expense of the excited modes, in qualitative agreement with the experimental outcome. By increasing the number of modes included in the calculation, this threshold becomes more and more evident. The experimental and theoretical results reported in Fig. 5 also show the same saturation effect occurring at large values of the energy input rate.

The observed spectra are certainly due to the light-
excited protein because the spectral feature at 0.314 THz is
not observed by illuminating (i) water alone, (ii) an aqueous
solution of the fluorochrome in the absence of the protein,
or (iii) an aqueous solution of the BSA protein without the
bound fluorochrome (see Sec. IV). On the other hand, the
observed spectral line at 0.314 THz immediately disappears
by switching off the laser. Remarkably, the spectra obtained
with two independent and different experimental setups,
based on two different methods of detection of the THz
radiation, operated in two different laboratories, show a
strikingly good overlap of the respective absorption line
profiles at 0.314 THz. This result is in excellent agreement
with the frequency of 0.308 THz predicted for the spher-
oidal (collective) vibrational mode computed on the basis
of recent experimental measurements of the relevant
parameters of the BSA protein. This triple concordance
among the results obtained so far should be enough
to exclude experimental artifacts. Nevertheless, various
sources of the artifact for the observed phenomenology
were considered. One objection is that the observed
phenomenology is just a trivial heating effect due to the
laser light. This would be true in the absence of a
frequency-dependent response to the injection of energy
into the proteins. Heating indeed means increasing the
kinetic energy of the atoms and groups of atoms of the
protein, entailing energy equipartition among the vibra-
tional modes of the protein. Hence, neither a frequency-
dependent effect nor a threshold effect for the energy input
rate would be measured, but a new thermal equilibrium
state would be achieved. To the contrary, each protein—
submitted to continuous energy feeding and energy dis-
sipation—behaves as an open system undergoing a non-
equilibrium phase transition: When the ratio between
energy gain and loss exceeds a critical value, a collective
behavior sets in, producing the phonon condensation, as

<figure><img src="image_7.png" /><figcaption>FIG. 5. Thresholdlike behavior 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 in qualitative agreement.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. Differential transmission spectra of the microwire-based absorption spectra. (a) Without blue light illumination (laser OFF), obtained with solvated AF488 in water (red circles), solvated nonlabeled BSA in water (green triangles), and BSA labeled with AF488 solvated in water (blue squares). (b) With blue light illumination (laser ON), obtained with solvated AF488 in water, at a concentration of 5 mg/mL (red circles); solvated nonlabeled BSA in water (green triangles); and BSA labeled with AF488 solvated in water (blue squares). Labeled and nonlabeled BSA concentrations are kept equal to 1 mg/mL.</figcaption></figure>

has already been discussed in Sec. II. Then, an important
point concerns the conversion mechanism of the visible
light energy absorbed by the electrons of the complex
protein or dyes into the vibrational modes of the proteins.
This mechanism of rapid intramolecular dissipation of
energy through quakelike structural motions as a conse-
quence of a perturbation (such as the breaking of a chemical
bond or the absorption of photons through electronic
transitions) is receiving increasing experimental attention
[32,33,49,50] and is referred to as the “protein quake.”
Similar to an earthquake, this effect describes how a protein
strain is released at a focus or “hot point” (in our case, the
fluorochromes) and then rapidly spreads as a structural
deformation through waves, thus exciting protein vibra-
tional modes. Another source of artifact could be the
apparition of standing waves and related interferences,
but these would be easily identified. Furthermore, there is
no reason for such interferences—if they exist—to manifest
themselves as a consequence of the excitation of fluoro-
chromes through the blue light [see Fig. 6(b)].

IV. MATERIALS AND METHODS

A. Sample preparation

BSA was purchased from Sigma-Aldrich (A7030)
(St. Louis, MO) in lyophilized powder—protease free, fatty
acid free (≤ 0.02%), essentially globulin free; since the degree
of purity of the BSA was higher than 98%, it was used without
further purification and was dissolved directly into bidistilled
water. The BSA molecules were labeled with the fluorochrome
AF488 5-TFP (Alexa Fluor 488 Carboxylic Acid, 2, 3, 5, 6—
Tetrafluorophenyl Ester), 5-isomer (A30005), which was
purchased from Molecular Probes Invitrogen.

The dye has excitation/emission of 495/515 nm and a
molar extinction coefficient of ε<sub>495</sub> = 71 000 M cm<sup>-1</sup>.
Both the protein and dye concentrations have been

determined by measuring their absorbance with a
Nanodrop 1000 Spectrophotometer (ThermoScientific)
at 280 nm with a molar extinction coefficient of
36 000 M cm⁻¹ and at 495 nm with a molar extinction
coefficient of 71 000 M cm⁻¹, respectively. The chemically
labeled BSA molecules were obtained by homemade
labeling using a 2 mg/mL aqueous solution concentration
of proteins with an initial ratio of concentrations [A488
5-TFP]/[BSA] equal to 8 in sodium bicarbonate at a
pH ≈ 8.5, during one hour. Unconjugated dye was
removed using a PD-mini-trap G25 (GE Health Care)
according to the instructions of the manufacturer, using
gravity protocol, and the degree of labeling was determined
spectroscopically. After purification, each BSA molecule
has between 5 and 6 fluorochromes attached.

B. Experimental detection of the collective mode

Two separated and different experiments, performed in Montpellier (France) and in Rome (Italy), respectively, have been performed to obtain the terahertz nonequilibrium spectra of the model protein chosen. The former used a microwire as a local probe, whereas the latter used a nanorectenna.

1. Microwire-based THz spectroscopy

A constitutive element of the experimental setup is a
tunable and continuous-wave primary source from Virginia
Diodes Inc. [51], thus emitting in the 0.22–0.33-THz
frequency range with an average power of 1 mW. The
high spectral resolution (< 300 Hz) of the continuous-
wave source allows an accurate detection of resonances.
The emitted radiation beam is focused on the samples of
protein solutions on which a 488-nm light beam—pro-
duced by an argon laser—can be focused too and switched
on and off at willing. The blue light provides the proteins
with the necessary energy to activate a collective vibra-
tional state. The latter is an out-of-equilibrium state because
it is kept by a nonthermal energy supply.

The THz near-field scanning spectroscopy technique in
aqueous medium is performed by resorting to a homemade
microcoaxial (i.e., subwavelength) near-field probe put
inside a metallic rectangular waveguide connected to a
heterodyne head and an electrical spectrum analyzer. The
subwavelength diameter of the wire allows an extremely
focused enhancement of the longitudinal component of the
electric field at its end over a volume of 4 pL. In 1995,
Keilmann [52] highlighted the advantages of introducing a
metal wire in a circular metal waveguide to produce probes
for near-field microscopy in the far-infrared and microwave
frequency domains. The advantage of this method is to
avoid the frequency cutoff when the diameter of the guide is
a subwavelength one. The circular waveguide is trans-
formed into a coaxial waveguide, which does not have a
low-frequency cutoff and which makes superfocalization
and high-resolution imaging possible.
However, for the experiments reported in the present work, a rectangular waveguide is used, and this entails some different phenomena since the microwire is soldered along the long axis of the waveguide and bent to exit it. This has two main consequences: On the one hand, the microwire enables a modal transition, and on the other hand, it serves as a waveguide. The bent portion of the microwire allows the conversion of the fundamental mode $TE_{01}$ inside the rectangular waveguide into a $TM_{01}$ mode along the wire. To optimize the coupling efficiency between the near field and the probe, special care is taken for the positioning of the microwire inside the guide. More precisely, the maximum of the near-field signal is theoretically attained for a wire positioned at $L_{\text{inside}} = p[(\lambda_{TE_{01}})/4]$ from the open side of the waveguide, where $\lambda_{TE_{01}}$ is the wavelength of the $TE_{01}$ mode at the considered frequency. At 0.3 THz, the microwire must be fixed at an entire multiple $p$ of 250 µm. The best compromise between efficiency and technical possibilities is found for $p = 4$, which gives $L_{\text{inside}} = 1$ mm. The second parameter to take into account is the angular positioning of the bent portion of the microwire inside the rectangular waveguide. The coupling is at a maximum when the wire is parallel to the orientation of the electric field in the guide, that is, along the long axis of the rectangular waveguide. Finally, we note two relevant parameters, which are the total length and the diameter of the wire. The intensity of the electric field is roughly sinusoidal, and its maximum is attained when the total length $L$ is a multiple of the half-wavelength. The best compromise between theory and technological possibilities has given a total length of $L = 2$ mm. Since the electric-field intensity at the microwire extremity increases when the diameter decreases, we use a wire of 12 µm diameter. All the previously mentioned parameters have also been simulated using CSTmicrowaves studio® (https://www.cst.com/products/cstmws) to ensure the best coupling efficiency possible.

All the measurements have been performed at room temperature. A drop of the protein solution sample is placed under the near-field probe, which is directly immersed inside the solution. A typical experiment consists in measuring the near-field electric-field intensity through a reference medium (water) and the protein solution sample. A sweep through the frequency domain accessible to the THz source is performed both when the blue light is switched on and off to yield a difference spectrum showing the absorption features that are attributed to collective vibrations of each protein molecule.

Figure 7 displays the spectra, stabilized in time, obtained by varying the power of the laser illuminating the solution of labeled proteins at a constant concentration of 1 mg/mL. Notice that the secondary peaks appear and develop together with the main peak at 0.314 THz; therefore, these secondary peaks are also vibrational modes of the entire protein, giving some support to the hypothesis attributing

<figure><img src="image_9.png" /><figcaption>FIG. 7. Differential transmission spectra of the microwire-based absorption spectra obtained with different excitation power of the laser. The labeled BSA concentration is kept fixed at 1 mg/mL. The first three spectra are exactly superposed.</figcaption></figure>

them to torsion modes. However, these kinds of modes cannot be accounted for by our theoretical model because it contains no explicit information about the geometry of the molecule.

## 2. Rectenna-based THz spectroscopy

In order to minimize the optical depth of water, in the second experiment (Rome), the probe domain was reduced to a volume of $10 \times 10$ microns in xy (horizontal plane) and to about 2 microns in z (vertical axis). To confine the THz radiation (wavelength $\lambda$ around 1 mm) to such a deeply subwavelength region, a plasmonic antenna is used [35]. This device is based on two main components: a planar metal antenna with length close to $\lambda/2$ (bow-tie, broadband type) that produces a high THz field region with an antenna feed gap of $10 \times 10$ microns; and a plasma-wave FET transistor, which is a THz nonlinear electronic device integrated in the feed gap of the antenna that provides an electric signal proportional to the THz field strength in the feed gap only. The plasma-wave transistor was first introduced by Dyakonov and Shur in 1993 [36] and further developed by many authors [37], the main advantage being that standard microwave transistor technology can be employed for fabricating a device that is sensitive to THz radiation. The device is mounted in a package with a silicon lens pressed on the back of the semiconductor substrate and illuminated from below with a tunable THz oscillator (by Virginia Diodes Inc., 0.18–0.4 THz) through a set of off-axis parabolic mirrors. The resolution of the free-running oscillator is 2 GHz. After acquiring the empty-channel response spectrum of the device, a drop of protein solution was cast on the top (air) side of the device with a micropipette (a volume of 1 microliter was drop-casted). The drop extends over the entire antenna; however, the radiation comes from below (i.e., from inside the substrate), and it is not attenuated by the whole drop. In this
experiment, absorption from molecules outside the antenna feed gap in the *xy* directions can be disregarded. In the *z* direction, the field also extends for less than 2 microns because of the plasma-wave properties (the extension of the plasmonic field in the *z* direction is calculated by finite element modeling simulations). Therefore, the number of molecules probed in this setup is that present in a 10 × 10 × 2 micron volume, 10⁻⁶ times less than the number of molecules present in a 1-mm³ diffraction-limited focus.

### C. Spectroscopic detection of the collective mode

In both experiments, the collective oscillation of the BSA protein is seen as a spectroscopic absorption feature. At variance with standard absorption spectroscopy, where the radiation entering the absorbing medium is responsible for the creation of atomic or molecular excited states, in the experiments reported here, the THz radiation is used to detect an already-excited state of the molecules. Actually, the collective oscillation of the proteins makes them behave as very small antennas (Hertzian dipoles), with the characteristic property of antennas simultaneously absorbing and emitting electromagnetic radiation. However, the absorption along the THz optical path cannot be compensated by the radiation emitted by the oscillating dipoles because it spreads over all the directions in space. The net result is an absorption feature. If we denote by $\vec{\mu}(t)$ the dipole moment of a protein and by $\vec{E}(t) = \vec{E}_0 \cos(\omega t)$ the electric field of the THz radiation, the attenuation of the electric energy density within the drop of protein solution is proportional to the work done by the electric field, that is, $L = -\vec{\mu}(t) \cdot \vec{E}(t)$. The oscillation of the dipole moment is necessarily damped, predominantly because of bremsstrahlung emission, so that, denoting by $\tau$ the lifetime of the activated collective oscillation and by $\omega_c$ its frequency, we can set $\vec{\mu}(t) = \vec{\mu}_0 e^{-t/\tau} \cos(\omega_c t)$. Thus, after averaging over all the relative orientations and all the phase differences $\phi$ between $\vec{\mu}(t)$ and $\vec{E}(t)$ such that the electric field does a positive work, we obtain

$$L(\omega) = 2 \int_{0}^{\pi} d\phi \int_{0}^{\infty} dt \mu_{0} E_{0} e^{-t/\tau} \omega_{c} \sin(\omega_{c} t) \cos(\omega t + \phi). \quad (52)$$

This is the elementary attenuation process of the THz radiation. This process is repeated in time for each molecule at a rate proportional to the intensity of the drop illumination with the blue light. Moreover, the total attenuation is proportional to the concentration of absorbing molecules in the protein solution. Equation (52) gives, for $L(\omega)$, a Lorentzian shape centered at $\omega_c$, the resonance frequency of the collective oscillation of the BSA protein. Figure 10 shows three different shapes of the function $L(\omega)$ obtained for different values of $\tau$ (in arbitrary units).

<figure><img src="image_10.png" /><figcaption>FIG. 8. Raw data of the rectenna-based absorption spectra. Black squares correspond to the THz emitting source spectrum without the protein solution. Red circles correspond to the absorption spectrum of the BSA protein solution without blue light illumination. Green triangles correspond to the absorption spectrum of the BSA protein solution recorded after 10 min of blue light illumination. Blue triangles correspond to the absorption spectrum of the BSA protein solution recorded after 15 min of blue light illumination. BSA concentration is equal to 1 mg/mL.</figcaption></figure>

These line shapes show that an already-vibrating underdamped dipole absorbs the weak terahertz-radiation probe with the same frequency pattern as the 0.314-THz absorption line reported in Fig. 3. The latter is well fitted by the $L(\omega)$ function by using a $Q$ factor of 50 ($Q = \Delta\nu/\nu$, that is, the ratio between the line width and the line frequency).

In Figs. 6(a) and 6(b), some control spectra are reported, which have been obtained with the microwire antenna. The two groups of spectra refer to the blue light emitting laser switched off and on, respectively. The observed absorption line at 0.314 THz is clearly found only when the fluorochrome AF488 is bound to the BSA molecules and in the presence of 488-nm laser light illumination of the protein solution. These spectra rule out any other origin of the observed absorption feature besides the proposed one in the main text of the present work.

In Fig. 8, the photoresponse spectra obtained with the rectenna probe highlight the same phenomenology: The absorption line at 0.314 THz is present only when the aqueous solution of AF488-labeled BSA molecules is illuminated with 488-nm laser light. A clear difference is again observed between the spectra when the blue light illumination is switched off and on. Artifacts possibly due to blue-light-illuminated water alone, or possibly due to the blue-light-illuminated AF488 dye in aqueous solution, are ruled out by the control spectra reported in Fig. 9 and are in agreement with the analogous controls performed with the microwire antenna reported in Figs. 6(a) and 6(b).
<figure><img src="image_11.png" /><figcaption>FIG. 9. Control data of the rectenna-based absorption spectra. We show the difference spectra obtained with and without blue laser illumination of water (red triangles) and of an aqueous solution of AF488 fluorochrome at a concentration of 5 mg/mL (blue triangles). The spectra have been recorded after 10 min of blue light illumination.</figcaption></figure>

D. Experimental activation timescale of the collective mode

In both Montpellier and Rome experiments, the absorp-
tion feature, which is attributed to the activation of an out-
of-equilibrium collective oscillation of the protein,
appeared after several minutes of blue light irradiation:
about 4 min in the Montpellier experiment and after 10 min
in the Rome experiment (though, in this latter case, no
measurement was performed before 10 min). In what
follows, we provide a consistency estimate of this long
activation time with the physical conditions of the experi-
ments. Going beyond a qualitative estimate is hardly
feasible and, in any case, would not affect the meaning

<figure><img src="image_12.png" /><figcaption>FIG. 10. The function $L(\omega)$ of Eq. (52) is plotted in arbitrary units for three different values of $\tau$. The blue line corresponds to $\tau = 1$, the red line to $\tau = 5$, and the black line to $\tau = 10$, with $\omega = 10$.</figcaption></figure>

of the experimental outcomes reported in the present work.
Thus, we proceed by estimating the balance between the
energy input rate and the energy loss rate of each molecule.
As we discuss at the end of this section, this long time is not
a hindrance to the biological relevance of the phenomenon
reported in the present work.

An elementary account of the balance between energy
gain and loss for each protein can be given by the equation

$$
\frac{dE}{dt} = -\frac{2}{3} \frac{(Ze)^2 |\ddot{\mathbf{x}}|^2}{c^3} - \Gamma + W, \quad (53)
$$

where *E* in the lhs is the numerical value of the energy of
the system described by the Hamiltonian (2) of the previous
section, the first term in the rhs represents the radiative
losses (bremsstrahlung) of the accelerated dipole of
moment *Zex* (where *Z* is the number of elementary charges
*e* and **x** is the separation between the center of positive and
negative charges), and *W* is the energy input rate due to the
488-nm light irradiation. Note that Γ stands for other kinds
of energy losses leading to thermalization of the protein
with its aqueous environment.

For what concerns the energy input rate *W*, the energy
difference between the entering photons, of wavelength
4.88 × 10⁻⁵ cm, and the fluorescently emitted photons, of
wavelength 5.3 × 10⁻⁵ cm, amounts to 3.18 × 10⁻¹³ erg =
0.19 eV. When the argon laser operates at 500 µW, the
number of emitted photons per second is 1.2 × 10¹⁵.
The cross section of the fluorochrome Alexa488 is
3.2 × 10⁻¹⁶ cm² (free dye in water) [53], so assuming a
Gaussian beam section (at 1/e²) of 7.8 × 10⁻³ cm² (cor-
responding to a diameter of about 0.1 cm), we find that each
Alexa488 molecule receives approximately 50 photons per
second. Each protein molecule has an average number of
5–6 fluorochromes bound to it so that, if we consider that
some of them could be partly shadowed by the protein
itself, we can assume that the average number of photons
received by each protein per second is in the interval
120–300. Hence, the upper bound for *W* is estimated to be
between 3.8 × 10⁻¹¹ erg s⁻¹ and 9.5 × 10⁻¹¹ erg s⁻¹. This
estimate is obtained under the assumption that all the
photons hitting the Alexa488 molecules are absorbed; the
quantum yield of these dye molecules is high (92%), so
practically all of the harvested photons are converted into
fluorescence, thus providing an upper bound of the energy
available for the protein excitation. Of course, part of this
energy can be dissipated in the surrounding environment in
the form of heat.

For what concerns the radiation losses, we proceed to
estimate the time average value of |**x**|². Since 63% of the
N-isoform of the BSA protein consists of α helices, we can
assume that the so-called accordionlike vibrational modes
of these helices provide the largest contribution to the
protein extension and thus to the variation of the protein's
dipole moment. At thermal equilibrium, we can assume that
the dipole elongation $\mathbf{x}(t)$ is the result of an incoherent superposition like $\mathbf{x}(t) = \sum_{i=1}^{30} A_i \cos(\omega_i t + \xi_i)$, where the $\xi_i$ are random phases and the sum runs over the 30$\alpha$ helices of which the protein is composed. To compute the time average of $|\ddot{\mathbf{x}}|^2$, we have to estimate the quantity $\sum_{i=1}^{30} \frac{1}{2} A_i^2 \omega_i^4$ (since the mixed terms average to zero). Using the formula $\omega_i = 2\pi(1/2L_i)(E/\rho)^{1/2}$ for the frequencies of the accordionlike modes of the $\alpha$ helices [54], where $1/\rho = 0.74$, $E = 2.31 \times 10^{11}$ dyne cm$^{-2}$ for $\alpha$ helices, $L_i = 1.5N$ Å, where $N$ is the number of amino acid residues, and assuming that $A_i$ belongs to the interval $0.05L_i-0.1L_i$, and $Ze = -13 \times 4.803 \times 10^{-10}$ Fr for the BSA protein at neutral pH [55], we finally find

$$ \frac{2}{3} \frac{(Ze)^2}{c^3} |\ddot{\mathbf{x}}|^2 \approx 3.25 \times 10^{-12} - 1.7 \times 10^{-11} \text{ erg s}^{-1}. \quad (54) $$

The range of values of the energy input rate and that of the radiative losses are thus almost overlapping. Making these estimates more precise is hardly feasible and is beyond the aim of the present work. What matters here is that, since we experimentally observe the activation of the collective mode of the BSA molecules, the energy input rate $W$ must exceed the rate of all the losses, radiative and nonradiative; if this happens only by a small amount, then $dE/dt$ can be so small as to require some minutes in order to accumulate enough energy in each molecule [this can be thought of as a steady increase of the value of the third term of the unperturbed part of the Hamiltonian in Eq. (2)]. When the energy accumulated exceeds a critical threshold value, the system undergoes a condensation transition, channeling the largest fraction of the input energy into the lowest-frequency mode(s). An elementary estimate of the rate of energy losses due to the bremsstrahlung in the condensed phase is obtained by entering into the Larmor formula (54) the dipole acceleration $|\ddot{\mathbf{x}}|$ computed from $\mathbf{x}(t) = A \cos(\omega_c t)$, putting $\omega_c = 0.314$ THz (which is the frequency of the collective mode), and assuming that the dipole elongation $A$ is about 10 Å. The break-even point between the energy input and the radiative losses is then found for a dipole moment in the range 14 500–23 000 Debye, corresponding to an effective number of charges $Z$ approximately in the range 290–460. In other words, in the condensed phase, the total dipole of the protein oscillates at a low frequency (with respect to the accordionlike modes of the $\alpha$ helices), so the collective oscillation can be stable with respect to radiative losses up to the activation of large values of the protein dipole moment.

We have neglected any estimate of the collisional losses—represented by $\Gamma$ in Eq. (53)—because their estimate as viscous losses—on the basis of the Stokes formula—gives an unreasonably high energy dissipation rate with respect to the energy input rate. If this were the case, no collective vibration would be observed at all. There is a vast amount of literature about relaxation phenomena, and thus

thermalization, of proteins subject to different kinds of radiative excitations. The typical relaxation timescales are in the range of picoseconds to nanoseconds, with a remarkable exception—recently reported in Ref. [31]— where the thermalization of terahertz photons in a protein crystal happens on a microsecond to millisecond timescale, which is interpreted as being due to Fröhlich condensation (either quantum or classical). On the contrary, typical picosecond to nanosecond relaxations usually pertain to side-chain relaxations or small groups of atoms. Also, in the case of relaxations of collective modes detected with FIR spectroscopy or with THz time domain spectroscopy [56,57], the physical conditions are very different with respect to our present experiments. In the mentioned literature, the collective modes of proteins are probed at thermal equilibrium—thus, in the presence of all the other vibrational modes and under the action of weak excitation fields. In our experiments, the protein molecules are out of thermal equilibrium and are strongly excited through an internal cascade of mode-mode couplings stemming from the fluorescence decay of the attached fluorochromes, and the THz radiation is only used to see the presence of the collective oscillation (see the preceding subsection).

On the other hand, it has recently been found that the hydration shell of the BSA protein is 25 Å thick [40], and this seems to be a generic property of solvated proteins [7]. The microrheology [58] of this kind of water-protein system—and, in particular, its high-frequency viscoelasticity—is still an open research subject, making a quantitative estimate of the term $\Gamma$ in Eq. (53) hardly feasible.

A comment about the prospective biological relevance of the observed phenomenology is in order. The main energy source within living cells is provided by ATP hydrolysis. The typical intracellular concentration of ATP molecules is given around 1 mM, implying that a protein molecule in the cell undergoes around $10^6$ collisions with ATP molecules per second [59]. Given the standard free energy obtained from ATP hydrolysis estimated around $50 \text{ kJ.mol}^{-1} = 8.306 \times 10^{-13} \text{ erg}$, if we assume that only 1% or 2% of the collisions with ATP provide energy, a power supply between $8.306 \times 10^{-9} \text{ erg s}^{-1}$ and $1.6 \times 10^{-8} \text{ erg s}^{-1}$ is potentially available. This could be at least 2 orders of magnitude larger than the power supplied to each protein in our experiments, but reasonably, it could be much more than 2 orders of magnitude because we have assumed 100% conversion efficiency of the energy supplied by the laser into internal vibrations of the protein. However, this cannot be the case; thus, the condensation mechanism in vivo can be considerably faster.

## V. CONCLUDING REMARKS

In the present paper, we have studied a classical version of a quantum model put forward many years ago by Fröhlich [9,10,12]. The classical model displays the same phenomenon of a Bose-like phonon condensation of the
normal vibrational modes of a macromolecule that was predicted by the original quantum model. Even though our classical model (just like its quantum predecessor) is very idealized, that is, even if our model cannot be quantitatively predictive, it can nevertheless qualitatively observe the presence of a robust and generic phenomenon undergone by open systems with internal nonlinear interactions: a nonequilibrium phase transition when a control parameter exceeds a critical threshold value. This kind of phase transition corresponds to a channeling into the lowest-frequency mode of almost all the energy pumped into the system, that is, to a collective global molecular vibration. This theoretical result has motivated an experimental—and successful—attempt at confirming it. So, in the second part of the present paper, we have reported—for a model protein—the unprecedented observation of an out-of-thermal-equilibrium collective oscillation, which is in qualitative agreement with the theoretical model and in excellent quantitative agreement with the theoretically expected value of a spheroidal collective vibration of the whole protein.

The novelty is not the collective oscillation itself because several terahertz spectroscopic studies have reported collective modes of proteins; however, all of these studies were performed at thermal equilibrium and mainly carried out using dry or low-hydrated powders because of the very strong absorption of water [2–5], even though more recent studies also addressed solvated proteins [6,7]. We stress at this point that the novelty of both our theoretical and experimental contributions consists of considering out-of-thermal-equilibrium conditions. On the other hand, recent studies on solvated BSA in the THz [60] and sub-THz frequency range have shown [40] broad resonances due to an efficient coupling of low-frequency modes of the protein with the surrounding water; though all of these works were performed at thermal equilibrium, as in these previous studies, the experimental part of our present work confirms the relevance of the coupling of the protein with the surrounding water molecules. In fact, the strong absorption feature that we observed in an aqueous solution of the BSA protein placed out of thermal equilibrium reveals that the protein vibrating in its collective mode has to be dressed by ordered layers of water molecules in order to attain an effective dipole moment sufficiently large enough to overcome the strong absorption of bulk water.

We anticipate that the theoretical and experimental sides of the work presented in this paper could open a broad domain of systematic investigations about out-of-equilibrium activation mechanisms and properties of collective oscillations of different kinds of biomolecules. Furthermore, as already mentioned in the Introduction, the possibility of exciting out-of-thermal-equilibrium collective oscillations of macromolecules is specifically interesting as a necessary condition to activate resonant long-distance electrodynamic intermolecular interactions

[8]. Thus, our results explain why electrodynamic interactions between biomolecules have hitherto eluded detection; in fact, no attempt has ever been done to detect them by involving biomolecules vibrating out of equilibrium. Consequently, our work also motivates new efforts to detect these electrodynamic intermolecular interactions [61,62].

## ACKNOWLEDGMENTS

We thank Jack A. Tuszynski, Michal Cifra, and Anirban Bandyopadhyay for useful discussions, and A. Penarier, T. Cohen, and F. Cano for their help in the simulation and the realization of the microwire local probes. The project leading to this publication received funding from the Excellence Initiative of Aix-Marseille University—A*Midex, a French “Investissements d’Avenir” programme. This work was also financially supported by the Seventh Framework Programme for Research of the European Commission under FET-Proactive Grant TOPDRIM (FP7-ICT-318121), by the projects SIDERANT and NEBULA financed by the French CNRS, by the Occitanie Region, and by Montpellier University through its TOP platform and NUMEV LabeX.

I. N., M. L., and J. T. performed the experiments in Montpellier. J. T. designed and built the experimental setup in Montpellier with the support of P. N., V. G. and M. O. performed the experiments in Rome. The classical formulation of Fröhlich condensation stems from the Ph.D thesis of M. G. J. P. contributed the theoretical part. I. D. did some numerical work. J. S. contributed to the biophysical and biochemical aspects of the work. J. T. and L. V. contributed to conceiving the Montpellier experiments and participated in the supervision of all the development of the project. M. P. in quality of project leader intervened in all aspects of the work. All the authors contributed to the discussion and to the analysis of the results. M. P. wrote the paper with the help of M. G., M. O., J. S., J. T., and L. V.; J. T. made all the figures.

[1] P. Painter, L. Mosher, and C. Rhoads, *Low-Frequency Modes in the Raman Spectra of Proteins*, Biopolymers **21**, 1469 (1982).

[2] D. Turton, H. M. Senn, T. Harwood, A. J. Lapthorn, E. M. Ellis, and K. Wynne, *Terahertz Underdamped Vibrational Motion Governs Protein-Ligand Binding in Solution*, Nat. Commun. **5**, 3999 (2014).

[3] G. Acbas, K. A. Niessen, E. H. Snell, and A. G. Markelz, *Optical Measurements of Long-Range Protein Vibrations*, Nat. Commun. **5**, 3076 (2014).

[4] R. J. Falconer and A. G. Markelz, *Terahertz Spectroscopic Analysis of Peptides and Proteins*, J. Infrared Milli Terahz Waves **33**, 973 (2012).

[5] A. G. Markelz, A. Roitberg, and E. Heilweil, *Pulsed Terahertz Spectroscopy of DNA*, Bovine Serum Albumin
and Collagen between 0.1 and 2.0 THz, Chem. Phys. Lett. **320**, 42 (2000).

[6] J. Xu, K. W. Plaxco, and S. J. Allen, Probing the Collective Vibrational Dynamics of a Protein in Liquid Water by Terahertz Absorption Spectroscopy, Protein Sci. **15**, 1175 (2006).

[7] S. Ebbinghaus, S. J. Kim, M. Heyden, X. Yu, M. Gruebele, D. M. Leitner, and M. Havenith, Protein Sequence- and pH-Dependent Hydration Probed by Terahertz Spectroscopy, J. Am. Chem. Soc. **130**, 2374 (2008).

[8] J. Preto, M. Pettini, and J. A. Tuszynski, Possible Role of Electrodynamic Interactions in Long-Distance Biomolecular Recognition, Phys. Rev. E **91**, 052710 (2015).

[9] H. Fröhlich, Long-Range Coherence and Energy Storage in Biological Systems, Int. J. Quantum Chem. **2**, 641 (1968).

[10] H. Fröhlich, Long Range Coherence and the Action of Enzymes, Nature (London) **228**, 1093 (1970).

[11] H. Fröhlich, Selective Long Range Dispersion Forces between Large Systems, Phys. Lett. A **39**, 153 (1972).

[12] H. Fröhlich, Long-Range Coherence in Biological Systems, Riv. Nuovo Cimento **7**, 399 (1977).

[13] J. R. Reimers, L. K. McKemmish, R. H. McKenzie, A. E. Mark, and N. S. Hush, Weak, Strong, and Coherent Regimes of Fröhlich Condensation and Their Applications to Terahertz Medicine and Quantum Consciousness, Proc. Natl. Acad. Sci. U.S.A. **106**, 4219 (2009).

[14] I. Altfeder, A. A. Voevodin, M. H. Check, S. M. Eichfeld, J. A. Robinson, and A. V. Balatsky, Scanning Tunneling Microscopy Observation of Phonon Condensate, Sci. Rep. **7**, 43214 (2017).

[15] Y. Sun, P. Wen, Y. Yoon, G. Liu, M. Steger, L. N. Pfeiffer, K. West, D. W. Snoke, and K. A. Nelson, Bose-Einstein Condensation of Long-Lifetime Polaritons in Thermal Equilibrium, Phys. Rev. Lett. **118**, 016602 (2017).

[16] In a classical context, the expression “Bose-like condensation” is used to signify that the largest part of the energy supplied to a macromolecule is channeled into its lowest-frequency normal mode.

[17] H. Haken, Cooperative Phenomena in Systems Far from Thermal Equilibrium and in Nonphysical Systems, Rev. Mod. Phys. **47**, 67 (1975).

[18] H. Haken, Synergetics—An Introduction: Nonequilibrium Phase Transitions and Self-Organization in Physics, Chemistry, and Biology (Springer, Berlin, 1977).

[19] T. M. Wu and S. Austin, Bose Condensation in Biosystems, Phys. Lett. A **64**, 151 (1977).

[20] T. M. Wu and S. Austin, Bose-Einstein Condensation in Biological Systems, J. Theor. Biol. **71**, 209 (1978).

[21] T. M. Wu and S. Austin, Cooperative Behavior in Biological Systems, Phys. Lett. A **65**, 74 (1978).

[22] T. M. Wu and S. Austin, Fröhlich's Model of Bose Condensation in Biological Systems, J. Biol. Phys. **9**, 97 (1981).

[23] J. Pokorny and T. M. Wu, Biophysical Aspects of Coherence and Biological Order (Springer, Berlin, Heidelberg, 2013).

[24] H. Bolterauer, Elementary Arguments that the Wu–Austin Hamiltonian has no Finite Ground State (The Search for a Microscopic Foundation of Fröhlichs Theory), Bioelectrochem. Bioenerg. **48**, 301 (1999).

[25] P. Kramer and M. Saraceno, Geometry of the Time-Dependent Variational Principle in Quantum Mechanics,

in Group Theoretical Methods in Physics (Springer, New York, 1980), pp. 112–121.

[26] P. Kramer, A Review of the Time-Dependent Variational Principle, J. Phys. Conf. Ser. **99**, 012009 (2008).

[27] H.-R. Jauslin and D. Sugny, Dynamics of Mixed Classical-Quantum Systems, Geometric Quantization and Coherent States, in Mathematical Horizons for Quantum Physics, edited by H. Araki, B.-G. Englert, L.-C. Kwek, and J. Suzuki, IMS-NUS Lecture Notes Series, Vol. 20, (World Scientific, Singapore, 2010).

[28] This formulation is equivalent to the Schrödinger representation in quantum mechanics.

[29] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.8.031061 for detailed derivation of Eq. (47) from Eq. (46).

[30] F. Piazza and Y.-H. Sanejouand, Energy Transfer in Nonlinear Network Models of Proteins, Eur. Phys. Lett. **88**, 68001 (2009).

[31] I. V. Lundholm, H. Rodilla, W. Y. Wahlgren, A. Duelli, G. Bourenkov, J. Vukusic, R. Friedman, J. Stake, T. Schneider, and G. Katona, Terahertz Radiation Induces Nonthermal Structural Changes Associated with Fröhlich Condensation in a Protein Crystal, Struct. Dyn. **2**, 054702 (2015).

[32] A. Ansari, J. Berendzen, S. F. Bowne, H. Frauenfelder, I. E. Iben, T. B. Sauke, E. Shyamsunder, and R. D. Young, Protein States and Protein Quakes, Proc. Natl. Acad. Sci. U.S.A. **82**, 5000 (1985).

[33] M. Levantino, G. Schirò, H. T. Lemke, G. Cottone, J. M. Glownia, D. Zhu, M. Chollet, H. Ihee, A. Cupane, and M. Cammarata, Ultrafast Myoglobin Structural Dynamics Observed with an X-ray Free-Electron Laser, Nat. Commun. **6**, 6772 (2015).

[34] T. Globus, D. Woolard, T. W Crowe, T. Khromova, B. Gelmont, and J. Hesler, Terahertz Fourier Transform Characterization of Biological Materials in a Liquid Phase, J. Phys. D **39**, 3405 (2006).

[35] V. Giliberti, A. Di Gaspare, E. Giovine, M. Ortolani, L. Sorba, G. Biasiol, V. V. Popov, D. V. Fateev, and F. Evangelisti, Downconversion of Terahertz Radiation Due to Intrinsic Hydrodynamic Nonlinearity of a Two-Dimensional Electron Plasma, Phys. Rev. B **91**, 165313 (2015).

[36] M. Dyakonov and M. Shur, Shallow Water Analogy for a Ballistic Field Effect Transistor: New Mechanism of Plasma Wave Generation by dc Current, Phys. Rev. Lett. **71**, 2465 (1993).

[37] W. Knap et al., Nonresonant Detection of Terahertz Radiation in Field Effect Transistors, J. Appl. Phys. **91**, 9346 (2002).

[38] P. Nouvel, H. Marinchio, J. Torres, C. Palermo, D. Gasquet, L. Chusseau, L. Varani, P. Shiktorov, E. Starikov, and V Gružinskis, Terahertz Spectroscopy of Plasma Waves in High Electron Mobility Transistors, J. Appl. Phys. **106**, 013717 (2009).

[39] A. Di Gaspare, R. Casini, V. Foglietti, V. Giliberti, E. Giovine, and M. Ortolani, Terahertz Current Oscillations in a Gated Two-Dimensional Electron Gas with Antenna Integrated at the Channel Ends, Appl. Phys. Lett. **100**, 203504 (2012).

[40] O. Sushko, R. Dubrovka, and R. Donnan, Sub-terahertz Spectroscopy Reveals that Proteins Influence the Properties
of Water at Greater Distances than Previously Detected, J. Chem. Phys. **142**, 055101 (2015).

[41] S. I. Bastrukov, Low-Frequency Elastic Response of a Spherical Particle, Phys. Rev. E **49**, 3166 (1994).

[42] S. Perticaroli, J. D. Nickels, G. Ehlers, H. O'Neill, Q. Zhang, and A. P. Sokolov, Secondary Structure and Rigidity in Model Proteins, Soft Matter **9**, 9548 (2013).

[43] E. Mylonas and D. I. Svergun, Accuracy of Molecular Mass Determination of Proteins in Solution by Small-Angle X-ray Scattering, J. Appl. Crystallogr. **40**, s245 (2007).

[44] K. Suhre and Y. Sanejouand, A Normal Mode Web-Server for Protein Movement Analysis and the Generation of Templates for Molecular Replacement, Nucleic Acids Res. **32**, W610 (2004).

[45] In Ref. [42], the authors seem to place serious limitations on the possibility of using the elastic global model [Y. Uezaki and N.Gō, Int. J. Peptide Protein Res. **7**, 333 (1975), where globular proteins in the native state are assumed to behave as continuous elastic spheres in the low-frequency breathing motions] because the resonance frequency observed with neutron scattering and with light scattering is independent of the linear dimensions of the protein, thus contradicting the model. However, the authors refer to the boson peak of the studied proteins. The boson peak is an equilibrium glasslike feature of proteins due to an excess of low-frequency states. The boson peak of the BSA protein is a broad feature peaked at 30 cm⁻¹, far from the narrow feature reported in the present work, which is peaked at 10 cm⁻¹.

[46] S. Boubanga-Tombet et al., Room Temperature Coherent and Voltage Tunable Terahertz Emission from Nanometer-Sized Field Effect Transistors, Appl. Phys. Lett. **97**, 262108 (2010).

[47] M. Heyden, D. J. Tobias, and D. V. Matyushov, Terahertz Absorption of Dilute Aqueous Solutions, J. Chem. Phys. **137**, 235103 (2012).

[48] J. Sun, G. Niehues, H. Forbert, D. Decka, G. Schwaab, D. Marx, and M. Havenith, Understanding THz Spectra of Aqueous Solutions: Glycine in Light and Heavy Water, J. Am. Chem. Soc. **136**, 5031 (2014).

[49] D. Arnlund et al., Visualizing a Protein Quake with Time-Resolved X-ray Scattering at a Free-Electron Laser, Nat. Methods **11**, 923 (2014).

[50] L. U. L. Brinkmann and J. S. Hub, Ultrafast Anisotropic Protein Quake Propagation after CO Photodissociation in Myoglobin, Proc. Natl. Acad. Sci. U.S.A. **113**, 10565 (2016).

[51] https://www.vadiodes.com/en/.

[52] F. Keilmann, FIR Microscopy, Infrared Phys. Technol. **36**, 217 (1995).

[53] L. Wang, A. K. Gaigalas, and V. Reipa, Optical Properties of Alexa™488 and Cy™5 Immobilized on a Glass Surface, Biotechniques **38**, 127 (2005).

[54] K. Itoh and T. Shimanouchi, Vibrational Frequencies and Modes of α-helix, Biopolymers **9**, 383 (1970).

[55] L. R. Barbosa, M. G. Ortore, F. Spinozzi, P. Mariani, S. Bernstorff, and R. Itri, The Importance of Protein-Protein Interactions on the pH-induced Conformational Changes of Bovine Serum Albumin: A Small-Angle X-ray Scattering Study, Biophys. J. **98**, 147 (2010).

[56] A. Xie, A. F. G. van der Meer, and R. H. Austin, Excited-State Lifetimes of Far-Infrared Collective Modes in Proteins, Phys. Rev. Lett. **88**, 018102 (2001).

[57] G. R. Kneller, Quasielastic Neutron Scattering and Relaxation Processes in Proteins: Analytical and Simulation-Based Models, Phys. Chem. Chem. Phys. **7**, 2641 (2005).

[58] T. A. Waigh, Microrheology of Complex Fluids, Rep. Prog. Phys. **68**, 685 (2005).

[59] B. Alberts et al., Molecular Biology of the Cell, 4th ed. (Garland Science, New York, 2002).

[60] J. W. Bye, St. Meliga, D. Feradom, G. Cinquel, J. A. Zeitler, and J. Falcones, Analysis of the Hydration Water around Bovine Serum Albumin Using Terahertz Coherent Synchrotron Radiation, J. Phys. Chem. A **118**, 83 (2014).

[61] I. Nardecchia, L. Spinelli, J. Preto, M. Gori, E. Floriani, S. Jaeger, P. Ferrier, and M. Pettini, Experimental Detection of Long-Distance Interactions between Biomolecules through Their Diffusion Behavior: Numerical Study, Phys. Rev. E **90**, 022703 (2014).

[62] I. Nardecchia et al., Detection of Long-Range Electrostatic Interactions between Charged Molecules by Means of Fluorescence Correlation Spectroscopy, Phys. Rev. E **96**, 022403 (2017).