# A Superfluid Hydrodynamic Model for the Enhanced Moments of Inertia of Molecules in Liquid 4He

## Abstract

We present a superfluid hydrodynamic model for the increase in moment of inertia, $`\mathrm{\Delta }I`$, of molecules rotating in liquid <sup>4</sup>He. The static inhomogeneous He density around each molecule (calculated using the Orsay–Paris liquid <sup>4</sup>He density functional) is assumed to adiabatically follow the rotation of the molecule. We find that the $`\mathrm{\Delta }I`$ values created by the viscousless and irrotational flow are in good agreement with the observed increases for several molecules \[OCS, (HCN)<sub>2</sub>, HCCCN, and HCCCH<sub>3</sub>\]. For HCN and HCCH, our model substantially overestimates $`\mathrm{\Delta }I`$. This is likely to result from a (partial) breakdown of the adiabatic following approximation.

The spectroscopy of atoms and molecules dissolved in He nanodroplets provides both a new way to study microscopic dynamics of this unique quantum fluid , and a very cold matrix (0.4 K ) to create and study novel species . Recent experiments have demonstrated that even heavy and anisotropic molecules display rotationally resolved vibrational spectra with a structure reflecting the gas phase symmetry of the molecule. However, the rotational constants required to reproduce the spectra are often substantially reduced from those of the isolated molecule. For example, the $`\nu _3`$ vibrational band of SF<sub>6</sub> dissolved in He nanodroplets (first observed by Goyal et al.  and later rotationally resolved and analyzed by Hartmann et al. ) indicates that the effective moment of inertia, $`I_{\mathrm{eff}}`$, in liquid <sup>4</sup>He is 2.8 times that of the isolated molecule. The same qualitative behavior has been found for a wide range of other molecules . In an elegant recent experiment, it has been demonstrated that the rotational structure of OCS broadens and collapses in pure <sup>3</sup>He droplets, and is recovered when $`60`$ <sup>4</sup>He atoms are co-dissolved in the <sup>3</sup>He . The association of the weakly damped, unhindered rotation with the Bose symmetry of <sup>4</sup>He suggests that this phenomenon is a manifestation of superfluidity, and has been called the microscopic Andronikashvili experiment .

A theory able to reproduce the observed increase, $`\mathrm{\Delta }I`$, in molecular moments of inertia would be of interest for at least two reasons. First, the enhanced inertia provides a window into the dynamics of the liquid. Second, the ability to predict the rotational constants would further improve the utility of He nanodroplet isolation spectroscopy for the characterization of novel chemical species.

The first model proposed to explain the observed $`\mathrm{\Delta }I`$ assumed that a certain number of He atoms, trapped in the interaction potential of the solute, rotate rigidly with the latter . In the case of SF<sub>6</sub>, 8 He atoms trapped in the octahedral global potential minima would create a rigidly rotating ‘supermolecule’ that would have approximately the observed $`I_{\mathrm{eff}}`$. In the case of OCS, putting a six He atom ‘donut’ in the potential well around the molecule also reproduces the observed $`I_{\mathrm{eff}}`$ . Recent Diffusion Monte Carlo (DMC) calculations  have predicted that the effective rotational constant of SF<sub>6</sub>–He<sub>N</sub> monotonically decreases from that of the isolated molecule to the large cluster limit, reached at $`N=8`$, and remains essentially constant for $`N=8`$$`20`$. The supermolecule model has been recently extended to consider the rigid rotation of a ‘normal fluid fraction’ of the He density which is claimed to be significant only in the first solvation layer , based on Path Integral Monte Carlo calculations of Kwon et al.  which show a molecule-induced reduction of the superfluid fraction. These calculations have been recently used to propose a definition of a spatially dependent normal fluid fraction which reproduces the observed $`I_{\mathrm{eff}}`$ of solvated SF<sub>6</sub> .

The limitations of the supermolecule model are made clear by the $`\mathrm{\Delta }I`$ observed for HCN in He droplets \[8, a\] which is only $`5\%`$ of the $`\mathrm{\Delta }I`$ observed upon formation of a gas phase He$``$HCN van der Waals complex . Furthermore, it has been previously recognized that in principle there is also a superfluid hydrodynamic contribution, $`I_\mathrm{h}`$, to $`\mathrm{\Delta }I`$ . Previous estimates, based upon a classic treatment of the rotation of an ellipsoid in a fluid of uniform density, found that $`I_\mathrm{h}`$ is only a small fraction of the observed $`\mathrm{\Delta }I`$, at least for heavy rotors such as OCS . In this report, we show that if the spatial variation of the He solvation density around the solute molecule is taken into account, the calculated $`I_\mathrm{h}`$ is instead rather large and agrees well with experimental data. We compare our calculations with the experimental results available in the literature (OCS , HCN \[8, a\]) and with results recently obtained in our laboratory for HCCCH<sub>3</sub> and HCCCN, and in the laboratory of R. E. Miller for (HCN)<sub>2</sub> \[8, b\] and HCCH \[8, c\].

We first calculate the ground state He density, $`\rho `$, around a static solute molecule. The molecule is then considered to undergo classical rotation, slowly enough that the helium ground state density adiabatically follows the molecular rotation. The kinetic energy associated with the He flow (assumed viscousless and irrotational) is used to calculate $`I_\mathrm{h}`$.

The main input of our hydrodynamic model is the ground state density of He around the solute. DMC calculations can provide this density with a minimum of assumptions beyond the interaction potentials , but are computationally expensive. The Density Functional method, which is a good compromise between accuracy and computational cost, consists in numerically minimizing the total energy of the many-body system in the form of a semi-empirical functional of the He density: $`E=𝑑𝐫[\rho (𝐫)]`$. The energy density $``$ contains an effective non-local interaction with a few parameters fixed to reproduce known properties of bulk liquid He. The functional used here is the one termed Orsay-Paris , which was shown to accurately reproduce the static properties of pure and doped He clusters . The need to treat axially symmetric molecules implies moving from one to two-dimensional equations. The new routines have been extensively tested against previously calculated spherically symmetric systems. The minimization of energy is carried out by mapping the density distribution on a grid of points and propagating it in imaginary time, starting from a trial distribution.

The density functional also contains the interaction between the He and the impurity molecule. The interaction, assumed pairwise, is treated as a static external potential, since the molecules considered here are expected to have negligible zero point motion. Existing potentials for He-HCN , He-HCCH , and He-OCS  have been used without modifications. The He-(HCN)<sub>2</sub> potential was generated as the superposition of the potential due to two HCN molecules whose centers of mass are separated by 4.44 $`\mathrm{\AA }`$ (the equilibrium distance for the gas phase dimer ). The repulsive part of the He-HCCCN potential has been taken from ; the attractive part from the He-HCN and He-HCCH potentials, using the concepts of distributed interaction and transferability . The He-HCCH  and He-CH<sub>4</sub>  interactions were used to generate the potential between He and HCCCH<sub>3</sub>, treating the latter molecule as cylindrically symmetric. Full detail on all potentials used are available from the authors, and will be published separately .

Once the helium density profiles are calculated, the molecules are assumed to rotate perpendicularly to their symmetry axis with angular velocity $`\omega `$. We assume that the He density adiabatically follows this rotation, which allows us to calculate the laboratory-frame time-dependent density at each point in the liquid. This assumption is only valid if at each point the velocity of the fluid, $`v(𝐫)`$, is less than a critical velocity, $`v_c`$. If $`v_c`$ is taken to be the velocity of sound, this is true for all our molecules, at the temperature of the droplet: 0.4 K. A further justification to our assumption is also the fact that no critical value of angular momentum is experimentally observed for a wide class of molecules (i.e. for a wide range of fluid velocities).

The second assumption that we make is that the He behaves entirely as a superfluid undergoing irrotational flow. The assumption that the motion is irrotational implies that $`𝐯(𝐫)`$ can be written as the gradient of a scalar potential: $`𝐯=`$$`\mathbf{}`$$`\mathit{\varphi }`$ (the dependence of $`𝝆\mathbf{,}𝐯\mathbf{,}\mathit{\varphi }`$ on $`𝐫`$ will be implicit from now on), where $`\mathit{\varphi }`$ is known as the velocity potential. These assumptions lead to the following hydrodynamic equation for the velocity potential :

$$\mathbf{}\mathbf{}\mathbf{(}𝝆\mathbf{}\mathit{\varphi }\mathbf{)}\mathbf{=}\frac{\mathbf{}𝝆}{\mathbf{}𝒕}\mathbf{=}\mathbf{}\mathbf{(}\mathbf{}𝝆\mathbf{)}\mathbf{}\mathbf{(}𝝎\mathbf{\times }𝐫\mathbf{)}\mathbf{.}$$
(1)   
The first equality is just the continuity equation, while the second reflects the statement that the density is time-independent in the rotating frame. We select our axis system with $`𝒛`$ along the symmetry axis of the molecule, and assume that rotation takes place round the $`𝒙`$ axis with angular velocity $`𝝎\mathbf{=}𝝎\widehat{𝐱}`$. In order to better exploit the symmetry of the problem, we have used elliptical coordinates $`𝝃\mathbf{,}𝜽\mathbf{,}𝝋`$, where $`𝒙\mathbf{=}𝒇\sqrt{𝝃^\mathrm{𝟐}\mathbf{}\mathrm{𝟏}}\mathrm{𝐬𝐢𝐧}\mathbf{(}𝜽\mathbf{)}\mathrm{𝐜𝐨𝐬}\mathbf{(}𝝋\mathbf{)}\mathbf{,}𝒚\mathbf{=}𝒇\sqrt{𝝃^\mathrm{𝟐}\mathbf{}\mathrm{𝟏}}\mathrm{𝐬𝐢𝐧}\mathbf{(}𝜽\mathbf{)}\mathrm{𝐬𝐢𝐧}\mathbf{(}𝝋\mathbf{)}`$, and $`𝒛\mathbf{=}𝒇𝝃\mathrm{𝐜𝐨𝐬}\mathbf{(}𝜽\mathbf{)}`$. The surfaces of constant $`𝝃`$ are ellipses of rotation with foci at $`𝒛\mathbf{=}\mathbf{\pm }𝒇`$. Two such surfaces limit the region where Eq. (1) is solved. The inner boundary excludes the volume occupied by the impurity, and is chosen as the largest ellipse contained in the region where $`𝝆\mathbf{<}\mathbf{0.005}𝝆_\mathrm{𝟎}`$ ($`𝝆_\mathrm{𝟎}\mathbf{=}\mathbf{0.0218}\mathit{\AA }^\mathbf{}\mathrm{𝟑}`$ is the bulk liquid density). Von Neumann boundary conditions $`\widehat{𝒏}\mathbf{}\mathbf{}\mathit{\varphi }\mathbf{=}\mathbf{}\widehat{𝒏}\mathbf{}\mathbf{(}𝝎\mathbf{\times }𝐫\mathbf{)}`$ insure that the normal component of velocity matches the normal component of motion of the boundary . For the outer boundary, any ellipse large enough that the motion of the outside fluid is negligible can be chosen (with Dirichlet boundary conditions $`\mathit{\varphi }\mathbf{=}\mathrm{𝟎}`$). These boundary conditions result in a unique solution to the hydrodynamic equations. Other solutions exist if we do not require the fluid to be irrotational, but it is known that these are higher in energy , and will include any solutions that have some portion (a “normal component” or a He “snowball”) of the He density that rigidly rotates with the molecule.

Given the solution, $`\mathit{\varphi }`$, to the hydrodynamic equation, we can calculate the kinetic energy, $`𝑲_𝐡`$, in the motion of the fluid by the following:

$`𝑲_𝐡`$ $`\mathbf{=}`$ $`{\displaystyle \frac{\mathrm{𝟏}}{\mathrm{𝟐}}}𝑰_𝐡𝝎^\mathrm{𝟐}\mathbf{=}{\displaystyle \frac{\mathrm{𝟏}}{\mathrm{𝟐}}}𝒎_{\mathrm{𝐇𝐞}}{\displaystyle \mathbf{}𝝆\mathbf{(}\mathbf{}\mathit{\varphi }\mathbf{)}\mathbf{}\mathbf{(}\mathbf{}\mathit{\varphi }\mathbf{)}𝒅𝑽}`$ (2)   
$`𝑲_𝐡`$ $`\mathbf{=}`$ $`{\displaystyle \frac{\mathrm{𝟏}}{\mathrm{𝟐}}}𝒎_{\mathrm{𝐇𝐞}}\mathbf{\left[}\mathbf{}{\displaystyle \mathbf{}\mathit{\varphi }\mathbf{\left(}\frac{\mathbf{}𝝆}{\mathbf{}𝒕}\mathbf{\right)}𝒅𝑽}\mathbf{+}{\displaystyle \mathbf{}𝝆\mathit{\varphi }\mathbf{(}\mathbf{}\mathit{\varphi }\mathbf{)}\mathbf{}𝒅𝐒}\mathbf{\right]}\mathbf{.}`$ (3)   
Eq. (2) follows directly from the definition of kinetic energy; Eq. (3) is derived from Eq. (2) using standard vector identities and assuming that $`\mathit{\varphi }`$ is a solution of Eq. (1). $`𝒅𝐒`$ is defined as positive when pointing out of the region of the fluid. $`𝑰_𝐡`$ is the hydrodynamic contribution to the moment of inertia for rotation about the $`𝒙`$ axis, and $`𝒎_{\mathrm{𝐇𝐞}}`$ is the atomic mass of helium. Both $`\mathbf{(}\mathbf{}𝝆\mathbf{/}\mathbf{}𝒕\mathbf{)}`$ and $`\mathit{\varphi }`$, are proportional to $`𝝎`$, thus the above definition of $`𝑰_𝐡`$ is independent of $`𝝎`$. The total kinetic energy of rotation will include the contribution from the molecule, $`𝑲_𝐦\mathbf{=}\frac{\mathrm{𝟏}}{\mathrm{𝟐}}𝑰_𝐦𝝎^\mathrm{𝟐}`$, where $`𝑰_𝐦`$ is the moment of inertia of the free molecule. We can also calculate the net angular momentum created by the motion of the He fluid: $`𝐉_𝐡\mathbf{=}𝒎_{\mathrm{𝐇𝐞}}\mathbf{}𝝆𝐫\mathbf{\times }\mathbf{(}\mathbf{}\mathbf{}\mathit{\varphi }\mathbf{)}𝒅𝑽`$. By use of standard vector identities and Eq. (1), this definition can be shown to lead to $`𝐉_𝐡\mathbf{=}𝑰_𝐡𝝎`$. The total angular momentum is the sum of that of the rotating molecule and the total moment of inetia the sum of the moment of inertia of the molecule and that due to hydrodynamic motion of the superfluid. The local shape of the velocity field $`𝐯\mathbf{(}𝐫\mathbf{)}`$ can be rather complex due to the presence of strong inhomogeneities in the density distribution.

We calculate $`𝑰_𝐡`$ by solving the hydrodynamic equation, Eq. (1) for $`\mathit{\varphi }`$, assuming unit angular velocity rotation around the $`𝒙`$ axis. It is computationally convenient to solve a slightly transformed version of Eq. (1), where the smoother function $`\mathrm{𝐥𝐧}𝝆\mathbf{(}𝐫\mathbf{)}`$ appears instead of $`𝝆\mathbf{(}𝐫\mathbf{)}`$:

$$\mathbf{}^\mathrm{𝟐}\mathit{\varphi }\mathbf{+}\mathbf{\left(}\mathbf{}\mathrm{𝐥𝐧}𝝆\mathbf{\right)}\mathbf{}\mathbf{\left(}\mathbf{}\mathit{\varphi }\mathbf{+}\widehat{𝒙}\mathbf{\times }𝐫\mathbf{\right)}\mathbf{=}\mathrm{𝟎}\mathbf{.}$$
(4)   
Eq. (4) is solved, subject to the boundary conditions, by converting it to a set of finite difference equations on a grid of points in our elliptical coordinate system and using the Gauss-Seidel relaxation method . Both Eq. (2) and Eq. (3) are then evaluated by simple numerical quadrature, and are found to give the same value of $`𝑰_𝐡`$ within a few percent. We also carefully tested the convergence of $`𝑰_𝐡`$ with grid size.

As an example of density distribution and velocity field, Fig. 1 shows our results for the OCS molecule in a cluster of 300 He atoms. On the left we give the contour plot of the He density near the molecule. One clearly sees the complex structure which results from the tendency to have He atoms near local minima of the impurity-He potential. The highest peak, at $`\mathbf{(}𝒚\mathbf{,}𝒛\mathbf{)}\mathbf{=}\mathbf{(}\mathbf{}\mathbf{3.6}\mathbf{,}\mathbf{}\mathbf{1.2}\mathbf{)}`$, corresponds to a ring of atoms perpendicular to the axis of the molecule. The integral of the density within this structure gives 6.5 atoms, and indeed 7–8 is the number of He atoms one expects to fit into such a ring by close-packing. On the right side of the same figure we plot the current density, $`𝝆𝐯`$. We find that most of the kinetic energy density, $`\frac{\mathrm{𝟏}}{\mathrm{𝟐}}𝝆𝐯^\mathrm{𝟐}`$, comes from the first solvation layer, the outer part of the cluster giving a negligible effect.

In Table I our results are compared with existing experimental values for several molecules in He nanodroplets. There is an overall good agreement between the predicted and observed enhancements of the effective moment of inertia. From a quantitative viewpoint, one notices that the predicted moments of inertia tend to overestimate the experimental values. In the case of the lightest rotors (HCN and HCCH) the large discrepancy suggests the breakdown of the assumption of adiabatic following as recently predicted . In that paper the importance of He exchange is pointed out; it is also shown that the interplay of the rotational constant with the potential anisotropy determines the extent to which the anisotropic He solvation density can adiabatically follow the rotation of the molecule. When the rotational constant of SF<sub>6</sub> is arbitrarily increased in the calculation by a factor of 10, the He density in the molecule-fixed frame becomes much more isotropic and the solvation-induced $`𝚫𝑰`$ decreases by a factor of 20 . We have recently obtained experimental evidence that $`𝚫𝑰`$ is larger for DCN than for HCN, which we believe to be direct experimental evidence for this effect . It is interesting to note that for these light (i.e. fast spinning) rotors the maximum of $`𝐯\mathbf{(}𝐫\mathbf{)}`$ approaches the bulk <sup>4</sup>He sound velocity.

The overestimate of the moments of inertia for the other molecules likely reflects the uncertainties in the calculated $`𝝆\mathbf{(}𝐫\mathbf{)}`$. We should remark here that while, by construction, the Orsay-Paris functional prevents $`\overline{𝝆}`$ (the density averaged over an atomic volume) from becoming much larger than $`𝝆_\mathrm{𝟎}`$, the functional was not constructed to deal with density gradients as high as those found in the first solvation layer. We observed that small changes in the form of the He density within the deep potential well of those molecules produce significant variations of the predicted moments of inertia, limiting the accuracy of the final results to 20% – 30%. This uncertainty does not affect the main result emerging from Table I that the hydrodynamic contribution to the moment of inertia of these systems, instead of being negligible, is rather large and can explain the observed rotational constants.

One could object that the density values found at the minima of the He-molecule interaction potential (e.g. $`\mathbf{}\mathrm{𝟏𝟏}𝝆_\mathrm{𝟎}`$ for OCS) are too high to be treated as those of a liquid, and should be interpreted as localized He atoms rigidly rotating with the molecule; it has been proposed that the He density distribution around the OCS-He<sub>6</sub> supermolecule is only weakly anisotropic and thus can rotate without generating a significant hydrodynamic contribution . We have calculated the above density distribution, and found that it is still strongly anisotropic, leading to a hydrodynamic moment of inertia of over 400 $`𝐮\mathbf{}\mathbf{\AA }^\mathrm{𝟐}`$. When combined with the moment of inertia of the OCS-He<sub>6</sub> supermolecule, this gives a total effective moment of inertia of over 650 $`𝐮\mathbf{}\mathbf{\AA }^\mathrm{𝟐}`$, dramatically larger than the experimental value (230 $`𝐮\mathbf{}\mathbf{\AA }^\mathrm{𝟐}`$).

In summary, the spatial dependence of the He density, which is caused by the molecule-He interaction, results in a hydrodynamic contribution to the moment of inertia more than an order of magnitude larger (in the case of the heavier rotors) than that predicted for the rotation of a reasonably sized ellipsoid in He of uniform bulk liquid density. Furthermore, the present calculations suggest that the effective moments of inertia of molecules in He nanodroplets (and likely also bulk He) can be quantitatively predicted by assuming irrotational flow of a spatially inhomogeneous superfluid.

We are pleased to acknowledge useful discussions and/or the sharing of unpublished information with D. Farrelly, Y. Kwon, E. Lee, R. E. Miller, K. Nauta, L. Pitaevskii, and K. B. Whaley. The work was supported by the National Science Foundation.