Absorption spectrum of Ca atoms attached to $^4$He nanodroplets

Within density functional theory, we have obtained the structure of $^4$He droplets doped with neutral calcium atoms. These results have been used, in conjunction with newly determined {\it ab-initio} $^1\Sigma$ and $^1\Pi$ Ca-He pair potentials, to address the $4s4p$ $^1$P$_1 \leftarrow 4s^2$ $^1$S$_0$ transition of the attached Ca atom, finding a fairly good agreement with absorption experimental data. We have studied the drop structure as a function of the position of the Ca atom with respect of the center of mass of the helium moiety. The interplay between the density oscillations arising from the helium intrinsic structure and the density oscillations produced by the impurity in its neighborhood plays a role in the determination of the equilibrium state, and hence in the solvation properties of alkaline earth atoms. In a case of study, the thermal motion of the impurity within the drop surface region has been analyzed in a semi-quantitative way. We have found that, although the atomic shift shows a sizeable dependence on the impurity location, the thermal effect is statistically small, contributing by about a 10% to the line broadening. The structure of vortices attached to the calcium atom has been also addressed, and its effect on the calcium absorption spectrum discussed. At variance with previous theoretical predictions, we conclude that spectroscopic experiments on Ca atoms attached to $^4$He drops will be likely unable to detect the presence of quantized vortices in helium nanodrops.


I. INTRODUCTION
Optical investigations of atomic impurities in superfluid helium nanodroplets have drawn considerable attention in recent years. 1,2 In particular, the shifts of the electronic transition lines with respect to the gas-phase transition lines (atomic shifts) are a very useful observable to determine the location of the foreign atom attached to the drop. Alkaline earth atoms appear to play a unique role in this context. While e.g., all alkali atoms reside in surface "dimple" states, and more attractive impurities like all noble gas atoms reside in the bulk of drops made of either isotope, 3 the absorption spectra of heavy alkaline earth atoms attached to 4 He drops clearly support an outside location of Ca, Sr, and Ba, 4,5 whereas for the lighter Mg atom the experimental evidence is that it resides in the bulk of the 4 He droplets. 6,7 We have recently presented Density Functional Theory (DFT) results for the structure and energetics of large 3 He and 4 He doped nanodroplets, showing that alkaline earth atoms from Mg to Ba go to the bulk of 3 He drops, whereas Ca, Sr and Ba reside in a dimple at the surface of 4 He drops, and Mg is in their interior. 8 This is in agreement with the analysis of available experimental data, although the case of Mg has been questioned very recently. 9 Moreover, according to the magnitude of the observed shifts, the dimple for alkaline earth atoms was thought to be more pronounced than for alkali atoms, indicating that the former reside deeper inside the drop than the later. This has been also confirmed by the calculations. In addition, the 5s5p ← 5s 2 experimental transition of Sr atoms attached to helium nanodroplets of either isotope has shown that strontium is solvated inside 3 He nanodroplets, also in agreement with the calculations. 8 Calcium atoms are barely stable on the surface of the drop, and the difference between the energy of the surface dimple state and that of the solvated state in the bulk of the drop is rather small and depends very sensitively on the Ca-He interatomic potential. 10 The aim of this work is to obtain the atomic shifts for Ca attached to large 4 He drops, and to compare them with the experimental data. To this end, we have improved our DFT approach, 8 treating the atomic impurity as a quantal particle instead of as an external field.
Laser Induced Fluorescence (LIF) experiments for Ca atoms in liquid 3 He and 4 He have been reported 11 and analyzed within a vibrating bubble model, which involves the formation of a bubble around the impurity, using Ca-He pair potentials based on pseudopotential SCF/CI calculations. 12 This work is organized as follows. In Sec. II we discuss the Ca-He interaction potentials we have used. In Sec. III we briefly present our density functional approach, as well as some illustrative results for the structure of Ca@ 4 He N drops. The method we have employed to obtain the atomic shifts is discussed in Sec. IV. In Sec. V we present the results obtained for calcium, discuss how thermal motion may affect the line shapes, and investigate how the presence of a quantized vortex line may change the Ca absorption spectrum. Finally, a summary is presented in Sec. VI. Figure 1 shows the X 1 Σ Ca-He adiabatic potential obtained by different authors. 13,14,15,16 This potential determines the dimple structure described in the next Section. It can be seen that apart from the unpublished potential by Meyer, the others are quite similar. As in our previous work, 8 we shall use the one obtained in Ref. 15. This will allow us to ascertain the effect of treating Ca as a quantal particle. Since the X 1 Σ potential seems to be fairly well determined, we have turned our attention to the excited adiabatic potentials.

II. CALCIUM-HELIUM INTERACTION POTENTIALS
In a previous work 17 the excited state potentials were calculated in a valence ab-initio scheme. The core electrons of calcium and helium were replaced by scalar-relativistic energyconsistent pseudopotentials, and the energy curves were calculated on the complete-activespace multiconfiguration selfconsisted field (CASSCF) 18 The calculations were performed within the MOLPRO suite of ab-initio programs. 24 The molecular orbitals used for the excited states calculations were optimized in the state averaged CASSCF method for all singlet states correlating to (4s 2 ) 1 S, (4s3d) 1 D, and (4s4p) 1 P atomic asymptotes. The active space was formed by distributing the two valence electrons of the Ca atom into 4s3d4p valence orbitals. The 1s2s2p3s3p orbitals of calcium and 1s orbital of helium were kept doubly occupied in all configuration state functions, but they were optimized in the CASSCF calculations. The resulting wave functions were used as references in the following ICMRCI calculations. At the ICMRCI level, the 1s2s2p3s3p orbitals of calcium and 1s orbital of helium were kept as doubly occupied in all reference configuration state functions, but these orbitals were correlated through single and double excitations. The 4s, 4p, and 3d calcium orbitals had not restricted occupation patterns. We show in Fig. 2 the excited adiabatic potentials; to have better insight into the potential minima of the 1 Σ and 1 Π potentials, they have been plotted correlating to the (4p) 1 P Ca term.

III. DFT DESCRIPTION OF HELIUM NANODROPLETS
In recent years, static and time-dependent density functional methods 25,26,27 have become increasingly popular to study inhomogeneous liquid helium systems because they provide an excellent compromise between accuracy and computational effort, allowing to address problems inaccessible to more fundamental approaches, see e.g. Ref. 3 for a recent review.
Obviously, DFT cannot take into account the atomic, discrete nature of these systems, but can address inhomogeneous helium systems at the nanoscale 28 and take into account the anisotropic deformations induced by some dopants in helium drops. Both properties are essential to properly describe these systems.
Our starting point is the Orsay-Trento density functional, 25 together with the Ca-He adiabatic potential X 1 Σ of Ref. 15, here denoted as V Ca−He . This allows us to write the energy of the Ca-drop system as a functional of the Ca wave function Φ(r) and the 4 He "order parameter" Ψ(r): The order parameter is defined as Ψ(r) = ρ(r) exp[ıS(r)], where ρ(r) is the particle density and v(r) =h∇S(r)/m 4 is the velocity field of the superfluid. In Eq. (1), E(ρ) is the 4 He "potential energy density". 25 In the absence of vortex lines, we set S to zero and E becomes a functional of ρ and Φ. Otherwise, we have used the complex order parameter Ψ(r) to describe the superfluid.
We have solved the Euler-Lagrange equations which result from the variations with respect to Ψ * and Φ * of the energy E[Ψ, Φ] under the constrain of a given number of helium atoms in the drop, and a normalized Ca wave function, namely: where µ is the helium chemical potential and ε is the lowest eigenvalue of the Schrödinger equation obeyed by the Ca atom. The effective potentials U He and U Ca are defined as The coupled Eqs.
On the figure are shown also the results obtained treating calcium as an external field. 8 It can be seen that for large drops the energy of the calcium atom is about 10 K less negative due to its zero point motion. We want to stress again how barely stable is the calcium atom on the surface of 4 He N drops. For instance, we have found that the total energy of the The dimple depth ξ, defined as the difference between the position of the dividing surface at ρ = ρ 0 /2, where ρ 0 = 0.0218Å −3 is the bulk liquid density, with and without impurity, respectively, is shown in Fig. 4. Due to the zero point motion that pushes the impurity towards lower helium densities, for large drops the dimple depth is about 0.8Å smaller when the zero point motion is included than when it is not. 8 This change in the depth is large enough to produce observable effects in the calculated absorption spectrum, as discussed below.
It can be seen that the dimple depth curve ξ(N) has some structure. This is not a numerical artifact, but a genuine effect due to the interplay between the Ca atom and the drop, whose density, even for pure drops, shows conspicuous oscillations all over the drop volume, extending up to the surface region irrespective of whether the drop is described within DFT or Diffusion Monte Carlo methods. 3,25,32 The interplay of these oscillations with those arising from the presence of the impurity little affects the total energy of the system and hence, the Ca energy, but yields some visible structure in the density distributions that shows up in related quantities, like the dimple depth. To illustrate it, we show in Fig where Z is the average distance in the z direction between the impurity and the geometrical center of the helium moiety and λ C is an arbitrary constant, large enough to guarantee that upon minimization, Z equals the desired Z 0 value. We have also applied this method to Mg doped helium drops, and will present the details of the calculation elsewhere. 33 We show in the bottom panel of Fig This behavior has not been disclosed before, and appears in the course of "pushing" the impurity inside the droplet from its equilibrium position (a similar structure shows up for Ca@ 4 He 500 , but at more negative ∆Z 0 values). While it affects rather little the equilibrium location of the Ca atom because of its clear surface location, and hence the atomic shift -see however the inset in Fig. 8, it plays a substantial role in the solvation of magnesium atoms in small helium drops. 33 We recall that this problem has recently drawn the attention of experimentalists and theoreticians as well. 6,7,8,9,35,36 How the position of the Ca atom affects the absorption spectrum throughout the change in the dimple structure will be discussed in Sec. V.

IV. EXCITATION SPECTRUM OF AN ATOMIC IMPURITY IN A 4 HE DROP
Lax method 37 offers a realistic way to study the absorption spectrum I(ω) of a foreign atom embedded in helium drops. It makes use of the Franck-Condon principle within a semi-classical approach, and it has been employed to study the absorption spectrum of several atomic dopants attached to fairly small 4 He drops. 35,38,39,40 The case of alkali atoms attached to large drops described within DFT has been also considered, see Refs. 41,42 and references therein. Lax method is usually applied in conjunction with the diatomics-in-molecules theory, 43 which means that the atom-drop complex is treated as a diatomic molecule, where the helium moiety plays a role of the other atom.
In the original formalism, to obtain I(ω) one has to carry out an average on the possible initial states of the system that may be thermally populated.

A. Line shapes
The line shapes for electronic transitions from the ground state (gs) to the excited state (ex) in a condensed phase system can be written as where D ge is the matrix element of the electric dipole operator, H gs and H ex are the Hamiltonians which describe the ground and excited states of the system respectively, and |Ψ gs represents the ground state. I(ω) can be evaluated using the Born-Oppenheimer approximation, which makes a separation of the electronic and nuclear wave functions |Ψ i = |e i |ψ i , and the Franck-Condon principle, whereby the heavy nuclei do not change their positions or momenta during the electronic transition. If the excited electron belongs to the impurity, the helium cluster remains frozen, so that the relevant coordinate is the relative position r between the cluster and the impurity. That principle amounts to assuming that D ge is independent of the nuclear coordinates. Taking into account that e − it h Hgs |Ψ gs = e −itωgs |Ψ gs and projecting on eigenstates of the orbital angular momentum of the excited electron |m , wherehω gs X and ψ gs X (r) are the energy and the wave function of the ro-vibrational ground state of the frozen helium-impurity system, and H ex m (r) is the ro-vibrational excited Hamiltonian with potential energy V ex m (r) determined by the electronic energy eigenvalue, as obtained in the next subsection for a p ← s transition. Eq. (9) will be referred to as the Fourier Formula, and it is the Fourier transform of the time-correlation function. It is nothing but a sum at the resonant energies weighted with the well-known Franck-Condon factors: where ω m ν and |ψ m ν are the ro-vibrational eigenvalues and eigenstates of the Hamiltonian If the relevant excited states for the transition have large quantum numbers, they can be treated as approximately classical 37,38,39 using the averaged energyhω m ν ≈ V ex m (r) which is independent of ν. In this case we obtain the expression where Ω m (ω) is the surface defined by the equation ω + ω gs X − V ex m (r)/h = 0. We will refer to this equation as the Semi-Classical Formula.
If the atom is in bulk liquid helium, or at the center of the drop, the problem has spherical symmetry and the above equation reduces to where r m (ω) is the root of the equation ω + ω gs X − V ex m (r)/h = 0. In the non-spherical case, we have evaluated I(ω) from the first expression in Eq. (10) using the discretization where Θ is the step function and ∆ω is a frequency step small enough so that the above discretization represents the delta function. We also take advantage that only points near the impurity contribute to the integral, by writing ∆r ijk +r 0 = (i∆x+x 0 , j∆y+y 0 , k∆z+z 0 ), with ∆x, ∆y, ∆z being the spatial mesh steps used in the discretization, and being r 0 an arbitrary point in the neighborhood of the impurity. Finally, we recall that I(ω) needs to be evaluated only in a narrow frequency range starting from an arbitrary ω 0 which can be, e.g., the free atom frequency. This range defines the maximum n value in the above equation.

B. Excited ro-vibrational potential for a p ← s transition
We next determine the potential energy surfaces V ex m (r) needed to carry out the calculation of the atomic shifts.

Pairwise sum aproximation
The pair-interaction between an atom in a s-state and an atom in a p-state can be expressed, in the cartesian eigenbasis (|x , |y , |z ) as where V Π (r) and V Σ (r) are the adiabatic potentials neglecting the spin-orbit interaction and r is the distance between atoms. For a system of N helium atoms and an excited impurity in a p-state, the total potential is approximated by the pairwise sum where r n is the distance between the n th helium atom and the impurity, and R n is the rotation matrix that transform the unity vectorẑ into ther n vector. It can be shown that, in cartesian coordinates, where x 1 = x, x 2 = y, x 3 = z, and r 2 n = x 2 n + y 2 n + z 2 n . Thus, the matrix elements of the total potential are Using the continuous density approach inherent to DFT [ n → d 3 r ′ ρ(r ′ )], this expression can be written as The eigenvalues of this symmetric matrix are the sought-after V ex m (r) which define the potential energy surfaces (PES) as a function of the distance between the centers of mass of the droplet and of the impurity, and are given by the three real roots λ i (r) of the equation It can be shown that for spherical geometry, Eq. (17) is diagonal with matrix elements (in spherical coordinates)

Spin-Orbit coupling
For atomic impurities in which the spin-orbit (SO) interaction is prominent and comparable to the splitting of the P -states due to the interaction with the droplet, it has to be taken into account in the calculation of the PESs. This is usually done considering that the SO splitting of the dopant is that of the isolated atom irrespective of the impurity-drop distance. 39,46,47 Given the atomic structure of alkaline earth atoms, the SO interaction can be safely neglected in the PES calculation. However, we discuss it here for the sake of completeness and future reference.
When the spin-orbit interaction is taken into account, the total potential can be written as V T = U + V SO , where V SO has the form, in the spin-cartesian orbit basis (|x, 1/2 , |x, −1/2 , |y, 1/2 , |y, −1/2 , |z, 1/2 , |z, −1/2 ): where A ℓs is 2/3 of the experimental SO splitting of an isolated atom. Kramers' theorem states that there is a two-fold degenerate manyfold of systems with a total half-integer spin value that cannot be broken by electrostatic interactions, 48 so that the two-fold degenerate eigenvalues that define the PES's are the roots of the equation with A, B and C defined in Eq. (19). In the case of spherical geometry, the eigenvalues adopt a simple expression:    13) shows. This explains why the main peak for N = 100 is the narrower one, and is the reason why the main peak is fairly apart from the Σ shoulder. As N increases, so it does the splitting, while the three components of the peak become broader. Eventually, if N is large enough, it is not possible to distinguish the components of the absorption line. It is also obvious that line broadening due to effects not considered here may wash out the blueshifted shoulder found for small N values.
The inset in Fig. 8 shows the calculated shifts relative to the gas-phase transition, compared with the experimental values. 4 One can appreciate a small oscillation for the largest drops; it is a genuine effect produced by the dimple structure. Further insight can be gained from the study of the absorption spectrum as a function of Z 0 . To this end, we display in the top panel of Fig. 9 the absorption spectra for the N = 500 droplet corresponding to seeing how the shift decreases as the distance between the centers of mass increases, and at the same time the absorption peak becomes more asymmetric as the Ca environment does (see also Subsection C).

B. Thermal broadening
As we have indicated, for very attractive impurities, the foreign atom is fully solvated, and its delocalization in the bulk of the drop due to thermal motion hardly introduces a significant change in the line shape. When the dopant is in a dimple state, it has to be checked whether thermal motion may have observable effects on the absorption spectrum, as Fig. 6 seems to indicate for calcium.
To ascertain this effect, we have carried out a thermal average of the spectrum using an approximate expression for the probability density. Referring the energy E i of a given Z 0 i configuration to the equilibrium value, ∆E i = E i − E gs , and neglecting the kinetic energy of the impurity and the displaced fluid, we write the probability density for the position Z 0 i of the Ca atom as where k B is the Boltzmann constant, T = 0.4 K, and the sum -actually integral-runs on the selected Z 0 i configurations. 55 The Z 2 0 i factor takes care of the relative volume available to each configuration. With this definition the probability of finding the Ca atom between Z 0 i and Z 0 i + ∆Z 0 i is w i ∆Z 0 i . We show in the top panel of Fig. 6 the probability densities w i corresponding to the configurations displayed in the bottom panel. We see that there is a non-negligible probability of finding the Ca atom in a broad region of the drop surface, and consequently we have addressed, as a case of study, the statistical properties of the Ca@ 4 He 500 system at this temperature.
We have found that the mean position, calculated as Z 0 = Σ i Z 0 i w i ∆Z 0 i , and the standard deviation, calculated as σ(Z 0 ) = Z 2 0 − Z 0 2 are 17.38Å and 0.76Å, respectively. This dispersion in the position generates a dispersion in the value of the shift. To quantify this effect we have evaluated the mean value of the shift, calculated as ∆ω = Σ i ∆ω i w i ∆Z 0 i , and its standard deviation, calculated as σ(∆ω) = ∆ω 2 − ∆ω 2 , using the values shown in the inset of the top panel of Fig. 9, which correspond to the seven lower energy configurations in Fig. 6. We have obtained ∆ω ± σ(∆ω) = 63.8 ± 11.5 cm −1 .
The thermally averaged Ca absorption spectrum is shown in the bottom panel of Fig. 9 (solid line), as well as that corresponding to the equilibrium configuration (dashed line). To carry out the average, we have used the I i (ω) in the top panel of Fig. 9, and averaged them as I(ω) = Σ i I i (ω)w i ∆Z 0 i . This procedure, consistent with the Franck-Condon principle, assumes that absorption proceeds instantaneously on any of the frozen drop-Ca configurations characterized by a Z 0 i value.
We are led to conclude that the thermal motion effect is rather small. This has prompted us to re-analyze the structure of a large N = 1000 drop hosting a calcium atom attached to a vortex line along the symmetry axis. Within DFT, a robust method to generate vortex configurations in liquid helium is described in Ref. 59. We adapt it here to the case of helium drops. For a n = 1 quantum circulation vortex line, we start the imaginary time evolution to solve Eq. (2) from the initial state if x and y are nonvanishing, and zero otherwise, where ρ(r) is the vortex-free Ca@ 4 He 1000 helium density. After the minimization procedure is converged, we have checked that the obtained final configuration is indeed a n = 1 vortex state. and without a vortex line along its symmetry axis. It can be seen that the vortex line draws the impurity towards the bulk of the droplet, but it still resides in a deeper surface dimple. Figure 11 shows the absorption spectrum for the two configurations displayed in Fig. 10.
The effect of the presence of the vortex on the absorption spectrum of calcium is twofold. On the one hand, the maximum of the absorption peak is shifted towards the bulk value because of the deeper dimple. On the other hand, the FWHM increases by about a factor of two.
The reason is the spreading of the Ca wave function within the stretched U Ca (r) well, that allows the atom to "probe" a wider region in the excited PES, thus increasing the width.
Notice also the larger splitting of the peaks that form the line due to the more anisotropic helium environment. This also contributes to increasing the width of the absorption peak.
Unfortunately, the experimental absorption line is so broad and asymmetric that the extra shift caused by the vortex is not enough to displace the line to a region where it could be distinguishable on top of the vortex-free absorption line.

VI. SUMMARY
Within density functional theory, we have carried out a detailed study of the absorption spectrum of calcium atoms attached to 4 He N drops in the vicinity of the 4s4p 1 P 1 ← 4s 2 1 S 0 transition, finding a semi-quantitative agreement with experiment. To this end, we have improved our previous implementation of the DF method by incorporating the zero point motion of the impurity, and have carried out ab-initio calculations to obtain the excited 1 Σ and 1 Π Ca-He potentials needed to obtain the potential energy surfaces.
We have studied the drop structure, finding that the "interference" between the density oscillations of the helium moiety arising from its intrinsic structure and those arising from the presence of the impurity plays a role in the determination of the position of the impurity.
This may be relevant for the solvation of alkaline earth atoms, especially for magnesium. 35 In a case of study, we have systematically addressed the dependence of the relative shift