Out-of-equilibrium collective oscillation as phonon condensation in a model protein

In the first part of the present paper (theoretical), the activation of out-of-equilibrium collective oscillations of a macromolecule is described 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 non-equilibrium 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, to 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 long time ago by H. Fr\"ohlich in the attempt to explain the huge speed of enzymatic reactions. In the second part of the present paper (experimental), we show that such a phenomenon is actually possible. Two different and complementary THz near-field spectroscopic techniques, a plasmonic rectenna, and a micro-wire near-field probe, have been used in two different labs to get rid of artefacts. By considering a aqueous solution of a model protein, the BSA (Bovine Serum Albumin), we found 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, and in excellent quantitative agreement with a theoretical result allowing to identify the observed spectral feature with a collective oscillation of the entire molecule.


I. INTRODUCTION
In the impressively complex molecular machinery at work in living cells, the biochemical players seem "to know" where to go and when. According to a longstanding working hypothesis, this could be explained if the encounters of distant cognate partners of biomolecular reactions would be actively driven by selective (resonant) attractive forces of electrodynamic nature. This was mainly put forward by H. Fröhlich [1][2][3][4], originally to account for the huge speed of enzyme reactions. In order to excite sufficiently intense and resonant electrodynamic attractive interactions, each of the cognate partners of a biochemical reaction should undergo collective vibrational oscillations at the same or almost the same frequency. And this would happen because molecular collective oscillations would bring about giant electric dipole oscillations.
In particular, this could explain why even in the context of low molecular concentration of one of the cognate partners of a reaction there is still a high efficiency and rapidity of mutual encounters of the reagents which can hardly be the result of random encounters only (Brownian diffusion). This is particularly evident in some ligand-receptor reaction or in the transcription machinery of the DNA. The molecular collective oscillations, driven by metabolic activity, in Fröhlich's proposal [1] are seen as a Bose-like condensation of the normal vibrational modes of a biomolecule. However, the condensation phenomenon originally surmised by Fröhlich has been criticized and marginalized because its quantum formulation can be hardly maintained to be a realistic one for biomolecules. In fact, the frequency of collective oscillations is expected in the sub THz domain, around 10 11 Hz, hence at room temperature it is k B T / ω 1 (where k B and are the Boltzmann and Planck's constants respectively) so that the thermal wavelength of the atoms, of which a biomolecule is composed, is so short as to make a quantum description unnecessary. And in fact, for example, the vibrational properties of proteins are very well described by Molecular Dynamics simulations performed in a classical context. Of course this does not apply to molecular orbitals and light induced electronic transitions, as bound electrons always behave quantum mechanically. Therefore, we have worked out a classical version of the original Fröhlich model, finding that -remarkably -also in a classical context a Bose-like phonon condensation phenomenon is possible. In this context a biomolecule is considered as an open system through which energy flows under the simultaneous action of an external supply and of dissipation due to radiative, dielectric, and viscous energy losses. The classical Bose-like condensation in the lowest vibrational mode occurs when the energy input rate exceeds some threshold value. Moreover, an important theoretical prediction [5] shows interactions. Materials and methods are described in Section IV. In Section V we draw some conclusion.

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 having in common the fact of being composed of many non-linearly 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 non-equilibrium phase transition [6,7]. The physically relevant fact is the genericity of this phenomenology once the mentioned basic ingredients are there. This means that even idealised 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 '60s by H. 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's rate equations can be derived -was later given by Wu and Austin [8][9][10][11]  In what follows we indicate with a † ω i , a ω i the quantum creation/annihilation operator for the vibrational normal modes of the main system (i.e. a biomolecule) with frequency ω i ∈ I sys .
Such a system is put in contact with a thermal bath which represents the degrees of freedom of the environment surrounding the protein and, possibly, other normal modes of the protein which 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 Ω j ∈ I bth whose annihilation/creation operators are b Ω j and b † Ω j , 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 T B . Also in this case the corresponding thermal bath is described by a collection of harmonic oscillators with frequencies Ω k ∈ I src , the quantum annihilation/creation operators of which are c Ω k and c † Ω k . These three sets of harmonic oscillators can be regarded as three subsystems of a larger isolated system S (we coherently indicate with I S the set of all the normal modes of the system) whose quantum dynamics is described by the Hamiltonian where H 0 is the free Hamiltonian of the three sets of harmonic oscillators representing the molecular normal modes and the two heat baths The interactions among normal modes are described byĤ Int ; in the original formulation by Wu and Austin the interaction term has the form: where 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. [8,9] and in the reference book in Ref. [12].
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 [13]. This problem can be easily fixed by adding a term with a quartic dependence on the creation/annihilation operators of the form so that the lower bound to the ground energy level does not go to −∞ for large values of This quartic interaction stands for a anharmonic interaction among the normal modes of the biomolecule, a broadly studied topic of relevance to energy transport in biomolecules like proteins.  [14,15]. The same technique, more deepened from both the mathematical and conceptual points of view, has been proposed as a "dequantization" technique in [16] as a kind 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 details, one begins with the ansatz that the wavefunction depends on N parameters where the parameters x i = x i (t) are in general functions of time. For a quantum system with Hamiltonian H T ot the equations of motion of the x i can be derived using the following variational principle (equivalent to the Least Action principle) where L(ψ, ) is the Lagrangian associated to the system 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 H tot over the state |ψ(x 1 , ..., x N ) , that is The Poisson brackets {·, ·} depend only on the chosen parametrization for the wavefunction.
Starting from the variables one defines the antisymmetric tensor so that the equations of motion are implicitly given by If the condition DetW ij = 0 holds, then the matrix W ij = W −1 ij defines the Poisson brackets for the classical Hamiltonian system This formalism can be applied to the quantum system described by the quantum Hamiltonian of Eq.(1) to associate to it a classical Hamiltonian system. The choice of the parametrization for the wavefunction 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/annihilation operators of the quantum harmonic oscillators describing the system, the wavefunction is been chosen as a product of the corresponding coherent states. In particular where |z A i (t) sys , |z B k (t) bth , z Cp (t) src are normalized coherent states for the normal modes of the main system, of the thermal bath and of the external source, respectively: their general form is given by From the definition of coherent states in Eq. (14) it follows that the expectation value for the occupation number n i is given by the squared norm of z so, as we are interested in writing rate equations for these quantities, we parametrize the wavefunction ψ(t) with the set of real parameters {(n i , θ i )} i∈I S such that Using Eq.(9) it is possible to derive the Poisson brackets associated with the variables {(n i , θ i )} i∈I S : and, consequently, using the definition (10) the entries of the matrix W are and its inverse has the form Consequently it follows that the variables J ω = n ω and θ ω are canonically conjugated variables. The classical Hamiltonian H = H 0 +H IntW A +H IntQuad for the variables {(θ ω , J ω )} ω∈I S is given by a free classical part by a semi-classical Wu and Austin interaction part and by the quartic term where each complex coupling constant is given in polar representation. In what follows, the coupling constants are considered real and rescaled s.t.
with these choices, the total Hamiltonian of the system reads In order to derive Fröhlich-like rate equations, the dynamics of the action variables J ω i of the system has to be studied. We could choose to investigate the dynamics of the system by In the following Section classical Fröhlich-like rate equations are derived by resorting to the time evolution of the distribution function ρ({(J ω , θ ω )} ω∈I S , t) satisfying the Liouville equation.
B. Derivation of classical rate equations using the Koopman-Von Neumann formalism Let ρ({(J ω , θ ω )} ω∈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 ρ associated with this Hamiltonian is given by where ı 2 = −1 and {·, ·} are the canonical Poisson brackets and L H (·) = ı {H, ·} is the Liouville operator acting on functions defined on the phase space In our case the Hilbert space of complex square integrable functions in phase space is L 2 (Λ {(Jω,θω)} ω∈I S ) with the inner product defined by with f, g ∈ L 2 (Λ {(Jω,θω)} ω∈I S ). On this space we can define the action of the Liouville operator and consider the domain DL Let ψ({(J, θ)} I S ; t) ∈ DL H be a normalized time-dependent function [17] such that then it can be proved that ρ = ψ L 2 (Λ (J,θ) ) = ψ * ψ is a normalized function for which (29) holds. Moreover asL H is a self-adjoint operator, it represents the unitary time evolution of the "wave function" as in analogy with quantum mechanics. With this formalism the rate equations for the average values of the actions (which are the counterpart of quantum occupation numbers) associated with the normal modes ω i of the main system are given by whereM Jω i is a multiplicative operator acting on L 2 (Λ {(Jω,θω)} ω∈I S ) as followŝ The Liouville operator can be decomposed aŝ and since the eigenfunctions of the operatorL H 0 are known, the action of the operator L H Int can be treated as a time dependent perturbation, which is adiabatically turned on and off from t 0 = 0 to t = +∞. Then a classical analogous of the interaction representation formalism is used. Thus, if |ψ(t) S is the "wave function" in the Schrödinger representation, then in interaction representation |ψ(t) I reads and, given a generic operatorÂ S in the Schrödinger picture, its expression in the interaction pictureÂ I risÂ With this formalism the time evolution of |ψ(t) I can be written through the unitary evolution operatorÛ (t; t 0 ) satisfying and the formal solution of (41) is given by At first order inL H int (t ), the unitary evolution operatorÛ (t; t 0 ) in the right hand side of eq. is substituted by the identity operator meaning that the state |ψ(t 0 ) I , if the perturbation is turned on at t 0 , at the lowest order can be approximated by |ψ(t 0 ) I |ψ(t 0 ) S and assumed to be coincident with Schrödinger picture (i.e. |ψ(0) = |ψ 0 ), so that and, asL H int (t ) is self-adjoint, the time evolution for the "bra" has the form The time derivative of the multiplicative operatorM Jω i in interaction picture is derived This means that the average (35) can be entirely rewritten using the interaction picture as d dt After lengthy computations (see SM for the 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  [1], 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 Eqs. (47) in non-dimensional form after introducing the following variables so that Eqs. (47) reaḋ By inspection, one can notice the following properties of the equations above • if the normal modes thermalize at the heat bath temperature T B , so that J ω i = k B T B /ω i , it follows that the variables y ω i are equal and take the value y ω i = 1; • by switching-off the external source of energy, i.e. putting S ω i = 0, the thermal solution above, that is y ω i = 1 for all the ω i , is a stationary solution of the system, namelẏ In order to understand whether the phonon condensation phenomenon exists also in the above defined classical framework, one has to work out non trivial out of equilibrium stationary states of the model equations. Unfortunately, doing this analytically appears 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 borrowed from Ref. The classical condensation phenomenon found numerically, as is detailed below, belongs to the family of non-equilibrium phase transitions mentioned at the beginning of Section II: by slightly varying the energy injection rate (the control parameter ) across a critical value, the system undergoes a major change in the energy distribution among its normal modes resulting in a more "organized" phase (the energy is mainly channelled into the lowest frequency mode). To characterize this major change of dynamical state, the equivalent of an order parameter for a finite set of normal modes is defined as where y ω i is the energy in the mode of frequency ω i in k B T B units. When the system is at thermal equilibrium, from the energy equipartition among the normal modes it follows that ε ω i = 0; whereas, for a non-equilibrium energy distribution among the normal modes it is The range of each parameter ε ω i is [−1/N sys ; 1 − 1/N sys ]; the lower bound corresponds to a system at an effective temperature T ω i = 0, while the upper bound corresponds to the condition where all the energy is concentrated in one normal mode ω i . We look for the asymptotic stationary solutions of the parameters ε ω i (t), corresponding to the fixed pointṡ y ω i = 0. As initial condition, the system was taken at equilibrium with the heath bath, that is, all the y ω i (0) were set equal to 1.

D. Results of numerical simulations
In Fig.1 the results of numerical integration of the rate equations (49)  has been reported in the literature [12] concerning quantum Fröhlich condensation. The integration has been performed using a fourth-order Runge-Kutta algorithm, in each case for a time interval τ int > 0 sufficiently long so as to guarantee that Three qualitatively different stationary states can be observed according to the value of S (the energy injection rate in adimensional units): • an equipartition regime (blue dots in Fig.1) for S 0.5 × 10 −1 there is energy equipartition among all the oscillators which attain an equilibrium stationary state at a somewhat higher temperature than the thermal bath, that is, T stat > T B ; • a partial condensation regime (green and purple dots in Fig.1) for 0.5×10 −1 S 100 the energy equipartition among the normal modes is violated, and some of the lowest frequencies normal modes contain much more energy than the other modes; • a saturated condensation regime (red dots in Fig.1) for S 100 the vibrational energy is entirely condensed in the lowest normal mode and there is practically no energy left in the higher frequencies normal modes.
By using also C = 0.01 and C = 0.001, these results have been found to be qualitatively stable. For the higher value C = 1, the results are close to those found for C = 0.1, what suggests a saturation effect at increasing values of C. Actually, this saturation effect is clearly evident in Figure 2, where the phenomenon of condensation is displayed by the deviation from energy equipartition of some low frequency modes as a function of the energy input parameter. Likewise, the saturation effect at increasing value of the a-dimensional energy input rate is well evident in Figure 5 reported in the following Section.
The effect of quartic direct interactions among normal modes on the order parameters ω i as functions of the energy injection rate S has been investigated in classical framework for the two frequency sets I sys1 and I sys2 . Some different tests have been done for C = 0.1 and and different values of Υ (1) and Υ (2) . To compare the results for the two systems I sys1 and I sys2 , the coefficients of the terms derived from quartic interactions have been chosen so that Υ (1) δω sys =Υ (1) and Υ (2) δω sys =Υ (2) , whereΥ (1) andΥ (2) are supposed to be constant.
The presence in the rate equations of terms proportional to Υ (1) does not significantly af-

MODE
Even though experimental evidence of the existence of collective modes of vibration of biomolecules has been provided at thermal equilibrium by means of Raman spectroscopy [19] already many years ago, and is still being the object of many investigations [20][21][22][23][24][25][26] The first setup (Fig. 3 (a), (b)) used a micro-coaxial near-field probe put inside a metal-  Since the BSA molecule can be modeled to first order as a threedimensional elastic nanoparticle [24], 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 [36] which holds for l ≥ 2. 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 [37], hydrodynamic (Stokes) radius R H = 35Å, and specific volume 1/ρ = 0.74 derived from small-angle X-ray scattering (SAXS) experiments [38], we find for the lowest mode (l = 2) the frequency ν 0 = 0.308 T Hz which agrees within an error of about 2% with the observed peak value at ν = 0.314 T Hz. Though such a modeling is unrealistic in what it does not take into account the details of the protein structure and the associated normal modes [39], 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. Secondary resonances are also present in both spectra. A possible explanation could be tentatively given considering torsional modes. These could be activated at the frequencies given by the relation [36] ν t = ν 0 (2l + 3) 2(2l + 1) where ν 0 is given by equation (51) is well known and studied in the literature [40] so that 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 [39] which have also been proved nearly uniformly optically active. Moreover, the coupling of these vibrational modes with water results in broad absorption features [41,42]. But this is true at thermal equilibrium, whereas µ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  is being given increasing experimental attention [27,28,43,44] and is referred to as "protein quake". Similarly 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 vibrational modes. Another source of artifact could be the apparition of standing waves and related interferences, but these would have been easily identified. Furthermore, there would be no reason for such interferences - The THz near-field scanning spectroscopy technique in aqueous medium is performed by resorting to a homemade micro-coaxial (i.e. subwavelength) near-field probe put inside a metallic rectangular waveguide connected to an heterodyne head and an electrical spectrum analyser. 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, F. Keilmann [45] 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 microwaves frequency domains. The advantage of this method is to avoid the frequency cut-off when the diameter of the guide is a subwavelength one. The circular waveguide is transformed into a coaxial waveguide which does not have low frequency cut-off and which makes superfocalisation and high-resolution imaging possible.
However, for the experiments reported in the present work, a rectangular waveguide was used, and this entails some different phenomena since the micro-wire is soldered along the long axis of the waveguide and bent to exit it. This has two main consequences: on the one hand, the micro-wire enables a modal transition and, on the other hand, it serves as a waveguide. The bent portion of the micro-wire allows the conversion of the fundamental mode T E 01 inside the rectangular waveguide into a T M 01 mode along the wire. To optimize the coupling efficiency between the near-field and the probe, a special care has been first paid to the positioning of the micro-wire inside the guide. More precisely, the maximum of the near field signal is theoretically attained for a wire positioned at L inside = p λ T E 01 4 from the open side of the waveguide, where λ T E 01 is the wavelength of the T E 01 mode at the considered frequency. At 0.3 THz, the micro-wire must be fixed at an entire multiple p of 250 µm. The best compromise between efficiency and technical possibilities was found for p = 4 that gives L inside = 1 mm. The second parameter to take into account is the angular positioning of the bent portion of the micro-wire inside the rectangular waveguide.
The coupling is 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 also paid attention to two relevant parameters that 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 micro-wire extremity increases when the diameter decreases we used a wire of 12 µm diameter. All the previously mentioned parameters have also been simulated using CSTmicrowaves studio R (https://www.cst.com/products/cstmws) to ensure the better coupling efficiency as possible.

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 × 10 microns in xy (horizontal plane), and to about 2 microns in z (vertical axis). To confine the THz radiation (wavelength λ around 1 mm) to such a deeply sub-wavelength region, a plasmonic antenna is used [30]. This device is based on two main components: a planar metal antenna with length close to λ/2 (bow-tie, broadband type) that produces a high THz field region with antenna feed  the balance between the energy input rate and the energy loss rate of each molecule. As we shall discuss at the end of this Section, this long time is not an 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 where E in the l.h.s is the numerical value of the energy of the system described by the Hamiltonian (2) 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 aims 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 On the other hand, it has been recently found that the hydration shell of the BSA protein is 25Å thick [35], and this seem to be a generic property of solvated proteins [25]. The microrheology [51] 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 Γ in equation (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 [52]. Given the standard free-energy obtained from ATP hydrolysis estimated around 50 kJ.mol -1 = 8.306 × 10 −13 erg, if we assume that only 1% or 2% of the collisions with ATP will provide energy, a power supply between 8.306 × 10 −9 erg s -1 and 1.6 × 10 −8 erg s -1 is potentially available. This could be at least two orders of magnitude larger than the power supplied to each protein in our experiments, but reasonably even much more than two orders of magnitude because we have assumed a hundred percent conversion efficiency of the energy supplied by the laser into internal vibrations of the protein. But this can be hardly 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 H. Fröhlich [1,2,4]. The novelty is not the collective oscillation in itself, because several terahertz spectroscopic studies have reported about collective modes of proteins, but all of these studies were performed at thermal equilibrium and mainly carried on using dry or low-hydrated powders because of the very strong absorption of water [20][21][22][23], even though more recent studies also addressed solvated proteins [24,25]. The novelty of both our theoretical and experi-mental contributions, let us stress this point again, consists of considering out of thermal equilibrium conditions. On the other hand, recent studies on solvated BSA in THz [53] and sub THz frequency range have shown [35] broad resonances due to an efficient coupling of low frequency modes of the protein with the surrounding water, and though all of these works were performed at thermal equilibrium, in common with 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 a watery solution of the BSA protein put 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 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-ofthermal-equilibrium collective oscillations of macromolecules is specifically interesting as a necessary condition to activate resonant long-distance electrodynamic intermolecular interactions [5]. 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 [54,55].