text:
# Dynamics of liquid 4He in confined geometries from Time-Dependent Density Functional calculations

## I Introduction

Density Functional (DF) methods evans  have become increasingly popular in recent years as a useful computational tool to study the properties of classical and quantum inhomogeneous fluids, especially for large systems where the DF methods provides a good compromise between accuracy and computational cost. In particular, a quite accurate description of the $`T=0`$ properties of liquid 4He has been obtained within a DF approach by using the energy functional proposed in Ref.dupont  and later improved in Ref.prica . In general, the Density Functionals for 4He depend on a number of phenomenological parameters which are adjusted to reproduce bulk experimental properties of liquid 4He at saturated vapor pressure. The resulting functionals are found to describe accurately also the properties of the inhomogeneous systems (e.g. the liquid-vapor interface of the free 4He surface) (for a thorough comparison between these and other functionals used to describe liquid 4He at $`T=0`$, see Ref.szybis ).

This phenomenological functional has been widely used in a variety of problems involving inhomogeneous 4He systems like pure and doped clusters dalfovo2 ; BarHer ; dalfovo3 , electron ”bubbles” in bulk 4Heanci , alkali atom adsorption on the surface of liquid 4He anci1 , vortices in 4He clustersbarranco1 , adsorption of liquid 4He on solid surfacessartori ; faccin , etc.

According to Ref.prica  the static properties of liquid 4He at T=0 are described by the following energy density functional (atomic units will be used throughout our paper):

$$E_0[\rho ]=\frac{1}{2M}𝑑𝐫\left(\sqrt{\rho }\right)^2+\frac{1}{2}𝑑𝐫𝑑𝐫^{}\rho (𝐫)\rho (𝐫^{})V_{HeHe}(|𝐫𝐫^{}|)+𝑑𝐫E_c(𝐫)$$
(1)
The first term on the right hand side is the quantum pressure (corresponding to the kinetic energy of a Bose gas of nonuniform density). The second term contains a two-body pair potential $`V_{HeHe}`$ screened at short distances, while the last term accounts for correlation effects due to the short-range part of the He-He interaction. The explicit form of the functional is given in Ref.prica , to which the reader is referred for more details (see also the Appendix in the present paper).

The minimization of the energy functional (1) with respect to density variations, subject to the constraint of a constant number of 4He atoms, $`\rho (𝐫)𝑑𝐫=N`$, leads to the equilibrium profile $`\rho (𝐫)`$, thus allowing to study the static properties of the 4He system.

In order to study the dynamics of 4He an additional term (which has no effect on the static properties derived from the minimization of the energy functional (1)) has been added prica  to the functional (1) in the form of a phenomenological current term:

$$H_J=\frac{\rho (𝐫)}{2}M𝐯(𝐫)^2\frac{M}{4}𝑑𝐫^{}V_J(|𝐫𝐫^{}|)\rho (𝐫)\rho (𝐫^{})[𝐯(𝐫)𝐯(𝐫^{})]^2$$
(2)
The first term is the usual hydrodynamic current density, while the second term accounts in a phenomenological way for non-local effects due to the ”backflow” current densityprica . The resulting DF (the so-called Orsay-Trento Functional), which will be used in our calculations, has thus the following form:

$$E[\rho ,𝐯]=E_0[\rho ]+𝑑𝐫H_J$$
(3)
One appealing feature of this functional, which will prove to be essential in the time-dependent calculations presented in the following Sections, is that it reproduces not only the static structure factor but also the bulk dispersion relations of sound excitations in liquid 4He, i.e. the effective current-current interaction $`V_J`$ in (2) is empirically adjusted in order to reproduce the phonon-roton dispersion curve of bulk 4He and the dynamic structure function as well.

The extended functional (3) has been applied in the past to the study of density oscillations in 4He clusterscasas , of surface excitations in 4He films prica ; dalfovo2  and of other related dynamical phenomena dalfovo . In all of the above applications a linear response approximation was used, as well as suitable symmetries were assumed in order to make the calculations feasible.

There are many different situations, however, where a more general approach to investigate the microscopic dynamics of liquid 4He in arbitrary environments, often characterized by geometrical confinement, is highly desirable. For instance, modifications of the excitation spectra due to geometrical confinement are expected in the case of 4He films adsorbed on solid surfaces. A second example involves impurity-doped 4He clusters. Atomic impurities are currently used as experimental probes of the superconducting behavior of 4He clusters: the coupling of the impurity motion with the cluster dynamics is crucial for a comprehension of the observed spectra. As a third example, we mention the problem of the dynamics of wetting phenomena: weakly attractive surfaces like the heavy alkali metal surfaces are not wet by liquid 4He at very low temperature, and thus finite amounts of 4He on these surfaces will form a droplet characterized by a non-zero contact angle. Even on a highly uniform surface, however, the contact angle between the 4He droplet and the surface is extremely hysteretic and its value strongly depends on whether the contact line is advancing or receding. Informations from microscopic calculations of the dynamic behavior of a liquid droplet/film spreading on a weak solid surface can represent a valuable help to understand these effects.

We present here calculations, done in the framework of Time-Dependent Density Functional (TD-DFT) phenomenological theory, of the dynamical properties of 4He (at T=0) in confined geometries. Our methods allows to calculate the direct real-time evolution of a 4He system in arbitrary $`3D`$ geometries. Some applications of the method, showing its efficiency and capabilities, will be presented in the following Sections.

## II Time-Dependent Density Functional Theory

In order to extend the DFT scheme described in the previous Section to the domain of time evolution (Time Dependent Density Functional Theory, TD-DFT), one can use a variational principle as follows. One starts from a suitable action integral $`A`$ written in terms of a complex ”wave function” $`\mathrm{\Psi }(𝐫,t)`$

$$A[\mathrm{\Psi }]=𝑑t𝑑𝐫\left\{[\mathrm{\Psi },\mathrm{\Psi }^{}]i\mathrm{\Psi }^{}\frac{\mathrm{\Psi }}{t}\right\}$$
(4)
where the energy density $``$ is defined by $`E[\rho ,𝐯]=𝑑𝐫`$. From the condition of stationarity

$$\frac{\delta A}{\delta \mathrm{\Psi }^{}}=0$$
(5)
a time-dependent Euler-Lagrange equation follows:

$$\left[\frac{^2}{2M}+U[\rho ,𝐯]\right]\mathrm{\Psi }H\mathrm{\Psi }=i\frac{\mathrm{\Psi }}{t}$$
(6)
The ”effective potential” $`U`$ appearing in the above non-linear time-dependent Schrodinger equation is defined as the variational derivative of the energy functional, $`U[\rho ,𝐯]\delta E[\rho (𝐫)]/\delta \mathrm{\Psi }^{}`$ and its explicit expression is given in the Appendix. From the solution $`\mathrm{\Psi }\varphi e^{i\mathrm{\Theta }}`$ of Eq.(6) one can get the time-dependent density $`\rho (𝐫,t)=\varphi ^2`$ and the fluid velocity field $`𝐯(𝐫,t)=\frac{1}{M}\mathrm{\Theta }`$.

We describe in the following Section the numerical procedure that we used in order to integrate numerically the Eq.(6)nota .

## III Computational scheme

The numerical solution of eq.(6), i.e. the wave function at an arbitrary time $`t`$, $`\mathrm{\Psi }(𝐫,t)`$, is obtained by using the Crank-Nicholson’s scheme (CN), which enables to evolve $`\mathrm{\Psi }(t)`$ to a later time $`t+\mathrm{\Delta }t`$:

$$\left[1+i\frac{\mathrm{\Delta }t}{2}H(t+\mathrm{\Delta }t)\right]\mathrm{\Psi }(t+\mathrm{\Delta }t)=\left[1i\frac{\mathrm{\Delta }t}{2}H(t)\right]\mathrm{\Psi }(t)$$
(7)
The CN algorithm is one of the many approximate methods that can be used to integrate numerically a time-dependent Schrodinger equation like the one shown in Eq.(6). It is $`2^{nd}`$ order accurate in $`\mathrm{\Delta }t`$, as can be easily seen by expanding the two sides of Eq.(7) in powers of $`\mathrm{\Delta }t`$. Among the advantages of such scheme is its numerical stability for long-time evolutions, and the unitariness (i.e. the number of particles is conserved during the time evolution) AC . The CN recursion can be easily solved iteratively for each time step evolution. We find that three iterations are usually enough to guarantee a stable evolution even during long simulation times.

The application of the CN formula requires the calculation of the action of the ”Hamiltonian” $`H\frac{1}{2M}^2+U[\rho ,𝐯]`$ to the complex wavefunction $`\mathrm{\Psi }`$. A direct evaluation of the potential term $`U`$ (whose explicit expression is reported in the Appendix) would be prohibitively costly, due to the presence of convolution integrals of the form $`I(𝐫)=𝑑𝐫^{}f(𝐫𝐫^{})g(𝐫^{})`$, whose evaluation in real space requires a huge number of floating point operations. For such reason we decided to employ a mixed real-reciprocal space representation of the basic quantities $`\mathrm{\Psi }`$ and of the density $`\rho `$. Since the Fourier transform of a convolution is the product of the Fourier transforms of the two terms entering the integrals, one can perform efficiently the calculations of a convolution integral on a computer by using Fast Fourier Transforms (FFT). The use of FFT implies Periodic Boundary Conditions (PBC), and thus our calculations have been performed using a periodically repeated supercell containing $`N`$ Helium atoms. Both $`\mathrm{\Psi }`$ and $`\rho `$ are expanded in plane waves: $`\mathrm{\Psi }(𝐫)=_𝐆\mathrm{\Psi }_𝐆e^{i𝐆𝐫}`$, and similarly for $`\rho `$. The $`𝐆`$’s are the supercell wave-vectors, $`𝐆=\pi (n/L_x,m/L_y,p/L_z)`$, with $`n,m,p=0,\pm 1,\pm 2,\mathrm{}`$ and $`L_x,L_y,L_z`$ are the sides of the supercell (a primitive orthorhombic supercell is used in all our calculations). The number of plane waves in the above expansions is chosen such as to give converged values for the total energy $`E`$ and for the structural parameters of the 4He system under investigation. In the various applications that we report in the following, the supercell size was always chosen in such a way to make negligible the interaction of the 4He system with its repeated images.

The action of the kinetic energy operator $`\widehat{T}\frac{1}{2M}^2`$ on the wavefunction $`\mathrm{\Psi }`$ is also evaluated in reciprocal space, i.e.

$$\widehat{T}\mathrm{\Psi }=\frac{1}{2M}FFT\{G^2\mathrm{\Psi }_G\}$$
(8)
Finally, the real-space expression of $`H\mathrm{\Psi }`$, computed as described above, is used in the iteration formula of the CN scheme to give the evolved wave-function $`\mathrm{\Psi }(t+\mathrm{\Delta }t)`$. A new potential $`U`$ is thus computed from the updated density and velocity, and the procedure is repeated for the next time-step. Typical time steps $`\mathrm{\Delta }t`$ which are found to give stable time evolutions are of the order of 1-3 femtoseconds.

We stress the fact that no geometric or simmetry constraints are imposed during the time evolution to reduce the computational effort. Our computational scheme is thus general and can be applied in particular to those situations characterized by geometrical confinement of 4He. We will give in the following Sections few examples of such case.

## IV 4He bulk excitations

One important feature of the most recent implementation of DF for 4He used here is the inclusion of a ”backflow” current term (the second term in Equation (2)) explicitly depending on the fluid velocity as well as on the fluid density, such that the resulting functional is able to correctly reproduce the bulk excitation spectrum of liquid 4He. This can be immediately verified within our computational scheme as follows.

We setup an initial state corresponding to a modulation of the bulk uniform density with a wave vector $`𝐪`$:

$$\mathrm{\Psi }(𝐫,t=0)=\left\{\rho _0[1+ϵsin(𝐪𝐫)]\right\}^{1/2}$$
(9)
where $`ϵ`$ is a small amplitude and $`\rho _0=0.02184\AA ^3`$ is the saturation density at zero temperature and pressure. The system is then allowed to evolve in time according to the CN scheme described above. By monitoring the periodic oscillations of $`\rho (𝐫,t)`$ during the simulation, we can immediately calculate the mode frequency:

$$\omega =\frac{2\pi }{T}=\omega (𝐪).$$
(10)
Here $`T`$ is the period of the observed oscillation. By varying the wave vector $`𝐪`$ we get the whole dispersion curve of 4He, which we compare in Figure 1 with the experimental data. As it appears from Fig.1, the agreement with the measured data is very good in a wide range of values of the wavevector $`q`$.

We finally note that if the bulk dispersion relation is calculated by neglecting the ”backflow” term in eq.(2), results deviating by as much as 50$`\%`$ from the experimental curva are obtained for $`q>0.5\AA ^1`$.

## V Dynamics of 4He clusters

The dynamics of liquid 4He clusters have been studied theoretically in recent years by means of a number of different methods and thus represents a convenient benchmark for our numerical scheme.

We have considered the monopole and quadrupole vibrational modes of cluster of different sizes, characterized by the angular momentum quantum number $`l=0`$ and $`l=2`$, respectively. The dipole mode ($`l=1`$) corresponds to a translation of the cluster as a whole, and should thus lie at zero energy. To study the monopole mode, which is a radial breathing oscillation of the cluster density, we compute the main frequencies by monitoring the periodic variations with time of the mean square radius of the cluster:

$$R(t)=\sqrt{\frac{1}{N}r^2\rho (𝐫,t)𝑑𝐫}.$$
(11)
To study the quadrupole oscillations we analyze instead the time evolution of the quadrupole moment tensor $`Q_{ij}`$:

$$Q_{ij}(t)=\rho (𝐫,t)[3r_ir_jr^2\delta _{ij}]d^3r.$$
(12)
Suitable initial states for both modes are constructed by applying to the ground state wave function $`\mathrm{\Psi }_0=\sqrt{\rho _{eq}}`$ (where $`\rho _{eq}`$ represents the equilibrium density profile for the cluster, calculated as described in Section 1 by direct minimization of the energy functional (1)) an excitation operator $`e^{i\xi F}`$, where $`F=r^2`$ to excite monopole oscillations, while $`F=r^2Y_{20}`$ to excite quadrupole oscillations:

$$\mathrm{\Psi }(𝐫,t=0)=e^{i\xi F}\mathrm{\Psi }_0(𝐫).$$
(13)
This initial state is thus allowed to evolve according to the equation (6). From a frequency analysis of the calculated $`M(t)`$ and of $`Q_{ij}(t)`$ we then calculate both the amplitudes and the frequencies of the observed oscillations.

We report in Fig.2 the calculated frequencies of the main oscillations observed in the discrete region of the spectrum (i.e. below the evaporation threshold $`\omega <|\mu |`$, $`\mu `$ being the chemical potential of a single 4He atom).

For comparison, we show in the same Figure the results taken from Ref.CaStr , obtained within the Random Phase Approximation (RPA) and using a zero-range Density Functional to describe the 4He system. The approach used in Ref.CaStr , which neglects completely the short-range He-He correlation effects, is accurate only when applied to the study of excitations involving variations in space larger than the average distance between the 4He atoms: this is not, however, a major limitation for the long-wavelength modes like the ones investigated here (external field proportional to $`r^lY_{l0}(\widehat{𝐫})`$.

In the same plot we compare also our calculated frequencies with those obtained from more recent calculations done by using a finite-range DF including non-local effects, very similar to the one used here, and again within RPABarHer . The results are not much different from those obtained using a zero-range functionalCaStr , showing that the effects of non-localities are not important in this case due to the $`q=0`$ character of the multipole excitations studied here.

Our calculated values, as it appears from Fig.2, do not differ much from the two previous theoretical results, even if they are obtained by explicitly including the ”backflow” term (see Eq.(2)) in the 4He density functional. The reason is that this term becomes important in the region where the bulk dispersion relation deviates from its phonon-like behavior, whereas in the present case only low-q excitations are considered. In the next sections we will present cases where the inclusion of the backflow term is instead crucial in order to get accurate results.

## VI Excitations of a 4He monolayer adsorbed on a weakly attractive surface

The study of the structure, growth and excitations of liquid 4He films adsorbed on solid surfaces is an important field of current research on superfluid properties. Reliable theories of inhomogeneous liquid 4He may represent a useful complementary tool for studying the properties of liquid 4He interacting with a solid surfaceKrot  .

Because of the extremely weak He-He interaction, liquid 4He in contact with almost any substrate spreads to form a continuous film over the surface, so that vapor and substrate are never in contact. A remarkable exception to this behavior is found when 4He is adsorbed on heavy alkali metal surfaces: due to the large electron spill-out at these surfaces, alkali metals provide the weakest adsorption potentials in nature for He atoms. Based on this observation, it was suggested cheng2 ; saam  that 4He might not wet some heavy alkali metal surfaces. Subsequent experiments have confirmed this remarkable prediction rutledge ; nacher ; ketola ; klier .

We study here the excitations of a single monolayer of liquid 4He, adsorbed on a weakly attractive surface. We choose here the surface of Rb: this is because Rb is the weakest surface that it is actually wet by 4He at $`T0`$, i.e. the stable state of liquid 4He on Rb is represented by a thin liquid film covering the whole surface, rather than being represented by a droplet as in the case of Cs. The extremely weak He-surface interaction on the other hand is expected to act as a small perturbation on the film. Thus the excitations studied here should be representative to some extent of those for an isolated 4He film.

To model the interaction with a Rb substrate, we use a binding potential $`V_s(z)`$ which describes the interaction between Rb, occupying the half space $`z0`$, and one 4He atom located at a distance $`z`$ above the surface, which is taken to be ideally flat. The explicit form for $`V_s(z)`$, originally proposed in Ref.chizme , has been later revised in Ref.faccin  in order to correctly reproduce the experimental wetting properties of the 4He/Rb system. We do not include any corrugation of the surface on the atomic scale to mimic its actual microscopic structure. This is a good approximation for the case of 4He adsorption on alkali metal surfaces, because the experiments diehl  indeed indicate that a surprisingly smooth surface is seen by adsorbed 4He atoms.

The total energy functional of liquid 4He interacting with a Rb surface is thus described by the energy functional (3) augmented by the additional term $`𝑑𝐫\rho (𝐫)V_s(𝐫)`$, which describes the interaction energy due to the presence of the surface. We have preliminarly studied the equilibrium density of liquid 4He on the Rb surface, at different areal coverages. Our results are reported in Fig.3, where the density profiles in the direction perpendicular to the surface (which is located at $`z=0`$) are shown. Note that as the coverage increases, a second layer starts to form above the first layer. Quite arbitrarily, we take as representative of a ”monolayer” the profile shown in Fig.3 with a thicker line, corresponding to $`n=4.2\times 10^2\AA ^2`$.

We thus prepare our monolayer in a suitable non-equilibrium initial state, characterized by a density modulation with a wave-vector $`𝐪`$ parallel to the surface. We note that, due to the presence of a symmetry-breaking external potential (the He-surface potential) which depends on $`z`$ (the coordinate normal to the surface), such initial state will be immediately coupled, during its time-evolution, to excitations of the film polarized perpendicularly to the surface plane.

From a Fourier analysis of the density variations during the time evolution from this initial state, we are able to get the main excitations of the 4He film. Fig.4 shows our results. The filled and empty squares show the main features observed in the frequency spectrum of the time-dependent density $`\rho (𝐫,t)`$. The filled squares indicate the most intense peaks, whereas the open squares indicate the somewhat weaker features in the spectrum.

For comparison we also show with solid line the experimental ”ripplon” dispersion relationlauter1 . It appears from Fig.4 that at low $`q`$ values the lower disperion curve becomes softer that the ripplon curve. This is not surprising since the ripplon dispersion, which relates to the surface excitations of a infinite half-space occupied by the liquid, is altered by the finite thickness of the film. A simple model calculation kala  gives in fact for the oscillation frequency of a thin 4He film of thickness $`a`$ in the low-$`q`$ limit:

$$\omega (q)\omega _{ripplon}(q)(qa)^{1/2}$$
(14)
Since $`\omega _{ripplon}q^{3/2}`$, the expected dispersion scales as $`q^2`$, in qualitative agreement with our findings. Although a definite character cannot be assigned unambiguously to the observed mode shown in Fig.4 due to the strong mixing between the two bands, we suggest that the upper branch should have a strong component of density modulations perpendicular to the surface. To verify this we have induced in the monolayer a simple $`q=0`$ mode where the 4He film is initially rigidly shifted (by a very small amount) towards the surface. Due to the presence of the holding He-surface potential, the film once let free to evolve, execute a motion where it executes oscillations mainly perpendicul to the surface, with a main frequency at about 5 K: this value is shown with a filled dot at $`q=0`$ in Fig.4.

We can compare our results for the film excitations with those obtained in Ref.clements1  by means of a microscopic statistical theory, where the dynamic excitations of thick 4He films on a Cs surface were considered. In that work it is found that ripplon modes appear in the calculated spectra, whose dispersion follows the $`q^{3/2}`$ law for surface excitations, as well as ”interfacial” ripplons localized mainly at the 4He/surface interface. In particular, the lowest mode calculated in Ref.clements1  (corresponding the the lower branch in Fig.4) has a ripplon character. The second mode has a ”longitudinal” character (i.e. parallel to the surface) on stronger substrates, only on very weak surfaces (i.e. Cs) it has a visible perpendicular component. A linear (third sound) dispersion has been also observed in Ref.clements1  as long as the wavelength is large compared with the film thickness. We do not observe any such modes, probably because of the finite thickness of our system.

We believe that the upper branch in Fig.4 is the counterpart of a similar excitations observed in neutron diffraction experiments on 4He adsorbed on a graphite substratelauter , where, among other excitations, almost dispersionless modes at low frequencies have been observed. These mode have been interpreted Apaja , by means of microscopic calculations, as standing waves of the 4He liquid perpendicular to the substrate. We cannot compare directly our results with the experimental measurements of lauter  and with the microscopic calculations of Ref.Apaja  because of the different substrate used (graphite is a much ”stronger” substrate than alkali metal surfaces) and because multylayer films were considered in these references. We hope that our results will stimulate further experimental measurements, in particular to study the excitations of 4He monolayers on weakly attractive surfaces.

## VII Dynamics of alkali-doped 4He clusters

The spectroscopy of doped liquid 4He clusters has become a current tool for the study of superfluidity in 4He dropletstoennies . As a consequence of the ultra-weak alkali-He interaction, an alkali atom picked up by a 4He cluster has its stable state in a ”dimple” on the surface of the clusteranci1 ; ernst1 , rather than being solvated into its interior like most of the other atomic impurities. Thus alkali could provide an excellent probe for the surface excitation of superfluid nanodropletsernst1 .

The analysis and interpretation of experimental studies of the formation of 4He-alkali complexes in doped clusters would benefit from theoretical work of the excitations of impurities attached to 4He clusters.

Although for light alkali the motion of the impurity can be considered to a first approximation as occurring in a static 4He environment (i.e. He does not adjust appreciably to the istantaneous atom position), however, for heavy alkali (Rb,Cs) a full dynamical approach is required, where the 4He atoms in the cluster are allowed to adjust dynamically upon the motion of the adsorbed impurity.

We address here the problem of finding the spectrum of vibrational excitations of a ”dimple” state. We consider in particular the system composed of a single Rb atom (in its electronic ground-state) attached to a 300-atom 4He cluster. The Rb-He interaction is taken from Ref.patil .

Fig. 5 shows the calculated equilibrium density profile (in a plane containing the center of the cluster and the Rb atom) for the ”dimple” state representing the lowest energy configuration of the Rb-4He system.

By applying to the Rb atom a small initial momentum towards the surface of the cluster, we observe that the impurity starts oscillating around its equilibrium position shown in Fig. 5, and at the same time the 4He density changes in time to adjust dynamically to the istantaneous Rb position. Our results for the Rb dynamics are shown in Fig.6. In the upper panel we show the calculated istantaneous position of the Rb atom as a function of time (measured with respect to the Rb-cluster center-of-mass). A straightforward Fourier transform of the signal shown in the upper panel allows to compute the frequency spectrum shown in the lower panel.

We compare our results with a previous estimate ernst2  based on static DF calculations, where the frequency of the oscillations of a Rb atom attached to a 300-atom cluster was estimated approximately by assuming that the He atoms adiabatically adjusted to the istantaneous position of the impurity. In that case the value of the frequency was extracted from the shape of the impurity-He potential energy surface, giving $`\omega 0.7K`$. It appears from our results that a more complex dynamics results, due to the coupling of the Rb atom motion with the surface oscillations of the cluster. The frequency of a surface mode of an N-atom 4He cluster can be written as $`\omega =\sqrt{l(l1)(l+2)}\omega _0/\sqrt{N}`$ where $`\omega _0=3.45K`$. This is approximately $`0.6K`$ for the lowest $`l=2`$ mode, i.e. close to the natural Rb frequency calculated in Ref.ernst2 . A large amount of mixing is thus expected, which determines the appearance of the complex vibrational spectrum shown in Fig.6.

Experimental measurements of such oscillations are under way ernst\_private  to confirm our predictions. A more complex and challenging application of time-dependent methods, i.e. the study of the oscillations of an alkali impurity in an excited electronic state, is currently in progress.

## VIII Spreading of a 4He nanodroplet on a weakly attractive substrate

The mechanism underlying the spreading of a liquid droplet on a solid substrate is a long-standing problem. Experiments done on a number of differente (classical) molecular fluids have revealed universal laws such as the Tanner’ spreading law tanner , where the droplet radius grows as $`r(t)t^{1/10}`$, as well as fascinating phenomena on the atomistic scale (such as the appearence of a monolayer film advancing in front of the dropheslot , two different regimes of the growth of the radius with timeconinck1 , dynamical ”layering”, etc.). Other experimentsheslot  have shown that one or more monomolecular precursor layers spread ahead of the droplet cap with an average radius $`R(t)\sqrt{t}`$.

Numerical studies of droplet spreading have been performed blake ; koplik ; nieminen ; coninck1  in recent years. Despite many efforts, however, many microscopic details of the spreading process are however still to be known.

Even more challenging is the behavior of a superfluid 4He droplet on a weakly attractive surface: it has been found that 4He droplets adsorbed on a Cs substrate have spreading and flow properties that are not simple consequences of bulk superfluid behaviorscience . When a 4He droplet is deposited on a Cs surface (the only material known that is not wetted by superfluid 4He at very low temperature), even on a highly uniform surface the contact angle between the droplet and the surface is extremely hysteretic and its value strongly depends on whether the contact line is advancing or receding. Superfluid 4He droplets on Cs are also remarkable because they can resist flow against a substantial chemical potential gradientscience .

We have studied the spreading of a 200-atom He cluster on a weakly attractive surface. We choose for the latter the case of Cs (but very similar results are obtained for a Rb surface). To simulate the Cs surface, we use the ”ab initio” potential developed in Ref.chizme , which has applied successfully to the wetting properties of 4He an a Cs surface faccin . The initial state is represented by a spherical cluster (whose density profile is obtained by minimization of the static functional (2) in the absence of the attractive He-surface potential), placed in the vicinity of the surface. The downward force responsible for the spreading of the droplet is provided solely by the long-range Van der Waals forces exerted by the substrate (due to the microscopic size of the droplet, the effect of gravity on spreading is negligible).

The spreading will continue, until the stable state is eventually reached. The nature of the latter depends on the strength of the He-surface potential: in the case of Cs, which is not wet by He at T=0, the final configuration will be that of an almost spherical cap characterized by a contact angle of $`40^{}`$ (partial wetting behavior)science ; sartori . More attractive substrate, on the other hand, will eventually be covered by a uniform microscopic film (complete wetting behavior). Due to the finiteness of our systems and our use of periodic boundary conditions our simulations are necessarily interrupted when the droplet hits its repeated images and starts coalescing with them.

The sequence characterized the spreading process is shown in Fig.7. The overall duration of the simulation shown in the Figure is about $`2ns`$. Interestingly enough, a precursor layer seems to form under the cap of the spreading cluster. We have computed from the cluster density profile $`\rho (𝐫,t)`$ the average radius $`R`$ of this precursor layer. The dependence of $`R`$ from the time, shown in Fig.8, shows a ”fast” regime where the precursor layer expands linearly in time, as found in similar simulations for classical fluidsnieminen . Then a slowing down occurs, in correspondence of the formation of the second layer (see panel (f) in Fig.7). It seems that after a transient, an almost linear spreading occurs again, although with a somewhat lower velocity. The speed from the linear portion of Fig.8 gives a velocity of about $`50m/s`$, which is close to the critical Landau velocity. The slowing down is perhaps due to the spontaneous formation of vortical rings during the spreading. One such ring is shown in Fig.9 by means of a vector plot showing the current density $`𝐉\rho (𝐫)𝐯(𝐫)`$ in a plane containing the center of the cluster and perpendicular to the Cs surface. Fig.9 also clearly shows the flow of atoms from the droplet cap feeding the underlying precursor film.

Although transient precursor films spreading linearly with time has been obtained from MD simulations of classical fluidsnieminen ; blake , they have not been seen in any experiment. This is probably due to the fact that they represent the very early stage of the spreading process: as the equilibrium shape of the droplet is reached, the initial fast spreading is found to slow down to a power law $`R(t)t^{1/7}`$ruijter , when a regime dominated by dissipation due to the viscous flow starts to dominate. It would be extremely interesting to study how this slowing down (which is only suggested from our short-time results in Fig.8) occurs in the case of superfluid flow.

The presence of surface disorder is expected to alter the dynamics of liquid 4He films spreading on a solid surface, because the three-phase contact line of the spreading film/droplet can easily be pinned to these defects, modifying the velocity of the contact line as it advances (or recedes) on the surface. We are currently investigating the effect of simple defects on the atomic-scale (adatoms or vacancies) on the droplet spreading process discussed above.

## IX Summary

We have presented TD-DFT calculations of the dynamical properties of 4He in situations characterized by geometrical confinement. We have shown the feasibility and efficiency of fully 3-D TD-DFT calculations to study the dynamics of liquid He-4 in a number of representative situations: the dynamics of pure and alkali-doped 4He clusters, the excitations of a liquid 4He monolayer adsorbed on a weakly attractive surface and the dynamics of a 4He nano-droplet spreading on a Cs surface.

Acknowledgments: We thank M. Pi, M. Barranco, E. Krotscheck, W. Ernst, F. Dalfovo, V.A. Apkarian and M.W. Cole for useful comments and discussions. We thank V.A. Apkarian for giving us a copy of the manuscript Ara  prior to publication. We acknowledge funding from MIUR-COFIN 2001.

Appendix

The effective potential entering Eq.(6) can be readily evaluated by functional differentiation of the energy functional (3) (see Ref.prica  for its detailed expression):

$`U[\rho ,𝐯]`$ $`=`$ $`V_s(𝐫)+{\displaystyle 𝑑𝐫^{}\rho (𝐫^{})V_{\mathrm{}}(|𝐫𝐫^{}|)}+{\displaystyle \frac{c_2}{2}}\overline{\rho }(𝐫)^2+{\displaystyle \frac{c_3}{3}}\overline{\rho }(𝐫)^3+{\displaystyle 𝑑𝐫^{}\rho (𝐫^{})[c_2\mathrm{\Pi }_h(|𝐫𝐫^{}|)\overline{\rho }(𝐫^{})+c_3\mathrm{\Pi }_h(|𝐫𝐫^{}|)\overline{\rho }^2(𝐫^{})]}`$ (15)
$`+{\displaystyle \frac{\alpha _s}{2M}}\left(1{\displaystyle \frac{\rho (𝐫)}{\rho _0}}\right){\displaystyle 𝑑𝐫^{}\left(1\frac{\rho (𝐫^{})}{\rho _0}\right)_r^{}\rho (𝐫^{})_rF(|𝐫^{}𝐫|)}{\displaystyle \frac{M}{2}}{\displaystyle 𝑑𝐫^{}V_J(|𝐫𝐫^{}|)\rho (𝐫^{},t)(v(𝐫)v(𝐫^{}))^2}`$
$`+{\displaystyle \frac{i}{2\rho (𝐫)}}{\displaystyle 𝑑𝐫^{}𝐠(𝐫,𝐫^{})}`$
$`V_s(𝐫)`$ represents an additional external potential acting on the 4He system. The second term contains a two-body 4He-4He pair potential $`V_{\mathrm{}}(r)`$ screened at distances shorter than a characteristic length $`h_{\mathrm{}}`$, while the third and the fourth term (correlation terms), which contains the average of the density over a sphere of radius $`h_{\mathrm{}}`$, $`\overline{\rho }_r`$, accounts for the internal kinetic energy and for the increasing contribution of the hard-core He-He repulsion when the density is increased. The last two terms represent the contribution to the potential $`U`$ of the ”backflow” term (2) introduced in Ref.prica  in order to reproduce the bulk excitation spectrum of liquid 4He. In particular, the last term contains the ”backflow” current $`𝐉_B𝑑𝐫^{}𝐠(𝐫,𝐫^{})=𝑑𝐫^{}V_J(|𝐫𝐫^{}|)\rho (𝐫)\rho (𝐫^{})(𝐯(𝐫)𝐯(𝐫^{}))`$

The free parameters $`h_{\mathrm{}}`$, $`c_2`$, $`c_3`$, $`\alpha _s`$ are adjusted in order to reproduce a number of experimental properties of bulk liquid 4He. For a detailed description of the various terms and the numerical values of the adjustable parameters, we refer the interested reader to Ref.prica .

Figure Captions

Figure 1: Bulk 4He excitation spectrum. The points show the calculated spectrum, while the solid line shows the experimental dispersion relation.

Figure 2: Calculated monopole (upper lines) and quadrupole (lower lines) oscillation frequencies of liquid 4He clusters, shown as a function of the cluster sizes. Dots: this work; squares: results from Ref.CaStr ; triangles: results from Ref.BarHer .

Figure 3: Equilibrium density profiles for 4He films adsorbed on a Rb surface. The z direction is perpendicular to the surface plane. Different curves represent increasing values of the surface coverages (number of 4He atoms per unit surface area), from $`n=1.6\times 10^2\AA ^2`$ (lower profile) to $`n=7.2\times 10^2\AA ^2`$ (upper profile). The thicker line (defining our ”monolayer”) is drawn at $`n=4.2\times 10^2\AA ^2`$.

Figure 4: Calculated dispersion curves for a 4He monolayer adsorbed on a Rb surface (squares). The dotted line is the bulk dispersion relation, while the solid line shows the experimental ”ripplon” dispersion.

Figure 5: Density contour plot for a 300-atom 4He cluster with a Rb atom attached on its surface: the dot shows the equilibrium position of the Rb impurity.

Figure 6: Upper panel: Rb displacement with respect to the center-of mass of the Rb-cluster system, measured along the radial direction, as a function of the simulation time. Lower panel: Fourier spectrum (in arbitrary units) of the time series shown in the upper panel.

Figure 7: Sequence of 4He density contour plots showing the spreading of a 200-atom cluster on a Cs surface.

Figure 8: Radius of the precursor layer as a function of time.

Figure 9: Vector plot of the current density $`\rho 𝐯`$ in a plane perpendicular to the surface and passing through the center of the cluster. The configuration shown corresponds to the density contours in panel (e) of Fig.7.


format: text

text:
# Electronic structure of overstretched DNA

## Abstract

Minuscule molecular forces can transform DNA into a structure that is elongated by more than half its original length. We demonstrate that this pronounced conformational transition is of relevance to ongoing experimental and theoretical efforts to characterize the conducting properties of DNA wires. We present quantum mechanical calculations for acidic, dry, poly(CG)$``$poly(CG) DNA which has undergone elongation of up to 90 % relative to its natural length, along with a method for visualizing the effects of stretching on the electronic eigenstates. We find that overstretching leads to a drastic drop of the hopping matrix elements between localized occupied electronic states suggesting a dramatic decrease in the conductivity through holes.

The possibility of efficient charge transfer along the stacked $`\pi `$ orbitals in DNA was proposed 40 years ago ELEY1961 . In recent years its conducting properties were tested in several contexts: DNA has been used as a template on top of which nanowires are built BRAUN1998  and there is hope that self-assembled DNA junctions between conductors will function soon. There has also been renewed interest in building and characterizing molecular wires of plain DNA between two metal leads FINK1999 ; PORATH2000 ; DEPABLO2000 ; KASUMOV2001 ; YOO2001 , possibly combined with changes of the major constituents of the molecule to improve its properties RAKITIN2001 . There is ongoing debate as to whether DNA is conducting, semi-conducting, insulating or all of the above FINK2001 . There are several reasons for this controversy. There is great diversity of the DNA forms in terms of its composition, length, and structure. Counterions and impurities that can attach to the phosphate groups and the grooves of the molecule need to be taken into account. Finally, environmental factors such as temperature and contact resistance can play a major role.

Experiments done 50 years ago by Wilkins et al. WILKINS1951 , suggested that overstretched DNA (that is, substantially longer than its natural length) undergoes a transition to a structure that can accommodate elongation up to twice the length of relaxed DNA. A major breakthrough in understanding the mechanical properties of DNA was achieved through single molecule stretching experiments SMITH1996 ; CLUZEL1996 ; STRICK2000 . It was shown that by gradually increasing the axial stretching force $`F`$ along the DNA helix, the molecule first uncoils, then, at $`F`$ 5 pN, it reaches its natural length, and then exhibits a stiff elastic response for forces of up to 50 pN. At $`F`$ 50 to 80 pN, DNA undergoes a pronounced and abrupt structural transformation to a yet unknown structure that is elongated by more than 50%. Further increase in force leads to breaking after it has stretched to twice its natural length. In this letter we argue that these overstretched forms of DNA are relevant in the quest for novel device components, we develop a visualization method for understanding the nature of electronic states in these forms, and we present the first detailed discussion of their electronic properties.

We use an efficient quantum mechanical electronic structure method ELSTNER1998a  that treats all the valence electrons and that accounts for the charge transfer involved in biological molecules and semiconducting surfaces. This approach, called self-consistent-charge density-functional tight-binding method, is based on an expansion of the Kohn-Sham energy functional HOHENBERG1964 ; KOHN1965  in terms of the charge density. The zeroth order term results in a tight-binding Hamiltonian, and the second order term incorporates the charge redistribution. The method is able to deal with the very large numbers of atoms involved in biologically relevant systems. As a test, we performed calculations on a periodic homopolymer of poly(C)$``$poly(G) with 10 base pairs per turn and reproduced the results of dePablo et al. DEPABLO2000 , which to our knowledge is the only ab initio electronic structure calculation of a periodic DNA structure LADIK . Our results are in excellent agreement to their calculations: we find a 2.1 eV bandgap between the highest occupied and lowest unoccupied eigenstates, compared to 2.0 eV reported in Ref. \[DEPABLO2000, \]. The wavefunctions of the highest occupied valence band are located on the guanines while the wavefunctions of the lowest unoccupied band are located on the cytosines, also in agreement with Ref. \[DEPABLO2000, \].

For the purpose of the present study, calculations were performed on a variety of DNA polymers with different compositions employing both infinite periodic and finite structures. We limit our discussion to poly(CG)$``$poly(CG) as a representative structure for which well established overstretched configurations have been published LEBRUN1996 . Our results on finite DNA polymers show the existence of states localized at the opposite ends of the molecule. The energy levels of these states lie between the regular bands that appear in the periodic structures. We removed these states before analyzing the electronic properties, and will refer to the remaining states as the “restricted band structure” of the non-periodic configuration. As a test, we compared 18 base pair restricted band sturcture with periodic structures with 10 base pairs per unit cell and found the band structures virtually indistinguishable. A boundary region of 3 base-pair layers at each end of the molecule was sufficient to contain the spurious end states. Fig. 1 displays the band structure of a periodic sequence with a unit cell consisting of 10 base pairs. The sharp features of the density of states are due to the one dimensional nature of DNA. In all periodic sequences there is an extra symmetry, which for the poly(CG)$``$poly(CG) structure is the translation by two base pairs along the stack, together with a rotation of 72 degrees around the helical axis: this symmetry is generated by the screw operator that is isomorphic to a pure translation operator. Structures obtained by short (50 ps) molecular dynamics runs at 50 K using classical force fields CORNELL1995 , indicate that minor deformations of the ideal configuration destroy the sharp features in the electronic structure.

We now move on to a discussion of overstretched DNA. First we elaborate on how this is relevant in the construction of DNA nanowires. One way to attach DNA molecules to substrates is the technique of “molecular combing” BENSIMON1994 . This method uses the meniscus forces that develop between the solution surface and the device surface template to uncoil the DNA and stretch it between two electrodes as the device is pulled out of the solution. A description of recent experiments estimates that the forces that are exerted on the DNA molecules are of order 100 pN KLEIN2001 , more than enough to stretch them past their natural length. This is indicated by recently built devices WU2001  where the distance between the electrodes was about 30% longer than the natural length of the DNA molecules. The overstretching can have a variety of consequences on the electronic states.

Elongated DNA structures have been determined in the pioneering study of Lebrun and Lavery LEBRUN1996 , which model the adiabatic elongation of selected DNA molecules. There are two distinct modes of overstretching corresponding to pulling the opposite 3’-3’ ends and 5’-5’ ends. The former mode leads to DNA unwinding, consistent with a recent estimate of the helicity of overstretched DNA LEGER1999 , whereas the latter leads to a contraction of the diameter of the helix. Both overstretching modes can accommodate elongations of up to 90% without breaking the molecule.

We next present the inherent electronic states of the overstretched structures and then consider electron transport aspects. For all the calculations we used acidic, dry, 18 base-pair strings of DNA, and we performed the projection of the end states in the manner stated before. The quasi one-dimensional structure of all the forms leads to the formation of clear minibands in the energy spectrum, similar to those shown in Fig. 1. The band structure of the overstretched forms have an overall resemblance to the band structure of the original unstretched form, and many of the minibands near the Fermi energy can be continuously mapped to their counterpart in the unstretched configuration. This is mostly because of the molecular structure of the overstretched forms as will become clear in the following discussion. For some of the most deformed structures, especially the 90% overstretching in the 5’-5’ mode, several of the minibands are mixed and small but clear gaps present in less deformed structures disappear.

We have devised a new scheme for visualizing the electronic states of DNA, which naturally accounts for the symmetry of this molecule. The basic concept of this scheme is to unfold contour plots of the band electron density evaluated at cylindrical surfaces centered on the helical axis, as shown in Fig. 2. In this figure we show the structure at its natural length together with one that is elongated by 30% in the 3’-3’ overstretching mode. With the help of these plots we can understand most of the band structure features of the stretched forms. The elongation to the overstretched form is achieved by changing the dihedral angle configuration of the backbone of the molecule, thus leaving the local part of the orbitals essentially intact.

In Fig. 3 (a) we present the valence band electron density of selected stretched poly(CG)$``$poly(CG) structures using our visualization scheme. We choose the radii of the cylinders to follow the atom with the largest electron density for the band that we present. In the unstretched structure, the $`\pi `$-band is along the guanines. Upon elongation, orbitals from the picture frame are removed while the remaining ones rotate partially as dictated by the new orientation of base pairs. The poly(CG)$``$poly(CG) structure contains pairs of guanines that are closer to each other. These guanines display an overlap of the electronic densities as seen in the top left frame of Fig. 3 (b). Elongation in the 3’-3’ mode keeps such paired guanines close to each other which is reflected in the overlap in their electronic densities; the overlap of the electronic densities between guanines from different pairs decreases with elongation. In Fig. 3 (b) we show the conduction band electron density. Note again how the orbitals rotate as the structure is being overstretched, following the rotation of the bases. Note that in the extreme 90% stretching case, the orbitals become perpendicular to the helical axis.

The changes in the structural features should be of profound consequence to any model dealing with electron transfer since, to a first approximation, the overlap matrix elements depend exponentially on the distance between orbitals. Because of the $`\pi `$ character of the states near the Fermi level, a non-parallel alignment will strongly affect the values of the overlap integrals. In order to quantify these effects, we considered the nearly degenerate eigenstates in the first valence miniband. We used the linear combination of these states that gives normalized, orthogonal states $`|\mathrm{\Psi }_i`$, which are maximally localized along the helical axis within the subspace spanned by the delocalized states. The resulting valence states are localized on guanine bases. With these states, we calculated the hopping matrix elements:

$$t_{ij}=\frac{1}{2}\mathrm{\Psi }_i|^2|\mathrm{\Psi }_j,$$
(1)
between nearest-neighbor sites. Note that for the poly(CG)$``$poly(CG) structure every guanine has 2 guanine neighbors, one above and one below as indicated in Fig. 3, separated by distances of 5.3 Å and 4.2 Å. The larger of the two matrix elements corresponds to hopping between sites that are close together; its value remains roughly constant for the 30$`\%`$ overstretched structures. The smaller of the two matrix elements, which corresponds to sites that are further apart, is dramatically reduced for all the overstretched forms. Since the two different hopping terms are connected in a 1-D series the smallest term will determine the bottleneck for electron transport. The absolute square values of these small hopping matrix elements between these neighbors are given in Fig. 4 (a). We infer from this figure that if the conduction mechanism is through holes there should be a dramatic drop in conductivity with 30% overstretching, which in the 3’-3’ mode drops further with elongation. In Fig. 4 (b) we give the hopping matrix elements between states in the conduction band. In the cases where the pairing between nearest neighbor pairs have very different matrix elements we only list the one that would be the bottleneck for this channel. Due to the stronger localization of the conduction band and because of the rather large distance between the occupied sites, the hopping matrix elements are substantially smaller than the ones between the valence states in the structure of natural length. In fact, it is the site that is two bases away which has the strongest overlap in this structure: the C6 atoms between these bases are actually closer to each other (9.0 Å) than the C6 atoms between consecutive base pairs (10.0 versus 10.3 Å). The first and second neighbor hopping terms are indicated in Fig. 3 (b). Since these terms are connected in parallel, this time it is the maximum of the two that will play the dominant role in the conduction mechanism and it is this one that is included in Fig. 4 (b). We note that if conduction takes place by electrons the conductivity drops with overstretching in the 3’-3’ mode but not necessarily in the 5’-5’ mode.

The concept of overstretched DNA in relation to electronic device components appears relevant to a variety of experimental and theoretical issues. Experimentally it is possible, and sometimes even unavoidable, to create overstretched DNA molecules. The window of allowable forces that can be applied to stretch DNA to its natural length, but not overstretch, is small (from 5 pN to 50 pN) compared to the forces applied during molecular combing. Another important aspect is owed to the nature of the extreme stretched forms, it seems possible to expose the bases to the exterior of the molecule allowing impurities to enter the crucial base-pair stack region; this could result in controlling the electronic properties with doping. This might also enable direct visualization of parts of the electronic wavefunctions which was is not possible in the natural DNA forms because of the insulating backbone layer. On the theoretical level, the major existing models (superexchange JORTNER1998 , hopping JORTNER1998 , and polaron HENDERSON1999 ), of pure DNA conductivity will be drastically affected by the consequences of overstretching, including changes in the hopping matrix elements and in the low-energy vibrational modes. For long fibers of homogeneous DNA (either a homopolymer, or a short repeated sequence) we expect stretching to appear homogeneously. Regions of varying richness in AT or GC mixture would induce different lengthenings at different parts of the molecule, permitting conduction mechanisms to change along the molecule.

We thank Leo Kouwenhoven and Nina Markovic for sharing their experimental insight with us. We are grateful to Richard Lavery for providing the overstretched structures. We thank Hatem Mehrez for his critical reading of the manuscript. This work was supported in part by Harvard’s Materials Research Science and Engineering Center, which is funded by the National Science Foundation. This work was partially supported by National Computational Science Alliance under MCB010005N and utilized the Boston University Scientific Computing and Visualization facilities. RLB was supported by a national science foundation graduate research fellowship.


text:
# 1 Introduction

## 1 Introduction

We describe a new cosmological model based on matching a critically expanding FRW metric to a metric which we call the Tolman-Oppenheimer-Volkoff (TOV) metric inside the Black Hole, across a shock wave that lies beyond one Hubble length from the center of the FRW spacetime. This implies that the spacetime beyond the shock wave must lie inside a Black Hole, thus extending the shock matching limit of one Hubble length that we identified in .

In the exact solutions constructed in this paper, the expanding FRW universe emerges behind a subluminous blast wave that explodes outward from the origin $`\overline{r}=0`$ at the instant of the Big Bang $`t=0,`$ at a distance beyond one Hubble length333We let $`(t,r)`$ denote standard FRW coordinates, so that $`\overline{r}=rR\left(t\right)`$ measures arclength distance at each fixed value of the FRW time $`t,`$ where $`R`$ denotes the cosmological scale factor. Barred coordinates also refer to TOV standard coordinates, in which case $`\overline{r}=rR\left(t\right)`$ also holds as a consequence of shock matching.. The shock wave then continues to weaken as it expands outward until the Hubble length eventually catches up to the shock, and this marks the event horizon of a Black Hole in the TOV metric beyond the shock. From this time onward, the shock wave is approximated by a zero pressure ($`k=0`$) Oppenheimer-Snyder (OS) interface, and thus the OS solution gives the large time asymptotics of these solutions. Surprisingly, the equation of state $`p=\frac{1}{3}\rho `$ of early Big Bang physics, is distinguished by the differential equations, and only for this equation of state does the shock wave emerge from the Big Bang at a finite nonzero speed, the speed of light. This is surprising because the equation of state $`p=\frac{1}{3}\rho `$ played no special role in shock matching outside the Black Hole, . We find it interesting that such a shock wave emerging from the Big Bang beyond the Hubble length, would thermalize the radiation in a region well beyond the light cone of an observer positioned at the FRW center, even though the model does not invoke inflation. Details will appear in our forthcoming paper ; we wish here to summarize this work, and describe its physical interpretation.

## 2 Statement of the Problem

If there is a shock wave at the leading edge of the expansion of the galaxies, then we can ask what is the critical radius $`\overline{r}_{crit}`$ at each fixed time $`t`$ in a $`k=0`$ FRW metric such that the total mass inside a shock wave positioned beyond that radius puts the universe inside a Black Hole? (There must be such a critical radius because the total mass $`M(\overline{r},t)`$ inside radius $`\overline{r}`$ in the FRW metric at fixed time $`t`$ increases like $`\overline{r}^3,`$ and so at each fixed time $`t`$ we must have $`\overline{r}>2M(\overline{r},t)`$ for small enough $`\overline{r},`$ while the reverse inequality holds for large $`\overline{r}.`$ We let $`\overline{r}=\overline{r}_{crit}`$ denote the smallest radius at which $`\overline{r}_{crit}=2M(\overline{r}_{crit},t).`$) We show that when $`k=0,`$ $`\overline{r}_{crit}`$ equals the Hubble length. Thus we cannot match a critically expanding FRW metric to a classical TOV metric beyond one Hubble length without continuing the TOV solution into a Black Hole—and we showed in  that the standard TOV metric cannot be continued into a Black Hole. Thus to do shock matching with a $`k=0`$ FRW metric beyond one Hubble length, we introduce the TOV metric inside the Black Hole.

## 3 The TOV Metric Inside the Black Hole

When the metric anzatz is taken to be of the TOV form

$`ds^2=B\left(\overline{r}\right)d\overline{t}^2+A^1\left(\overline{r}\right)d\overline{r}^2+\overline{r}^2d\mathrm{\Omega }^2,`$ (3.1)
and the stress tensor $`T`$ is taken to be that of a perfect fluid co-moving with the metric, the Einstein equations $`G=\kappa T,`$ inside the Black Hole, take the form

$`\overline{p}^{}`$ $`=`$ $`{\displaystyle \frac{\overline{p}+\overline{\rho }}{2}}{\displaystyle \frac{N^{}}{N1}},`$ (3.2)
$`N^{}`$ $`=`$ $`\left\{{\displaystyle \frac{N}{\overline{r}}}+\kappa \overline{p}\overline{r}\right\},`$ (3.3)
$`{\displaystyle \frac{B^{}}{B}}`$ $`=`$ $`{\displaystyle \frac{1}{N1}}\left\{{\displaystyle \frac{N}{\overline{r}}}+\kappa \overline{\rho }\right\}.`$ (3.4)
We let $`\overline{\rho },`$ $`\overline{p}`$ denote the density and pressure, and $`\overline{r}`$ is taken to be the timelike variable because we assume

$`A\left(\overline{r}\right)=1{\displaystyle \frac{2M\left(\overline{r}\right)}{\overline{r}}}1N\left(\overline{r}\right)<0.`$ (3.5)
Here, $`M\left(\overline{r}\right)`$ has the interpretation as the total mass inside the ball of radius $`\overline{r}`$ when $`\overline{r}>2M,`$ but $`M`$ does not have the same interpretation inside the Black Hole444System (3.2)-(3.4) for $`A<0`$ differs substantially from the TOV equations for $`A>0`$ because, for example, the energy density $`T^{00}`$ is equated with the timelike component $`G^{rr}`$ when $`A<0,`$ but with $`G^{tt}`$ when $`A>0.`$ In particular, this implies that $`M^{}=4\pi \overline{p}\overline{r}^2`$ when $`A<0,`$ versus $`M^{}=4\pi \overline{\rho }\overline{r}^2`$ when $`A>0,`$ the latter being the equation that gives the mass function its physical interpretation outside the Black Hole. because $`\overline{r}<2M.`$  The system (3.2)-(3.4) defines the simplest class of gravitational metrics that contain matter, and evolve inside the Black Hole.

## 4 Shock Matching Inside the Black Hole

Matching a given $`k=0`$ FRW metric to a TOV metric inside the Black Hole across a shock interface, leads to the system of ODE’s

$`{\displaystyle \frac{du}{dN}}=\left\{{\displaystyle \frac{\left(1+u\right)}{2\left(1+3u\right)N}}\right\}\left\{{\displaystyle \frac{\left(3u1\right)\left(\sigma u\right)N+6u\left(1+u\right)}{\left(\sigma u\right)N+\left(1+u\right)}}\right\},`$ (4.1)
$`{\displaystyle \frac{d\overline{r}}{dN}}`$ $`=`$ $`{\displaystyle \frac{1}{1+3u}}{\displaystyle \frac{\overline{r}}{N}},`$ (4.2)
with conservation constraint

$`v={\displaystyle \frac{\sigma \left(1+u\right)+\left(\sigma u\right)N}{\left(1+u\right)+\left(\sigma u\right)N}}.`$ (4.3)
Here

$`u={\displaystyle \frac{\overline{p}}{\rho }},v={\displaystyle \frac{\overline{\rho }}{\rho }},\sigma ={\displaystyle \frac{p}{\rho }},`$ (4.4)
$`\rho `$ and $`p`$ denote the (known) FRW density and pressure, $`\sigma =p/\rho ,`$ and all variables are evaluated at the shock. Solutions of (4.1)-(4.3) determine the (unknown) TOV metrics that match the given FRW metric Lipschitz continuously across a shock interface, such that conservation of energy and momemtum hold across the shock, and such that there are no delta function sources at the shock, .

For such solutions, the speed of the shock interface relative to the fluid comoving on the FRW side of the shock, is given by

$`s=R\dot{r}=\sqrt{N}\left({\displaystyle \frac{\sigma u}{1+u}}\right).`$ (4.5)
Note that the dependence of (4.1)-(4.3) on the FRW metric is only through the variable $`\sigma .`$ Since solutions of (4.1)-(4.3) are formally time-reversible but shock waves are not, system (4.1)-(4.3) must be augmented by entropy condition for shocks that breaks the time symmetry. Since we are interested in solutions that model the “Big Bang” as a localized explosion with an outgoing blast wave emanating from $`\overline{r}=0`$ at time $`t=0,`$ we impose the entropy conditions,

$`0<\overline{p}<p,0<\overline{\rho }<\rho .`$ (4.6)
Condition (4.6) for outgoing shock waves implies that the shock wave is compressive, and is sufficient to rule out unstable rarefaction shocks in classical gas dynamics.

## 5 Exact Shock Wave Solutions Inside the Black Hole

In the case when the FRW pressure is given by the equation of state

$`p=\sigma \rho ,`$ (5.1)
$`\sigma `$ assumed constant, $`0<\sigma <1,`$ the FRW equations have the exact solutions

$`\rho ={\displaystyle \frac{4}{3\kappa \left(1+\sigma \right)^2}}{\displaystyle \frac{1}{t^2}},`$ (5.2)
$`R=\left({\displaystyle \frac{t}{t_0}}\right)^{\frac{2}{3\left(1+\sigma \right)}},`$ (5.3)
which assumes an expanding universe, ($`\dot{R}>0`$), and initial conditions $`R=0`$ at $`t=0,`$ and $`R=1,`$ at $`t=t_0.`$ Since $`\sigma `$ is constant, equation (4.1) uncouples from (4.2), and thus solutions of system (4.1)-(4.3) are determined by the scalar non-autonomous equation (4.1). Making the change of variable $`S=1/N,`$ which transforms the “Big Bang” $`N\mathrm{}`$ over to rest point at $`S0,`$ c.f. , we obtain,

$`{\displaystyle \frac{du}{dS}}=\left\{{\displaystyle \frac{\left(1+u\right)}{2\left(1+3u\right)S}}\right\}\left\{{\displaystyle \frac{\left(3u1\right)\left(\sigma u\right)+6u\left(1+u\right)S}{\left(\sigma u\right)+\left(1+u\right)S}}\right\},`$ (5.4)
We take as entropy condition (4.6), and in addition, we require that the TOV equation of state meet the physical bound

$`0<\overline{p}<\overline{\rho }.`$ (5.5)
Note that the conditions $`N>1`$ and $`0<\overline{p}<p`$ restrict the domain of (5.4) to the region $`0<u<\sigma <1,`$ $`0<S<1.`$ Inequalities (4.6) and (5.5) are both implied by the single condition

$`S<\left({\displaystyle \frac{1u}{1+u}}\right)\left({\displaystyle \frac{\sigma u}{\sigma +u}}\right).`$ (5.6)
We prove the following theorem:

###### Theorem 1

For every $`\sigma ,`$ $`0<\sigma <1,`$ there exists a unique solution $`u_\sigma (S)`$ of (5.4), such that (5.6) holds on the solution for all $`S,`$ $`0<S<1,`$ and on this solution, $`0<u_\sigma (S)<\overline{u},`$ $`lim_{S0}u_\sigma (S)=\overline{u},`$ where $`\overline{u}=Min\{1/3,\sigma \}.`$ and

$`\underset{S1}{lim}\overline{p}=0=\underset{S1}{lim}\overline{\rho }.`$ (5.7)
Concerning the the shock speed, we have:

###### Theorem 2

Let $`0<\sigma <1.`$ Then the shock wave is everywhere subluminous, that is, the shock speed $`s_\sigma (S)s(u_\sigma (S))<1`$ for all $`0<S1,`$ if and only if $`\sigma 1/3.`$

By a careful analysis of the asymptotics of the solution near $`S=0,`$ we can prove

###### Theorem 3

The shock speed at the Big Bang $`S=0`$ is given by:

$`\underset{S0}{lim}s_\sigma \left(S\right)=0,\sigma <1/3,`$ (5.8)
$`\underset{S0}{lim}s_\sigma \left(S\right)=\mathrm{},\sigma >1/3,`$ (5.9)
$`\underset{S0}{lim}s_\sigma \left(S\right)=1,\sigma =1/3,`$ (5.10)
## 6 Bounds on the Shock Position

Let $`r_{}`$ be the FRW radial position of the shock wave at the instant of the Big Bang, (the arclength distance $`\overline{r}_{}=r_{}R\left(0\right)=0`$ when $`R\left(0\right)=0`$). The analysis implies that the shock wave will first become visible at the FRW center $`\overline{r}=0`$ at the moment $`t=t_0,`$ ($`R\left(t_0\right)=1`$), when the Hubble length $`H_0^1=H^1\left(t_0\right)`$ satisfies

$`{\displaystyle \frac{1}{H_0}}={\displaystyle \frac{1+3\sigma }{2}}r_{},`$ (6.1)
where $`r_{}`$ is the FRW position of the shock at the instant of the Big Bang. At this time, we prove that the number of Hubble lengths $`\sqrt{N}_0`$ from the FRW center to the shock wave at time $`t=t_0`$ can be estimated by

$`1{\displaystyle \frac{2}{1+3\sigma }}\sqrt{N}_0{\displaystyle \frac{2}{1+3\sigma }}e^{\sqrt{3\sigma }\left(\frac{1+3\sigma }{1+\sigma }\right)}.`$
Thus, in particular, the shock wave will still lie beyond the Hubble length $`1/H_0`$ at the FRW time $`t_0`$ when it first becomes visible. Furthermore, we prove that the time $`t_{crit}>t_0`$ at which the shock wave will emerge from the White Hole given that $`t_0`$ is the first instant at which the shock becomes visible at the FRW center, can be estimated by

$`{\displaystyle \frac{2}{1+3\sigma }}e^{\frac{1}{4}\sigma }{\displaystyle \frac{t_{crit}}{t_0}}{\displaystyle \frac{2}{1+3\sigma }}e^{\frac{2\sqrt{3\sigma }}{1+\sigma }},`$ (6.2)
for $`0<\sigma 1/3,`$ and by the better estimate

$`e^{\frac{\sqrt{6}}{4}}{\displaystyle \frac{t_{crit}}{t_0}}e^{\frac{3}{2}},`$ (6.3)
in the case $`\sigma =1/3.`$ For example, (6.2), (6.3) imply that at the OS limit $`\sigma =0,`$

$$\sqrt{N_0}=2,\frac{t_{crit}}{t_0}=2,$$

and in the limit $`\sigma =1/3,`$

$$1.8\frac{t_{crit}}{t_0}4.5,1<\sqrt{N_0}4.5,$$

We conclude in these shock wave cosmological models, that the moment $`t_0`$ when the shock wave first becomes visible at the FRW center, it must lie within 4.5 Hubble lengths of the center. Throughout the expansion up until this time, the expanding universe must lie entirely within a Black Hole—the universe will eventually emerge from this Black Hole, but not until some later time $`t_{crit},`$ where $`t_{crit}`$ does not exceed $`4.5t_0.`$

## 7 Concluding Remarks

We have constructed global exact solutions of the Einstein equations in which the expanding FRW universe extends out to a shock wave that lies arbitrarily far beyond the Hubble length. The distance to the shock wave at any given value of the Hubble constant is determined by one free paramenter which can be taken to be the FRW position of the shock wave at the instant of the Big Bang.

The critical OS solution inside the Black Hole is obtained in the limit of zero pressure, and provides the large time asymptotics, but the shock wave solutions differ qualitatively from the OS solution. For example, the shock wave models contain matter, and thus do not rely on any portion of the empty space Schwarzschild metric inside the Black Hole for their construction. In both models, the interface propagates outward from the FRW center $`\overline{r}=0`$ at the instant of the Big Bang, but in the shock wave model the mass function $`M(\overline{r},t)`$ is infinite at the instant of the Big Bang, and immediately becomes a finite decreasing function of FRW time, for all future times $`t>0.`$ And although the OS solution is time reversible, the directional orientation of the shock interface relative to the various observers is determined by an entropy condition, . The entropy condition selects the explosion over the implosion, and the condition that the entropy condition be satisfied globally, determines a unique solution. Since the TOV radial coordinate $`\overline{r}`$ is timelike inside the Black Hole, we can also say that the density $`\overline{\rho }\left(\overline{r}\right)`$ and mass $`M\left(\overline{r}\right)`$ are both constant at each fixed time in the TOV spacetime beyond the shock wave.

The shock interface that marks the boundary of the FRW expansion continues out through the White Hole event horizon of an ambient Schwarzschild metric at the instant when the wave is exactly one Hubble length from the FRW center $`\overline{r}=0,`$ and then continues on out to infinity along a geodesic of the Schwarzschild metric outside the Black Hole. These solutions thus indicate a scenario for the Big Bang in which the expanding universe emerges from an explosion emanating from a White Hole singularity inside the event horizon of an asymptotically flat Schwarzschild spacetime of finite mass. The model does not require the physically implausible assumption that the uniformly expanding portion of the universe is of infinite mass and extent at every fixed time, and it has the nice feature that it embeds the Big Bang singularity of cosmology within the event horizon of a larger spacetime, the Schwarzschild spacetime. Moreover, the model also allows for the uniform expansion of arbitrarily large densities within an arbitrarily large mass extended over an arbitrary number of Hubble lengths early on in the Big Bang, a prerequisite for the standard physics of the Big Bang at early times.

We conclude that these shock wave solutions give the global dynamics of strong gravitational fields in an exact solution, the dynamics is qualitatively different from the dynamics of solutions when the pressure $`p0,`$ the solution suggests a Big Bang cosmological model in which the expanding universe is bounded throughout its expansion, and the equation of state most relevant at the Big Bang, $`p=\frac{1}{3}\rho ,`$ is distinguished by the differential equations. But these solutions are only rough qualitative models because the equation of state on the TOV side is determined by the equations, and therefore cannot be imposed. That is, the TOV density $`\overline{\rho }`$ and pressure $`\overline{p}`$ only satisfy the entropy conditions (4.6), together with the loose physical bounds (5.5); and on the FRW side, the equation of state is taken to be $`p=\sigma \rho ,`$ $`\sigma const.,`$ $`0<\sigma <1.`$ Nevertheless, these bounds imply that the equations of state are qualitatively reasonable, and we expect that these solutions will capture the gross dynamics of solutions when more general equations of state are imposed. For more general equations of state, other waves, (e.g. rarefaction waves), would need to be present to meet the conservation constraint, and thereby mediate the transition across the shock wave. Such transitional waves would be pretty much impossible to model in an exact solution, but the fact that we can find global solutions that meet our physical bounds, and that are qualitatively the same, for all values of $`\sigma (0,1),`$ and all initial shock positions, leads us to expect that such a shock wave would be the dominant wave in a large class of problems.

As a final remark, we note that because Einstein’s theory by itself does not choose an orientation for time, it follows that if we believe that a Black Hole can exist in the forward time collapse of a mass through an event horizon as $`t\mathrm{},`$ (the time $`t`$ as observed in the far field), then we must also allow for the possibility of the time reversal of this process, a White Hole explosion of matter out through an event horizon coming from $`t\mathrm{}.`$ That is, if we are willing to accept Black Holes in the forward time dynamics of astrophysical objects whose collapse appears so great as to form an event horizon in the future, then by symmetry, we may well also be forced to accept White Holes in the backward time dynamics of astrophysical objects which, like the expanding universe, appear to have expansions so great as to have emerged from an event horizon in the past. Given this, it is natural to wonder if there is a connection between the mass that disappears into Black Hole singularities, and the mass that emerges from White Hole singularities.