text:
# References

Preprint SPhT-T04/010
$`1/N^2`$ correction to free energy in hermitian two-matrix model
B.Eynard1, A.Kokotov2, D.Korotkin2
1 Service de Physique Théorique, CEA/Saclay, Orme des Merisier F-91191 Gif-sur-Yvette Cedex, France

2 Concordia University, 7141 Sherbrooke West, Montreal H4B1R6, Montreal, Quebec, Canada

Abstract. Using the loop equations we find an explicit expression for genus 1 correction in hermitian two-matrix model in terms of holomorphic objects associated to spectral curve arising in large N limit. Our result generalises known expression for $`F^1`$ in hermitian one-matrix model. We discuss the relationship between $`F^1`$, Bergmann tau-function on Hurwitz spaces, G-function of Frobenius manifolds and determinant of Laplacian over spectral curve.

In this letter we derive an explicit formula for the $`1/N^2`$ correction to free energy $`F`$ of hermitian two-matrix model:

$$e^{N^2F}:=𝑑M_1𝑑M_2e^{N\mathrm{tr}\{V_1(M_1)+V_2(M_2)M_1M_2\}}.$$
(1)
It is hard to overestimate the interest to random matrix models in modern physics and mathematics; we just mention their appearance in statistical physics, condensed matter and 2d quantum gravity (see e.g. ) . The expansion $`F=_{G=0}^{\mathrm{}}N^{2G}F^G`$ ($`N`$ is the matrix size) in hermitian matrix models is one of the cornerstones of the theory, due to its clear physical interpretation as topological expansion of the functional integral, which appears in $`N\mathrm{}`$ limit; in statistical physics the term $`F^G`$ plays the role of free energy for statistical physics model on genus $`G`$ Riemann surface. From the whole zoo if the random matrices one of the simplest is the hermitian one-matrix model with partition function $`e^{N^2F}=𝑑Me^{N\mathrm{tr}V(M)}`$ (V is a polynomial), which can be used as testing ground for the methods applied in more general situations of two- and multi- matrix models. The most rigorous way to compute the $`1/N^2`$ expansion for both one-matrix and two-matrix models is based on the loop equations. The loop equations follow from the reparametrization invariance of matrix integrals; for one-matrix case the loop equations were used to compute $`F^1`$ (see ). Later the loop equations were written down for the case of two-matrix model  and $`F^1`$ was found for the case when the spectral curve has genus zero and one ; for arbitrary genus of spectral curve of two-matrix model only the leading term $`F^0`$ is known (see ).

Let us write down polynomials $`V_1`$ and $`V_2`$ in the form $`V_1(x)=_{k=1}^{d_1+1}\frac{u_k}{k}x^k`$ and $`V_2(y)=_{k=1}^{d_2+1}\frac{v_k}{k}y^k`$. It is sometimes convenient to think of $`V_1`$ and $`V_2`$ as infinite formal power expansions: $`V_1(x)=_{k=1}^{\mathrm{}}\frac{u_k}{k}x^k`$, $`V_2(y)=_{k=1}^{\mathrm{}}\frac{v_k}{k}y^k`$, where coefficients $`u_k`$ vanish for $`kd_1+2`$, and $`v_k`$ vanish for $`kd_2+2`$. According to this point of view the operators of differentiation with respect to coefficients of $`V_1`$ and $`V_2`$ have the following meaning (see ):

$$\frac{\delta }{\delta V_1(x)}|_x:=\left\{\underset{k=1}{\overset{\mathrm{}}{}}x^{k1}k_{u_k}\right\}|_{u_k=0,kd_1+2},\frac{\delta }{\delta V_2(y)}|_y:=\left\{\underset{k=1}{\overset{\mathrm{}}{}}y^{k1}k_{v_k}\right\}|_{v_k=0,kd_1+2}.$$
(2)
As it was discussed in detail in , (2) is a formal notation which makes sense only order by order in the infinite power series expansion; it allows to write an infinite number of equations at once. Consider the resolvents $`𝒲(x)=\frac{1}{N}\mathrm{tr}\frac{1}{xM_1}`$ and $`\stackrel{~}{𝒲}(y)=\frac{1}{N}\mathrm{tr}\frac{1}{yM_2}`$ The free energy of two-matrix model (1) satisfies the following equations with respect to coefficients of polynomials $`V_1`$ and $`V_2`$:

$$\frac{\delta F}{\delta V_1(x)}=𝒲(x),\frac{\delta F}{\delta V_2(y)}=\stackrel{~}{𝒲}(y).$$
(3)
The equations (3) were solved in  in the zeroth order assuming the finite-gap structure of distribution of eigenvalues of $`M_1`$ (and, a posteriori, also of $`M_2`$) as $`N\mathrm{}`$. Here we find the next coefficient $`F^1`$, using the loop equations. The spectral curve $``$ is defined by the following equation, which arises in the zeroth order approximation:

$$^0(x,y):=(V_{1}^{}{}_{}{}^{}(x)y)(V_{2}^{}{}_{}{}^{}(y)x)𝒫^0(x,y)+1=0$$
(4)
where the polynomial of two variables $`𝒫^0(x,y)`$ is the zeroth order term in $`1/N^2`$ expansion of the polynomial

$$𝒫(x,y):=\frac{1}{N}\mathrm{tr}\frac{V_{1}^{}{}_{}{}^{}(x)V_{1}^{}{}_{}{}^{}(M_1)}{xM_1}\frac{V_{2}^{}{}_{}{}^{}(y)V_{2}^{}{}_{}{}^{}(M_2)}{yM_2};$$
(5)
the point $`P`$ of the curve is the pair of complex numbers $`(x,y)`$ satisfying (4) (on the “physical” sheet the equation of spectral curve (4) defines an implicit function $`y(x)`$, which gives the zeroth order approximation to $`V_{2}^{}{}_{}{}^{}(x)𝒲(x)`$). The spectral curve (4) comes together with two meromorphic functions $`f(P)=x`$ and $`g(P)=y`$, which project it down to $`x`$ and $`y`$-planes, respectively. These functions have poles only at two points of $``$, called $`\mathrm{}_f`$ and $`\mathrm{}_g`$: at $`\mathrm{}_f`$ function $`f(P)`$ has simple pole, and function $`g(P)`$ \- pole of order $`d_1`$ with singular part equal to $`V_{1}^{}{}_{}{}^{}(f(P))`$. At $`\mathrm{}_g`$ the function $`g(P)`$ has simple pole, and function $`f(P)`$ \- pole of order $`d_2`$ with singular part equal to $`V_{2}^{}{}_{}{}^{}(g(P))`$. In addition, in our normalization of partition function (1) we have the asymptotics $`𝒲(x)_x\mathrm{}1/x+\mathrm{}`$ and $`\stackrel{~}{𝒲}(y)_y\mathrm{}1/y+\mathrm{}`$, which imply  $`\mathrm{Res}_\mathrm{}_fgdf=1`$ and $`\mathrm{Res}_\mathrm{}_gfdg=1`$, respectively. Therefore, one gets the moduli space $``$ of triples $`(,f,g)`$, where functions $`f`$ and $`g`$ have this pole structure. The natural coordinates on this moduli space are coefficients of polynomials $`V_1`$ and $`V_2`$ and $`g`$ numbers, called “filling fractions” $`ϵ_\alpha =\frac{1}{2\pi i}_{a_\alpha }g𝑑f`$, where $`(a_\alpha ,b_\alpha )`$ is some basis of canonical cycles on $``$. The additional constraints which should be imposed a posteriori to make the “filling fractions” dependent on coefficients of polynomials $`V_1`$ and $`V_2`$ are (according to one-matrix model experience, these conditions correspond to non-tunneling between different intervals of eigenvalues support): $`_{b_a}g𝑑f=0`$. Denote the zeros of differential $`df`$ by $`P_1,\mathrm{},P_{m_1}`$ ($`m_1=d_2+2g+1`$) (these points play the role of ramification points if we realize $``$ as branched covering by function $`f(P)`$); their projections on $`x`$-plane are the branch points, which we denote we denote by $`\lambda _j:=f(P_j)`$ . The zeros of the differential $`dg`$ (the ramification points if we consider $``$ as covering defined by function $`g(P)`$) we denote by $`Q_1,\mathrm{},Q_{m_2}`$ ($`m_2=d_1+2g+1`$); there projections on $`y`$-plane (the branch points) we denote by $`\mu _j:=g(Q_j)`$. We shall assume hat our potentials $`V_1`$ and $`V_2`$ are generic i.e. all zeros of differentials $`df`$ and $`dg`$ are simple, and none of the zeros of $`df`$ coincides with a zero of $`dg`$. If is well-known (see for instance ) how to express all standard algebro-geometrical objects on $``$ in terms of the previous data. In particular, the Bergmann bidifferential $`B(P,Q)=d_Pd_Q\mathrm{ln}E(P,Q)`$ ($`E(P,Q)`$ is the prime-form) can be represented as follows:

$$B(P,Q)=\frac{\delta g(P)}{\delta V_1(f(Q))}|_{f(Q)}df(P)df(Q)$$
(6)
The Bergmann bidifferential has the following behaviour near diagonal $`PQ`$: $`B(P,Q)=\{(z(P)z(Q))^2+\frac{1}{6}S_B(P)+o(1)\}dz(P)dz(Q)`$, where $`z(P)`$ is some local coordinate; $`S_B(P)`$ is the Bergmann projective connection . Consider also the four-differential $`D(P,Q)=d_Pd_Q^3\mathrm{ln}E(P,Q)`$, which has on the diagonal the pole of 4th degree: $`D(P,Q)=\{6(z(P)z(Q))^4+O(1)\}dz(P)(dz(Q))^3`$. From $`B(P,Q)`$ and $`D(P,Q)`$ it is easy to construct meromorphic normalized (all $`a`$-periods vanish) $`1`$-forms on $``$ with single pole; in particular, if the pole coincides with ramification point $`P_k`$, the natural local parameter near $`P_k`$ is $`x_k(P)=\sqrt{f(P)\lambda _k}`$; then $`B(P,P_k):=\frac{B(P,Q)}{dx_k(Q)}|_{Q=P_k}`$ and $`D(P,P_k):=\frac{D(P,Q)}{(dx_k(Q))^3}|_{Q=P_k}`$ are meromorphic normalized 1-forms on $``$ with single pole at $`P_k`$ and the following singular parts:

$$B(P,P_k)=\left\{\frac{1}{[x_k(P)]^2}+\frac{1}{6}S_B(P_k)+o(1)\right\}dx_k(P);D(P,P_k)=\left\{\frac{6}{[x_k(P)]^4}+O(1)\right\}dx_k(P)$$
(7)
as $`PP_k`$, where $`S_B(P_k)`$ is the Bergmann projective connection computed at the branch point $`P_k`$ with respect to the local parameter $`x_k(P)`$.

Equations (3) in order $`1/N^2`$ look as follows (we write only equations with respect to $`V_1`$):

$$\frac{\delta F^1}{\delta V_1(f(P))}=Y^1(P)$$
(8)
where the $`Y^1`$ is the (taken with minus sign) $`1/N^2`$ contribution to the resolvent $`𝒲`$. The function $`Y^1`$ can be computed using the loop equations  and the “normalization conditions”

$$_{a_\alpha }Y^1(P)𝑑f(P)=0$$
(9)
over all basic $`a`$-cycles (these conditions mean that the “filling fractions” do not have the $`1/N^2`$ correction). We introduce also the polynomial

$$(x,y):=(V_1(x)y)(V_2(y)x)𝒫(x,y)+1,$$
(10)
the function $`𝒰(x,y)`$, which is a polynomial in $`y`$ and rational function in $`x`$:

$$𝒰(x,y):=\frac{1}{N}\mathrm{tr}\frac{1}{xM_1}\frac{V_{2}^{}{}_{}{}^{}(y)V_{2}^{}{}_{}{}^{}(M_2)}{yM_2},$$
(11)
and rational function $`𝒰(x,y,z)`$:

$$𝒰(x,y,z):=\frac{\delta 𝒰(x,y)}{\delta V_1(z)}=\mathrm{tr}\frac{1}{xM_1}\frac{V_{2}^{}{}_{}{}^{}(y)V_{2}^{}{}_{}{}^{}(M_2)}{yM_2}\mathrm{tr}\frac{1}{zM_1}N^2𝒰(x,y)𝒲(z).$$
(12)
Then the loop equation looks as follows:

$$𝒰(x,y)=xV_{2}^{}{}_{}{}^{}(y)+\frac{(x,y)}{yY(x)}\frac{1}{N^2}\frac{𝒰(x,y,x)}{yY(x)},$$
(13)
where $`Y(x):=V_{1}^{}{}_{}{}^{}(x)𝒲(x)`$; it arises as a corollary of reparametrization invariance of the matrix integral (1) . To use the loop equation effectively we need to consider the $`1/N^2`$ expansion of all of their ingredients. In this way we get the following expression for $`Y^1`$:

$$Y^1(P)df(P)=\frac{𝒫^1(f(P),g(P))df(P)}{_y^0(f(P),g(P))}+\underset{QP:f(Q)=f(P)}{}\frac{B(P,Q)}{df(Q)}\frac{1}{g(P)g(Q)},$$
(14)
where $`_y^0(x,y)`$ means partial derivative with respect to the second argument. All ingredients of (14) arise already in the leading term, except $`𝒫^1`$. However, from (5) we see that $`𝒫(x,y)`$ is a polynomial of degree $`d_11`$ with respect to $`x`$ and $`d_21`$ with respect to $`y`$; moreover, the coefficient in front of $`x^{d_11}y^{d_21}`$ does not have $`1/N^2`$ correction. Thus we can conclude that the one-form $`Y^1(P)df(P)`$ is non-singular on the spectral curve outside of the branch points $`P_m`$ (where it has poles of order 4); moreover, the first term in (14) is non-singular on $``$ (the first order zeros of $`_y^0`$ at the branch points are cancelled by first order zeros of $`df(P)`$ at these points). The form of singular parts at $`P_m`$ allows to determine $`Y^1(P)df(P)`$ completely in terms of differentials $`B(P,P_k)`$ and $`D(P,P_k)`$ if we take into account the absence of $`1/N^2`$ correction to the “filling fractions” (9); the result looks as follows:

$$Y^{(1)}(P)df(P)=\underset{k=1}{\overset{m_1}{}}\left\{\frac{1}{96g^{}(P_k)}D(P,P_k)+\left[\frac{g^{\prime \prime \prime }(P_k)}{96g^2(P_k)}\frac{S_B(P_k)}{24g^{}(P_k)}\right]B(P,P_k)\right\}$$
(15)
Then the solution of (8), (15), which is symmetric with respect to the projection change (and, therefore, satisfies also equations (3) with respect to $`V_2`$), looks as follows:

$$F^1=\frac{1}{24}\mathrm{ln}\left\{\tau _f^{12}(v_{d_2+1})^{1\frac{1}{d_2}}\underset{k=1}{\overset{m_1}{}}dg(P_k)\right\}+\frac{d_2+3}{24}\mathrm{ln}d_2$$
(16)
where $`\tau _f`$ is the so-called Bergmann tau-function on Hurwitz space , which satisfies the following system of equations with respect to the branch points $`\lambda _k`$:

$$\frac{}{\lambda _k}\mathrm{ln}\tau _f=\frac{1}{12}S_B(P_k);$$
(17)
In derivation of (15) we have used the following variational formulas, which can be easily proved in analogy to Rauch variational formulas:

$$\frac{\delta \lambda _k}{\delta V_1(f(P))}g^{}(P_k)df(P)=B(P,P_k),$$
(18)
$$\frac{\delta \{g^{}(P_k)\}}{\delta V_1(f(P))}|_{f(P)}df(P)=\frac{1}{4}\left\{D(P,P_k)\frac{g^{\prime \prime \prime }(P_k)}{g^{}(P_k)}B(P,P_k)\right\}$$
(19)
The Bergmann tau-function (17) appears in many important problems: it coincides with isomonodromic tau-function of Hurwitz Frobenius manifolds , and gives the main contribution to $`G`$-function (solution of Getzler equation) of these Frobenius manifolds; it gives the most non-trivial term in isomonodromic tau-function of Riemann-Hilbert problem with quasi-permutation monodromies. Finally, its modulus square essentially coincides with determinants of Laplace operator in metrics with conical singularities over Riemann surfaces . The solution of the system (17) looks as follows . Define the divisor $`(df)=2\mathrm{}_f(d_2+1)\mathrm{}_g+_{k=1}^{m_2}P_k:=_{k=1}^{m_2+2}r_kD_k`$. Choose some initial point $`P\widehat{}`$ and introduce corresponding vector of Riemann constants $`K^P`$ and Abel map $`𝒜_\alpha (Q)=_P^Qw_\alpha `$ ($`w_\alpha `$ form the basis of normalized holomorphic 1-forms on $``$). Since some points of divisor $`(df)`$ have multiplicity 1, we can always choose the fundamental cell $`\widehat{}`$ of the universal covering of $``$ in such a way that $`𝒜((df))=2K^P`$ (for an arbitrary choice of fundamental domain these two vectors coincide only up to an integer combination of periods of holomorphic differentials), where the Abel map is computed along the path which does not intersect the boundary of $`\widehat{}`$. The main ingredient of the Bergmann tau-function is the following holomorphic multivalued $`(1g)g/2`$-differential $`𝒞(P)`$ (the higher genus analog of Dedekind eta-function) on $``$:

$$𝒞(P):=\frac{1}{W(P)}\underset{\alpha _1,\mathrm{},\alpha _g=1}{\overset{g}{}}\frac{^g\mathrm{\Theta }(K^P)}{z_{\alpha _1}\mathrm{}z_{\alpha _g}}w_{\alpha _1}(P)\mathrm{}w_{\alpha _g}(P).$$
(20)
where $`W(P):=\mathrm{det}_{1\alpha ,\beta g}w_\beta ^{(\alpha 1)}(P)`$ denotes the Wronskian determinant of holomorphic differentials at point $`P`$; $`K^P`$ is the vector of Riemann constants with basepoint $`P`$. Introduce the quantity $`𝒬`$ defined by

$$𝒬=[df(P)]^{\frac{g1}{2}}𝒞(P)\underset{k=1}{\overset{m+N}{}}[E(P,D_k)]^{\frac{(1g)r_k}{2}};$$
(21)
which is independent of the point $`P`$. Then the Bergmann tau-function (17) of Hurwitz space is given by the following expression :

$$\tau _f=𝒬^{2/3}\underset{k,l=1k<l}{\overset{m+n}{}}[E(D_k,D_l)]^{\frac{r_kr_l}{6}};$$
(22)
together with (16) this gives the answer for $`1/N^2`$ correction in two-matrix model. If $`\tau _f`$ and $`\tau _g`$ are Bergmann tau-functions (17) corresponding to divisors $`(df)`$ and $`(dg)`$, respectively, then

$$\left(\frac{\tau _f}{\tau _g}\right)^{12}=C\frac{(u_{d_1+1})^{1\frac{1}{d_1}}}{(v_{d_2+1})^{1\frac{1}{d_2}}}\frac{_kdf(Q_k)}{_kdg(P_k)}$$
(23)
where $`C=d_1^{d_1+3}/d_2^{d_2+3}`$ is a constant independent of moduli parameters. Using the transformation (23) of the Bergmann tau-function under projection change, we find that the solution expression (16) for $`F^1`$ satisfies also the necessary equations with respect to $`V_2`$. This could be considered as a confirmation of consistency of our computation. Derivatives of function $`F^1`$ (16) with respect to the filling fractions look as follows:

$$\frac{F^1}{ϵ_\alpha }=_{b_\alpha }Y^1(P)𝑑f(P);$$
(24)
these equations are $`1/N^2`$ counterparts of Seiberg-Witten type equations for $`F^0`$ (see for example ).

One-matrix model. If potential $`V_2`$ is quadratic, integration with respect to $`M_2`$ in (1) can be performed explicitly, and the free energy (16) gives rise to the free energy of one-matrix model. The spectral curve $``$ in this case becomes hyperelliptic, and the formula (16) gives, using the expression for $`\tau _f`$ obtained in :

$$F^1=\frac{1}{24}\mathrm{ln}\left\{\mathrm{\Delta }^3(\mathrm{det}𝐀)^{12}\underset{k=1}{\overset{2g+2}{}}g^{}(\lambda _k)\right\}$$
(25)
where $`\lambda _k`$, $`k=1,\mathrm{},2g+2`$ are branch points of $``$; $`\mathrm{\Delta }`$ is their Wronskian determinant; $`𝐀`$ is the matrix of $`a`$-periods of non-normalized holomorphic differentials on $``$; this agrees with previous results .

$`F^1`$, isomonodromic tau-function and $`G`$-function of Frobenius manifolds. The genus 1 correction to free energy in topological field theories is given by so-called $`G`$-function, which for an arbitrary $`m`$-dimensional Frobenius manifold related to Hurwitz space looks as follows : $`\mathrm{exp}\{G\}=\tau _f^{1/2}_{k=1}^m\{\mathrm{res}|_{P_k}\frac{\phi ^2}{df}\}^{1/48}`$, where $`\tau _f`$ is the Bergmann tau-function, $`\phi `$ is an “admissible” one-form on underlying Riemann surface. If, trying to build an analogy with our formula (16) for $`F^1`$, we formally choose $`\phi (P)=dg(P)`$, the formula for G-function coincides with (16) up to small details like sign, additive constant, and the highest coefficient of polynomial $`V_2`$. However, the differential $`dg`$ is not admissible, and, therefore, does not really correspond to any Frobenius manifold; therefore, the true origin of his analogy is still unclear

$`F^1`$ and determinant of Laplace operator. Existence of close relationship between $`F^1`$ and determinant of Laplace operator was suggested by several authors (see e.g.  for hermitial one-matrix model,  for hermitian two-matrix model and, finally,  for normal two-matrix model with simply-connected support of eigenvalues). After an appropriate regularization the (formal) determinant of Laplace operator $`\mathrm{\Delta }^f`$ over $``$ in the singular metric $`|df(P)|^2`$ is given by the expression  $`\mathrm{det}\mathrm{\Delta }^f=C𝒜\{\mathrm{det}\mathrm{}𝐁\}|\tau _f|^2,`$ where $`𝒜`$ is a regularized area of $``$, $`𝐁`$ is the matrix of $`b`$-periods of $``$, $`C`$ is a constant. In the “physical” case of real coefficients of $`V_1`$ and $`V_2`$ and real filling fractions the empirical expression for $`\mathrm{ln}\{\mathrm{det}\mathrm{\Delta }^f\}`$ differs from our expression (16) by several explicit terms. Therefore, the relationship between hermitian and normal two-matrix models on the level of $`F^1`$ seems to be not as straightforward as on the level of functions $`F^0`$ .

Acknowledgements We thank M.Bertola, L.Chekhov, B.Dubrovin, T.Grava, V.Kazakov, I.Kostov, M.Staudacher and S.Theisen for important discussions. The work of BE was partially supported by the EC ITH Network HPRN-CT-1999-000161. The work of DK was partially supported by NSERC, NATEQ and Humboldt foundation. AK and DK thank Max-Planck Institute for Mathematics in Bonn for support and nice working conditions. DK thanks also SISSA and CEA for support and hospitality.


format: text

text:
# Multiperiodic RR Lyrae stars in the center of the Galaxy

## 1 Introduction

It is well known that RR Lyrae stars play very important role in modern astrophysics. They are used as distance indicators, their spatial distribution can provide us with useful information about the structure of the Galaxy and they help us understand stellar evolution and pulsation. Especially pulsation theory benefits from investigation of large number of RR Lyrae stars after the discovery of nonradial oscillations in these stars. Bulge catalog presented in this paper contains over 700 RR Lyrae stars showing multiperiodic behavior. In many cases such behavior can only be explained by excitation of nonradial pulsation modes.

## 2 Observational data

In the years 1997 - 1999, during the second phase of the Optical Gravitational Lensing Experiment - OGLE-II, 49 fields in the center of the Galaxy were monitored. The fields covered approximately 11 squared degrees with most of them being located in regions of high stellar density. With the aid of Warsaw 1.3-m telescope at Las Campanas Observatory of the Carnegie Institution of Washington about 200 observations in the *I* band were made for each field, mostly in driftscan mode. Technical details of the telescope are given in Udalski, Kubiak & Szymański (1997). The photometry was obtained using Difference Image Analysis software written by Woźniak (2000).

## 3 Data analysis

Image subtraction method preclassified as variable stars about 4500 objects per field. Since there are 49 fields in OGLE-II database, selection of RR Lyrae stars had to be automated. There are approximately 200 measurements from three observational seasons for each star and the time span is about 1000 days. Observations are sampled unevenly. As the first step Lomb-Scargle periodograms, with frequency resolution (5$`\times `$$`10^4`$ c/d) for all stars were computed. Such a resolution is high enough to classify a star as a periodic variable and to examine if it is a RR Lyr star. At later stages periodograms for all selected RR Lyr stars were computed with much higher resolution. After a periodogram was calculated for every star, the most probable period was selected out of 25 periods corresponding to the highest peaks in its periodogram. For all of those 25 periods my algorithm fitted the data with a Fourier series with various number of harmonic terms. This number ran from 1 to 6 for each period which was done because of the diversity of RR Lyrae light curve shapes. Strongly skewed curves of RRab stars require more harmonic terms than sinusoidal curves of RRc stars. The criterion of selection was the smallest value of $`\chi ^2`$ per degree of freedom. The program always chose the pair (period, number of harmonics) that minimized that quantity. As expected in the case of RRc stars selected number of harmonic terms was mainly 1 or 2, but for RRab stars it was typically 4 or 5 in most cases.

### 3.1 Identification of RR Lyrae stars

When the most probable period and the best number of harmonic terms were found for every star the algorithm shifted the zero point of every light curve in the time domain to the maximum of the curve. The maximum was estimated from the Fourier fit. To separate RR Lyrae stars from all other variables the code computed the following parameters:

$$A_{ij}=\frac{A_i}{A_j}$$
(1)
$$\mathrm{\Phi }_{ij}=j\varphi _ii\varphi _j$$
(2)
Quantities $`A_i`$ and $`\varphi _i`$ are respectively the amplitude and the phase of a cosine term with $`\omega =i\omega _0`$ where $`\omega _0`$ is the frequency of variability. Parameters defined above, when combined with the value of the period, allow not only to distinguish RR Lyrae stars from other variables, but also to separate RRab and RRc classes. In the search for RRabs my code used parameters $`A_{21}`$, $`A_{32}`$, $`A_{43}`$ and $`\mathrm{\Phi }_{21}`$. For RRcs only $`A_{21}`$, $`A_{32}`$ and $`\mathrm{\Phi }_{21}`$ were used. Output RRab and RRc samples were contaminated by small amount of variables of different types, most frequently High Amplitude $`\delta `$ Sct stars. These were eliminated with the aid of color-magnitude diagrams and by inspection of all the light curves by eye. This method of selecting RR Lyrae stars was also used by Mizerski and Bejger (2001).

### 3.2 Detection of additional periodicities

After RR Lyrae stars were separated from all other variables by the algorithm described in subsection 3.1, the search for multiperiodic behavior was performed. The whole process was fully automated and for each RR Lyrae star consisted of three main phases:

* Calculation of Lomb-Scargle periodogram with high resolution in frequency
(5$`\times `$$`10^6`$ c/d) and determination of the most probable period of variability.
* Prewhitening of the light curve with the determined period with dynamically adjustable number of harmonic terms.
* Calculation of Lomb-Scargle periodogram of the prewhitened light curve.

In the first phase my code calculated Lomb-Scargle periodogram of the original light curve. This method of periodogram calculation works well for unevenly sampled data. Setting frequency resolution to 5$`\times `$$`10^6`$ c/d allowed me to determine the first, main period with high accuracy, considering that the time span of the observations was about 1000 days. Out of the periodogram 25 frequencies corresponding to the highest peaks were selected. For all those frequencies my code fitted the data by means of least squares method with a Fourier series. Number of harmonic terms was adjusted dynamically in a range between 1 and 6. The program searched for a pair (frequency,number of harmonics) which minimized the value of $`\chi ^2`$ of the fit. After such a pair was found the inverse of the selected frequency was taken for the primary period of pulsation. The light curve was prewhitened with this model, which means that the Fourier series with frequency and number of harmonics that minimized $`\chi ^2`$ was subtracted from the data. Then periodogram was calculated for the prewhitened data. If the periodogram of the prewhitened light curve showed presence of peaks 15 times greater than the median of the peaks the star was considered a multiperiodic candidate. The most probable value of the second period of variability was then determined in exactly the same way as in the case of the primary period. In a case of a star with two detected periods a two period model was fitted to the original light curve. The formula for such a fit is given in equation 3. During this phase not only the amplitudes and phases were computed. Both periods were adjusted as well. The values of two periods of pulsation found in the earlier phases were the input values for the procedure. In the last step the model with two periods was subtracted from the data and search for additional periodicities was performed.

$$F(t_k)=A_{00}+\mathrm{\Sigma }_{j=1}^mA_{0j}\mathrm{cos}(j\omega _2t_k+\varphi _{0j})+\mathrm{\Sigma }_{i=1}^m\mathrm{\Sigma }_{j=im}^{mi}A_{ij}\mathrm{cos}[|i\omega _1+j\omega _2|t_k+\varphi _{ij}]$$
(3)
Number m is the rank of the fit. In a case where *m*=2 the series includes the constant term and the terms with frequencies $`\omega _1`$, $`\omega _2`$, 2$`\omega _1`$, 2$`\omega _2`$, $`\omega _1+\omega _2`$, $`|\omega _1\omega _2|`$. Amplitude errors are computed as well, and only the terms with amplitudes 3 times greater than their errors are included in the final fit. Value of $`\chi ^2`$ of a two period model strongly depends on the rank of the fit. As an example Table 1 demonstrates how $`\chi ^2`$ values change with $`m`$ for a two period model of an RRd star. It this star a model with two periods of variability is satisfactory, and fits the data very well for $`m4`$. In the case shown in Table 1 setting $`m>4`$ does not reduce $`\chi ^2`$, as the high order terms do not improve the fit, while consuming additional degrees of freedom.

## 4 Results. Statistics of multiperiodic variables.

As the result I obtained an RRab sample containing 1942 stars and an RRc sample of 771 stars. This gives the total number of RR Lyraes equal to 2713, which is much larger than 215 stars in OGLE-I RR Lyrae database. Although my current sample is not as numerous as LMC sample of Alcock et al. (1996) it reveals a larger diversity of pulsation behaviors. Both RRab and RRc stars cover a wide range of periods with $`P_{ab}[0.340,0.922]d`$ and $`P_c[0.201,0.446]d`$. Average values of pulsation periods in both classes are respectively $`P_{ab}^{avr}=0.554d`$ and $`P_c^{avr}=0.308d`$. Figure 1 shows period histogram for all RR Lyrae stars. Presence of two peaks reflects the fact of existence of two distinct populations of RR Lyrae stars. These populations are fundamental mode pulsators RRabs and first overtone pulsators RRcs.

Another interesting observation is a very high number of RR Lyrae stars that exhibit multiperiodic behavior. I discovered about 800 such stars. Most of them are Blazhko stars but there are also a few RRd stars, that is double mode pulsators oscillating simultaneously in fundamental mode and first overtone. Multiperiodic RR Lyrae stars will be discussed more thoroughly in the next subsections.

### 4.1 Double mode pulsators

RRd stars are very important in astrophysics because they provide us with one additional observable, that is the second period of variability. A widely known and intriguing fact concerning these variables is that they are very numerous in some stellar systems while almost completely absent in others. This fact is not understood. So far there is no clear evidence that any specific physical quantity is responsible for this difference. Moskalik (2000) suggests that different incidences of RRd stars in stellar systems could result from different metalicities in those systems. Many RRd variables were discovered in the Galactic disk and in some globular clusters in the Galaxy. Magellanic Clouds are also known to have quite numerous populations of RRd stars. Alcock et al. (2000) report 181 such variables in their LMC database. The fundamental question which arises from the analysis of OGLE data is why do we observe so very few RRd stars in the Galactic bulge. Previous investigations of OGLE-I database performed by Moskalik (2000) and Pigulski et al. (2003) resulted in only two RRd variables found. Those were the only RRd stars known in this system until now. Analysis of OGLE-II database, which includes the fields of OGLE-I, led to the discovery of three additional RRd stars in the center of the Galaxy. Total number of only 5 known RRd variables in the Galactic bulge is in sharp contrast to the large number of objects of this kind observed in globular clusters and in the Magellanic Clouds. Incidence rate, defined as the ratio of the number of RRd to the number of RRc variables, is $`0.007`$ for the bulge whereas for the LMC Alcock et al. (2000) give the value of $`0.134`$. This is obviously a striking difference: there are proportionally 20 times more RRd stars in the LMC than in the Galactic bulge. Unfortunately we still do not know why RRds are so rare in the center of the Galaxy.

Two of three RRd stars in OGLE-II database are very typical representatives of their class. The ratio of the periods of the first overtone and the fundamental mode, hereafter referred to as $`\zeta `$, is $`0.745`$ in both cases. The ratio of the amplitudes of the second periodicity to the dominant periodicity, although substantially different for both stars, is in good agreement with the statistics given by Alcock et al. (2000). This ratio is denoted as $`\xi `$. Table 2 gives the basic characteristics of the three RRd stars I discovered, that is both periods $`P_0`$ and $`P_1`$, values of parameters $`\zeta `$ and $`\xi `$ and two values of $`\chi ^2`$. The first value $`\chi _1^2`$ corresponds to a single period model with only the dominant period present. The second value $`\chi _2^2`$ is the $`\chi ^2`$ of a double period model with both periods of variability. These numbers are given for comparison to demonstrate that taking into account the second periodicity significantly reduces $`\chi ^2`$ of the fit. As it can be seen in table 2, the last star, namely bul\_sc39.1568, differs from the first two objects. It has a shorter period and a greater $`\xi `$ value than the two other stars. It means that this star is separated from the other two on the Petersen diagram, and therefore, according to Popielski et al. (2000), has a different metallicity.

### 4.2 Blazhko stars

Blazhko effect manifests itself in two different ways. Prewhitened periodograms show either one or two peaks in the close vicinity of the primary period of pulsation. As for now it is thought that these two behaviors divide Blazhko stars into two classes: BL1 class with frequency doublet and BL2 class with frequency triplet. Members of BL1 class have their additional frequency very close to the frequency of the radial mode, typically closer than 0.01 \[c/d\]. In the case of BL2s the dominant frequency and the two additional frequencies create an equidistant triplet centered on the dominant frequency. This triplet is typically asymmetric, that is the heights of the additional peaks are different. So far we do not have a satisfactory and complete explanation of Blazhko effect. Most widely accepted theoretical approach explains the presence of additional peaks by excitation of nonradial modes. Those modes are very close in the frequency domain to the frequency of the radial mode. Existence of equidistant, almost symmetric frequency triplets could be explained by other means but presence of strongly asymmetric triplets and especially doublets seem to require some kind of nonradial oscillations.

It was thought for a long time that Blazhko effect exists only among RRab stars. Recent studies indicate that it also concerns RRc variables. Moskalik & Poretti (2002) report the discovery of two RRc Blazhko stars. My research confirms that Blazhko effect is present in some RRc stars, although is not as frequent as in the RRab class. Two next subsections cover consecutively RRab and RRc Blazhko stars found in the OGLE-II database.

#### 4.2.1 RRab stars

My code identified 243 BL1 and 143 BL2 type stars among RRab stars.
BL1s constitute $`12.5\%`$ of all RRab in the OGLE-II database, while BL2s only $`7.4\%`$. Figure 4 shows two examples of such variables, one for BL1 and one for BL2 classes.
A purely radial pulsator, with no additional frequencies, was included in the figure for comparison. Average frequency separation between the additional peak and the radial mode is |$`\mathrm{\Delta }f=0.01`$|. In the case of BL2s, that is stars with equidistant triplets, $`\mathrm{\Delta }f`$ is understood as the difference between the frequency of the higher of the two subsidiary peaks and the frequency of the radial mode. In the literature one can also find a quantity called modulation period, which is the inverse of the frequency separation: $`P^m=\frac{1}{|\mathrm{\Delta }f|}`$.

Basic properties of RRab Blazhko stars are as follows:

* BL1

+ $`P_{min}`$=0.370d $`P_{max}`$=0.779d $`P_a`$=0.543d
    + $`P_{min}^m`$=5d $`P_{max}^m`$=968d $`P_a^m`$=91d
    + $`\mathrm{\Psi }_{min}`$=0.039 $`\mathrm{\Psi }_{max}`$=0.537 $`\mathrm{\Psi }_a`$=0.160
    + $`\mathrm{\Delta }`$f\>0 in 81$`\%`$ cases
* BL2

+ $`P_{min}`$=0.362d $`P_{max}`$=0.680d $`P_a`$=0.521d
    + $`P_{min}^m`$=15d $`P_{max}^m`$=977d $`P_a^m`$=96d
    + $`\mathrm{\Psi }_{min}`$=0.067 $`\mathrm{\Psi }_{max}`$=0.595 $`\mathrm{\Psi }_a`$=0.174
    + $`\psi _{min}`$=0.429 $`\psi _{max}`$=0.996 $`\psi _a`$=0.795
    + $`\mathrm{\Delta }`$f\>0 in 73$`\%`$ cases

where $`P`$ represents radial pulsation period, $`\mathrm{\Psi }`$ is the amplitude ratio of the additional peak to the peak representing the radial mode, and $`\psi `$ is the amplitude ratio of the two additional peaks in the case of BL2 stars. Subscripts min, max, a denote minimal, maximal and average values respectively. It can be seen that the amplitudes of additional peaks can be comparable with the amplitude of the dominant peak, and that the triplets range from highly symmetric to highly asymmetric. The fact that $`\mathrm{\Delta }f`$ is positive in most cases means that the additional peak in BL1s, or higher of the two additional peaks in BL2s, corresponds to a higher frequency than the frequency of the radial mode. This finding is in agreement with the results of Moskalik & Poretti (2002). As it will be shown in the next paragraph, RRc Blazhko stars behave differently.

Period distribution of RRab Blazhko stars is consistent with period distribution of all RRab stars. Figure 3 demonstrates period distributions of all RRab and RRab BL stars.

#### 4.2.2 RRc stars

RRc Blazhko stars are less numerous than RRab Blazhko stars. I detected 22 BL1 and 30 BL2 RRc stars, and it is easy to see that the proportions are different than those for RRabs. RRab Blazhko stars seem to prefer BL1 subclass, such effect is not observable among RRc stars. The statistical properties of RRc Blazhko stars differ from those of RRabs and are as follows:

* BL1

+ $`P_{min}`$=0.246d $`P_{max}`$=0.390d $`P_a`$=0.288d
    + $`P_{min}^m`$=7d $`P_{max}^m`$=986d $`P_a^m`$=247d
    + $`\mathrm{\Psi }_{min}`$=0.049 $`\mathrm{\Psi }_{max}`$=0.820 $`\mathrm{\Psi }_a`$=0.279
    + $`\mathrm{\Delta }`$f\>0 in 54$`\%`$ cases
* BL2

+ $`P_{min}`$=0.256d $`P_{max}`$=0.439d $`P_a`$=0.316d
    + $`P_{min}^m`$=2d $`P_{max}^m`$=913d $`P_a^m`$=397d
    + $`\mathrm{\Psi }_{min}`$=0.054 $`\mathrm{\Psi }_{max}`$=0.494 $`\mathrm{\Psi }_a`$=0.178
    + $`\psi _{min}`$=0.426 $`\psi _{max}`$=0.999 $`\psi _a`$=0.844
    + $`\mathrm{\Delta }`$f\>0 in 33$`\%`$ cases

with the symbols having exactly the same meaning as in the case of RRab stars. As one can see there are substantial differences between RRab and RRc Blazhko stars. Average modulation periods of RRc BLs are much longer and there is no strong preference for $`\mathrm{\Delta }f>0`$ among RRc BL1s. In fact for RRc BL2s we can observe an opposite preference, $`\mathrm{\Delta }f`$ is mostly negative, which means that the higher of the additional peaks corresponds to a lower frequency than the radial mode. Also the value of $`\mathrm{\Psi }_{max}=0.820`$ for RRc BL1s seems to be an interesting result, as it shows, that in the case of RRc Blazhko stars the nonradial mode can have almost the same amplitude as the radial mode. Examples of RRc BL1s and BL2s are shown in figure 5.

### 4.3 Stars with multiple periods

My sample of multiperiodic RR Lyrae stars contains some objects that can not be classified as either Blazhko type 1 stars BL1s or type 2 BL2s on the grounds of current definitions of these classes. Although their light curves look very similar to light curves of Blazhko stars, with the characteristic scatter around maximum and minimum of brightness, prewhitened periodograms of these objects reveal presence of more than two additional peaks. There are also a few objects with frequency triplets, but unlike in the case of BL2 stars, those triplets are not equidistant. Existence of RR Lyr stars with such signatures in the prewhitened spectra was reported by Moskalik & Poretti (2002), who found a few stars of this type in OGLE-I database. OGLE-II database contains much more examples. I found 86 such objects among RRab stars and 41 among RRc stars. This gives respectively 4.4$`\%`$ and 5.3$`\%`$ incidence rates in both classes. Examples are shown in figure 6.

It is not clear yet if these types of objects could also be referred to as Blazhko stars. It is possible that in the future they will be called BL3 stars. They certainly share some features with Blazhko stars, but on the other hand constitute a distinct group. Their nature is uncertain, but they seem to be another strong evidence for excitation of nonradial modes.

### 4.4 Objects of uncertain nature

The last class of my multiperiodic RR Lyrae catalog consists of objects that could not be classified reliably as members of any of the classes described above. Most of these objects are probably BL1 or BL2 stars, but computed modulation periods are longer than the time span of the data. Additional observations in the future are needed to verify the hypothesis that these objects are BL1 or BL2 stars with very long modulation periods. My catalog contains 14 such objects in the RRab class and 2 objects in the RRc class.

## 5 Summary

I presented the results of an automated search and analysis of RR Lyrae stars in the Galactic Bulge. My RR Lyrae samples are very numerous, about 12 times richer than OGLE-I RR Lyrae catalogs, and are therefore well suitable for statistical analysis. Out of more than 2700 RR Lyrae stars, that I identified in the OGLE-II database, about 600 proved to be multiperiodic. Most of the multiperiodic RR Lyrae stars are Blazhko stars. They show a large variety of pulsation periods and frequency patterns. The incidence rate of Blazhko stars in the Galactic Bulge, determined from my analysis, is very high comparing to some other systems like the LMC and the SMC. With about 25$`\%`$ in RRab and 10$`\%`$ in RRc classes, incidence rates are almost 2.5 times greater than those for the LMC found by Soszynski et al. (2003). Nature of this phenomenon is unknown and is yet to be understood. High number of about 100 RRc Blazhko stars is especially interesting because until now not many variables of this kind were known. They are very different from RRab Blazhko stars in certain aspects.

On the other hand I have found only 3 RRd stars in addition to the 2 previously known. This means that their incidence rate, with respect to RRc stars, is $`0.007`$, and is about 20 smaller than their incidence rate in the LMC given by Alcock et al. (2000).

As for now we still do not understand why there are so many more Blazhko stars and so many less RRd stars in the Galactic Bulge comparing to the Magellanic Clouds. These striking differences call for an explanation.

## 6 Acknowledgments

I would like to thank Dr. Wojciech Dziembowski and Dr. Andrzej Udalski for their helpful comments on the script and valuable discussions. This work was supported by Polish KBN grant 5 P03D 030 20.


text:
# Short range repulsive interatomic interactions in energetic processes in solids

## I Introduction

A high-energy ion (tens to hundreds of keV) propagating through the system can approach other atoms to very close distances, typically between 1 Å and 0.1 Å. The strong repulsive potential determines the scattering process. Ion implantation, laser ablation and displacement cascades produced by $`\alpha `$-decay events in materials used for nuclear waste immobilization, are examples in which these energetic collisions play an important role. The energetic ion loses energy due to inelastic collisions with the electrons and other atoms in the materialBohr . The latter, called nuclear stopping, dominates for relatively low energies ($``$ 1 keV/amu), while electronic stopping dominates for higher energies.

Nuclear stopping processes represent an enormous challenge because the kinetic energies involved allow for displacements of the atoms in large regions, and the energy dissipation takes place over long time scales (at least tens of picoseconds). Complicated structural configurations are obtained, with atoms breaking and making new bonds, and ionisation processes affecting the interactions between the particles. The complexity in physical and chemical interactions, the large sizes, and the long times involved in the problem demands high-efficiency first-principles calculations.

Many Density Functional Theory (DFT)DFT  implementations used by the condensed matter community are based in the pseudopotential approximation, in which the core electrons are considered as frozen in their atomic configuration, and only the valence electrons are responsible of the chemical and physical properties of the solid, resulting in a substantial reduction of the degrees of freedom to be solved and hence in much faster calculations than the equivalent all-electron. For the very challenging problems in this area, it would be highly desirable to asses the usefulness and the range of validity of such approximation. It has to be remembered, however, that conventional all-electron approaches are not much better suited for very short distances. If augmented plane waves (APW) are used to solve numerically the all-electron problem in the core region, problems arise when the augmentation spheres overlap. The need for an accurate basis set for the core electrons might prevent the use of available gaussian basis sets, normally devised for unperturbed core electrons.

Reference calculations are thus needed to test the accuracy of pseudopotential-based first principles calculations for describing electronic screening at short distances. Experimental uncertainties in the determination of the potentials are usually of the order of 10% or more. We have chosen to follow the approach of Nordlund et. al.Nordlund , and use as a reference calculations at the Hartree-Fock limit, where all these problems are eluded and the only remaining approximations are the absence of electronic correlation and of relativistic effects. In this work we have studied the relative importance of each of the effects that could be relevant to the problem, and assessed the validity of the pseudopotential approximation.

## II Methodology

For high energies the interatomic separation in a head-on collision can be small, making the binary collision (BC) an appropriate approximation. In this case, the complexities of the repulsive interaction between atoms can be modelled by a simple two-body interatomic potential for short distances. The interatomic potential $`V(R)`$ is defined as the difference between the total energy of the two atoms at a distance $`R`$, and the energy of the isolated atoms:

$$V(R)=E(R)E(\mathrm{})$$
(1)
We use two descriptions to characterize this potential: $`\delta V(R)`$, and the screening function $`\mathrm{\Phi }(R)`$. The former is defined as

$$\delta V(R)=V(R)\frac{Z_1Z_2e^2}{R}$$
(2)
where $`Z_1`$ and $`Z_2`$ are the atomic numbers of the two atoms. This can be used to characterize the potentials at $`r0`$. In the limit, we have the joint-atom solution, with just one atom with atomic number $`Z_1+Z_2`$. Notice that $`\delta V(0)=E^{atom}(Z_1+Z_2)E^{atom}(Z_1)E^{atom}(Z_2)`$, and we can use all-electron atomic calculations (without the problems related to basis commented before) to obtain information of this limit of the repulsive potential.

A different recasting of the interatomic potential that is also very usefulSmith , is obtained if we consider the effect of screening due to the electronic cloud around the nuclei:

$$V(R)=\frac{Z_1Z_2e^2}{R}\mathrm{\Phi }(R)$$
(3)
where $`\mathrm{\Phi }(R)`$ is defined as the screening function. For very short distances, there is almost no screening of the internuclear repulsion because the probability of finding electrons in the interatomic region is very small, and $`\mathrm{\Phi }(0)1`$. For long distances the potential goes to zero faster than the Coulomb potential and $`\mathrm{\Phi }(\mathrm{})0`$ (for neutral atoms). Some popular parametrizations of this potentialTorrens  are the ones due to Molière, Lenz-Jenson, Born-Mayer, and the one by Ziegler, Biersack and LittmarkZBL  (ZBL). We will compare our results with the latter, that takes the form:

$$\mathrm{\Phi }(R)=\underset{i=1}{\overset{4}{}}a_ie^{b_ir}$$
(4)
where the numerical coefficients $`a_i`$ and $`b_i`$ are fitted to a variety of interatomic repulsive potentials for different pairs. This gives a relatively good agreement with experiment, specially for light elements. Nevertheless, the deviation with experimental ion scatteringexpr  data can be as large as 20%.

The DFT energy $`E(R)`$ is obtained with the Siesta methodsiesta  using separablekb  norm-conserving Troullier-Martins pseudopotentialstm2 . The valence wave functions are expanded in pseudoatomic numerical orbitals including multiple-$`\zeta `$ and polarization function. Periodic boundary conditions require the definition of a supercell large enough for the interaction between replicas of the atoms to be sufficiently small (L$``$20Å). An uniform real-space grid defined by an equivalent plane-wave cutoff of 350 Ry was used for numerical integration.

In our simulations, scalar relativistic pseudopotentials are generated using different electronic configurations, including semicore states, partial core corrections, and different exchange-correlation functionals. The use of semicore shells, that is, core electrons considered as part of the valence, has the benefit of increasing the accuracy as compared to all-electron calculations, making the pseudopotential harder, and more demanding computationally. In table 1 the different configurations used for core/valence electrons are shown. Partial core correctionsPCC  (PCC) account at least partially for nonlinearity of the exchange-correlation functional due to the charge-density overlap between core and valence electrons. Notice that including scalar relativistic effects, as well as the correlation, represents, at least in some respects, an improvement to HF.

We compare the potentials obtained with Siesta with those obtained using the same numerical Hartree-Fock method as Nordlund et. al.Nordlund  and developed by Laaksonen et. al.2dhf . This HF approach, uses wave-functions of the general form:

$$\psi (\xi ,\eta ,\varphi )=\frac{e^{im\varphi }}{(2\pi )^{1/2}}u(\xi ,\eta )$$
(5)
and the Hartree-Fock equations are reduced to partial differential equations for the function $`u(\xi ,\eta )`$ in prolate spheroidal coordinates. The equations are then discretized using finite differences.

This numerical method is used as a reference in our electronic structure calculation. The effect of electronic correlation was shown not to be relevant, at least for the light element potential studied by Nordlund et. al.Nordlund . We confirm this result on other systems. We also checked that relativistic effects not included in the Hartree-Fock method, are important but not relevant for the purpose of this work. Thus, albeit the computational cost of HF restricts its application to relatively simple systems, it appears to be a good tool to gauge our approximation and test the accuracy achieved with other methods.

## III Results

The pairs considered in this work were C-C, O-O, Si-Si, Ca-O, and Ca-Ca with both HF and DFT, and U-O with only DFT. The screening functions $`\mathrm{\Phi }(R)`$ are shown in figures 1 and 2. Deviations from HF become relevant for distances $`R0.7`$Å, when the attractive bonding interaction begins to be important. It is in this region of bonding character in the interatomic potential where the different descriptions of exchange and correlation give different results. It is, however, not important for this study. We used the Ceperley-Alderca  parametrization of exchange-correlation in the Local Density Approximation (LDA), and the Perdew, Burke and Ernzerhofpbe  scheme for Gradient-corrected functional (GGA), both with and without spin polarization. We observed (Fig. 1) that any of these approximations for exchange and correlation gives essentially the same repulsive potential. The agreement with HF is good in general, even for short distances, showing that electronic correlation plays only a minor role in the repulsive potential. Note that ZBL gives a good description of the repulsive HF potential for light elements (C, O). For larger elements this agreement is not so good (see below).

### III.1 Analysis of the different approximations

#### III.1.1 Effect of semicore states

In Fig. 2 we analyze the influence of semicore states in the screening. When the core electrons of one atom begin to overlap with the core electrons of the other atom, the Pauli repulsion between them is not properly taken into account in the pseudopotential approximation and consequently the result deviates from the Hartree-Fock description. The use of semicore states in the pseudopotential improves this, with part of the core included as valence. This can clearly be seen in the top Fig.2a, comparing the HF results with those obtained with pseudopotentials. In 2b, we observe the influence of these semicore states in the hardness/softness of the repulsive potential. The absence of semicore electrons in the valence, makes $`V(R)`$ softer. Fig.2c shows the deviation of the potential with respect to the HF value. Although this result was expected, it is important to quantify it with respect to other effects, and to obtain an easy rule to assess the limit of validity of the pseudopotential approximation.

When the interatomic distance is comparable to bonding distances (R$``$1Å), the valence electrons are spread around the atoms and in the bonding region, screening the Coulomb interaction between the nuclei. As the distance between the nuclei is reduced, the repulsive electronic interactions (Coulomb and Pauli) in the bonding region increase, and the electrons are more unlikely to be in the confined area between the nuclei. This is shown in Fig. 3, where the valence charge distribution around the atoms is displaced from the interatomic center, to the surrounding area as the distance decreases. In a simulation within the pseudopotential approximation, the valence electrons have a smaller participation in the nuclear screening, but the core electrons are assumed to be fixed in their initial configuration, and thus, their contribution to the screening is roughly the same. For this reason $`\mathrm{\Phi }(R)`$ is always smaller than the expected screening function with an all electron approach.

Core electrons deep in the inner shells (see Fig.2d), are more strongly bonded to the nuclei and their participation in the collision would only be relevant for very high collision energies when the interatomic distance becomes really small. For that energy scale, the electronic stopping power is more important than the nuclear stopping power. Hence, an accurate description of the repulsive potential for shorter (higher) distances (energies), is beyond the scope of this study.

Only the outer core-shell electrons are included in the semicore. How localized these electrons are can be roughly estimated using the spread of the corresponding atomic Kohn-Sham orbitals in a DFT all-electron calculation. In table 2 we show for each atomic shell, the radii $`r_{90}`$ and $`r_{99}`$, in which 90% or 99% of the atomic wavefunction’s norm is localized. Two shells corresponding to different atoms will likely interact, if the interatomic distance is smaller than sum of the radii at which the charge is more localized. According to this, and considering $`r_{90}`$ as a reference, in a Si-Si collision, we can say that 2s electrons will interact with 2s or 2p valence electrons in the other atom for distances R$``$0.95Å, which is an appreciable distance (see arrows in Fig.2a). So, if we want to have a good description of the collision for shorter ranges, we would need to include 2s electrons into the semicore. With that, we would ensure a reasonable description of the interactions up to distances of the order of 0.2Å, when the 1s electrons start to play a role. This would explain why the pseudopotential without any semicore state gives such a soft repulsive potential.

The same argument may be used for 1s electrons in a C-C collision for distances smaller than 0.4Å, and for the 2s-2p electrons in a Ca-Ca collision below $``$0.6Å. In the Ca-Ca case, both the 3s and 3p shells have to be included in the semicore for reasonable screening down to 2$`r_{90}`$, in accordance to the model. In Si-Si, however, the 2p electrons seem to be enough, and the 2s apparently do not affect the screening for distances smaller than $``$1Å.

#### III.1.2 Effect of the core radius

In the pseudopotential approximation, it is assumed that the core electrons of different atoms do not overlap. For the interatomic distances of interest in this work, this is not the case. Although the core radii of the pseudopotentials ($`r_c`$) play an important role in the description of ground state properties of real materials, the effects for the repulsive potential at short distances is not as dramatic. Changes of up to 40% in $`r_c`$ give differences in the repulsive potentials smaller than 2% when the interatomic distance is smaller than 2$`r_c`$, whereas differences can be as large as 20% when the distance is close to 2$`r_c`$.

#### III.1.3 Effect of charge state, correlation, and relativistic corrections

The charge state of the shooting ion can play an important role in the description of the electronic energy loss for intermediate energies. For the short range interatomic repulsion, however, these effects are not so relevant: since electrons move out of the internuclear region for short interatomic distances (Fig. 3), there will be no important effect on the repulsive potential if the ionization involved not too deep electrons.

Correlation and relativistic corrections are known to be relevant for very short distances or very heavy nuclei. In the pseudopotential method, scalar relativistic corrections are included, but not in the HF approach that we use as a reference. With an all-electron DFT-based atomic calculation which includes relativistic effects (scalar and spin-orbit), we use $`\delta V(0)`$ to estimate an upper limit of the relativistic corrections. As an example a limit of $``$1.3 keV is obtained in a Ca-Ca collision. Although this is a considerable correction, the energy scales we are interested in (see Fig. 4), are of the order of tens of keV. Furthermore, the relativistic effects are expected to be smaller for distances $``$0.4Å. The difference between DFT and HF (both relativistic or both non-relativistic) is of the order of $``$5 eV and gives an idea of the magnitude of correlation effects. It is evident from the figure that the approximations involved in the pseudopotential approach are much more relevant than correlation or relativistic effects.

### III.2 Practical realization: Heavy ions

We can now try to obtain some qualitative information about the interatomic potential for heavy ions. The use of pseudopotentials for equilibrium (low energy) situations of lanthanide and actinide compounds has been shown to be accuratePickard ; Louie  when the effects coming from overlap between core and valence are properly treated (either with the inclusion of semicore states, or with partial core correction), even if only scalar relativistic corrections are considered. Relativistic spin-orbit terms can also be important in these elements, and the results obtained within this method should be taken with caution.

We have computed the uranium-oxygen interatomic potential $`V(R)`$ with different pseudopotentials, and compared it to ZBL (see Fig. 5). The electronic configuration is shown in table 1, and the localization radii $`r_{90}`$ and $`r_{99}`$, in table 3. Using the same arguments as before, we could expect that using a semicore with 6s and 6p electrons in the valence, the U-O interatomic potential obtained would be appropriate up to distances of the order of 0.7Å. Note that the ZBL potential would give a lower screening than expected (at least for distances larger than 0.7Å), giving a harder description of the binary collision.

### III.3 Correction of the screening function

A systematic feature of the pseudopotential calculations is that in the limit of small distances, the screening functions obtained do not approach unity. This is due to the frozen core approximation and the fact that the inner shells of the electronic cloud are not contributing to the screening. We propose a method improve $`\mathrm{\Phi }`$ by adding a continuous function that gives the proper limit ($`\mathrm{\Phi }(0)1`$) for values below the radius $`r_0`$ at which the core states begin to overlap. By inspection of figures 2 and 5 we could use the simplest function $`\mathrm{\Phi }(R)=e^{\alpha R^2+\beta R}`$, with parameters $`\alpha `$ and $`\beta `$ fixed by the conditions over $`\mathrm{\Phi }(r_0)`$ and $`\frac{d\mathrm{\Phi }}{dR}(r_0)`$, giving continuity of the screening and its derivative (and thus, the interatomic force) at $`r_0`$. Notice that the matching radii would depend on the particular configuration used for the generation of the pseudopotential. We propose to use the sums of the $`r_{90}`$ for the relevant core states of both atoms considered: $`r_0=r_{90}^{atom1}+r_{90}^{atom2}`$. In figure 6 we compare the previous screening functions with the ones corrected with this scheme. This will allow realistic molecular dynamic simulations with pseudopotentials, even if the atoms get extremely close in a high energy collision. Numerical fits of these potentials have already been used to do semiempirical simulations of collision cascades for studies of radiation damage in CaTiO3Kostya .

## IV Conclusions

The validity of the pseudopotential approximation has been checked at close distances as found in high energetic collisions in materials. We have compared the interatomic repulsive potential computed with an all-electron Hartree-Fock method, with that obtained with pseudopotentials. In this HF reference method, no approximation is assumed with respect to the core shells, or the basis used to represent the electronic states. The validity of the HF approximation itself has been assessed in order to use it as a reference for the pseudopotentials. The effect of different approaches to treat the electronic correlation was previously shownNordlund  not to give substantial differences in the screening function, a result confirmed in this study. The relativistic corrections, though important, can be considered small for the energy scales of interest in this work.

It has been shown that as the nuclear distance decreases, the electrons are expelled from the internuclear region to the surrounding area, and the screening is completely dominated by the deeper electrons, bound at short length scales. If the core is frozen at length scales longer or comparable to the internuclear distance, the screening is larger than expected and the pseudopotential description fails. We have shown this effect to be the most important one in the problem, actually the only one of relevance at the energy scales involved.

We have quantified to what extent the use of deeper electrons as part of the valence configuration is essential to describe the repulsive potential for relatively short internuclear distances. The localization of the core electrons, $`r_{90}`$, has been found as the most appropriate criterion used to decide which states have to be included for a particular range of collision energies. The pseudopotential description would then be valid until the more external core shells not included in the calculation start to overlap. Finally, we have proposed a scheme to correct the screening given in the pseudopotential approximation, that recovers the right trend for very short distances, and improves the use of pseudopotentials for the description of energetic collisions with first principles simulations.

## V Acknowledgments

This work was supported by British Nuclear Fuels (BNFL) and NERC. We would like to thank M. Dove, K. Trachenko for helpful discussions and sharing their experience with MD methods in radiation damage simulations.


text:
# The ATLAS Pixel Detector

## 1 Introduction

The ATLAS Inner Detector  is designed for precision tracking of charged particles with 40 MHz bunch crossing identification. It combines tracking straw tubes in the outer transition-radiation tracker (TRT) and microstrip detectors of the semiconductor tracker (SCT) in the middle with the Pixel Detector, the crucial part for vertex reconstruction, as the innermost component.

The Pixel Detector  is subdivided into three barrel layers in its centre, one of them around the beam pipe ($`r=5`$ cm), and three disks on either side for the forward direction. With a total length of approx. 1.3 m this results in a three-hit system for particles with $`|\eta |<2.5`$.

The main components are approx. 1700 identical sensor-chip-hybrid modules, corresponding to a total of $`810^7`$ pixels. The modules have to be radiation hard to an ATLAS life time dose of 50 MRad or $`10^{15}`$ neutron-equivalent.

## 2 Module Layout

A pixel module consists of a single n-on-n silicon sensor, approx. 2$`\times `$6 cm2 in size. The sensor is subdivided into 47,268 pixels which are connected individually to 16 front-end (FE) chips via “bumps” . These chips are connected to a module-control chip (MCC)  mounted on a kapton-flex-hybrid glued onto the back-side of the sensor. The MCC communicates with the off-detector electronics via opto-links, whereas power is fed into the chips via cables connected to the flex-hybrid.

To provide a high space-point resolution of approx. 12$`\mu `$m in azimuth ($`r\varphi `$), and approx. 90$`\mu `$m parallel to the LHC beam ($`z`$), the sensor is subdivided into 41,984 “standard” pixels of 50 $`\mu `$m in $`r\varphi `$ by 400 $`\mu `$m in $`z`$, and 5284 “long” pixels of $`50\times 600`$ $`\mu `$m2. The long pixels are necessary to cover the gaps between adjacent front-end chips. The module has 46,080 read-out channels, which is smaller than the number of pixels because there is a 200 $`\mu `$m gap in between FE chips on opposite sides of the module, and to get full coverage the last eight pixels at the gap must be connected to only four channels (“ganged” pixels). Thus on 5% of the surface the information has a two-fold ambiguity that will be resolved off-line.

The FE chips  contain 2880 individual charge sensitive analogue circuits with a digital read-out that operates at 40 MHz clock. The analogue part consists of a high-gain, fast preamplifier followed by a DC-coupled second stage and a differential discriminator. The threshold of the discriminator ranges up to 1 fC, its nominal value being 0.5 fC. When a hit is detected by the discriminator the pixel address is provided together with the time over threshold (ToT) information which allows reconstruction of the charge seen by the preamplifier.

## 3 Module Prototypes

Several ten prototype modules have been built with a first generation of radiation-hard chips in $`0.25\mu `$m-technology. In order to assure full functionality of the modules in the later experiment, measurements at the production sites, after irradiation, and in a test beam are performed.

### 3.1 Laboratory measurements

An important test that allows a large range of in-laboratory measurements is the threshold scan. Signals are created with on-module charge injection and scanning the number of hits versus the so injected charge yields the physical value of the threshold of the discriminator and the equivalent noise charge as seen by the preamplifier. A set of such scans is used to reduce the threshold dispersion by adjusting a DAC-parameter individually for each channel. The resulting threshold and noise after threshold tuning is shown in figure 1. Typically approx. 100 e threshold dispersion across a module and a noise value of below 200 e for standard pixels is achieved, as is needed for good performance.

In a similar fashion, the cross-talk is measured to a few per cent for standard $`50\times 400`$ $`\mu `$m2 pixels.

A measurement of the timewalk, i.e. the variation in the time when the discriminator input goes above threshold, is an issue since hits with a low deposited charge have an arrival time later than those with high charges, in particular for ganged pixels.

Data taken when illuminating the sensor with a radioactive source allows in-laboratory detection of defective channels. The source-spectrum reconstructed from the ToT-readings is in agreement with expectations.

### 3.2 Irradiation

Some of the prototype modules have been irradiated to a dose of 50 MRad, approx. the dose expected after 10 years of ATLAS operation. The radiation damage is monitored reading the leakage current individually for each pixel. The single event upset rate is measured during irradiation.

The threshold dispersion and the noise after irradiation as shown in figure 2 is only modestly increased and still well in agreement with requirements for operation in ATLAS.

### 3.3 Test beam measurements

Tests have been performed in the beam line of the SPS at CERN using 180 GeV/c hadrons. The setup consists of a beam telescope for the position measurement , trigger scintillators for timing measurement to 36 ps, and up to four pixel modules. The number of defective channels is observed to less than $`10^3`$ and for standard $`50\times 400`$ $`\mu `$m2 pixels the efficiency for normal incidence particles is 99.57$`\pm `$0.15%. The timewalk is measured to values similar to those from lab tests (see above).

Modules irradiated as described in section 3.2 have been tested in the beam line and the bias voltage needed for full depletion is measured to be between 500 and 600 V, see figure 3. The deposited charge is measured via the ToT readings and no significant changes in the uniformity w.r.t. unirradiated modules are observed.

### 3.4 Next generation of chips

The results of the prototype modules show that the chips fulfil largely the requirements of electrical performance and radiation hardness and that the production process has a high yield. Problems such as timewalk and too little SEU tolerance have been addressed and solved in a new design of FE and MCC chips .

## 4 Off-detector electronics

The off-detector readout electronics is designed to process data at a rate of up to 100 kHz level-1 triggers. The main data-processing component is the “read-out driver” (ROD), of which first prototypes have been built to pixel specifications and are being evaluated. The first-step event-building and error flagging is done via FPGAs. The communication to the rest of the DAQ-system is run through a 1.6 Gbit/s opto-link. The communication to modules, online monitoring and calibration runs are performed with SRAMs and DSPs; their programming is ongoing and modules have already been configured and operated successfully with a ROD.

## 5 System aspects

### 5.1 Support structures

The mechanics of the system has to guarantee good positional stability of the modules during operation while the amount of material has to be kept to a minimum. At the same time it has to provide cooling to remove the heat load from the modules and maintain the sensors at a temperature of -$`6^{}`$C to keep the radiation damage low.

Barrel-modules are glued to “staves”, long, flat carbon-structures with cooling pipes embedded. The staves are mounted inside halfshells, which themselves are assembled into frames to form the barrel system.

The disks are assembled from carbon-sectors with embedded cooling covering 1/12 of a wheel. The modules are glued directly to either side of the disk sectors.

### 5.2 Systemtests

First mini-systemtests have been performed with six modules on a disk sector and three modules on a barrel-stave. The noise behaviour on the disks or staves shows no significant differences compared to similar measurements with the same modules individually, see figure 4. Larger systemtests are already in preparation and will include realistic powering and read-out.

## 6 Conclusions

Prototype modules built with a first generation of radiation hard chips show largely satisfying performance in laboratory-test, in test beam studies and after irradiation. Remaining problems have been solved in a new generation of chips which is now ready for production.

Work on the off-detector electronics and the support structures have been going on in parallel and are well on track. First systemtest results are promising.


text:
# Transiciones de fase topológicas en cristales líquidos

## I Introducción

Mermin y WagnerMermin-Wagner  establecieron que para sistemas con grupo de simetría continuo no existe fase ferromagnética o de orden a largo alcance (LRO), a temperatura diferente de cero, para dimensiones menores o iguales que dos. Sin embargo, sistemas de este tipo pueden presentar un tipo de transición que está determinado por la aparición de defectos topológicos en pares a bajas temperaturas que justo en la transición, $`T_{KT}`$, se desliganBerezinskii ; Kost-Thou ; Kosterlitz . A este fenómeno se le llama transición de fase topológica o transición de Berezinskii, Kosterlitz y Thouless (BKT).

El modelo XY, con grupo de simetría $`O(2)`$, presenta este tipo transición en dos dimensionesKosterlitz . Entre sus características resaltantes tenemos que la fase LRO es sustituida, a bajas temperaturas, por un orden a cuasi-largo-alcance (QLRO). A muy bajas temperaturas, las correlaciones son dominadas por las ondas de espín obteniéndose una dependencia del exponente de la función de correlación, $`\eta `$, con la temperatura, $`T`$, de la forma: $`\eta =T/2\pi `$, donde la constante de Boltzmann $`k_B=1`$ y el factor de acoplamiento $`J=1`$. Debido a este comportamiento QLRO, a temperaturas menores que $`T_{KT}`$, la susceptibilidad, que mide las fluctuaciones de la magnetización, diverge. Otra característica de este tipo de transición es que a temperaturas justo por encima de la temperatura de transición $`t=(TT_{KT})/T_{KT}0`$, la longitud de correlación, $`\xi `$, diverge de una manera mucho más fuerte que la típica ley de potencias, $`\xi t^\nu `$, encontrada en las transiciones de segundo orden. Dicha divergencia es del tipo de una singularidad escencial: $`\xi \mathrm{exp}(bt^{1/2})`$.

Por otro lado, otro modelo de simetría continua como el de Heisenberg ferromagnético, con grupo de simetría $`O(3)`$, no presenta ningún tipo de transición de fase en $`d=2`$ mientras que si presenta la típica transición para-ferromagnética en $`d=3`$Polyakov ; Wolff2 . Por esta razón surgió la pregunta de que si sistemas con grupo de simetría no abelianos podían presentar transiciones del tipo BKT. Se ha reportado que el modelo de Heisenberg antiferromagnético completamente frustrado presenta una transición del tipo BKT pero sin correlaciones del tipo de ondas de espín a bajas temperaturaskawamura ; wintel .

Kunz y ZumbachKunzZumbach  realizaron un estudio intensivo en $`d=2`$, mediante simulaciones, del modelo $`RP^2`$ el cual consiste en un sistema con grupo de simetría global $`O(3)`$ pero con grupo de simetría local $`Z_2`$. El modelo $`RP^2`$ describe la transición nemática isotrópica de cistales líquidos en $`d=3`$. Ellos determinaron la longitud de correlación para $`t0`$ y encontraron un buen ajuste para una singularidad escencial. Por otro lado, mediante escalamiento de tamaño finito, estimaron que también era válido un ajuste del tipo ley de potencia. Sin embargo, en base a cálculos de energía y calor específico así como de ciertas cantidades que estiman el número de defectos topológicos, argumentaron que la transición de fase debería ser del tipo BKT.

En el 2003 se retomó el estudio pero ahora utilizando el modelo de Lebwohl-Lasher(LL)LL  para describir a los cristales líquidos en $`d=2`$. Este modelo ha sido muy exitoso para detectar la transición de fase discontínua débil que se observa en los experimentos de cristales líquidos en $`d=3`$. En este modelo se representan las moléculas mediante vectores unitarios $`\stackrel{}{\sigma }_w`$ colocados en los sitios de una red hipercúbica $`\mathrm{\Lambda }`$ de longitud $`L`$. Para este sistema el Hamiltoniano viene dado por:

$$\frac{H}{k_BT}=\frac{J}{k_BT}\underset{w}{}\underset{\mu }{}P_2(\stackrel{}{\sigma }_w\stackrel{}{\sigma }_{w+\mu }),$$
(1)
donde $`P_2`$ es el segundo polinomio de Legendre y la interacción es de primeros vecinos. En $`d=2`$, se determinóPLA  la presencia de un orden del tipo QLRO a bajas temperaturas con soluciones del tipo de ondas de espín cuando $`T`$ tiende a cero. Esto se hizo mediante el novedoso método de las transformaciones conformes. De esta manera se concluyó, que sistemas con grupo de simetría no abeliano tenían transición BKT con correlaciones del tipo de ondas de espín a bajas temperaturas.

En el presente artículo revisará el método de transformaciones conformes y su utilidad para determinar exponentes de la función correlación para sistemas con invarianza de escalas. Luego estudiando el parámetro de orden nemático bajo $`T_{KT}`$, mediante escalamiento de tamaño finito, se obtendrá de nuevo el exponente de la correlación y se mostrará la excelente concordancia entre los dos procedimientos. Finalmente se discutirá el efecto de la presencia de impurezas sobre este tipo de transición.

## II Método de las transformaciones conformes

Uno de los grandes problemas que existe al simular un sistema físico es el de tener evidencia de lo que ocurre en el límite termodinámico. Por lo general, en los estudios de fenómenos críticos, se simulan sistemas a diferentes tamaños y se estudian como escalan las cantidades termodinámicas en función de $`L`$. El exponente de escalamiento de dichas cantidades tiene relación con los exponentes críticos asociados a dichas cantidades. El costo computacional es demasiado alto. Recientemente, para sistemas en $`d=2`$ se ha comenzado a utilizar las transformaciones conformes (TC). Está técnica consiste en realizar simulaciones en sistemas finitos y conectar los resultados con los de un sistema infinito vía una TC. Esto lo podemos realizar en sistemas que presenten invarianza de escalas y para $`d=2`$.

Sistemas con invarianza de escala cumplen con la hipótesis de homogeneidad para la función de correlación de dos puntos para cualquier densidad $`\varphi `$, tal como el parámetro de orden, la energía, etc., $`\varphi (b\stackrel{}{r}_2)\varphi (b\stackrel{}{r}_1)=b^\eta \varphi (\stackrel{}{r}_2)\varphi (\stackrel{}{r}_1)`$, donde $`b`$ es un factor de escala y $`\eta `$ es la dimensión anómala o exponente de la correlación. De la misma manera, si existe invarianza de escala, podemos relacionar mediante la hipótesis de homogeneidad la función de correlación entre sistemas que se conectan bajo una transformación conforme. Supongamos que las simulaciones se realizan en un sistema cuadrado ($`w=u+iv`$) de tamaño $`L\times L`$ ($`L/2uL/2,0vL`$) y mediante una transformación de Schwarz-Christoffel lo mapeamos en el plano semi-infinito ($`z=x+iy,0y<\mathrm{}`$). Para esta transformación el cambio de escala es local y las funciones de correlación se relacionan mediante $`\varphi (w_1)\varphi (w_2)=|w^{}(z_1)|^{x_\varphi }|w^{}(z_2)|^{x_\varphi }\varphi (z_1)\varphi (z_2)`$, donde $`x_\varphi =\frac{1}{2}\eta `$.

Si las simulaciones se realizan en una red cuadrada con condiciones de borde fijas se obtiene que el perfil de densidad (función de correlación de un punto) se comporta como una ley de potencia de la forma: $`\varphi (w)_{\text{sq}}\kappa (w)^{\eta /2}`$ con $`\kappa (w)=\mathrm{}[z(w)](|1z(w)^2||1k^2z(w)^2|)^{1/2}`$, donde $`z(w)=\text{sn}(2Kw/L)`$, siendo $`K`$ la integral elíptica completa de primera clase, sn el seno elíptico de Jacobi y $`k`$ una constante. Esta expresión se obtiene ya que el perfil para una red semi-infinita con condiciones de borde fijas se conoce exactamente $`\varphi (z)_{\frac{1}{2}\mathrm{}}y^{x_\varphi }`$Cardy .

Gracias a la invarianza de escalas presente en las fases del tipo QLRO, la metodología de las TC se aplicó con exito en la determinación de los exponentes de la función de correlación del modelo XY en $`d=2`$ a $`t<0`$EPL ; JPA . Más recientemente, utilizando las TC en el modelo LL para los cristales líquidos en $`d=2`$, se concluyó la existencia de un orden QLRO en este sistemaPLA . Para este problema la función de correlación de un punto utilizada es: $`m(w)=P_2(\stackrel{}{\sigma }_w\stackrel{}{h}_{\mathrm{\Lambda }(w)})_{\text{sq}}`$ donde $`\stackrel{}{h}_{\mathrm{\Lambda }(w)}`$ indica que todas las moléculas en el borde $`\mathrm{\Lambda }(w)`$ tienen la orientación fija $`\stackrel{}{h}`$. En la Fig. 1(a) se muestran los perfiles del parámetro de orden como función de $`\kappa `$. Se nota, como a temperaturas menores que $`T_{KT}`$, se presenta un comportamiento del tipo ley de potencias mientras que muy por encima de esta temperatura dicho comportamiento es abandonado indicando que el sistema se encuentra en una región donde no existe invarianza de escalas. La determinación del cambio de comportamiento se puede hacer de manera cuantitativa calculando el $`\chi ^2`$PLA . En la Fig. 1(b) se grafican los exponentes de la correlación obtenidos a partir de los ajustes de la Fig. 1(a) como función de $`T`$. A muy bajas temperaturas se observa el comportamiento lineal típico de las ondas de espín.

## III Escalamiento de tamaño finito

En el presente trabajo se realizaron simulaciones del modelo LL en dos dimensiones con condiciones de borde periódicas. Las simulaciones se hicieron para temperaturas comprendidas entre $`T=0.4`$ y $`T=0.58`$. Se simularon sistemas con longitudes comprendidas entre $`L=16`$ y $`L=512`$. Se empleo el algoritmo de WolffWolff . Se tomaron del orden de $`10^5`$ pasos de equilibración y $`10^6`$ pasos de promediación. Los tiempos de relajación para todos los tamaños y longitudes no sobrepasaron los $`200`$ pasos de Monte Carlo. En particular, se reportan estimados del exponente de la función de correlación del parámetro de orden utilizando la técnica de escalamiento de tamaño finito para el parámetro de orden.

Un buen parámetro de orden para describir la transición de fase nemática-isotrópica, en $`d=3`$, para el modelo LL es:

$$M=L^2\underset{w}{}P_2(\stackrel{}{\sigma }_w\widehat{n})=L^2\underset{w}{}P_2(\mathrm{cos}\theta _w),$$
(2)
donde $`\widehat{n}`$ es un vector unitario que indica la dirección preferencial y se le denomina director y $`\theta _w`$ es el ángulo entre $`\stackrel{}{\sigma }_w`$ y el director. $`M`$ tiende a $`1`$ a muy bajas temperaturas y a $`0`$ a altas temperaturas. En $`d=2`$, $`M`$ debería ser $`0`$ para todo $`T`$ ya que no puede haber rompimiento de la simetría continuaMermin-Wagner .

En la Fig 2.(a) se tiene una gráfica del parámetro como función de la temperatura. Se observa un valor finito del parámetro de orden a temperaturas bajas. Sin embargo, la tendencia es que disminuya al aumentar $`L`$. En la Fig. 2(b) se muestra el comportamiento del parámetro de orden con la longitud del sistema. Obsérvese que para temperaturas bajas el comportamiento es del tipo ley de potencia. Esto indica la presencia de invarianza de escalas para $`T<T_{KT}`$ o lo que es lo mismo un orden del tipo QLRO. Sobre $`T_{KT}`$ la caida es exponencial, comportamiento característico de una fase paramagnética. El exponente de la ley de potencia obtenido a bajas temperauras es el correspondiente a la función de correlación de un punto, es decir, $`\frac{1}{2}\eta `$. Colocando los exponentes provenientes de los ajustes realizados con los datos de la Fig. 2(b) en la Fig. 1(b) tenemos que la correspondencia entre las dos metodologías reportadas es excelente.

## IV Discusión

En el presente trabajo se reportan evidencias de que a bajas temperaturas existe un orden del tipo QLRO. La manera más eficiente de realizar estos estimados es mediante el uso de la técnica de las transformaciones conformes ya que solo se necesita hacer simulaciones en un solo tamaño de red, inclusive para los primeros resultados bastó el uso de una red pequeña ($`L=48`$). El costo computacional del método de escalamiento de tamaño finito es mucho mayor al tener que simular para muchas longitudes de red y tamaños mayores. Realizando estudios del escalamiento de la susceptibilidad se tiene un porcedimiento alternativo para el cálculo del exponente $`\eta `$ al igual de una forma muy precisa para determinar que en $`T_{KT}`$ la singularidad es escencialPRB . Por lo tanto, definitivamente se concluye que, para los cristales líquidos en $`d=2`$, se tiene un sistema con grupo de simetría global no-abeliano que presenta una transición BKT con orden QLRO a bajas temperaturas.

En la literatura se ha reportado que para sistemas con transicón de tipo BKT, la introducción de desorden en los enlaces es totalmente irrelevanteEPJB . Esto se debe a que debido a la singularidad escencial para $`t>0`$ el exponente del calor específico $`\alpha \mathrm{}`$ ($`2d\nu =\alpha `$) y debido al criterio de Harrisharris  los exponentes críticos del sistema puro no deberían cambiar. Sin embargo, debido a cambios que aparecen en la coordinación el valor de $`T_{KT}`$ debería depender de la intensidad, $`c`$, del desorden. Si se realizaran simulaciones del LL con enlaces aleatorios obtendríamos que el exponente de la correlación sería idéntico para todos los valores $`T_{KT}(c)`$.

Donde resulta muy interesante discutir el problema del desorden sería en $`d=3`$ con campo magnético aleatorio. Según el criterio de Imry y Ma ImryMa  la fase nemática debería desaparecer y podría ser sustituida por un orden del tipo QLROchakrabarti ; feldman .

## V Agradecimientos

Este trabajo cuenta con el apoyo del programa PCP Venezolano-Francés titulado “Fluidos Petroleros”.