On the existence of soliton-like collective modes in liquid water at the viscoelastic crossover

The problem of large-density variations in supercooled and ambient water has been widely discussed in the past years. Recent studies have indicated the possibility of nanometer-sized density variations on the subpicosecond and picosecond time scales. The nature of fluctuating density heterogeneities remains a highly debated issue. In the present work, we address the problem of possible association of such density variations with the dynamics of terahertz longitudinal acoustic-like modes in liquid water. Our study is based on the fact that the subpicosecond dynamics of liquid water are essentially governed by the structural relaxation. Using a mode coupling theory approach, we found that for typical values of parameters of liquid water, the dynamic mechanism coming from the combination of the structural relaxation process and the finiteness of the amplitude of terahertz longitudinal acoustic-like mode gives rise to a soliton-like collective mode on a temperature-dependent nanometer length scale. The characteristics of this mode are consistent with the estimates of the amplitudes and temperature-dependent correlation lengths of density fluctuations in liquid water obtained in experiments and simulations. Thus, the fully dynamic mechanism could contribute to the formation and dynamics of fluctuating density heterogeneities. The soliton-like collective excitations suggested by our analysis may be relevant to different phenomena connected with supercooled water and can be expected to be associated with some ultrafast biological processes.

The existence of density variations in both supercooled and ambient water on the subnanometer/nanometer length scales and subpicosecond/picosecond time scales has been widely discussed in recent years  . Bulk water under ambient conditions is traditionally considered in the framework of a homogeneous model as a mainly tetrahedrally hydrogen-bonded, continuous liquid [2][3][4][5][7][8][9]13,14,21,22,30,31 . Within this model, the temporal existence of low-and high-density regions in liquid water can be explained as a result of natural density fluctuations of an equilibrium system. For instance, molecular dynamics (MD) simulations based on the SPC/E model showed the presence of transient low-and high-density (< 0.9 or > 1.1 kg/m 3 ) regions in water for both supercooled and ambient conditions, interpreted within a continuum model. The correlation length of the corresponding density fluctuations was estimated to be of 0.7 nm, and the timescale for low-and high-density regions was estimated as < 4 ps 8 . Likewise, numerical simulations by Soper 3 showed that on the length scale of 0.9 nm, density fluctuations can, in principle, be up to ~ ± 30% compared to the average water density.
Different local and global order parameters have been used in order to assess the presence of low-and highdensity regions in liquid water 26,27 . Recently, coexisting low-and high-density regions in ambient water have been identified using Voronoi voids on a larger length scale of 1-2 nm than those predicted by local order parameters and those reported previously 19 . Atomistic simulations showed that the low-and high-density regions are formed around the regions of empty space within the hydrogen-bond network, occurring as a result of thermal fluctuations 18 . The regions of empty space are frequently formed as asymmetric, fractal voids that are highly delocalized over the hydrogen-bond network, and the corresponding nanometer-sized density fluctuations can be associated with these voids. The results of simulations 18 suggested a picosecond time scale for voids with volumes larger than 0.15 nm 3 , while small voids can be characterized by a subpicosecond timescale. Similar results were obtained for supercooled water 29 .
Large-density fluctuations in liquid water can be explained also in the framework of heterogeneous models based on fluctuations between two classes of different structures 12 , known as low-density liquid (LDL) and high-density liquid (HDL) structures 1 , locally favoured and normal structures 32  www.nature.com/scientificreports/ In the simulations 27 the effects connected with LDL-like and HDL-like species were quite local since a significant part of water molecules was sufficiently close to the average; thus the observed LDL-like and HDLlike species represented rather extreme and rare fluctuations. However, Iwashita et al. 17 showed that the local molecular dynamics of ambient water can possibly reflect local fluctuations into LDL 17 . This conclusion is based on the results of direct experimental measurements of real-space, real-time motions of water molecules via the van Hove function determined by inelastic X-ray scattering. Iwashita et al. 17 showed that water molecules are strongly correlated in space and time with coupling between the first and second nearest-neighbour molecules. Water seems to be different from some other liquids in this respect 17 . The observed local dynamics were found to be crucial to a fundamental understanding of the origin of the physical properties of water, including viscosity and molecular transport 17 .
Starting structures classified by the LSI, seeding in the simulations 27 , may be treated as both extreme, rare spontaneous fluctuations or externally induced perturbations. As described above, these provide formation of LDL-like and HDL-like regions that are well-defined on nanometer length and picosecond time scales and strongly influence the translational dynamics of water 27 . Such density variations may be relevant in various cases, such as dynamical heterogeneity in supercooled water, a phenomenon where spatially separated, extended regions of relatively mobile and immobile molecules coexist 66 . Specifically, recent simulations showed that in supercooled water low-mobility regions have a lower local density (and higher degree of tetrahedrality), while high-mobility regions have a higher local density 61,62 . This finding has significant implications for the understanding of diverse phenomena associated with supercooled water; for example, it was shown that ice nucleation occurs in the low-density regions 67 .
Large-density variations may also be relevant in the following context. Inelastic neutron scattering measurements allow to estimate that density fluctuations in biomolecular hydration water can be affected by the biomolecule surface at a distance up to 1.2 nm from the surface 68 . The results of THz spectroscopy suggest that this distance may be significantly larger. It was found that the perturbation induced by some proteins to the collective dynamics of hydration water can extend to about 2 nm or more from the protein-water interface 69 . Also, MD simulations showed the similarity of the acoustic-like modes in the protein and protein hydration water 64,70 . Such long-range protein-water interactions associated with heterogeneous hydration dynamics were proven to be able to contribute to the activity of proteins 64 . The possibility of similar long-range interactions between bilayer lipid membranes and solvent in the terahertz frequency range was indicated in Ref. 71 .
Large-density variations can also arise in such natural processes as vibrations of nanoparticles in water in the terahertz frequency range which produce large-amplitude perturbations of adjacent water, regions of rarefication 72,73 . We finally note that recent simulations revealed the existence of propagating nanometer-sized (~ 5 nm) density pulses in liquid water, induced by a large-amplitude initial perturbation 74 .
In the present work, we address the problem of possible association of the discussed nanometer-sized density variations with the dynamics of terahertz longitudinal acoustic-like modes in liquid water. Our analysis is based on the fact that the subpicosecond dynamics of liquid water are essentially governed by a (α-) structural relaxation process, i.e., the time decay of density fluctuations induced by viscous dissipation. Within the viscoelastic crossover region, the structural relaxation causes a positive sound dispersion, i.e., upward bending of the wave number-dispersion of the longitudinal acoustic-like mode from the low-frequency sound dispersion to the highfrequency dispersion [75][76][77][78][79][80][81] . Using a mode coupling theory approach, we show that for typical values of parameters of both bulk water and protein/DNA hydration water, the dynamic mechanism coming from the combination of the structural relaxation process and the finiteness of the amplitude of terahertz longitudinal acoustic-like mode gives rise to the existence of soliton-like collective modes on a temperature-dependent nanometer length scale. We also show that the characteristics of these modes are consistent with the estimates of the amplitudes and temperature-dependent length scales of density fluctuations in liquid water obtained in experiments and simulations.

Model
We consider dynamics of longitudinal acoustic-like mode in liquid water on the length scale corresponding to the viscoelastic transition region in the framework of a mode-coupling theory (MCT) approach. Note that a model based on MCT reproduces correctly the experimental results of the time-resolved optical Kerr effect measurements of low-frequency vibrational dynamics in water over the temperature range of 247-353 K 11,82 . MCT can be formulated in the framework of a generalized hydrodynamic approach 83,84 . Several generalized hydrodynamic models have been developed. They are typically based on retaining the formal structure of the classical hydrodynamic theory and replacing macroscopic thermodynamic and transport coefficients with appropriate wave number-dependent memory functions. We follow the method proposed by Zwanzig et al. 85 .
Based on the data of Liao et al. shown in Fig. 5 of Ref. 86 , we may approximately consider the heat mode to be decoupled from the sound mode for the parameters under study. Thus we neglect thermal relaxation and consider an isothermal case.
The dynamic structure factor of the system can be defined as 87 where δρ(q, t) = dre iq·r δρ(r, t) is the Fourier transform of the mass density fluctuation from the equilibrium value ρ 0 , q is the wave vector, ω is the angular frequency, t is the time, m is the molecular mass, and the bracket <> implies a thermal ensemble average. www.nature.com/scientificreports/ In the linear approximation, the generalized hydrodynamic equations for an isotropic system are reduced to the dynamic structure factor determined by the relation 87 Here ω 0 (q) is the excitation frequency, Ŵ(q) is the damping factor, S(q) is the static structure factor, q = q ; with k B and T being the Boltzmann constant and the absolute temperature, respectively.
Due to the viscoelastic transition, on increasing q , ω 0 (q) gradually bends upward from the linear adiabatic sound dispersion to the infinite frequency dispersion [75][76][77][78][79][80][81] . The onset of a substantial positive sound dispersion induced by the viscoelasticity accompanies also the decrease of temperature under supercooling conditions 80,88 . Accordingly, we use the following approximation for ω 2 0 (q) up to wave numbers corresponding to the upper boundary of the viscoelastic transition region q u : with the coefficients k 1 and k 2 being constant. It should be emphasized that the approximation given by Eq. (2) is valid only for wave numbers less than ~ q u ; this equation does not have a physical meaning for larger q . We estimate the values of k 1 and k 2 from a least-squares fit to the dispersion curves for ω 0 obtained in inelastic X-ray scattering experiments 81 for water at ≈ 278 K (case I) and in MD simulations for SPC/E model 89 for D 2 O molecules at the temperature of maximum density for the SPC/E potential (case II). We get and from the positions of ω 0 (q) shown, respectively, in Fig. 2 of Ref. 81 for q up to q u ≈ 3.0 nm −1 and in Fig. 2a of Ref. 89 for q up to q u ≈ 3.1 nm −1 . Note that we also complemented the data of Ref. 81 by those from the inelastic ultraviolet scattering and Brillouin light scattering measurements 90 for the wave number region from 0.02 to 0.1 nm −1 ; however, this did not noticeably change the values of k 1 and k 2 .
(2) and Ŵ(q) given by Eq. (3), with η L being the generalized longitudinal viscosity, corresponds to a linearized version of the following MCT model: Here r = (x 1 , x 2 , x 3 ) are the rectangular Cartesian coordinates, u = (u 1 , u 2 , u 3 ) is the velocity vector, F[ρ] is the free energy of the system, δF[ρ]/δρ is the functional derivative of F[ρ] with respect to the density ρ (the same symbols are used for the functional derivative and for the density fluctuation δρ which can lead to a confusion, www.nature.com/scientificreports/ thus we must conclude from the context which quantity is meant), S is the deviatoric stress tensor that is represented as a sum of the mean and fluctuating components, η L and η S describe, respectively, the response of longitudinal and shear stresses to a change in the strain rate, δ αβ is the Kronecker delta, δ(t) is the delta function, and the asterisk denotes complex conjugation. Equations (4)- (7) are treated as the equations of motion of hydrodynamic variables 84 .
Assuming that the relative amplitudes of the density fluctuations are small but finite, we set for the finite amplitude density fluctuations in the first approximation, instead of Eq. (6): where α is a constant. Evidently, the coefficient α in Eq. (8) corresponds to the second derivative of pressure with respect to density, and thus to 91,92 where β is the acoustic nonlinear coefficient and c the speed of sound. For water at room temperature for low frequency 92 β ≈ 3.5. To obtain estimates for α , we first estimate the apparent sound velocity, i.e., c(q) = ω 0 (q)/q , in the following way. The center of the viscoelastic crossover corresponds to the inflection point of the sound dispersion curve. Accordingly, for the dispersion curves given in Refs. 81,89 , the characteristic wave number for the region where the viscoelastic transition predominantly takes place is chosen to be q * = 2.0 nm −1 (case I) and q * = 2.5 nm −1 (case II). Then we get, respectively, ω 0 (q * ) ≈ 4.9 ps −1 and ω 0 (q * ) ≈ 6.8 ps −1 , and consequently, c(q * ) ≈ 2.5 × 10 3 m/s and c(q * ) ≈ 2.7 × 10 3 m/s. Then, assuming that β ≈ 3.5 and using Eq. (9), we have the following estimates in the case I and II, respectively. Note that the viscoelastic transition signifies a transformation from the lowfrequency liquid-like behavior to the high-frequency solid-like behavior, while typically β is higher for solids than for liquids 92 . Thus the above estimates for α may be considered as the estimates from below.
In order to estimate the characteristic magnitude of density fluctuation, we note that the root mean square deviation of the number of molecules N in a specified volume V , expressed as a fraction of the average number of molecules in the volume, is determined by 21 where κ T is the isothermal compressibility. The characteristic magnitude A * of density fluctuation δρ can be determined by this relation for a volume of water of size 1 nm 3 . (Note that the length scale of 1 nm is, probably, the smallest that can be examined outside the range of the local structure around individual molecules 21 .) We thus estimate A * ≈ 43 kg/m 3 for water at ≈ 277 K ( ρ 0 ≈ 1.0 × 10 3 kg/m 3 ).

Results and discussion
It is reasonable to focus on the dynamics of density variations on the time interval 0.5-1 ps, as shown by the time scales relevant to the considered system, listed below. In the context of our model, we note that it is well-known that the two factors we consider, the nonlinear dispersion and nonlinear amplitude effects, can provide soliton 94 (cf. Ref. 95 ). In the idealized limit case η L = 0 , δS ≡ 0 , Eq. (13) is known as the Boussinesq equation, a nonlinear wave equation that indeed possesses exact n -soliton solutions 96 . A one-soliton solution of the form 97 corresponds to a relatively stable localized density pulse. Here the parameters A , v , and δ characterize, respectively, the amplitude, velocity, and width of a pulse. These parameters are connected with each other through Eqs. www.nature.com/scientificreports/ (16a) and (16b), and δ can be considered as a free parameter, given that the condition of reality for v is fulfilled. In view of Eq. (16b), the propagation velocity will be less than the sound velocity and will decrease with decreasing δ. Under the condition Equation (16) determines the exact solution of Eq. (14) in the case η L = 0 , δS ≡ 0: Formally, this solution requires the fulfillment of the condition δS ≡ 0 . However, numerical simulations (see "Methods") of the corresponding non-dimensional problem (Eqs. (15) and (15a)) with the indicated above values of the parameters showed that the contribution of the stochastic term ∂ 2 δS/∂x 2 into the profile of δρ is negligibly small on the time interval up to 1 ps, for the initial profile given by Eqs. (17). Thus the exact solution (17) can be considered as a good approximation to the solution of Eqs. (14) and (14a) on this time interval.
This result suggests the existence of a stationary soliton-like mode in ambient water on the time scale of 1 ps relevant to the considered system. In order to test this hypothesis, we compare the values of the amplitude and length scale of a localized density variation, following from our model, with available experimental and numerical data. At the beginning, it should be emphasized that in the viscoelastic transition region the dynamics of density fluctuations are strongly coupled with the formation and break-up of the hydrogen-bond network 78 . Thus the considered soliton solution may be interpreted as a stationary nanosized region of rarefication associated with rearrangement of hydrogen-bond network at the peak of the pulse.
From Eqs. (17a) and (17c), we obtain the following estimates for A , δ , and the full width at half maximum (FWHM) of a pulse. In the case I: A ≈ 3.9 × 10 2 kg/m 3 , δ ≈ 0.70 nm, and the FWHM of a pulse is ≈ 1.2 nm. In the case II: A ≈ 2.6 × 10 2 kg/m 3 , δ ≈ 0.97 nm, and the FWHM of a pulse is ≈ 1.7 nm. As indicated above, one may expect that α would be somewhat higher than the value given by Eq. (10) and, consequently, the amplitude A would be lower.
The obtained values of the amplitude A of the soliton density pulse turn out to be comparable with an expected difference between the densities of LDL-like and HDL-like species, which is of the order of 200-300 kg/ m 3 (Refs. 21,[34][35][36]. It should be emphasized that these LDL-like and HDL-like species were numerically observed as extreme and rare fluctuations 27 , as discussed in the introductory part of the paper in detail. The considered density variations can also be treated as those induced by externally applied perturbations. For instance, starting structures classified by the LSI, seeding in the simulations of Camisasca et al. 27 , may be treated, in fact, as initial perturbations which produce LDL-like and HDL-like regions well-defined on nanometer length and picosecond time scales. Considering external perturbations, we note that nanometer-sized density variations with amplitudes of the order of 2-4 × 10 2 kg/m 3 could be induced by mechanical forces in the biological range. Indeed, for a surface of size ~ 1 nm 2 , the effective pressure p induced by external force of the order of 2 × 10 -9 N can be estimated as ~ 2 × 10 9 Pa. Then from the relation ∂p/∂ρ = c 2 with c = 2.5 × 10 3 m/s as the characteristic value of the apparent sound velocity c(q * ) (see earlier), one could estimate δρ ~ 3.2 × 10 2 kg/m 3 , in the very first approximation.
We next compare the obtained values of the length scale with available experimental and numerical results. It can be seen that the length scale of 1-2 nm, suggested by the soliton solution given by Eqs. (17), is consistent with the following results. (i) The molecular arrangement around the LDL-like and HDL-like species, determined by the extreme values of the LSI, in the simulations 27 is well-defined over time at the distances ~ 2.2 nm. (ii) X-ray scattering studies 1,16 determined the correlation length ξ of density fluctuations to be of 0.3-0.4 nm. This correlation length can be related to an average length scale of density fluctuations, assuming a spherical shape with the diameter ~ 4.5 ξ , which gives ~ 1.5 nm at 280 K 1,16 . It should be noted that in many cases density variations will have sizes both smaller and larger than the indicated average size 8,16 . (iii) The length scale of low-and high-density regions in ambient water identified using Voronoi voids was found to be 1-2 nm 19 .
We now examine the temperature dependence of the length scales of the considered density variations. It is known that the low-frequency sound velocity is significantly decreased with decreasing temperature 90,98 , while the high-frequency sound velocity is increased 99 . Therefore one can expect that the ratio k 2 /k 1 (and, consequently, δ ) will be increased with lowering temperature. Note that the average correlation length of density fluctuations in water is also increased with decreasing temperature 16 . In the case of supercooled water at 250 K (case III), one has the following data: c(q) ≈ 1.3 × 10 3 m/s at q = 0.023 nm -1 (Ref. 90 ) and c(q) ≈ 3.5 × 10 3 m/s at q = q u = 3 nm -1 (Ref. 99 ). Then one can estimate ω 0 (q) using the relation ω 0 (q) = c(q)q . In order to estimate k 1 and k 2 from these data, we first note that for the case I (the data of Refs. 81,90 ), a least-squares fit to the dispersion curve for ω 0 with only two reference points, q = 0.023 nm −1 and q = q u = 3.0 nm −1 , gives k 1 ≈ 1.8 × 10 6 m 2 /s 2 and k 2 ≈ 7.3 × 10 -13 m 4 /s 2 . These values are comparable to those obtained in Eq. (2a). Therefore we use this fit with the indicated above two reference points to estimate the values of k 1 and k 2 for the cases I and III by one and the same procedure. Then we obtain, in the very first approximation, k 1 ≈ 1.8 × 10 6 m 2 /s 2 , k 2 ≈ 1.2 × 10 -12 m 4 /s 2 , and δ ≈ 1.6 nm for supercooled water at 250 K (case III), and δ ≈ 1.3 nm for water at 278 K (case I). The ratio of these values of δ , approximately equal to 1.3, is consistent with the ratio of the average correlation lengths of density fluctuations for 250 K and 277 K ( ≈ 1.2) 16 . Thus it seems that the temperature-dependent length scales for density fluctuations in water and the stationary soliton solution to our model are in agreement with each other.
(17a) www.nature.com/scientificreports/ The soliton ansatz (16) suggests the existence of propagating nonlinear eigenmodes of the system. In order to consider this possibility, we choose 21 A/ρ 0 = 0.043, in accordance with the estimation based on the root mean square deviation of the number of molecules in a volume of water of size 1 nm 3 . This corresponds to the characteristic amplitude of thermal fluctuations. Then from Eqs. (16a) and (16b), we have the following estimates: δ ≈ 2.1 nm (the FWHM of a pulse is ≈ 3.7 nm) and v ≈ 1.9 × 10 3 m/s in the case I; δ ≈ 2.4 nm (the FWHM of a pulse is ≈ 4.2 nm) and v ≈ 1.6 × 10 3 m/s in the case II. Our numerical simulations showed that on the time scale of 0.5 ps indicated in the beginning of this section, the exact soliton solution (16) can be considered as an approximation to the numerical solution of Eqs. (14) and (14a).
We now consider the case of biomolecules hydration water. Within the low-wave number region q ≤ 4 nm −1 , the damping factors for hydration water of DNA and some proteins, such as the ribonuclease protein, are less than those for bulk water (in contrast to the high-wave number region) 68,100 . For example, for coherent density fluctuations propagating through DNA hydration water, 2Ŵ(q)/ω 0 (q) is about 0.38 for q = 3 nm −1 and 0.32 for q = 2 nm −1 , at the temperature 300 K and at the hydration level corresponding to DNA molecules surrounded by approximately two hydration shells 68 . In this case, since approximately ω 0 (q) ∝ q and Ŵ(q) ∝ q 2 (Refs. 68,100 ), one can anticipate that the ratio Ŵ(q)/ω 0 (q) will be further decreased with decreasing q . Note that both the dispersion curves and the damping factors of the ribonuclease protein and DNA hydration water are very similar to each other in this wave number region 68,101 . Thus one can assume for this characteristic case As a prototypical case of DNA/protein hydration water, we consider the case when ν ′ L = 0.32 and other parameters remain the same as earlier. Our numerical simulations showed that the exact solution (16) can be regarded as an approximation to the numerical solution of Eqs. (14) and (14a) on approximately two times longer time scale in comparison with the case of bulk water ( ν ′ L = 0.60). We note that for hydration water of DNA and proteins, the high-frequency sound velocity is typically larger than that for bulk water 68,100 . For example, the propagation velocity of density fluctuations within DNA hydration water is of about 3500 m/s, which is significantly higher than for bulk water (3040 m/s) 68 . In this case, the ratio k 2 /k 1 would be somewhat higher, and consequently, δ would be larger.
To summarize, the results obtained predict the existence of stable, nearly stationary LDL-like states in ambient and supercooled water on the nanometer length scale and on the time scale up to ~ 0.5-1 ps. These states are associated with nonlinear eigenmodes of the system. Thus, the fully dynamic mechanism coming from the combination of the structural relaxation process and the finiteness of the amplitudes of terahertz longitudinal acoustic-like modes could contribute to the formation and dynamics of fluctuating nanosized density heterogeneities. The obtained results also suggest that for smaller amplitudes of density variations, determined by the average isothermal compressibility 21 , a perturbed soliton-like regime is also possible in liquid water on the time scale of ~ 0.5-1 ps.
The soliton-like collective modes proposed by the present analysis may be connected, in particular, with dynamical heterogeneity in supercooled water. Indeed, spatially separated, extended low-mobility low-density regions and high-mobility high-density regions were observed in supercooled water in recent simulations 61,62 . Subsequent simulations 67 revealed that ice nucleation occurs in the low-density regions. Thus one can suggest that the soliton-like collective excitations might be associated with rare, collective rearrangements in supercooled water from which ice is born 67,102 . The considered soliton-like collective modes can also be expected to be associated with some ultrafast biological processes such as those involving the energy transfer. For instance, they may be connected with protein-induced long-range perturbation of water dynamics extending up to ~ 2 nm or more from the protein surface 64,69 , discussed in the introductory part of the paper. Such long-range perturbation of the solvent's collective dynamics was shown to be able to assist 64,103,104 and even provide a necessary condition for activity of some proteins 103 . The dynamical heterogeneous distribution of water molecules in the vicinity of the interface of a biomolecule may also be important in such biological processes as protein-DNA binding which may occur in positions of local peaks of water density 68,105,106 . Thus, the mechanism suggested by the present analysis may support the corresponding long-range interactions 71,107 .

Methods
For numerical simulations, we rewrite Eq. (15) as two first-order equations: where and h ′ is related to the fluctuating force term ∂ 2 δS ′ /∂x ′2 in Eq. (15). Equations (18) with non-reflecting boundary conditions were solved numerically using a variant of the MacCormack method 108 with spatial step x ′ = 10 -2 and time step t ′ = 2.5 × 10 -5 . The length of the spatial lattice was chosen to be 120. We assume that at the grid node x ′ i , the following relations are satisfied:

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.