text:
# On the construction of 𝒩=4 SYM effective action beyond leading low-energy approximation 11footnote 1Talk given at International Workshop ”Supersymmetries and Quantum Symmetries” - SQS’03 (Dedicated to the 75th anniversary of the birth of V.I.Ogievetsky); JINR, Dubna; July 24-29, 2003

## 1 Introduction

Different aspects of low-energy string dynamics and of the AdS/CFT correspondence  can be studied in terms of the quantum field theory effective action. Due to the proposition put forth in (see discussed in greater detail ), the superconformal version of the Born-Infeld (BI) action is expected to be considered as the result of summing up leading and subleading terms in the quantum effective action of the $`𝒩=4`$ SYM theory in the Coulomb branch. Some terms of the effective action expansion should be constrained by new non-renormalization theorems.

In the off-shell superfield formulation of the SYM and supergravity theories, the supersymmetry must be realized linearly on physical fields and on an infinite set of auxiliary fields so that the supersymmetry transformations are independent of the form of the action. But in the on-shell formalism, the supersymmetry transformation is realized non-linearly. When obtaining higher derivative contributions to the effective action preserving extended supersymmetry, we must obtain self-consistently the deformations of the classical supersymmetry transformation rules order by order and, simultaneously, construct the supersymmetry-invariant higher-order terms in the action: $`(\delta _0+_n\delta _n)(S_0+_nS_n)=0`$. Here $`\delta _0`$ is classical supersymmetry transformation, $`S_0`$ is classical action and $`\delta _n,S_n`$ are the quantum deformations of transformation law and higher derivative corrections to the classical action respectively. It is hardly possible to compute in a closed form the full derivative dependence of the effective action (see for review ). Therefore, what we can do in this situation is, relying on particular results known, to list all the supersymmetry invariants as well as deformed transformation rules with a given number of derivatives. For the well-known leading potential $`F^4`$ in the vector field strength sector (see for references and a recent progress ), the problem of constructing the full $`𝒩=4`$ superinvariant has been solved in .

There exist many different approaches that are used for the construction of the supersymmetric higher-derivative string effective action . In the instructive paper , an off-shell $`𝒩=3`$ supersymmetric extension of the Abelian $`D=4`$ BI action was constructed starting from the action of supersymmetric Maxwell theory in $`𝒩=3`$ harmonic superspace . The $`𝒩=3`$ superfield strength contains combination of the auxiliary field and the gauge field strength. The nonlinearity in the ordinary gauge field strength arises in its full form as the result of elimination of these auxiliary fields with the aid of their nonlinear equations of motion. This is different from $`𝒩=1,2`$, where superextensions of the bosonic BI action are fulfilled in each order of expansion in powers of the Maxwell field strength.

The purpose of the present paper is to consider a possible self-consistent way to find hypermultiplet dependent complements and the correspondent deformed hidden supersymmetry transformations, which are needed for manifestly $`𝒩=4`$ supersymmetric next-to-leading terms in the $`𝒩=4`$ SYM theory effective action. The while such an approach is useful for the resolution of existence problem concerning higher-derivative invariants, obtaining fully explicit expressions is extremely cumbersome due to the enormous number of terms and the problem of dealing witch partial integration.

## 2 On hidden $`𝒩=2`$ invariance of the $`F^6`$ term in $`𝒩=4`$ SYM theory

In order to construct subleading $`𝒩=2`$ hidden invariant terms in the derivative expansion of the effective action, which depend on all the $`𝒩=4`$ SYM multiplet fields, one can apply the Noether procedure with classical supersymmetry transformation modulo to boundary terms and free equations of motions ,

$`\delta _0𝒲={\displaystyle \frac{1}{2}}\overline{\epsilon }^{\dot{\alpha }a}\overline{D}_{\dot{\alpha }}^{}q_a^+,\delta _0\overline{𝒲}={\displaystyle \frac{1}{2}}\epsilon ^{\alpha a}D_\alpha ^{}q_a^+,`$ (1)
$`\delta _0q_a^\pm ={\displaystyle \frac{1}{4}}(\epsilon _a^\alpha D_\alpha ^\pm 𝒲+\overline{\epsilon }_a^{\dot{\alpha }}\overline{D}_{\dot{\alpha }}^\pm \overline{𝒲}).`$
Let us try to find a possible hypermultiplet completion for the following two-loop term $`F^6`$ found in  with the supergraph technique in the harmonic superspace:

$$\mathrm{\Gamma }_{(6|0)}=c_2d^{12}z[\frac{1}{\overline{𝒲}^2}\mathrm{ln}𝒲D^4\mathrm{ln}𝒲+\frac{1}{𝒲^2}\mathrm{ln}\overline{𝒲}\overline{D}^4\mathrm{ln}\overline{𝒲}]=c_2d^{12}z_{(6|0)}+c.c.,$$
(2)
where $`c_2=N^2g_{YM}^2\frac{1}{48(4\pi )^4}`$. This term consists of two different parts, marked by the fourth powers of distinct derivatives: $`D^4`$ or $`\overline{D}^4`$. This parts should be studied separately, because the variation rules (1) do not mix them. The variation of first part (2) induced by (1) can be written in the form

$$\delta _0_{(6|0)}=\frac{2q^{+a}}{\overline{𝒲}^3𝒲}(\overline{\epsilon }_a^{\dot{\alpha }}\overline{D}_{\dot{\alpha }}^{}\overline{𝒲}+\epsilon _a^\alpha D_\alpha ^{}𝒲)D^4\mathrm{ln}𝒲+\frac{q^{+a}\epsilon _a^\alpha D_\alpha ^{}𝒲}{𝒲\overline{𝒲}^3}D^4\mathrm{ln}𝒲.$$
(3)
The heart of all problems is the fact that the classical variation (3) of $`_{(6|0)}`$ generates terms non-symmetric under the replacement $`\epsilon \overline{\epsilon }`$.

Further we consider classical transformations $`\delta _0`$ defined by (1) along with their deformations $`\delta _{(n|k)}`$. The full deformed transformations are considered as an expansion in powers of $`D,\overline{D}`$ as well as in powers of $`X=\frac{2q^{+a}q_a^{}}{𝒲\overline{𝒲}}`$, i.e. $`\delta =\delta _0+\delta _1(D^4)+\delta _2(D^8)+\mathrm{}`$ The subscript $`k`$ in deformations indicates the power of $`X`$, e.g. $`\delta _1=\delta _{(1|k)}`$. Let us introduce the first complement to (2) in the form

$$_{(6|1)}=d_1[X\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲+X\frac{1}{𝒲^2}\overline{D}^4\mathrm{ln}\overline{𝒲}].$$
(4)
Its variation $`\delta _0^{(q)}_{(6|1)}`$ in $`q^\pm `$ cancel the first term in (3) if $`d_1=2`$, but the other part of (3) is not cancelled. Since the structure either of the two functionals (2) is not symmetric in respect to $`𝒲\overline{𝒲}`$, while the $`\delta ^{(q)}`$-variation is symmetric, the discrepancy between variations like (2) and (3) will always appear in each step of the variational procedure. To resolve this discrepancy, we consider the one-loop $`\mathrm{\Gamma }_{(4)}F^4`$ term along with its well-known hypercomplement :

$$\mathrm{\Gamma }_{(4)}=c_1d^{12}z[\mathrm{ln}𝒲\mathrm{ln}\overline{𝒲}+\frac{1}{2}X+\frac{1}{43}X^2+\mathrm{}],$$
(5)
where $`c_1=N\frac{1}{(4\pi )^2}`$. We know that this term will be renormalizable by neither the higher loop nor the instanton contributions. Suppose that the classical hidden supersymmetry is deformed as follows

$$\delta _{(1|0)}𝒲=\frac{\overline{A}}{2}\overline{\epsilon }^{\dot{\alpha }a}\overline{D}_{\dot{\alpha }}^{}q_a^+\frac{1}{𝒲^2}\overline{D}^4\mathrm{ln}\overline{𝒲},\delta _{(1|0)}\overline{𝒲}=\frac{A}{2}\epsilon ^{\alpha a}D_\alpha ^{}q_a^+\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲$$
(6)
$$\delta _{(1|0)}q_a^\pm =\frac{1}{4}[B\epsilon _a^\alpha D_\alpha ^\pm 𝒲\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲+\overline{B}\overline{\epsilon }_a^{\dot{\alpha }}\overline{D}_{\dot{\alpha }}^\pm \overline{𝒲}\frac{1}{𝒲^2}\overline{D}^4\mathrm{ln}\overline{𝒲}].$$
(7)
If the following conditions

$$c_2+c_1\frac{(\overline{A}\overline{B})}{2}=0,c_2+c_1\frac{(AB)}{2}=0$$
(8)
for the coefficients introduced in (6 \- 7) are satisfied, then deformed variations of the first two terms in $`\mathrm{\Gamma }_{(4)}`$ can cancel the last term in (3). The variation of the first complement $`\delta _0^{(𝒲)}_{(6|1)}`$ (4) in $`𝒲`$ under classical transformation rules (1) is

$$\delta _0^{(𝒲)}_{(6|1)}=4\frac{5}{3}\frac{q^+q^{}}{𝒲^2\overline{𝒲}^4}q^{+a}(\overline{\epsilon }_a^{\dot{\alpha }}\overline{D}_{\dot{\alpha }}^{}\overline{𝒲}+\epsilon _a^\alpha D_\alpha ^{}𝒲)D^4\mathrm{ln}𝒲$$
(9)
$$+\frac{4}{3}\frac{2q^+q^{}}{𝒲^2\overline{𝒲}^4}q^{+a}\epsilon _a^\alpha D_\alpha ^{}𝒲D^4\mathrm{ln}𝒲+\frac{4}{3}\frac{\overline{\epsilon }_a^{\dot{\alpha }}\overline{D}_{\dot{\alpha }}^+q^a}{𝒲^2\overline{𝒲}^3}(q^+D^+q^{})\frac{1}{16}D_\alpha ^+D^2\mathrm{ln}𝒲.$$

Now introduce the second complement

$$_{(6|2)}=d_2[X^2\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲+X^2\frac{1}{𝒲^2}\overline{D}^4\mathrm{ln}\overline{𝒲}].$$
(10)
Its variation in $`q^\pm `$ is exactly the first term in (9), and if we choose $`d_2=\frac{5}{3}`$, it will cancel the variations involving $`\overline{\epsilon }`$. At the same time, the rest of variation (9) is saved. In order to cancel its first part, one can consider variation of (5) under the following deformation of the transformations

$$\delta _{(1|1)}\overline{𝒲}=\frac{A_1}{2}X\epsilon ^{\alpha a}D_\alpha ^{}q_a^+\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲,\delta _{(1|1)}𝒲=\frac{\overline{A}_1}{2}X\overline{\epsilon }^{\dot{\alpha }a}\overline{D}_{\dot{\alpha }}^{}q_a^+\frac{1}{𝒲^2}\overline{D}^4\mathrm{ln}\overline{𝒲},$$
(11)
$$\delta _{(1|1)}q_a^{}=\frac{B_1}{4}X\epsilon _a^\alpha D_\alpha ^{}𝒲\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲+\frac{\overline{B}_1}{4}X\overline{\epsilon }_a^{\dot{\alpha }}\overline{D}_{\dot{\alpha }}^{}\overline{𝒲}\frac{1}{𝒲^2}\overline{D}^4\mathrm{ln}\overline{𝒲}.$$
(12)
Using these deformations, one can find variation of (5): the first term under $`\delta _{(1|1)}^{(\overline{𝒲})}`$, the second term under $`\delta _{(1|0)}^{(\overline{𝒲})}`$ as well as $`\delta _{(1|1)}^{(q)}`$, and the third term under $`\delta _{(1|0)}^{(q)}`$. The part of the variations, we are interested in, is

$$\delta _{(4)}=c_1\left[\frac{2}{3}A_1+B_1+\frac{BA}{3}\right](q^+q^{})q^{+a}\epsilon _a^\alpha D_\alpha ^{}𝒲\frac{D^4\mathrm{ln}𝒲}{𝒲^2\overline{𝒲}^4}.$$
(13)
The requirement of the cancellation between the second term in $`\mathrm{\Delta }_{(6|1)}`$ (9) and the corresponding term $`\delta _{(4)}`$ (13) gives the following equation

$$4\frac{2}{3}c_2=c_1(\frac{BA}{3}+B_1\frac{2}{3}A_1),\mathrm{or}c_1(B_1\frac{2}{3}A_1)=2c_2.$$
(14)
These relations define coefficients in the decomposition of the modified supersymmetry transformations in powers of $`X`$. There is arbitrariness in choosing $`A,B`$, and, therefore, some additional information is needed to get rid off this arbitrariness. For cancelling the second part of (9), we introduce a complement of the new type $`_{(6|1)}^{}=\frac{1}{3}X(\frac{q^+D^{+\alpha }q^{}}{𝒲\overline{𝒲}})\frac{1}{\overline{𝒲}^2}D_\alpha ^+D^2\mathrm{ln}𝒲`$.

The example presented above shows that, as a matter of principle, a complement for $`_{(6|0)}D^4`$ (2) is defined by the classical transformations (1) generated by $`\overline{\epsilon }`$, while all contradictions arising in the $`\epsilon `$ sector variations can be eliminated by the hidden transformation modifications $`\delta _{(1|n)}X^n`$ order by order. Thus, the the problem is split into two separate tasks. This is our main idea how to overcome the difficulty in constructing of the hidden $`𝒩=2`$ invariants with derivatives of the vector strength of the $`𝒩=2`$ multiplet.

Let’s consider a particular task of obtaining the leading term in the full $`F^6`$ effective action, which can be solved with the above considerations. For this purposes, it is sufficient to consider a generic term of the series of complements to (2) in the form

$$\mathrm{\Gamma }_{(6|n)}=d_nd^{12}z\left(\frac{2q^+q^{}}{𝒲\overline{𝒲}}\right)^n\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲+c.c.$$
(15)
The classical variation $`\delta _0_{(6|n)}`$, generated to the parameter $`\overline{\epsilon }`$, for terms $`D^4\mathrm{ln}𝒲`$ is

$$\left[d_nn\frac{(2q^+q^{})^{n1}}{𝒲^n\overline{𝒲}^{n+2}}+d_n\frac{n(n+4)}{n+2}\frac{(2q^+q^{})^n}{𝒲^{n+1}\overline{𝒲}^{n+3}}\right](q^{+a}\overline{ϵ}_a^{\dot{\alpha }}\overline{D}_{\dot{\alpha }}^{}\overline{𝒲})D^4\mathrm{ln}𝒲.$$
(16)
The requirement of cancellation variations of $`\delta _0_{(6|n)}`$ and $`\delta _0_{(6|n+1)}`$ is fulfilled if

$$d_n=d\frac{(n+2)(n+3)}{n}.$$
(17)
Summing all complements $`_{(6)}^q=_{n=0}^{\mathrm{}}_{(6|n)}(X)`$, one obtains

$$\mathrm{\Gamma }_{(6)}^q=\frac{c_2}{6}d^{12}z[\frac{X}{(1X)^2}+\frac{5X}{1X}6\mathrm{ln}(1X)]\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲.$$
(18)
Of course, this result is not entirely full in the sense that it should be completed by contributions containing hypermultiplet derivatives. Generic terms with $`q^+D^+q^{}`$ derivatives can be found. We introduce a new type of complement

$$_{(6|n)}^{}=p_n\frac{(2q^+q^{})^n}{𝒲^{n+1}\overline{𝒲}^{n+3}}q^+D_\alpha ^+q^{}D^{+\alpha }D^3\mathrm{ln}𝒲.$$
(19)
Then the requirement of the cancellation of its variation $`\delta _0^q(\overline{\epsilon })`$ with an appropriate term in variation (15) leads to a inhomogeneous recurrent relation. The relation has the solution

$$p_n=\frac{1}{6}\frac{(n+2)(n+3)}{n+1}\frac{1}{12}\frac{(n+2)^2(n+3)}{n+7}H_{n+1}^{(2)},$$
(20)
where $`H_n^{(2)}=_{k=1}^n\frac{1}{k^6}`$ is the harmonic number $`[n,2]`$. For the coefficients of the next generic term of the series $`_{(6|n)}^{\prime \prime }=h_nX^n(q^+D^+q^{})^2D^2\mathrm{ln}𝒲`$, there must be the recurrent relation on $`p_n,d_n,h_n`$.

We see that for $`𝒩=4`$ supersymmetrization of next-to-leading terms, one has to consider transformations that mix different terms in the derivative expansion of the effective action! To make sure that the guessed transformations (6-7) are not unreasonable, one can consider a variation of the well-known classical action . The variation of the hypermultiplet action is proportional to the on-shell equation of motion $`D^{++}q^+=0`$, but variation of the vector strength action in (6) is

$$\delta _{(1|0)}\mathrm{\Gamma }_0=\frac{1}{8}d^{12}z\{\overline{A}\overline{\epsilon }^{\dot{\alpha }a}\overline{D}_{\dot{\alpha }}^{}q_a^+\frac{1}{𝒲}\mathrm{ln}\overline{𝒲}+A\epsilon ^{\alpha a}D_\alpha ^{}q_a^+\frac{1}{\overline{𝒲}}\mathrm{ln}𝒲\}.$$
(21)
This expression has the same structure as $`\delta _0_{(4|0)}`$. The fact that in order to cancel (21) one should take into account classical variations of $`_{(4|0)}\mathrm{ln}𝒲\mathrm{ln}\overline{𝒲}`$ serves as an additional cross-checking of the consistency of the proposed recipe.

It is also interesting to consider the first deformed variations of $`F^6`$ (2). We obtain

$$\delta _{(1|0)}_{(6|0)}=c_2[\frac{\overline{A}}{𝒲^3\overline{𝒲}^2}\overline{\epsilon }^{\dot{\alpha }a}\overline{D}_{\dot{\alpha }}^{}q_a^+D^4\mathrm{ln}𝒲\overline{D}^4\mathrm{ln}\overline{𝒲}\frac{A}{\overline{𝒲}^5}\epsilon ^{\alpha a}D_\alpha ^{}q_a^+\mathrm{ln}𝒲(D^4\mathrm{ln}𝒲)^2].$$
(22)
The first term in the brackets looks like classical variation of the one-loop $`F^8`$ structure $`_{(8|0)}=𝚿^2\overline{𝚿}^2=\frac{1}{𝒲^2}\overline{D}^2\mathrm{ln}\overline{𝒲}\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲`$ but with coefficient $`c_2`$. This means that the one-loop coefficient $`\frac{1}{2(24\pi )^2}`$ in front of $`F^8`$ structure , of the effective action will be renormalized by two-loop contributions. The second term in (22) is a term of new type. For its cancellation, one needs to add another structure

$$\mathrm{\Gamma }_{(8|0)}^{}=c_8d^{12}z\mathrm{ln}𝒲\left(\frac{1}{\overline{𝒲}^2}D^4\mathrm{ln}𝒲\right)^2,$$
(23)
with typical $`\frac{1}{(4\pi )^4}`$ two-loop coefficient $`c_8=c_2\frac{A}{2}`$. This allows us to conjecture that accurate two-loop calculations should give such $`F^8`$-structure in the effective action.

Thus, the self-consistent obtaining of the appropriate hypermultiplet dependent contributions in the effective action and modification of the hidden supersymmetry transformation allows one to obtain information about renormalizable terms as well as non-renormalizable higher corrections terms.

## 3 Acknowledgements

The authors thank the Organizing Committee of the SQS’03 conference for warm welcome and partial support. The work was supported in part by INTAS grant, INTAS-00-00254 and RFBR grant, project No 03-02-16193. The authors also would like to thank I.L.Buchbinder and E.A. Ivanov for numerous discussion and a critical remarks.


format: text

text:
# Reverse Monte Carlo modeling of amorphous silicon

## I introduction

The structure of amorphous semiconductors is well represented by continuous random network (CRN) model introduced by Zachariasen zach  70 years ago. The CRN model has the simplicity that each of the atoms should satisfy its local bonding requirements and should have as small strain as possible in the network which is generally characterized by having a narrow bond angle as well as bond length distribution. In spite of its apparent simplicity, the structural modeling of high quality tetrahedral amorphous semiconductors appears to be quite difficult. There have been many models of amorphous silicon polk ; gutt ; woot ; bark ; Car ; Thorpe ; Stillinger ; Kaxiras ; Keating ; Holender  proposed in the last 30 years which include from very simple hand-built model of Polk polk , computer generated periodic network model of Guttman gutt  to the complex model of Wooten, Winer and Weaire (WWW) woot , but most of these models have some limitations in one way or another in describing the true nature of the amorphous state. The last method, the so called sillium approach of WWW, is based on the strategy of randomizing and relaxing the network is so far the most successful method of producing minimally strained CRN. The algorithm, in its modified form developed by Djordjevi$`\stackrel{´}{c}`$ et. al. Thorpe  and Barkema and Mousseau bark , can produce a CRN which is comparable to experimental results and is capable of producing a clean band gap without any defect states in the gap.

In this paper we develop a different approach to model amorphous semiconductors known as reverse monte carlo (RMC) simulation gere ; laszlo ; walt ; rmc1 ; rmc2 ; rmc3 ; rmc4 . Our primary objective is to produce structural configurations that are consistent with experimental data but at the same time we go one step further to generate realistic configurations for comparison with models obtained via other routes. We emphasize that producing realistic models (meaning models which agree with all experiments) requires more than spatial pair correlations, and identify additional constraints which lead to realistic models.

The existing RMC models of amorphous semiconductors are found to be inadequate and fail to produce some of the basic experimental features of amorphous tetrahedral semiconductors. Gereben & Pusztai gere ; laszlo  have carried out RMC simulation of tetrahedral semiconductors using a number of models ranging from completely disordered configuration to randomized diamond structure. Although a certain degree of tetrahedral character in the bond angle distribution was reflected in their work, most of the models show an unphysical peak in bond angle distribution around $`60^{}`$. The work of Walters and Newport walt  on amorphous germanium made some progress toward getting the correct bond angle distribution, but the number of 3-fold coordinated atoms are quite high in their model and in absence of any discussion on local strain and electronic properties it is difficult to say how reliable their models are when it comes looking at the electronic properties.

A developing area where RMC may be applied successfully is for modeling amorphous materials exhibiting medium-range order (MRO). Such MRO is characterized by the existence of $`1020`$Å scale structure. Recent developments in fluctuation electron microscopy (FEM)fem  and its application on amorphous germanium and silicon have indicated that computer generated CRN model of these materials lack the characteristic signature of MRO. Since RMC is based on experimental data, it provides a promising scheme to model amorphous materials having medium-range order by including the experimentally measured FEM signal as input data to augment pair correlations.

The plan of the paper is as follows. In section II, we briefly mention the basic philosophy of reverse monte carlo modeling and some of its salient features. This is followed by role of constraints in RMC modeling in section III where we illustrate how a set of judiciously chosen constraints can be used to construct a reliable model of amorphous silicon. Finally we compare our results with those obtained from earlier RMC models and a model obtained via WWW algorithm.

## II Basics of RMC

The RMC method has been described in detail elsewhere rmc4 . Here we briefly outline the basic philosophy of RMC. At the very basic level, RMC is a technique for generating structural configurations based on experimental data. The logic is very appealing: any model of a complex material worthy of further study should, at a minimum, agree with what is known (that is, the experiments). By construction, the RMC scheme enforces this (and for contrast, a molecular dynamics simulation may not). In an ideal implementation, one should find a model agreeing with all known information, but this is not easy to accomplish, though we make some progress below. The approach was originally developed by McGreevy & Pusztai rmc4 ; gere  for liquid and glassy materials for lack of different routes to explore experimental data but in recent years progress has been made toward modeling crystalline systems as well keen . Starting with a suitable configuration, atoms are displaced randomly using the periodic boundary condition until the input experimental data (either the structure factor or the radial distribution function) match with the data obtained from the generated configuration. This is achieved by minimizing a cost function which consists of either structure factor or radial distribution function along with some appropriately chosen constraints to restrict the search space. Consider a system having $`N`$ number of atoms with periodic boundary condition. One can construct a generalized cost function for an arbitrary configuration by writing :

$$\xi =\underset{j=1}{\overset{K}{}}\underset{i=1}{\overset{M}{}}\eta _i^j\{F_E^j(Q_i)F_c^j(Q_i)\}^2+\underset{l=1}{\overset{L}{}}\lambda _lP_l$$
(1)
where $`\eta _i^j`$ is related to the uncertainty associated with the determination of experimental data points as well as the relative weight factor for each set of different experimental data. The quantity $`Q`$ is the appropriate generalized variable associated with experimental data $`F(Q)`$ and $`P_l`$ is the penalty function associated with each constraint. For example, in case of radial distribution function and structure factor, $`Q`$ has the dimension of length and inverse length respectively. In order to avoid the atoms getting too close to each other, a certain cut-off distance is also imposed which is typically of the order of interatomic spacing. In RMC modeling, this is usually obtained from the radial distribution function by Fourier transform of the measured structure factor. This is equivalent to adding a hard sphere potential cut-off in the system which prevents the catastrophic build up of potential energy.

In spite of the fact that RMC has been applied to many different types of systems – liquid, glasses, polymer and magnetic materials, questions are often raised about the reliability of results obtained from RMC simulation. The method has never been accepted without some degree of controversy and the most popular criticism is the lack of unique solution from RMC. RMC can produce multiple configurations having the same pair correlation function. This lack of uniqueness, however, is not surprising, since usually only the pair correlation function or structure factor is used in modeling the structure, while there exists an infinite hierarchy of higher order correlation functions carrying independent structural information are neglected. In other words, RMC samples from the space of all models consistent with some limited body of data – in its simplest form (analyzing a single experiment) RMC is an ideal gauge of how non-specific the data is with respect to identification of an atomistic model. If the modeler possesses a priori information independent of that implicit in the experiment being fit to, it is necessary to add this information to the modeling in some fashion to receive a model in joint agreement with the experiment and the additional information.

## III A new RMC model

We begin by including the minimal information that is necessary to model a configuration of $`a`$-Si. In so doing, we use the radial distribution function obtained from a high quality model of amorphous silicon. This latter model was generated by Barkema and Mousseau bark  using a modified form of WWW algorithm woot  having bond angle distribution close to $`10^{}`$ with 100% 4-fold coordination. In addition to this RDF, we also impose the conditions that the average bond angle of all the triplets Si-Si-Si should be near $`109.5^{}`$ and the corresponding root mean square deviation should be no less than $`10^{}`$. The number of 4-fold coordinated atoms is driven to a specified value during the simulation by including a constraint on the average coordination number. It is to be noted that while there is no limit to the number of constraints that can be included in the system, there is no guarantee that mere inclusion of more constraints will necessarily give better results. Forcing a completely random configuration with too many competing constraints may cause the configuration to be trapped in the local minimum of the function $`\xi `$ and may prevent the system from exploring a large part of the search space. By adding only the essential constraints that describe the chemical and geometrical nature of the bonding correctly, Eq. 1 can be written as :

$`\xi `$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{M}{}}}\lambda _1\{F_E(x_i)F_c(x_i)\}^2`$ (2)
$`+`$ $`\lambda _2(\theta _0\theta )^2`$
$`+`$ $`\lambda _3(\delta \theta _0\delta \theta )^2`$
$`+`$ $`\lambda _4\{1\mathrm{\Theta }(xx_c)\}`$
$`+`$ $`\lambda _5\{\varphi _0\varphi \}^2`$
where,

$`\theta `$ $`=`$ $`{\displaystyle \frac{1}{N_\theta }}{\displaystyle \underset{i\{j,k\}}{}}\theta _{ijk}`$
$`\delta \theta `$ $`=`$ $`\sqrt{(\theta \theta _0)^2}`$
In Eq. 2, $`\theta `$ and $`\delta \theta `$ are the average angle and the rms deviation while $`\varphi `$ and $`\varphi _0`$ are the current and proposed concentration of the 4-fold coordinated atoms. It is important to note that each of the terms in Eq. 2 is non-negative, and should decrease ideally to zero during the course of minimization.

Since the cost of energy associated with the bond length relaxation is more than the bond angle relaxation, atomic arrangements with large bond angle distribution but having correct RDF frequently result. The coefficients, $`\lambda _1`$ to $`\lambda _3`$, for the different terms in Eq. 2 can be chosen appropriately to minimize this effect. In general the coefficients $`\lambda `$ are constant during the course of simulation but the minimization procedure can be slightly accelerated by making them vary in such a way that the contribution from each of the term are of the same order during the course of simulation. The coefficient $`\lambda _4`$ is usually assigned a large value in order to include a hard sphere cut-off as mentioned earlier so that no two particles can come closer to $`x_c`$ while the coefficient $`\lambda _5`$ maintains the number of 4-fold coordinated atoms to a specified value. In RMC simulation of amorphous tetrahedral semiconductors one usually encounters the problem of having a pronounced peak at $`60^{}`$. This peak is a characteristic feature of unconstrained RMC simulation and is due to the formation of equilateral triangles by three atoms. In the work of Gereben and Pusztai gere ; laszlo , attempts were made to overcome this difficulty by constraining the bond angle distribution as well as by making an initial configuration which is 100% 4-fold obtained from a diamond lattice. The resulting structure is, however, found to be unstable and on relaxation using a suitable potential, the configuration tends to get back toward the starting structure, i.e., randomized diamond in this case note1 .

In the approach of Walter and Newport walt , the initial random configuration was examined and any “triples”, i.e. , three atoms forming an equilateral triangle was removed before the beginning of RMC fit. By selective removal of such unwanted triplets, they have been able to generate configuration of $`a`$-Ge without having a peak at $`60^{}`$. The approach that we have taken in our work is more general and starts with a completely random configuration. This eliminates, in the first place, any possible local ordering that may exist in the starting structure (e.g., randomized diamond structure retains the memory of tetrahedral ordering). Furthermore, we have not included or excluded any special configuration in our starting structure, e.g., three atoms forming an equilateral triangle. Based on experimental consideration, we have included only the key features of amorphous tetrahedral semiconductors – an average bond angle of $`109.5^{}`$ having rms deviation of $`10^{}`$ which is consistent with the RDF obtained from a WWW relaxed model used in our calculation. For the 500-atom model reported in this work, we have chosen a cubic box of length 21.18Å which corresponds to number density 0.0526 atom/Å3. The initial configuration is generated randomly so that no two atoms can come closer to 2.0Å. The configuration is then relaxed by moving the atoms to minimize the cost function $`\xi `$. In addition to applying standard monte carlo moves in which a single or a group of atoms is randomly displaced, a variety of monte carlo moves have been implemented in our work. For example, in one of such moves, a 3-fold or 5-fold atom is selected and the nearest neighbor distance is examined. If the distance is greater than 2.7Å, the neighboring atom is displaced in order to bring the distance within a radius of 2.7Å.

The maximum displacement of a monte carlo move is limited to 0.2-0.4Å throughout the simulation. Since we are interested in the electronic structure as well, we confine ourselves within a reasonable system size for studying the generated structure using a first principles density functional Hamiltonian. The density functional calculations were performed within the local density approximation (LDA) using the local basis first principles code Siesta siesta1 . We have used a non self-consistent version of density functional theory based on the linearization of the Kohn-Sham equation by Harris functional approximation Harris  along with the parameterization of Perdew and Zunger Perdew  for the exchange-correlation functional.

## IV Results

The results for the model including all the constraints are presented in Figs. 1-4. Since the structure factor is generally considered to be more sensitive to an arbitrary small change in the atomic positions than the radial distribution function, we have plotted the structure factors for the constrained RMC and WWW model in Fig. 1. It is evident from Fig. 1 that the agreement between the RMC and WWW model is very good both for small and large values of Q. In order to further justify the credibility of our model, we have plotted in Fig. 2 the structure factor from the experiment of Laaziri et. allaaziri  along with the same obtained from our RMC model. Once again we find that the agreement between the structure factor from RMC and the experimental results is quite good except for the few points near the first peak. It is very tempting to think this deviation as a finite size effect coming from the finiteness of our model. We have therefore calculated the structure factor for WWW models containing 300 to 4096 atoms of Si but the deviation continues to remain. Holender and Morgan Holender  also observed similar deviation near the first peak in their work with a much larger model containing 13824 atoms which was compared with the experimental data obtained by Fortner and Lannin Fortner .

In Fig. 3, we have plotted the bond angle distributions (BADF) for both the RMC and WWW model. As we have discussed in section II, the radial distribution function or structure factor can not alone provide all the necessary information that are needed to characterize an atomic configuration obtained from a reverse monte carlo simulation. A further characterization beyond pair correlation function is therefore vital and necessitates the need for getting some idea about the 3-body correlation function. It is clear from the Fig. 3 that the distribution obtained from the RMC model follows the tetrahedral character observed in amorphous semiconductors. The average bond angle in this case is found to be $`109.01^{}`$ with rms deviation of $`12.5^{}`$. An important aspect of the bond angle distribution in Fig. 3 is that most of the angles are lying between $`70^{}`$-$`150^{}`$ compared to $`80^{}`$-$`140^{}`$ in WWW case. We emphasize at this point that the earlier works on modeling amorphous tetrahedral semiconductors using RMC predicted a much wider bond angle distributions. Gereben and Pusztai gere  have observed a pronounced, unphysical peak at $`60^{}`$ except for the model starting with diamond structure while Walters and Newport walt  have reported a bond angle distribution of $`a`$-Ge which is as wide as $`60^{}`$-$`180^{}`$. It is an important development here that by adding three more constraints ($`\lambda _2,\lambda _3`$, and $`\lambda _4`$) we have achieved a significantly improved results. Both the radial and the bond angle distribution functions reported here are at par with the results obtained from molecular dynamics simulation and is comparable to those obtained from WWW model. The fact that the inclusion of these two constraints leads to a significant improvement is not surprising. For a large continuous random network (CRN) model of amorphous tetrahedral semiconductor, one can approximate the bond angle distribution as nearly Gaussian note3 . This approximated Gaussian distribution can defined by the first two moments of the distribution function. By specifying these two moments as constraints in Eq. 2, we correctly describe the tetrahedral bonding geometry of the atoms which along with the radial distribution function produces a configuration more realistic than those obtained from models based on RDF or structure factor only. This suggests that in addition to the radial distribution of the atoms, one needs to include some relevant information about the nature of 3-body correlation among the atoms to construct a realistic configuration.

Having studied the radial and bond angle distribution we now address the electronic density of states calculations. While the width of the bond angle distribution function (BADF) and the structure factor together indeed gives some idea about the quality of the model, some of the features e.g., the existence of spectral gaps and the position of defects states in the spectrum can be studied by looking at the electronic density of states only. The structure obtained from RMC simulation is first relaxed using the density functional code Siesta and is found to be close to an energy minimum in the local density approximation (LDA). This is an important test for determining the stability of the structure obtained from RMC simulation and as far as we are concerned almost all earlier works on RMC have completely neglected this issue. In Fig. 4, we have plotted the electronic density of states (EDOS) for the constrained model. The EDOS appears with all the characteristic features of $`a`$-Si with the exception of a clean gap in the spectrum. This behavior is not unexpected in view of the fact that 88% of the total atoms are found to be 4-fold coordinated with an average coordination number 3.85. The presence of the defect states makes the gap noisy and at the same time the use of LDA underestimates the size of the gap. This EDOS is in significantly better agreement with optical measurements than conventional RMC models with much higher defect concentrations and spurious bond angles. It is interesting to observe that the average coordination number from our model is very close to the experimental value of 3.88 reported by Laaziri et. al. laaziri .

## V conclusions

We have presented a model of amorphous silicon based on reverse monte carlo simulation. One of the novel features of our model is to start with a completely random structure and then to relax toward a realistic configuration by adding a number of physically relevant constraints. The characteristic features of the tetrahedral bonding are taken into account by adding constraints on average bond angle and its deviation from the mean while the number of 4-fold coordinated atoms is maintained at a specified value by further use of a constraint on average coordination number. The radial and the bond angle distribution obtained from our model is found to be in excellent agreement with a high quality CRN model produced by WWW algorithm. We have also compared the structure factor with the experimental data obtained by Laaziri et al. and observed a reasonably good agreement. By relaxing the model using the first principles density function code Siesta, we find that the model is close to the energy minimum for LDA and is stable. The electronic density of states (EDOS) obtained from our model contains all the essential feature of amorphous silicon including a signature of the band gap. Although the model does not produce a clean gap in the spectrum, the quality of the EDOS is at par with models obtained from molecular dynamics simulation. Our RMC algorithm presents a significant improvement on previous RMC studies and makes it possible to compare for the first time, albeit qualitatively, the structural and electronic properties of RMC models with its WWW counterpart. We expect that further developments toward this direction will eventually make RMC as an useful modeling tool incorporating experimental information and can be used effectively without any criticisms in modeling complex materials.

###### Acknowledgements.

We thank John Abelson (UIUC) for motivating us in this work. We thank Gerard Barkema (ITP, Utrecht) for providing us with WWW models of amorphous silicon. We thank Stephen Elliott (Trinity College, Cambridge) for very valuable discussions. We acknowledge the support of National Science Foundation under grant Nos. DMR-0310933 and DMR-0205858.


text:
# A Very Large Array 3.6 cm continuum survey of Galactic Wolf-Rayet stars

## 1 Introduction

Wolf-Rayet (WR) stars are the evolved descendants of massive O-type progenitors ($`M_\mathrm{i}`$$`>`$ 22 M) and are considered to be the immediate precursors of supernovae Type Ib and Type Ic (e.g., Maeder 1981; Meynet & Maeder 2003). They are characterized by powerful stellar winds with terminal velocities $`v_{\mathrm{}}`$$``$ 400 - 5500 km s-1 (e.g., Prinja et al. 1990; Eenens & Williams 1994; Kingsburgh et al. 1995; for an overview see van der Hucht 2001, table 15, column 14), and mass loss rates in the range (1 - 5)$`\times `$10-5 M yr-1 (e.g., Abbott et al. 1986 \[hereafter AB86\]; Leitherer et al. 1995, 1997 \[hereafter LC95, LC97\]; Chapman et al. 1999 \[hereafter CL99\]). Due to these properties WR stars are not only interesting in themselves, but also important species for chemical enrichment and galactic evolution (e.g., Dray et al. 2003; De Donder & Vanbeveren 2003). For general reviews of the WR phenomenon, see, e.g., van der Hucht (1992) and Maeder & Conti (1994).

The strong stellar winds of WR stars reveal themselves in the electromagnetic spectrum from the UV (P-Cygni profiles of resonance lines) to the IR/radio range (free-free radiation). Therefore, an important diagnostic tool to obtain information about stellar wind properties of massive stars is the study of the radio continuum emission from their close environs. Radio observations provide a reliable method for determining stellar mass loss rates if the emission is thermal in origin, since the measured radio flux density is related to the mass loss rate by an analytical expression (e.g., Wright & Barlow 1975; Panagia & Felli 1975). This emission is the free-free radiation originating in the ionized wind, expanding at constant velocity. It arises from the outer parts of the wind envelope and is partially optically thick in the radio region, showing positive spectral indices ($`\alpha `$$``$ +0.6, where $`S_\nu `$$``$$`\nu ^\alpha `$).

However, the last decade has seen an increasing number of discoveries of WR stars with non-thermal radio emission in addition to their free-free emission. Non-thermal emission is thought to be produced by electrons accelerated in strong shocks in the winds of single stars (White 1985) or in the shocked wind-wind collision region between massive WR+O binaries (e.g., Williams et al. 1990; Eichler & Usov 1993). Unlike thermal emission, non-thermal emission has a negative spectral index and can be variable. Van der Hucht et al. (1992) suggested that all cases of WR non-thermal radio emission correspond to WR+OB colliding wind binaries, which has been corroborated by Dougherty & Williams (2000).

Due to their relatively large distances (van der Hucht 2001), the number of WR radio continuum detections was limited to only about 44 cases (AB86, LC95, LC97 and CL99). The VLA survey by AB86, which also quotes previous observations, contains radio data for 40 WR stars at $`\delta `$$`>`$$`47`$ with 27 detections using the VLA, with rms noise values of 0.06 - 0.23 mJy beam-1. They found 77% probable or definite thermal wind sources and 12% probable or definite non-thermal sources. In the southern hemisphere, LC95, LC97 and CL99 surveyed 36 WR stars ($`\delta `$$`<`$$`45`$) at two frequencies using the ATCA, with rms noise values of 0.10 - 0.15 mJy beam-1. Other WR radio continuum studies have been performed for WR 140 by Williams et al. (1987, 1990, 1994) and White & Becker (1995); for WR 146 by Dougherty et al. (1996, 2000) and Setia Gunawan et al. (2000); for WR 147 by Moran et al. (1989), Churchwell et al. (1992), Davis et al. (1996), Williams et al. (1997), Skinner et al. (1999), Setia Gunawan et al. (2001), Watson et al. (2002) and Dougherty et al. (2003); for Galactic Center WR stars by Lang et al. (2001) and Lang (2003); for WR 98a, WR 104 and WR 112 by Monnier et al. (2002); and for massive stars in Cyg OB2 by Setia Gunawan et al. (2003b). In these studies, some 70 out of the $``$ 250 known Galactic WR stars (van der Hucht 2001, 2003) have been investigated for radio continuum emission, resulting into 44 positive detections. The ATCA survey also indicates that at least 40% of the WR stars with measured radio spectral indices contain a non-thermal component at centimeter wavelengths, a percentage that, coincidentally, is equal to the percentage of WR+OB binaries among the known Galactic WR stars (van der Hucht 2001).

In this paper we report on the results of a survey of radio continuum emission from Galactic WR stars at 3.6 cm with higher sensitivity than available before. Our aim is to increase the sample of WR stars with radio continuum observations in order to derive mass loss rates or upper limits, and to identify non-thermal cases. As a first step, we performed observations at one frequency only. Consequently, spectral indices of detected sources can only be determined here using radio data at other wavelengths from the literature and may be influenced by variability. Thus, new radio continuum observations are important to find new thermal and non-thermal WR radio sources, the thermal ones providing mass loss rates, crucial for testing existing wind models, and the non-thermal ones pointing to colliding winds in WR+OB binaries.

## 2 Selected Wolf-Rayet stars

In order to achieve these goals, we selected 33 WR stars from the catalog by van der Hucht (2001). To select the targets, we took into account all WR stars north of $`\delta `$ =$`46`$, and selected the ones with expected free-free flux densities $`S_{\mathrm{exp}}`$$`>`$ 0.075 mJy. $`S_{\mathrm{exp}}`$ was determined with the classical expressions by Panagia & Felli (1975) and Wright & Barlow (1975) for the radio emission from an optically thick stellar wind flowing at constant velocity. In order to estimate these values, we adopt for all sources the same mass loss rate $`\dot{M}`$ = (1.5 - 2.0)$`\times `$10-5 M yr-1, and stellar distances and terminal velocities from van der Hucht (2001).

With these assumptions, we found $`S_{\mathrm{exp}}`$$`>`$ 0.075 mJy for 64 WR stars. Only 27 of these stars had previous radio detections, and radio flux density upper limits were available for 12 stars. We were then left with 37 stars which either had not been observed yet or had not been detected. We selected 26 of these stars as our targets, listed in Table 1. In our sample we included some WR stars previously detected: WR 79a (HD 152804, detected with the VLA by Bieging et al. 1989), WR 81, WR 89, WR 98a, WR 104, WR 105, and WR 125. The detection of WR 81 by AB86 was doubtful. WR 98a, WR 104, WR 105 and WR 125 are known non-thermal emitters and, thus, may vary in flux density.

————–

TABLE 1

————–

In order to find new non-thermal candidates, we included binaries and suspected binaries listed by van der Hucht et al. (2001). About 70% of the stars are closer than 3.0 kpc, while no stars beyond 5.3 kpc are included.

## 3 VLA observations

We performed radio continuum observations at 8.46 GHz with the VLA in the direction of the selected WR stars. Thirty-three fields were observed on 15 September 2001 (during the move from the C- to the DnC-array), 5 and 8 October 2001 (in the DnC-array) and 12 November 2001 (in the D-array). With the exception of WR 89, the field centers coincide with the stellar positions. The bandwidth was 50 MHz. The on-source observing time was about 25 minutes for each field to achieve a rms noise of $``$0.04 mJy beam-1. The synthesized beam is about 6$`\stackrel{}{\mathrm{.}}`$1$`\times `$9$`\stackrel{}{\mathrm{.}}`$0 in all cases.

The data were edited, calibrated and imaged following a standard way using aips tasks. After editing, the data were calibrated both in amplitude and phase. The primary flux density calibrators were 3C 48 and 3C 286. We used phase referencing taking into account a selection of nearby point sources with precision positions (0$`\stackrel{}{\mathrm{.}}`$1 or better). The secondary calibrators were 0244+624, 0747–331, 0828–375, 1001–446, 1604–446, 1744–312, 1751–253, 1820–254, 1832–105, 1851+005, 1924+156, 2015+371, 2230+697 and 2250+558. The fields corresponding to WR 4, WR 96 and WR 142 were also self-calibrated in phase, since peak flux densities within the primary beam exceed 10 mJy beam-1. Primary beam correction was applied to the images corresponding to WR 89 and WR 105. WR 89 is in the same field as WR 87, 1$`\stackrel{}{\mathrm{.}}`$9 far from the field center. The primary beam correction factor for WR 89 is 1.5. In the case of WR 105, two additional point-like sources, with no relation to the WR star, are detected close to the center of the field. Diffuse emission from extended sources detected in the fields of WR 87, WR 104, WR 106, WR 113, WR 114, WR 121, WR 125, WR 142, WR 143 and WR 153ab was removed by applying an uv-range $``$ 5K$`\lambda `$ (rejecting spatial scales $``$ 40″).

## 4 Observed radio emission

We report 15 definite and 5 probable detections at 3.6 cm out of the 34 observed WR stars. For 10 out of the 15 definite detections and 3 out of the 5 probable detections this represents a first radio detection, while four definite detections had been previously identified at 3.6 cm and one had been detected previously at a different frequency. All 20 detections are unresolved. The images of the certain and probable detections are displayed in Figs. 1 to 20. The synthesized beam of each field is shown in a corner of the contour plots.

Radio detections with signal-to-noise S / N $``$ 5, as derived from the Gaussian fitting, and close agreement between radio and optical positions were classified as definite detections, while probable detections correspond to sources with S / N $``$ 4 (WR 81, WR 95, WR 114 and WR 155) or with relatively large difference between the optical and the radio positions (WR 79a and WR 95). Detections classified as probable need further confirmation.

————–

FIGURES 1 - 20

————–

Flux densities and radio positions were obtained by fitting a 2-D Gaussian with the aips-jmfit tool. In most of the detections, the uncertainty in the radio position is less than 1″. Except for WR 79a and WR 95, the optical and radio positions agree within 2″. The flux density uncertainties quoted in the table are equal to the rms noise (1$`\sigma `$) of the images. Flux density uncertainties were derived from the images by taking into account a relatively large region near the detected sources. For the 14 undetected sources, the upper limits in flux density correspond to 3 $`\sigma `$. Typically, the 1 $`\sigma `$ level is in the range 0.025 - 0.084 mJy beam-1. In the cases of WR 79a and WR 89, the 1 $`\sigma `$ level is 0.13 mJy beam-1. The large rms noise associated with the detection of WR 89 is due to the fact that this star is displaced about 1$`\stackrel{}{\mathrm{.}}`$9 from the field center. The typical uncertainty in the radio position is 0.5″.

Table 2 lists the observed 3.6 cm flux densities and radio coordinates of the definite and probable detections, along with the optical coordinates from the catalog of van der Hucht (2001). Non-detected sources are also indicated at the botton of the table. The latter have typical uncertainties of about 1″ (Leitherer et al. 1997).

————–

TABLE 2

————–

In Table 3, we compare our 3.6 cm detections with previously determined flux densities or upper limits at six radio wavelengths. The table also gives spectral indices $`\alpha `$ when calculable.

Since we performed observations at one frequency only, spectral indices can not be derived from our data. Thus, the nature of the emission was inferred by comparing the observed fluxes with earlier determinations, allowing us to identify variable sources. Variability in flux density, which is a hint for the presence of non-thermal emission, is also indicated in Table 3.

————–

TABLE 3

————–

## 5 Comments on individual stars

WR 79a, WN9ha (Fig. 4). This object is a known visual binary with $`\mathrm{\Delta }\varphi `$ = 6$`\stackrel{}{\mathrm{.}}`$9 (Mason et al. 1998; van der Hucht 2001, table 22). The radio source extends to the W. The earliest flux density determinations by BA89 (VLA) indicated thermal emission with a spectral index $`\alpha `$ = $`+`$0.8. Recent ATCA observations (November 2000) by Setia Gunawan et al. (2003a) show a flat non-thermal energy distribution with spectral index $`\alpha `$ = $``$0.0 $`\pm `$ 0.1. Our 3.6 cm observation, less than a year later, is consistent with the ATCA observation and, with respect to BA89 observation, with a variable non-thermal component.

WR 81, WC9 (Fig. 5). Our VLA image shows a source with a complex structure coincident with the optical position of the star. We identified the source to the N as the WR star because of its proximity (2$`\stackrel{}{\mathrm{.}}`$5) to the optical position. The radio source detected by AB86 and identified by these authors as the WR star (see their Fig. 1) is 4$`\stackrel{}{\mathrm{.}}`$6 south of the optical position and closer to the southern extension of the radio source detected in the new image. A spectral index cannot be derived.

WR 89, WN8h+OB (Fig. 7). This source was detected in the image corresponding to WR 87, 1$`\stackrel{}{\mathrm{.}}`$9 displaced from the field center. Previous flux density determinations by LC95 were consistent with a thermal wind. Judging from our observation, the 3.6 cm flux density dropped some 30 % in recent years. WR 89 is associated with the open cluster Havlen-Moffat No. 1 (C1715$``$387, Vázquez & Baume 2001).

WR 98a, WC8-9vd+? (Fig. 9). In the infrared this object shows a expanding dust spiral, rotating with a period $`P`$ = 565 $`\pm `$ 50 d (Monnier et al. 1999), which is considered to be the binary period. Monnier et al. (2002) observed this object with the VLA in 1999.7 and 2000.2 and found the source to be non-thermal. Our 3.6 cm flux density determination is $``$ 20 % below their result. Apparently, the non-thermal source is variable.

WR 104, WC9d+B0.5V (Fig. 12). In the infrared this object shows a expanding dust spiral, rotating with a period $`P`$ = 243 $`\pm `$ 3 d (Tuthill et al. 1999), which is considered to be the binary period. Monnier et al. (2002) found the source to be non-thermal. Our 3.6 cm flux density determination is $``$ 40 % below theirs, thus, the non-thermal source is apparently variable.

WR 105, WN9h (Fig. 13). This star shows variable non-thermal flux densities. Our 3.6 cm flux density is $``$ 50 % larger than the one measured by LC97. The double radio source to the S of the optical position of WR 105 was also detected by AB86 (see their fig. 1). These radio sources were not detected in the NVSS (Condon et al. 1998), nor were optical counterparts found. Zoonematkermani et al. (1990, radio continuum at 20 cm) detected the double source with flux densities of 22 mJy and 16 mJy, while from our VLA data we estimate 6.4 $`\pm `$ 0.2 mJy and 5.4 $`\pm `$ 0.2 mJy. These non-thermal sources have spectral indices $`\alpha `$ of $``$0.7 and $``$0.6, respectively.

WR 125, WC7ed+O9III (Fig. 17). The AB86 1982 observations at 6.3 cm and 2.0 cm showed a non-thermal radio source with $`\alpha `$ = $``$0.5. Our 2001 3.6 cm flux density observation agrees with the AB86 data for a non-thermal source with that spectral index. The source is variable in flux and turned from non-thermal to thermal in WH92 (see Table 3). We suggest here that the emission is again non-thermal because of the similarity in flux density levels between the AB86 and the present data. The detection of non-thermal emission in 1982 and in 2001 suggests a period of around 20 years. On the other hand, the last IR excess of the WC+O binary WR 125 started in 1991 (WH92), continued till 1993 and then faded away. IR monitoring performed during the last 18 years showed only one period of infrared excess (Williams 2002), implying an orbital period of at least 18 years. WR 125 has much in common with the archetype WR colliding wind binary WR 140 (HD 193793, WC7pd+O4-5, Williams et al. 1990), which has a thermal to non-thermal radio light curve with a period of 7.94 yr (Williams et al. 1994). If WR 125 is an analog of WR 140, then its 1992 IR excess occurred at periastron passage in an eccentric orbit, and an IR excess will happen again in $``$2010 - 2012.

WR 133, WN5+O9I (Fig. 18). We detect a double point-like source (separation $``$ 10″) with the southern component coincident with the optical position of the WN5+O9I binary (Fig. 18, upper panel), which is a known visual binary ($`\mathrm{\Delta }\varphi `$ = 5$`\stackrel{}{\mathrm{.}}`$4, Hartkopf et al. 1999; van der Hucht 2001, table 22). The source to the NE has a flux density of $`S_{3.6\mathrm{cm}}`$ = 0.39 $`\pm `$ 0.03 mJy beam-1. The lower panel of Fig. 18 displays an overlay of the radio continuum emission and the DSS2 red image, suggesting the presence of an additional star NE of the WR system, coincident in position with the northern radio source. Pollock et al. (1995) found a ROSAT X-ray excess for WR 133, indicative for a colliding wind binary.

Finally, we note that a radio source with a flux density $`S_{3.6cm}`$ = 0.15 $`\pm `$ 0.04 mJy is detected 2$`\stackrel{}{\mathrm{.}}`$5 far from the optical position of WR 4. We believe that this star has probably been barely detected, but, because of the low S / N ratio and the difference between the optical and radio positions (2.1$`\sigma `$), we have listed this star as non-detected.

In summary, we note that WR 79a, WR 89, WR 98a, WR 104, WR 105 and WR 125 are variable, which can be explained by the presence of a non-thermal component. This implies that our observations of these six stars cannot really be used for mass loss rate determinations, or will yield only upper limits. Of the above six cases all three WCd objects (WR 98a, WR 104 and WR 125) turn out to be variable non-thermal radio sources, suggesting a large number of non-thermal cases between these objects. Additional observations are necessary to investigate the nature of the other three WCd stars in our sample (WR 95, WR 103 and WR 113) since previous observations only provided upper limits to the flux densities.

## 6 Radio mass loss rates

Following Panagia & Felli (1975) and Wright & Barlow (1975), the flux density due to thermal Bremsstrahlung in an ionized expanding stellar wind can be expressed as

$$S_\nu =2.32\times 10^4(\frac{\dot{M}Z}{v_{\mathrm{}}\mu })^{4/3}(\frac{\gamma g_{\mathrm{ff}}\nu }{d^3})^{2/3}(1),$$

where $`S_\nu `$ is the flux density in mJy, $`\dot{M}`$ is the mass loss rates in M yr-1, $`v_{\mathrm{}}`$ is the terminal velocity in km s-1, $`\nu `$ is the frequency in Hz and $`d`$ is the distance in kpc. $`\mu `$, $`Z`$ and $`\gamma `$ are the mean molecular weight, the rms ionic charge and the mean number of electrons per ion, respectively. This expression corresponds to a spherically symmetric, stationary, isothermal wind flowing at constant velocity. The free-free Gaunt factor $`g_{\mathrm{ff}}`$ can be approximated, following Leitherer & Robert (1991) as

$$g_{\mathrm{ff}}=9.77(1+0.13log\frac{T_\mathrm{e}^{3/2}}{Z\nu })(2).$$

In this expression, $`T_\mathrm{e}`$ is the electron temperature of the wind in K.

Based on Equations (1) and (2), we have derived mass loss rates $`\dot{M}`$ corresponding to the observed flux densities at 3.6 cm, assuming that the observed radio emission is due to free-free emission. In the case of non-thermal contributions, our values represent upper limits to the true mass loss rate. All non-detections are indicated as upper limits to $`\dot{M}`$. Following LC97, we adopted $`T_\mathrm{e}`$ = 104 K in expression (2). As pointed out by these authors, this temperature has a very minor effect on the $`g_{\mathrm{ff}}`$. Values for the mean molecular weight, the rms ionic charge and the mean number of electrons per ion were taken from LC97. The values adopted for the WO star in our sample, WR 142, are those for the WC4 subtype (Crowther et al. 1998).

Table 4 lists the derived mass loss rates $`\dot{M}`$ and upper limits. Uncertainties in these values come from several sources. Errors in mass loss rate determinations were estimated following eq. (5) in LC97. We adopted similar criteria as LC97 to estimate the uncertainty in each parameter. For the error in the distances we assumed 20% in case of cluster/associations distances and 50% in case of field stars. The logarithmic uncertainty is then $`\pm `$0.09 and $`\pm `$0.25, respectively. Errors in $`g_{\mathrm{ff}}`$ are about 10% (0.05 in logarithmic scale). We adopted uncertainties in wind terminal velocities of 10% (i.e., a logarithmic error of $`\pm `$0.04). Finally, logarithmic errors in $`Z`$, $`\gamma `$ and $`\mu `$ were adopted as $`\pm `$0.08. Typical logarithmic errors in log $`\dot{M}`$ are in the range $`\pm `$ 0.20 to $`\pm `$ 0.24 in case of cluster/association distances and $`\pm `$ 0.40 in case of field stars. With the quoted uncertainties for the stellar distances and other wind parameters, uncertainties in mass-loss rates determinations are dominated by errors in distances, while uncertainties in flux densities have a minor impact.

————–

TABLE 4

————–

Other sources of error are non-thermal contributions to the flux at 3.6 cm. Inhomogeneities or clumping in the stellar winds lead to over estimates of the final values as emphasized by, e.g., Contreras et al. (1996) and Morris et al. (2000). The latter found volume filling factors of the range 0.04 - 0.25, leading to mass-loss rates a factor of 2.5 - 5 lower than expected for smooth, homogeneous winds. These effects have not been considered here.

## 7 Discussion

In order to compare mass loss rates derived for different WR subtypes from available radio surveys, we plotted the ATCA and VLA mass loss rates versus WR subtype in Fig. 21. The ATCA surveys have been performed by LC95, LC97 and CL99. Previous VLA results were published by Willis (1991) (who updated AB86 results by taking into account new criteria in terminal velocities and in the chemical composition of the region where radio emission originates) and by Monnier et al. (2002).

————–

FIGURE 21

————–

Monnier et al. observed WR 98a and WR 104 at several frequencies and found a non-thermal component in addition to the thermal one. They separated both components and derived a mass-loss rate from the best-fitting value to the thermal emission. The $`\dot{M}`$-values they derived are 0.8$`\times `$10-5 M yr-1 for WR 98a and 0.5$`\times `$10-5 M yr-1 for WR 104.

The upper panel of Fig. 21 displays the mass loss rate ($`\mathrm{log}\dot{M}`$) as a function of the spectral WR subtype for the stars included in our VLA survey. Bona-fide mass loss rate estimates are indicated by filled circles, while open circles indicate upper limits both for the variable detected sources and the undetected sources. Uncertainties in the derived values are also indicated. Note that we have plotted the mass loss rates derived by Monnier et al. instead of our estimates listed in Table 4.

The lower panel of Fig. 21 shows all radio mass loss rates from the ATCA and VLA samples. Although more than one mass loss rate determination is available for some of the stars, we plotted only one value for each star from (in order of preference) Monnier et al. (2002), new VLA data, ATCA data and previous VLA surveys. In the cases of non-thermal sources listed by LC97 and CL99, the $`\dot{M}`$-values plotted are upper limits. As for the upper panel, open circles indicate upper limits both for undetected sources and variable detected sources.

Excluding definite and suspected non-thermal cases (WR 79a, WR 89 and WR 105) the average mass loss rate for the WN stars detected in our survey is $`\dot{M}`$(WN) = \[4.3 $`\pm `$ 2.5\]$`\times `$10-5 M yr-1, equal to the average WN mass loss rate in ATCA and previous VLA data sets of $`\dot{M}`$(WN) = 4$`\times `$10-5 M yr-1.

Mass loss rates determined from emission lines studies by Nugis & Lamers (2000) differ from values derived from radio continuum data: in some cases (WR 104 and WR 113), $`\dot{M}`$ derived from emission lines are factors of two smaller than our radio estimates, while in other cases (WR 81, WR 95 and WR 103) factors of two larger than our results. Clumping and/or an ionization or temperature structure different from our assumptions may account for the differences.

We have thirteen WC8-9 stars included in our sample and detected seven of them. Four of these stars had no radio detections yet (WR 88, WR 95, WR 103 and WR 113), thus increasing to 12 the number of WC8-9 stars with radio detections. Most of these stars have heated circumstellar dust shells (see Williams 1995). WR 95, WR 103 and WR 113 are known to have persistent dust shells (see Table 1). The distribution of these stars in the diagram displays the same trend as in fig. 17 of LC97, with relatively low mass loss rates. The average mass loss rate for WC8-9 stars detected in our sample, excluding definite non-thermal cases (WR 98a and WR 104) is $`\dot{M}`$(WC8-9) = \[1.7 $`\pm `$ 1.0\]$`\times `$10-5 M yr-1. The average mass loss rate for WC5-7 stars, obtained from bona-fide values from the ATCA (5 stars) and the previous (2 stars) and present (1 star) VLA samples, is $`\dot{M}`$(WC5-7) = \[4.4 $`\pm `$ 1.6\]$`\times `$10-5 M yr-1, higher than the value derived for WC8-9 stars. Note, on the other hand, that average values for WN and WC5-7 stars are equal. Uncertainties in stellar distances clearly contribute to the observed scatter of mass loss rates.

The presence of dust in the wind envelopes might decrease the escaping UV photon flux, lowering the thermal radio emission, as suggested by Monnier et al. (2002), and thus the derived mass loss rates. LC97 suggest that this may explain the low mass loss rates obtained for WC9 stars from radio data. Note that the mass loss rates for the WC8 star WR 98a, which also displays a dust persistent shell, is also relatively low (see Fig. 21, lower panel).

The percentage of non-thermal emitters between the detected sources in our sample is 20 - 30 %. A similar value is derived for the 1986 VLA survey (which includes also previous results). Our percentage is lower than the one derived from ATCA data (40 %).

The case of WR 89 (suggested to be a binary by van der Hucht 2001), which shows time variations in flux density (see Table 3), hints to additional non-thermal emission (the radio spectrum was thermal at the time of the ATCA observations). WR 79a also appears to be variable in flux, suggesting the presence of a non-thermal component. However, a variable thermal component can not be discarded (e.g., Watson et al. 2002). It would be desirable to perform multifrequency radio observations of WR 79a and WR 89 to look for variability and to analyze the nature of the flux variations.

Finally, we would like to comment on the undetected stars. We repeat here that uncertainties in the expected flux densities $`S_{\mathrm{exp}}`$ are quite large, arising mainly from uncertainties in distances of WR field stars and in the expected range of mass loss rates (1 - 5$`\times `$10-5 M yr-1). $`S_{\mathrm{exp}}`$ values for WR 79b, WR 92, WR 96 and WR 119 are 0.12 - 0.26, 0.09 - 0.14, 0.10 - 0.16 and 0.11 - 0.16 mJy, respectively, which are below our detection limit (see Table 2). $`S_{\mathrm{exp}}`$ values were estimated adopting the stellar parameters of Tables 1 and 4, and mass loss rates in the range (1.5 - 2.0)$`\times `$10-5 M yr-1. For the cases of WR 80 and WR 121, however, $`S_{\mathrm{exp}}`$ = 0.5 - 0.8 and 0.4 - 0.6 mJy, respectively. The derived upper limit is 0.17 mJy (3$`\sigma `$) for both stars. Although a mass loss rate lower than typical can not be discarded, the objects may be at larger distances than adopted here, thus lowering the expected flux density.

## 8 Conclusions and suggestions for future work

We performed a survey of 34 Galactic Wolf-Rayet stars at 3.6 cm within the declination range of the Very Large Array. We report 15 definite and 5 probable detections. Of these 20 sources, 13 were detected for the first time at radio frequencies.

We confirm time variations in flux of WR 98a, WR 104, WR 105, and WR 125, which support the presence of a non-thermal component. For WR 125 (WC7ed+O9III) we suggest a binary period lower than 20 yr. WR 79a (HD 152408, WN9ha) and WR 89 are also variable in flux and we suspect they are also non-thermal emitters. Non-thermal radiation is indicative for wind-wind collision in WR+OB binaries. Thus, of our sample 20 - 30 % of the detected stars are non-thermal radio sources.

Averages of mass loss rate determinations yield similar values for WN (all subtypes) and WC5-7 stars ($`\dot{M}`$(WN) = \[4 $`\pm `$ 3\]$`\times `$10-5 M yr-1 and $`\dot{M}`$(WC5-7) = \[4 $`\pm `$ 2\]$`\times `$10-5 M yr-1), while a slightly lower value was derived for WC8-9 stars ($`\dot{M}`$(WC8-9) = \[2 $`\pm `$ 1\]$`\times `$10-5 M yr-1).

Future observations at several radio frequencies should monitor the objects WR 79a, WR 89, WR 98a, WR 104 and WR 125 during long time periods. These observations can help to determine the presence of non-thermal contributions to the flux density as a function of time. Radio observations at low and high frequencies of the sources first detected will be useful to investigate the nature of the radio emission.

Acknowledgements. We thank the anonymous referee for very helpful comments and suggestions that help to improve the presentation of the paper. C.C. is grateful to the kind hospitality during her stay at AOC, NRAO in Socorro, New Mexico. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research was partially supported by Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata; project 11/G049 (UNLP); and CONICET project PIP 607/98, Argentina. C.C. acknowledges a travel grant from Fundación Antorchas, Argentina, through project 13622/10.