text:
# Comment on “Gravitating magnetic monopole in the global monopole spacetime”

## Abstract

We point out a problem with the stability of composite (global–magnetic) monopoles recently proposed by J. Spinelly, U. de Freitas and E.R. Bezerra de Mello \[Phys. Rev. D66, 024018 (2002)\].

Spinelly et. al. SBF  have written an interesting paper on composite monopoles, i.e., a global monopole and a magnetic monopole bound together by their gravitational interaction. This paper was followed by two works considering more general cases BH ; BBH . Such composite configurations are of interest for several reasons

One is in connection with the (magnetic) monopole problem. Global monopoles have long range interactions, and their annihilation rate is very efficient. Composite monopoles could share this property, if they remained bound as the global monopoles move to annihilate each other. One could even imagine a scenario in which the global and magnetic monopoles are created independently, with the magnetic monopoles coming from an earlier phase transition. A second transition leads to global monopoles which then capture and drag the magnetic monopoles with them as they go towards a neighbouring defect. Therefore magnetic monopole annihilation could be enhanced with respect to the case of purely magnetic monopoles. If there is no correlation between the signs of the global and magnetic charges, the effect is quite mild, and gives rise to a remnant population of doubly charged, purely magnetic monopoles. On the other hand, if the global and magnetic charges are correlated the effect could be dramatic, strong enough to solve the monopole problem.

Another reason is the spacetime structure around such defects. Pure global monopoles cause a deficit solid angle in the surrounding spacetime111The area of a sphere of radius $`r`$ surrounding the monopole is $`4\pi (1\mathrm{\Delta })r^2`$ $`\mathrm{\Delta }`$, plus a small effect due to the core BV . The mass that can be ascribed to the core turns out to be negative HL , leading to repulsive interactions. The presence of a magnetic monopole component can in principle change the sign of the core mass to make it positive SBF . Finally, the critical phenomena associated with gravitational collapse of the monopole to form a black hole show an interesting intermediate behaviour between that of the global and magnetic cases BH .

Our comment consists basically of two points. First, some of the results of SBF  have a very simple physical interpretation that was somewhat obscured by their choice of units and that we want to clarify. Second, and more important, we want to emphasise that the composite monopole that is bound only gravitationally is almost certainly unstable to breaking up into its global and magnetic monopole constituents; therefore, if it is to be cosmologically relevant some stronger type of interaction must be introduced between the two sets of fields. Non-minimal coupling has been considered in BBH , but only in spherically symmetric situations, and this problem was not addressed. More work will be needed in order to establish the existence of stable composite monopoles.

We begin by clarifying the interpretation behind the apparent dependence of the global monopole field profile with the parameter $`\beta `$, shown in figures 3a) and 6c) of ref. SBF , which is just an artifact of their choice of units.

Considering first the case where there is only a gravitating global monopole, we see that there are two mass scales in the problem. One is set by the v.e.v of the scalar field, $`\eta `$, the other by the massive scalar excitations, which have mass $`\sqrt{\lambda }\eta `$ ($`\lambda `$ is the quartic coupling). We can use the second one to define a natural length scale, and the first one will define an overall scale for the energy. In other words, in the absence of a magnetic monopole, the rescaling

$$\chi ^a\eta \chi ^a,x_\mu \widehat{x}_\mu =\frac{x_\mu }{\eta \sqrt{\lambda }},$$
(1)
in Lagrangian (7),(10) in reference SBF  makes the parameter $`\lambda `$ disappear from the equations of motion, while the parameter $`\eta `$ appears only in the adimensional combination $`G\eta ^2`$, which measures the gravitational strength of the global monopole (and, in particular, the deficit solid angle $`\mathrm{\Delta }=8\pi G\eta ^2`$ in the spherically symmetric configuration). The spherical monopole profile $`f(\widehat{r})`$ is known to be fairly independent of $`\mathrm{\Delta }`$ HL , and we expect this feature to persist in the presence of the magnetic monopole; this independence is observed in ref. SBF .

In general, the magnetic monopole would introduce three new mass scales: the v.e.v. of its scalar field, the mass of the scalar excitations and the mass of the vector boson, $`e\eta `$ ($`e`$ is the gauge coupling constant). In ref. SBF  the first two were chosen to coincide with those of the global monopole, so the only new adimensional parameter is the ratio of the scalar mass to the vector mass. The parameter $`\beta =\lambda /e^2`$ measures this ratio (squared). There are no more free parameters in the problem.

We do not expect the presence of the magnetic monopole to alter the global monopole profile $`f(\widehat{r})`$ substantially, since there is no direct coupling between them. In particular $`f`$ should be quite insensitive to the value of $`\beta `$, if $`\mathrm{\Delta }`$ is kept fixed. The only reason why figs. 3a) and 6c) in ref. SBF  appear to give such a dependence is that in SBF  the radial coordinate was rescaled using the vector boson mass instead of the scalar mass (a desirable choice in other respects). The difference is just a factor of $`\sqrt{\beta }`$ in the abscissa and once this is taken into account the curves agree perfectly, as shown in figure 1.

We now turn to the composite monopole’s stability. To understand the nature of the problem, consider the simpler case of the interaction between a magnetic monopole and a conical singularity (such as would be produced, for example, by an infinitely thin, straight, idealised cosmic string). In the plane perpendicular to the string, the metric is flat and has no effect on test masses other than lensing. But this is not true of gravitating masses and charges, whose long–range electric or magnetic fields can detect the presence of the singularity. The net effect turns out to be that masses are attracted to the conical singularity and charges (electric and magnetic) are repelled Smith . In the case of a magnetic charge, the convergence of magnetic field lines behind the singularity is responsible for the repulsion.

Now consider the three-dimensional version, a magnetically charged point particle in a spherically symmetric spacetime with a deficit solid angle concentrated in a point-like singularity. Without loss of generality the charge can be assumed to move in the equatorial plane ($`\theta =\pi /2`$).

This is a reasonable approximation to the interaction between a global monopole and a magnetic monopole that are well separated. But then the arguments of Smith  apply, and we expect the monopoles to repel. Moreover, in this case, even the gravitational interaction is repulsive, as the “mass” of the core of the global monopole is negative. An explicit calculation of the interaction between a pointlike electrically charged particle and a global monopole (including the effect of the core) confirms that it is repulsive BF ; we expect this result to hold when considering a magnetic charge and a global monopole.

We conclude that the cosmological capture of a magnetic monopole by a global monopole seems quite improbable (we have neglected the effect of the finite–size cores in this argument, but since the monopoles only interact gravitationally we do not expect our conclusions to change qualitatively). While this is disappointing, it is not necessarily a problem for the composite monopole since in that case the cores are superposed, the magnetic field lines are evenly distributed with spherical symmetry and thus there is no net repulsion between the two cores.

The problem arises when one considers the stability of the composite monopole to angular perturbations. This is an important consideration because, in flat space, global monopoles have a zero mode which allows the redistribution of gradient energy density in an axisymmetric way, without any cost in energy G . Following the notation in SBF , where the spherical global monopole is given by:

$$\chi ^a=f(r)\frac{x^a}{r},$$
(2)
purely angular deformations with axial symmetry can be described using the following ansatz G

$`\chi ^1`$ $`=`$ $`f(r)\mathrm{sin}\overline{\theta }(r,\theta )\mathrm{cos}\phi ;`$
$`\chi ^2`$ $`=`$ $`f(r)\mathrm{sin}\overline{\theta }(r,\theta )\mathrm{sin}\phi ;`$
$`\chi ^3`$ $`=`$ $`f(r)\mathrm{cos}\overline{\theta }(r,\theta ).`$ (3)
Note that $`\overline{\theta }=\theta `$ corresponds to the unperturbed monopole.

As shown by Goldhaber G , in flat space, configurations of the form

$`\mathrm{tan}\left({\displaystyle \frac{\overline{\theta }}{2}}\right)=\mathrm{tan}\left({\displaystyle \frac{\theta }{2}}\right)e^\xi ,\xi =\mathrm{const},`$ (4)
where $`\xi `$ is an arbitrary constant, are degenerate in energy with the spherical monopole. For any given $`r`$, $`\xi 1`$ corresponds to concentrating the gradient energy in an arbitrarily small region around the north-pole with no cost in energy.

The extra tension created at the north pole drags the core of the global monopole upwards, restoring spherical symmetry BR ; P ; AU ; WT . But consider its effect on the magnetic monopole: a larger gradient energy density in the global monopole causes a redistribution of the deficit angle “density”, with the larger contribution now coming from the north pole. But then magnetic field lines are pushed together there, causing the magnetic core to move downwards, i.e. in the opposite direction to the global monopole core. Once the cores are separated, they are expected to repel, so an attractive interaction between the cores seems necessary for the composite object to survive.

To summarise, the stability of composite monopoles can by no means be taken for granted and is expected to depend quite sensitively on the details of the interaction between the two constituents. Any attempt to consider composite monopoles in an astrophysical or cosmological context must address this problem.

###### Acknowledgements.

We are grateful to E.R. Bezerra de Mello for pointing out ref. BF . J.U. thanks the University of Leiden, where part of this work was done, for hospitality. This work is partially supported by the ESF COSLAB programme, and by AEN99-0315, FPA 2002-02037 and 9/UPV00172.310-14497/2002 grants.


format: text

text:
# Multiplication Formulas for the Elliptic Gamma Function

## 1. Introduction

Special functions defined by infinite products often have duplication formulas. Here are some examples.

$`\text{sin}(2\pi z)`$ $`=`$ $`2\text{sin}(\pi z)\text{sin}(\pi (z+{\displaystyle \frac{1}{2}})),`$
(1) $`\mathrm{\Gamma }(2z)\sqrt{\pi }`$ $`=`$ $`2^{2z1}\mathrm{\Gamma }(z)\mathrm{\Gamma }(z+{\displaystyle \frac{1}{2}}),`$
(2) $`\mathrm{\Gamma }_q(2z)\mathrm{\Gamma }_{q^2}({\displaystyle \frac{1}{2}})`$ $`=`$ $`[2]_q^{2z1}\mathrm{\Gamma }_{q^2}(z)\mathrm{\Gamma }_{q^2}(z+{\displaystyle \frac{1}{2}}),`$
(3) $`\theta _0(2z,\tau )`$ $`=`$ $`\theta _0(z,\tau )\theta _0(z+{\displaystyle \frac{1}{2}},\tau )\theta _0(z+{\displaystyle \frac{\tau }{2}},\tau )\theta _0(z+{\displaystyle \frac{1+\tau }{2}},\tau ).`$
The function $`\mathrm{\Gamma }(z)=_0^{\mathrm{}}t^{z1}e^t𝑑t`$ is the Euler gamma function. It satisfies the functional equation $`\mathrm{\Gamma }(z+1)=z\mathrm{\Gamma }(z)`$. Formula (1) is Legendre’s duplication formula.

The function $`\mathrm{\Gamma }_q(z)`$ is Jackson’s $`q`$-gamma function. Set $`x=e^{2\pi iz}`$, $`q=e^{2\pi i\tau }`$, and denote

$$(x;q)=\underset{j=0}{\overset{\mathrm{}}{}}(1xq^j).$$

Then

$$\mathrm{\Gamma }_q(z)=\mathrm{\Gamma }_{\mathrm{trig}}(z,\tau )=(1q)^{1z}\frac{(q;q)}{(q^z;q)}.$$

The $`q`$-gamma function obeys the functional equation

$$\mathrm{\Gamma }_q(z+1)=[z]_q\mathrm{\Gamma }_q(z),$$

where $`[z]_q=\frac{1e^{2\pi i\tau z}}{1e^{2\pi i\tau }}`$ is the trigonometric analog of the number $`z`$. The $`q`$-gamma function degenerates to Euler’s gamma function,

$$\underset{\tau 0}{lim}\mathrm{\Gamma }_{\mathrm{trig}}(z,\tau )=\mathrm{\Gamma }(z).$$

Formula (2) is Askey’s duplication formula, see \[A\].

The function $`\theta _0(z,\tau )=(x,q)(q/x,q)`$ in (3) is one of Jacobi’s theta functions. Formula (3) see for instance in \[Ra\].

In this paper we give two duplication formulas for elliptic analogs of the gamma function,

$`\mathrm{\Gamma }(2z,\tau ,\sigma )=`$ $`\mathrm{\Gamma }(z,\tau ,\sigma )\mathrm{\Gamma }(z+{\displaystyle \frac{\tau }{2}},\tau ,\sigma )\mathrm{\Gamma }(z+{\displaystyle \frac{\sigma }{2}},\tau ,\sigma )\mathrm{\Gamma }(z+{\displaystyle \frac{\tau +\sigma }{2}},\tau ,\sigma )`$
$`\mathrm{\Gamma }(z+{\displaystyle \frac{1}{2}},\tau ,\sigma )\mathrm{\Gamma }(z+{\displaystyle \frac{1+\tau }{2}},\tau ,\sigma )\mathrm{\Gamma }(z+{\displaystyle \frac{1+\sigma }{2}},\tau ,\sigma )\mathrm{\Gamma }(z+{\displaystyle \frac{1+\tau +\sigma }{2}},\tau ,\sigma ),`$
$`\overline{\mathrm{\Gamma }}(2z,\tau ,\sigma )\overline{\mathrm{\Gamma }}({\displaystyle \frac{1}{2}},2\tau ,\sigma )=\left({\displaystyle \frac{\theta _0(2\tau ,\sigma )}{\theta _0(\tau ,\sigma )}}\right)^{2z1}\overline{\mathrm{\Gamma }}(z,2\tau ,\sigma )\overline{\mathrm{\Gamma }}(z+{\displaystyle \frac{1}{2}},2\tau ,\sigma ),`$
see definitions below. The expression $`\frac{\theta _0(z\tau ,\sigma )}{\theta _0(\tau ,\sigma )}`$ is an elliptic analog of the number $`z`$. We have the trigonometric and rational limits of the theta function:

$$\frac{\theta _0(z\tau ,\sigma )}{\theta _0(\tau ,\sigma )}\stackrel{\sigma i\mathrm{}}{}\frac{1e^{2\pi i\tau z}}{1e^{2\pi i\tau }}\stackrel{\tau 0}{}z.$$

## 2. Elliptic gamma function

The elliptic gamma function is an elliptic generalization of the Euler gamma function. It is the meromorphic function of three complex variables $`z,\tau ,\sigma `$, with $`\mathrm{Im}\tau ,\mathrm{Im}\sigma >0`$ defined by the convergent infinite product

$$\mathrm{\Gamma }(z,\tau ,\sigma )=\underset{j,k=0}{\overset{\mathrm{}}{}}\frac{1e^{2\pi i((j+1)\tau +(k+1)\sigma z)}}{1e^{2\pi i(j\tau +k\sigma +z)}}.$$

It is the unique solution of a functional equation involving the Jacobi theta function $`\theta _0`$.

Theorem. \[FV1\] Suppose that $`\tau ,\sigma `$ are complex numbers with positive imaginary part. Then $`u(z)=\mathrm{\Gamma }(z,\tau ,\sigma )`$ is the unique meromorphic solution of the difference equation

$$u(z+\sigma )=\theta _0(z,\tau )u(z)$$

such that:

1. $`u(z)`$ obeys $`u(z+1)=u(z)`$ and is holomorphic on the upper half plane $`\mathrm{Im}z>0`$,
2. $`u((\tau +\sigma )/2)=1`$.

The elliptic gamma function first appeared in \[R\]. The modular properties of the elliptic gamma function and their relations to $`SL(3,)`$ are discussed in \[FV1\], appearances and applications of the elliptic gamma function can be found in \[B, DP, JMN, JKKMW, FTV, FV2, FV3, FV4\].

Let $`\overline{\mathrm{\Gamma }}`$ be the function

$$\overline{\mathrm{\Gamma }}(z,\tau ,\sigma )=\frac{(q;q)}{(r;r)}\theta _0(\tau ,\sigma )^{1z}\mathrm{\Gamma }(z\tau ,\tau ,\sigma ),q=e^{2\pi i\tau },r=e^{2\pi i\sigma }.$$

Then $`u(z)=\overline{\mathrm{\Gamma }}(z,\tau ,\sigma )`$ is a solution of the functional equation

$$u(z+1)=\frac{\theta _0(\tau z,\sigma )}{\theta _0(\tau ,\sigma )}u(z).$$

The normalization was chosen here so that $`u(1)=1`$. As $`\sigma i\mathrm{}`$ we recover Jackson’s $`q`$-gamma function,

$$\mathrm{\Gamma }_{\mathrm{trig}}(z,\tau )=\underset{\sigma i\mathrm{}}{lim}\overline{\mathrm{\Gamma }}(z,\tau ,\sigma ).$$

## 3. Multiplication Formulas

### 3.1. The first multiplication formula

Theorem. For any natural $`n`$ we have

$$\mathrm{\Gamma }(nz,\tau ,\sigma )=\underset{k_1,k_2,k_3=\mathrm{\hspace{0.17em}0}}{\overset{n1}{}}\mathrm{\Gamma }(z+\frac{k_1+k_2\tau +k_3\sigma }{n},\tau ,\sigma ).$$

Proof. Let $`w=e^{2\pi i/n}`$. Then the right hand side of this formula is

$`{\displaystyle \underset{k_1,k_2,k_3=\mathrm{\hspace{0.17em}0}}{\overset{n1}{}}}{\displaystyle \underset{l,m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1w^{k_1}q^{l+\frac{nk_2}{n}}r^{m+\frac{nk_3}{n}}x^1}{1w^{k_1}q^{l+\frac{k_2}{n}}r^{m+\frac{k_3}{n}}x}}`$
$`={\displaystyle \underset{k_2,k_3=\mathrm{\hspace{0.17em}0}}{\overset{n1}{}}}{\displaystyle \underset{l,m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1q^{(l+1)nk_2}r^{(m+1)nk_3}x^n}{1q^{ln+k_2}r^{mn+k_3}x^n}}={\displaystyle \underset{l,m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1q^{l+1}r^{m+1}x^n}{1q^lr^mx^n}}=\mathrm{\Gamma }(nz,\tau ,\sigma ).`$
$`\mathrm{}`$

### 3.2. The second multiplication formula

Theorem. For any natural $`n`$ we have

$`\overline{\mathrm{\Gamma }}(nz,\tau ,\sigma )\overline{\mathrm{\Gamma }}({\displaystyle \frac{1}{n}},n\tau ,\sigma )\overline{\mathrm{\Gamma }}({\displaystyle \frac{2}{n}},n\tau ,\sigma )\mathrm{}\overline{\mathrm{\Gamma }}({\displaystyle \frac{n1}{n}},n\tau ,\sigma )=`$
$`\left({\displaystyle \frac{\theta _0(n\tau ,\sigma )}{\theta _0(\tau ,\sigma )}}\right)^{nz1}\overline{\mathrm{\Gamma }}(z,n\tau ,\sigma )\overline{\mathrm{\Gamma }}(z+{\displaystyle \frac{1}{n}},n\tau ,\sigma )\overline{\mathrm{\Gamma }}(z+{\displaystyle \frac{2}{n}},n\tau ,\sigma )\mathrm{}\overline{\mathrm{\Gamma }}(z+{\displaystyle \frac{n1}{n}},n\tau ,\sigma ).`$
The theorem is an easy consequence of the following two lemmas.

Lemma. For any natural $`m`$ and $`n`$, we have

$$\mathrm{\Gamma }(z,\tau ,\sigma )=\underset{a=0}{\overset{m1}{}}\underset{b=0}{\overset{n1}{}}\mathrm{\Gamma }(z+a\tau +b\sigma ,m\tau ,n\sigma ).$$

Lemma. For any natural $`n`$, we have

$$\mathrm{\Gamma }(\tau ,n\tau ,\sigma )\mathrm{\Gamma }(2\tau ,n\tau ,\sigma )\mathrm{}\mathrm{\Gamma }((n1)\tau ,n\tau ,\sigma )=\frac{1}{(q,q^n)(q^2,q^n)\mathrm{}(q^{n1},q^n)}.$$


text:
# 1 Introduction.

## 1 Introduction.

String theory provides a consistent framework for investigating the physics of very high energies in general, and that of additional compact dimensions in particular. The physics of extra dimensions can be studied with some success on effective field theory grounds as well. Establishing an exact link of the string results with those of the effective field theory is however a difficult task. The latter theory may reproduce some results from string theory in the limiting case of an infinite string scale $`M_S`$ (zero slope $`\alpha ^{}1/M_S^2`$). It is this point of view that we adopt in our effective field theory analysis below, while seeking the exact relationship with string results for the problem of radiative corrections to the gauge couplings. However, the results we present are also relevant for models with “large” extra dimensions, without any reference to string theory.

Upon compactification, provided some conditions are fulfilled, remnant effects of the extra dimensions may be present and affect the low energy physics, such as the gauge couplings of the theory. For our discussion below we consider the example of a two-torus compactification. The model we address is an $`𝒩=1`$ supersymmetric orbifold with one $`𝒩=2`$ sub-sector (e.g. $`𝐙_4`$ orbifold). This two dimensional sub-sector (“bulk”), if charged under the gauge group of the theory, may bring one-loop corrections to the gauge couplings of the theory. At the string level these will be due to Kaluza-Klein and winding modes excitations with respect to the additional compact dimensions. These corrections are well-known in string case and were investigated in , see also e.g. , .

An effective field theory approach to investigating such corrections to the gauge couplings due to the two extra dimensions is also possible. Such an approach can only account for the effects of Kaluza-Klein states: the winding states effects cannot be described by a field theory approach. However, in the limit of a large compactification radius (compared to the ultraviolet length scale $`1/\mathrm{\Lambda }`$), one would hope that the effects of winding states are minimised. The reason for this would be that the above condition corresponds on the string side to the limit $`\alpha ^{}0`$ in which case the mass of the winding states is significantly larger than that of the momentum modes. Given the infinite number of these states it is not actually clear that they decouple in the low energy limit. Indeed, the density of momentum and winding modes is $`\rho =R+1/R`$, thus winding states may still affect the gauge couplings even if $`R`$ is very large  in string units. Comparing the effective field theory result (due to effects of momentum modes only) with that of the string (in its $`\alpha ^{}0`$ limit) will allow us understand whether the winding states have any UV effect on the gauge couplings. Previous effective field theory approaches  to such radiative corrections were restricted to understanding the UV behaviour of the couplings in the presence of extra dimensions. The exact link with string theory was recently addressed in .

We thus consider an $`𝒩=1`$ orbifold with a $`𝒩=2`$ sub-sector, and attempt to keep a general approach, without unnecessary model dependence. To begin with we note that the radiative corrections to the “bare” or string coupling $`\alpha _u`$, induced at one-loop level in the two-torus compactification are usually written as

$$\alpha _i^1(Q)=\alpha _u^1+\frac{b_i}{2\pi }\mathrm{ln}\frac{M_S}{Q}+\mathrm{\Delta }_i+\mathrm{},$$
(1)
Here $`Q`$ is a low energy scale, $`b_i`$ the $`𝒩=1`$ beta function coefficient, and $`M_S`$ is the ultraviolet scale (string scale). The first logarithm in eq. (1) accounts for the one loop effects of the massless states of the model. These states can include for example “twisted” string states (candidates for the MSSM-like matter fields), and also gauge bosons’ contributions. Computing their overall radiative effect on the gauge couplings requires an infrared (IR) and ultraviolet (UV) regularisation. The second term $`\mathrm{\Delta }_i`$ is due to massive Kaluza-Klein and, in string case winding states as well. The dots stand for higher order corrections and for mixing effects between the massless and $`𝒩=2`$ massive states sectors.

## 2 One-loop corrections from the string.

We first review some details of the string calculation for the one-loop corrections $`\mathrm{\Delta }_i`$ to the gauge couplings . At string level there exists an additional, gauge independent (universal) correction to the gauge couplings, which was “absorbed” into the (re)definition of the “bare” or string coupling $`\alpha _u`$. This ensures that at the string level this coupling is invariant under the symmetry $`SL(2,𝐙)_T\times SL(2,𝐙)_U\times 𝐙_2^{TU}`$ of the (heterotic) string .

The separation of the radiative effects on the gauge couplings into massless and massive modes contributions in (1) in the string calculation is not imposed by a string symmetry or principle. In general this is done because string calculations only compute the corrections due to massive modes alone. The expression for $`\mathrm{\Delta }_i`$ in string case (hereafter denoted $`\mathrm{\Delta }_i^H`$) is

$$\mathrm{\Delta }_i^H=\frac{\overline{b}_i}{4\pi }_\mathrm{\Gamma }\frac{d\tau _1d\tau _2}{\tau _2}\left(Z_{torus}1\right)$$
(2)
Here $`\overline{b}_i`$ is the beta function coefficient associated with the $`𝒩=2`$ sub-sector of the $`𝒩=1`$ orbifold. $`Z_{torus}`$ is the two-torus string partition function. The modulus $`\tau =\tau _1+i\tau _2`$ of the world sheet torus is integrated over the fundamental domain $`\mathrm{\Gamma }=\{\tau _2>0,|\tau _1|<1/2,|\tau |>1\}`$. As the massless states contribution is added separately in eq. (1), the massless (“zero”) mode contribution has been subtracted out in (2) and accounted for by “-1”. One can show that the integral over $`Z_{torus}`$ alone is modular invariant but is (infrared) divergent. However, by subtracting the massless contribution the result is (made) finite, but is not modular invariant anymore due to this last contribution. Briefly, modular invariance does not ensure a finite string result. Finally, there seems to be no reason to integrate this massless modes contribution “-1” over the fundamental domain $`\mathrm{\Gamma }`$, except to make the string result finite. Quantitatively, $`\mathrm{\Delta }_i^H`$ can be written as

$$\mathrm{\Delta }_i^H=\frac{\overline{b}_i}{4\pi }_\mathrm{\Gamma }\frac{d\tau }{\tau _2}\left[\frac{T_2}{\tau _2}\underset{A}{}e^{2\pi iTdetA}\mathrm{exp}\left[\frac{\pi T_2}{\tau _2U_2}\left|(1U)A(\begin{array}{c}\tau \\ 1\end{array})\right|^2\right]1\right];A=(\begin{array}{cc}n_1& p_1\\ n_2& p_2\end{array})$$
(3)
where integers $`p_{1,2}`$ are (Poisson re-summed) Kaluza-Klein levels and $`n_{1,2}`$ are winding modes. The moduli fields $`T=T_1+iT_2`$, $`U=U_1+iU_2`$ can be expressed in terms of the anti-symmetric tensor background $`B_{ij}=Bϵ_{ij}`$ and the radii $`R_{1,2}`$ and angle $`\theta `$ of the two-torus as $`T=2[B+iR_1R_2\mathrm{sin}\theta /(2\alpha ^{})]`$, $`U=R_2/R_1\mathrm{exp}(i\theta )`$.

The integral in eq.(3) can be written as a sum over the “orbits” of the modular group $`SL(2,𝐙)`$. This is just a sum over classes of matrices $`A`$ with entries integer numbers labelling (Poisson re-summed) Kaluza-Klein and winding levels, giving

$$\mathrm{\Delta }_i^H=\frac{\overline{b}_i}{4\pi }\left[𝒥^{(A=0)}+𝒥^{(detA=0)}+𝒥^{(detA0)}+_\mathrm{\Gamma }\frac{d\tau _1d\tau _2}{\tau _2}(1)\right]_{reg}$$
(4)
The subscript “reg” indicates that the splitting of the integrals in eqs.(2), (4) only makes sense in the presence of an infrared regulator. For example, each of the separate terms can be multiplied by $`(\tau _2)=(1\mathrm{exp}(N/\tau _2))`$ with $`N\mathrm{}`$. This regularisation breaks modular invariance, but other IR regularisation schemes, which are modular invariant  can be used.

The three (regularised) contributions in eq.(4) are classified in function of the matrices $`A`$ and correspond to: the zero orbit $`𝒥^{(A=0)}`$, the degenerate orbit $`𝒥^{(detA=0)}`$, and the non-degenerate orbit $`𝒥^{(detA0)}`$. $`𝒥^{(A=0)}`$ is due to infinitely many original Kaluza-Klein modes. $`𝒥^{(detA=0)}`$ is due to a mixing of momentum and winding modes (if $`n_1=p_1=0`$), or momentum modes alone (if $`n_{1,2}=0`$) or winding modes alone (if $`p_{1,2}=0`$). In the limit of an infinite string scale or $`\alpha ^{}0`$ the momentum modes contribution is dominant. $`𝒥^{(detA0)}`$ is due to a mixing of momentum and winding modes and is vanishing in the limit of an infinite string scale or $`\alpha ^{}0`$ when winding modes’ effects are suppressed. Quantitatively, one has

$`𝒥^{(A=0)}`$ $`=`$ $`{\displaystyle _\mathrm{\Gamma }}{\displaystyle \frac{d\tau _1d\tau _2}{\tau _2^2}}T_2={\displaystyle \frac{\pi }{3}}T_2,`$ (5)
$`𝒥^{(detA=0)}`$ $`=`$ $`\mathrm{ln}N\mathrm{ln}\left[4\pi e^{2\gamma _E}T_2U_2|\eta (U)|^4\right]+𝒪\left((T_2U_2/N)^{\frac{1}{2}}\right)`$ (6)
$`𝒥^{(detA0)}`$ $`=`$ $`\mathrm{ln}{\displaystyle \underset{n_1=1}{\overset{\mathrm{}}{}}}\left|1e^{2\pi in_1T}\right|^4,𝒥^{(detA0)}0\text{ if }\alpha ^{}0`$ (7)
$`{\displaystyle _\mathrm{\Gamma }}{\displaystyle \frac{d\tau _1d\tau _2}{\tau _2}}(1)|_{reg.}`$ $`=`$ $`\mathrm{ln}N+\mathrm{ln}3\sqrt{3}e^{1\gamma }/2,`$ (8)
The term $`𝒪((T_2U_2/N)^{1/2})`$ is just a correction depending on the infrared regulator and vanishes in the limit of removing it, $`N\mathrm{}`$. Adding together the above equations, one obtains the final result

$$\mathrm{\Delta }_i^H=\frac{\overline{b}_i}{4\pi }\mathrm{ln}\left\{C_{reg}U_2|\eta (U)|^4T_2\left|\eta \left(T\right)\right|^4\right\}+𝒪\left((T_2U_2/N)^{\frac{1}{2}}\right),$$
(9)
where $`C_{reg}=8\pi e^{1\gamma }/(3\sqrt{3})`$ is a regularisation scheme dependent constant (for the DR scheme , $`C_{reg}=4\pi e^\gamma `$).

## 3 One-loop corrections from a field theory approach.

A field theory approach to computing the effects of the two extra dimensions of the two-torus on the gauge couplings can only sum the effects due to Kaluza-Klein states (no winding modes). The calculation we present is also relevant in the context of models with “large” compactification radii (relative to the UV cut-off length) and without any reference to the string. As the massless states’ effects are already accounted for in eq.(1) the one-loop radiative corrections to the gauge couplings due to massive modes alone is (hereafter denoted $`\mathrm{\Delta }_i^{}`$)

$$\mathrm{\Delta }_i^{}=\frac{1}{4\pi }\underset{i}{}T(R_i)\underset{m_{1,2}𝐙}{\overset{}{}}_\xi ^{\mathrm{}}\frac{dt}{t}e^{\pi tM_{m_1,m_2}^2/\mu ^2},$$
(10)
The result of summing a finite or infinite tower of Kaluza-Klein states to the gauge couplings is divergent. Since the integral above is divergent in the UV ($`t0`$), we introduced a UV dimensionless (proper time) regulator $`\xi 0`$ as the lower limit of the integral, and an arbitrary (finite) mass scale $`\mu `$. A “prime” on the sum indicates that it does not include the effects of the massless state $`(m_1,m_2)(0,0)`$, which is already present in (1). Finally, the mass of the Kaluza-Klein states is

$$M_{m_1,m_2}^2=\frac{1}{\mathrm{sin}^2\theta }\left[\frac{m_1^2}{R_1^2}+\frac{m_2^2}{R_2^2}\frac{2m_1m_2\mathrm{cos}\theta }{R_1R_2}\right]=\frac{|m_2Um_1|^2}{(\mu ^2T_2(\mu ))U_2},$$
(11)
where $`T(\mu )iT_2(\mu )=i\mu ^2R_1R_2\mathrm{sin}\theta `$ and the complex moduli $`U`$ is identical to that in the string case, $`U=R_2/R_1\mathrm{exp}(i\theta )`$. For the particular case of an orthogonal torus $`\theta =\pi /2`$ one recovers a more familiar mass formula: $`M_{m_1,m_2}^2=m_1^2/R_1^2+m_2^2/R_2^2`$, and if the radii are equal one obtains that $`M_{m_1,m_2}^2=(m_1^2+m_2^2)/R_2^2`$.

The calculation of the integral for $`\mathrm{\Delta }_i^{}`$ is rather technical, so we only quote the final result

$`\mathrm{\Delta }_i^{}`$ $`=`$ $`{\displaystyle \frac{\overline{b}_i}{4\pi }}\mathrm{ln}\left[4\pi e^{\gamma _E}e^{T_2^{}}T_2^{}U_2|\eta (U)|^4\right],`$ (12)
$`T_2^{}`$ $``$ $`\mathrm{\Lambda }^2R_1R_2\mathrm{sin}\theta ,\mathrm{\Lambda }^2\mu ^2/\xi ,\mathrm{max}\{1/(R_2\mathrm{sin}\theta ),1/R_1\}\mathrm{\Lambda }`$ (13)
where $`\overline{b}_i`$ is a sum over the beta function contributions of those states with a Kaluza-Klein tower and $`\mathrm{\Lambda }`$ is the UV scale associated with the UV regulator $`\xi `$. This result for $`\mathrm{\Delta }_i^{}`$ holds true only in the limit of “removing” the dependence on the regulator $`\xi 0`$, when additional corrections in $`\xi `$ vanish. This condition is converted into bounds on the mass scales of the theory presented in (13). Notice that $`R_2\mathrm{sin}\theta `$ plays the role of an effective radius.

## 4 A comparison with string theory corrections.

We can now compare the field theory result (12) with the heterotic string expression (9) for the gauge thresholds. For this comparison to make sense we note that one should actually compare the field theory results with the limit $`\alpha ^{}0`$ of the string results (9), when the effects of the winding modes, not included in the field theory calculation, are suppressed.

We observe that the $`U`$ dependence of the two equations is identical. As both the effective field theory and the string theory mass spectra are $`SL(2,𝐙)_U`$ symmetric, and $`U`$ is by definition $`\alpha ^{}`$ independent (i.e. scale independent), it may not be surprising that the $`U`$ dependence of the final string result can be entirely recovered on effective field theory grounds.

Regarding the $`T`$ dependent part of the two results (9), (12) we note the following. The field theory UV cut-off in (12) can be identified with the string scale: $`\mathrm{\Lambda }^2\mu ^2/\xi 1/\alpha ^{}`$. Removing the regulator, $`\xi 0`$ on the field theory side in (10) corresponds to an infinite string scale or $`\alpha ^{}0`$. In this limit eq. (9) (see also (5), (6)) has quadratic and logarithmic divergences in the string scale similar to those in eq.(12). However, the coefficient of the quadratic divergence is different: $`T_2^{}`$ in $`\mathrm{\Delta }_i^{}`$, and $`(\pi /3)T_2`$ in $`\mathrm{\Delta }_i^H`$. In the field theory case this coefficient is regulator dependent, thus it can be chosen such that this difference is avoided. Such choice provides the appropriate definition for the UV cut-off scale of an effective field theory, compatible with the modular invariance of its string embedding. We conclude that the full UV structure of the string thresholds in the limit $`\alpha ^{}0`$ can be obtained on pure field theory grounds, except that the coefficient of the quadratically divergent part is arbitrary in the field theory approach.

To explain how the factor $`\pi /3`$ arises in the heterotic string, and its link with winding modes effects, one should analyse its origin in eq. (5). This factor is consequence of the integration over the fundamental domain $`\mathrm{\Gamma }`$, and is thus an effect of modular invariance symmetry and indirectly, of the winding modes. As this symmetry is not present in the effective field theory, the factor $`\pi /3`$ cannot be recovered. Thus winding modes indirectly control the coefficient of the term quadratic in string scale ($`T_2M_S^2`$), therefore they do have an UV role, even in the limit $`\alpha ^{}0`$.

We would like to end with an additional remark on the field theory limit $`\alpha ^{}0`$ of the string result, eq. (9). A closer look at the final (infrared regularised) string result of eq. (9) shows that higher corrections due to the infrared regularisation $`𝒪((U_2T_2/N)^{1/2})`$ are discarded in the final result for $`\mathrm{\Delta }_i^H`$ when $`N\mathrm{}`$. This is true if $`T_2`$ is finite. However, in the field theory limit of $`\alpha ^{}0`$ or equivalently $`T_2\mathrm{}`$ (with $`T_2`$ expressed in string units) the term $`𝒪((U_2T_2/N)^{1/2})`$ is not vanishing. Therefore this term should be kept in the field theory limit of the string result $`\mathrm{\Delta }_i^H`$. This issue requires an investigation at the string level, to clarify if any string symmetry left after compactification can still impose the order to take the limits of removing the infrared regulator $`N\mathrm{}`$ and the field theory limit $`T_2\mathrm{}`$ on the string result. This is relevant because these two limits do not commute. For further discussions on this and its relationship to the field theory approach see recent results in .

## 5 Conclusions.

We provided an effective field theory approach to computing the corrections to the gauge couplings due to Kaluza-Klein states and a comparison with the string result which also includes the effects of the winding states.

While the heterotic string calculation is well known and studied in the literature, we attempted to emphasize some points which are usually overlooked: one is that that the finiteness of the string threshold correction $`\mathrm{\Delta }_i^H`$ is not due to modular invariance. It is the massless modes’ contribution which is introduced to “make” the result finite, by subtracting the (infrared) divergence due to the massive momentum and winding modes. The integration of the massless modes contribution (non-modular invariant) over the fundamental domain remains a procedure to make the string result finite, but is not required or supported by a string symmetry. One should clarify at the string level if any string symmetry can dictate the order to take the limits of removing the string infrared regulator and the field theory limit, since these two limits do not commute.

The effective field theory approach to computing corrections to the gauge couplings due to the two extra dimensions of the torus provides results similar to those of the string in the limit of an infinite string scale. The full dependence on the shape moduli ($`U`$) is recovered. In addition, the UV divergences of the string result in the limit $`\alpha ^{}0`$ were reproduced by the field theory calculation. The coefficient in front of the leading UV divergences in string theory (in the limit $`\alpha ^{}0`$) turns out to be equal to $`\pi /3`$ as a result of modular invariance symmetry, while at field theory level this coefficient is regularisation scheme dependent.

## Acknowledgements.

This work was supported by PPARC UK (D.G.), and in part by CITA and NSERC (S.G.N.).


text:
# Ac Susceptibility and Static Magnetization Measurements of CeRu2Si2 at Small Magnetic Fields and Ultra Low Temperatures

## Abstract

The magnetic properties of CeRu2Si2 at microkelvin temperatures (down to 170 $`\mu `$K) and ultra small magnetic fields ($`0.026.21`$ mT) are investigated experimentally for the first time. The simultaneously measured ac susceptibility and static magnetization show neither evidence of the magnetic ordering, superconductivity down to the lowest temperatures nor conventional Landau Fermi-Liquid behavior. The results imply the magnetic transition temperature in undoped CeRu2Si2 is very close to absolute 0 K. The possibility for proximity of CeRu2Si2 to the quantum critical point without any doping is discussed.

preprint: APS/123-QED

The unusual properties of heavy fermion (HF) systems are determined by the competition between intersite spin couplings, Ruderman-Kittel-Kasuya-Yosida interaction, and intrasite Kondo interaction. Doniach  In a system dominated by the Kondo effect, the Pauli paramagnetic (PP) state with massive quasiparticles is achieved through screening of the *f* electron’s magnetic moments by conduction electrons below the characteristic temperature $`T_\text{K}`$. The physical properties of the HF compounds below $`T_\text{K}`$ are well understood within the framework of the Landau Fermi-liquid (LFL) theory.

Recently, however, non-Fermi-liquid (NFL) behavior was observed in a large class of HF compounds near the quantum critical point (QCP). QCP1 ; QCP2  NFL systems exhibit anomalous temperature dependence of the physical quantities in contrast to the LFL theory, such as specific heat $`\mathrm{\Delta }C/T\mathrm{ln}T`$, resistivity $`\mathrm{\Delta }\rho T^ϵ(1ϵ<2)`$ and magnetic susceptibility $`\mathrm{\Delta }\chi `$ either $`1\sqrt{T}`$ or $`\mathrm{ln}T`$. In general, the quantum (zero-temperature) phase transition is driven by a control parameter other than temperature, for example, composition, pressure, or magnetic field, and is accompanied by a qualitative change in the correlations in the ground state. The second order quantum phase transitions and QCPs in HF systems can be classified into two types: (i) the long-wavelength fluctuations of the order parameter are the only critical degrees of freedom and the quantum criticality is developed as spin-density wave instability, SCRHF ; Millis  here the zero-temperature spin fluctuations are given by the Gaussian fluctuations of the order parameter; (ii) local critical modes coexist with long-wavelength fluctuations of the order parameter and there is non-Gaussian distribution of the fluctuations. Si  These are the so-called locally critical phase transitions where the quantum criticality of CeCu(6-x)Aux (Ref. CeCu6Au, ) and YbRh2Si2 (Ref. YbRhSi1, ) are regarded as type-(ii) QCP. Si

CeRu2Si2 with a ThCr2Si2-type crystal structure is well known to be a typical HF compound with an electronic specific-heat coefficient $`\gamma =350\text{mJ/K}^2\text{mol}`$ below $`T_\text{K}=20`$ K. Gamma ; KondoTemp  This compound exhibits the pseudo-metamagnetic transition into the ferromagnetically ordered state induced by the magnetic field at $`H_\text{M}=7.8`$ T below 10 K. meta1 ; meta2 ; meta3 ; meta4 ; meta5  The neutron-scattering measurements note short-range antiferromagnetic (AFM) correlations in CeRu2Si2 even below $`T_\text{K}`$. These time-fluctuating correlations are described by different incommensurate wave vectors. AFM1 ; AFM2  A $`\mu `$SR experiment shows an ultra small static moment of the order of 10$`{}_{}{}^{3}\mu _{\text{B}}^{}`$mSR  It is most remarkable that the alloying compound systems Ce(1-x)LaxRu2Si2 and Ce(Ru(1-x)Rhx)2Si2 show an incommensurate spin density wave (SDW) ground state in a concentration range $`x>0.08`$ (Ref. La1, ) and $`0.03<x<0.4`$ (Ref. Rh1, ). This long-range ordering has the form of a sine-wave modulated structure with the short-range correlation of CeRu2Si2, as described above. At the critical concentrations of $`x_c=0.075`$ for La and $`x_c=0.03`$ for Rh doping, the SDW transition vanishes but $`T_\text{K}`$ remains a finite temperature. Therefore, these small critical concentrations of La and Rh suggest that CeRu2Si2 might be located in the vicinity of the type-(i) QCP, and AFM spin fluctuations are expected to play a key role in this magnetic ground state. The $`C`$, $`\chi `$ and $`\rho `$ measurements show the conventional LFL ground state for CeRu2Si2 below $`T_\text{K}`$ down to 20 mK. meta1  All magnetization and susceptibility measurements have been performed at the magnetic field above 1 T and at low temperatures. In this letter, we report for the first time the results of ac susceptibility and static magnetization measurements at microkelvin temperatures (down to 170 $`\mu `$K) and ultra small magnetic fields ($`0.026.21`$ mT). The obtained magnetic field and temperature dependence of the susceptibility and magnetization provide evidence of NFL behavior and allow us to think about the proximity of CeRu2Si2 to QCP.

The single crystal of CeRu2Si2 was prepared by a Czochralski pulling method with starting materials Ce ($`99.99`$%), Si ($`>99.999`$%) and Ru ($`99.99`$%) and purified by a solid state transport method. The sample size was $`11\times 4.2\times 1.5`$ mm3. The sample was cooled with a copper nuclear demagnetization refrigerator and a 3He-4He dilution refrigerator. It was sandwiched between two silver plates which were parts of the thermal link to the copper nuclear stage. The temperature was measured by a Pt NMR thermometer, a 3He melting curve thermometer MeltingCurve  and a RuO2 resistance thermometer. All these thermometers and the thermal link were attached to the same flange of the nuclear stage. The temperature difference between sample and heat bath was estimated to be less than the order of 0.1 % at all temperatures.

The ac susceptibility and static magnetization of CeRu2Si2 were measured simultaneously in a static field $`0.02B6.21`$ mT by an ac impedance bridge using a SQUID magnetometer. The applied static field declined a few tens of degrees from the crystalline *c*-axis. All of the ac susceptibility measurements were performed at a frequency of 16 Hz with an excitation field below $`0.75`$ $`\mu `$T parallel to the static field. The primary coil, secondary coil, and static field coil were placed inside a Nb superconducting magnetic shield, which was surrounded by a $`\mu `$ metal magnetic shield to suppress any external stray field. The static magnetization was calibrated against the absolute value measured by another magnetometer in the temperature range from 4 to 2 K.

The ac susceptibility was measured during cooling and warming, and the results showed no appreciable hysteresis. The static magnetization was measured in the warming procedure. Figure 1 shows the temperature dependence of the inphase ($`\chi ^{}`$) and quadrature ($`\chi ^{\prime \prime }`$) components of the ac susceptibility ($`M/H`$) at different magnetic fields below 400 mK. The inset of Fig.1 shows the temperature dependence of the static susceptibility ($`M/H`$) derived from the static magnetization. We calibrated all data against the temperature independent PP susceptibility which observed above $``$50 mK. Below $``$50 mK, we observed an excess susceptibility obeying the Curie law.

The ac susceptibility shows a peak at the magnetic field between 0.20 and 0.94 mT. The peak temperatures $`T_\text{P}`$ shift to higher values and the ac susceptibility is suppressed with increase in the applied magnetic field. In particular, the ac susceptibility at 6.21 mT is suppressed almost to the level of the PP susceptibility. The ac and static susceptibilities deviate from the Curie law as they approach $`T_\text{P}`$. The Curie constant $`C`$ can be written in the form $`C=N_\text{A}\mu _\text{0}\mu _\text{p}^2/3k_\text{B}V_{\text{mol}}`$. The effective magnetic moment $`\mu _\text{p}`$ turns out to be $`0.020\pm \mathrm{\hspace{0.17em}0.003}\mu _\text{B}`$/unit cell from the static susceptibility and is independent of applied magnetic fields. The value of $`\mu _\text{p}`$ is in agreement with the ultra-small static moment observed in the $`\mu `$SR experiment. mSR  The static susceptibility, however, becomes flat with no peaks in fields higher than 0.20 mT.

Figure 2 shows the $`H/T`$ dependence of the static magnetization below 400 mK. The dashed line corresponds to the Curie law with $`\mu _\text{p}=0.02\mu _\text{B}`$. In the fields above 0.20 mT, the magnetization clearly shows the saturation. The saturated magnetic moment $`\mu _\text{s}`$ can be evaluated in each field using the following relation : $`M_\text{s}=N_\text{A}\mu _\text{s}/V_{\text{mol}}`$. The calculated $`\mu _\text{s}`$ are $`1.20\times 10^5,1.60\times 10^5,1.95\times 10^5`$ and $`2.98\times 10^5\mu _\text{B}`$/unit cell, and the ratio of $`\mu _\text{p}`$ to $`\mu _\text{s}`$ is derived as $`1.80\times 10^3,1.15\times 10^3,0.98\times 10^3`$ and $`0.73\times 10^3`$ at 0.20, 0.39, 0.94 and 6.21 mT, respectively. Figure 2 suggests that the magnetization cannot be attributed to an impurity effect. At very small concentrations, the magnetic impurities contribution to the total magnetization should behave according to a Brillouin function. With increasing concentration, a locally ordered state like spin-glass can be formed. The impurity effect on the static magnetization in CeCu$`_\text{6}`$ at low temperatures is one example for the first case. CeCu6static  It indicates that the ratio of $`\mu _\text{p}`$, as deduced by Curie law, to $`\mu _\text{s}`$ has to be of the order of 1. However, this does not agree with our results, $`\mu _\text{p}/\mu _\text{s}10^3`$. In the second case, many compounds with a spin glass transition show quite different magnetization behavior between zero field cooling (ZFC) and field cooling (FC) measurements through the transition temperature. The inset in Fig. 2 shows the magnetization at 0.94 mT measured by ZFC. The results of ZFC and FC do not indicate different behavior below $`T_\text{P}`$ in this field. Consequently, the possibility of spin glass transition is also strongly denied.

Below $`T_\text{P}`$, there is unique scaling behavior of the ac susceptibility and static magnetization from 0.20 to 0.94 mT, as shown in Fig. 3. We normalized the susceptibility by its peak height and the magnetization by the saturated value at each field. And the temperature is normalized by $`T_\text{P}`$. The scaling behavior provides evidence of the proximity of CeRu2Si2 at small magnetic fields and ultra low temperatures to some phase transition. S-K.Ma  The nature of such transition is determined by interplaying between ferromagnetic (FM) and AFM fluctuations observed by neutron scattering experiments. meta4 ; meta5 ; AFM1 ; AFM2  And this type of interplaying was evidently observed around $`H_\text{M}`$meta4  In order to shed light on the nature of this transition we analyzed this scaling behavior from the viewpoint of up-to-date theories for ordering in HF compounds.

According to the mean-field theory, the temperature dependence of the susceptibility with the SDW ground state below Néel temperature obeys the following expression : $`\chi (T)=\chi _0+B\mathrm{exp}(a/T)`$, where $`\chi _0`$ and $`a`$ are the $`\chi (T0)`$ and an energy gap, respectively. There, the ratio of the gap energy to $`T_\text{N}`$ should be above 1.76. In our case, the normalized temperature $`T^{}`$ dependence of the normalized ac susceptibility $`\chi ^{}(T^{})`$ is obeying above the expression with $`a/T^{}1.0`$ at the peak temperature, as shown in Fig. 3. This ratio is in contrast to the SDW state case and there is no indication for the AFM transition at $`T_\text{P}`$. The scaling behavior and the exponent type temperature dependence of $`\chi ^{}`$ below $`T_\text{P}`$, however, suggest that CeRu2Si2 is in some magnetic field arranged state between 0.20 and 0.94 mT. We speculate that the physical background of this scaling behavior is the quantum critical fluctuation effect of CeRu2Si2 which is in proximity to the QCP discussed below.

Further, we compare CeRu2Si2 with CeCu6 which is also a typical HF compound and very similar to CeRu2Si2 in its 4*f* electron behavior below $`T_\text{K}`$. Recently, it has been shown that magnetic ground state of CeCu6 is the SDW state and $`T_\text{N}`$ is in fair agreement with the estimated one by the self consistent renormalization (SCR) spin fluctuation theory. Tsujii  The SCR theory for the HF system predicts the value of $`T_\text{N}`$ as an equation : $`T_\text{N}=0.1376p_\text{Q}^{4/3}T_A^{2/3}T_0^{1/3}`$, where $`p_\text{Q}`$ is the staggered spontaneous moment in $`\mu _B`$ at $`T=0`$ K, $`T_\text{A}`$ and $`T_\text{0}`$ are the characteristic temperatures in the *q* and $`\omega `$ space, respectively. SCRHF  If we use the derived values of $`T_\text{A}=16`$ K, $`T_\text{0}=14.1`$ K (Ref. CharacteristicTemp, ) and $`p_\text{Q}=7\times 10^3`$ $`\mu _B`$ (Ref. mSR, ) for CeRu2Si2, the predicted $`T_\text{N}`$ is estimated as $`T_\text{N}2.8`$ mK. However, the magnetic properties of CeRu2Si2 do not show any indication of the magnetic ordering in the smallest applied magnetic field. The estimation of $`T_\text{N}`$ for SDW state in framework of SCR theory takes into account only AFM characteristic wave vectors while our data indicate an existence of both FM and AFM fluctuations. And these two types of fluctuations can lead to some disordered state in the smallest field. This means that the magnetic transition temperature in CeRu2Si2 is possibly close to $`T=0`$ K and the spin system of CeRu2Si2 under our conditions is in the vicinity of the QCP.

In the case of type-(i) QCP, we cannot explain the temperature dependence of the magnetic susceptibility in CeRu2Si2 only from AFM fluctuations which have been observed in neutron experiments. AFM1 ; AFM2  Based on the SCR theory, the uniform susceptibility of an itinerant weak AFM compound should not indicate the Curie-Weiss behavior without taking account of the FM fluctuations. SCRFerro  The large value of the $`\mu _\text{p}`$/$`\mu _\text{s}`$ also suggests the weak FM properties in this compound.

We consider two known examples of type-(ii) QCP for the 4*f* electron system. The well-known type-(ii) QCP doped compound is CeCu$`_{\text{5.9}}`$Au$`_{\text{0.1}}`$, which exhibits the critical scaling behavior for the differential susceptibility ($`M/H`$) in the form : $`(M/H)^1=\chi _0^1+CT^\alpha g(H/T)`$. Here $`\alpha `$ is the critical exponent and the universal scaling function $`g(H/T)`$ is given by Schröder *et al*.. CeCu6Au  This scaling function, however, does not lead to the peak for $`M/H`$ ; for this reason, it does not explain our results. On the other hand, our results are very similar to the case of undoped YbRh2Si2, which is classified as a type-(ii) QCP compound. A plateau in the Knight shift below 1 K is observed in YbRh2Si2, while the ac susceptibility shows a peak at the marginal temperature with the magnetic field dependence. YbRhSi1 ; YbRhSi2 ; YbRhSi3  The Curie-Weiss behavior of magnetic susceptibility $`\chi (T)`$ for YbRh2Si2, however, hints to large fluctuating localized Yb3+ moments, while that for our system indicates very tiny fluctuating moments with itinerant nature.

We speculate that the magnetic properties of CeRu2Si2 at small magnetic fields and ultra low temperatures are determined by competition between FM and AFM fluctuations. The narrow range of the applied magnetic fields above 0.20 mT can modulate the FM fluctuations and show a magnetic field arranged state below $`T_\text{P}`$. The FM fluctuations, however, decrease with increasing magnetic field because the magnetic ground state is recovered nearly to the LFL state at the field above 6.21 mT. In the field of 0.02 mT, the FM and AFM fluctuations compete strongly and show a non-analytical temperature dependence. This possible scenario of the QCP in our compound shows that CeRu2Si2 at small magnetic fields and ultra low temperatures can probably be considered one of the candidates for investigations of quantum phase transitions at ambient pressure and without any doping. Similar to the type-(ii) QCP compound YbRh2Si2, the external magnetic field is the control parameter for that transition. The NMR measurements for CeRu2Si2 under the conditions described above would be very useful in identifying of the nature of the itinerant 4*f* electron system and that QCP.

This work was partially supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan. The first author is grateful to JSPS for a Research Fellowship for Young Scientists. The authors are also grateful to H. Tsujii and Y. Tabata for valuable discussions and to K. Mukai, T. Tsunekawa and K. Nunomura for technical support.


text:
# Ab Initio Molecular Dynamics Study of Aqueous Solvation of Ethanol and Ethylene

## I Introduction

The study of the solvation of alcohols in aqueous solution is of fundamental interest in physics, biology and chemistry, but also of importance in technical applications FranksIves ; Eisenberg ; Franks ; FranksReid ; Blandamer ; Luck ; Franksv4 ; Alfsen ; Franks2 ; FranksDes ; Mohr . Among the alcohols, ethanol is one the most well studied compounds. Aqueous ethanol solutions are a common in chemical research and industry applications. Ethanol can be produced by aqueous hydration of ethylene, that is readily available from natural sources. This process that can be accelerated by acid catalysis. The reverse route of decomposing of ethanol into water and ethylene is also of importance. For example, for developing countries who do not have a large supply of fossil fuels, dehydration of ethanol obtained from biomass is often the most economical way to produce ethylene. As is well known, solvation structures play a crucial role in aqueous solution chemistry, where reactive events often require a significant reordering of the water molecules in the solvation shell. The solvation of ethanol and ethylene are therefore crucial in the course of their (acid-catalyzed) interconversion. The aqueous solvation of these molecules will be addressed in the present paper. Elsewhere we report on an ab initio molecular dynamics study of the acid-catalyzed conversion ownfut .

Ethanol and ethylene have distinct solvation properties in aqueous solution. Ethanol is easily soluble as its polar hydroxyl group can participate in the hydrogen bonded network and the hydrophobic ethyl group is relatively small. In contrast, the apolar ethylene molecule has a much weaker interaction with water and is generally considered to be hydrophobic. Mixtures of water and ethanol have been studied extensively, both experimentally as by molecular simulation. Experimental studies employing NMR Coccia ; Harris ; Mizuno ; Ludwig ; Lamanna , ultrasonic absorption Brai , infrared absorption spectroscopy Harris ; Mizuno ; Nishi ; Onori , mass spectroscopy Nishi ; Nishi0 , X-ray diffraction measurements Nishi ; Petrillo , neutron diffraction Petrillo ; Turner ; Sidhu , and dielectric relaxation measurements Mashimo ; Bao ; Sato ; Petong  have been performed to unravel the solvation properties of ethanol. Molecular simulation studies using empirical force fields have addressed the equation of state, thermodynamics, and structure and dynamics of solvation Mueller-Plathe ; Levchuk ; Fidler ; Tarek ; Tarek2  of aqueous ethanol solutions. The following general picture of the aqueous solvation in dilute solutions has emerged: the hydroxyl group participates in the hydrogen bonded network, while the hydrophobic alkyl group is accommodated in the hydrogen bonded network of water molecules. The nature of hydration structure around the hydrophobic part of alcohols is still a controversial subject. It has been suggested FranksIves ; Ludwig ; Lamanna ; Franks ; FranksReid ; Franksv4 ; Blandamer ; FranksDes  that hydrophobic solutes enforces the network of hydrogen bonded water molecules around it and decreases their mobility, a notion referred to as hydrophobic hydration Frank ; Ben-Naim . However, recent experimental and computational studies have shown that for small alcohols the structure and dynamics of the water molecules in the solvation shell is almost identical to that of bulk water Fidler ; Turner ; Sidhu .

Much less is known for the solvation of ethylene in water under ambient conditions. Experimental Sugahara  and theoretical Tanaka ; Kvamme  work addressed the clathrate hydrates of ethylene in water. The isolated ethylene-water complex has been a subject of various studies. In the lowest energy configuration, the ethylene molecule forms a weak bond with water. An early *ab initio* study of Del Bene in 1974 DelBene  has characterized this interaction as a $`\pi `$-H hydrogen bond of a water proton with the $`\pi `$ electrons of the C=C bond. Several experimental techniques have been applied to measure the strength of this interaction, such as the matrix isolation study of Engdahl and Nelander Engdahl ; Engdahl86 , the microwave spectra study of Andrews and Kuczkowski Andrews  and the molecular-beam measurements of Peterson and Klemperer Peterson . Recently Tarakeshwar et al. Tarakeshwar ; Tarakeshwar2 ; Tarakeshwar3  and Dupré and Yappert Dupre  have performed calculations on the ethylene-water complexes with high level *ab initio* methods. The interaction is weak compared to a hydrogen bonds such as in the water dimer. The role of the $`\pi `$-H bond in aqueous solvation under ambient conditions is still an open question.

Molecular simulation provides an approach to study the microscopic behavior of liquids complementary to experimental studies. All molecular simulations studies of aqueous ethanol solutions up to now are based on empirical force fields that are designed to reproduce a selection of experimental data. Obviously, molecular simulations based on these potentials do not provide a picture completely independent from experiment. Moreover, the reliability of the results at conditions that are significantly different from those where the potential was designed for, may be questionable. Density functional theory (DFT) based molecular dynamics (MD) simulation, such as the Car-Parrinello molecular dynamics method CaPa85 , where the interactions are calculated by accurate electronic structure calculations, provides a route to overcome these limitations. This has been demonstrated in studies of liquid water Laasonen ; Sprik ; SiBe97  and aqueous solvation Marx ; ErpMe01\_erratum ; Raugei . Important advantages of DFT-MD over force-field MD are that it intrinsically incorporates polarization, that it accounts for the intra-molecular motion and therefore allows for a direct comparison with spectroscopy of intra-molecular vibrations, and that it yields detailed information on the electronic properties, such as the energy levels of electronic states and the charge distribution. In a broader chemical perspective it is important to note that DFT-MD is capable to study chemical reactions in solution, where force-field MD would fail completely as it cannot account for the change in chemical bonding.

Here, we report on a DFT-based MD simulation of the solvation of ethanol and ethylene. First we describe the simulation methods. Then we show results of geometries and energetics of relevant gas-phase complexes, that will serve as a validation of the numerical methods employed. Subsequently, results of structure, dynamics, and polarization of the solvated species will be shown. We conclude with a discussion.

## II Computational Methods

Electronic structure calculations are performed using the Kohn-Sham KoSh65  formulation of DFT HoKo64 . We employed the gradient-corrected BLYP functional Beck88\_2 ; LeYa88 , that has proven to give a good description of the structure and dynamics of liquid water Sprik . The DFT-based MD simulations of aqueous ethanol and ethylene are performed using the Car-Parrinello method as implemented in the CPMD package CPMD30f . Norm-conserving Martins-Trouillier pseudopotentials TrMa91  are used to restrict the number of electronic states to those of the valence electrons. Cut-off radii for H, O and C atoms were chosen to be 0.50, 1.11 and 1.23 a.u., respectively both for the l=s and l=p terms. The Kohn-Sham orbitals are expanded in a plane wave basis, matching the periodicity of the periodic box with waves up to a kinetic energy of 70 Ry. With this basis set energies and geometries are converged within 0.25 kcal/mol and 0.01 Å , respectively. Vibrational frequencies are converged within 1 %, except for C-O and O-H stretch modes that are underestimated by 3 % and 5 % compared to the basis set limit values ErpMe01\_erratum . In the molecular dynamics calculations, the fictitious mass associated with the plane-wave coefficients is chosen at 900 a.u., which allowed for a time step in the numerical integration of the equations-of-motion of 0.145 fs.

To validate the computational approach we compared CPMD results for energies and structures of relevant gas-phase molecules and complexes with state-of-the-art atomic-orbital DFT-BLYP calculations performed with ADF ADF2002 ; Velde ; Fonseca ; ADF\_basisVII , and other high-level quantum chemical results taken from literature. The gas-phase calculations with CPMD were performed in a large periodic box of 10 Å using the screening technique of BaLa93  to eliminate the interactions among periodic images. We have not included zero-point energies in the energies of the gas-phase compounds. This also holds for computed energies taken from literature and referred to in the present paper. This, to ensure a proper comparison between our results and those from literature.

The MD simulations of the solutions were performed for a ’small’ and a ’large’ system to asses the finite size effects. For the small system a simulation of 10 ps was performed in a cubic periodic box of length 10.07 Å, both for an ethanol solution of 31 waters and 1 ethanol as for an ethylene solution of 32 water and 1 ethylene. For the large system 5 ps simulations were performed using a periodic box with bcc symmetry and a volume of 1977.6 Å3. This periodic cell, a truncated octahedron, is in shape closer to a sphere than a simple cube, and therefore better suited for liquid simulation. The large ethanol solution consisted of 63 waters and 1 ethanol, while the large ethylene solution consisted of 64 waters and 1 ethylene. The box sizes for both the ethanol and ethylene solutions were set to match the experimental densities of the ethanol solutions under ambient conditions. For the ethylene system this will be slightly larger than the experimental density, as the effective volume of ethanol is a bit smaller than the combined volume of ethylene + water. However, we do not expect this to give rise to observable changes in the calculated properties. For reference we performed 10 ps MD simulations of a single ethanol and ethylene in a periodic cubic box of 10 Å and a pure water systems of 32 water molecules in a cubic box of 9.85 Å for simulation times of 10 ps. For all simulations there was an initial equilibration trajectory of 1 ps. Temperature was controlled by a Nosé-Hoover thermostat Nose84\_1 ; Nose84\_2 ; Hoov85  and fixed at 300 K.

## III Gas-phase complexes

The ethanol monomer has two stable conformers very close in energy, the symmetric trans structure and the a-symmetric gauche structure (see fig. 1). The main distinction is the orientation of the OH bond with respect to the CCO plane. A microwave study Kakar  has shown that the trans form is slightly (0.12 kcal/mol) more stable than the gauche form. Table 1 lists the most important geometric data of the trans and gauche conformers, and compares the CPMD results with ADF calculations, B3LYP calculations Borowski , and experimental data Sasada . CPMD and ADF bond lengths differ at most 0.01 Å and angles are within 0.5o. Comparing CPMD with B3LYP and experimental values yields differences upto 0.03 Å and 0.5o, and 0.04 Å and 1o respectively. The calculated energy difference between the two conformers is listed in fig. 2. CPMD, ADF, B3LYP of Ref. Coussan , and MP4 Senent  predict the trans conformer to be stable by 0.07-0.10 kcal/mol, in good agreement with the experimental value of $`0.12`$ kcal/mol. In contrast, the B3LYP calculation of Borowski  yields the opposite, with the gauche conformer stable by $`0.23`$ kcal/mol. Vibrational frequencies obtained by a Fourier transform of the velocity auto correlation function (VACF) of a single ethanol at 300 K yields an OH stretch frequency of 3200 cm-1 and a CH stretch frequencies in the range of 2600-2800 cm-1. This should be compared to the experimental values of 3653 cm-1 Coussan , and 2800-3000 cm-1 Coussan ; Nishi , respectively. The tendency of BLYP to underestimate frequencies by $``$ 10% is a known feature, and also observed in BLYP calculations of water Sprik  and methanol ErpMe01\_erratum .

For the ethanol-water complex we distinguished four complexes, with both the ethanol trans and gauche conformer acting as proton donor or acceptor in the hydrogen bond with water. Table 2 lists the most important geometric data obtained with CPMD and ADF. Fig. 2 shows the energy differences for all complexes. As for the ethanol monomer there is excellent agreement between CPMD and ADF results, with deviations within 0.02 Å and 2o for bond lengths and angles, and the binding energies within 0.2 kcal/mol. This indicates a state-of-the art accuracy of the electronic-structure method employed in CPMD. The CPMD-BLYP calculation predicts the complex with the water molecule donating a proton to the ethanol gauche conformer to be the most stable, with a binding energy relative to the isolated water and trans ethanol conformer of 4.75 kcal/mol. The complex with the ethanol trans conformer donating a proton is less stable by 0.82 kcal/mol. Switching to the other ethanol conformer within either complex destabilizes the complex by 0.07 and 0.4 kcal/mol, respectively. The relative stability of the ethanol donor and acceptor complexes is similar to that for the methanol-water complex where CPMD-BLYP ErpMe01\_erratum  and the complete-basis-set MP2 estimate Kirschner  yields the methanol acceptor complex stable by 0.74 and 0.35 kcal/mol, respectively. There are no experimental data for the structure and energetics of the ethanol-water dimer. The MP2 result of Ref. Masella  provides the only high-level quantum mechanical study reported in literature. This study does not distinguish between the trans and gauche conformers. Therefore, a comparison with our results is somewhat limited. For the ethanol acceptor configuration the MP2 and CPMD-BLYP result for the complexation energy are similar. However, for the ethanol donor configuration the MP2 complexation energy is more than 1.3 kcal/mol larger than the CPMD-BLYP results. Consequently, the MP2 calculations yield an opposite relative stability of the two water-ethanol configurations, with an energy difference of 0.92 kcal/mol. It should be noted that in the same study the methanol-water complex with the methanol as hydrogen bond donor is found to be the most stable, in contradiction with the MP2 basis-set limit result of Ref. Kirschner . This suggests that in the Ref. Masella  a limited basis set or other factors could have lead to spurious reversal of the relative stability of the two water-ethanol configurations. The CPMD-BLYP hydrogen bond interaction energy for the ethanol-water dimer is of the same order of magnitude as the CPMD-BLYP result for the water-water Sprik , water-methanol ErpMe01\_erratum , and methanol-methanol Handgraaf  dimer. Comparison of these three dimers against high-level quantum chemical calculations and experimental values indicates that BLYP underestimates the binding energy by approximately 1 kcal/mol. For the ethanol-water dimer we could expect a similar difference. Here we should add that going from water, via methanol, to ethanol the dispersion forces become increasingly important. These are not accounted for in gradient-corrected functionals such as BLYP. Correlated methods such MP2 incorporate, to a good approximation, dispersion forces. In Ref. Handgraaf  we estimated, in a comparison of BLYP and high-level MP2 calculations Mooij  for the methanol-dimer in non-hydrogen bonded configurations, that the absence of dispersion interaction in BLYP amounted to an underestimate of the binding energy of $``$1 kcal/mol. For the ethanol-water dimer this number could serve as an underestimate. Yet, although by far no insignificant, the magnitude of the deviation is much smaller than the hydrogen-bond interaction. Therefore, it can be argued that for a study of aqueous ethanol neglecting the dispersion interaction is acceptable.

Next we discuss the ethylene-water complexes. We will consider complexes with a single water and a water dimer. The dominant interaction of water with ethylene is a $`\pi `$-H bond, where a proton of the water molecule binds to $`\pi `$ electrons of the double C=C bond. There exist two stable complexes between ethylene and a single water molecule, indicated as the EW1a and EW1b geometry (see fig. 3). Both geometries have Cs symmetry with the plane of the water molecule orthogonal to that of ethylene molecule. In the EW1a structure the water-plane is parallel to the C=C bond, whereas in EW1b it is orthogonal. In table 3 we compare the CPMD-BLYP geometries with ADF-BLYP and MP2Tarakeshwar ; Dupre  and CCSD Tarakeshwar3  results. Fig. 2 shows the energy differences, and compares these to values reported in literature. Again the CPMD and ADF results are in excellent agreement, with energies smaller than 0.02 kcal/mol, and geometries within 0.01 Å and 0.3, except for the non bonded OHnb distance that is related to water orientation, a coordinate along with the energy surface is relatively flat. The BLYP result shows a significant $`\pi `$-H binding energy of $``$ 1.4 kcal/mol. However, in MP2 Tarakeshwar ; Tarakeshwar3 ; Dupre  and CCSD Tarakeshwar3  calculations, binding energies are $``$ 1-2 kcal/mol stronger, accompanied by a shorter $`\pi `$-H bond with differences of 0.12 Å and 0.14 Å for EW1a and EW1b, respectively. The comparison with the experimental value gives a similar picture with the BLYP underestimating the binding energy $``$ 2 kcal/mol.

Matrix isolation studiesEngdahl ; Engdahl86  have revealed that the complex of ethylene with two water molecules consists of a water dimer that has one water molecule $`\pi `$-H bonded to the C=C double bond. Recent high-level MP2 calculations have indicated that the presence of the second water molecule enhances the strength of the $`\pi `$-H bondTarakeshwar2 ; Dupre . From fig. 2 we see that BLYP result are qualitatively in line with this observation, with a binding energy of the (water-dimer)-ethylene complex of 2.25 kcal/mol, up from 1.41 kcal/mol for the single-water binding. However, quantitatively MP2 and BLYP compare less well, with BLYP underestimating both the total binding energy as well as the increase from the single-water binding. Again, the experimental valueEngdahl ; Engdahl86 , that is closer to the MP2 result, shows a significantly stronger binding than BLYP. Both for the water-ethylene and the (water-dimer)-ethylene complex the differences of BLYP with MP2 and experiment are significant relative to the total binding energy. This absence of dispersion interaction in BLYP will be an important factor contribution to this discrepancy.

Overall, we conclude that the reference calculations provide confidence that BLYP is capable of a quantitative study of a the aqueous solvation of a single ethanol molecule where interactions are dominated by relatively strong hydrogen bonds. BLYP qualitatively accounts for the weaker $`\pi `$-H binding in water-ethylene systems. However, comparison with MP2 and experimental data suggest that its strength is significantly underestimated. Yet, for a single solvated ethylene the hydrogen bonds among the solvating water molecules will be the dominant interaction. Hence, we believe that BLYP will be able to quantitative described dilute aqueous ethylene solution. Both for ethanol and ethylene the absence of dispersion attraction in BLYP will have some impact, in particular for the coordination around the CH2 and CH3 groups, where BLYP only accounts for the steric repulsion.

## IV Solvation Structure

### IV.1 Ethanol Solvation

Fig. 4 shows ethanol-water radial distribution functions (RDFs) of the small and large ethanol solutions. The pronounced structure in the hydrogen bonding RDFs (HH, OH, HO, OO) are a clear indication of the presence of hydrogen bonds. All RDFs show that the small system gives a good description of the first solvation shell, while the large system also includes the second solvation shell. The peak positions of the force-field simulations of Fidler and Rodger Fidler , indicated by crosses, are close to our results. Note that the good agreement for the CO RDFs, that are potentially sensitive to a proper description of the dispersion attraction, suggests that the absence of the dispersion interaction in our CPMD-BLYP simulation is of limited importance for a proper description of the aqueous solvation of the hydrophobic group of ethanol. The position of the first peak of the OH- and HO-RDF indicates that the average hydrogen-bond length is 1.7 Å for both the ethanol-donor and -acceptor bond. Integration of the OO-RDF upto the first minimum $`r=3.3`$ Å yields on average 3 water molecules in the first solvation shell of the hydroxyl group, in good agreement with the neutron diffraction value of $``$ 3 Sidhu  . Integrating the methyl-oxygen (C2O) RDF up to $`r=5.7`$ Å indicates that the first solvation shell of the methyl group consists of approximately 21 water molecules. These coordination numbers are of the same order as the experimental estimation of Petrillo et al. Petrillo , who found that there are $`18\pm 2`$ water molecules within 4 Å from the center-of-mass of an ethanol in aqueous solution.

### IV.2 Ethylene solvation

In fig. 5 we show the ethylene-water RDFs of the small and large ethylene solution. Both the carbon-oxygen and the intermolecular carbon-hydrogen RDFs are shown. Experimental or simulation data for ethylene in aqueous solution at room temperature were not available. The comparison between the small and large system shows a larger deviation in the first solvation shell than was found for the ethanol solutions. Apparently, the small ethylene system is not able to accommodate properly the first solvation shell. Integration the CO-RDF of the large system up to the first minimum at 5.7 Å yields a hydration shell of 23 water molecules. More spatial information can be extracted from fig. 6 where, for the large system, the distribution of water H-atoms around the ethylene C=C axis is shown. $`X`$ and $`Y`$ are the components of the vector joining a water hydrogen and the midpoint of the C=C bond. Here, $`X`$ is the distance orthogonal to the C=C axis and $`Y`$ the parallel distance. Note that we used the four-fold symmetry to improve the statistical accuracy yielding four identical quadrants. The ethylene hydrogens are shown to visualize the size of the ethylene molecule, but do not indicate any angle dependence around the C=C axis. The figure shows a well defined elliptic region with no hydrogens present except for two weak, but clearly visible, peaks in the mid plane $`Y=0`$ at $`X=\pm 2.5`$. These should be attributed to the presence of $`\pi H`$ bonded configurations.

### IV.3 Water Structure

In fig. 7 the calculated water oxygen-oxygen RDFs of both the ethanol and the ethylene solutions are compared to those calculated for pure water. Note that the first solvation shell contains a large fraction of the total number of water molecules, especially for the small system. Structural changes due to the solute molecule should therefore be detectable by comparing these systems with the pure water system. For the ethanol solution we see a small drop of the first peak accompanied by a slight increase at the first minimum. This indicates some decrease in the hydrogen-bond structure, when compared to the pure water. As the small ethylene-water solution is too small to accommodate a fully relaxed water solvation shell fig. 7 only shows the RDFs of the large ethylene system. The RDFs show a similar behavior as for the ethanol solvation, with a slight decrease of the first peak and a small increase of the RDF in the region around the first minimum. The small changes in the structure of the solvating water shell around ethanol and ethylene indicate that the hydrogen-bonded water network is very flexible and can easily open up to accommodate small apolar solute groups without significantly changing its structure. Our findings are consistent with the neutron diffraction experiments of Turner and Soper Turner  and force-field molecular dynamics simulation of Fidler and Rodger Fidler  who also did not find any evidence of structural enhancement in the hydration shell of ethanol. Turner and Soper found that only for larger hydrophobic solutes this effect was experimentally detectable, but even then very small. The same trend was found for in a CPMD-BLYP molecular dynamics study of the solvation of methanol ErpMe01\_erratum .

### IV.4 Hydrogen Bonds

To examine the hydrogen bond statistics in the ethanol solution we adopted the definition of Ref. FeHa90 : two molecules are hydrogen bonded if simultaneously the inter-oxygen distance is less than 3.4 Å and the OHO angle is smaller than 30. With this definition we found that, for the large system, ethanol oxygen donates on average 0.9 hydrogen bonds and accepts 1.7. For the small system we found 1.0 and 1.6, respectively. This in consistent with the fact that approximately three water molecules occupy the first solvation shell of the hydroxyl group. For comparison, in a CPMD-BLYP simulation we found that methanol in a dilute aqueous solution donates on average 0.9 hydrogen bonds and accepts 1.5 ErpMe01\_erratum . ¿From the pure water simulation these numbers were measured to be 1.7.

For the ethylene molecule a well defined definition for the $`\pi `$-H hydrogen bond does not exist. To investigate the influence of the $`\pi `$-H hydrogen bond on the solvation structure we looked at the water hydrogen positions relative to the C=C axis (see fig 6). Inside the elliptic region around the C=C bond that is depleted of hydrogens we clearly detect near the center-of-mass of the ethylene molecule at a distance of $``$ 2.5 Å a region with an increased amount of water hydrogens. Integration over this region with $`0<X<2.6`$ and $`0<Y<1`$ yields a value of 0.42. This implies that approximately 40% of the time a water molecule is oriented towards the double bond, forming a $`\pi `$-H bond.

## V Solute Dynamics

The time scale of the present simulations (5-10 ps) allows for an analysis of the short-time dynamics. As we mentioned above BLYP underestimates most vibrational frequencies by more than 10 %. However, the frequency shifts upon solvation still provide valuable information and allow for a direct comparison to experiments. Fig. 8 shows the spectrum of the velocity auto correlation functions (VACFs) of the ethanol hydroxyl hydrogen for the large ethanol solution and for the single isolated ethanol. The spectrum is of limited accuracy due to the relative short trajectory. However, the calculations show that, upon solvation, there is a clear red shift of about 200 cm-1 of the OH stretch frequency. Experimental data for the OH frequency shift of ethanol in dilute aqueous solution is not available. However, the OH red shift is a typical characteristic of a hydrogen bonded liquid. A comparison with measured shifts in liquid ethanol, from $`3676cm^1`$Perchard  to 3330 cm-1Lock , shows a similar trend. In a CPMD-BLYP simulation of a dilute aqueous solution of methanol we found a similar red shift for the methanol OH frequency of about 200 cm-1

The vibrational spectrum of the VACF of the ethylene molecule is not shown. As for the ethanol system there was a limited accuracy due to the relatively short calculated trajectory. We observed that upon solvation the CH and CC peaks do broaden. However an estimate for peak shift falls outside the accuracy of the calculated spectra. The matrix isolation studies of Engdahl and Nelander Engdahl ; Engdahl86  show minor changes in frequencies when isolated ethylene is compared with the ethylene-water complex, with the largest shift being a blue shift of 12 cm-1 of the out-of-plane bending mode ($`\nu _7947`$ cm-1).

## VI Polarization

As the electronic structure is an intrinsic part of a CPMD simulation, detailed information about the electronic charge distribution can be obtained. To quantify the charge distribution we used the method of maximally localized Wannier functions that transforms the Kohn-Sham orbitals into Wannier functions, whose centers (WFC) can be assigned a chemical meaning such as being associated with an electron bonding- or lone-pair Marzari ; Silvestrelli4 . It also provides a unambiguous route to determine the dipole moment of individual molecules in a condensed phase by assuming the electronic charge to be distributed as point charges located on the WFCs.

Table 4 list the calculated dipole moments for the gas-phase systems and the average dipole moment for the solvated systems. The latter were obtained from 28 independent configurations of the small solvated systems. The distributions of these dipole moments are shown in fig. 9. For ethanol we observe an significant increase in going from the gas-phase molecule via the water-ethanol cluster to the fully solvated system. This trend is very similar to what is found is found for water Silvestrelli2 ; Silvestrelli3  and methanol Handgraaf  and may be considered typical feature of a strongly hydrogen bonding molecule. The calculated distribution of dipole moments of solvated ethanol shows that thermally driven fluctuations give rise to a significant variation ranging from 2.0 D to 4.0 D.

Also for ethylene we observed a significant change of the dipole moment upon solvation. Being apolar in the gas phase, the average dipole moment increases up to 0.5 D when complexed with the water dimer or when in solution. In solution it exhibits rare fluctuations where the dipole moment reaches values of up to 1.0 D. Fig. 10 shows a typical snapshot where such an extreme high value of the dipole moment is reached. The two WFCs of the ethylene double bond are located near the middle of the C=C bond, just below and above the plane of the molecule. In the figure we see clearly that the upper $`\pi `$-WFC is acting as proton acceptor. The $`\pi `$-H bond shifts this WFC further out of plane inducing an ethylene dipole moment orthogonal to the plane of the molecule. The snapshot suggests an enhancement of the ethylene polarization by the fact that the $`\pi `$-H bonded water molecule has in turn a rare four-fold hydrogen-bonded water coordination, with three of the waters donating a proton. This should be seen as a manifestation of the strengthening of the $`\pi `$-H bond upon complexing water molecules with the water molecule that is $`\pi `$-H bonded to the ethylene, a feature found in gas-phase MP2 calculations Tarakeshwar2 ; Dupre  and also in the present BLYP calculations, and discussed above in section III.

## VII Conclusions

We have studied the solvation of ethanol and ethylene in water by DFT-based (Car-Parrinello) Molecular Dynamics. Validation of the computational approach by comparing structure and energetics of relevant gas-phase complexes against experimental results and state of the art quantum chemical calculations showed that CPMD employing the DFT-BLYP is capable of qualitatively describing the aqueous solvation of a single ethanol of ethylene molecule.

The structural properties of the ethanol solvation were in good agreement with both neutron diffraction data Turner  and force-field MD simulations Fidler . We found that in aqueous solution ethanol accepts on average 1.7 hydrogen bonds and donates on average 0.9 hydrogen bonds. For ethylene we found it has approximately 0.4 $`\pi `$-H bonds with a water molecule. Both for ethanol and ethylene the simulations provided no structural evidence for hydrophobic hydration: the structure of the hydrating water shell was not enhanced compared to that of pure water. The calculation even indicated a slight decrease in the structure. For aqueous ethanol the calculated red shift of the hydroxyl vibration upon solvation was consistent with experimental findings.

Analysis of the electronic charge distribution by means of Wannier functions showed that the interaction with water can significantly increase the dipole moment of the ethanol and the ethylene molecule. The average dipole moment of ethanol increases from 1.8 D in the gas phase to 3.1 D in aqueous solution. Ethylene, that is apolar in the gas phase, attains an average dipole moment of 0.5 D in solution. We identified configurations, with a $`\pi `$-H bonded water molecule that has a rare four-fold hydrogen-bonded water coordination, where the instantaneous dipole moment of ethylene takes values of up to 1 D. Such configurations with large solute dipole moments may also play an important role in activating chemical reactions involving the solutes as we have seen in a CPMD-BLYP study of the acid-catalyzed hydration of ethene ownfut . The strong polarization effect raises questions towards the common consideration of thinking ethylene as a apolar and hydrophobic molecule. The electronic charge analysis also points out the necessity of polarizable force fields for both ethanol and ethylene when dissolved in water. The present results may be considered valuable for the design of such force-fields.

###### Acknowledgements.

T.S.v.E acknowledges NWO-CW (Nederlandse Organisatie voor Wetenschappelijk Onderzoek, Chemische Wetenschappen). E.J.M. acknowledges the Royal Netherlands Academy of Art and Sciences for financial support. We acknowledge support from the Stichting Nationale Computerfacileiten (NCF) and the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) for the use of supercomputer facilities.

TABLES

FIGURES