# Theoretical modeling of electron mobility in superfluid

Frédéric Aitken, Nelly Bonifaci, Klaus von Haeften, Jussi Eloranta

► To cite this version:

Frédéric Aitken, Nelly Bonifaci, Klaus von Haeften, Jussi Eloranta. Theoretical modeling of electron mobility in superfluid. The Journal of Chemical Physics, 2016, 10.1063/1.4959293. hal-01362267

HAL Id: hal-01362267

https://hal.science/hal-01362267

Submitted on 14 Jan 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Theoretical modeling of electron mobility in superfluid $^4$He

Frédéric Aitken,¹ Nelly Bonifaci,¹ Klaus von Haeften,² and Jussi Eloranta³,⁴

¹G2ELab-GreEn-ER, Equipe MDE, 21 Avenue des Martyrs, CS 90624, 38031 Grenoble Cedex 1, France

²Department of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, United Kingdom

³Department of Chemistry and Biochemistry, California State University at Northridge, 18111 Nordhoff St., Northridge, California 91330, USA

The Orsay-Trento bosonic density functional theory model is extended to include dissipation due to the viscous response of superfluid $^4$He present at finite temperatures. The viscous functional is derived from the Navier-Stokes equation by using the Madelung transformation and includes the contribution of interfacial viscous response present at the gas-liquid boundaries. This contribution was obtained by calibrating the model against the experimentally determined electron mobilities from 1.2 K to 2.1 K along the saturated vapor pressure line, where the viscous response is dominated by thermal rotons. The temperature dependence of ion mobility was calculated for several different solvation cavity sizes and the data are rationalized in the context of roton scattering and Stokes limited mobility models. Results are compared to the experimentally observed “exotic ion” data, which provides estimates for the corresponding bubble sizes in the liquid. Possible sources of such ions are briefly discussed.

## I. INTRODUCTION

An excess electron confined in superfluid $^4$He is an excellent prototype system for modeling the general aspects of nonpolar solvation and the interaction of nanoscale objects with liquids. For example, the static solvation structure of electrons in liquid helium has received a considerable amount of experimental and theoretical attention.¹ These studies established that electrons reside in spherical cavities (“bubbles”) in the liquid with a barycenter radius of ca. 18.5 Å at saturated vapor pressure as determined by earlier density functional theory (DFT) calculations.²,³ In contrast, such bubbles are considerably smaller in polar liquids, such as water, due to the strong electrostatic interactions.⁴ Since the bubble size in superfluid helium is determined by a delicate balance between the electron zero-point energy and the liquid surface energy, its radius is strongly dependent on the external pressure. The quantized energy levels of an electron trapped in such a spherical well have been studied experimentally (i.e., optical transitions between $1s-1p$ and $1s-2p$ states) as a function of external pressure.⁵⁻⁹ These observations are in excellent agreement with the earlier DFT calculations, which employ several different types of electron-helium pseudopotentials to include the correlation effects into the model.²,³ Due to the quantum analog of the Bernoulli force, electron bubbles may also become trapped on quantized vortex lines and rings present in superfluid helium.¹,¹⁰,¹¹ In fact, the first experimental visualization of quantized vortex lines was based on this principle.¹²,¹³

Dynamic phenomena of solvated electrons in superfluid helium have been mostly studied experimentally where they were subject to either an external static electric field

(i.e., electron mobility experiments) or ultrasound excitation to explode the electron bubbles by negative pressure waves.¹,¹⁴,¹⁵ The latter experiments established that electrons can be optically detected in the liquid through such expanded bubbles.¹⁶ The electron mobility experiments, on the other hand, have demonstrated that such measurements can yield detailed information about the dissipative interaction with the surrounding superfluid.¹⁷⁻²⁰ In the low velocity regime (a few m/s), the measurements have shown that the electron mobility is reduced by collisions with the thermal excitations present in the liquid, i.e., phonons and rotons. Three different temperature regimes have been approximately identified at saturated vapor pressure: (1) phonon limited mobility when $T < 0.6$ K, (2) $0.6$ K $< T < 1.2$ K where the dissipation takes place primarily by discrete roton - electron bubble collisions and phonons, and (3) $1.2$ K $< T < 2.17$ K where the roton collision frequency is sufficiently high such that continuum models become applicable for describing the effective interaction.¹ The transition temperature between cases 2 and 3 will also be discussed further in this paper as it appears slightly different than given in Ref. 1. At higher drift velocities (some tens of m/s), a critical threshold value has been observed after which the electron mobility suddenly begins to decrease.¹,²¹⁻²⁴ This transition point has been attributed to the creation of rotons or vortex rings when the electron reaches the Landau critical velocity corresponding to the temperature, pressure, and geometry used in the experiment.

Time-of-flight measurements of the solvated electrons in an external electric field have revealed that, in addition to the normal electron bubble, a sequence of additional negative ion signals were also present. Such faster than electron “exotic ion” signals have been reported over the years by many groups²⁵⁻²⁹ and various explanations for their origin have been proposed in the literature:³⁰,³¹ (1) fractional

<footnote>a)E-mail: Jussi.Eloranta@csun.edu</footnote>
electrons where the electron wavefunction is split into different solvation bubbles, i.e., a superposition state, (2) electronically excited negatively charged He⁻ and He⁻,³²⁻³⁷ and (3) solvated impurities (e.g., O⁻, O⁻₂ H₂). However, none of the previously given explanations fully account for the experimental observations and hence the problem of identifying these exotic ions remains open.

To date, theoretical modeling of the time-dependent liquid helium response in the presence of an excess electron has received relative little attention in the literature.¹ For example, the early work of Rosenblit and Jortner used a hydrodynamic model to estimate the free electron localization time in liquid helium.³⁸ In the superfluid phase, time-dependent DFT calculations have been used to predict the variation of the energy level separations between 1s, 2s, 1p, and 2p states following an optical excitation.³ This calculation demonstrated the feasibility of performing picosecond-times scale optical pump-probe experiments on electron bubbles to monitor their evolution. Recently, a similar DFT model has been applied to evaluate the hydrodynamic added mass of negative and positive charge carriers in superfluid helium.³⁹ As the current DFT models do not include the thermal excitations that are responsible for viscous dissipation, Aitken et al. have developed a phenomenological model where the hydrodynamic radius is related to the free volume and thermodynamic state equations that describe its variation with phase, pressure, and temperature.⁴⁰,⁴¹ The state equations are calibrated to the known experimental data for the particular ion and medium under consideration. This model has been successful in describing electron mobility in supercritical and liquid helium,⁴⁰ helium gas at low temperature,⁴¹ and positive ion mobility in liquid helium⁴² and supercritical helium.⁴³ Although the model itself is less computationally demanding than DFT and not limited to a specific system, it obscures many of the microscopic details of ion solvation and the role of dissipative processes.

In this work, we extend the time-dependent Orsay-Trento DFT (OT-DFT) model to include the viscous response of superfluid helium present at finite temperatures. This is accomplished by transforming the classical viscous response terms from the Navier-Stokes equation to DFT by using the Madelung transformation. The overall viscous response in this model is not only determined by the bulk liquid shear viscosity but also by the variation of viscous response at the inhomogeneous bubble interface. The latter contribution is obtained by calibrating the model against the experimentally known electron mobilities in superfluid helium along the saturated vapor pressure line. Finally, mobilities for several different bubble sizes in the bulk liquid are predicted as a function of temperature, which can be applied, for example, to identify unknown ionic species based on their experimentally observed mobility.

## II. THEORY

### A. Density functional theory model

Previously, time-independent OT-DFT has been successfully applied to model electron bubbles in superfluid helium

as a function of both pressure and temperature.²,³,⁴⁴⁻⁴⁶ In this model, the coupled non-linear Schrödinger equations corresponding to superfluid helium OT-DFT and the solvated electron are expressed as:²,³,⁴⁶

$$ \left\{ \begin{aligned} i\hbar \frac{\partial \psi}{\partial t} &= -\frac{\hbar^2}{2m_{He}} \Delta \psi + V_{OT}[\rho] \psi + V_{ext}[\rho_e] \psi \\ i\hbar \frac{\partial \psi_e}{\partial t} &= -\frac{\hbar^2}{2m_e} \Delta \psi_e + V_{ext}[\rho] \psi_e, \end{aligned} \right. \qquad (1) $$

where $m_{He}$ is the helium atom mass, $m_e$ is the electron mass, $\psi$ is the liquid order parameter with the corresponding one-particle density given as $\rho = |\psi|^2$, $\psi_e$ is the electron wavefunction with $\rho_e = |\psi_e|^2$, $V_{OT}$ is the non-linear OT-DFT potential (excluding the kinetic correlation and backflow terms),⁴⁶ and the external potential functional is defined as:

$$ V_{ext}[\rho] = \int \rho(r')V_{ps}(|r-r'|)d^3r' \qquad (2) $$

where $V_{ps}$ represents the electron-helium pseudopotential⁴⁷ and the argument $\rho$ corresponds to either liquid helium or electron density. The previously developed finite-temperature OT-DFT includes only the static contributions of the liquid thermal response whereas the dynamic dissipative processes (e.g., viscosity) are neglected.⁴⁵ To include the dissipation due to viscous liquid into Eq. (1), the connection between the OT-DFT equation (first line of Eq. (1)) and the hydrodynamic Euler equation must first be briefly considered. By using the Madelung transformation (i.e., $\psi = \sqrt{\rho} \exp(iS/\hbar)$),⁴⁸ the OT-DFT equation can be written as (the corresponding continuity equation not shown),

$$ \rho_m \left( \frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \vec{\nabla}) \vec{v} \right) = -\rho \vec{\nabla} \xi \qquad (3) $$

where the helium mass density is defined as $\rho_m = m_{He} \rho$, the liquid velocity field $\vec{v} = \vec{\nabla} S/m_{He}$, and the effective pressure is given by:

$$ \xi = V_{OT} + V_{ext} - \frac{\Delta\sqrt{\rho}}{2m_{He}\sqrt{\rho}}. \qquad (4) $$

The last term on the right hand side of Eq. (4) is termed quantum pressure. In the Navier-Stokes equation, additional terms are included on the right hand side of Eq. (3) to account for viscous response. Following the Stokes hypothesis (i.e., monoatomic gas), we set the bulk viscosity to zero and write the Navier-Stokes equation as (corresponding continuity equation not shown),

$$ \rho_m \left( \frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \vec{\nabla}) \vec{v} \right) = -\rho \vec{\nabla} \xi + \vec{\nabla} \cdot \left[ \eta \left( \vec{\nabla} \vec{v} + (\vec{\nabla} \vec{v})^T - \frac{2}{3} (\vec{\nabla} \cdot \vec{v}) \vec{1} \right) \right], \qquad (5) $$

where $\eta = \eta(r)$ represents the shear viscosity and $\vec{1}$ is the unit tensor. The term contained between the brackets corresponds to the symmetrized linear stress tensor. Note that this equation is valid in the presence of both liquid compression ($\vec{\nabla} \cdot \vec{v} \neq 0$) and rotational flow ($\vec{\nabla} \times \vec{v} \neq 0$). To recast the viscous terms of Eq. (5) into the appropriate form for OT-DFT, the viscous
response functional, $\xi_{NS}$, must be solved from

$$ - \rho \nabla \xi_{NS} = \vec{\nabla} \cdot \left[ \eta \left( \vec{\nabla} \vec{v} + \left( \vec{\nabla} \vec{v} \right)^T - \frac{2}{3} \left( \vec{\nabla} \cdot \vec{v} \right) \vec{1} \right) \right]. \quad (6) $$

This is equivalent to the following Poisson equation for $\xi_{NS}$:

$$ \Delta\xi_{NS} = -\vec{\nabla} \cdot \left( \frac{1}{\rho} \vec{\nabla} \cdot \left[ \eta \left( \vec{\nabla} \vec{v} + \left( \vec{\nabla} \vec{v} \right)^T - \frac{2}{3} \left( \vec{\nabla} \cdot \vec{v} \right) \vec{1} \right) \right] \right). \quad (7) $$

We assume that both $\rho$ and $\vec{v}$ vary slowly around the electron bubble such that the viscosity becomes a functional of density only. An often employed form for this functional dependence at a given temperature is $\eta[\rho] \propto \rho^\alpha$ where $\alpha$ is a constant defining the variation of viscosity across the interface ("interfacial viscosity coefficient"). While no rigorous justification for employing this form for liquids exists, it may be considered to bear similarities with the polynomial-type expansions used for expressing viscosity of gases in terms of density.49,50 The following ansatz is employed to represent the interfacial viscosity:

$$ \eta_T [\rho, T] = \left( \frac{\rho}{\rho_0(T)} \right)^{\alpha(T)} \eta_0(T), \quad (8) $$

where $\eta_0$ is the shear viscosity of the bulk liquid, which is determined by the product of the normal fluid viscosity and the normal fluid fraction (i.e., $\eta_0 = (\rho_n/\rho_0)\eta_n$). Note that this power form can also approximately describe the density dependence of viscosity for classical liquids.51,52

Since there is no known analytic form for function $\alpha(T)$, it must be obtained by calibration against the experimentally available electron mobility data.14 The following empirical form was used to approximate the temperature dependence:

$$ \alpha(T) = a_0 + a_1 \exp(a_2 T), \quad (9) $$

where the empirical constants $a_0$, $a_1$, and $a_2$ were found to be 1.73, 2.32 × 10⁻¹⁰, and 11.15 K⁻¹, respectively, by minimizing the difference between the viscous OT-DFT calculations and experimental mobility data.14

To model a moving electron in superfluid helium, the first equation in Eq. (1) is solved in the coordinate frame moving with the electron by introducing the following kinetic energy constraint operator, $\hat{T}_c$, to the Hamiltonian:53

$$ \hat{T}_c \psi = -\vec{v}_0 \cdot \vec{p} \psi, \quad (10) $$

where $\vec{v}_0$ is the velocity vector corresponding to the electron movement (liquid flowing in opposite direction) and $\vec{p}$ is the momentum vector operator for the liquid. For simplicity, the electron motion is restricted along the x-axis and hence $\vec{v}_0 = (v_0, 0, 0)$. Note that when periodic boundary condition is imposed on the simulation box, $v_0$ must be an integer multiple of $\frac{2\pi}{N_x\Delta x}$, where $N_x$ is the number of discretization points along the x-axis and $\Delta x$ is the spatial step length. Therefore, the lowest applicable constrained velocity is dictated by the spatial extent of the grid. For the present calculations, $v_0 = 2.3$ m/s, which is well within the experimentally observed linear mobility regime and especially much lower than the Landau critical velocity for roton/vortex creation.1

The steady-state flow solution to Eq. (1) can be conveniently achieved by using the imaginary time method.46,54 Since the calculation is carried out in the moving electron

frame of reference, it converges to the steady-flow solution around the electron with no net acceleration. Thus the viscous drag force and the force due to the external electric field must be equal in strength but in opposite directions, $\vec{F}_{drag} = -\vec{F}_{field}$. The drag force can be evaluated from the steady-state solution of Eq. (1) by

$$ \vec{F}_{drag} = - \int \int \vec{\nabla}_r (\rho_e(r)V_{ps}(|r-r'|)) \rho(r')d^3r'd^3r. \quad (11) $$

Since the viscous response term in Eq. (7) introduces dissipation near the gas-liquid interface, asymmetric variations in both the liquid velocity and density appear between the front and the back of the bubble relative to the incident flow. According to Eq. (11), such an asymmetric density profile is responsible for producing the drag force. The incident field strength at the electron, $V$, is directly proportional to the external force, $\vec{F}_{field}$, acting on it: $V = |\vec{F}_{field}|/e$, where $e$ is the elementary charge. In the linear regime, the electron mobility, $\mu$, can now be calculated from (in units of m² V⁻¹ s⁻¹),

$$ \mu \equiv \frac{v_0}{V} = \frac{ev_0}{|\vec{F}_{field}|} = \frac{ev_0}{|\vec{F}_{drag}|}. \quad (12) $$

Furthermore, the following derived quantities are useful for analyzing the results:39,55

$$ m_{add} = \int \rho(r)v_x(r)d^3r/v_0, \quad (13) $$

$$ R = \frac{e}{4\pi\mu\eta}, \quad (14) $$

where $m_{add}$ is the hydrodynamic added mass of the electron bubble and $R$ is the Stokes radius of the bubble based on the Stokes law. Note that the viscosity used above corresponds to the viscosity of the whole fluid, which is determined by the product of the normal fluid fraction and the normal fluid viscosity.

To predict the dependence of mobility on the bubble radius, the following artificial external potential, $V_e$, was employed:

$$ V_e(r) = c_0 \exp(-c_1(r - c_2)), \quad (15) $$

where constants $c_0$ and $c_1$ where fixed to 3.8003 × 10⁵ K and 1.6245 Å⁻¹, respectively. The bubble radius corresponding to this exponentially repulsive potential can be varied by changing the value of $c_2$. For example, choosing $c_2 = 10.1$ Å, results in a density profile that resembles the electron bubble structure very closely. In such calculations, only the liquid portion of Eq. (1) was considered with $V_{ext}$ given by Eq. (15) and integration over the (non-existent) electronic degrees of freedom was omitted in Eq. (11).

The model defined by Eq. (1) consists of two coupled non-linear Schrödinger equations (i.e., superfluid helium and electron), which were solved numerically in imaginary time as described previously.2,3,46,54,56 The applied imaginary time step for the liquid degrees of freedom was 100 fs whereas a 0.1 fs time step was used for the electron due to its lower mass. The spatial grid applied in the calculations consisted of 1024 × 512 × 512 grid points with a step size of 0.8 bohr. Since the interaction between the electron and the liquid involves exchange of long wavelength phonons,
the simulation box must be sufficiently large to accommodate
such excitations. Note that the energy of these phonons is
very low and therefore the mobility converges rather slowly
in imaginary time, particularly near 2.1 K. The convergence
criterion in the present calculations was based on the observed
difference between the consequent mobilities calculated at
every 500 imaginary time iterations. When this difference fell
below 0.1% of the mobility, the calculation was considered
to be converged. Typically 10⁴ imaginary iterations were
required to reach this threshold. However, by using a solution
obtained previously at a nearby temperature or bubble size, the
convergence time could be reduced significantly. Furthermore,
the applied predictor-corrector scheme used for solving Eq. (1)
allowed the imaginary time step length to be increased up
to 100 fs without introducing any time step bias from the
applied operator splitting scheme into the solution.⁴⁶ The
bulk liquid density was maintained during the imaginary time
iterations by introducing a constant chemical potential to
Eq. (1), which corresponds to the applied OT-DFT functional.
Because the calculations were carried out in imaginary time, it
was possible to employ regular periodic boundary conditions
in the simulation box. On the contrary, real time calculations
would require the use of absorbing boundaries to eliminate
back reflections of waves that originate from the accelerating
object. To compute the liquid viscous response functional, the
Poisson equation given by Eq. (7) was solved numerically
in the Fourier space by using the centered finite difference
formula⁵⁷ where a constant of 10⁻⁹ a.u. was added to ρ to
avoid numerical instabilities arising from the division by ρ.
Superfluid helium shear viscosity values at the considered
temperatures along the saturated vapor pressure curve were
taken from the literature.¹⁴

Finally, we note that the OT-DFT part of Eq. (1) could
also be written separately for the superfluid and normal
fluid order parameters (“two-fluid DFT”). In this case the
OT-DFT potential should be evaluated by using the total
liquid density, which is the sum of the normal and super
fractions, ρ = ρ<sub>s</sub> + ρ<sub>n</sub>. Only the normal fluid equation should
include the viscous response determined by η<sub>n</sub>. However, this
model allows unrestricted spatial separation of the normal
and superfluid components, which is not compatible with the
original thermal OT-DFT functional. For the solvated electron
considered in this work, both Eq. (1) and this two-fluid OT-
DFT model produced very similar results for the calculated
electron mobilities where in the latter case only a small spatial
separation between the two liquid fractions was observed
(<2% of bulk density). Within this model, the normal fluid
accumulates in the front of the electron bubble whereas the
superfluid tends to fill the void behind the bubble.

### III. RESULTS AND DISCUSSION

Plotting the inverse of bulk liquid shear viscosity against
the experimentally observed mobility indicates nearly linear
behavior down to ca. 1.2 K as demonstrated in Fig. 1. Such
behavior is predicted by the Stokes law, Eq. (14), assuming
that the bubble radius is independent of temperature. However,
careful examination shows that this result holds only very
approximately and, thus, for the Stokes law to hold exactly,

<figure><img src="image_2.png" /><figcaption>FIG. 1. Experimental electron mobility ($\mu$) is plotted vs. the inverse bulk shear viscosity ($\eta_0^{-1}$) as black circles.¹⁴ The blue continuous line demonstrates the near linear relationship above 1.2 K and the red dashed line is provided as a guide to the eye to highlight the non-linear behavior below 1.2 K. A red dotted vertical line marks the approximate transition point, which separates the hydrodynamic continuum and discrete roton collision regimes (see text for discussion).</figcaption></figure>

$R$ must be temperature dependent, $R = R(T)$. Below 1.2 K temperature the behavior becomes notably non-linear (cf. “transition point” indicated in Fig. 1), which can be interpreted as the point where the hydrodynamic continuum model breaks down and discrete roton scattering events at the electron bubble dominate the dissipative response (“roton gas”).

An example of the superfluid helium density and the
velocity field obtained from the imaginary time solution of
Eq. (1) is shown in Fig. 2. Since the electron drift velocity
is low (≈2 m/s), the velocity field appears nearly symmetric
between the front and the back of the bubble. The presence of
the viscous term given by Eq. (7) is responsible for producing
the liquid flow pattern that moves around the bubble. Note
that the behavior of the hydrodynamic variables $\rho$ and $\vec{v}$
may become problematic at the bubble interface and inside

<figure><img src="image_3.png" /><figcaption>FIG. 2. Liquid density contours ($\rho$) and the associated velocity field ($v$) for an electron moving in superfluid helium at 2.1 K obtained by viscous OT-DFT model. The velocity field shown corresponds to ($v_x - v_0, v_y, v_z$) where the liquid background moves at 2.3 m/s.</figcaption></figure>
<table><caption>TABLE I. Viscous OT-DFT optimized values for the exponent $\alpha(T)$ in Eq. (8), the resulting electron mobility $\mu_{DFT}$ from Eq. (12), the experimental electron mobility $\mu_{exp}$,¹⁴ the calculated hydrodynamic added mass $m_{add}$ from Eq. (13) in units of $m_{He}$, the Stokes radius $R_{Stokes}$ based on $\mu_{DFT}$ (Eq. (14)), and the bubble radius at half bulk density value $R_{\rho0/2}$.</caption><thead><tr><th>T (K)</th><th>α</th><th>μ<sub>DFT</sub> (cm<sup>2</sup> V<sup>-1</sup> s<sup>-1</sup>)</th><th>μ<sub>exp</sub> (cm<sup>2</sup> V<sup>-1</sup> s<sup>-1</sup>)</th><th>m<sub>add</sub></th><th>R<sub>Stokes</sub> (Å)</th><th>R<sub>ρ0/2</sub> (Å)</th></tr></thead><tbody><tr><td>1.2</td><td>2.00</td><td>0.984</td><td>0.988</td><td>273</td><td>27.9</td><td>18.2</td></tr><tr><td>1.4</td><td>1.58</td><td>0.363</td><td>0.364</td><td>256</td><td>33.7</td><td>18.4</td></tr><tr><td>1.6</td><td>1.62</td><td>0.176</td><td>0.177</td><td>279</td><td>34.3</td><td>18.5</td></tr><tr><td>1.8</td><td>1.95</td><td>0.108</td><td>0.109</td><td>281</td><td>29.4</td><td>18.6</td></tr><tr><td>2.0</td><td>2.80</td><td>0.0655</td><td>0.0686</td><td>286</td><td>23.9</td><td>18.7</td></tr><tr><td>2.1</td><td>5.20</td><td>0.0503</td><td>0.0505</td><td>304</td><td>19.0</td><td>18.8</td></tr></tbody></table>

the bubble itself where $\rho \approx 0$ but $\vec{v} \neq 0$. On the other hand, the order parameter representation handles such situations naturally because both its real and imaginary parts approach zero when $\rho \to 0$ and the velocity does not need to be considered. This may be important in situations where high liquid velocities are present in low-density liquid regions (e.g., rapid bubble expansion dynamics). Ideally, for this reason, functionals that do not depend on liquid velocity or involve division by liquid density should be employed.

Results from the viscous OT-DFT model calibration calculations are shown in Table I. At each temperature listed, minimization of the difference between the calculated (Eq. (12)) and the experimental mobilities with respect to parameter $\alpha(T)$ was carried out. The optimized values of $\alpha$ are plotted in Fig. 3 as a function of temperature. It can be seen that $\alpha$ varies smoothly from 5.20 at 2.1 K down to 1.58 at 1.4 K but begins to increase rapidly below 1.2 K. Note that this takes place near the same temperature where the mean free path of rotons becomes larger than the size of the electron bubble and, consequently, the Knudsen number (Kn) exceeds 1 (i.e., gaseous in terms of roton density).¹ For parameter $\alpha$ this effect is very pronounced as its value essentially approaches

infinity below 1.0 K. To extend the present model to this temperature regime, a Cunningham-Millikan-type correction factor should be included to provide an effective value for viscosity.⁵⁸ An approximate value for $\alpha(T)$ between 1.4–2.1 K can be obtained from Eq. (9), but for accurate calculations, the exact values provided in Table I should be employed instead. The viscosity values predicted by Eq. (8) are plotted in Fig. 4 as a function of the liquid density. At 1.4 K, the limited cold helium gas viscosity data⁵⁹ appears comparable to the values obtained from the viscous OT-DFT optimization calculations. However, the gas phase values should be somewhat lower because the cold gas can only carry thermal phonons and, at this temperature, thermal rotons also contribute to liquid viscosity. Unfortunately, such data is not available at other temperatures to allow for more exact comparison. For the electron bubble, such spatial density variations occur at the gas-liquid interface, which, in turn, contribute to the interfacial viscous drag through Eq. (7). Finally, a comparison between the experimental electron mobility data obtained along the saturated vapor pressure line and the results from the viscous OT-DFT model employing Eq. (8) is shown in Fig. 5. Note that this agreement would naturally be exact if the optimized

<figure><img src="image_4.png" /><figcaption>FIG. 3. The exponent $\alpha$, which defines the dependence of viscosity on liquid density, is plotted as a function of temperature. An exponential fit (see Eq. (9)) is shown for comparison by a red line. Note that the continuum approximation begins to break down below 1.2 K where the normal fluid Knudsen number (Kn) for rotons becomes greater than one (dotted black line). The dashed black line highlights the rapid increase below 1.2 K.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Dependence of liquid helium viscosity ($\eta$) on density ($\rho$) at selected temperatures as predicted by Eq. (8) with the values of $\alpha(T)$ specified in Table I. The gas phase data were taken from Ref. 59 with the experimental viscosity values scaled by the normal fluid fraction to account for the two-fluid nature of superfluid helium. Note the logarithmic scale for $\eta$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Comparison between the experimental and viscous OT-DFT electron mobilities (μ) as a function of temperature (T) along the saturated vapor pressure line.14,60 The interfacial viscosity coefficient α(T) was taken from Eq. (8). Note the logarithmic scale for mobility.</figcaption></figure>

values given in Table I were used. Despite the approximate form used for expressing α(T) with Eq. (8), the model still reproduces the experimental electron mobility values with a good accuracy.

The small deviations of the experimental data from the Stokes law in the linear regime (see Fig. 1) imply that the Stokes radius must include some additional temperature dependent contribution. This is demonstrated in Fig. 6 where the calculated Stokes radius reaches a maximum value of ca. 35 Å between 1.4 and 1.6 K. In general, below 2.1 K, the Stokes radius appears much larger than the ρ₀/2 radius (≈18.5 Å) obtained from the OT-DFT density profiles. This means that the Stokes radius should be interpreted as an effective radius that incorporates additional factors such as temperature dependent changes in the interfacial viscosity. Coincidentally, when the normal liquid fraction approaches one (T ≈ 2.1 K), the Stokes and ρ₀/2 radii become equal.

<figure><img src="image_7.png" /><figcaption>FIG. 6. Temperature dependence of electron bubble Stokes radius and radius defined by ρ₀/2 density obtained from OT-DFT calculations. The lines connecting the discrete points are drawn as guides to the eye.</figcaption></figure>

Ion mobilities predicted by the developed viscous OT-DFT model for selected bubble sizes at specified temperatures are shown in Fig. 7. The calculated values are observed to follow the general power-law:

$$ \mu(R) = \mu_e \left( \frac{R_e}{R} \right)^{\gamma}, \qquad (16) $$

where μ_e is the electron mobility at the given temperature with bubble radius R_e (reference system) and μ and R are the corresponding quantities for the (unknown) species. The exponent γ = 2 when the mobility is limited by roton scattering (see, for example, Refs. 31 and 61) and γ = 1 for the continuum Stokes limit. Least squares fits of Eq. (16) to the calculated OT-DFT mobilities with respect to the exponent γ are indicated as continuous lines in Fig. 7. In general, the obtained exponent values vary between the two limits but appear closer to the Stokes limit within the temperature range considered. While the temperature variation of γ appears systematic, it is not clear what microscopic phenomenon is responsible for this behavior. Nevertheless, the form given by Eq. (16) can be used to estimate ion bubble ρ₀/2 radii (as opposed to an effective hydrodynamic radius) based on their experimentally observed mobilities. Note, however, that such an estimate assumes that (1) ion is solvated in a bubble rather than in a snowball,62 (2) the interfacial viscosity (see Eq. (8)) is independent of the bubble radius, and (3) the bubble radius is large enough that the continuum model is applicable. To avoid complications from the last assumption, experiments measuring ion mobilities in superfluid helium should be carried out close to the lambda temperature. On the contrary, however, ion mobilities can be related to the

<figure><img src="image_8.png" /><figcaption>FIG. 7. Predicted ion mobilities (μ) for spherical bubbles as a function of ρ₀/2 radius (Rρ₀/2) at temperatures indicated in the legend. The bubbles were defined by the repulsive potential of Eq. (15) where c₂ (a.u.) was set to 6.0, 10.0, 12.0, 16.0, or 19.0 to yield bubbles with the corresponding R values of 11.1 Å, 13.4 Å, 14.6 Å, 16.9 Å, and 18.6 Å at 2.1 K. The continuous lines and exponents γ shown in the legend correspond to Eq. (16) with the applied (μ_e, Rρ₀/2,ε) values of (0.051, 18.6), (0.066, 18.6), (0.109, 18.5), (0.187, 18.4), (0.379, 18.4), and (0.970, 18.4) in cm² V⁻¹ s⁻¹ and Å units at temperatures 2.1, 2.0, 1.8, 1.6, 1.4, and 1.2 K, respectively.</figcaption></figure>
bubble radii more accurately at low temperatures (cf. Fig. 7) where the variation of mobility vs. bubble radius is more pronounced.

An example application of the data given in Fig. 7 is to estimate the bubble radii for unknown ions, such as the exotic ions, in superfluid helium. The ion mobilities of exotic ions extrapolated to 1.2 K vary between 1.4 cm<fnref n="2" /> V⁻<fnref n="1" /> s⁻<fnref n="1" /> (Ion #8) and 3.5 cm<fnref n="2" /> V⁻<fnref n="1" /> s⁻<fnref n="1" /> (Ion #1).⁶³ Furthermore, the observed “Fast Ion” species has a very high mobility, 6.0 cm<fnref n="2" /> V⁻<fnref n="1" /> s⁻<fnref n="1" />. Assuming a spherical geometry for the corresponding bubble, the data shown in Fig. 7 predicts that the radius for Ion #8 is ca. 13.6 Å. Although the present calculations do not extend to radii below 11 Å, Eq. (16) can be used to extrapolate the OT-DFT results such that Ion #4 radius is estimated to be around 10 Å. The negative helium ions, such as He⁺⁻ and He⁺⁻, are expected to have solvation cavities in this size regime but their limited observed lifetime in the gas phase as well as the existence of only two such ionic species are solid arguments against assigning them to these exotic ion signals (see discussion, for example, in Ref. 31; see also Ref. 32). However, they could contribute to the experimentally observed continuous background due to ionization during their flight through the cell.⁶³ Our preliminary *ab initio* electronic structure calculations at the coupled clusters level of theory (CCSD(T))⁶⁴ with the basis set specified in Ref. 65 indicate that He⁺⁻ species is also stable in the symmetric linear geometry (bond lengths 1.22 Å). However, at present, it is not clear if this type of negative ion series, He⁺⁻, extends any further than *n* = 3. Note that typical impurities in the ion mobility experiments include oxygen from air, nitrogen from air, water, and possible residues of hydrocarbons used for pumping the sample cell. The high energy process used to produce the electrons (e.g., discharge) can also dissociate and ionize these species, producing many fragments with various electron affinities. Examples of possible ions that can form include O⁻, O₂⁻, OH⁻, C⁻, C₂⁻, H⁻, H₂, and possibly even N⁻ and N₂⁻. The neutral parents of the latter two species have, however, very low electron affinities. If the previously mentioned extrapolation procedure is applied to fast ion,⁶³ a radius of 4 Å is obtained. Such a small bubble could be formed by a high electron affinity atom or diatomic molecule (e.g., O⁻ or O₂⁻).

Finally, we note that the present viscous OT-DFT method is not just limited to describing electron mobility in superfluid helium, but it can also model ions with a complex solvation structure. Our preliminary calculations show that, for example, the predicted ion mobilities for Ca⁺ and Sr⁺ at 1.2 K are 1.06 cm<fnref n="2" /> V⁻<fnref n="1" /> s⁻<fnref n="1" /> and 0.97 cm<fnref n="2" /> V⁻<fnref n="1" /> s⁻<fnref n="1" />, respectively, which can be compared with the corresponding experimental values of 1.1 cm<fnref n="2" /> V⁻<fnref n="1" /> s⁻<fnref n="1" /> and 1.2 cm<fnref n="2" /> V⁻<fnref n="1" /> s⁻<fnref n="1" />.⁶⁶ Both ions have a partial snowball solvation structure (i.e., hollow core with a radius of approximately 4 Å)⁶⁷ where the first rigidly bound solvation layer follows the ion. This, in turn, reduces the ion mobility and increases its apparent size in the liquid. Note that, based on the calculated hollow core sizes and the position of the first solvation shell,⁶⁷ Ca⁺ should have a slightly higher mobility than Sr⁺, which is opposite to the experimental findings.⁶⁶ The only input to these calculations is the ion–helium pair potential that was obtained

from earlier electronic structure calculations.⁶⁷ Application of the developed method to describe mobilities of positive ion snowballs and the behavior of the model outside the saturated vapor pressure line will be considered in a future publication.

## ACKNOWLEDGMENTS

Financial support from the National Science Foundation Grant Nos. CHE-0949057 and DMR-1205734 to J.E. and the Leverhulme Trust (Res. Grant F00212AH), the Royal Society (Int. Exchange Grant IE130173), and COST action CM 1405 (Molim) to K.V.H. are gratefully acknowledged.

<footnote n="1">¹A. F. Borghesani, *Ions and Electrons in Liquid Helium* (Oxford Science Publications, Oxford, 2007).</footnote>

<footnote n="2">²V. Grau, M. Barranco, R. Mayol, and M. Pi, *Phys. Rev. B* **73**, 064502 (2006).</footnote>

<footnote n="3">³J. Eloranta and V. A. Apkarian, *J. Chem. Phys.* **117**, 10139 (2002).</footnote>

<footnote n="4">⁴J. R. Casey, A. Kahros, and B. J. Schwartz, *J. Phys. Chem. B* **117**, 14173 (2013).</footnote>

<footnote n="5">⁵C. C. Grimes and G. Adams, *Phys. Rev. B* **41**, 6366 (1990).</footnote>

<footnote n="6">⁶C. C. Grimes and G. Adams, *Phys. Rev. B* **45**, 2305 (1992).</footnote>

<footnote n="7">⁷A. Y. Parshin and S. V. Pereverzev, *JETP Lett.* **52**, 282 (1990).</footnote>

<footnote n="8">⁸A. Y. Parshin and S. V. Pereverzev, *JETP* **74**, 68 (1992).</footnote>

<footnote n="9">⁹S. V. Pereverzev and A. Y. Parshin, *Physica B* **197**, 347 (1994).</footnote>

<footnote n="10">¹⁰D. Mateo, J. Eloranta, and G. A. Williams, *J. Chem. Phys.* **142**, 064510 (2015).</footnote>

<footnote n="11">¹¹R. J. Donnelly, *Quantized Vortices in Helium II*, Studies in Low Temperature Physics Vol. 3 (Cambridge University Press, Cambridge, 1991).</footnote>

<footnote n="12">¹²G. A. Williams and R. E. Packard, *Phys. Rev. Lett.* **33**, 280 (1974).</footnote>

<footnote n="13">¹³G. A. Williams and R. E. Packard, *J. Low Temp. Phys.* **39**, 353 (1980).</footnote>

<footnote n="14">¹⁴R. J. Donnelly and C. F. Barenghi, *J. Phys. Chem. Ref. Data* **27**, 1217 (1998).</footnote>

<footnote n="15">¹⁵A. Ghosh and H. J. Maris, *Phys. Rev. Lett.* **95**, 265301 (2005).</footnote>

<footnote n="16">¹⁶W. Guo, D. Jin, G. M. Seidel, and H. J. Maris, *Phys. Rev. B* **79**, 054515 (2009).</footnote>

<footnote n="17">¹⁷L. Meyer and F. Reif, *Phys. Rev. B* **110**, 279 (1958).</footnote>

<footnote n="18">¹⁸F. Reif and L. Meyer, *Phys. Rev. B* **119**, 1164 (1960).</footnote>

<footnote n="19">¹⁹L. Meyer and F. Reif, *Phys. Rev. B* **123**, 727 (1961).</footnote>

<footnote n="20">²⁰L. Meyer and F. Reif, *Phys. Rev. Lett.* **5**, 1 (1960).</footnote>

<footnote n="21">²¹G. W. Rayfield and F. Reif, *Phys. Rev. Lett.* **11**, 305 (1963).</footnote>

<footnote n="22">²²G. Careri, S. Cunsolo, P. Mazzoldi, and M. Santini, *Phys. Rev. Lett.* **15**, 392 (1965).</footnote>

<footnote n="23">²³G. W. Rayfield, *Phys. Rev. Lett.* **20**, 1467 (1968).</footnote>

<footnote n="24">²⁴D. M. Strayer and R. J. Donnelly, *Phys. Rev. Lett.* **26**, 1420 (1971).</footnote>

<footnote n="25">²⁵M. P. Seah and D. P. Woodruff, *Phys. Lett. A* **30**, 250 (1969).</footnote>

<footnote n="26">²⁶G. G. Ihas and J. T. M. Sanders, *Phys. Rev. Lett.* **27**, 383 (1971).</footnote>

<footnote n="27">²⁷V. L. Eden and P. V. E. McClintock, *Phys. Lett. A* **102**, 197 (1984).</footnote>

<footnote n="28">²⁸C. D. H. Williams, P. C. Hendry, and P. V. E. McClintock, *Jpn. J. Appl. Phys.* **26**, 105 (1987).</footnote>

<footnote n="29">²⁹W. Wei, Z.-L. Xie, G. M. Seidel, and H. J. Maris, *J. Low Temp. Phys.* **171**, 178 (2013).</footnote>

<footnote n="30">³⁰A. G. Khrapak and W. F. Schmidt, *Low Temp. Phys.* **37**, 494 (2011).</footnote>

<footnote n="31">³¹H. J. Maris, *J. Low Temp. Phys.* **120**, 173 (2000).</footnote>

<footnote n="32">³²T. Andersen, *Phys. Scr.* **T59**, 230 (1995).</footnote>

<footnote n="33">³³H. T. Schmidt, P. Reinhed, A. Orbán, S. Rosén, R. D. Thomas, H. A. B. Johansson, J. Werner, D. Misra, M. Björkhage, L. Brännholm *et al.*, *J. Phys. Conf. Ser.* **388**, 012006 (2012).</footnote>

<footnote n="34">³⁴H. H. Michels, *Chem. Phys. Lett.* **126**, 537 (1986).</footnote>

<footnote n="35">³⁵H. H. Michels, *Phys. Rev. Lett.* **52**, 1413 (1984).</footnote>

<footnote n="36">³⁶Y. K. Bae, M. J. Coggiola, and J. R. Peterson, *Phys. Rev. Lett.* **52**, 747 (1984).</footnote>

<footnote n="37">³⁷A. Mauracher, M. Daxner, J. Postler, S. E. Huber, S. Denifl, P. Scheier, and J. P. Toennies, *J. Phys. Chem. Lett.* **5**, 2444 (2014).</footnote>

<footnote n="38">³⁸M. Rosenblit and J. Jortner, *J. Phys. Chem. A* **101**, 751 (1997).</footnote>

<footnote n="39">³⁹D. Mateo and J. Eloranta, *J. Phys. Chem. A* **118**, 6407 (2014).</footnote>

<footnote n="40">⁴⁰F. Aitken, Z.-L. Li, N. Bonifaci, A. Denat, and K. von Haeften, *Phys. Chem. Chem. Phys.* **13**, 719 (2011).</footnote>

<footnote n="41">⁴¹F. Aitken, N. Bonifaci, A. Denat, and K. von Haeften, *J. Low Temp. Phys.* **162**, 702 (2011).</footnote>
⁴²F. Aitken, N. Bonifaci, L. G. Mendoza-Luna, and K. von Haeften, Phys. Chem. Chem. Phys. **17**, 18535 (2015).

⁴³H. G. Tarchouna, N. Bonifaci, F. Aitken, L. G. Mendoza-Luna, and K. von Haeften, J. Phys. Chem. Lett. **6**, 3036 (2015).

⁴⁴F. Dalfovo, A. Lastri, L. Pricaupenko, S. Stringari, and J. Treiner, Phys. Rev. B **52**, 1193 (1995).

⁴⁵F. Ancilotto, F. Faccin, and F. Toigo, Phys. Rev. B **62**, 17035 (2000).

⁴⁶L. Lehtovaara, T. Kiljunen, and J. Eloranta, J. Comput. Phys. **194**, 78 (2004).

⁴⁷J. Jortner, N. R. Kestner, S. A. Rice, and M. H. Cohen, J. Chem. Phys. **43**, 2614 (1965).

⁴⁸E. Madelung, Z. Phys. **40**, 322 (1927).

⁴⁹J. C. Rainwater, J. Chem. Phys. **81**, 495 (1984).

⁵⁰D. G. Friend and J. C. Rainwater, Chem. Phys. Lett. **107**, 590 (1984).

⁵¹K. R. Harris, P. J. Newitt, and L. A. Woolf, J. Chem. Eng. Data **49**, 138 (2004).

⁵²K. Rah and B. C. Eu, J. Chem. Phys. **114**, 10436 (2001).

⁵³F. Ancilotto, M. Barranco, and M. Pi, Phys. Rev. B **82**, 014517 (2010).

⁵⁴L. Lehtovaara, J. Toivanen, and J. Eloranta, J. Comput. Phys. **221**, 148 (2007).

⁵⁵G. Batchelor, *Introduction to Fluid Dynamics* (Cambridge University Press, Cambridge, 1970).

⁵⁶L. Lehtovaara, D. Mateo, and J. Eloranta, Implementation of the Orsay-Trento density functional for describing superfluid helium, 2015, available on http://sourceforge.net/projects/libdft/.

⁵⁷W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes in C: The Art of Scientific Computing*, 2nd ed. (Cambridge University Press, Cambridge, 2002).

⁵⁸E. Cunningham, Proc. R. Soc. A **83**, 357 (1910).

⁵⁹A. van Itterbeek, F. W. Schapink, G. J. van den Berg, and H. J. M. van Beek, Physica **19**, 1158 (1953).

⁶⁰K. W. Schwarz, Phys. Rev. A **6**, 837 (1972).

⁶¹G. I. Hlas, Ph.D. thesis, University of Michigan, 1971.

⁶²K. R. Atkins, Phys. Rev. **116**, 1339 (1959).

⁶³W. Wei, Z. Xie, L. N. Cooper, G. M. Seidel, and H. J. Maris, J. Low Temp. Phys. **178**, 78 (2015).

⁶⁴H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz et al., *MOLPRO*, version 2012.1, a package of *ab initio* programs, 2012, see http://www.molpro.net.

⁶⁵J. Eloranta and V. A. Apkarian, J. Chem. Phys. **115**, 752 (2001).

⁶⁶W. I. Glaberson and W. W. Johnson, J. Low Temp. Phys. **20**, 313 (1973).

⁶⁷S. L. Fiedler, D. Mateo, T. Aleksanyan, and J. Eloranta, Phys. Rev. B **86**, 144522 (2012).