Advances in Quantum Monte Carlo - American Chemical Society

Chapter 12

Rotations and Exchange in Doped Helium Clusters: Insight from Imaginary-Time Correlation Functions

Downloaded by UNIV OF GUELPH LIBRARY on July 3, 2012 | http://pubs.acs.org Publication Date: December 31, 2006 | doi: 10.1021/bk-2007-0953.ch012

Nicholas Blinov and Pierre-Nicholas Roy Department of Chemistry, University of Alberta, Edmonton, Alberta T6G 2G2, Canada

We analyze the behavior of doped helium clusters in terms of imaginary time correlation functions. Illustrative examples are presented in order to highlight the importance of the inclusion of quantum exchange effects in the path integral Monte Carlo calculation of imaginary time orientational correlation functions. We relate the oscillatory behavior of the rotational constant to the variation of the minimum value of the orientational correlation functions for Bose-Einstein statistics.

Introduction Doped helium clusters provide a feature-rich environment for the study of quantum statistical effects over a broad range of size scales. In the nano-scale regime, rovibrational spectra of doped helium nanodroplets show sharp spectral lines and a renormalized moment of inertia of the dopant species (1,2,3). In the case of small to medium size doped helium clusters, a non-monotonic size evolution of the effective rotational constant has been observed experimentally © 2007 American Chemical Society

In Advances in Quantum Monte Carlo; Anderson, J., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2006.

165

166 for N 0 (4,5) and C 0 (6) impurities, for instance, and predicted via theoretical simulations for O C S (7,8) and N 0 (5,9,10), and C 0 (6,11,12) dopant species. The observed features have been attributed to quantum exchange effects associated with the bosonic character of the helium environment. Theoretically, the importance of exchange effects on the rotational dynamics of a doped cluster can be assessed with the help of path-integral Monte Carlo (PIMC) simulations (5,9,12,14). One of the most interesting results which so far can only be obtained from P I M C simulations is a qualitatively correct prediction of size evolution of the rotational constant for small clusters (in the range of cluster sizes corresponding to the completion of the first solvation shell) based on a two-fluid hydrodynamic model (5,11). The approach allows one to relate the superfluid response of helium to the rotational dynamics of the doped cluster but relies on some assumptions (such as a rigid coupling between the normal fraction of helium and the dopant). It is therefore interesting to provide evidence for the importance of exchange effects without invoking these assumptions. The ultimate test would be based on the energy spectrum of the doped cluster, but such an approach is impractical even for small (more than three helium atoms) clusters due to exponential scaling of the complexity of the problem. We discuss here some aspects of the extraction of rotational constants from imaginary-time orientational correlation functions. 2

2


2

2

Formal aspects and practical considerations The most practical (numerically exact) approach to treat the rotational dynamics of relatively large clusters is based on the analysis of imaginary time dipole-dipole correlation functions. These functions can be evaluated using ground state quantum Monte Carlo (QMC) (7,8) or P I M C approaches (9,13,14). In the rigid rotor approximation, the dipole-dipole correlation function is proportional to an orientational correlation function defined as

0) where the unit vector η represents the orientation of the molecule, β is the reciprocal temperature, Ζ is the partition function, and Η is the Hamiltonian of the system. In the path-integral picture the difference between the ground state and the finite temperature (canonical) average in the above equation appears as a difference in topology of the Feynman paths. In the ground state, open paths are


167 involved, while at finite temperature, only closed paths contribute to the physical properties (other than the Green's functions). There are also important differences in the eigenstate (spectrum) representation of the correlation functions. In the ground state, the correlation functions are represented by the T

sum of decaying exponential functions C ( )

=

g

2-d

energy eigenstate of the Hamiltonian: Ηψ = ε^,

ε

s

\ *

a n

I

2

labels the eigenstates, and A = n J . Here n is a matrix element of the operator i i taken between the ground state, g, and the fth eigenstate. For long enough projection times (τ-»αο) the contribution of the highly excited states to the correlation function can be projected out. The excitation energies corresponding to the transitions involving the ground and few (usually, one or two) excited sates can then be extracted from the imaginary time propagation via the use of a multi-exponential fit or from maximum entropy approaches (7,5). For finite temperatures, the spectral representation of the correlation i


e r e

i is a composite index that

ί

ι

w

**' ' ^

gi

function of Eq. (1) is: C^(r) = — ^

Α^

τ < ΰ ή

, where ω - ε - ε is the ί}

]

ι

transition frequency and τ varies in the interval [0, β]. The important difference with the ground state is that each decaying exponential function is matched here by a growing exponential function associated with the corresponding reverse transition. This guarantees the periodicity of the correlation function in imaginary time. This is made especially clear using the following representation of the correlation function β( +ε;)/2

Γ

C/,(r) = X ( 2 - ^ y 0

εί

Ay COSh(fA

(2)

where τ e [-β/2,β/2]. In contrast to ground state simulations, the correlation function cannot be propagated long enough in imaginary time in order to project out excited states for an arbitrary temperature. Still, in the low temperature regime (β-»οο), transitions involving only the ground and the first excited states will dominate. In this case, a time interval, τ e [0,/?/2], should exist where the evolution of the correlation function will be determined by the low lying energies (just as in the ground state case). A multi-exponential fitting procedure can then be applied on such an interval in order to extract the excitation energies, again, as in ground state applications. It should be advantageous to take into account the explicit spectral form of the correlation function in Eq. (2) in the fitting procedure at finite temperature (13). If the temperature is low enough, and the amplitudes of transitions are treated as fitting parameters along with the


168 transition frequencies, such a procedure should give exact excitation energies (and intensities).


The challenge of fitting The low temperature regime required to obtain accuracies on par with ground state calculations demands a large number of discretizations in the pathintegral representation for the density matrix. Such an implementation is currently impractical for the bosonic case. The actual temperatures used in previous P I M C calculations (0.15-0.4K) are too high to single-out the contribution of a single excited state to the multi-exponential sum in the expression for the correlation function. For relatively high temperatures, i.e. βΒ ~ 1, (where Β is the effective rotational constant of the complex) the above fitting may become extremely unreliable because many exponential functions contribute simultaneously to the correlation function for any τ. In this case a model fit can be used as a qualitative way to analyze the rotational dynamics of doped helium clusters. The difference here, when comparing to the above multiexponential fit, is that one uses a specific model of the transition amplitudes and excitation energies. Rigid top and ro-vibrational Hamiltonian models are possible. If the actual rotational dynamics can be accurately described by a rigid top model, this procedure provides quantitative results for the rotational constant. The accuracy of the fit also provides some insight into the accuracy of the model, and so provides useful information (72). We analyze this type of approach in the following section.

Illustrative examples To circumvent the technical difficulties associated with the contribution of many transitions to the Monte Carlo sampling in the low temperature regime, we first choose to illustrate the above considerations using basis-set calculations for the He-OCS dimer. The approach and results have been partially discussed in Ref. (13). The method allows us to obtain (exact) bound energy levels and related eigenfunctions of the He-OCS Hamiltonian. We can then evaluate canonical averages for any physical observable including imaginary time correlation functions. We present the ro-vibrational energies obtained using parity-adapted basis functions by Lanczos diagonalization (75,76) in Table I. The energy levels are labeled by the total angular momentum J, the parity, p, and the vibrational quantum number, n. These energies are degenerate in M


169 1

4

Table I. Energy levels (cm ) for He-OCS dimer


η 0 1 2 3

JoPo -18.6 -10.3 -9.2 -6.5

JiPi

JiPo -18.0 -7.8 -6.5 -4.0

J2P1 -17.3 -16.6 -7.1 -6.7

J2P0 -17.7 -17.5 -16.6 -9.6

-18.3 -18.1 -10.1 -9.0

(the projection of the total angular momentum on the z-axis in the laboratory fixed frame). The imaginary time orientational correlation function can now be calculated as follows,

(«(r)«z(0))4 Σ·"**"*"·** Σ*-*" ){^ Ρ\*ζ\»'' 4 > :Ρ

Μ

Μ

2

z

(3)

where the matrix elements of the direction cosines are, (nJMp\h \n'J'Mp')= z

Σ&"α™ '(^Μρ\ή [ΓkMp'), Ρ

(4)

ζ

ojjc

where c ^

p

is an element of the eigenvectors for state \nJMp) in the parity

adapted ro-vibrational basis (See Refs. (13) and (77) for details). The non-zero matrix elements of the direction cosine in the parity adapted rotational basis are, {jkMp\n \JkMp>) =

(1 - δ .),

z

( ^ k

+

l ^ )

=

V ^ )

'

2

/

2

(5)

ρρ

^

+

1

)

2

- ^

(

G/+1)^(27+1)^7+3)

l - V ) , p p

(6)

2

(jkMp\ ,\j-XkMp') n

= V-/ -Wg^ , _ JiAJ -\

(

δA

(?)

2

Note that the above matrix elements are diagonal in M, and that the contributions from the other two direction cosines to the correlation function are the same. The


170 correlation function from Eq. (3) is plotted in Figure 1 for some representative temperatures. We can observe that the lower the temperature, the lower the amplitude of the correlation function at τ=β/2. To illustrate the contribution of individual energy levels to the canonical averages in Eq. (3), we show in Figure 2 the Boltzmann factors (including degeneracy), ( 2 J + l ) e ~ ^ , for some of the lowest energies.

0.4

T = 0.04K Τ = 0.1 Κ α Τ = 0.38Κ 1 w Τ= 1 κ


ο

Figure 1. Orientational imaginary time correlation functions for the He-OCS dimer. Results for different temperatures are presented.

It is clear from the above picture (Figure 2) that at 7M).04 Κ (and even at 7M). 1 K ) , the ground state contribution dominates. A multi-exponential fit would therefore yield results with accuracies comparable to that of a ground state analysis in this low temperature regime. On the other hand, at 7M).38 K , the temperature of actual nanodroplets experiments, the populations of the excited states are relatively high. Under these conditions, a multi-exponential fit could be unreliable, but one can use a model fit to obtain the information regarding the rotational dynamics. Because many excited states (with temperature dependent populations) contribute to the correlation function, the resulting effective rotational constant obtained from the model fit will depend on temperature. For larger clusters, the rotational constant will be renormalized and the Boltzmann weight distribution will be different.



171

Figure 2. Boltzmann weights for the lowest states of the He-OCS dimmer at Τ = 0.04 Κ O.IK 0.38 Κ and 1.0 Κ.

Path integral Monte Carlo calculations We illustrate the above considerations in Figure 3, where we plot the effective rotational constants obtained from fits to linear rotor and symmetric top models (75). The temperature dependence of the effective rotational constants for the He-OCS complex indicates a deviation from the rigid rotor models. This fact should not be surprising since the presence of the helium atom hinders the free rotation of the O C S monomer and this floppy complex should not be expected to behave as a simple free linear rotor. We present orientational correlation functions for ΗβΝ-Ν 0 clusters of varying sizes in Figure 4. The calculations are related to those presented in Ref. (5) and include Bose-Einstein statistics (exchange). The figure shows that the size evolution of the minimum value of the correlation functions at τ=β/2 exhibits an oscillatory behavior. This pattern is very similar to the size evolution of the Β rotational constant. To explore this further, we show in Figure 5 the size evolution of the reciprocal value of the correlation function at τ=β/2 for BoseEinstein and Boltzmann statistics. We observe that Bose-Einstein statistics are required in order to capture the experimentally observed two turnarounds in the size evolution. For comparison, we also show in the inset of Figure 5 the actual experimental and PIMC Β constants as reported in Réf. (J). The l/C(p/2) values and the Β rotational constants clearly have a similar behavior. We explore this connection in Figure 6 and observe that the result of the fit is directly correlated 2


172 0.3

1 ι

- - - -O

He-OCS dimer

S 0.2

o-Ô^symmtop - D - /r Q -o- ocs


B

111

w

"0

φ

0.5

1

1.5 T(K)

B

e

,lnear

H-OCS

2

e

*P

e r i m e n t

2.5

I

3

Figure 3. Effective rotational constants obtained from symmetric top (open circles) and linear molecule (open squares) modelfitsof the PIMC orientational correlation function. The experimental value (18) of the Β constant (filled diamond) is also presented along with the free OCS Β value (open triangles). (Reproduced with permission from Path Integral Monte Carlo Approach for Weakly Bound van der Waals Complexes with Rotations: Algorithms and Benchmark Calculations by Nicholas Bilnov. Copyright American Institute of Physics.)

Figure 4. Orientational correlation functions for various cluster sizes.



173

•

Boltzmann statistics

0

5

10

15

20

Ν

Figure 5. Size evolution of l/C(fi/2) for two statistics. Inset: experimental (filled triangles) and PIMC (open symbols) Β versus Ν (Réf. (5)).

•

•

ο0

•ο

ο Φ

0.05

φ * 0

0.1

0

Bose-Einstein statistics Boltzmann statistics Bose-Einstein statistics Boltzmann statistics

·

00

0.15

0.2

0.25

0.3

Β, Κ Figure

6. Correlation between l/C(fi/2) and the PIMC (filled symbols) and experimental (open symbols) Β constants.


174 to the \/C(fi/2) value. We also present the correlation between l/C(p/2) and the experimental Β values. The analysis shows that the Ι/Οφ/2) values obtained for Bose-Einstein statistics correlate much better with experiment than results where exchange has been neglected (i.e., Boltzmann statistics). The bosonic results still deviate from experiment since the height of the correlation function cannot simply be connected to the Β constant because of thermal excitations.


Concluding Remarks We have shown how the analysis of orientational imaginary time correlation functions obtained from PIMC simulations can provide insight into the size evolution of the Β rotational constants in doped helium clusters. We noted that for the case of an N 0 dopant, exchange effects (Bose-Einstein staitistics) must be included in order to reproduce the experimentally observed oscillatory behavior of the Β constants. This oscillatory behaviour has been connected to the onset of superfluidity where decoupling between a rotor and its environment is observed. The oscillation of the value of the correlation function at β/2 directly follows the trend in Β constants and therefore provides a measure of decoupling that is free from a fit. We note that a connection between C(p/2) and Β has recently been used by Miura (19) for a particular cluster size. 2

Acknowledgments We thank the Natural Sciences and Engineering Research Council of Canada and the Canada Foundation for Innovation for financial support.

References 1. 2. 3. 4. 5. 6.

Grebenev, S; Toennies, J.P.; Vilesov, A . F . Science 1998, 279, 2083. Toennies, J.P.; Vilesov, A.F., Annu. Rev. Phys. Chem. 1998, 49, 1. Reinhard, I.; Callegari, C.; Conjusteau, C.; Lehmann, K . K . ; Scoles, G . Phys. Rev. Lett. 1999, 82, 5036. X u , Y . ; Jäger, W.; Tang, J.; McKellar, A . R. W. Phys. Rev. Lett. 2003, 91, 163401. X u , Y . ; Blinov, N . ; Jäger, W.; Roy, P.-N.; J. Chem. Phys. 2006, 124, 081101. Tang, J.; McKellar, A . R. W.; Mezzacapo, F.; Moroni, S. Phys. Rev. Lett. 2004, 92, 145503.


175 7. 8.


9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19.

Moroni, S.; Sarsa, Α.; Fantoni, S.; Schmidt, S. Κ. Ε.; Baroni, S. Phys. Rev. Lett. 2003, 90, 143401. Paesani, F.; Viel, Α.; Gianturco, F.A.; Whaley, K . B . Phys. Rev. Lett. 2003, 90, 073401. Moroni, S.; Blinov, N . ; Roy, P.-N. J. Chem. Phys. 2004, 121, 3577. Paesani F.; Whaley, Κ. B . J. Chem. Phys. 2004, 121, 5293. Paesani, F.; Kwon, Y . ; Whaley, Κ. B . Phys. Rev. Lett. 2005, 94, 153401. Blinov, N . ; Roy, P.-N. Low Temp. Phys. 2005,140, 255. Blinov, N . ; Song, X . ; Roy, P.-N. J. Chem. Phys. 2004, 120, 5916. Zillich, R.; Paesani, F.; Kwon, Y . ; Whaley, K . B . J. Chem. Phys. 2005 123, 114301. Cullum, J.K.; Willoughby, R . A . Lanczos Algorithms for Large Symmetric Eigenvalue Computations; Birkhäuser: Boston, M A , 1985. Carrington, T. Jr., Encyclopedia of Computational Chemistry, edited by P. von Ragué Schleyer; Wiley: New York, N Y , 1998; V o l . 5. Song, X . ; X u , Y . ; Roy, P.-N.; Jäger, W. J. Chem. Phys. 2004, 121, 12308. Higgins, K . ; Klemperer, W. J. Chem. Phys. 1999, 110, 1383. Miura, S.; J. Phys.: Condens. Matter 2005, 17, S3259.


Advances in Quantum Monte Carlo - American Chemical Society

Recommend Documents