text:
# In-situ diagnostics of the crystalline nature of single organic nanocrystals by nonlinear microscopy

## Abstract

We elucidate the crystalline nature and the three-dimensional orientation of isolated organic nanocrystals embedded in a sol-gel matrix, using a polarized nonlinear microscopy technique that combines two-photon fluorescence and second harmonic generation. This technique allows the distinction between mono-crystalline structures and nano-scale poly-crystalline aggregates responsible for incoherent second harmonic signals.

preprint: APS/123-QED

The optical properties of nanoparticles have recently attracted much attention. In addition to metallic and semiconductor nanoparticles which are now used as bio-markers and as the building blocks of nanostructured materials NanopartBio , their organic counterparts constitute an interesting alternative. Advances in molecular engineering have enabled the design of molecular structures of various resonances and symmetries with optimized one- and two-photon absorption cross sections 2PhotonDyes , or combining different optical properties such as luminescence and second harmonic generation (SHG) Blanchard . In addition, macroscopic molecular arrangements have been optimized using the tensorial oriented gas model NLOcrystals , which predicts that an enhancement of the SHG efficiency is expected from the non-centrosymmetric crystalline arrangement of efficient nonlinear molecules. Molecular nanocrystals can be therefore envisioned as a new class of multi-functional nano-scale materials. In the case of organic nanocrystals however, the traditional crystalline characterization techniques have raised many practical barriers due to their low concentration and fragility. Consequently, the elucidation of their crystalline nature has been so far indirect or averaged over a large number of nanocrystals Vallee .

In this letter, we show that two-photon nonlinear microscopy permits in-situ characterization of isolated organic nanocrystals grown in an amorphous sol-gel matrix. The diagnostic is based on polarization resolved two-photon excited fluorescence (TPF) and SHG. TPF is an incoherent process allowed in centrosymmetric media, which exhibits a specific anisotropy depending on the medium symmetry. On the other hand, SHG is the signature of a crystalline non-centrosymmetric phase in the sample, with a sensitivity down to the nanometric scale NLO-LB ; PrasadNanoCrystals . We show that the polarization analysis of both TPF and SHG from nanocrystals allows the unambiguous discrimination between isolated mono-crystalline and poly-crystalline systems. Moreover, once a nanocrystal has been identified as mono-crystalline, a detailed model for both TPF and SHG polarization responses accounting for the unit-cell symmetry allows the determination of its three-dimensional orientation within the host matrix.

The organic nanocrystals that we investigate are based upon the $`\alpha `$-((4’-methoxyphenyl)methylene)-4-nitro-benzene-acetonitrile) molecule (CMONS), which exhibits efficient luminescence and quadratic nonlinearity under two-photon excitation Wang ; Treussart ; Sanz . The bulk crystalline phases of such crystals have three possible polymorphic forms, two being non-centrosymmetric with a very similar unit-cell crystal geometry (forms $`(a)`$ and $`(b)`$)Sanz ; CMONScristallo ; Oliver . The preparation of the organic nanocrystals in sol-gel glasses relies on the control of the nucleation and growth kinetics of the dye confined in the pores of the gel Sanz . CMONS nanocrystals were grown in 1:1 tetramethoxysilane(TMOS):methyltrimethoxysilane(MTMOS) matrices with a controlled CMONS:alkoxydes molar fraction. The sols containing the alkoxyde solvent, water and organic phases are directly spin-coated on a glass substrate. The presence of nanocrystals in the sol-gel film after spin-coating was ascertained by measurement of the material melting point using differential scanning calorimetry. Previous works on such materials have shown through bulk spectroscopy characterizations that the nanocrystals are in the predominant thermodynamically stable $`(a)`$ form, with a luminescence maximum peaking around 500-580 nm Sanz . The particle mean size, which depends on the preparation parameters such as temperature and matrix porosity, ranges from 20 to 100 nm. In this work, we focus on the smallest particles, which are obtained using a CMONS:alkoxydes molar fraction of $`4\times 10^3`$.

CMONS nanocrystals immobilized in the sol-gel film (0.5 - 1 $`\mu `$m thickness) are imaged using an inverted two-photon microscopy set-up with an excitation wavelength of 987 nm from a Ti:Sa laser (150 fs pulses duration, 82 MHz repetition rate) LeFloch . The reflection geometry, which uses a dichroic mirror to provide the incident IR light onto the sample, is well suited for the observation of SHG from structures of sub-wavelength size, this optical process being in this case insensitive to nonlinear phase-mismatch effects. The SHG signal, at half of the incident wavelength (493.5 nm), is spectrally shifted from the TPF emission maximum (at around 570 nm in the studied films), and can therefore be detected separately using appropriate spectral filters 111Note that for all quantitative analysis, we account for residual TPF emission in the SHG detection channel.. TPF and SHG scanning microscopy images exhibit isolated spots (about $`1`$ per 10 $`\mu `$m2) with diffraction limited resolution (Fig. 1(a)). The polarized two-photon microscopy technique consists of rotating the incident IR linear polarization while recording the SHG or TPF signals along two perpendicular ($`X`$ and $`Y`$) polarization directions for each isolated nanocrystal. We will show that this simple scheme allows the distinction between mono-crystalline and poly-crystalline forms.

In order to model the polarization responses of both optical processes, we first assume a mono-crystalline CMONS nanocrystal whose unit-cell orientation is defined by the Euler set of angles $`\mathrm{\Omega }=(\theta ,\varphi ,\psi )`$, as illustrated in Fig. 1(b).

The CMONS crystalline unit-cell in its non-centrosymmetric form $`(a)`$ is represented in Fig. 2. The unit-cell is composed of four CMONS molecules arranged according to the monoclinic $`C_c`$ crystalline space group CMONScristallo . The macroscopic polarization responses from mono-crystalline nanocrystals are dependent on their unit-cell orientation $`\mathrm{\Omega }`$. The TPF optical response involves the molecular fourth order $`\gamma `$ and second order $`\alpha `$ susceptibility tensors, which participate respectively to the two-photon absorption and the one-photon emission processes. The TPF intensity analyzed in the $`I`$ polarization direction is proportional to LeFloch :

$`_I^{\text{TPF}}(\mathrm{\Omega })=\alpha _{II}(\mathrm{\Omega }){\displaystyle \underset{J,K,L,M}{}}\gamma _{JKLM}(\mathrm{\Omega })\overline{E_JE_KE_LE_M}`$ (1)
where the $`I,J,K,L,M`$ indexes span the macroscopic frame ($`X,Y`$), and $`\overline{(\mathrm{})}`$ denotes photodetection temporal averaging. The $`E_I`$ components are the incident field polarization projections on each macroscopic axis. The SHG optical response stems from the molecular third order $`\beta `$ susceptibility tensor. The resulting intensity analyzed in the $`I`$ polarization direction for a mono-crystalline structure originates from the coherent addition of induced nonlinear polarizations from each molecule in the crystal, and is therefore proportional to:

$`_I^{\text{SHG}}(\mathrm{\Omega })={\displaystyle \underset{J,K,L,M}{}}\beta _{IJK}(\mathrm{\Omega })\beta _{ILM}(\mathrm{\Omega })\overline{E_JE_KE_LE_M}`$ (2)
Using Eqs. (1), (2) and in-plane rotation of the incident polarization, each of the macroscopic coefficients $`\alpha _{II}`$, $`\beta _{JKL}`$ and $`\gamma _{IJKL}`$ can be deduced. These coefficients depend on the orientation angle $`\mathrm{\Omega }`$ and on the unit-cell susceptibility components $`\alpha _{ij}`$, $`\gamma _{jklm}`$ and $`\beta _{ijk}`$ with $`(i,j,k)=`$ (a,b,c). For any $`t`$ tensor with components $`t_{I\mathrm{}N}`$ expressed in the macroscopic frame, the $`\mathrm{\Omega }`$ dependence of these coefficients can be readily calculated from a projection law $`t_{I\mathrm{}N}(\mathrm{\Omega })=_{i,\mathrm{},n}t_{i\mathrm{}n}\mathrm{cos}(i,I)(\mathrm{\Omega })\mathrm{}\mathrm{cos}(n,N)(\mathrm{\Omega })`$, where the $`\mathrm{cos}(i,I)`$ functions of ($`\theta ,\varphi ,\psi `$) express the unit-cell frame projection in the macroscopic frame. Consequently, the exploration of polarization responses finally relies on the knowledge of the microscopic $`\alpha `$ and $`\gamma `$ tensor coefficients for TPF detection, and of the $`\beta `$ tensor coefficients for SHG. At the molecular scale, we modelled a single CMONS molecule by a rod-like system of susceptibility components denoted $`\alpha _{uu}^{\text{cmons}}`$, $`\beta _{uuu}^{\text{cmons}}`$ and $`\gamma _{uuuu}^{\text{cmons}}`$, $`u`$ defining the CMONS fundamental dipole direction in the molecular frame Wang . The unit-cell tensor coefficients are constructed from the addition of each molecular contribution accounting for their respective orientation in the unit cell, according to the oriented gas model. This model neglects possible contributions from intermolecular interactions in the susceptibility calculations NLOcrystals . Numerical values of the resulting $`\alpha `$, $`\beta `$ and $`\gamma `$ tensor coefficients in the unit-cell (a,b,c) frame are given in the table of Fig. 2, the few non-vanishing remaining coefficients being consistent with the unit-cell symmetry. The off-diagonal tensor coefficients are the signature of the multipolar symmetry of the CMONS nanocrystals, which strongly influences the polarization dependence of the macroscopic optical responses.

Fig. 3 shows a typical example of signals from two isolated nanocrystals, exhibiting various TPF and SHG polarization responses analyzed in the $`X`$ and $`Y`$ directions. These signals can be simultaneously fitted with the previous model, which is a clear evidence of mono-crystalline structures. Moreover, the fit of both TPF and SHG features permits the retrieval of the orientation parameters of such nanocrystals, using Eqs (1) and (2) 222The fits account for the large aperture of the microscope objective, as well as the slight ellipticity and dichroism resulting from reflection of the incident IR beam on the dichroic mirror, characterized separately by ellipsometry.. This example shows a typical situation where two nanocrystals of different orientations exhibit similar TPF polarization responses, with on the contrary very different SHG responses. This originates from the distinct tensorial nature of these two optical processes: TPF involves even order tensors and is sensitive to the axial order terms of the molecular distribution (symmetric contributions), whereas SHG is sensitive to odd order parameters (non-centrosymmetric contributions). This complementarity is therefore the key point in the determination of the 3D orientation of isolated nanocrystals 333The 3D orientation determination would however not be possible in the rare case of rod-like symmetry ($`C_\mathrm{}v`$) crystals, for which the projection of the unit-cell in the sample plane is of the same symmetry as the system itself.. In comparison, traditional fluorescence microscopy only measures the in-plane orientation of the emitters, while the access to out-of plane contributions requires the use of complex optical field geometries or variable incident angle measurements SingleMolecule ; Wild .

In the present measurements, the determination of the three ($`\theta ,\varphi ,\psi `$) angles is furthermore unambiguous. Indeed, a slight change of one of the three Euler angles (by typically $`\pm \mathrm{\hspace{0.17em}5}10^{}`$) would clearly modify differently the TPF and SHG polarization responses. Moreover another type of analysis, which consists of using the experimental TPF and SHG anisotropy ratios in order to determine the unit-cell orientation, gives similar results by ruling-out any solution that does not correspond to polar plots such as those shown in Fig. 3.

It is visible from the measurements on mono-crystalline structures of Fig. 3 that the TPF polarization responses in the $`X`$ and $`Y`$ analysis directions have the same shape, although they are of different amplitudes. This effect is also predicted mathematically by the model. A departure from this main characteristic is therefore evidence of a non mono-crystalline phase which contains at least two or more nanocrystals. Such behavior, as represented in Fig. 4, is typically observed for about 30% of the measured nanocrystals. The TPF and SHG features of Fig. 4(a) are identical and have a ($`X,Y`$) axial symmetry for crossed polarization excitations. Such symmetry, fixed by the detection polarization directions, is a signature of incoherent emission from an isotropic orientational distribution of emitters. This effect, which has been observed in Rayleigh scattering experiments in solutions HLS , has been only moderately investigated so far in the study of nanoparticles NLOgold . In the case of the CMONS nanoparticles, an incoherent SHG emission is made possible by their sub-wavelength size, contrary to micrometer-scale particles for which the SHG efficiency is expected to depend strongly on the particle size when it approaches the coherence length. Such an effect can occur from the nano-scale assembly of randomly oriented SHG-active nanocrystals formed during the crystallization process. An incoherent summation of the responses from $`N>1`$ nanocrystals of different orientations contained in a sub-wavelength size aggregate can be modelled by the following SHG or TPF intensities:

$`_I^{\text{incoherent}}`$ $`=`$ $`{\displaystyle \underset{N}{}}_I(\mathrm{\Omega }_N)`$ (3)
with $`(\mathrm{\Omega }_N)`$ being the SHG or TPF intensities from a nanocrystal of given orientation $`\mathrm{\Omega }_N`$ in the aggregate.

Computing the SHG and TPF responses with an increasing amount of randomly generated $`\mathrm{\Omega }_N`$ orientations, the completely symmetric features of Fig. 4(a) start appearing from an incoherent summation over $`N>10`$ nanocrystals. Less-symmetric features are the signature of an intermediate situation between mono-crystalline and isotropic assemblies of nanocrystals originating from a sum over 2 to 10 nanocrystals. In Fig. 4(b), a reasonable fit is obtained for $`N`$=$`2`$ nanocrystals in the aggregate, having different relative orientations. Our observations therefore allow the distinction between mono-crystalline and more complex poly-crystalline structures formed by either a nearly-isotropic aggregate of nanocrystals, or by a small number ($`N<10`$ for the current analysis) of nanocrystals with sub-wavelength sizes.

In conclusion, the combination of TPF and SHG in two-photon nonlinear polarized microscopy has allowed the examination of organic nanoparticles, with the possibility to clearly distinguish mono-crystalline structures from poly-crystalline arrangements. We have also shown that a quantitative model of the TPF and SHG polarization responses enables the 3D orientation of mono-crystalline nanocrystals imbedded in an amorphous host matrix. This method can be extended to a broad variety of complex structures. The study of both polarized SHG and TPF at the nanoscale level can be applied to investigate local field effects and local orientation of molecules under electrical or optical fields perturbation. It can furthermore be extended to orientational tracking of optical markers in complex environments.

The authors thank Pr. J.-F. Nicoud (IPCMS Strasbourg) for his assistance in the selection and preparation of the chromophore, V. De Beaucoudrey for experimental contributions, and R. Pansu for fruitful discussions. This work is partially supported by an ACI ”Jeune Chercheur” grant from Ministère de la Recherche, France.


format: text

text:
# The deconfinement phase transition in Yang-Mills theory with general Lie group 𝐺

## 1 Introduction and Overview

Let us consider Yang-Mills theory with gauge symmetry group $`G`$ on a periodic lattice in $`(d+1)`$-D with the Wilson action. This theory is invariant under the multiplication by an element of the center subgroup $`C(G)`$ of $`G`$ of all time-like links in a given time-slice. This global center symmetry is unbroken at low temperatures and it gets spontaneously broken at the deconfinement phase transition. The Polyakov loop $`\mathrm{\Phi }`$ transforms non-trivially under this symmetry and thus is an order parameter for the deconfinement transition.

Integrating out the spatial degrees of freedom of the Yang-Mills theory, one can write down an effective action for $`\mathrm{\Phi }`$. It describes a scalar model with global symmetry $`C(G)`$ in $`d`$-D: the $`C(G)`$-symmetric confined phase of the gauge theory corresponds to the disordered phase of the scalar model, while the $`C(G)`$-broken deconfined phase has its counterpart in the ordered phase. Svetitsky and Yaffe  argued that the interactions in the effective description are short ranged. Hence, if the deconfinement phase transition in the Yang-Mills theory is second order, approaching criticality, the details of the complicated short-ranged effective scalar action become progressively unimportant. Only the dimensionality $`d`$ and the symmetry $`C(G)`$ of the scalar model are relevant. Thus, one can exploit the universality of the critical behavior to use a simple scalar model to obtain information about the much more complicated Yang-Mills theory. If the phase transition is first order the correlation length does not diverge: there are no universal features and the $`G`$-symmetric Yang-Mills theory in $`(d+1)`$-D does not share the critical behaviour with a $`C(G)`$-symmetric $`d`$-D scalar model.

Yang-Mills theory on the lattice is naturally formulated in terms of group elements while in the continuum the fundamental field is the gauge potential, living in the algebra. An algebra can generate different groups, however it is natural to expect that lattice Yang-Mills theories whose gauge groups correspond to the same algebra have the same continuum limit. Hence, instead of $`SO(N)`$ we consider its covering group $`Spin(N)`$. Keeping this in mind, we look at the center subgroups $`C(G)`$ of the various simple Lie groups $`G`$

$`C(SU(N))=𝖹𝖹(N);C(Sp(N))=𝖹𝖹(2)`$ (1)
$`C(SO(N))C(Spin(N))=\{\begin{array}{c}𝖹𝖹(2);N\text{odd}\hfill \\ 𝖹𝖹(2)^2;N=4k\hfill \\ 𝖹𝖹(4);N=4k+2\hfill \end{array}`$ (5)
$`C(G(2))=C(F(4))=C(E(8))=\{\mathrm{𝟣}\mathrm{𝟣}\}`$ (6)
$`C(E(6))=𝖹𝖹(3);C(E(7))=𝖹𝖹(2)`$ (7)
Many numerical simulations in $`(2+1)`$-D and $`(3+1)`$-D have been performed for $`SU(N)`$ Yang-Mills theory in order to investigate the order of the deconfinement transition and – if it is second order – to check the validity of the Svetitsky-Yaffe conjecture. A recent study on this subject can be found in . The currently known results are:
$``$ $`(3+1)`$-D: for $`N=2`$ – consistent with the Svetitsky-Yaffe conjecture – the deconfinement transition is second order, in the universality class of the $`3`$-D Ising model; for $`N=3,4,6,8`$ the phase transition is first order with no universal features.
$``$ $`(2+1)`$-D: fluctuations are stronger and the deconfinement phase transition turns out to be second order for $`N=2,3,4`$. In agreement with the Svetitsky-Yaffe conjecture, the $`2`$-D universality classes are, respectively, those of the Ising, 3-state Potts, and Ashkin-Teller models.

The $`SU(N)`$ branch is not a good choice to study the relation between the order of the deconfinement transition and the size of the group. In fact, when increasing the size $`N^21`$, also the center $`𝖹𝖹(N)`$ changes. In order to disentangle these two features we have considered $`Sp(N)`$ Yang-Mills theory. In this case, the available universality class is fixed and the relevance of the size of the group on the order of the deconfinement transition can be directly addressed.

## 2 $`Sp(N)`$ Yang-Mills Theory

The matrices $`USU(2N)`$ with the property

$$U^{}=JUJ^+,\text{where}J=i\sigma _2\mathrm{𝟣}\mathrm{𝟣}$$
(8)
form the symplectic group $`Sp(N)`$; $`\sigma _2`$ is the imaginary Pauli matrix and $`\mathrm{𝟣}\mathrm{𝟣}`$ is the $`N\times N`$ unit matrix. From the definition (8), it follows that the Hermitean generators $`H`$ of $`Sp(N)`$ satisfy the constraint $`H^{}=JHJ^+=JHJ`$. Then we can write the following general forms for $`U`$ and $`H`$

$$U=\left(\begin{array}{cc}W& X\\ X^{}& W^{}\end{array}\right),H=\left(\begin{array}{cc}A& B\\ B^{}& A^{}\end{array}\right)$$
(9)
where $`A,B,C,D`$ are $`N\times N`$ complex matrices satisfying: $`A=A^+`$, $`B=B^T`$, $`WW^++XX^+=\mathrm{𝟣}\mathrm{𝟣}`$ and $`WX^T=XW^T`$. Counting the number of degrees of freedom, it follows that $`Sp(N)`$ has $`N(2N+1)`$ generators and its rank is $`N`$. For the center elements – since they are multiples of the unit matrix – it holds that $`X=0`$ and $`W=W^{}`$: hence $`C(Sp(N))=𝖹𝖹(2)`$. We note that (8) implies that $`Sp(N)`$ is a pseudo-real group with the special cases $`Sp(1)=SU(2)`$ and $`Sp(2)=Spin(5)`$. The $`Sp(N)`$ Yang-Mills theory on the lattice can be formulated in the usual way in terms of group-valued link variables. We have carried out our simulations using the Wilson plaquette action.

## 3 Numerical Results

In the next two subsections we report on the results of numerical simulations in $`Sp(2)`$ and $`Sp(3)`$ Yang-Mills theories in $`(2+1)`$-D and $`(3+1)`$-D.

### 3.1 $`Sp(2)`$ Yang-Mills Theory

As a first step we have scanned the expectation value of the plaquette from the strong to the weak coupling regime. We find no bulk phase transition that might interfere with the study of the deconfinement transition. In $`(2+1)`$-D we observe a second order deconfinement transition, signalled by the broadening of the probability distribution of $`\mathrm{\Phi }`$ and, hence, by the increase of the Polyakov loop susceptibility $`\chi `$ at criticality. A finite size scaling analysis confirms the expectation that the universality class is that of the $`2`$-D Ising model. Fig. 1 shows the collapse of $`\chi `$ data – collected on lattices of different sizes $`L^2\times 2`$ and at various gauge couplings $`\beta `$ – on a single curve. We also plot rescaled $`\chi `$ data for $`SU(2)`$ Yang-Mills theory in $`(2+1)`$-D: it again deconfines with a second order transition in the 2-D Ising universality class. The two sets agree excellently.

In $`(3+1)`$-D – contrary to what one might have expected – the probability distribution of $`\mathrm{\Phi }`$ in the critical region shows the coexistence of the symmetric and of the broken phases, indicating that the deconfinement transition is first order. Fig. 2 indeed shows that $`\chi `$ scales with the spatial volume $`L^3`$ at criticality.

### 3.2 $`Sp(3)`$ Yang-Mills Theory

The results of $`Sp(2)`$ Yang-Mills theory show that in $`(2+1)`$-D fluctuations are stronger than in $`(3+1)`$-D and the deconfinement transition is second order. Expecting that the larger the group the weaker the fluctuations, we have also investigated the deconfinement transition in $`Sp(3)`$ Yang-Mills theory. Consistent with this picture, we find that $`(2+1)`$-D $`Sp(3)`$ Yang-Mills theory has a first order deconfinement transition. The probability distribution of $`\mathrm{\Phi }`$ in the critical region, indeed displays the coexistence of the broken and of the symmetric phases. In $`(3+1)`$-D – similar to the $`Sp(2)`$ case – $`Sp(3)`$ Yang-Mills theory deconfines with a first order transition. Fig. 3 shows tunneling events between the coexisting symmetric and broken phases as a function of Monte Carlo time $`t_{\text{MC}}`$.

$`Sp(3)`$ Yang-Mills theory also has no bulk phase transition in $`(2+1)`$-D and $`(3+1)`$-D.

## 4 Conjecture

Our numerical results show that $`(2+1)`$-D $`Sp(2)`$ Yang-Mills theory has a second order deconfinement transition in the $`2`$-D Ising universality class. However, $`(3+1)`$-D $`Sp(2)`$, $`(2+1)`$-D and $`(3+1)`$-D $`Sp(3)`$ Yang-Mills theories deconfine with a non-universal first order transition. Hence, despite the fact that a universality class is available, Yang-Mills theory can have a non-universal first order deconfinement transition. A non-trivial center plays no role in determining the order of this transition. Instead our $`Sp(N)`$ and the $`SU(N)`$  results indicate that the order of the deconfinement transition is a dynamical issue related to the size of the gauge group. We conjecture that the difference in the number of the relevant degrees of freedom between the confined phase (color singlet glueballs) and the deconfined phase (gluon plasma) determines the order of the deconfinement transition. Thus, we expect that in $`(3+1)`$-D only $`SU(2)`$ Yang-Mills theory has a second order deconfinement transition; in $`(2+1)`$-D, due to stronger fluctuations, only $`SU(N)`$, $`N=2,3,4`$, and $`Sp(2)`$ Yang-Mills theories should have second order transitions. According to this picture, $`E(6)`$ and $`E(8)`$ Yang-Mills theories should also have a first order transition due to the large size of the groups: 78 and 133 generators, respectively. For Yang-Mills theories with trivial center gauge groups $`G(2)`$, $`F(4)`$, $`E(8)`$ there is no compelling reason for a finite temperature phase transition . Then, although a first order transition can not be ruled out, we expect a crossover.


text:
# Vertex Intrinsic Fitness: How to Produce Arbitrary Scale-Free Networks

(September 29, 2003)

## Abstract

We study a recent model of random networks based on the presence of an intrinsic character of the vertices called fitness. The vertices fitnesses are drawn from a given probability distribution density. The edges between pair of vertices are drawn according to a linking probability function depending on the fitnesses of the two vertices involved. We study here different choices for the probability distribution densities and the linking functions. We find that, irrespective of the particular choices, the generation of scale-free networks is straightforward. We then derive the general conditions under which scale-free behavior appears. This model could then represent a possible explanation for the ubiquity and robustness of such structures.

In the last few years, much attention has focused on the study of complex networks. A network is a mathematical object consisting of a collection of vertices (nodes) connected by edges (links)Bollobas01 ; DM03book . Networks arise in many areas of science: biology KKKR02 ; bio ; bio2 , social sciences ASBS00 ; social1 ; social2 , Internet int1 ; int2 ; int3 , WWWAH00 , etc., where vertices and links can be for example, proteins and their mutual interaction, individuals and sexual relationship Liljeros01 , computers and cable connections. Very interestingly the same non trivial statistical properties appear ubiquitously in all the above situations. A more traditional view, indeed, is represented by the binomial model inspired to the random graph model of Erdős-Rényi ER61 . Here, each vertex has the same probability to connect to any other, resulting in a network whose degree, i.e. number of edges per vertex, has a binomial distribution. This is not the case of the above real data, where instead, the structure is self similar resulting in a scale-free (SF) probability distribution for the degree. More specifically, the degree $`k`$ of the vertices, i.e. the number of links entering them, is distributed according to a power law $`P(k)k^\alpha `$ with usually $`3<\alpha <2`$.

In order to explain the occurrence of SF networks the ingredients of growth preferential attachment have been introduced bara . The network increases the number of vertices with time, the newcomers tend to be connected with old vertices with large degree. Nevertheless, in some cases, we have the same SF properties without either growth of the system or preferential attachment mechanism. As an example, the finite set of protein interactions in a cell forms a self-similar network. This is done without growth of the system size and ignoring their reciprocal degree. Possibly, some external influence on intrinsic properties like chemical affinity is instead driving the phenomenon.

To take into account this new mechanism, the varying fitness model has been introduced by Caldarelli et al. gcalda . In this model, considering e.g. only undirected graphs, one extracts a real non-negative variable $`x`$ (the hidden variable) for each vertex of the graph from a probability distribution density $`\rho (x)`$. This variable $`x`$ is the fitness of the vertex. Links between vertices are successively formed with a probability function $`f(x,y)`$, a symmetric function of its arguments.

A static simplified form of the vertex hidden variable model has been considered for only one particular case by Goh et al. goh01 . Bianconi et al. BB01b  introduced a fitness mechanism coupled to the preferential attachment. In the paper of Caldarelli et al. gcalda , the onset of SF behavior is instead directly related only to the fitness presence of any kind. This behavior is also checked for different fitness probability distribution densities. Following Ref. gcalda , we present here an exhaustive study on the conditions needed in order to produce a SF network.

The aim of this work is to provide some ingredients to generate SF networks with the vertex hidden variable model and to provide the analytic expressions for the functions $`\rho (x)`$ and $`f(x,y)`$ that define SF networks in three special cases.

The fitness model can be easily generalized in order to have more than one fitness variable per vertex PScondmat . In the following, we consider a single real variable $`x`$ per vertex, with $`x0`$. As a probability distribution density function, $`\rho `$ satisfies $`\{\rho (x)0|_0^{\mathrm{}}\rho (z)𝑑z=1\}`$, while the linking probability $`0f(x,y)1`$. We define the primitive function of $`\rho (x)`$, the probability distribution $`R(x)=_0^x\rho (z)𝑑z`$. Indicating the number of vertices in the graph with $`N`$, one has the vertex degree

$$k(x)=N_0^{\mathrm{}}f(x,z)\rho (z)𝑑z.$$
(1)
Other quantities of interest are the average nearest neighbor connectivity (vertex degree correlation),

$$K_{\mathrm{nn}}(x)=\frac{N}{k(x)}_0^{\mathrm{}}f(x,z)k(z)\rho (z)𝑑z,$$
(2)
expressing the average degree of vertices that are nearest neighbors of vertices with fitness $`x`$, and the clustering coefficient (vertex transitivity),

$$C(x)=N^2\frac{_0^{\mathrm{}}_0^{\mathrm{}}f(x,y)f(y,z)f(z,x)\rho (y)\rho (z)𝑑y𝑑z}{k(x)^2},$$
(3)
that counts the fraction of nearest neighbors of vertices with fitness $`x`$ that are also nearest neighbors each other. Eqs. (1), (2), (3) are valid asymptotically when $`N`$ approaches infinity. Eqs. (2), (3) were first derived in Ref. BP03 , and expressed in a different form.
If $`k(x)`$ is an invertible and increasing function of $`x`$ then the probability distribution $`P(k)`$ is given by

$$P(k)=\rho (x(k))x^{}(k)$$
(4)
or, as a function of $`x`$,

$$P(k(x))=\frac{\rho (x)}{k^{}(x)}.$$
(5)
Since the degree probability is power-law distributed in most of the physical situations, we impose in Eq. (5) $`P(k)=ck^\alpha `$ with $`\alpha `$. The constant $`c`$ is fixed by the the normalization condition $`_{k_0}^NP(k)𝑑k=1`$:

$$c=\{\begin{array}{cc}\frac{\alpha +1}{N^{\alpha +1}k_0^{\alpha +1}}\hfill & \text{if }\alpha 1\hfill \\ \left(\mathrm{log}\frac{N}{k_0}\right)^1\hfill & \text{if }\alpha =1\hfill \end{array}$$
(6)
with $`k_0=\underset{x+0}{lim}k(x)`$. Note that $`k_0=\beta N`$ for some $`0<\beta <1`$, and $`cN^{(\alpha +1)}`$. Eq. (5) becomes:

$$ck^{}(x)(k(x))^\alpha =\rho (x).$$
(7)
By integrating Eq. (7) from $`0`$ to $`x`$ we get the following non linear integral equation:

$$k(x)=\{\begin{array}{cc}\left(k_0^{\alpha +1}+\frac{\alpha +1}{c}R(x)\right)^{\frac{1}{\alpha +1}}\hfill & \mathrm{if}\alpha 1\hfill \\ k_0e^{R(x)/c}\hfill & \mathrm{if}\alpha =1\hfill \end{array}$$
(8)
with $`k(x)`$ given by Eq. (1).

By multiplying both sides of Eq. (8) by $`\rho (x)`$ and integrating from $`0`$ to $`\mathrm{}`$ we get an analytic expression for the average vertex degree $`k`$. This expression can be used to express $`k_0`$ as a function of $`k`$, so that the final expressions do depend on the physical quantity $`k`$ only. For this purpose, the integral on the rhs is simply solved using the relation $`\rho (x)dx=dR(x)`$.

In the following we show an application of the model in three special cases of interest, comparing the analytic results with numerical simulations. It has to be noticed that once fixed $`N`$, in order to compute the quantities $`P(k)`$, $`K_{\mathrm{nn}}(k)`$, $`C(k)`$, from ensemble statistics, we need to perform two different average procedures. Firstly, we should extract a $`\{x_i\}_{i=1\mathrm{}N}`$ configuration with the distribution density $`\rho (x)`$ and keep it fixed, while creating ensemble elements using the linking probability $`f(x,y)`$ and averaging at the end. Secondly, we should repeat the above procedure a sufficient number of times. We assume that for large enough $`N`$ and ensemble elements, the procedure of first averaging with respect to the $`f`$ can be skipped.
Here we focus on two different problems: firstly, what we call direct problem, one assigns a distribution density function $`\rho (x)`$ and tries to find the linking probability function $`f(x,y)`$; secondly, what we call inverse problem, one assigns the linking probability function and tries to determine the fitness probability distribution density $`\rho (x)`$. The inverse problem is by far more complex and interesting than the direct one. For instance, in the case of protein SF network by assuming a reasonable linking function, we can retrieve the probability density distribution of fitness (e.g. some basic property of the macro-molecules).

We start with the special case of $`f(x,y)=g(x)h(y)`$ where both the direct and inverse problems can be analytically solved. Because of the symmetry of $`f(x,y)`$ with respect to its arguments one has $`g(x)h(x)`$, so that $`f(x,y)=g(x)g(y)`$. Eq. (1) becomes:

$$k(x)=Ng(x)_0^{\mathrm{}}g(z)\rho (z)𝑑z$$
(9)
that substituted into Eq. (8) gives equations in $`g`$ and $`\rho `$. If one fixes a given function $`\rho (x)`$, the equations in $`g(x)`$ can be easily solved. Take for instance the second equation corresponding to $`\alpha =1`$. One gets:

$$Ng(x)g=k_0e^{R(x)/c},g=_0^{\mathrm{}}g(z)\rho (z)𝑑z.$$

By multiplying left and right hand side by $`\rho (x)`$ and integrating from $`0`$ to $`\mathrm{}`$, considering that $`\rho (x)dx=dR(x)`$, we get:

$$g=\sqrt{k_0c(e^{1/c}1)/N}.$$

Finally, the solution reads:

$$g(x)=\sqrt{\frac{k_0}{Nc(e^{1/c}1)}}e^{R(x)/c}.$$
(10)
This procedure is applicable for any value of $`\alpha `$. Eq. (10) generates random networks with degree probability distribution $`P(k)1/k`$.
In order to test the result we take the choice reported in the caption of Fig. 1.

We conclude that for any given $`\rho (x)`$ there exists a function $`g(x)`$ such that the network generated by $`\rho (x)`$ and $`f(x,y)=g(x)g(y)`$ is scale free with arbitrary real exponent.

In this case both the average nearest neighbors connectivity and clustering coefficient are constant BP03 . Respectively:

$$K_{\mathrm{nn}}=Ng^2,C=\frac{g^2^2}{g^2},$$
(11)
as it can be derived from Eq. (2) and Eq. (3).

The inverse problem for $`f(x,y)=g(x)g(y)`$ is solved by substituting Eq. (9) into Eq. (7):

$$\rho (x)=cg^{}(x)g(x)^\alpha (Ng)^{\alpha +1}.$$

Let us remark that the assumptions on $`k(x)`$ forces $`g(x)`$ to be non-decreasing with $`g(\mathrm{})>g(0)>0`$.

In the case $`\alpha 1`$, the normalization condition $`R(\mathrm{})=1`$ results in

$$\rho (x)=\frac{\alpha +1}{g(\mathrm{})^{\alpha +1}g(0)^{\alpha +1}}g^{}(x)g(x)^\alpha .$$
(12)
By multiplying both sides of the previous equation by $`g(x)`$ and integrating from $`0`$ to $`\mathrm{}`$, one gets an expression for $`g`$ that does not explicitly contain the function $`\rho (x)`$. This expression can be used to get the allowed value of the constant $`c`$:

$$c=\frac{1}{N^{\alpha +1}}\left[\frac{g_{\mathrm{}}^{\alpha +1}g_0^{\alpha +1}}{\alpha +1}\right]^\alpha \left[\frac{\alpha +2}{g_{\mathrm{}}^{\alpha +2}g_0^{\alpha +2}}\right]^{\alpha +1}$$

if $`\alpha 2`$, with $`g_{\mathrm{}}=g(\mathrm{})`$ and $`g_0=g(0)`$. The particular cases $`\alpha =\{1,2\}`$ can be similarly treated.

The case $`f(x,y)=f(xy)`$ is more complicated. In this case both the nearest neighbor connectivity and clustering coefficient depend on the fitness $`x`$ and conversely on the degree $`k`$. We managed to solve this case in the particular case of an exponentially distributed fitness. We indicate with $`F(x)`$ the rhs of Eq. (8). Thus Eq. (8) becomes:

$$_0^{\mathrm{}}f(xu)\rho (u)𝑑x=F(x)/N.$$

By changing the integration variable into $`z=xu`$ we get:

$$_{\mathrm{}}^x\rho (xz)f(z)𝑑z=F(x)/N$$

that in the special case $`\rho (x)=e^x`$ becomes:

$$_{\mathrm{}}^xe^zf(z)𝑑z=e^xF(x)/N.$$

By differentiating with respect to the variable $`x`$ we finally obtain:

$$f(x,y)=\frac{F(xy)+F^{}(xy)}{N}.$$
(13)
In order to test the result we take the function and parameter choice of Fig. 2 caption.

The case $`f(x,y)=f(x+y)`$ is analogous.

Again, we consider the special case $`\rho (x)=e^x`$, getting now:

$$f(x,y)=\frac{F(x+y)F^{}(x+y)}{N}.$$
(14)
The solution of Eq. (8) for $`\alpha =2`$ obtained via Eq. (14) reads, recalling that $`k_0=\beta N`$ and using Eq. (6),

$$f(x,y)=\frac{1}{[1(\beta ^11)e^{(x+y)}]^2}.$$
(15)
Through Eq. (15) we clarify the assumption made in the original paper by Caldarelli et al. gcalda , where $`f(x,y)=\mathrm{\Theta }(x+yz)`$ with $`z=z(N)`$. Note that now with the latter choice of $`f(x,y)`$ one gets $`P(k)=Ne^zk^2`$ that forces $`z`$ to depend upon $`N`$ in order to get the correct normalization. The functional form of the $`z(N)`$ was numerically guessed by Ref. millozzi . To test the result we take the parameters reported in the caption of Fig. 3.

In these last two cases, both the nearest neighbor connectivity and clustering coefficient show non trivial $`k`$ dependence.

In conclusion we present a general procedure to reproduce real SF networks with arbitrary vertex degree distribution densities. More specifically, we found that, given a fitness distribution density $`\rho (x)`$, it is always possible to find a symmetric linking probability function $`f(x,y)`$ such that the resulting random network is scale-free with a given real exponent. We give the recipe to find these linking functions, in three cases of interest. In order to allow the generation of networks even closer to the real data, it would be desirable to have control not only on the vertex degree distribution, but also on the vertex transitivity and vertex degree correlation, by solving simultaneously Eqs. (2), (3), (5). As a first step, the compatibility of these three equations should be addressed, once the functions $`P(k)`$, $`K_{\mathrm{nn}}(k)`$, $`C(k)`$ are given. The solution of this problem is certainly very hard and is left open for the future. The relative ease with which we obtain SF structures seems to be the key ingredient in order to explain the ubiquitous presence and robustness of the real data.

We acknowledge R. Pastor-Satorras, P. De Los Rios, D. Garlaschelli, F. Squartini, S. Millozzi, S. Leonardi for useful discussions and support from EU FET Open Project IST-2001-33555 COSIN (www.cosin.org). VDPS and GC thank IMFM PAIS PA\_G02\_4 for support.


text:
# Young star clusters: Metallicity tracers in external galaxies

## 1. Models & tests

We use our evolutionary synthesis code galev to study the basic physical parameters of star clusters in interacting galaxies and violent starbursts. Star clusters can easily be approximated as simple stellar populations (SSPs), since all stars have the same age, metallicity and extinction. However, more complicated star formation / metallicity evolution histories can be studied by superimposing appropriate SSPs. The main ingredients of our code are: stellar isochrones (ensemble of stars with different masses at a given age) from the Padova group for metallicities in the range of -1.7 $``$ \[Fe/H\] $``$ +0.4, a stellar initial mass function (usually assumed to be Salpeter-like), the library of stellar spectra by Lejeune et al. (1997,1998), and gaseous emission (both lines and continuum). From the integrated spectra we derive magnitudes in a large number of filter systems. The gaseous line and continuum emission contributes significantly to the integrated light of stellar populations younger than $`3\times 10^7`$ yr both in terms of absolute magnitudes and derived colors (Anders & Fritze - v. Alvensleben 2003). The updated models are available from

http://www.uni-sw.gwdg.de/$``$galev/panders/

We developed a tool to compare our evolutionary synthesis models with observed cluster SEDs to determine the basic cluster parameters age, metallicity, internal extinction and mass independently. Using artificial clusters with various input parameters (with cluster SEDs taken directly from the evolutionary synthesis models) we systematically studied the impact of the choice of passbands, of finite observational magnitude uncertainties, and a priori assumptions. Two examples of the artificial cluster tests are shown in Fig. 1 and 2. Additional tests were performed using broad-band observations of star clusters in NGC 3310 (de Grijs et al. 2003a,b), confirming the results from the artificial cluster tests. Due to the young age of this cluster system these additional tests are restricted to ages younger than approx. 200 Myr. From these tests we conclude that:

1. At least 4 passbands are necessary to determine the 3 free parameters age, metallicity and extinction, and the mass by scaling the SED, independently.
2. The most important passbands are the U and B bands; for systems older than roughly 1 Gyr the V band is equally important.
3. NIR bands significantly improve the results by constraining the metallicity efficiently.
4. A wavelength coverage as long as possible is desirable. Best is UV through to NIR, thus tracing the pronounced kink/hook in the SEDs around the B band.
5. Large observational errors and/or wrong a priori assumptions may lead to completely wrong results.

## 2. The case of NGC 1569

As a first application the dwarf starburst galaxy NGC 1569 was chosen. Our sample enlarges the number of star clusters studied in this galaxy by a factor of 3. Our results for the derived physical parameters are in agreement with previous results of the starburst history in NGC 1569 in general, and of the two prominent “super star clusters” in particular, regarding age, mass and metallicity. In addition, we find a surprising change in the cluster mass function with age: The clusters formed during the onset of the burst (approx. 25 Myr ago) seem to exhibit an excess of massive clusters as compared to clusters formed more recently. Using various statistical methods we show the robustness of this result (Anders et al. 2003b).

## 3. References

Anders, P., Fritze - v. Alvensleben, U. 2003, A&A, 401, 1063
Anders, P., Bissantz, N., Fritze - v. Alvensleben, U., de Grijs, R. 2003a,
MNRAS, submitted
Anders, P., de Grijs, R., Fritze - v. Alvensleben, U., Bissantz, N. 2003b,
MNRAS, submitted
de Grijs, R., Fritze - v. Alvensleben, U., Anders, P., Gallagher, J. S., Bastian, N., Taylor, V. A., Windhorst, R. A. 2003a, MNRAS, 342, 259
de Grijs, R., Anders, P., Bastian, N., Lynds, R., Lamers, H.J.G.L.M., O’Neil, Jr., E.J. 2003b, MNRAS, 343, 1285
Lejeune, T., Cuisinier, F., Buser, R. 1997, A&AS, 125, 229
Lejeune, T., Cuisinier, F., Buser, R. 1998, A&AS, 130, 65


text:
# 1 Introduction

## 1 Introduction

In spite of the great success in the quantum field theory, there is still no consistent scenario in which the quantum gravity is included. The most promising scenario in this direction would be the string theory, where the consistency is ensured by the requirement of additional dimensions and supersymmetry. Supersymmetry is sometimes used to explain the large hierarchy between the scales of Grand Unified Theory(GUT) and the Standard Model(SM). Unlike the standard model, the Supersymmetric Standard Model(SSM) contains many flat directions in its scalar sector. The flat directions might appear from the moduli fields that parametrizes the compactified space of the string theory. Moreover, considering the brane extensions of the string theory, flat potential might appear from the potential for the distances between branes, or simply for the positions of branes in the compactified space.222Some cosmological implications and defect configurations of the flat directions are discussed in ref.. In the string theory, initially the sizes of extra dimensions had been assumed to be as small as the Planck mass, but later it has been observed that there is no reason to believe such a tiny compactification radius. The idea of the large extra dimension may solve the hierarchy problem. Denoting the volume of the $`n`$-dimensional compact space by $`V_n`$, the observed Planck mass is obtained by the relation $`M_p^2=M_{}^{n+2}V_n`$, where $`M_{}`$ denotes the fundamental scale of gravity. If one assumes more than two extra dimensions, $`M_{}`$ may be close to the TeV scale without conflicting any observable bounds. In this scenario the standard model fields are expected to be localized on a wall-like structure and the graviton propagates in the bulk. The most natural embedding of this picture in the string theory context is realized by a brane construction.333Constructing successful models for inflation with a low fundamental scale is still an interesting problem. Baryogenesis and inflation in models with a low fundamental scale are discussed in . Affleck-Dine baryogenesis in such scenarios are discussed in ref.. We think constructing models of particle cosmology with large extra dimensions is very important since we are expecting that future cosmological observations would determine the fundamental scale of the underlying theory.

On the other hand, we know historically that the characteristic features of the phenomonological models are revealed by discussing their cosmological evolutions. For example, so called the Polonyi problem for the flat directions of the supersymmetric models has been discussed by many authors. If the historical knowledge still applies to the above brane-extended flat directions, it should be important to seek other options, since the alternative approaches might reveal the underlying problems of the model or show us the novel approach to a new scenario of cosmology.

From the above point of view, we will consider Q-balls and show that inflation may start inside a Q-ball. It is well known that a variety of models possess flat directions that may lead to non-topological solitons of the Q-ball type. On the other hand, it is sometimes discussed that topological defects other than the Q-balls, which might be produced in the early stage of the Universe, would play important roles in particle cosmology.444Inflation from topological defect is discussed in ref. for conventional gravity and in ref. for models of the brane world. Constraints on hidden-sector walls from weak inflation are discussed in ref.. Here we briefly summarize the interesting differences between conventional defect inflation and the Q-ball inflation.

* In conventional models of defect inflation, the typical parameters such as the mass scales or the width of the defect are the constants that do not change with time. In this case, the conditions for defect inflation to start are determined solely by the form of the potential that induces defect configuration. On the other hand, the typical parameters of the Q-balls will be determined by their charges that may evolve with time by absorbing other Q-balls.
* Q-balls will be produced during the oscillation after inflation. If the Q-balls decay safely, the reheating may be dominated by the “surface reheating” that was discussed in ref.. On the other hand, if a produced Q-ball triggers inflation in its core, the second stage of inflation will start in the local area. Of course the second inflation might produce another Q-balls again, which might induce problematic eternal inflation.
* The Q-ball inflation might start as the secondary weak inflation that induces the temporal swmall expansion of a local area. Thus, unlike conventional defect inflation, the Q-ball inflation might be used to explain the bubble structure of the Universe.

In section 2 we show a simple example of the Q-ball that appears on a hybrid potential. Although there had been no explicit argument on such an extension of the flat potential of the Q-balls, it is easy to see that the extension of this type is quite natural in many phenomonological models. In section 3 we examine the conditions for the Q-ball inflation. In section 4, we apply the results obtained in section 3 to a specific model of surface reheating. We stress here that the ideas we will show in this paper are quite simple and generic.

## 2 Q-balls from hybrid potential

Here we consider a typical hybrid potential that has been used in models of D-term inflation. As usual, the trigger field $`\varphi `$ is assumed to have a steep potential, while the field $`\sigma `$ parameterizes a flat potential. The explicit form of the effective potential is

$$V(\sigma )=M_0^4\mathrm{log}\left(1+K_1\frac{|\sigma |^2}{M_1^2}\right)+m_{3/2}^2|\sigma |^2\left[1+K_2\mathrm{log}\left(\frac{|\sigma |^2}{M_2^2}\right)\right]$$
(2.1)
where $`m_{3/2}`$ is the gravitino mass, which is obtained by using the supersymmetry breaking scale $`\mathrm{\Lambda }_{SUSY}`$, as $`m_{3/2}=\frac{\mathrm{\Lambda }_{SUSY}^2}{M_p}`$. The constant $`K_1,K_2`$ represent the renormalization factors at one-loop, and $`M_1,M_2`$ are the renormalization scales. We introduce a dimensionless constant $`\eta M_0/\mathrm{\Lambda }_{SUSY}`$ for later convenience. The second term dominates when $`\sigma >\sigma _c`$, where $`\sigma _c`$ is defined as

$`\sigma _c`$ $``$ $`{\displaystyle \frac{M_0^2}{m_{3/2}}}`$ (2.2)
$`=`$ $`M_p\left({\displaystyle \frac{M_0}{\mathrm{\Lambda }_{SUSY}}}\right)^2`$
$`=`$ $`M_p\times \eta ^2,`$
where the numerical factor is neglected. For $`\eta 1`$, $`\sigma _c`$ becomes much larger than the Planck scale ($`\sigma _cM_p`$), and the second term in (2.1) is always negligible. On the other hand, for $`\eta 1`$, one should consider two types of Q-balls.

Let us first consider the case where $`\eta <1`$. As far as $`\sigma <\sigma _c`$, the potential is dominated by the first term in eq.(2.1) and the Q-balls will have the following well-known properties;

$`r_Q{\displaystyle \frac{Q^{1/4}}{M_0}},`$ $`\omega {\displaystyle \frac{M_0}{Q^{1/4}}}`$
$`\sigma M_0Q^{1/4}.`$ (2.3)
On the other hand, when $`\sigma >\sigma _c`$, the second term will dominate and the Q-balls will have the following properties,

$`r_Q{\displaystyle \frac{1}{\sqrt{|K_2|}m_{3/2}}},`$ $`\omega m_{3/2}`$
$`\sigma |K|^{3/4}m_{3/2}Q^{1/2}.`$ (2.4)
As these results are well established, we use (2) and (2) hereafter to examine the conditions for the Q-ball inflation. Because we are considering hybrid potential, we may assume $`M_0\mathrm{\Lambda }_{SUSY}`$ as well as $`M_0\mathrm{\Lambda }_{SUSY}`$. If one assumes $`M_0\mathrm{\Lambda }_{SUSY}`$, $`\eta `$ becomes much larger than $`O(1)`$. In this case, $`\sigma `$ is always much smaller than $`\sigma _c`$, which means that one may safely assume that (2) is always satisfied for any $`Q`$.

In the next section we use the above results and examine the conditions for the Q-ball inflation. We then analyze the properties of the Q-ball inflation.

## 3 Inflation from a Q-ball

As in the conventional models of topological inflation, the following conditions will be required.

* The radius of the Q-ball must become larger than the Hubble radius. Since the field $`\sigma `$ is trapped at a false vacuum within the Q-ball, inflation will start when the boundary of the Q-ball exits the horizon.
* If we consider the scenario where the charges of the Q-balls evolve with time, the expectation value of the $`\sigma `$ in the core will also change with time. In this case we should examine when the transition from $`\sigma <\sigma _c`$ to $`\sigma >\sigma _c`$ occurs. The properties of the Q-ball will be modified at $`\sigma _c`$. 555Of course one might consider the scenario where a huge Q-ball is produced at the earliest stage of the Universe and inflation starts abruptly.
* Naively, Q-balls might become a black hole before it induces inflation. This condition is rather trivial as we will show later.

Let us examine the above conditions in more detail. First we assume $`\eta 1`$, where both phases $`\sigma <\sigma _c`$ and $`\sigma >\sigma _c`$ might appear. When the charge $`Q`$ is small, the first term in (2.1) dominates and then (2) is reliable.666We are not compelling the scenario where small Q-balls evolves with time. Q-balls might be huge when it is produced before inflation. Since the vacuum energy $`\rho `$ inside Q-balls is as large as $`M_0^4`$, the Hubble constant when the boundary exits the horizon will be $`H\frac{M_0^2}{M_p}`$. Then the condition $`r_Q>H^1`$ is represented by

$$Q^{1/4}\frac{M_p}{M_0}.$$
(3.1)
The above result can be applied as far as $`\sigma <\sigma _c`$. However, the condition for $`\sigma <\sigma _c`$ is

$`Q^{1/4}`$ $`<`$ $`\left({\displaystyle \frac{M_0}{m_{3/2}}}\right)`$ (3.2)
$`=\left({\displaystyle \frac{M_p}{\mathrm{\Lambda }_{SUSY}}}\right)\left({\displaystyle \frac{M_0}{\mathrm{\Lambda }_{SUSY}}}\right)`$
$`={\displaystyle \frac{M_p}{M_0}}\eta ^2,`$
which is much smaller than (3.1). Thus in this case, we must conclude that evolving Q-balls will alter their properties from (2) to (2) before it induces inflation. Of course the evolving Q-balls must not decay into black holes before it induces inflation. Denoting the Schwarzschild radius by $`r_g`$, one can easily find the condition

$$Q^{1/4}\frac{M_p}{M_0}$$
(3.3)
which suggests that the critical charge for the black hole formation is much larger than the criteria (3.2).777To understand (3.3), we think it is helpful to consider a simplest toy model. Imagine a spherical region with the radius $`r_b`$. The vacuum energy is $`\rho =\rho _b>0`$ inside the “ball”, while $`\rho =0`$ outside. Then the mass of the “ball” is given by the formula $`M_b=\frac{4\pi }{3}r_b^3\rho _b`$. The Schwarzschild radius of the “ball” is $`r_g^2\rho _b/M_p^2`$. On the other hand, the Hubble parameter inside the “ball” when the boundary exits horizon is $`H^2\rho _b/M_p^2`$. Thus the naive black hole condition will give a rather trivial result.

Now the above results are suggesting that we should consider the Q-balls of $`\sigma >\sigma _c`$ in the case when $`\eta 1`$. When the Q-ball become large and the expectation value of the field $`\sigma `$ inside the Q-ball becomes larger than $`\sigma _c`$, the above criteria are represented by the following conditions. Here the vacuum energy density inside Q-balls is estimated by $`\rho m_{3/2}^2\left[(|K|^{3/4})m_{3/2}Q^{1/2}\right]^2`$. The radius of the Q-ball will exit the horizon when $`r_Q>H^1`$, which happens when

$`Q`$ $`>`$ $`|K_2|^{1/2}\left({\displaystyle \frac{M_p}{m_{3/2}}}\right)^2`$ (3.4)
$`=|K_2|^{1/2}\left({\displaystyle \frac{M_p}{\mathrm{\Lambda }_{SUSY}}}\right)^4`$
$`=|K_2|^{1/2}\left({\displaystyle \frac{M_p}{M_0}}\right)^4\times \eta ^4,`$
To verify our argument, we should examine if $`\sigma >\sigma _c`$ is satisfied when inflation starts. $`\sigma >\sigma _c`$ is satisfied if

$`Q`$ $`>`$ $`|K_2|^{3/2}\left({\displaystyle \frac{M_p}{m_{3/2}}}\right)^2\times \eta ^4`$ (3.5)
$`=|K_2|^{3/2}\left({\displaystyle \frac{M_p}{M_0}}\right)^4\times \eta ^8.`$
Because we are considering the case $`\eta 1`$, the condition (3.5) is safely satisfied when inflation starts.

Now we will examine the case $`\eta 1`$. As we are considering hybrid potential in this paper, $`\eta 1`$ is not unnatural. As we have discussed above, we should always use (2), because in this case $`\sigma _c`$ is much larger than $`M_p`$. The conditions for $`r_Q>H^1`$ and $`r_Q>r_g`$ are the same as the result obtained above. The properties of the Q-balls do not alter as the charge evolves, which is the only difference from the above result. In this case, $`\sigma \sigma _c`$ is always satisfied for any realistic value of $`Q`$. Thus our conclusion for $`\eta 1`$ is that the Q-ball inflation will start when $`Q(M_p/M_0)^4`$.

In the above discussions we have examined the conditions for inflation to start within Q-balls. Our second task is to examine the evolution of the field $`\sigma `$ during the Q-ball inflation. At the earliest stage of the Q-ball inflation, the field $`\sigma `$ is trapped at the false vacuum because of the large $`\omega `$. When inflation starts inside the Q-ball, the friction term dominates the equation of motion. Then $`\dot{\sigma }`$ decays as $`\dot{\sigma }e^{Ht}`$. Assuming that the change in $`|\sigma |`$ is much slower than that in $`\omega `$, one can obtain

$$\omega \omega _0e^{Ht},$$
(3.6)
where $`\omega _0`$ denotes the initial value of $`\omega `$ when inflation starts. In this case, one may expect two kinds of inflation that may start subsequently. The first inflation occurs during the period of $`\omega ^2>V^{\prime \prime }`$ when the field $`\sigma `$ is trapped at the false vacuum because of the large $`\omega `$. The e-foldings of the first inflation is

$$N_e\frac{1}{2}\mathrm{log}\left(\frac{w_0^2}{V^{\prime \prime }}\right),$$
(3.7)
where $`V^{\prime \prime }`$ is the effective mass of the field $`\sigma `$ at the false vacuum. Then after this period, the conventional slow-roll inflation will start if some conditions are satisfied. If $`\eta >1`$, the Hubble parameter might be smaller than the effective mass, where the situation is the same as the conventional D-term inflation. The e-foldings of the second inflation is determined by the charge of the Q-ball, which determines the initial value of $`\sigma `$. Even if the slow-roll inflation does not start, the expansion of the local area induced by the small inflation will affect the later cosmological structure formation. In any case, the cosmological observations of the present Universe might put some bound on the Q-ball inflation, which might put a bound on the phenomonological models.

Our last example is the Randall-Sundrum Type 2(RS-2) model. In this case the Friedman equation receives an additional term that is quadratic in the density. The Hubble parameter is related to the energy density by

$$H^2=\frac{8\pi }{3M_p^2}\rho \left(1+\frac{\rho }{2\lambda }\right),$$
(3.8)
where $`\lambda `$ is the brane tension. Denoting the ratio between $`\rho `$ and $`\lambda `$ by $`ϵ\lambda /\rho `$, the modified condition for $`H^1<r_Q`$ is relaxed when $`ϵ<1`$. Assuming that $`ϵ1`$, one can obtain for $`\sigma <\sigma _c`$,

$$Q\left(\frac{M_p}{M_0}\right)^4\times ϵ^2.$$
(3.9)
For $`\sigma >\sigma _c`$, it becomes

$$Q>|K_2|^{1/2}\left(\frac{M_p}{M_0}\right)^4\times \eta ^4ϵ.$$
(3.10)
From (3.9) and (3.10), one can see that the condition for Q-ball inflation is rather relaxed in the Randall-Sundrum Type 2 scenario.

## 4 Surface reheating and Q-ball inflation

In this section we will examine whether the Q-balls produced just after inflation leads to the surface reheating or to the Q-ball inflation. The surface reheating in models of running mass inflation is already examined in ref.. Here we mainly follow the setups of , and discuss if the Q-ball inflation takes place. The potential which can lead to a Q-ball formation is

$$V(\sigma )=m^2\left(1+|K|\mathrm{log}\left[\frac{\sigma ^2}{M_p^2}\right]\right)$$
(4.1)
We will assume that the maximum charge of a Q-ball might be as large as the maximum charge within the Hubble radius when Q-balls are produced. Denoting the Hubble parameter and the expectation value of the field $`\sigma `$ at the time of the Q-ball formation by $`H_Q`$ and $`\sigma _Q`$, one can obtain the maximum charge $`Q_{MAX}`$,

$$Q_{MAX}\frac{\omega _Q\sigma _Q^2}{H_Q^3}.$$
(4.2)
In generic situations, $`\omega _Q`$ is expected to be about the same order as the Hubble parameter $`H_Q`$. Finally, one can obtain the simple result,

$$Q_{MAX}\sqrt{K}\left(\frac{M_p}{m}\right)^2.$$
(4.3)
which is nearly the same order as (3.4).888Note that $`m`$ in (4.2) corresponds to $`m_{3/2}`$ in (3.4). However, it becomes much smaller than (3.4) if $`K1`$. Although it is quite difficult to calculate the exact values of the properties of the cosmological Q-balls, it seems fair to conclude from the above rough estimations that the surface reheating will be reliable if $`K`$ is much smaller than $`O(1)`$.999This result does not exclude the possibility that a few Q-balls might become abnormally large absorbing other Q-balls.

## 5 Conclusions and Discussions

In this paper we have examined the idea that a modified version of the conventional topological inflation might start within Q-balls. An extension of the Q-balls is discussed by using an explicit hybrid potential, which is useful for our discussion. As the flat direction might appear in the hybrid potential of the brane distance, Q-balls might appear for the brane rotation. We have obtained the criteria for the Q-ball inflation, which is comparable to the phenomenological values of the conventional Q-ball formation. Since the Q-ball formation is quite general in many models where flat directions are contained in the low-energy effective Lagrangian, we think the cosmological bound that will be obtained by considering Q-ball inflation is quite important.

## 6 Acknowledgment

We wish to thank K.Shima for encouragement, and our colleagues in Tokyo University for their kind hospitality.


text:
# Heteroepitaxial growth of high-K gate oxides on silicon: insights from first-principles calculations on Zr on Si(001)

## I Introduction

The scaling of the CMOS transistor has been the driving force behind the tremendous increase in microprocessor performance observed during the last decades. While the problems of the past were dominated by manufactural aspects, one now faces the first fundamental physical limitations as structures in logical devices approach atomic dimensions.

As a rule of thumb, the thickness of the gate-oxide has to be directly proportional to the channel length in a MOSFET (metal oxide silicon field effect transistor) device. In the course of the ongoing miniaturization, also the thickness of this insulating layer is being continually reduced. By the year 2007, gate-oxides in a transistor will approach a thickness of 1.5 nmroadmap , which corresponds to about ten typical atomic distances. The quantum mechanical tunneling currents through such a thin oxide are intolerable, and cause increased power consumption and deteriorated switching characteristics of the transistor. Replacing SiO2-based oxides, nowadays employed as gate dielectric, is one of the “key challenges” to the semiconductor industryroadmap , which has to find a solution within the next 4 to 5 years.

Employing high-k ($``$ large dielectric constant) oxides would allow for a greater physical thickness while preserving the electrical properties. Possible candidates have to meet an extensive list of requirementswil01  such as sufficiently large band offsets, thermodynamical stability or interface quality.

The Si–SiO2 system meets all these requirements in an unparallelized waybuc00 . First attempts to grow alternative oxides on Si(001) did not yield satisfactory results for a variety of reasons – often originating in the interface region.

In the course of changing to a new gate material, one also considers to now grow crystalline oxide layers. It has been a paradigm that the gate oxide must be amorphous in order to avoid dislocations and grain boundaries, which provide natural pathways for leakage currents and atom diffusion. However, in the very small devices used in the near future, the probability for such a defect may be sufficiently small. On the other hand, epitaxial oxides bear the promise of very low defect concentrations at the interface to silicon, which is crucial for efficient device operation.

The first account of heteroepitaxial growth of a transition metal oxide has been presented by a seminal work of McKeemck01  who succeeded in growing SrTiO3 on Si(001). For the divalent elements detailed studies have been performedyao99 ; oji01 ; her00 ; bak96 ; fan90 ; wan99 . As similar understanding does not yet exist for the technological relevant early transition metals such as Zirconium. Theoretical studies on Zr are limited bulk interface calculations of silicateskaw01 .

The deposition and formation of an interface between one of the major contenders for high-k oxides, namely Zirconia (ZrO2), and silicon has not yet been investigated with ab-initio simulations. As a first step we performed state-of-the-art electronic structure calculations of the deposition of Zr on Si(001) up to a coverage of two monolayers (ML).

## II Computational Aspects

We performed first-principles calculations within the framework of density functional theorydft1 ; dft2  using the gradient corrected density functional of Perdew, Burke and Enzerhofpbe . The electronic structure problem has been solved with the projector augmented waves methodblo94 , which uses augmented plane waves to describe the full wave-functions and densities without shape approximation. The PAW method as implemented in the CP-PAW code employs the Car-Parrinello approachcpmd  to minimize the total energy functional.

The core electrons are described within the frozen core approximation with the semi-core (4$`s`$ and 4$`p`$) shells of the Zr atoms treated as valence electrons. Plane wave cutoffs of 30 and 60 Ry for the wave functions and the density have been used.

For metallic system we minimized the Mermin functional with respect to the occupation numbers which yields the Fermi-Dirac distribution for the electrons. An electronic temperature of 1000 K was used. The zero-Kelvin result has been extrapolated using the method suggested by sch97 .

The reconstructed Si(001) surface is the template for our growth studies. It is modeled by a slab consisting of 5 silicon layers, where the bottom layer is saturated with two hydrogen atoms per silicon. The position of the hydrogen atoms and the lowest silicon layer has been frozen. Supercells with 16 atoms per layer have been used for the calculations of Zr in the dilute limit, 8 atoms per layer for all higher coverages. In the first case, the distance between two periodic images of an adatom is 15.36 Å. For all surface calculations the k-mesh density corresponds to 64 lateral k-points per $`p(1\times 1)`$ surface unit cell.

## III Dilute Limit

Fig. 1. shows the structure of the reconstructed Si(001) surface. Excellent accounts covering the full complexity of the surface reconstruction can be found in literaturecha79 ; sch01 . The resulting structure can, however, be explained by the following considerations: While every bulk silicon is 4-fold coordinated, the surface atoms lack their upper bonding partners which leaves 2 dangling bonds per atom, each occupied by one electron. In a first step, two neighboring silicon atoms form a dimer bond – initially parallel to the surface – which leaves one dangling bond per atom. Now both electrons are transferred into one dangling bond to reduce the number of unpaired electrons to zero. This results in the buckling as the silicon with the filled dangling bond prefers a tetrahedral bonding arrangement, while the other prefers a planar sp2 configuration. From the electrostatic point of view, the alternating buckling behavior is favorable. Our calculations readily reproduce previously reported results such as the $`c(4\times 2)`$ reconstruction and the difference in $`z`$-coordinate of two silicon atoms within a dimernor93 .

The energy surface of a Zr adatom on Si(001), has been obtained by freezing the lateral position of the adsorbed atom relative to the slab backplane, while all other degrees of freedom were fully relaxed. We used a grid of twelve Zr positions in the irreducible zone of the reconstructed $`p(2\times 1)`$ silicon surface.

The resulting total energy surface for an isolated Zr on top of Si(001) is shown in Fig. 2. Two nearly degenerate positions can be identified. One is located in the valley right in the middle of two dimers of neighboring rows and the other on top of a dimer row between to adjacent dimers. A third, local minimum, is located in the valley and has an energy 0.30 eV higher relative to the most favorable positions. The diffusion is quasi one-dimensional with barriers of 0.70 eV for diffusion parallel to the dimer rows and 1.63 eV from the valley to the row.

## IV Coverage of 0.25 and 0.5 Monolayers

In order to investigate the formation of a continuous film, we increased the coverage. The structures with coverages of 0.25 ML and 0.5 ML have similar energies per Zr atom as the dilute limit (see Fig. 3). A wealth of complex structures has been found. Here, however, we only summarize the main trends of the chemical binding.

To first approximation Zr prefers a formal 4+ charge state – the projected density of states of the Zr $`d`$-states is located well above the Fermi level. It is well known that for transition metal cations the $`s`$-and $`p`$-electrons are located above the $`d`$-states, and have only a minor effect on the occupied states. The silicon dimers at the surface accept up to two electrons in their dangling bonds. A clear structural indication that the dangling bonds are filled is the disappearance of the dimer buckling.

Further electrons that are supplied at increased coverage of Zr, occupy the antibonding states of the dimer bond. As a consequence the silicon dimers break up. This happens at a coverage of half a ML and above. It should be noted, however, that metallic states accept some of the electrons so that the number of broken dimers does not directly correlate with the number of Zr adatoms in the ratio one to one. As a result of the interplay between breaking dimers and metallic states we find fairly complex reconstructions of the surface structure for intermediate coverages.

## V Silicide Formation

Upon increase of the coverage to a full ML, the energy drops by 0.64 eV per adatom (see Fig. 3). The stable structure with ML coverage is shown in Fig. 4. All dimer bonds of the surface layer are broken, and the Zr atoms occupy the centers of the square array of surface silicon atoms in the resulting $`p(1\times 1)`$ reconstruction. Such a layer is one structural element of bulk ZrSi2.

The energy gain is, however, not due to the surface geometry but can be attributed to a dimer reconstruction of the silicon subsurface. A metastable state without this reconstruction is even higher in energy than adsorbed Zr at lower coverages.

During the reconstruction, the silicon atoms in the layer underneath the ZrSi surface layer form dimers analogous to the bare silicon surface. In contrast to the dimer row reconstruction of the silicon surface, these subsurface dimers are not buckled and are arranged in a checkerboard instead of a row pattern.

At a coverage of two MLs a second ZrSi layer is formed. This configuration is shown in Fig. 5. The second layer is nearly identical to the first, but is shifted laterally. The double layer is again a structural element of bulk ZrSi2, the structure of which is shown in Fig. 6. The ZrSi double layer consists of one-dimensional Si zig-zag chains separated by Zr atoms. The Zr atoms lie approximately in the plane of the upper and lower atoms of the Si chain. As for ML coverage the first silicon layer underneath the ZrSi double layer exhibits a dimer reconstruction.

Despite an additional energy gain from the ML coverage to a coverage of two MLs of 0.25 eV per adatom, even this structure is still 0.67 eV higher in energy than bulk ZrSi2. While the atomic process has not yet been resolved in every detail, our findings give strong indications for the nucleation of the silicide.

## VI Discussion and Conclusion

We presented the results of state of the art ab-inito electronic structure calculations aiming at understanding the deposition of Zr atoms on a Si(001) surface as it is the case in an MBE reaction chamber.

Our results are summarized in Fig. 3, which shows the energy per Zr adatom as a function of the coverage. For coverages below one monolayer the the energy is nearly independent of the coverage. At a coverage of 1 ML we observe a sharp drop in energy by 0.64 eV, followed by further drops in energy for higher coverages. All structures are less stable than bulk silicide.

Our findings suggest that islands with a local coverage of 1 ML or higher are formed even at low coverages. The islands contain structural elements of bulk ZrSi2, which is more stable than any surface structure. Therefore a likely scenario is the formation of bulk silicide grains, that disrupt the surface morphology, and are detrimental for epitaxial growth. Modification of growth conditions, such as exposing the surface to a oxygen containing ambient, may bypass silicide formation during the first growth steps.

## Acknowledgments

This work has been funded by the European Commission in the project ”INVEST” (Integration of Very High-k Dielectrics with CMOS Technology) IST-2000-28495 and by the AURORA project (SBF F011) of the Austrian Science Fond. This work has benefited from the collaborations within the ESF Programme on ’Electronic Structure Calculations for Elucidating the Complex Atomistic Behavior of Solids and Surfaces’.


text:
# Generalized Verlinde formulae for some Riemann surface automorphisms

## Acknowledgements

I would like to thank P. Bantay for helpful discussions.


text:
# Quark mass dependence of pseudo-scalar masses and coupling constants

## 1 Introduction

Monte Carlo simulations of QCD with dynamical quarks are done in most cases at relatively large quark masses (typically two quark flavours with $`m_{ud}>m_s/2`$). This makes the extrapolation to the physical point $`m_{ud}m_s/20`$ rather uncertain. The extrapolation is done by using (PQ)ChPT – typically to NLO (1-loop) order. Estimates  show that one should perform simulations in the range $`m_{ud}<m_s/4,m_s/3`$ in order to see the expected logarithmic dependence. Matching the predicted functional dependence is a crucial test for Lattice QCD.

Until recently the comparison between lattice data and ChPT was not satisfactory. In a recent paper  we suggested that the agreement can be found when going to light enough dynamical quarks. Here we supplement some integration of the analysis of  and provide some further comments.

## 2 Simulations

For our simulation we used the algorithm described in  and further improved as in  and references therein. In our range of parameters we found that the cost for producing one independent gauge field configuration roughly goes as:

$$\text{Cost}F(am_q)^2\mathrm{\Omega },$$

where $`am_q`$ is the quark mass in lattice unit, $`\mathrm{\Omega }`$ is the number of lattice points and the factor $`F10^7`$flop for the plaquette, but one order of magnitude smaller for $`m_\pi `$ and $`f_\pi `$.

We produced three sets of $`O(1000)`$ thermalized configurations for $`N_f=2`$ unimproved Wilson fermions at Vol$`=16^4`$, $`\beta =5.1`$ and $`\kappa =0.176,0.1765,0.177`$. From these points we extrapolated a value of $`r_0/a`$ at $`\kappa _{cr}=0.1773(1)`$ which is $`r_0(\kappa _{cr})/a=2.65(7)`$. This corresponds to an UV cutoff $`a=0.189(5)`$fm $`(1.04\text{ GeV})^1`$, and to a physical volume $`L3`$fm. In those three points we found a ratio $`m_q/m_s`$ equal to $`0.68`$, $`0.59`$ and $`0.35`$ respectively. Here $`m_q`$ represents the sea quark mass defined as $`M_r=(m_\pi r_0)^2`$, when the strange quark mass corresponds to $`M_r=3.1`$ . For the lightest valence quarks the ratios $`m_q^{val}/m_s`$ become respectively $`0.37`$, $`0.34`$ and $`0.17`$. Even with these small quark masses finite volume effects are expected to be under control since we have always $`Lm_\pi ^{val}4.8`$. Of course this is payed with a very low UV cutoff, and one expects large lattice artifacts. These are taken into account in the analysis.

## 3 Comparison with ChPT

Our first goal is to confront the results of numerical simulations with the (PQ)ChPT formulas . In order to cancel the $`Z`$-factors of multiplicative renormalization, which in the case of a mass-independent renormalization scheme only depend on the gauge coupling and not on the quark mass, we considered ratios of quark masses ($`m_{qA}`$), pion masses ($`m_{AB}`$) and pion decay constants ($`f_{AB}`$). Here $`A,B`$ stand for the flavor indices of valence ($`V`$) or sea ($`S`$) quarks.

If we assume that there are no lattice artifacts, no NNLO corrections and we take for $`\chi _S=\frac{2B_0m_{qS}}{f_0^2}`$ its tree level estimate, $`\chi _S^{est}\frac{(r_0m_\pi )^2}{(r_0f_0)^2}`$, then the ratios (here $`\xi \frac{\chi _V}{\chi _S}`$)

$`RRf`$ $``$ $`{\displaystyle \frac{f_{VS}^2}{f_{VV}f_{SS}}}=1+{\displaystyle \frac{\chi _S^{est}}{32N_s\pi ^2}}[\xi 1\mathrm{log}(\xi )],`$
$`RRn`$ $``$ $`{\displaystyle \frac{4\xi m_{VS}^4}{(\xi +1)^2m_{VV}^2m_{SS}^2}}=1+{\displaystyle \frac{\chi _S^{est}[\mathrm{log}(\xi )\xi +1]}{16N_s\pi ^2}}`$
are non trivial and parameterless predictions of ChPT. This provides a strong check of how far we are from the NLO ChPT regime. In figure 1 and 2 we superimpose the predicted functions to the data. The agreement is of course not complete (there are indeed $`O(a)`$ and NNLO effects), but the corrections are sub-dominant contributions. Moreover the agreement improves when the sea quark masses decrease.

Encouraged by this results, we systematically compared our data with those ratios of pion masses and coupling constants which determine the Gasser-Leutwyler  coefficients $`L_4`$, $`L_5`$, $`L_6`$ and $`L_8`$. Since we expected large lattice artifacts in the data, the comparison was done with W(PQ)ChPT , including $`O(a)`$ lattice artifacts in the effective continuum theory. Besides that, we also included the relevant contributions from NNLO ChPT . Although this involve many parameters, one can obtain enough constraints from Partially Quenched simulations. The details of the fitting procedure and the results are described in . Here we simply add that, following the analysis of , we can now include also $`O(a^2)`$ lattice artifacts in our fit of the pion mass ratios. It turns out that this does not add new parameters to the fit, but it amounts to a redefinition of the Wilson-ChPT coefficients $`W_i`$.

In order to have a consistency check of the surprisingly small lattice artifacts that we found, we compared the results of  with those obtained in an older simulation with larger lattice spacing ($`\beta =4.68`$, $`a0.27`$fm) . The comparison (fig. 3 and 4) shows quite small scale breaking. Larger scale breaking ($`8\%`$) are observed for the ratio $`m_\pi /m_\rho `$, which at fixed $`M_r=1.088`$ goes from $`0.613`$ (at $`\beta =5.1`$) to $`0.664`$ (at $`\beta =4.68`$).

To summarize: we showed that it is possible to simulate, with reasonable costs, very light dynamical quarks. Compensating $`O(a)`$ effects by introducing $`O(a)`$ terms in the PQCh-Lagrangian is, in this case, a viable alternative to $`O(a)`$ improvement of the action. The observed $`O(a)`$ contributions in the light Goldstone boson sector are surprisingly small, while NNLO are still important. The expected behavior predicted by PQChPT is already visible, although a quantitative determination of the LEC’s still needs further simulations at smaller masses and lattice spacing. Most of the numerical calculations presented here have been done at the computers of NIC-Juelich and Zeuthen.


text:
# Auger Electron Cascades in Water and Ice

## 1 Introduction

High energy X-rays damage matter, mainly through the photoelectric effect. A photon is absorbed by an atom and a photoelectron is emitted from the inner (K) shell of the atom. In biological molecules, the remaining hollow ion might relax through Auger decay, in which an electron from a higher level falls down and a further electron is emitted with an energy between 250 eV and 2 keV . The basic atomic physics for these processes is well understood and cross sections for photoionization have been tabulated . Furthermore, Auger line widths have been measured giving K-hole life times, of e.g. 11.1 fs (C), 9.3 fs (N), 6.6 fs (O) and 1.3 fs (S) . During the photoionization process, the photoelectron may interact with electrons in the valence shell. For elements of biological significance this may lead to the so called shake-up and shake-off effects  in which a further low energy electron (10-100 eV) is emitted. Semi-empirical quantum calculations were used to estimate that such shake-off ionizations occur in 10 - 30% events .

The advent of X-ray free electron lasers (FEL) will enable a whole range of new experiments in physics, chemistry and biology within a few years . In the category of biological applications, we envision performing structural studies on large biomolecules, biomolecular aggregates or nanocrystals. Molecular dynamics simulations of a protein molecule in a FEL beam are encouraging as to the feasibility of such experiments . Preliminary studies with a simple hydrodynamic model have suggested that  a layer of water or ice may delay the eventual Coulomb explosion of a biomolecule, such as a protein, in a FEL beam . The reason for this is that water molecules at the surface supply electrons to neutralize the charging core of the droplet, where the protein is expected to be situated, while protons, carrying the positive charge, are expelled from the water shell. This may be advantageous in future experiments if the sample is held in water droplets.

In the aftermath of a photoionization event in a liquid or solid, further electrons will be liberated through collisions between the Auger electron and the atoms, leading to a cascade of ionizations. Recent studies  have made a critical assessment of the total amount of electrons due to an Auger electron of 250 eV (carbon)111Although the actual Auger energy for carbon in diamond is 265 eV , we have used the value of 250 eV for comparison with earlier work.. It was predicted that in diamond, between 6 and 13 electrons will be released, while for amorphous carbon the number is somewhat higher due to the absence of a band gap . Here we present results on water, both in liquid and solid state. The understanding of Auger electron cascades in water is essential for modeling the shielding effects water molecules will have on e.g. a protein sample. Furthermore, water is less dense than diamond and therefore more comparable to biomolecules. The Auger electron energy for oxygen in a water molecule is 507.9 eV , providing twice as much energy to the electron than carbon.

## 2 Theory

To calculate the cross sections for scattering of electrons on atoms used in the molecular dynamics model, we have used a solid state approach, assuming a neutral crystal, or, in the case of liquid water, a neutral cluster with periodic boundary conditions. The same approach has been used earlier by Ziaja et al. . Under this assumption, there are three different interactions that a free electron in the medium can undergo: (i) elastic collision with atoms, (ii) inelastic collision with atoms or (iii) recombination with atoms or ions. The probability of electron-ion recombination is low on the time scale considered here ($`t100`$ fs) and therefore neglected . Inelastic collision of electrons with ions are not considered, since the number of ionizations caused by a photoelectron or single primary Auger electron is much smaller than the total number of atoms in the sample. Scattering of electrons on phonons and the vibrational cross section are not considered here, though they become relevant for low energies ($`E10`$ eV) .

### Elastic scattering

The cross section for the elastic scattering of electrons on atoms is calculated using programs from the Barbieri/van Hove Phase Shift package . The calculations are based on the partial wave expansion technique, using the muffin-tin potential approximation. The expression for the elastic cross section as a function of the electron energy $`E`$, is

$$\sigma _{\mathrm{elastic}}(E)=\frac{2\pi \mathrm{}^2}{Em_\mathrm{e}}\underset{l}{}(2l+1)\mathrm{sin}^2\delta _l,$$
(1)
where $`\delta _l`$ is the phase shift of each partial wave $`l`$. For further details on the elastic cross section calculations see reference  and references therein.

The energy of an elastically scattered electron is conserved, and the only effect the collision has on the system is that the direction of the electron velocity changes. For the calculations of the cross sections in liquid water we have used a unit cell with hydrogen bonded water molecules at liquid density (0.9965 $`g/cm^3`$), while for ice calculations we have used hexagonal ice, $`I_\mathrm{h}`$ (density 0.9197 $`g/cm^3`$). The results of our calculations are shown in Fig. 1 together with results for diamond . The results compare very well with the experimental data for electron scattering on water molecules.

### Inelastic scattering

A correct treatment of the inelastic scattering is more complicated and, to our knowledge, there is no complete method to describe this. We followed the approach presented in our previous work . We have used two different optical data models, the Ashley model  and the Tanuma, Powell and Penn model (TPP2) , together with the Lindhard dielectric function $`\epsilon (q,\omega )`$  for the calculations of the inelastic cross section

$$\sigma _{\mathrm{inelastic}}(E)\frac{1}{E}_0^{\mathrm{}}d\omega _q_{}^{q_+}\frac{\mathrm{d}q}{q}\mathrm{Im}\left(\frac{1}{\epsilon (q,\omega )}\right),$$
(2)
were $`\mathrm{}q_\pm =\sqrt{2m_\mathrm{e}}\left(\sqrt{E}\pm \sqrt{E\mathrm{}\omega }\right)`$. The above cross section depends on the knowledge of the dielectric function, which is only known for photons, i.e. for momentum transfer equal to zero. The two different models, Ashley and TPP2, present different ways to empirically compensate for that. The Ashley approach includes the exchange between the incoming electron and the electrons in the crystal, in accordance with the non relativistic Møller cross section . The TPP2 model is developed for calculating the differential inelastic mean free path of electrons in a solid. A more thorough discussion of the two models can be found in ref. . From the absorption coefficient for liquid water , we deduced the dielectric function $`\epsilon (0,\omega )`$ for zero momentum transfer. For ice we used experimental data  for energies up to 30 eV, and for higher energies we have extrapolated the values assuming the shape of the curve to be the same as for liquid water (Fig. 2).

An electron that scatters inelastically will change its path and loose energy to the lattice and, if the lost energy is high enough to excite a bound electron from the valence band to the conduction band, to the excited electron. For simplicity we assumed that the energy lost to the lattice is absorbed into the system and in the present study we have not included any treatment for holes. The band gaps of water and ice we have used for the calculations were $`E_{\mathrm{gap}}^{\mathrm{water}}=8.7`$ eV  and $`E_{\mathrm{gap}}^{\mathrm{ice}}=7.8`$ eV . For water the Fermi level $`E_\mathrm{F}^{\mathrm{water}}=14.85`$ eV is taken from experimental data , and for ice $`E_\mathrm{F}^{\mathrm{ice}}=11.74`$ eV was calculated using the free electron approximation with the density taken from reference . The calculations were performed at temperature $`T_{\mathrm{ice}}=270`$ K and $`T_{\mathrm{water}}=300`$ K. The inelastic cross sections for ice and water are shown in Fig. 3 together with those from diamond, taken from reference .

From the wealth of experimental data for electron impact ionization of atoms and molecules, models and recommended cross sections have emerged for different targets. We have compared the theoretical estimates for inelastic scattering of electrons in solids (diamond, ice), with the recommended cross sections for carbon  and a successful model that describes ionization of water molecules  (Fig. 4). The former represent a fit to the atomic data for several light atoms and ions , while the latter, known as the Binary-Encounter-Bethe (BEB) model, provides analytical formulae for the ionization cross section of atomic and molecular orbitals by summing up all contributions. The BEB model was found to reproduce experimental data for small atoms and a variety of molecules very well .

We find that, in the case of diamond, the inelastic cross section is lower than the ionization of single atoms, reflecting the nonlocalized nature of the electrons in the solid. This difference could prove significant for the estimation of electron interaction in solids, if the solid loses its crystalline structure (due to Coulomb explosion while exposed to a FEL beam). For ice and water, however, we find a smaller difference between the ionization cross sections  and our theoretical estimates (Fig. 4), and with the opposite sign.

## 3 Methods

Classical simulations of the electron trajectories were performed, where Newton’s equations of motion were solved using a Leap-Frog integration scheme  with a time step of $`\mathrm{\Delta }t=1`$ attoseconds. The lattice, be it diamond, water or ice, was modeled implicitly through the cross sections for collisions. Elastic as well as inelastic collisions were modeled as stochastic processes, where the probability $`P(E)`$ for an event was determined from:

$$P(E)=\frac{\rho N_{\mathrm{av}}}{m}\sigma (E)|v_\mathrm{e}(E)|\mathrm{\Delta }t,$$
(3)
where $`\rho `$ is the density of the material, $`N_{\mathrm{av}}`$ is Avogadro’s number, $`m`$ is the mass of the atom (C in diamond), respectively water molecule, and $`v_\mathrm{e}`$ is the velocity of the electron. $`\mathrm{\Delta }t`$ was chosen such that $`P1`$ in all cases. At each time step in the simulation the following processes occur:
(1) An electron with a certain kinetic energy $`E_\mathrm{k}`$ travels through the sample. The scattering probability (both inelastic and elastic) $`P(E_\mathrm{k})`$ is compared to a random number $`0r<1`$. If no scattering occurs, the velocity is kept unchanged.
(2a) If the electron scatters elastically it changes its momentum according to a probability function derived from the differential cross section (the formulae in  were applied for water and ice). No momentum transfer to the lattice has been considered at this point.
(2b) In the case of inelastic scattering, the electron will lose energy $`E_{\mathrm{lost}}`$ and change its path. If $`E_{\mathrm{lost}}>E_{\mathrm{gap}}`$ an additional electron is liberated from the valence band with $`E_\mathrm{k}=E_{\mathrm{lost}}E_{\mathrm{gap}}`$. If $`E_{\mathrm{lost}}<E_{\mathrm{gap}}`$, the energy is transfered to the lattice. In both cases the momentum is changed based on the differential cross section calculated for water and ice using formulae from .
(3) The electron position is integrated.
The spatial electron dynamics program is part of the GROMACS software package .

## 4 Results and Discussion

The classical dynamics simulations of the electrons produced the same results for diamond as our earlier Monte Carlo simulations , showing that the implementation of the algorithm is correct, and that the integration time step was small enough. In principle, the model could be enhanced by adding electrostatic interactions between electrons and ions . It is however not trivial to model the interactions between electrons and ions classically without severe approximations. Since we have only simulated single Auger electron trajectories, and the cross sections were determined for neutral systems, we have deemed it unnecessary to add the Coulomb term. Finally, we note that correct treatment of electron/ion or electron/hole systems can be performed using e.g. time-dependent density functional calculations, but it is complicated, and not computationally feasible for macroscopic systems .

### Ionization rates

We find that an Auger electron of 500 eV will produce typically between 5 and 10 secondary electrons in the first femtosecond in ice and water, a result which is comparable with the average number of ionization events in diamond from a 250 eV Auger electron. In ice and water, however, saturation appears much later, after about 100 fs, when the average number of ionizations is around 20 for liquid water and 25 for ice. These results are shown in Fig. 5, which not only compares the ionization rates for ice, water and diamond, but also shows that the simulation of diamond using classical dynamics reproduces previous results obtained with a Monte Carlo method well . We note that, although the cross sections for ice and water are rather similar, the difference in the number of ionizations in the Auger electron cascade is determined mainly by the the approximately 1 eV difference in the band gap between liquid water and ice.

### Energy distribution

The total energy of the electrons shows a substantial drop during the first 10 femtoseconds, with the rest of the energy being transfered to the lattice (Fig. 6). The electron gas is cooling down very rapidly, due to the large number of ionizations, and the average temperature of the electrons $`kT_{\mathrm{el}}`$ is similar for ice and water. The energy lost to the lattice is more than 50% of the initial Auger electron energy. For comparison, the results obtained for diamond are also shown, consistent with the previous calculations .

The electron energy distribution as a function of time in the cascade (Fig. 7) is particularly interesting. In the early stages of the cascade (of the order of 0.1 fs), the electron population is divided in two sets, initial Auger electrons which are highly energetic (450-500 eV) and the first secondary electrons of lower energy. At the critical time of about 1 fs, the population of high energy electrons fades away, leaving a spread of electrons with intermediate and low energies. Subsequent collisions will cool down the electron gas very fast, and after 10 fs all the electrons will have energies less than 20 eV. At such low energies (below 10 eV) and long time scales (100-1000 fs), new processes become relevant, such as electron-ion recombination or electron-phonon scattering, which are not yet included in our model.

### Spatial distribution

We have performed an analysis of spatial distribution of electrons in the cascade in our classical dynamics model for ice, water and diamond. Fig. 8 shows the average gyration radius of the electron cloud for ice, water and diamond. The spatial extension of the cloud in ice and water is three times larger than the one in diamond, due to the considerably lower density and the corresponding longer mean free path an electrons in water/ice.

## 5 Conclusion

An understanding of the behaviour of water and ice exposed to an intense X-ray beam is important for planned biology experiments with FELs. Biological macromolecules carry about 30-60% water with them as structural water and water bound as a surface layer. In the light of results from a hydrodynamic model , it was suggested that embedding a sample in water or ice may be beneficial in the sense that the eventual Coulomb explosion  could be delayed by the presence of the water or ice layer. Our results show that in hexagonal ice, $`I_\mathrm{h}`$, a slightly higher number of secondary electrons are liberated as compared to water (Fig. 5), despite the slightly lower temperature. Furthermore, the electrons are spread over a larger volume (Fig. 8). Although the differences are small, it seems therefore that water is the preferred medium for encapsulating biomolecules when studied in FEL beams. If the sample is embedded in vitrious ice, which could be true in case of a rapid freezing, we expect the behaviour to be unaffected compared to water.

### Acknowledgments

We would like to thank Beata Ziaja, Magnus Bergh, Jorge Llano, Lars Nordström, Richard London, Gyula Faigel and Zoltán Jurek for stimulating discussions. The Swedish research foundation is acknowledged for financial support.

## Figure captions

1. Elastic cross sections for scattering of electrons on ice, water and diamond. Lines represent the scattering cross section per scattering center, derived using the Barbieri/van Hove Phase Shift method , for ice (full line) and water (dashed line). For comparison, the cross section for diamond taken from  is shown (dotted line). The experimental points represent the elastic scattering of electrons on water molecules as measured in .
2. The imaginary part of the inverse dielectric function, $`\mathrm{Im}(1/\epsilon (0,E))`$, for water and ice plotted as a function of the incoming photon energy.
3. Inelastic cross section for electron scattering on ice and water as a function of the electron energy. Lines represent the cross section per scattering center, calculated using two optical data models, the Ashley model  and the Tanuma, Powell and Penn model (TPP2) . For comparison, the cross sections for diamond taken from  are also shown.
4. Comparison between the inelastic cross section for electron scattering on solids (ice and diamond) and the ionization cross section for electron scattering on single atoms and molecules. Lines represent the calculated cross section using TPP2 model  for diamond (solid line) and ice (dashed line). Points represent the recommended cross section for carbon atoms  (solid points), and the calculated cross section Binary-Encounter-Bethe (BEB) model  for ionization of water molecules (circles).
5. Average number of secondary electrons as a function of time, from an initial Auger electron (at time t=0 fs) in ice (solid line), water (dashed line) and diamond (dotted line). Lines correspond to calculations using Ashley’s model .
6. Average energy of the free electrons in an electron cascade coming from an initial Auger electron in ice and water as a function of time (thin lines). The rest of the Auger electron energy is transfered to the lattice. Average temperature $`kT_{\mathrm{el}}`$ of the electron gas versus time (thick lines) for water and diamond.
7. Histogram of the energy distribution $`N(E,t)/N(t)`$ among electrons at different time steps in ice. At a fixed time, the histogram is normalized to the number of electrons in the cascade $`N(t)`$. Results are shown for an average cascade produced by an Auger electron in ice.
8. The gyration radius (average distance of electrons from the cloud center of mass) of the electron cascade as a function of time, for ice, water and diamond.