text:
# 1 Introduction

## 1 Introduction

The principles of equilibrium statistical mechanics are well established. But, thermal equilibrium is a special case, and little is known about the properties of systems not in equilibrium, for example about the relaxation toward the stationary state. Some interesting problems in non-equilibrium systems are non-equilibrium phase transitions described by phenomenological rate equations, and the way the system relaxes to its steady state. As mean-field techniques, generally, do not give correct results for low-dimensional systems, people are motivated to study exactly-solvable stochastic models in low dimensions. Moreover, solving one-dimensional systems should in principle be easier. Exact results for some models on a one-dimensional lattice have been obtained, for example in \[1–14\]. Different methods have been used to study these models, including analytical and asymptotic methods, mean field methods, and large-scale numerical methods.

The Glauber dynamics was originally proposed to study the relaxation of the Ising model near equilibrium states. It was also shown that, there is a relation between the kinetic Ising model at zero temperature and the diffusion annihilation model in one dimension. There is an equivalence between domain walls in the Ising model and particles in the diffusion annihilation model. Kinetic generalizations of the Ising model, for example the Glauber model or the Kawasaki model, are phenomenological models and have been studied extensively \[15–20 \]. Combination of the Glauber and the Kawasaki dynamics has been also considered \[21–23\].

In , an asymmetric generalization of the zero-temperature Glauber model on a lattice with boundaries was introduced. It was shown there that, in the thermodynamic limit, when the lattice becomes infinite, the system shows two kinds of phase transitions. One of these is a static phase transition, the other a dynamic one. The static phase transition is controlled by the reaction rates, and is a discontinuous change of the behavior of the derivative of the stationary magnetization at the end points, with respect to the reaction rates. The dynamic phase transition is controlled by the spin flip rates of the particles at the end points, and is a discontinuous change of the relaxation time towards the stationary configuration. Other generalizations of the Glauber model consist of, for example, alternating-isotopic chains and alternating-bound chains (see , for example). People have also considered phase transitions induced by boundary conditions (see \[26–28\], for example).

The scheme of the paper is as follows. In section 2, the model is introduced, the rates are determined using the detailed balance, and the steady state configuration of the magnetization is obtained. In section 3, the dynamical phase transition of the system is investigated, and it is shown that it does show a dynamical phase transition, provided the temperature of the system is not zero.

## 2 Kinetic Ising model on a one-dimensional lattice with boundaries

The model being addressed, is the Glauber model on a one-dimensional lattice with boundaries. In the Glauber model, the interaction is between three neighboring sites. Spin flip brings the system into equilibrium with a heat bath at temperature $`T`$. A spin is flipped with a rate $`\mu =1\mathrm{tanh}(2\beta J)`$ whenever the spin of both of its neighboring sites are the same as it and is flipped with a rate $`\lambda =1+\mathrm{tanh}(2\beta J)`$ whenever the spin of both of its neighboring sites are in the opposite direction. (Here $`\beta :=1/(KT)`$.) At domain boundaries, spins are flipped with unit rate. So the interactions can be written as,

$`\text{ and }`$ with rate $`\mu `$
$`\text{ and }`$ with rate $`\lambda `$
$`\text{ and }`$ with rate $`1`$
Consider a lattice with $`L`$ sites and the Glauber dynamics as the interaction. The spin of the first site may flip with the following rates

$``$ with rate $`g_1`$
$``$ with rate $`g_2`$
$``$ with rate $`g_3`$
$``$ with rate $`g_4,`$
and the spin of the last site may flip with the following rates

$``$ with rate $`h_1`$
$``$ with rate $`h_2`$
$``$ with rate $`h_3`$
$``$ with rate $`h_4.`$
It is known that the time evolution equation for the one-point functions in the bulk are expressed in terms of only the one-point functions . To make this true for the boundaries as well, the following relations should hold.

$`g_1+g_4=`$ $`g_2+g_3`$
$`h_1+h_4=`$ $`h_2+h_3.`$ (1)
One may give a physical meaning to the parameters $`g_i`$ and $`h_i`$, by demanding the detailed balance to hold. Consider the energy $``$ of the system to be

$$=(B_1s_1+B_Ls_L+J\underset{i=1}{\overset{L1}{}}s_is_{i+1}),$$
(2)
then, the detailed balance demands

$`(s_1s_2s_1^{}s_2)\mathrm{exp}\{\beta (B_1s_1+Js_1s_2+\mathrm{})\}`$
$`=`$ $`(s_1^{}s_2s_1s_2)\mathrm{exp}\{\beta (B_1s_1^{}+Js_1^{}s_2+\mathrm{})\},`$ (3)
where $`(s_1s_2s_1^{}s_2)`$ is the rate of the spin flip of the first site from $`s_1`$ to $`s_1^{}`$. Equation (2) shows that

$$(s_1s_2s_1^{}s_2)=f(s_2)\mathrm{exp}\{\beta s_1(B_1+Js_2)\}.$$
(4)
The exponential term in the above equation is at most linear in terms of $`s_1`$. So,

$$(s_1s_2s_1^{}s_2)=\overline{f}(s_2)[1s_1\mathrm{tanh}\beta (B_1+Js_2)].$$
(5)
Then,

$`g_1`$ $`=\overline{f}(1)[1\mathrm{tanh}\beta (B_1J)]`$
$`g_2`$ $`=\overline{f}(1)[1\mathrm{tanh}\beta (B_1+J)]`$
$`g_3`$ $`=\overline{f}(1)[1+\mathrm{tanh}\beta (B_1+J)]`$
$`g_4`$ $`=\overline{f}(1)[1+\mathrm{tanh}\beta (B_1J)].`$ (6)
The condition of exact solvability (2) (the closure of time evolution equation of one-point functions) leads to

$$\overline{f}(1)=\overline{f}(1).$$
(7)
This means that the inertia of the first spin against spin-flip does not depend on the second spin. A similar expression can be written for the rate of the spin-flip of the last site.

For the infinite lattice, the Glauber model has a particle reaction-diffusion interpretation. If the spins of the neighboring sites are different (at a domain wall), one may consider the link between that sites as a particle. When the spins of the neighboring sites are the same (no domain wall), one may consider the link between the sites as a vacancy. Then the Glauber model turns into a reaction-diffusion model:

$``$ with rate $`1+\mathrm{tanh}(2\beta J)`$
$``$ with rate $`1\mathrm{tanh}(2\beta J)`$
$``$ with rate $`1,`$
where a particle (a vacancy) is denoted by $``$ ($``$). For the Glauber model with boundaries, to have a consistent particle model, one has to impose

$`g_1=g_3`$ $`g_2=g_4`$
$`h_1=h_3`$ $`h_2=h_4.`$ (8)
Then, the injection and extraction of particle at the first site are

$``$ with rate $`g_1=g_3`$
$``$ with rate $`g_2=g_4,`$ (9)
and the injection and extraction of particles at the last site are

$``$ with rate $`h_1=h_3`$
$``$ with rate $`h_2=h_4,`$ (10)
Now, consider the general case where only the conditions (2), which guarantee the closure of the time evolution, are satisfied. We have

$`\dot{s}_k=`$ $`2s_k+(s_{k+1}+s_{k1})\mathrm{tanh}(2\beta J)1<k<L`$
$`\dot{s}_1=`$ $`(g_2+g_3)s_1+(g_1g_2)s_2+(g_3g_1)`$
$`\dot{s}_L=`$ $`(h_2+h_3)s_L+(h_1h_2)s_{L1}+(h_3h_1).`$ (11)
The steady-state solution to (2) is

$$s_k=D_1z_1^k+D_2z_2^{kL1},$$
(12)
where

$$z_1=z_2^1=\mathrm{tanh}(\beta J).$$
(13)
It can be shown that in the thermodynamic limit ($`L\mathrm{}`$),

$`D_1={\displaystyle \frac{g_1g_3}{(g_1g_2)z_1^2(g_2+g_3)z_1}}`$
$`D_2={\displaystyle \frac{h_1h_3}{(h_1h_2)z_1^2(h_2+h_3)z_1}}.`$
$`D_1`$, and $`D_2`$ are continuous functions of the rates. So the behavior of $`s_k`$ near the ends of the lattice varies continuously with rates, and there is no phase transition.

## 3 The dynamical phase transition of the system

The average magnetization per site $`m(t)`$ is

$$m(t)=\frac{1}{L}\underset{k=1}{\overset{L}{}}s_k(t).$$
(14)
In the thermodynamic limit, the boundary terms are negligible, and

$$\frac{\mathrm{d}}{\mathrm{d}t}m(t)=2[\mathrm{tanh}(2\beta J)1]m(t).$$
(15)
Then, the same as Glauber model on an infinite lattice, the average magnetization do not show any phase transition. But, as it will be shown, the system does exhibit dynamical phase transition.

The homogeneous part of (2) can be written as

$$\dot{s}_k=h_k^ls_l.$$
(16)
The eigenvalues and eigenvectors of the operator $`h`$ satisfy

$`Ex_k`$ $`=2x_k+\mathrm{tanh}(2\beta J)(x_{k+1}+x_{k1}),k1,L`$
$`Ex_1`$ $`=(g_2+g_3)x_1+(g_1g_2)x_2,`$
$`Ex_L`$ $`=(h_2+h_3)x_L+(h_1h_2)x_{L1},`$ (17)
where the eigenvalue and the eigenvector have been denoted by $`E`$ and $`x`$, respectively. The solution to these is

$$x_k=az_1^k+bz_2^k,$$
(18)
where

$`(E+g_2+g_3)(az_1+bz_2)+(g_1g_2)(az_1^2+bz_2^2)`$ $`=0`$
$`(E+h_2+h_3)(az_1^L+bz_2^L)+(h_1h_2)(az_1^{L1}+bz_2^{L1})`$ $`=0,`$ (19)
and $`z_j`$’s satisfy

$$E=2+\mathrm{tanh}(2\beta J)(z+z^1).$$
(20)
So, $`z_1z_2=1`$. Two cases may occur,

a) $`z_1`$ and $`z_2`$ are phases.

b) $`z_1`$ and $`z_2`$ are both real, not equal to $`\pm 1`$. Then the modulus of one of them is less than one, that of the other is greater than one.

Using (20) and $`z_1z_2=1`$, one can eliminate $`E`$, and arrives at

$`z^{1L}[2g_2g_3+z(g_1g_2\mathrm{tanh}(2\beta J))z^1\mathrm{tanh}(2\beta J)]`$
$`\times [2h_2h_3+z(h_1h_2\mathrm{tanh}(2\beta J))z^1\mathrm{tanh}(2\beta J)]`$
$`z^{L1}[2g_2g_3+z^1(g_1g_2\mathrm{tanh}(2\beta J))z\mathrm{tanh}(2\beta J)]`$
$`\times [2h_2h_3+z^1(h_1h_2\mathrm{tanh}(2\beta J))z\mathrm{tanh}(2\beta J)]=0`$ . (21)
Obviously, $`z_j=\pm 1`$ satisfies (3). But these solutions lead to

$$x_k=z^k(a+bk).$$
(22)
And this form for $`x_k`$, generally does not satisfy the boundary conditions at $`k=1,L`$. Equation (3) can be written in the form

$$G(z):=F(z)F(z^1)=0,$$
(23)
where

$`F(z):=z^{1L}\{2g_2g_3+z[g_1g_2\mathrm{tanh}(2\beta J)]z^1\mathrm{tanh}(2\beta J)\}`$
$`\times \{2h_2h_3+z[h_1h_2\mathrm{tanh}(2\beta J)]z^1\mathrm{tanh}(2\beta J)\}.`$ (24)
For a phase solution to (3), $`z=e^{i\vartheta }`$, we have

$$E=2+2\mathrm{tanh}(2\beta J)\mathrm{cos}\vartheta .$$
(25)
In the thermodynamic limit $`(L\mathrm{})`$, in any neighborhood of $`z=1`$ there exist a phase solution to (3). The supermum of the eigenvalues determines the relaxation time toward the stationary average-density profile. So, if all of the solutions are phase,

$$\tau =[2+2\mathrm{tanh}(2\beta J)]^1$$
(26)
But, if there exist solutions which are not phases, they should be real. Consider $`z>1`$. Then for $`L\mathrm{}`$, (3) becomes

$`[2g_2g_3+`$ $`z^1(g_1g_2\mathrm{tanh}(2\beta J))z\mathrm{tanh}(2\beta J)]`$
$`\times [2h_2h_3+`$ $`z^1(h_1h_2\mathrm{tanh}(2\beta J))z\mathrm{tanh}(2\beta J)]=0.`$ (27)
Changing the rates, one may arrive at a situation where the above equation has a real solution greater than one. The transition occurs at the point that this equation has a solution equal to one. When the system has passed this point, the relaxation time becomes

$$\tau =[2+2(\mathrm{\Lambda }+\mathrm{\Lambda }^1)\mathrm{tanh}(2\beta J)]^1,$$
(28)
where $`\mathrm{\Lambda }`$ is that solution to (3), which is greater than one. (here we have assumed $`J>0`$, ferromagnetism. If $`J<0`$, anti-ferromagnetism, $`\mathrm{\Lambda }`$ is that solution to (3) which is less than $`1`$.) Putting $`z=1`$ in (3), at least one of the following equations should hold

$`2[1\mathrm{tanh}(2\beta J)]g_2g_4=0`$
$`2[1\mathrm{tanh}(2\beta J)]h_2h_4=0`$ (29)
If the temperature is zero, (3) for example gives $`g_2+g_4=0`$. Remembering that these parameters are rates, one arrives at $`g_2=g_4=0`$. So, at zero temperature, the solution cannot pass $`z=1`$. But at any other temperature, $`1\mathrm{tanh}(2\beta J)`$ is positive, and changing the parameters, $`g_2+g_4`$ can be made more than or less than $`1\mathrm{tanh}(2\beta J)`$.

If one uses the expressions (2) and (7) for $`g_i`$’s, then (3) becomes

$$2\mathrm{tanh}(2\beta J)+\overline{f}[\mathrm{tanh}\beta (JB_1)+\mathrm{tanh}\beta (J+B_1)]=0.$$
(30)
Putting $`\overline{f}=1`$, means that the inertia of the first spin against the spin flip is the same as those of the bulk spins. In this case, however, (30) has no solution. That is, there is no phase transition. In fact, (30) has no solution for $`\overline{f}1`$. For $`\overline{f}>1`$, however, it may have a solution.

It is seen that the parameters $`g_2`$ and $`g_4`$ (or $`h_2`$ and $`h_4`$) are control parameters of the dynamical phase transition. The parameters $`g_1`$ and $`g_3`$ (or $`h_1`$ and $`h_3`$) do not have any contribution in the dynamical phase transition. The rates $`g_1`$ and $`g_3`$ are the rates of the disappearance of the domain walls. But we note that the eigenvector corresponding to $`z=1`$ is a configuration where all the spins are the same ($`s_kz^k=1`$.) It is this configuration which corresponds to the largest value of $`E`$, which determines the relaxation time, and in this configuration, there is no domain wall. The disappearance rate of this configuration determines the relaxation time towards the steady state, and $`g_1`$ and $`g_3`$ (or $`h_1`$ and $`h_3`$) are irrelevant to this rate. In the particle-vacancy picture, this means that the rate of change of vacancy to particle is important, since the configuration corresponding to the maximum value of $`E`$ is the empty lattice.

This arguments are true for $`J>0`$, the ferromagnetic case. If $`J<0`$, then the relaxation time is determined by the value of $`E`$ at the smallest possible value of $`z`$ (which is less than $`1`$), and the transition occurs as $`z=1`$ becomes a solution to (3). It is not difficult to see that in this case $`g_1+g_3`$ (or $`h_1+h_3`$) determine the phase transition. The reasoning is the same as above, except that here the configuration determining the relaxation time is that corresponding to $`z=1`$, which means that the spins are alternating. So, in this configuration there are no $``$ or $``$ configurations and $`g_2`$ and $`g_4`$ (or $`h_2`$ and $`h_4`$) are irrelevant.

Acknowledgement The authors would like to thank Institute for Studies in Theoretical Physics and Mathematics for partial support.


format: text

text:
# An improved Parton Shower algorithm in QED

## 1 Introduction

The precise luminosity monitoring and data analysis at $`e^+e^{}`$ flavour factories (DA$`\mathrm{\Phi }`$NE, PEP-II, KEKB, CESR) require that cross section calculations for QED final state processes ($`e^+e^{}e^+e^{}`$, $`\mu ^+\mu ^{}`$, $`\gamma \gamma `$) have theoretical error below the 1% level. Also Monte Carlo simulations with the same accuracy are strongly demanded for data analysis. This implies to include the relevant radiative corrections, in particular the large effects due to multiple photon emission.

Aiming at providing a useful tool for data analysis, a Monte Carlo event generator (BABAYAGA) for QED processes at flavour factories has been developed. An exhaustive description of its main features and its basic theoretical ingredient, the QED Parton Shower (PS), can be found in . Here, we remind only of the basics of the theoretical approach which the generator is based on, in order to describe recent improvements introduced in the code to give a more precise simulation also for exclusive radiative events (i.e. with observed photons in the final state). They concern the inclusion of the effects due to interference between initial state and final state radiation, which are not present in a pure PS framework.

The outline of the paper is the following: in the next section, the main theoretical framework is summarized and the Parton Shower in QED is presented. In section 3, the generation of photons kinematics is revised and a new generation, inspired to the Yennie-Frautchi-Suura formula, is introduced. In section 4, the results obtained by the PS are compared to known perturbative results and the improvements due to the new generation are shown. The accuracy of the approach is also estimated. In the concluding section, the present study is summarized and the open issues and possible developments are discussed.

## 2 The Parton Shower algorithm in QED

A widely used theoretical tool, to take into account the bulk of the photonic radiative corrections, consists in the calculation of the corrected cross section following the master formula

$`\sigma (s)=`$ $`{\displaystyle }dx_{}dx_+dy_{}dy_+{\displaystyle }d\mathrm{\Omega }_{cm}D(x_{},Q^2)D(x_+,Q^2)\times `$ (1)
$`D(y_{},Q^2)D(y_+,Q^2){\displaystyle \frac{d\sigma _0}{d\mathrm{\Omega }_{cm}}}(x_{}x_+s,\vartheta _{cm})\mathrm{\Theta }(cuts),`$
which is based on the factorization theorems of universal infrared and collinear singularities. In the previous formula, $`d\sigma _0/d\mathrm{\Omega }`$ represents the Born-like differential cross section for the process under consideration and $`D(x,Q^2)`$ are the electron structure functions (SF) for initial and final state radiation. Equation (1) is also suited for the Monte Carlo generation of unweighted events.

The SF’s describe the effects of multiple emission of soft and hard photons in the collinear limit to all orders of perturbation theory. The SF is the solution of the QED Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation  in the non-singlet channel

$$Q^2\frac{}{Q^2}D(x,Q^2)=\frac{\alpha }{2\pi }_x^1\frac{dy}{y}P_+(y)D(\frac{x}{y},Q^2),$$
(2)
where $`P_+(x)`$ is the regularized $`ee+\gamma `$ splitting function. The electron SF has a clear and intuitive meaning: it represents the probability density of finding “inside” a parent electron an electron with momentum fraction $`x`$ and virtuality $`Q^2`$.

In the literature (see for example  and references therein) several solutions of DGLAP equation are available, but we focus our attention on its Monte Carlo solution, known as the Parton Shower algorithm . The advantages of the PS are mainly two: it is an exact numerical solution of the DGLAP equation and it allows an approximate generation of the photons momenta. The steps of the algorithm can be extracted introducing the Sudakov form factor  and, using it, getting an iterative solution of DGLAP equation. Namely, the solution is

$`D(x,Q^2)`$ $`=`$ $`\mathrm{\Pi }(Q^2,m^2)\delta (1x)+`$ (3)
$`+`$ $`{\displaystyle _{m^2}^{Q^2}}\mathrm{\Pi }(Q^2,s^{}){\displaystyle \frac{ds^{}}{s^{}}}\mathrm{\Pi }(s^{},m^2){\displaystyle \frac{\alpha }{2\pi }}{\displaystyle _0^{x_+}}𝑑yP(y)\delta (xy)+`$
$`+`$ $`{\displaystyle _{m^2}^{Q^2}}\mathrm{\Pi }(Q^2,s^{}){\displaystyle \frac{ds^{}}{s^{}}}{\displaystyle _{m^2}^s^{}}\mathrm{\Pi }(s^{},s^{\prime \prime }){\displaystyle \frac{ds^{\prime \prime }}{s^{\prime \prime }}}\mathrm{\Pi }(s^{\prime \prime },m^2)\times `$
$`(`$ $`{\displaystyle \frac{\alpha }{2\pi }})^2{\displaystyle }_0^{x_+}dx_1{\displaystyle }_0^{x_+}dx_2P(x_1)P(x_2)\delta (xx_1x_2)+({\displaystyle \frac{\alpha }{2\pi }})^3\mathrm{}`$
where $`\mathrm{\Pi }(s_1,s_2)=\mathrm{exp}[\alpha /2\pi _{s_2}^{s_1}𝑑s^{}/s^{}_0^{x_+}𝑑zP(z)]`$ is the Sudakov form factor, which represents the probability that the electron evolves from virtuality $`s_2`$ to virtuality $`s_1`$ with no emission of photons of energy fraction greater than (an infrared regulator) $`ϵ=1x_+`$. Note that the eq. (3) accounts for soft + virtual and real photons radiation up to all order of perturbation theory, in the leading logarithmic approximation. Concerning the soft + virtual cross section, it is worth noticing that, by setting the scale $`Q^2`$ to be equal to $`st/u`$, the Sudakov form factor exponentiates the leading logarithmic contribution of the $`𝒪(\alpha )`$ soft + virtual cross section as well as the dominant contribution coming from the infrared cancellation between the virtual box and the initial-final state interference of the bremsstrahlung diagrams. For a more detailed discussion on this topic, we refer to .

For the steps of our implementation of the PS algorithm, they can be found in . For the scope of this paper, it is sufficient to notice that the algorithm simulates a shower of photons emitted by the electron following the solution of eq. (3). When the algorithm stops, we are left with the energy fraction $`z_i`$ of each photon (distributed according to the Altarelli-Parisi splitting function), the virtualities of the electron at each branching and the remaining energy fraction $`x`$ of the electron after the showering. The $`x`$ variable is distributed according to $`D(x,Q^2)`$. By means of these quantities, an approximate branching kinematics can be obtained, as discussed below.

It is worth stressing that, when an event sample has been generated following equations (1) and (3), the sharing of the sample in elastic events (which correspond to soft + virtual cross section), one-photon events, two-photons events and so on, is automatic and built-in in the method itself. On the event sample, any experimental cut can be imposed, by means of the $`\mathrm{\Theta }`$ function which is present in eq. (1).

The reliability of the PS method to simulate and generate events can be checked by a systematic comparison of its results with known perturbative results, in particular with exact $`𝒪(\alpha )`$ matrix element (as discussed in section 4). The comparison allows to keep under control the approximations intrinsic into the PS approach and to estimate its theoretical accuracy. In order to perform such tests, an up to $`𝒪(\alpha )`$ PS algorithm has been developed as well. It allows to calculate the corrected cross section of eq. (1) up to $`𝒪(\alpha )`$. Such a calculation is strongly required for fully consistent comparisons between the PS predictions and the exact $`𝒪(\alpha )`$ ones. The steps required for the $`𝒪(\alpha )`$ PS, which are described in , can be obtained by using eq. (3) and expanding up to $`𝒪(\alpha )`$ the product $`D(x_{},Q^2)D(x_+,Q^2)D(y_+,Q^2)D(y_{},Q^2)`$ present in eq. (1).

## 3 Branching kinematics

The main advantage of the PS algorithm with respect to the collinear treatment of the electron evolution is the possibility of going beyond the strictly collinear approximation and generating transverse momentum $`p_{}`$ of electrons and photons at each branching. In fact, the kinematics of the branching process $`e(p)e^{}(p^{})+\gamma (q)`$ can be written as

$$p=(E,\stackrel{}{0},p_z),p^{}=(zE,\stackrel{}{p}_{},p_z^{}),q=((1z)E,\stackrel{}{p}_{},q_z).$$
(4)
Once the variables $`p^2`$, $`p_{}^{}{}_{}{}^{2}`$ and $`z`$ are generated by the PS algorithm, the on-shell condition $`q^2=0`$, together with the longitudinal momentum conservation, allows to obtain an expression for the $`p_{}`$ variable:

$$p_{}^2=(1z)(zp^2p^2)$$
(5)
at first order in $`p^2/E^21`$, $`p_{}^2/E^21`$. From now on, we will refer to eq. (5) as the “PS prescription”.

Some not correct behaviours of the exclusive photon kinematics reconstruction are connected with this PS picture, due to the approximations intrinsic in eq. (5). First of all, since within the PS algorithm the generation of $`p^2`$ and $`z`$ are independent, it can happen that in some branching the $`p_{}^2`$ as given by eq. (5) is negative. In order to avoid this problem, the introduction of any kinematical cut on $`p^2`$ or $`z`$ generation (or the regeneration of the whole event) would mean a not correct reconstruction of the SF $`x`$ distribution, which is important for a precise cross section calculation. Furthermore, in the PS scheme, each fermion produces its photons cascade independently from the other ones, missing the effects due to the interference of radiation coming from different charged particles. As far as inclusive cross sections (i.e., no cuts are imposed on the generated photons) are considered, these effects are integrated out and we do not expect them to be large, but when exclusive photon variables distributions are looked at, they may be important.

Concerning the first problem, it can be overcome choosing a generation of photons $`p_{}`$ different from eq. (5). For example, we can choose to extract the photon $`\mathrm{cos}\vartheta _\gamma `$ according to the universal leading poles $`1/pk`$ present in the matrix element for photon emission (in the following referred to as the “LL prescription”), as done for example in . Namely, we can generate $`\mathrm{cos}\vartheta _\gamma `$ as

$$\mathrm{cos}\vartheta _\gamma \frac{1}{1\beta \mathrm{cos}\vartheta _\gamma },$$
(6)
where $`\beta `$ is the speed of the emitting particle. In this way, photon energy and angle are generated independently, differently from eq. (5). The nice feature of this prescription is that $`p_{}^2=E_\gamma ^2\mathrm{sin}^2\vartheta _\gamma `$ is always well defined and the $`x`$ distribution reproduces exactly the SF, because we do not need to impose further kinematical cuts to avoid unphysical events. At this stage, the PS is used only to generate photons energies and photons multiplicity. The problem of including the radiation interference is still unsolved, because the variables of photons emitted by a fermion are still not correlated to the other charged particles. The phenomenological comparison between PS and LL prescriptions will be shown in the next section.

The issue of including photon interference can be successfully worked out looking at the Yennie-Frautschi-Suura (YFS) formula :

$$d\sigma _nd\sigma _0\frac{e^{2n}}{n!}\underset{l=1}{\overset{n}{}}\frac{d^\mathit{3}𝐤_l}{(2\pi )^32k_l^0}\underset{i,j=1}{\overset{N}{}}\eta _i\eta _j\frac{p_ip_j}{(p_ik_l)(p_jk_l)}.$$
(7)
It is a very general formula which gives the differential cross section $`d\sigma _n`$ for the emission of $`n`$ photons, whose momenta are $`k_1,\mathrm{},k_n`$, from a kernel process described by $`d\sigma _0`$ and involving $`N`$ fermions, whose momenta are $`p_1,\mathrm{},p_N`$. In eq. (7), $`\eta _i`$ is a charge factor, which is $`+1`$ for incoming $`e^{}`$ or outgoing $`e^+`$ and $`1`$ for incoming $`e^+`$ or outgoing $`e^{}`$. Note that eq. (7) is valid in the soft limit ($`k_i0`$). The important point is that it also accounts for coherence effects. From YFS formula, it is straightforward to read out the angular spectrum of the $`l^{th}`$ photon (“YFS prescription”):

$$\mathrm{cos}\vartheta _l\underset{i,j=1}{\overset{N}{}}\eta _i\eta _j\frac{1\beta _i\beta _j\mathrm{cos}\vartheta _{ij}}{(1\beta _i\mathrm{cos}\vartheta _{il})(1\beta _j\mathrm{cos}\vartheta _{jl})}.$$
(8)
It is worth noticing that in the LL prescription, the same quantity writes as

$$\mathrm{cos}\vartheta _l\underset{i=1}{\overset{N}{}}\frac{1}{1\beta _i\mathrm{cos}\vartheta _{il}}$$
(9)
whose terms are of course contained in eq. (8).

In order to consider also coherence effects in the angular distribution of the photons, we may generate $`\mathrm{cos}\vartheta _\gamma `$ according to eq. (8), rather than eq. (9). We expect that the improvements on photon angular distribution will be sizable, as discussed in the next section.

For the sake of clearness, it is worth stressing that we are not exploiting the so called “YFS exclusive exponentiation” procedure; we are simply considering the photons angular spectrum implied by equation (7) and we are generating the photons angular variables according to the distribution (8). Following such a procedure, we are embedding the more accurate angular spectrum (8) in the framework of the PS algorithm. If the reader is interested in the different approach of the YFS exclusive exponentiation, it is extensively and exhaustively treated in the papers of reference .

## 4 Comparisons and results

In this section, we deal with the application of the PS, mainly in its $`𝒪(\alpha )`$ realization, to the Bhabha and radiative Bhabha process at large angle, as implemented in BABAYAGA. To test the theoretical accuracy of the PS approach in total cross section calculation as well as in event generation, an exact $`𝒪(\alpha )`$ calculation has been carried out in the program LABSPV, used as a benchmark calculation. The code is described in  and it is based on the formulae given in .

As first comparison, we consider the results of PS and LL prescriptions on photon distributions. They are obtained by means of the PS algorithm truncated at $`𝒪(\alpha )`$. We generated a sample of radiative Bhabha events at a typical DA$`\mathrm{\Phi }`$NE energy ($`\sqrt{s}=M_\mathrm{\Phi }`$) and requiring the final state $`e^+`$ and $`e^{}`$ to lie within $`20^{}`$ and $`160^{}`$, with at least $`0.4`$ GeV of energy and a maximum acollinearity between them of $`10^{}`$. Furthermore, we require that the photon has at least $`10`$ MeV of energy. For those unphysical $`p_{}^2<0`$ events generated in the pure PS prescription, we arranged to set the photon collinear to the emitting particle ($`p_{}^2=0`$), without regenerating the event. In fig. 1, the energy and angular distributions of the photon are plotted. We note that, except some details, the two angular generation schemes give roughly the same answer: this is consistent with the PS picture, based on the leading logarithmic treatment of the branching kinematics.

Next, we consider how the distributions change if we weight photon angle by eq. (8) and we compare the results to the exact $`𝒪(\alpha )`$ distributions.

On the left, fig. 2 shows the angular distribution of the photon, imposing the same event selection criteria of the previous plots: the open circles represent the LL prescription, the stars represent the YFS prescription and the solid histogram is the exact $`𝒪(\alpha )`$ distribution. The improvement due to the inclusion of radiation coherence by means of YFS formula is evident: the YFS-PS spectrum fits perfectly to the exact $`𝒪(\alpha )`$ one. On the right, the photon energy distributions are compared: we note that the result of LL description is very good by itself (except close to the distribution end, where anyway the statistics is very poor), as a consequence of the good approximation of the Altarelli-Parisi splitting function to reproduce the correct energy spectrum. Nevertheless, the YFS-PS fits better to the exact distribution, even at its end. The same comments apply to fig. 3, where variables referring to final state fermions are plotted for the same sample of events: electron energy (on the left) and acollinearity (on the right) distributions are shown.

Up to now, we presented the results of the $`𝒪(\alpha )`$ PS, but the YFS prescription of eq. (8) has been applied to the all orders PS as well. As an example, we considered a sample of Bhabha events generated by means of the PS to all orders, imposing the same cuts as before and requiring that the most energetic photon of each event has at least 10 MeV of energy. The differences (due to coherence) on the angular distribution of the most energetic photon as given by the LL PS and the YFS PS are evident in fig. 4 (on the left). In the same figure, on the right, the effects of higher order photon emission are shown. The distribution of the energy lost because of radiation is plotted, comparing the $`𝒪(\alpha )`$ and the up to all order PS distribution. In this case, no cut is applied to the photon(s).

Finally, we give an estimate of the theoretical accuracy of the PS approach in inclusive cross section calculation. The bulk of the terms missed in the PS scheme in its all orders implementation can be estimated, up to very negligible uncertainty at next to leading $`𝒪(\alpha ^2)`$ level , by means of the relative difference between the exact $`𝒪(\alpha )`$ corrected cross section and the $`𝒪(\alpha )`$ PS one. In figure (5), such a difference is plotted, as a function of the acollinearity cut, for inclusive Bhabha cross section at $`\sqrt{s}=M_\mathrm{\Phi }`$, requiring the final state $`e^+`$ and $`e^{}`$ to have at least $`0.4`$ GeV of energy and considering $`20^{}<\vartheta _\pm <160^{}`$ and $`50^{}<\vartheta _\pm <130^{}`$ acceptance region. For the “best” prescription, the YFS one, the amount of $`𝒪(\alpha )`$ missing contributions is contained within the 0.4% and they could be further reduced, by an ad hoc and tuned choice of the $`Q^2`$ scale present in the electron SF. This topic is widely discussed in . It is worth noticing that the exact versus PS difference is less dependent on the acollinearity cut when the YFS-improved PS is considered. The next step toward the reduction of the theoretical error (and the complete independence from ad hoc recipes) would be the real merging of the exact $`𝒪(\alpha )`$ matrix element into the PS. Work is in progress in this direction.

## 5 Conclusions

The QED Parton Shower algorithm has been applied to QED processes event generation at $`e^+e^{}`$ flavour factories, in order to simulate the phenomenologically relevant radiative corrections to all orders of perturbation theory and to provide a useful Monte Carlo event generator (BABAYAGA) for data analysis. The theoretical accuracy of the event generator is estimated to be within the 0.5% level for Bhabha scattering, in typical event selection criteria for data analysis at flavour factories. After a brief sketch of the theoretical framework which the PS algorithm is based on, in this paper the photons kinematics generation is critically revised. We introduced three different schemes for generating photons angles, named the pure PS, the LL and the YFS prescription. In particular, the last one is inspired to the YFS formula  of eq. (7). By means of an $`𝒪(\alpha )`$ PS, the results from the three prescriptions have been systematically compared to the exact $`𝒪(\alpha )`$ ones. We pointed out and showed that the YFS-PS has the best behaviour to simulate and generate also radiative events, fitting very well to all the exact $`𝒪(\alpha )`$ distributions. The new release of the BABAYAGA event generator is built on this theoretical background. The code is available and can be downloaded from the web site http://www.pv.infn.it/~nicrosi/programs.html.

The open issue for the next future is the merging of the exact $`𝒪(\alpha )`$ matrix element with the PS, which is an interesting task also for QCD simulations at high energies, as documented in . The merging would make the PS independent from ad hoc recipes and would actually shift the theoretical uncertainty to the next to leading order $`\alpha ^2`$ corrections.

## 6 Acknowledgments

My work is done in collaboration with Guido Montagna, Oreste Nicrosini and Fulvio Piccinini, whom I wish to thank for their constant support to my work. I’m grateful also to Achim Denig, of the KLOE Collaboration, for his interest in our work and very useful discussions. Finally, I thank the INFN, Sezione di Pavia, for the use of the computer facilities. The author acknowledges partial support from the EEC-TMR Program, Contract N. CT98-0169.


text:
# Theory of resonant Raman scattering of tetrahedral amorphous carbon

## I Introduction

Tetrahedral amorphous carbon is a very studied material because of its physical properties, which are close to those of diamond. It is composed of carbon atoms with $`sp^3`$ and $`sp^2`$ hybridization, with a fraction of $`sp^3`$ sites larger than 40 %. The properties of tetrahedral amorphous carbon are very dependent on its microscopic structure. Thus, the knowledge of the microscopic structure can be used to improve the synthesis of the material. Raman spectroscopy have often been applied to the study of tetrahedral amorphous carbon. However visible Raman spectroscopy is much more sensitive to $`sp^2`$ than to $`sp^3`$ carbons. Recently, it has been shown that it is possible to probe the vibrations of $`sp^3`$ carbons with UV Raman, but few UV spectra have been presented in literature up to now. In absence of further experimental data, a theoretical study can elucidate the detailed dependence of the Raman spectra on the laser frequency.

A theory to compute the Raman spectra under non-resonant conditions is well established. Indeed, the Placzek approximation links the scattered Raman intensity with the electronic susceptibility of the sample. However, for resonant Raman the Placzek approximation is not justified. In this work, we develop a practical method for the calculation of resonant Raman spectra in solids with delocalized excitations. We apply this theory to amorphous carbon, using density functional theory in the local density approximation to compute the vibrational modes and frequencies and a tight binding approximation to obtain the resonant Raman intensities. We compare our results to the experimental visible and UV Raman spectra and we analyze our theoretical spectra in terms of the motion of microscopic units.

## II Theoretical results

In Raman spectroscopy a sample of matter is irradiated by a monochromatic laser beam of pulsation $`\omega _\mathrm{L}`$, and the intensity of the scattered light $`I(\omega )`$ is measured, as a function of the pulsation $`\omega `$. The intensity $`I(\omega )`$ has a main elastic peak at the incident light pulsation, $`\omega =\omega _\mathrm{L}`$, and a smaller inelastic contribution at $`\omega \omega _\mathrm{L}`$ associated to energy transfers between the light and the sample. In this paper, we will focus on first-order vibrational Stokes spectroscopy, in which the light excites a single vibrational mode.

In order to describe the vibrational Raman scattering and to fix the notation, we need to introduce the Born-Oppenheimer (BO) approximation. Within the BO approximation the eigenstates of the Hamiltonian with eigenvalues $`E_k^e`$ can be written as $`|\psi _{e,𝐑}|\nu _{e,k}`$, where the kets $`|\psi _{e,𝐑}`$ and $`|\nu _{e,k}`$ are defined on the Hilbert space of the electronic and nuclear coordinates, respectively. The kets $`|\psi _{e,𝐑}`$ are the eigenvectors with eigenvalues $`E_e(𝐑)`$ of the electronic Hamiltonian $`H^{el}(𝐑)`$, which depends parametrically on the positions $`𝐑`$, of the N nuclei ($`𝐑`$ represents a vector of 3N coordinates). For every $`e`$, each $`E_e(𝐑)`$ defines a BO potential energy surface. Finally the wave-functions $`|\nu _{e,k}`$ satisfy the eigenvalue equation:

$$[T+E_e(𝐑)]|\nu _{e,k}=E_k^e|\nu _{e,k},$$
(1)
where $`T`$ is the kinetic energy of the nuclei. It is possible to define a vibrational energy as:

$$ϵ_k^e=E_k^eE_e(𝐑_e^{\mathrm{eq}}),$$
(2)
where $`𝐑_e^{\mathrm{eq}}`$ is the equilibrium position of the nuclei on the BO surface $`E_e(𝐑)`$.

In vibrational Raman spectroscopy, before the scattering, the system is on the ground state BO surface, $`E_0(𝐑)`$, and in the vibrational state $`ϵ_i^0`$, after the scattering the system is still on the ground state BO surface but in a different vibrational state $`ϵ_f^0`$. In first-order Stokes Raman the harmonic approximation is assumed and the final state differs from the initial one by the creation of a phonon with energy $`\omega _s`$.

In non-resonant Raman scattering the energy of incident light is very far from any electronic excitation. Under this condition, a simple and well established theory to compute the scattering intensity $`I(\omega )`$, have been developed by Placzek. Within the Placzek’s approximation the intensity of scattered light is expressed in terms of the electronic susceptibility at the laser frequency $`\chi _{𝐮_\mathrm{L}𝐮}(\omega _\mathrm{L})`$:

$$I(𝐮_\mathrm{L},𝐮,\omega )\underset{s}{}\frac{\omega ^4}{\omega _s}\left|\frac{\chi _{𝐮_\mathrm{L}𝐮}(\omega _\mathrm{L})}{Q_s}\right|^2(n_s+1)\delta (\omega \omega _\mathrm{L}+\omega _s),$$
(3)
where $`𝐮_\mathrm{L}`$ and $`𝐮`$ are the polarization of the incident and scattered light, $`n_s=[\mathrm{exp}(\omega _s/(k_BT))1]^1`$ is the average thermal occupation of the vibrational state $`s`$, and the derivative of the susceptibility with respect to the vibrational normal mode $`Q_s`$ is:

$$\frac{\chi _{𝐮_\mathrm{L}𝐮}(\omega _\mathrm{L})}{Q_s}=\underset{J}{}\underset{\alpha =x,y,z}{}\frac{\chi _{𝐮_\mathrm{L}𝐮}(\omega _\mathrm{L})}{R_J^\alpha }q_J^{\alpha ,s}=\frac{\chi _{𝐮_\mathrm{L}𝐮}(\omega _\mathrm{L})}{𝐑}𝐪^s.$$
(4)
Here $`q_J^{\alpha ,s}`$ is a component of the unit vector of normal mode $`Q_s`$ for vibration $`s`$ on the atom J with normalization $`𝐪^s𝐪^s^{}=\delta _{s,s^{}}`$. The Placzek’s approximation have been used successfully for the first principles computation of the non-resonant vibrational Raman intensity, (see e.g. Refs. ).

In resonant Raman spectroscopy the energy of incident light is in resonance with an electronic transition. Under this condition the Placzek’s approximation is, in principle, not justified. However Eq.(3) has also been applied to resonant Raman scattering in solids in order to interpret experimental measurements on semiconductors.

In this paper, we will derive an expression for the scattering intensity in solids which is valid both in the non-resonant and resonant case. This derivation will be used in the following section to analyze the Raman spectra of amorphous carbon.

We start our derivation from the general expression for the vibrational Raman intensity, which can be obtained from theory of second quantization of light:

$`I(𝐮_\mathrm{L},𝐮,\omega )`$ $``$ $`{\displaystyle \underset{i,f}{}}\rho (i)\omega ^4|𝒜_{i,f}|^2\delta (\omega \omega _\mathrm{L}ϵ_i^0+ϵ_f^0),`$ (5)
with

$`𝒜_{i,f}={\displaystyle \underset{e,k}{}}({\displaystyle \frac{\nu _{0,i}|f_{𝐮_\mathrm{L}}^e(𝐑)|\nu _{e,k}\nu _{e,k}|f_𝐮^e(𝐑)|\nu _{0,f}}{E_e(𝐑_e^{\mathrm{eq}})+ϵ_k^eE_0(𝐑_0^{\mathrm{eq}})+ϵ_i^0\omega _\mathrm{L}i\gamma }}`$ (7)
$`+{\displaystyle \frac{\nu _{0,i}|f_𝐮^e(𝐑)|\nu _{e,k}\nu _{e,k}|f_{𝐮_\mathrm{L}}^e(𝐑)|\nu _{0,f}}{E_e(𝐑_e^{\mathrm{eq}})+ϵ_k^eE_0(𝐑_0^{\mathrm{eq}})+ϵ_i^0+\omega +i\gamma }}).`$
Here, $`\rho (i)`$ is the probability to find the system in the initial state $`i`$, $`\gamma `$ is a small real number added for calculation in order to treat correctly the poles of the expression and

$$f_𝐮^e(𝐑)=\frac{1}{\sqrt{V}}\psi _{0,𝐑}|𝐮𝐃(𝐫,𝐑)|\psi _{e,𝐑},$$
(8)
where $`𝐃`$ is the ionic and electronic dipole moment, $`V`$ is the sample volume, and $`𝐫`$ is the position vector for the electrons. Placzek’s expression, Eq. (3), has been derived in the non resonant case starting from Eq. (7). Notice that Eq. (7) is much more complicated to be evaluated than the Placzek’s expression, Eq. (3), since it requires the knowledge of the vibrational eigenfunctions, $`|\nu _{e,k}`$, on the excited BO surfaces.

Before presenting our derivation we recall, for comparison, that of the Placzek’s approximation. In the Placzek derivation Eq. (7) is approximated by:

$`𝒜_{i,f}={\displaystyle \underset{e}{}}({\displaystyle \frac{_k\nu _{0,i}|f_{𝐮_\mathrm{L}}^e(𝐑)|\nu _{e,k}\nu _{e,k}|f_𝐮^e(𝐑)|\nu _{0,f}}{E_e(𝐑)E_0(𝐑)\omega _\mathrm{L}i\gamma }}`$ (10)
$`+{\displaystyle \frac{_k\nu _{0,i}|f_𝐮^e(𝐑)|\nu _{e,k}\nu _{e,k}|f_{𝐮_\mathrm{L}}^e(𝐑)|\nu _{0,f}}{E_e(𝐑)E_0(𝐑)+\omega +i\gamma }}).`$
Using the completeness relation $`_k|\nu _{e,k}\nu _{e,k}|=1`$, Eq. (10) becomes:

$$𝒜_{i,f}=\underset{e}{}\left(\frac{\nu _{0,i}|f_{𝐮_\mathrm{L}}^e(𝐑)f_𝐮^e(𝐑)|\nu _{0,f}}{E_e(𝐑)E_0(𝐑)\omega _\mathrm{L}i\gamma }+\frac{\nu _{0,i}|f_𝐮^e(𝐑)f_{𝐮_\mathrm{L}}^e(𝐑)|\nu _{0,f}}{E_e(𝐑)E_0(𝐑)+\omega +i\gamma }\right),$$
(11)
which finally leads to Eq. (3). In order to obtain Eq. (10), we drop the dependence of the denominators of Eq. (7) on the excited vibrational states $`|\nu _{e,k}`$, by replacing $`[E_e(𝐑_e^{\mathrm{eq}})+ϵ_k^eE_0(𝐑_0^{\mathrm{eq}})+ϵ_i^0]`$ with $`[E_e(𝐑)E_0(𝐑)]`$. In this way, we neglect in the denominators an energy term

$$\mathrm{\Delta }_e^k=E_e(𝐑_e^{\mathrm{eq}})E_e(𝐑)+ϵ_k^e[E_0(𝐑_0^{\mathrm{eq}})E_0(𝐑)+ϵ_i^0].$$
(12)
This approximation is justified if $`E_e(𝐑)E_0(𝐑)\omega _\mathrm{L}\mathrm{\Delta }_e^k`$. Since we are using the completeness relation, the above inequality should, in principle, held for every vibrational states $`k`$, and not just for the one phonon excitations. In particular, the relation must be verified for excitations, and energy differences $`\mathrm{\Delta }_e^k`$, involving an arbitrary number of phonons and even for unbounded states in the energy continuum. This condition can obviously never be satisfied. The use of Placzek approximation relies on the hope that the inclusion of excitations involving a limited amount of phonons are sufficient to fulfill the completeness relation. When the laser frequency is close to an electronic transition, the Placzek’s approximation could not be applied in general since, in this case, the dependence on $`\mathrm{\Delta }_e^k`$ could not be ignored and the completeness relation could not be used.

Here, to compute the Raman spectra in a periodic solid, we simplify Eq. (7) using an approximation different from Placzek’s. In a periodic solid with dispersive bands, the electronic excitations are generally delocalized. In this case, the excited BO surfaces are obtained from the ground-state one by a constant vertical displacement, and we can assume that $`|\nu _{e,k}|\nu _{0,k}`$, $`𝐑_e^{\mathrm{eq}}𝐑_0^{\mathrm{eq}}`$, and $`ϵ_k^eϵ_k^0`$. This leads to:

$`𝒜_{i,f}(𝐮_\mathrm{L},𝐮)={\displaystyle \underset{e,k}{}}({\displaystyle \frac{\nu _{0,i}|f_{𝐮_\mathrm{L}}^e(𝐑)|\nu _{0,k}\nu _{0,k}|f_𝐮^e(𝐑)|\nu _{0,f}}{E_e(𝐑_0^{\mathrm{eq}})E_0(𝐑_0^{\mathrm{eq}})+ϵ_k^0ϵ_i^0\omega _\mathrm{L}i\gamma }}`$ (14)
$`+{\displaystyle \frac{\nu _{0,i}|f_𝐮^e(𝐑)|\nu _{0,k}\nu _{0,k}|f_{𝐮_\mathrm{L}}^e(𝐑)|\nu _{0,f}}{E_e(𝐑_0^{\mathrm{eq}})E_0(𝐑_0^{\mathrm{eq}})+ϵ_k^0ϵ_i^0+\omega +i\gamma }}).`$
We expand the expression of $`f_𝐮^e(𝐑)`$ at first order around $`𝐑_0^{\mathrm{eq}}`$, $`f_𝐮^e(𝐑)f_𝐮^e(𝐑_0^{\mathrm{eq}})+(𝐑𝐑_0^{\mathrm{eq}})f_𝐮^e(𝐑_0^{\mathrm{eq}})/𝐑`$ and we extract all the terms of degree one in $`(𝐑𝐑_0^{\mathrm{eq}})`$ in Eq. (14). The term $`𝐑𝐑_0^{\mathrm{eq}}`$ can be written using the operators of creation and annihilation of one phonon of energy $`\omega _s`$. If we extract Stokes terms, we can write the amplitude as a function of a vibration $`s`$,

$`I(𝐮_\mathrm{L},𝐮,\omega )`$ $``$ $`{\displaystyle \underset{s,n_s}{}}\rho (n_s)\omega ^4|𝒜_s|^2\delta (\omega \omega _\mathrm{L}+\omega _s),`$ (15)
where $`n_s`$ is the occupation number of the mode $`s`$, $`\omega _s=ϵ_f^0ϵ_i^0`$, and

$`𝒜_s`$ $`=`$ $`{\displaystyle \underset{e}{}}\sqrt{{\displaystyle \frac{n_s+1}{2\omega _s}}}({\displaystyle \frac{𝐪_s(\frac{f_{𝐮_\mathrm{L}}^e}{𝐑})_{𝐑_0^{\mathrm{eq}}}f_𝐮^e(𝐑_0^{\mathrm{eq}})}{E_e(𝐑_0^{\mathrm{eq}})E_0(𝐑_0^{\mathrm{eq}})+\omega _s\omega _\mathrm{L}i\gamma }}+{\displaystyle \frac{f_{𝐮_\mathrm{L}}^e(𝐑_0^{\mathrm{eq}})𝐪_s(\frac{f_𝐮^e}{𝐑})_{𝐑_0^{\mathrm{eq}}}}{E_e(𝐑_0^{\mathrm{eq}})E_0(𝐑_0^{\mathrm{eq}})\omega _\mathrm{L}i\gamma }}`$ (16)
$`+`$ $`{\displaystyle \frac{𝐪_s(\frac{f_𝐮^e}{𝐑})_{𝐑_0^{\mathrm{eq}}}f_{𝐮_\mathrm{L}}^e(𝐑_0^{\mathrm{eq}})}{E_e(𝐑_0^{\mathrm{eq}})E_0(𝐑_0^{\mathrm{eq}})+\omega _\mathrm{L}+i\gamma }}+{\displaystyle \frac{f_𝐮^e(𝐑_0^{\mathrm{eq}})𝐪_s(\frac{f_{𝐮_\mathrm{L}}^e}{𝐑})_{𝐑_0^{\mathrm{eq}}}}{E_e(𝐑_0^{\mathrm{eq}})E_0(𝐑_0^{\mathrm{eq}})+\omega _\mathrm{L}+i\gamma \omega _s}}).`$ (17)
Finally, we neglect the dependence of the denominator over $`\omega _s`$. Indeed, we expect that $`𝒜_s`$ has the same behavior with respect to $`\omega _L`$ as the electronic susceptibility, $`\chi (\omega _L)`$. In a solid with delocalized excitations, away from the van Hove singularities, $`\chi (\omega _L)`$ is a smooth function such that $`\chi (\omega _L\pm \omega _s)\chi (\omega _L)`$. Therefore, we expect that in a solid with delocalized excitations we can drop the term $`\pm \omega _s`$ in the denominator of Eq. (16). Notice that the validity requirements of our approximation are less stringent than those of the Placzek’s approximation. Our approximation neglects the dependence of the denominators on the one phonon excitation energies $`\omega _s`$, whereas the Placzek’s approximation neglects the multi-phonons excitation energies $`\mathrm{\Delta }_e^k`$. Therefore our approximation requires that $`\chi (\omega )`$ is flat on the energy scale of one phonon excitations whereas the Placzek’s approximation requires that $`\chi (\omega )`$ is flat on the energy scale of multiple phonon excitations.

$$𝒜_s=\sqrt{\frac{n_s+1}{2\omega _s}}𝐪_s\frac{\kappa (𝐑,𝐑^{},\omega _\mathrm{L})}{𝐑}|_{𝐑=𝐑^{}=𝐑_0^{\mathrm{eq}}},$$
(18)
where

$$\kappa _{𝐮_\mathrm{L},𝐮}(𝐑,𝐑^{},\omega _\mathrm{L})=\underset{e}{}\frac{f_{𝐮_\mathrm{L}}^e(𝐑)f_𝐮^e(𝐑)}{E_e(𝐑^{})E_0(𝐑^{})\omega _\mathrm{L}i\gamma }+\frac{f_𝐮^e(𝐑)f_{𝐮_\mathrm{L}}^e(𝐑)}{E_e(𝐑^{})E_0(𝐑^{})+\omega _\mathrm{L}+i\gamma }.$$
(19)
Eqs. (15), (18), and (19) are the working expressions that we will use in the next section to compute the Raman intensity.

It is interesting to compare our final expression with that of the Placzek’s approximation. To this purpose we notice that the electronic susceptibility $`\chi `$ can be obtained as a special case of the function $`\kappa `$, indeed $`\chi _{𝐮_\mathrm{L},𝐮}(𝐑,\omega _\mathrm{L})=\kappa _{𝐮_\mathrm{L},𝐮}(𝐑,𝐑,\omega _\mathrm{L})`$. Thus, the Placzek’s expression differs from the expression we derived just by an additional partial derivative. Indeed, if we assume that

$$𝒜_s^{\mathrm{Placzek}}=\sqrt{\frac{n_s+1}{2\omega _s}}𝐪_s\left[\frac{\kappa (𝐑,𝐑^{},\omega _\mathrm{L})}{𝐑}+\frac{\kappa (𝐑,𝐑^{},\omega _\mathrm{L})}{𝐑^{}}\right]_{𝐑=𝐑^{}=𝐑_0^{\mathrm{eq}}},$$
(20)
we obtain Eq. (3) from Eq. (15).

## III Calculation on amorphous carbon

We compute the vibrational spectra and the Raman intensities of tetrahedral amorphous carbon using a model generated in Ref.  by tight binding (TB) molecular dynamics. The model contains 64 atoms in a periodic cubic cell of 7.57 Å. Among the 64 atoms, 18 atoms are three-fold coordinated ($`sp^2`$ hybridized) and the others are four-fold coordinated ($`sp^3`$ hybridized). In this model all $`sp^2`$ atoms are inserted into a double bonded link, i.e. there are no electronic defects associated to dangling $`p`$ orbitals. This kind of defects has been identified in other computed generated models of amorphous carbon. Fig. 1 shows a ball and stick model of the $`sp^2`$ carbon atoms in the cell. The $`sp^2`$ atoms are arranged in four pairs and two chains of four and six atoms respectively. No aromatic ring are present in the sample. We compute the dynamical matrix and the vibrational eigenvalues and eigenvectors using density functional theory in the local density approximation. The C atoms are described by a Troullier Martins pseudo-potential with $`s`$ non-locality. The wave-functions are expanded in a plane wave basis set with a kinetic energy cut-off of 40 Ry. The dynamical matrix is obtained by displacing each atom and computing the resulting forces. The Brillouin zone is sampled with the $`\mathrm{\Gamma }`$ $`𝐤`$-point only. We compute the intensities of the resonant Raman spectra using Eqs. (15), (18), and (19). The evaluation of the intensity at a resonant frequency requires a very fine Brillouin zone k-point mesh, for which a fully ab-initio calculation is not affordable. Therefore we evaluate the function $`\kappa _{𝐮_\mathrm{L},𝐮}(𝐑,𝐑^{},\omega _\mathrm{L})`$ and the electronic susceptibility $`\chi _{𝐮_\mathrm{L},𝐮}(𝐑,\omega _\mathrm{L})`$ within a TB approximation. We use a TB model with s and p orbitals. We verified that this TB Hamiltonian well reproduces the ab-initio random-phase-approximation (RPA) dielectric constant of diamond for laser frequencies smaller than 10 eV, which is the range of energies used in this work for the calculation of the Raman spectra. The correct description of the dielectric constant for frequencies larger than 10 eV would require the inclusion in the TB model of d or s orbitals. The Brillouin zone of the 64 atoms supercell is sampled with 108 special k-points. We use a value of 0.05 eV for the constant $`\gamma `$ in Eq. (19). The derivative in Eq. (18) is computed by numerical differentiation.

In Fig. 2 we present our theoretical vibrational density of state (DOS). The total DOS is decomposed in its partial contributions coming from the $`sp^2`$ and $`sp^3`$ hybridized carbons and from the stretching and bending modes. The $`sp^2`$ vibrations cover the entire range of frequencies of the total DOS, whereas the $`sp^3`$ vibrations form a large peak centered at 1000 $`\text{cm}^1`$, which dies over 1400 $`\text{cm}^1`$. Thus, the modes above 1400 $`\text{cm}^1`$ involve only $`sp^2`$ carbons. The decomposition between bending and stretching contributions shows two separated peaks, centered at 600 $`\text{cm}^1`$ and 1000 $`\text{cm}^1`$, respectively. Above 1000 $`\text{cm}^1`$ the stretching modes contribute alone to the DOS. Similar theoretical results have been presented in Ref.  for a model of hydrogenated amorphous carbon.

Before discussing the results for the Raman spectra, we present in Fig. 3 the real and imaginary part of the electronic dielectric constant, $`ϵ(\omega _L)=1+4\pi \chi (\omega _L)`$. Our theoretical value of 5.8 for the static dielectric constant is in good agreement with the value of 6.2, measured experimentally for a tetrahedral amorphous carbon sample. Notice that, for $`\omega >5`$ eV, the function $`ϵ(\omega _L)`$ is quite smooth, whereas, under 5 eV, three resonant peaks at 1.3, 3.0 and 3.9 eV show up. The presence of such peaks, which involve transitions between localized $`\pi `$ states, is an artifact due to the limited size of our model which contains only 9 $`\pi `$-bonds. This kind of features should disappear in larger models. The electronic structure has a gap of about 0.9 eV, above this value the incident light can be absorbed by the sample, as shown by the non zero value of the imaginary part of the dielectric constant. Thus, Raman spectra both in the visible and ultra-violet (UV) regions are collected under resonant conditions.

Fig. 4 and 5 present the calculated Raman spectra of amorphous carbon for an incident laser beam of 1.8 eV and 4.3 eV, respectively. The laser beam energies are chosen to be away from the three resonant peaks of $`ϵ(\omega _L)`$. For each spectrum, we show a decomposition between $`sp^2`$, $`sp^3`$, bending, and stretching vibrations. We also plot the overlap between the different contributions. The decompositions are meaningful only if the overlap is small. Again, the details of the shape of the Raman spectra coming from $`sp^2`$ atoms are affected by the limited number of $`sp^2`$-$`sp^2`$ bonds in our model. For example the division of the contribution above 1300 $`\text{cm}^1`$ in two peaks is an artifact due to the lack of statistic. In good agreement with the experimental observation, the visible Raman spectrum is essentially due to $`sp^2`$ carbons which give a large contribution above 1200 $`\text{cm}^1`$, centered at 1500 $`\text{cm}^1`$. Only a little $`sp^3`$ contribution, centered at 1000 $`\text{cm}^1`$, can be noticed. In the UV Raman spectrum, on the contrary, a contribution of $`sp^3`$ hybridized atoms can clearly be noticed through the large peak centered at 1000 $`\text{cm}^1`$. Due to a larger number of $`sp^3`$ hybridized atoms in our model, the shape of this peak is more reliable, and its position and shape closely match the experimental results.  A peak centered at 1600 $`\text{cm}^1`$, due to $`sp^2`$ sites, can still be observed, as in experiments. Both visible and UV Raman spectra are in a large majority due to the stretching vibrations, only a little contribution of bending vibrations can be noticed at about 600 $`\text{cm}^1`$, which has sometimes been observed in the experimental spectra (see e.g. Fig. 1 of Ref. ).

For a given laser pulsation, $`\omega _\mathrm{L}`$, we compute the integrated Raman Stokes intensities as the total area of the Stokes spectra. The total integrated Raman intensity and its decomposition in $`sp^2`$, $`sp^3`$, bending, and stretching contributions are presented in Figs. 6 and 7 as a function of $`\omega _\mathrm{L}`$. These data show that whatever the energy of the incident light, the intensity of the Raman spectrum is always dominated by stretching modes, whose contribution is always at least an order of magnitude larger than that of bending modes. Regarding the relative decomposition in the $`sp^2`$ and $`sp^3`$ contributions, Fig. 6 and 7 clearly show an inversion at about 5 eV. Under 4 eV almost 90% of the intensity is due to $`sp^2`$ carbons, whereas over 6 eV almost 90% of the intensity is due to $`sp^3`$ carbons. Notice that this inversion happens very quickly, thus, in an experimental measurement a little increase in the energy of the incident light in the UV region can lead to a major modification of the Raman spectra.

Finally, in Fig. 8, we compare the total Raman intensities computed using our method and the Placzek’s approximation. The Raman spectra computed with the Placzek’s approximation show stronger singularities in correspondence to the $`sp^2sp^2`$ electronic transitions. Instead, the intensity of the spectra computed with our approach is similar for laser frequency in resonance with $`sp^2sp^2`$ transitions and with $`sp^3sp^3`$ transitions. We did not find any experimental data for the total intensity as a function of the laser frequency, but according to the discussion of Section II the use of our approximation is more justified in the resonant case.

## IV conclusion

In this paper, we have presented a practical method to compute the vibrational resonant Raman spectra in solids. Our approach is justified if the electronic excitations are delocalized. We use our method to compute the vibrational resonant Raman spectra of tetrahedral amorphous carbon. The computed spectra are in good agreement with the experimental ones measured with visible and UV lasers. In particular, the theoretical visible Raman spectra present a broad feature between 1200 and 1700 cm<sup>-1</sup> associated with $`sp^2`$ stretching as found in experiments. In the theoretical UV Raman spectra a second broad peak around 1000 cm<sup>-1</sup> associated to $`sp^3`$ stretching shows up. The location and the shape of this peak agree very well with the experimental UV spectra. We have analyzed in detail the evolution of the Raman spectra as a function of the laser frequency. We have shown that, at any frequency, the spectra are dominated by the stretching vibrations. We have identified a very rapid inversion in the relative Raman intensities of the $`sp^2`$ and $`sp^3`$ sites with the frequency of the incident laser beam. In particular, the spectra are dominated by $`sp^2`$ atoms below 4 eV and by $`sp^3`$ atoms above 6 eV. According to our results, it would be interesting to collect a UV spectrum of tetrahedral amorphous carbon with a laser frequency larger than those used in the actual experiments, to further emphasize the contribution of the tetrahedral $`sp^3`$ carbon sites.

We acknowledge D. A. Drabold for providing us with the model of tetrahedral amorphous carbon. We thank Dr. M. Marangolo for a critical reading of the manuscript. The calculations have been performed at the Idris computer center of the CNRS.