text:
# 1 Introduction

## 1 Introduction

The determination of the critical temperature and the corresponding critical exponents for the systems exhibiting second-order phase transitions is the major task in many areas of physics. In order to determine the critical parameters of the system there exist two approaches. Traditionally, one uses universality and the scaling relations in the equilibrium stages of the statistical mechanical system. In other words, one studies the long-time regime of the system. The second approach is to study dynamic scaling which seems to exist in the early stages of the quenching process in the system exhibiting second-order phase transition. Jansen, Schaub and Schmittmann  have shown that for a dynamic relaxation process, in which a system is evolving according to a dynamics of Model A  and is quenched from a very high temperature to the critical temperature, a universal dynamic scaling behaviour within the short-time regime exists .

Monte Carlo simulations are the most commonly used tools to study statistical mechanical systems in thermal equilibrium around the phase transition point. Creating statistically independent configurations is controlled by the autocorrelations in the system which depend on the system as well as the simulation algorithm. Around the phase transition point, the correlation length grows exponentially and this results in the effect known as the critical slowing down. On finite lattices the autocorrelation time diverges as $`\tau L^z`$, where $`L`$ is the linear size of the lattice and $`z`$ is the dynamic critical exponent.

Dynamic finite-size scaling relation is a generalized form of the finite-size scaling relation valid for the equilibrium stages. For the $`k^{th}`$ moment of the magnetization of a system with spatial size $`L`$ this relation can be written as

$$M^{(k)}(t,ϵ,m_0,L)=L^{(k\beta /\nu )}^{(k)}(tL^z,ϵL^{1/\nu },m_0L^{x_0})$$
(1)
where $`\beta `$ and $`\nu `$ are the well-known critical exponents, $`t`$ is the simulation time and $`ϵ=(TT_c)/T_c`$ is the reduced temperature. In the study of the short-time dynamic scaling of a system, in addition to the critical exponents, the dynamic critical exponent $`z`$ and a new and independent exponent $`x_0`$ appear, where $`x_0`$ is the anomalous dimension of the initial magnetization $`m_0`$.

Eq.(1) can be reformulated for the time evolution of the magnetization at early stages of the simulation. For very small initial magnetization, time evolution of the magnetization undergoes a critical initial increase  in the form

$$M(t)m_0t^\theta $$
(2)
where the exponent $`\theta `$ is related to the exponent $`x_0`$ by $`\theta =(x_0\beta /\nu )/z`$. The power-law behaviour relation (Eq.(2)) of the dynamic finite-size scaling is used to estimate the critical exponents $`z`$, $`\theta `$ and $`\beta /\nu `$. In literature, it has been shown that the dynamic behaviour of various spin models obeys the dynamic scaling relations (Eqs. (1) and (2)) ; the universality of $`x_0`$ has also been studied and numerical evidence for the existence of the universality of the dynamic scaling has been given .

In this work we aim to reconsider dynamic finite-size scaling by approaching the problem from a different point of view. Instead of the critical initial increase (Eq. (2)) of the magnetization used in the literature, we have used the dynamic scaling relation (Eq.(1)) by considering the Taylor series expansion of $`M^{(k)}(t,ϵ,m_0,L)`$ with respect to small initial magnetization $`m_0`$. We have seen that the new exponent ($`x_0`$) exists not only in the initial stages, but also until the effects of the initial magnetization of the system diminish. This approach seems to give the critical exponents and the dynamic exponents with very high accuracy. In order to discuss the proposed method we have considered $`2`$\- and $`3`$-dimensional Ising models which are well-studied by finite-size scaling techniques for both equilibrium and non-equilibruium stages.

In the next section we present the model, the method and the scaling relation for the magnetization. In section 3 numerical results for $`2`$\- and $`3`$-dimensional Ising models will be reported. The last section contains the conclusions.

## 2 Model and the Method

In this work we have employed $`2`$\- and $`3`$-dimensional Ising models which are described by the Hamiltonian

$$\beta H=K\underset{<ij>}{}S_iS_j.$$
(3)
Here $`\beta =1/kT`$ and $`K=J/kT`$, where $`k`$ is the Boltzmann constant, $`T`$ is the temperature and $`J`$ is the magnetic interaction between the spins. In the Ising Model the spin variables take the values $`S_i=\pm 1`$. In this work we have performed Monte Carlo simulations by using Metropolis algorithm.

The measurements of the thermodynamic quantities which are obtained using thermalized configurations in Monte Carlo simulations are independent of the initial configuration. Here, before the thermalization, the convergence path (to the thermalization) may be different depending on which initial configuration is used. Although the system reaches the same thermalized configurations in the long-time regime, the convergence path is highly sensitive to the initial configuration. In this work we propose to consider the dynamics of the thermalization for a system with two different initial configurations, in order to obtain the new dynamic index $`x_0`$. One initial configuration can be considered as the configuration with vanishing initial magnetization. The second initial configuration can be taken with arbitrary, but small initial magnetization. The dynamic finite-size scaling form of the $`k^{th}`$ moment of the magnetization is given by Eq.(1). At the critical point $`(ϵ=0)`$, for small initial magnetization ($`m_0`$) values, the function given in Eq.(1) can be expanded into Taylor series around $`m_0=0`$ resulting in

$`M^{(k)}(t,0,m_0,L)`$ $`=`$ $`L^{k\beta /\nu }^{(k)}(tL^z,0,0)+`$
$`{\displaystyle \frac{^{(k)}(tL^z,0,m_0)}{m_0}}_{m_0=0}m_0L^{x_0k\beta /\nu }+O(m_{0}^{}{}_{}{}^{2}).`$
The first term in this expansion has no $`m_0`$ dependence and it can be calculated through Monte Carlo simulations by considering the initial configuration with vanishing initial magnetization. The $`m_0`$-dependent terms (second and higher order terms in the expansion) can be represented by a function $`M_s(t,0,m_0,L)`$ which is basically

$$M_s(t,0,m_0,L)=\frac{^{(k)}(tL^z,0,m_0)}{m_0}_{m_0=0}m_0L^{x_0k\beta /\nu }+O(m_{0}^{}{}_{}{}^{2})$$
(5)
and it can be obtained by calculating the difference of $`M^{(k)}(t,0,m_0,L)`$ and $`M^{(k)}(t,0,0,L)`$. This means it has the form,

$$M_s(t,0,m_0,L)=M^{(k)}(t,0,m_0,L)M^{(k)}(t,0,0,L).$$
(6)
Here $`M^{(k)}(t,0,0,L)`$ corresponds to the first term in the expansion. The terms $`M^{(k)}(t,0,m_0,L)`$ and $`M^{(k)}(t,0,0,L)`$ can be measured through Monte Carlo simulations using two different initial configurations, one with initial magnetization $`m_0`$ and the other one with vanishing initial magnetization.

The function given in Eqs.(5) and (6) contains all $`m_0`$ dependence. The thermalized configurations are independent of the initial magnetization, hence $`M_s(t,0,m_0,L)`$ is a decaying function. This function is expected to be a function of the autocorrelation time and the initial magnetization. In this work, we have chosen an exponentially decaying function in the form

$$M_s(t,0,m_0,L)=f(m_0,L)e^{(t/\tau (L))^{\varphi (L)}}.$$
(7)
Here $`t`$ is the simulation time, $`\tau (L)`$ is the autocorrelation time and $`\varphi (L)`$ is an exponent which is related to the initial slip. If this form fits the data, two dynamical exponents can be calculated separately. $`f(m_0,L)`$ is a $`m_0`$-dependent function, but it is has no autocorrelation time dependence. Hence the new dynamic scaling exponent $`x_0`$ can be calculated from the size dependence of this function for small values of $`m_0`$. The autocorrelation time and hence $`z`$ can be calculated from the decay rate of $`M_s(t,0,m_0,L)`$. The validity of our assumptions has been tested by fitting the function given in Eq(7) to our data and the results are presented in Section 3.

For very small initial magnetization values, it is expected that the leading term in the expansion (Eq. (2)) is dominant and the rest of the terms (represented by the function $`M_s(t,0,m_0,L)`$) is in the same order of magnitude with the simulation errors. Increasing initial magnetization results in a sizeable $`m_0`$-dependent term in Eq. (2). In this range of $`m_0`$ one can expect that $`M_s(t,0,m_0,L)`$ scales like $`L^{x_0\beta /\nu }`$. Further increase in $`m_0`$ reduces the reliability of such scaling since the contribution of higher order terms in the expansion (Eq. (2)) becomes effective. In our simulations we have used different $`m_0`$ values to investigate the proposed scaling behaviour of $`M_s(t,0,m_0,L)`$. To test the scaling we have used the literature values of $`\beta /\nu `$ and $`x_0`$ for $`2`$\- and $`3`$-dimensional Ising models and also calculated the dynamic exponents $`z`$ and $`x_0`$ for $`3`$-dimensional model through a fit to our Monte Carlo data.

## 3 Results and Discussions

We have studied $`2`$-and $`3`$-dimensional Ising models evolving in time according to the dynamics of the Model A. We have prepared the lattices with varying initial magnetization values starting from $`m_0=0`$. Totally random initial configurations are quenched at the corresponding infinite lattice critical temperature. Simulations are performed by using Metropolis algorithm on lattices with linear sizes $`L=512,640,768,896`$ and $`L=32,64,96`$ for $`2`$\- and $`3`$-dimensional Ising models respectively. For error calculations we have binned data points obtained from $`4000`$ to $`15000`$ independent runs for $`3`$-dimensional model, and $`50000`$ to $`500000`$ independent runs for $`2`$-dimensional model, starting from randomly prepared configurations. For $`3`$-dimensional model, $`3000`$ iterations are performed for each independent run. For $`2`$-dimensional model there are two different sets of runs. One set of data is obtained after $`500`$ iterations for each independent run and for the other set of data, each independent run includes $`10000`$ iterations.

The dynamic critical exponent $`z`$ has a little dependence on dimension for local algorithms. For a given numbers of spins, the saturation in magnetization can be observed after relatively small number of iterations in $`3`$-dimensional lattice compared two $`2`$-dimensional one, since the saturation is observed when the correlation length reaches the linear size of the lattice. In order to observe the behaviour of the $`m_0`$-dependent terms, one needs to study the time development of the system until both $`M^{(k)}(t,0,m_0,L)`$ and $`M^{(k)}(t,0,0,L)`$ reach approximately the same value. In our simulation studies, we have observed that for $`3`$-dimensional lattice, configurations with both vanishing and non-vanishing initial magnetization reach the same value in (approximately) $`3000`$ iterations for the largest lattice $`(L=96)`$. For $`2`$-dimensional model, since much larger linear sizes are used, the initial configuration dependence remains after $`30000`$ iterations and observing the decaying behaviour of $`M_s(t,0,m_0,L)`$ becomes extremely expensive. For this reason, we have tested our ideas for $`3`$-dimensional model. In $`2`$-dimensional lattice, we have performed simulations up to $`10000`$ iterations, but even this many iterations is extremely short to observe the decaying behaviour of $`M_s(t,0,m_0,L)`$. Nevertheless, it is sufficient to test the scaling behaviour and to obtain the new dynamic critical exponent $`x_0`$.

In Figure 1, the magnetization data are presented for vanishing initial magnetization and for $`m_00.1`$ in $`3`$-dimensional Ising Model with linear sizes $`L=32,64,96`$. The magnetization data converge approximately to the same value for both initial configurations after approximately $`500(L=32),1500(L=64),3000(L=96)`$ iterations. The differences between the data obtained by starting from different initial magnetizations are large initially for large lattices and vanishes after above mentioned number of iterations are reached. The final magnetization value decreases as the lattice size increases, and the magnetization is smaller initially for vanishing $`m_0`$. The function $`M_s(t,0,m_0,L)`$ which represents the difference between $`M^{(k)}(t,0,m_0,L)`$ and $`M^{(k)}(t,0,0,L)`$ is given in Figure 2, for $`m_00.1`$ and for linear sizes, $`L=32,64,96`$ in $`3`$-dimensional lattice. The function $`M_s(t,0,m_0,L)`$ is a decreasing function with an initial value proportional to the lattice size $`L`$. Despite the fact that our main aim is not to obtain the dynamical critical exponent $`(z)`$, since we are interested in the functional form of the $`m_0`$-dependent term and the speed of its decay, we have fitted our data to the functional form $`(\tau L^z)`$ and obtained $`z`$ from the fit results. We have fitted the data to the exponential function given in Eq.(7). This function fits the data reasonably well, as it can be seen from Figure 2. The correlation times $`(\tau )`$ obtained from the fit are $`\tau (L)=80,\mathrm{\hspace{0.33em}372},\mathrm{\hspace{0.33em}839}`$ for $`L=32,64,96`$ respectively. The errors are in the order of one. We have fitted the functional form $`(\tau L^z)`$ to our data and we have obtained the value for $`z`$ as $`z=2.04\pm 0.07`$ for $`3`$-dimensional model. Figure 3 shows the fitted function and the data points for three lattice sizes used.

As it is mentioned in Section 2, the function $`M_s(t,0,m_0,L)`$ representing $`m_0`$-dependent terms scales like $`L^{z\theta }`$. Here we have used the literature values of the exponent $`\theta =(x_0\beta /\nu )/z`$, $`\theta =0.191(2)`$  and $`\theta =0.108(2)`$  for $`2`$\- and $`3`$-dimensional Ising models respectively, to test the scaling. The exponent $`z\theta =x_0\beta /\nu `$ is obtained from these literature values of both $`\theta `$ and $`z`$. For $`2`$-dimensional Ising Model, the value of the dynamic critical exponent $`z`$ for Metropolis algorithm varies from $`2.13`$ to $`2.18`$ . For the $`3`$-dimensional case, the equilibrium case, the value of the dynamic critical exponent is $`z=2.04(3)`$ , while from the dynamic scaling, $`z`$ ranges from $`2.032(11)`$ to $`2.052(17)`$ . We have used $`z=2.18`$ and $`z=2.04`$ for $`2`$\- and $`3`$-dimensional Ising models respectively. The scaling of the function $`M_s(t,0,m_0,L)`$ obtained from the Monte Carlo data for $`3`$-dimensional model is given in Figure 4. As it can be seen from this figure, these terms are scaled reasonably well. We have also fitted the data to a function to obtain the best fit and to obtain the value for $`z\theta `$. Our value $`z\theta =0.27\pm 0.03`$ (for $`3`$-dimensional model) is slightly larger than the value in the literature.

For $`2`$-dimensional Ising Model, the autocorrelation time $`(\tau )`$ is extremely long for the lattice sizes we use. Consequently, we didn’t aim to see the decaying form of $`M_s(t,0,m_0,L)`$. Instead we have tested the scaling behaviour of $`M_s(t,0,m_0,L)`$. In Figure 5, we have presented the data for $`M_s(t,0,m_0,L)`$ for after $`500`$ iterations and Figure 6 shows the scaling of the data when the literature values are used.

## 4 Conclusions

The short-time dynamic behaviour of the magnetization is studied in $`2`$\- and $`3`$-dimensional Ising models. The time-dependent magnetization profiles with different initial magnetization ($`m_0`$) values exhibit an initial increase during the time scale set by the initial magnetization, the dynamic critical exponent ($`z`$) and the new dynamic exponent ($`x_0`$). After a long-time ($`tL^z`$), all profiles gradually join together, since in the long-time regime the effects of the initial conditions are ignorable. In this work we aimed to study the effects of the initial configurations on magnetization by using Monte Carlo simulation. This can be observed by expanding the magnetization around $`m_0=0`$. This expansion (Eq. (1)) gives a leading term which is independent of $`m_0`$ and the terms containing the information of how the initial configuration is prepared.

In this approach, the dynamic critical behaviour can be studied by using the data from the early stages of the Monte Carlo simulation. In order to obtain a relation between the size and the initial magnetization dependence of the Monte Carlo data, we have considered both vanishing and non-vanishing initial magnetization runs. The iteration by iteration differences of these two sets of runs yield the initial magnetization dependence of the system. Using different lattice sizes we have shown that the proposed scaling holds. We have also shown that the dynamic finite-size scaling is valid not only for the very early stages (within $`100`$ iterations) of the Monte Carlo runs, but it is valid until the effects of the the initial magnetization diminish.

In order to test our proposed approach, we have employed the literature values of $`\beta /\nu `$, $`z`$, $`x_0`$ and we have observed a good scaling behaviour for $`m_0`$-dependent terms for the lattice sizes we used. Such a good scaling gives hopes for the applicability of this approach to a wide range of statistical mechanical systems for calculations of both dynamic and equilibrium critical exponents.

## Acknowledgements

We greatfully acknowledge Hacettepe University Research Fund (Project no : 01 01 602 019) and Hewlett-Packard’s Philanthropy Programme.

## Figure Captions

Figure 1. Magnetization data for $`3`$-dimensional Ising Model with linear lattice sizes $`L=32,64,96`$, with vanishing and nonzero initial magnetization ($`m_00.1`$) for each size, as a function of simulation time $`t`$. The curves for vanishing initial magnetization start from the origin.

Figure 2. The function $`M_s(t,0,m_0,L)`$ corresponding to $`m_0`$-dependent terms in the expansion as a function of simulation time $`t`$, for $`3`$-dimensional Ising Model. The plot of the function fitted to the form given in Eq.(7) is also shown.

Figure 3. The autocorrelation time $`\tau (L)`$ as a function of the linear lattice size $`L=32,64,96`$ in $`3`$-dimensional Ising Model. The function fitted to the form $`\tau L^z`$ is also shown.

Figure 4. Scaled form of $`M_s(t,0,m_0,L)`$ obtained from the simulation data (given in Figure 2), for $`3`$-dimensional Ising Model.

Figure 5. $`M_s(t,0,m_0,L)`$ obtained from the simulation data as a function of simulation time $`t`$ for $`2`$-dimensional Ising Model ($`m_00.07`$).

Figure 6. Scaled form of the data given in Figure 5 for function $`M_s(t,0,m_0,L)`$ for $`2`$-dimensional Ising Model.


format: text

text:
# Massive or Massless Scalar Field and Confinement

## 1 The model

Recently, Dick and Fulcher  have proposed an interesting model where the lowest glueball, represented by a massive scalar field $`\varphi `$, has been in a non-minimal way coupled to the $`SU(2)`$ gauge fields. They have chosen, in analogy with the chiral quark models , one–glueball coupling: $`\varphi F^{a\mu \nu }F_{\mu \nu }^a`$. Their model provides confinement of quarks and gives a quite reasonable interquark potential $`U_{q\overline{q}}r^{1/3}`$.
However, as it was pointed out by Zalewski and Motyka, the most probable potential, which gives the best fit to the quarkonium system takes a different form , namely

$$U_{MZ}(r)=C_1\left(\sqrt{r}\frac{C_2}{r}\right),$$
(1)
where $`C_10.71`$ Gev$`^{\frac{1}{2}}`$ and $`C_20.46`$ Gev$`^{\frac{3}{2}}`$. Due to that the model should be modified. It can be done by using slightly more complicated, effective coupling. It has to be underlined that in spite of that modification the scalar field still represents the lowest scalar glueball.
In the present paper we focus on the following scalar–gauge action

$$S=d^4x\left[\frac{1}{4}\frac{(\frac{\varphi }{\mathrm{\Lambda }})^{8\delta }}{1+(\frac{\varphi }{\mathrm{\Lambda }})^{8\delta }}F_{\mu \nu }^aF^{a\mu \nu }+\frac{1}{2}_\mu \varphi ^\mu \varphi \frac{m^2}{2}\varphi ^2\right].$$
(2)
Here, $`m`$ is a mass of the scalar field, $`\mathrm{\Lambda }`$ is a dimensional and $`\delta `$ a dimensionless constant. In order to have the standard Maxwell–like behavior of the fields in the neighborhood of sources the particular form of the denominator has been added. It does not inflect the long range behavior of the fields but provides the asymptotic freedom in the limit of the strong scalar field .
The field equations corresponding to the action (2) take the form

$$D_\mu \left(\frac{(\frac{\varphi }{\mathrm{\Lambda }})^{8\delta }}{1+(\frac{\varphi }{\mathrm{\Lambda }})^{8\delta }}F^{a\mu \nu }\right)=j^{a\nu },$$
(3)
and

$$_\mu ^\mu \varphi =2\delta F_{\mu \nu }^aF^{a\mu \nu }\frac{\varphi ^{8\delta 1}}{\mathrm{\Lambda }^{8\delta }\left(1+\frac{\varphi ^{8\delta }}{\mathrm{\Lambda }^{8\delta }}\right)^2}m^2\varphi ,$$
(4)
where $`j^{a\mu }`$ is an external current.
In the next section it will be shown that the action (2) possesses confining solutions in the massive sector. Moreover, these solutions describe a wide family of confining potentials which can be compared with various potentials found in fits to the phenomenological data , , .
Interestingly enough, there is also a confining solution in the massless sector, where the scalar field cannot be identified with any glueball. Finally we prove that model with the non-dynamical scalar field provides confinement of external sources as well. Problem of glueball states in these models is discussed in the last section.

## 2 Confining solutions

### 2.1 Massive case

Let us now discuss the Coulomb problem. We find solutions of the field equations generated by an external, static, point-like source:

$$j^{a\mu }=4\pi q\delta (r)\delta ^{a3}\delta ^{\mu 0}.$$
(5)
The restriction to the Abelian source is not essential. It is possible to investigate more general, static non-Abelian source : $`j_\mu ^a=4\pi q\delta (r)C^a\delta ^{\mu 0}`$. Here $`C^a`$ is the expectation value of the $`su(N_c)`$ generator for a normalized spinor in color space. However, fields generated by above sources have identical dependence on the spatial coordinates. They differ only by a multiplicative color factor. Due to that it is sufficient to analyze only Abelian source (5). Then the field equations read

$$\left[r^2\frac{\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta }}{1+\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta }}E\right]^{}=4\pi q\delta (r),$$
(6)
$$_r^2\varphi =4\delta E^2\frac{\varphi ^{8\delta 1}}{\mathrm{\Lambda }^{8\delta }\left(1+\frac{\varphi ^{8\delta }}{\mathrm{\Lambda }^{8\delta }}\right)^2}+m^2\varphi .$$
(7)
We use the standard definition $`E^{ai}=F^{a0i}`$. The electric field is chosen in the same color direction as the source $`\stackrel{}{E}^a=E(r)\delta ^{3a}\widehat{r}`$. One can easily derive the electric field from equation (6) and express it in terms of the scalar field:

$$E(r)=\frac{q}{r^2}\left(1+\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta }\right).$$
(8)
This field can be understood as usual Coulomb field $`E(r)=\frac{\alpha _{eff}}{r^2}`$ in unusual medium. The medium is characterized by the scalar field which modifies the effective coupling constant $`\alpha _{eff}`$:

$$\alpha _{eff}=q\left(1+\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta }\right).$$
(9)
Using relation (8) we are able to rewrite the equation (7) in the following form:

$$_r^2\varphi =4\delta \frac{q^2}{\mathrm{\Lambda }r^4}\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta 1}+m^2\varphi .$$
(10)
Unfortunately, this equation is still too complicated to integrate it analytically. However, on account of the fact that we are mainly interested in the long range behavior of the fields we can analyze (10) in the asymptotic regime i.e. for $`r\mathrm{}`$. The asymptotic solution is found to be:

$$\varphi \left(\frac{2\sqrt{\delta }q\mathrm{\Lambda }^{4\delta }}{m}\right)^{\frac{1}{1+4\delta }}\left(\frac{1}{r}\right)^{\frac{2}{1+4\delta }}.$$
(11)
Then the electric field:

$$E\frac{q}{r^2}+\left(\frac{4\delta q}{m^2}\right)^{\frac{4\delta }{1+4\delta }}\mathrm{\Lambda }^{\frac{8\delta }{1+4\delta }}\left(\frac{1}{r}\right)^{\frac{28\delta }{1+4\delta }}.$$
(12)
Finally, the corresponding potential has the form:

$$U\frac{q}{r}+\frac{1+4\delta }{12\delta 1}\left(\frac{4\delta q}{m^2}\right)^{\frac{4\delta }{1+4\delta }}\mathrm{\Lambda }^{\frac{8\delta }{1+4\delta }}\left(\frac{1}{r}\right)^{\frac{112\delta }{1+4\delta }},$$
(13)
for $`\delta \frac{1}{4}`$ or $`\frac{1}{12}`$. One can observe that $`\delta =\frac{1}{12}`$ corresponds to logarithmic behavior of the electric potential:

$$U\frac{q}{r}+\left(\frac{q}{3m^2}\right)^{\frac{1}{4}}\mathrm{\Lambda }^{\frac{1}{2}}\mathrm{ln}r.$$
(14)
In the framework of such classical models the confinement is understood as appearance of singularity of the electric potential (energy density) at the spatial infinity. Energy becomes infinite due to the long range behavior of fields. This effect takes place for $`\delta \frac{1}{12}`$. However, it was shown by Seiler  that confining potentials cannot grow faster than linearly for large $`r`$. This gives us the upper bound for parameter $`\delta `$. Concluding, the model given by the action (2) can be used to modelling confinement of the external sources for the following parameter $`\delta `$:

$$\delta [\frac{1}{12},\frac{1}{4}]$$
(15)
The standard linear potential is obtained for $`\delta =\frac{1}{4}`$. Then we achieve a model discussed previously by Dick . In this model a global, confining solution was found. As one could expect the asymptotic behavior of the solution is linear.

### 2.2 Massless case

In this section we are going to focus on the problem of the confinement in case of the massless scalar field. Apparently, such a field can be no longer identified with the scalar glueball. Nonetheless, one can still use it to modify properties of the vacuum. The scalar field becomes an effective field describing dynamical permittivity of the medium i.e. vacuum. One can notice that the model belongs to the so–called color dielectric field theories.
Equation (10) for the scalar field takes simpler form now:

$$_r^2\varphi =4\delta \frac{q^2}{\mathrm{\Lambda }r^4}\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta 1},$$
(16)
which can be integrated analytically. We have found a family of solutions regular at the spatial infinity and labelled by a positive parameter $`\beta _0`$:

$$\varphi =A\mathrm{\Lambda }\left(\frac{1}{r\mathrm{\Lambda }}+\frac{1}{\beta _0}\right)^{\frac{1}{1+4\delta }},$$
(17)
where the constant $`A`$ is given by

$$A=[q(1+4\delta )]^{\frac{1}{1+4\delta }}.$$

The corresponding electric field has the form:

$$E=\frac{q}{r^2}+A^{8\delta }\frac{q}{r^2}\left(\frac{1}{r\mathrm{\Lambda }}+\frac{1}{\beta _0}\right)^{\frac{8\delta }{1+4\delta }}.$$
(18)
We see that there exists infinite set of electric solutions, generated by external point source, which falls as $`r^2`$ for large $`r`$. The behavior of the electric field near a source strongly depends on the value of the parameter $`\delta `$. If $`\delta 0`$ the singularity at $`r=0`$ is identical to that generated by point source in the standard Maxwell electrodynamics. If $`\delta <0`$ then the singularity can be arbitrarily large.
Besides the family of solutions mentioned above, there is a single field configuration, which solves the equations of motion, and gives confinement. The solution is:

$$\varphi =A\mathrm{\Lambda }\left(\frac{1}{r\mathrm{\Lambda }}\right)^{\frac{1}{1+4\delta }}.$$
(19)
The electric field is then given by the formula:

$$E=\frac{q}{r^2}+A^{8\delta }q\mathrm{\Lambda }^2\left(\frac{1}{r\mathrm{\Lambda }}\right)^{\frac{2}{1+4\delta }}.$$
(20)
Thus the electric potential for $`\delta \frac{1}{4}`$ takes the form:

$$U=\frac{q}{r}+\frac{4\delta +1}{4\delta 1}A^{8\delta }q\mathrm{\Lambda }^{\frac{8\delta }{1+4\delta }}\left(\frac{1}{r}\right)^{\frac{14\delta }{1+4\delta }}$$
(21)
and for $`\delta =\frac{1}{4}`$

$$U=\frac{q}{r}+qA^{8\delta }\mathrm{\Lambda }\mathrm{ln}(\mathrm{\Lambda }r).$$
(22)
In fact, we observe confining behavior (in the same sense as for the massive field) for the parameter $`\delta `$:

$$\delta [\frac{1}{4},\mathrm{}).$$
(23)
As it was recently shown , it is possible to get rid of the non-confining solutions (17), (18) adding the following potential term

$$V(\varphi )=\alpha \varphi ^4\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta },$$
(24)
where $`\alpha `$ is a positive constant. Then only the confining solution survives. Quite interesting, if we omit the denominator in the dielectric function in the action (2) (that means we are no longer interested in restoration of the standard $`r^2`$ behavior of the electric field in small neighborhood of a source) then the non-confining solutions describe the so–called screening phenomena known from quantum Yang–Mills theory . The field configuration generated by a fixed charge can have arbitrary small but positive energy .

### 2.3 Non-dynamical permittivity

Let us now proceed and discuss the third possibility that is to neglect the kinetic term of the scalar field in (2). In this case $`\varphi `$ is no longer a dynamical field . In the other words it is possible to treat it as an additional field which, after expressing it in terms of the gauge fields, can be completely removed from the action. One can easily observe that in order to deal with a non-trivial theory we are forced to add a potential term for the scalar field. In our calculation the potential is chosen as in the massless case (24).
Then the action takes the following form:

$$S=d^4x\left[\frac{1}{4}\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta }F_{\mu \nu }^aF^{a\mu \nu }\alpha \varphi ^4\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta }\right].$$
(25)
After variation with respect to the scalar field one gets

$$8\delta F\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta 1}4(2\delta +1)\alpha \varphi ^4\left(\frac{\varphi }{\mathrm{\Lambda }}\right)^{8\delta 1}=0,$$
(26)
where $`F=\frac{1}{4}F^{a\mu \nu }F_{\mu \nu }^a`$. As previously we will investigate only the electric part of the theory. Then one can trivially solve it and find

$$\varphi ^4=\frac{E^2}{a},$$
(27)
where for simplicity a new constant $`a=(2\delta +1)\alpha `$ was defined. Substituting it into the Gauss law we obtain

$$\left[r^2\left(\frac{E^2}{a\mathrm{\Lambda }^4}\right)^{2\delta }E\right]^{}=4\pi q\delta (r).$$
(28)
The solution reads

$$E=a^{\frac{2\delta }{1+4\delta }}\left(\frac{|q|}{\mathrm{\Lambda }^2r^2}\right)^{\frac{1}{1+4\delta }}\mathrm{\Lambda }^2.$$
(29)
This gives us the electric potential

$$U=a^{\frac{2\delta }{1+4\delta }}|q|^{\frac{1}{1+4\delta }}\frac{4\delta +1}{4\delta 1}\left(\frac{1}{\mathrm{\Lambda }r}\right)^{\frac{14\delta }{1+4\delta }}\mathrm{\Lambda }$$
(30)
for $`\delta \frac{1}{4}`$ and

$$U=\sqrt[4]{aq^2}\mathrm{\Lambda }\mathrm{ln}(\mathrm{\Lambda }r)$$
(31)
for $`\delta =\frac{1}{4}`$. Identically as for the massless scalar field the confinement of the point sources is found for all parameter $`\delta `$ equal or larger than $`\frac{1}{4}`$. From the confinement point of view these models are equivalent.
As it was said before the action (25) can be expressed only by gauge fields. After a simple calculation one obtains the Pagels–Tomboulis lagrangian :

$$L=\frac{1}{4}F^2\left[\left(\frac{F^2}{\mathrm{\Lambda }^4}\right)^2\right]^\delta .$$
(32)
Of course, the pertinent equations of motion, in the electric sector, derived for the Pagels–Tomboulis model are the the same as (28) (up to unimportant multiplicative constant).
Models, where the permittivity depends on the strength tensor $`F`$ has been very successfully analyzed in the context of the effective theory for the low energy QCD , . Here no additional, effective scalar field is needed. Confining solutions emerge due to the modification of the vacuum by the gauge fields. In fact, in the Pagels–Tomboulis model a source charge generates infinite energy field configuration whereas a configuration originating from dipol-like source has finite energy , .

## 3 Conclusions

In the present paper we have analyzed the $`SU(2)`$ gauge field coupled to the scalar field. For the massive and massless scalar field as well as for the non-dynamical scalar field solutions describing confinement have been found. It is relatively easy to model confinement of external sources using the scalar dilaton/glueball field. Our models seem to be very similar from the confinement point of view. Of course, there are also some differences between them.
Firstly, the occurrence of confinement strongly depends on the value of the parameter $`\delta `$. The massive and massless models require different values of $`\delta `$. For example the best phenomenological quark–antiquark potential, found by Zalewski and Motyka in the fits to the spectrum of quarkonia, is obtained for $`\delta =\frac{3}{20}`$ for the massive and for $`\delta =\frac{3}{4}`$ for the massless case respectively. In fact, then $`U\sqrt{R}`$, where $`R`$ is the distance between sources. Another known interquark potential $`UR^{0.1}`$, used by Martin , is realized for $`\delta =\frac{11}{116}`$ (massive field) or for $`\delta =\frac{11}{36}`$ (massless field). The difference between these models is even more apparently seen in case of the linear confinement. The massive model describe it for $`\delta =\frac{1}{4}`$. It is unlikely for the massless model, where the linear potential is realized only in the limit $`\delta \mathrm{}`$, which cannot be implemented on the action level. The linear confinement can be only approximated with arbitrary accuracy by taking sufficiently large $`\delta `$.
To conclude, in spite of the differences mentioned above all the models considered here can serve very well to model confinement of quarks. Nevertheless these models differ essentially. One of the most profound differences is visible if we look at the way how glueball states appear in the models.
If we would like to treat these models as candidates for description of low energy QCD, problem of glueballs must be taken into consideration. It is obvious that a good effective model is expected not only to provide confinement of quarks (and give a reasonable quark–antiquark potential) but is also expected to have solutions which could be interpreted as glueballs. In the massive case the situation is clear. The scalar field represents the scalar glueball $`0^{++}`$ which has mass $`m`$. The dielectric function determines coupling between the glueball and gluons. For instance, $`\delta =\frac{1}{8}`$ gives one glueball coupling. For different values of $`\delta `$ we have more complicated, effective coupling. Unfortunately, this model is unable to describe other glueball states. In order to do it one has to enlarge the number of fields. For example the pseudo-scalar glueball $`0^+`$, represented by a new field $`\psi `$, requires axion-like coupling .
Situation changes drastically if we consider remaining models (with the dynamical or non-dynamical scalar field). Obviously, the interpretation that the scalar field corresponds to the glueball is no longer correct. Because of that, glueballs should appear in a different, non-trivial way. It is believed that in the framework of these models glueballs could be described as (topological) solitons i.e. static, localized and finite energy solutions of the sourcesless field equations. This is a very attractive way of studying glueballs. Different glueballs would be given by different solitons characterized by a topological index (cf. e.g. the Faddeev–Niemi model ).
As a result, existence of topologically non-trivial solutions in the framework of scalar–gauge models is crucial from the QCD point of view. It would let us justify which model can be really relevant as the effective model for the low energy quantum chromodynamics. We plan to address this problem in our proceeding paper.


text:
# Edge dislocations in crystal structures considered as traveling waves of discrete models

## Abstract

The static stress needed to depin a 2D edge dislocation, the lower dynamic stress needed to keep it moving, its velocity and displacement vector profile are calculated from first principles. We use a simplified discrete model whose far field distortion tensor decays algebraically with distance as in the usual elasticity. Dislocation depinning in the strongly overdamped case (including the effect of fluctuations) is analytically described. $`N`$ parallel edge dislocations whose average inter-dislocation distance divided by the Burgers vector of a single dislocation is $`L1`$ can depin a given one if $`N=O(L)`$. Then a limiting dislocation density can be defined and calculated in simple cases.

In many fields, genuinely microscopic phenomena affect macroscopic behavior in a way that is difficult to quantify precisely. Typical cases are the motions of dislocations nab67 ; hir82 , cracks fre90 , vortices in Josephson arrays zan95  or other defects subject to pinning due to the underlying crystal microstructure. Emerging behavior due to motion and interaction of defects might explain common but poorly understood phenomena such as friction ger01 . Macroscopic theories consider the continuum mechanics of these solids subject to forces due to the defects and additional equations for the densities of defects and properties of their motion ll7 . The latter are usually postulated by phenomenological considerations. An important problem is to derive a consistent macroscopic description taking into account the microstructure.

Here we tackle a simplified problem containing all the ingredients of the previous description: the pinning and motion of edge dislocations. Firstly, we study a two-dimensional (2D) discrete model kkl  describing the damped displacement of atoms subject to the field generated by a 2D edge dislocation and a constant applied shear stress of strength $`F`$. If $`|F|<F_{cs}`$ ($`F_{cs}`$ is related to the static Peierls stress), the stable displacement field is stationary, whereas the dislocation core and its surrounding displacement field move if $`|F|>F_{cs}`$. A crucial observation is that there exists a stable uniformly moving dislocation with both core and far field advancing at the same constant velocity. This suggests that a moving edge dislocation is a traveling wave of the discrete model.

Our self-consistent calculation self  based on this picture predicts the following magnitudes: (i) the critical static stress needed to depin a static dislocation, (ii) the dynamic stress $`F_{cd}<F_{cs}`$ below which a moving dislocation stops (in the strongly overdamped case, $`F_{cd}=F_{cs}`$), and (iii) the dislocation velocity as a function of applied stress. The latter information has to be taken from experiments in the standard treatment hul01 . Then a macroscopic quantity, the dislocation velocity, is obtained from analysis of a microscopic model.

Secondly, we consider a distribution of many parallel edge dislocations separated by macroscopic distances comprising many lattice periods. A dislocation cannot move under the influence of other dislocations far away unless the latter have finite density (there are $`N`$ such dislocations and the average distance between them is $`L`$ with $`N=O(L)`$ as $`L\mathrm{}`$). Under the influence of such a distribution, one dislocation may be pinned or move depending on the dislocation density. The latter is calculated at the critical stress in a simple configuration.

A simplified discrete model of edge dislocations. Consider an infinite three dimensional cubic lattice with symmetry axes $`x,y,z`$. We insert an extra half plane of atoms parallel to the plane $`yz`$. The border of this extra half plane is a line (parallel to the $`z`$ axis), which is called an edge dislocation in the crystal. The Burgers vector of the dislocation points in the $`x`$ direction and the plane $`xz`$ is the glide plane of the dislocation (see ll7 ). If we apply an external stress $`\sigma `$, the dislocation moves on its glide plane and in the direction of its Burgers vector in response to just one component of $`\sigma `$: the stress $`\tau `$ resolved on the glide plane in the glide direction hir82 ; ll7 . All the sections of the lattice by planes parallel to $`xy`$ look alike. Thus, we reduce the problem to a 2D lattice with an extra half line of atoms, see Fig. 22 of Ref. ll7, . Assuming that glide is only possible in the $`x`$-direction, the dynamics of an edge dislocation in a 2D lattice can be described by kkl :

$`m{\displaystyle \frac{d^2u_{i,j}}{dt^2}}`$ $`+`$ $`{\displaystyle \frac{du_{i,j}}{dt}}=u_{i+1,j}2u_{i,j}+u_{i1,j}`$ (1)
$`+`$ $`A[\mathrm{sin}(u_{i,j+1}u_{i,j})+\mathrm{sin}(u_{i,j1}u_{i,j})].`$
The lattice is a collection of chains in the $`x`$ direction with elastic interaction between nearest neighbors within the same chain and sinusoidal interaction between chains. $`u_{i,j}/(2\pi )`$ is the dimensionless displacement of atom $`(i,j)`$ in the $`x`$ direction measured in units of the Burgers vector length $`b`$. $`A>0`$ measures the relative strengths of the nonlinear forces exerted by atoms on different planes $`y=k`$ (constant) and the linear forces exerted within any plane $`y=k`$. The dimensionless parameter $`A`$ also determines the width of the dislocation core ($`1/\sqrt{A}`$). Lastly, the time unit is the ratio between the friction coefficient and the spring constant in the $`x`$ direction. Then $`m`$ is the dimensionless ratio between the atomic mass times the spring coefficient and the square of the friction coefficient. In dislocation dynamics, an important case is that of overdamped dynamics, $`m=0`$ gro00 . Eq. (1) can be generalized to a vector model having a displacement vector $`(u_{ij},v_{ij})`$ and a continuum limit yielding the 2D Navier equations with cubic symmetry car03 . Such model has among its solutions edge dislocations with Burgers vectors in the $`x`$ or $`y`$ directions gliding in the direction thereof (which does not have to be assumed as in the present simple model). This model can also be solved using our methods at the expense of technical complications and high computational cost.

In this geometry, the far field of a static 2D edge dislocation is approximately given by the corresponding continuum elastic displacement $`u_{ij}=u(x,y)`$ with $`x=ϵi`$, $`y=ϵj`$ (where $`ϵ=b/L1`$, $`i`$, $`j`$ are large and $`L`$ is the appropriate mesoscopic length). Then the stationary solutions of Eq. (1) satisfy the equations of anisotropic linear elasticity, $`u_{xx}+Au_{yy}=u_{xx}+u_{YY}=0`$, ($`Y=y/\sqrt{A}`$) far away from singularities and jumps kkl . The solution corresponding to the edge dislocation is the polar angle $`\theta (x,Y)[0,2\pi )`$, measured from the positive $`x`$ axis. Continuum approximations break down near the dislocation core, which should be described by the discrete model car97 . The advantadge of Eq. (1) compared to other 2D generalizations lom86  of the Frenkel-Kontorova model is that it yields the correct decay for strains and stresses: $`r^1`$ as $`r^2=x^2+Y^2\mathrm{}`$, instead of exponential decay.

Overdamped dynamics and static Peierls stress. We shall now study the structure of a static edge dislocation of Eq. (1), the critical stress needed to set it in motion and its subsequent speed. We solve numerically Eq. (1) with $`m=0`$ on a large lattice $`|i|,|j|N`$ using Neumann boundary conditions (NBC) corresponding to applying a shear stress of strength $`F`$ in the $`x`$ direction. The (far field) continuum elastic displacement for a static 2D edge dislocation subject to such a shear stress is $`\theta (x,y/\sqrt{A})+Fy`$. Then the NBC are $`u_{\pm (N+1),j}u_{\pm N,j}=\pm [\theta _{\pm (N+1),j}^A\theta _{\pm N,j}^A]`$ and $`u_{i,\pm (N+1)}u_{i,\pm N}=\pm [\theta _{i,\pm (N+1)}^A\theta _{i,\pm N}^A]\pm F`$, where $`\theta _{i,j}^A=\theta (i,j/\sqrt{A})`$ with $`\theta (0,0)=\pi /2`$. If $`F=0`$ and the initial condition is the elastic far field $`\theta _{i,j}^A`$, the system relaxes to a stationary configuration $`u_{i,j}`$. The dislocation is expected to remain stationary for $`F0`$ unless $`|F|`$ is larger than a critical value $`F_{cs}(A)`$, related to the so called static Peierls stress nab67 ; hir82 . Nonlinear stability of the stationary edge dislocation for $`|F|<F_{cs}`$ was proven in Ref. car02, . To test this picture, we solve numerically Eq. (1) in a large lattice, using NBC and the static dislocation obtained for $`F=0`$ as initial data. For large times and $`|F|`$ small, the system relaxes to a steady configuration $`u_{i,j}`$ which provides the structure of the core, see Fig. 1. When $`|F|`$ is large enough, the dislocation is observed to glide in the $`x`$ direction: to the right if $`F>0`$, and to the left if $`F<0`$.

To calculate $`F_{cs}`$, we extend the depinning calculations of Ref. carPRL01,  to 2D systems. We redefine $`u_{i,j}=U_{i,j}+Fj`$, insert $`U_{i,j}=U_{i,j}(F,A)+v_{i,j}(t)`$ in Eq. (1) with $`m=0`$ and expand the resulting equation in powers of $`v_{i,j}`$, about the stationary state $`U_{i,j}(A,F)`$ up to cubic terms. Subscripts in the resulting equation can be numbered with a single one starting from the point $`i=j=N`$: $`U_{i,j}=U_k`$ and $`v_{i,j}=v_k`$, $`k=i+(j+N)(2N+1)`$ for $`i,j=N,\mathrm{},N`$. The resulting equation can be written formally as $`d𝐯/dt=(F)𝐯+(𝐯,𝐯;F)`$, where the vector $`𝐯`$ has components $`v_k`$. The linear stability of the stationary state $`U_k(A,F)`$ depends on the eigenvalues of the matrix $`(F)`$. These eigenvalues are all real negative for $`|F|<F_{cs}`$ whereas one of them vanishes at $`|F|=F_{cs}`$. This criterion allows us to numerically determine $`F_{cs}`$ as a function of $`A`$; see Fig. 2(a). Notice that the critical stress increases with $`A`$. Thus narrow core dislocations ($`A`$ large) are harder to move.

Dislocation velocity. Let us assume that $`F>F_{cs}`$ (the case $`F<F_{cs}`$ is similar). Then $`v_{ij}=(FF_c)j+\varphi (t)r_{ij}`$ (plus terms that decay exponentially fast in time). The procedure sketched in Refs. carPRL01, ; carPRE01,  for discrete 1D systems yields the amplitude equation $`d\varphi /dt=\alpha +\beta \varphi ^2`$. Here $`\alpha =𝐥𝒩(F_{cs})(FF_{cs})/(𝐥𝐫)`$, $`𝒩_{ij}=A[\mathrm{cos}(U_{i,j+1}U_{i,j}+F_c)\mathrm{cos}(U_{i,j}U_{i,j1}+F_c)]`$ and $`\beta =𝐥(𝐫,𝐫;F_{cs})/(𝐥𝐫)`$, $`𝐥`$ and $`𝐫`$ are the left and right eigenvectors of the matrix $`(F_{cs})`$ corresponding to its zero eigenvalue (its largest one). From the amplitude equation, the approximate dislocation velocity is carPRL01 : $`c\sqrt{\alpha \beta }/\pi =O(|FF_{cs}|^{\frac{1}{2}})`$. Numerically measured and theoretically predicted dislocation velocities are compared in Fig. 2(b). Calculations in lattices of different sizes yield similar results.

How do we calculate numerically the dislocation velocity? This is an important point for using the calculated dislocation velocity as a function of stress in mesoscopic theories and a few comments are in order. If we solve numerically Eq. (1) with static NBC for $`|F|>F_{cs}`$, the velocity of the dislocation increases as it moves towards the boundary. The dislocation accelerates because we are using the far field of a steady dislocation as boundary condition, instead of the (more sensible) far field of a moving dislocation. However, the latter is in principle unknown because we do not know the dislocation speed. We will assume nevertheless that the dislocation moves at constant speed $`c`$ once it starts moving, as it would in a stressed infinite system. Then the correct dislocation far field is $`\theta (ict,j/\sqrt{A})+Fj`$. With this far field in the NBC, Eq. (1) has traveling wave solutions $`u_{i,j}(t)`$ whose velocity can be calculated self-consistently. How? By an iterative procedure that adopts as initial trial velocity that of a dislocation subject to static NBC as it starts moving. Near threshold, step-like profiles are observed (see Fig. 3), that become smoother as $`F`$ increases. The profiles have been calculated by following the trajectories of points with the same value of $`i`$ and different values of $`j`$ for $`F>F_{cs}(A)`$, according to the formula $`u_{i,j}(t)=u(\zeta ,j)`$, $`\zeta =ict`$. Notice that the wave front profiles are kinks for $`j<0`$ and antikinks for $`j0`$.

Influence of fluctuations. The original discrete model contains both damping and fluctuation terms kkl . Fluctuation terms are appreciable only near $`F_{cs}`$, and contribute an additive white noise term to the amplitude equation. Due to this term, there is a small probability for the dislocation to move even if $`|F|<F_{cs}`$ and $`m=0`$. The resulting average velocity can be estimated by observing that the potential of the corresponding Fokker-Planck equation is cubic and it has a small barrier of height proportional to $`(F_{cs}|F|)^{\frac{3}{2}}`$. Then the exponentially small velocity of the dislocation under the critical stress is the reciprocal of the mean escape time from the barrier vankampen . Provided $`(F_{cs}|F|)D`$ (where $`D`$ measures the noise strength), we have $`\mathrm{ln}c(F_{cs}|F|)^{\frac{3}{2}}/D`$.

Inertia and dynamic Peierls stress. Inertia changes the previous picture of dislocation motion in one important aspect: the dislocations keep moving for an interval of stresses below the static Peierls stress, $`F_{cd}<|F|<F_{cs}`$. On this stress interval, stable solutions representing static and moving dislocations coexist: to depin a static dislocation, we need $`|F|>F_{cs}`$. However if $`|F|`$ decreases below $`F_{cs}`$, a moving dislocation keeps moving until $`|F|<F_{cd}`$; see Fig. 2. Thus $`F_{cd}`$ represents the dynamic Peierls stress of the dislocation nab67 . Our theory therefore yields the static and the dynamic Peierls stresses and the velocity of a dislocation.

Interaction between edge dislocations. Let us assume that there are $`N`$ static edge dislocations at the points $`(x_n,y_n)`$ parallel to one dislocation at $`(x_0,y_0)`$, and that all dislocations are separated from each other by distances of order $`L1`$ (measured in units of the Burgers vector length, $`b`$). We want to analyze whether the collective influence of the $`N`$ distant dislocations can move that at $`(x_0,y_0)`$. This problem is similar to that of deriving a reduced dynamics for the centers of 2D vortices of Ginzburg-Landau equations subject to their mutual influence neu90 . In the case of dislocations, the existence of a pinning threshold implies that the reduced dynamics is that of a single dislocation subject to the mean field created by the others. We thus have a reduced field dynamics, not particle dynamics as in the case of the Ginzburg-Landau vortices.

Displacement vectors pose the problem of defining branch cuts in the continuum (elastic) limit, which leads us to consider instead the distortion tensor as a primary quantity ll7 . For our discrete system, the distortion tensor has nonzero components $`w_{i,j}^{(1)}=u_{i+1,j}u_{i,j}`$ and $`w_{i,j}^{(2)}=\mathrm{sin}(u_{i,j+1}u_{i,j})`$ that become $`ϵ(u/x)`$ and $`ϵ(u/y)`$, respectively, in the continuum limit $`ϵ=1/L0`$, $`xx_0=ϵi`$, $`yy_0=ϵj`$ finite. In the continuum limit, the distortion tensor of an edge dislocation centered at the origin has nonzero components $`w^{(1)}=ϵ\sqrt{A}y/(Ax^2+y^2)`$ and $`w^{(2)}=ϵ\sqrt{A}x/(Ax^2+y^2)`$. If we have $`N`$ edge dislocations at $`(x_n,y_n)`$, $`1nN`$, far from one at $`(x_0,y_0)`$, the distortion tensor is sum of individual contributions. Then the far field distortion tensor seen by the dislocation at $`(x_0,y_0)`$ is:

$`w_{i,j}^{(1)}={\displaystyle \frac{\sqrt{A}j}{Ai^2+j^2}}+F_1,w_{i,j}^{(2)}={\displaystyle \frac{\sqrt{A}i}{Ai^2+j^2}}+F_2,`$ (2)
$`F_1=ϵ{\displaystyle \underset{n=1}{\overset{N}{}}}{\displaystyle \frac{\sqrt{A}(y_0y_n)}{A(x_0x_n)^2+(y_0y_n)^2}}+\mathrm{},`$ (3)
$`F_2=ϵ{\displaystyle \underset{n=1}{\overset{N}{}}}{\displaystyle \frac{\sqrt{A}(x_0x_n)}{A(x_0x_n)^2+(y_0y_n)^2}}+\mathrm{}.`$ (4)
The dislocation at $`(x_0,y_0)`$ moves if $`F_2>F_{cs}(A)`$. This cannot be achieved as $`ϵ0`$ unless $`N=O(1/ϵ)`$. Then the sums in Eqs. (3) and (4) become integrals. We define a static dislocation density $`\rho (x,y)`$ as the limit of $`N^1_{n=1}^N\delta (xx_n)\delta (yy_n)`$ as $`N\mathrm{}`$. Then

$`F_1=\alpha {\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\sqrt{A}(y_0y)\rho (x,y)}{A(x_0x)^2+(y_0y)^2}}𝑑x𝑑y,`$ (5)
$`F_2=\alpha {\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\sqrt{A}(x_0x)\rho (x,y)}{A(x_0x)^2+(y_0y)^2}}𝑑x𝑑y,`$ (6)
where $`\alpha =ϵN`$ is the ratio of the total Burgers vector to the mesoscopic length measuring average interdislocation distance. As an example, let us assume that $`y_0=y_n=0`$. Then $`\rho =\rho (x)\delta (y)`$ and $`F_1=0`$. Let us assume that the dislocations are constrained by two obstacles at $`x_0=\pm l`$ and subject to the same critical stress. Then the critical dislocation density is $`\rho (x)=[1(\sqrt{A}F_{cs}/\alpha )x]/[\pi \sqrt{l^2x^2}]`$, provided $`\alpha >\sqrt{A}F_{cs}l`$ (cf. ll7, , page 127).

In conclusion, edge dislocations can be characterized as traveling waves of discrete models. The dislocation far field moves at a constant velocity equal to that of the dislocation core. Static and dynamic Peierls stresses and the dislocation velocity as a function of applied stress can be found numerically (or analytically near critical stress in the overdamped case) and adopted as the basis of a mesoscopic theory gro00 . We have also shown that the interaction between distant edge dislocations can be described in terms of a continuous dislocation density. This field-theoretical reduced description greatly contrasts with the case of interacting point vortices, which is completely described by the particle dynamics of the vortex centers neu90 . Extension to fully vectorial models and to other types of dislocation should follow along similar lines.

We thank Joe Keller for helpful discussions. This work has been supported by the MCyT grant BFM2002-04127-C02, by the Third Regional Research Program of the Autonomous Region of Madrid (Strategic Groups Action), and by the European Union under grant HPRN-CT-2002-00282.


text:
# 1 Introduction

## 1 Introduction

PSR B1929$`+`$10 is one of the closest and oldest known ordinary pulsars that has been detected in X-rays. The two latest investigations of this object in X-rays were based on ROSAT PSPC (Yancopoulos et al. 1994, hereafter Y94) and ASCA (Wang & Halpern 1997, hereafter WH97) observations with exposure durations of $`45\mathrm{ks}`$ and $`54\mathrm{ks}`$, respectively. Y94 obtained, after background subtraction, a total of $`420\pm 25`$ photons in the 0.1 - 2.0 keV band, corresponding to a luminosity of $`1.210^{30}\mathrm{ergs}\mathrm{s}^1`$ for a source distance of $`250\mathrm{pc}`$, or $`310^4`$ of the pulsar’s spin-down luminosity. After folding the barycentered arrival times of the selected events Y94 fitted the pulse shape using a sinusoid yielding a pulsed fraction of $`0.28\pm 0.10`$. The spectrum was fitted with a black-body with temperature $`T_{\mathrm{}}3.210^6\mathrm{K}`$, which indicated that the implied emitting area has a radius of less than $`50\mathrm{m}`$. WH97 used ASCA data from October 1994. The spectrum of PSR B1929$`+`$10 over the 0.5 - 5.0 keV range was fitted with a single black-body. The neutral Hydrogen column density was not included as a free parameter in the fit. It was fixed at the level $`110^{20}\mathrm{cm}^2`$, and the luminosity was calculated for a distance of $`250\mathrm{pc}`$ used by Y94. Their results yield $`T=5.14\pm 0.5310^6\mathrm{K}`$ and $`L1.2810^{30}\mathrm{ergs}\mathrm{s}^1`$. In the timing analysis WH97 extracted events from a circle of radius 2.5 arcmin for both GIS 2 and 3 and obtained 747 counts within this radius, 413 of which were estimated to belong to the background. The lightcurve was characterized by a single broad maximum, similar to the ROSAT soft X-ray profile, with a modulation significance of $`3.1\sigma `$. The pulse fraction was estimated to be $`0.35\pm 0.15`$. Both groups were unable to discriminate conclusively between a thermal and a non-thermal spectrum. While a black-body is a better overall fit, the inferred emitting area is surprisingly small.

In this study we revisited the ROSAT PSPC and ASCA analyses. Moreover, we also included ROSAT HRI data from a combination of 3 observations lasting 346 ks in total, together with data from a 30 ks RXTE exposure. The main goal was to characterize the timing and spectral properties of PSR B1929$`+`$10 in more detail using all available X-ray data.

## 2 Observations

In our study we made use of archival ROSAT PSPC/HRI, ASCA SIS/GIS and RXTE PCA data. An overview of the X-ray observations of PSR B1929$`+`$10 is presented in Table I.

## 3 Timing analysis

### 3.1 ROSAT PSPC

Point-like X-ray emission was clearly detected in the vicinity of the radio position of the pulsar on the ROSAT PSPC Maximum Likelihood Ratio map (MLR map, see e.g. Kuiper et al. 1998). The number of source and background counts are derived simultaneously by this method. In the timing analysis we extracted events within an optimum radius (the radius where the signal-to-noise ratio S/N maximizes) of $`50^{\prime \prime }`$ from the X-ray centroid. The 604 events (0.1 - 2.5 keV) falling within this radius have subsequently been barycentered and folded through appropriate pulsar parameters (see e.g. Table II for general characteristics of PSR B1929$`+`$10) to obtain the pulse-phase distribution shown in left panel of Fig. 1. The profile is characterized by one broad maximum and deviates from a statistically flat distribution at a $`4.2\sigma `$ level when applying a $`Z_1^1`$ method (Buccheri et al. 1983). Due to the limited accuracy of the ROSAT absolute timing we could not relate the X-ray phase to the radio phase. In order to quantify the pulsed fraction we have applied two methods: one based on bootstrapping (Swanepoel et al. 1996) and another based on sinusoid fitting. For broad sinusoidal profiles it is expected that the bootstrap method will systematically underestimate the genuine pulsed fraction because parts of the leading and trailing wings of the broad profile are considered to be part of the interval used to construct the unpulsed level. Applying the bootstrap method we obtain a pulsed fraction, defined as $`f_{pulsed}=\frac{N_P}{N_P+N_{DC}}`$, of $`0.24\pm 0.07`$. The number of pulsed counts $`N_P`$ follows directly from the bootstrap method, while the number of DC counts $`N_{DC}`$ requires an estimate of the number of background counts within the extraction radius of $`50^{\prime \prime }`$. From the MLR map we could estimate that $`140.5\pm 2.2`$ background counts are expected within the extraction radius.

Because the pulse profile shows one broad maximum sinusoid fitting using one harmonic provides an adequate description of the measured distribution. Applying sinusoid fitting to the 0.1 - 2.5 keV PSPC pulse profile (see Fig. 2 *left panel*) we find a pulse fraction of $`0.36\pm 0.08`$ slightly higher than the bootstrap value, as expected, and consistent with the value found earlier by Y94.

By selecting the events further on energy we applied the same sinusoid fitting method to the 0.1 - 1.28 keV and 1.28 - 2.5 keV pulse profiles. These two profiles are also shown in Fig. 1: 0.1 - 1.28 keV, *middle panel*; modulation significance $`3.0\sigma `$, and 1.28 - 2.5 keV, *right panel*, modulation significance also $`3.0\sigma `$. The pulsed fractions applying sinusoid fitting for these energy ranges are $`0.33\pm 0.09`$ and $`0.42\pm 0.13`$, respectively.

### 3.2 ROSAT HRI

For the 3 long ROSAT HRI observations (see Table I) the timing analysis consisted in extracting the events from a $`11^{\prime \prime }`$ aperture around the centroid of the PSR B1929$`+`$10 counterpart and subsequently folding the barycentered event times using an appropriate pulsar ephemeris. The obtained pulse profiles for each observation are all of low significance (between 2 - 3 sigma), although the broad bump can be discerned easily. Cross-correlating each of the three HRI profiles with the 0.1 - 2.5 keV PSPC profile and applying the shifts in the combination of the three separate aligned profiles, yielded a phase distribution which deviates from a flat distribution at the $`3.7\sigma `$ level (see Fig. 2 *middle panel*). In this case the sinusoid fitting method yielded a pulsed fraction of $`0.21\pm 0.05`$, lower but consistent with the PSPC value. The combined ROSAT PSPC and HRI profile is shown in the *right panel* of Fig. 2 and it deviates from being flat at a $`5.6\sigma `$ level. The pulsed fraction of the combined profile is $`0.25\pm 0.04`$.

### 3.3 ASCA GIS

Of the two detector systems aboard ASCA only the GIS detectors have sufficient timing resolution to support timing analyses at millisecond scales. Data from the high and medium rate telemetry mode were used. The 3.9 ms timing resolution valid for these rates was sufficiently high to sample the pulse period of 227 ms (see Table II). In the spatial analysis a Maximum Likelihood method similar to the one used for the ROSAT PSPC/HRI data, provided the optimum position for the X-ray counterpart of PSR B1929$`+`$10. An extraction radius of $`130^{\prime \prime }`$ from the optimum position was found to optimize the S/N ratio. In the 0.5 - 5 keV energy interval 493 event timetags were barycentered and subsequently folded through appropriate timing parameters to yield the corresponding pulse profile. This 0.5 - 5 keV profile is shown in Fig. 3 *left panel* and deviates from uniformity at a $`2.7\sigma `$ level. The same broad enhancement is visible as in the ROSAT profiles at softer X-ray energies. Constraining the energy to the 0.5 - 2.0 keV interval, fully overlapping the ROSAT energy window, yielded a significance of $`3.1\sigma `$ (see Fig. 3 *right panel*). The pulsed fractions derived from sinusoid fitting for the (integral) 0.5 - 5 keV and (ROSAT overlapping) 0.5 - 2.0 keV energy ranges are $`0.36\pm 0.11`$ and $`0.69\pm 0.18`$, respectively. The former value is consistent with the one found by W97, while the latter value derived for the ROSAT overlapping interval seems to confirm the increasing trend seen at soft X-ray energies for the pulsed fraction as a function of energy.

### 3.4 RXTE PCA

Finally, we analysed RXTE PCA (2 - 60 keV) data from a 30 ks observation of PSR B1929$`+`$10 (Obs. id. 20156; 21-Nov-1997) obtained in Good Xenon mode, time tagging each event with a $`0.9\mu \mathrm{s}`$ time resolution. For the timing analysis we used ftools program fasebin to select events from the top layer of each involved PCU, to barycenter the times of the selected events and finally to fold the barycentered event times with an appropriate ephemeris. The ephemeris we used in the folding procedure within fasebin is based on radio timing observations of PSR B1929$`+`$10 with the 32m TCfA radio-telescope at Toruń, Poland since the beginning of June 1997. A dual-channel, circular polarization L-band receiving system at frequencies around 1.7 GHz and a 2 x 64 x 3 MHz pulsar back-end, the Penn Stata Pulsar Machine - 2 (Konacki et al. 1999) are used. For analyzing the radio pulsar timing data the TEMPO software package (http://pulsar.princeton.edu/tempo) has been used.

The folding procedure embedded in fasebin resulted in pulse phase histograms for each of the 256 PHA channels. Unfortunately, in none of the PHA channels or channel intervals significant pulsed emission has been detected.

## 4 Spectral analysis

In the spectral analysis we first derived the number of counts assigned to the X-ray counterpart of PSR B1929$`+`$10 in several narrow energy slices using a spatial (MLR) analysis which took into account the presence of nearby sources. We applied this approach to the ROSAT PSPC-B (20 energy bins in the interval 0.1 \- 2.5 keV ) and ASCA GIS 2+3 (10 energy bins in the interval 0.5 - 10 keV) data. We chose for ASCA GIS data (contrary to W97) because the GIS sensitivity is higher than the SIS sensitivity for energies beyond $`5\mathrm{keV}`$ and because of uncertainties in the SIS efficiency for energies below 1 keV. In the combined ROSAT PSPC-B and ASCA GIS 2+3 spectral fits we used the latest upgrades for the response matrices and we took into account the vignetting corrections for the off-axis (3.7 arcmin) ASCA observation.

The spectral models we fitted to the combined data set consist of a single absorbed power-law, a single absorbed black-body, an absorbed black-body plus power-law and an absorbed double black-body. The best fits were obtained for a single absorbed power-law (Fig. 4 *left panel*) and an absorbed double black-body model (Fig. 4 *right panel*). The fit results for the latter two cases, each with a free Hydrogen column density, are shown in Table III.

In case of the double black-body model we are dealing with two different thermal X-ray components with temperatures $`(2.0_{0.05}^{+0.05})10^6\mathrm{K}`$ and $`(6.9_{0.35}^{+0.23})10^6\mathrm{K}`$ for $`T_1`$ and $`T_2`$, respectively. For the characteristic age of PSR B1929$`+`$10 ($`\tau 3.1\mathrm{Myr}`$) cooling models predict a surface temperature of the neutron star to be $`\text{ }<10^5\mathrm{K}`$, much too low to be responsible for the observed X-ray spectrum of PSR B1929$`+`$10. However, a plausible explanation could be that the observed X-ray radiation originates from a heated polar cap characterized by two different thermal components (see e.g. Cheng & Zhang 1999). The N$`_\text{H}`$ value obtained from the fit is about two times higher than that used by Y94 and WH97. Moreover, the (non-thermal) power-law fit resulting in a photon index of $`2.54\pm 0.12`$ yielded an estimate for the neutral Hydrogen density towards the source of $`9.8_{1.0}^{+1.4}10^{20}\mathrm{cm}^2`$. This value is six times larger than the adopted one in earlier works (Y94, WH97). In both cases the derived column density is consistent with the updated (radio) distance determined by Brisken et al. (2002), and the Hydrogen distribution measured in the direction of the pulsar (Frisch & York, 1983). Additionally, as Pavlov et al. (1996) mentioned, values of N$`_\text{H}`$ obtained from direct measurements of stars in the neighbourhood of the pulsar are rather controversial. They differ from $`<10^{19}`$ up to $`10^{21}\mathrm{cm}^2`$. This suggests that interstellar medium is very patchy in that direction.

## 5 Conclusion

Our revisited analysis of X-ray data of PSR 1929$`+`$10 from ROSAT, ASCA and RXTE observations yielded an improved timing signal at soft X-rays, $`5.6\sigma `$ for the combined ROSAT PSPC/HRI pulse profile, and a non-detection for energies beyond $`2`$ keV (RXTE PCA). Analyzing the pulse profile at soft X-rays (ROSAT PSPC) in two different energy bands suggested that the pulsed fraction increases from $`0.33\pm 0.09`$ in the 0.1 - 1.28 keV band to $`0.42\pm 0.13`$ in the 1.28 - 2.5 keV band (even $`0.69\pm 0.18`$ for the 0.5 - 2 keV ASCA GIS band). This tendency must stop abruptly near $`2`$ keV to be consistent with the non-detection of pulsed emission above $`2`$ keV. A spectral analysis of the combined ROSAT/ASCA data in the 0.1 - 10 keV energy range indicated that both an absorbed power-law (non-thermal) model and an absorbed double black-body model provide an adequate description of the observed X-ray spectrum. The derived Hydrogen column density N$`_\text{H}`$ is in both cases much larger than the value used in previous X-ray studies. However, the newly determined N$`_\text{H}`$ values are more consistent with the updated (increased) distance to PSR B1929$`+`$10 in combination with the measured Hydrogen distribution in the direction of the pulsar.

A scheduled observation with the XMM-Newton satellite could confirm our results and the much better statistics would help to fix the column density N$`_\text{H}`$ and to discriminate conclusively between thermal and non-thermal models for the X-ray emission from this pulsar.

## Acknowledgements

This research was supported by KBN grant 2P03D02117 and Polish Foundation of Astronomy. Aga Woźna acknowledges people from SRON for their hospitality during her stay there, when most of this work was done.


text:
# Connectivity Distribution of Spatial Networks

## I Introduction

In contrast to abstract graphs, many real networks are embedded in a metric space: The interactions between the nodes depend on their spatial distance and usually take place between nearest neighbors. Examples of such networks are transportation and communication networks, friendship or contact networks. An especially important example is the Internet, which is a set of routers linked by physical cables with different lengths. Several recent studies have investigated networks whose nodes are embedded in a metric space, and where the probability of connecting two nodes with an edge depends on their distance .

On the other hand, the concept of scale-free network has emerged in the last few years as a powerful unifying paradigm in the study of complex systems of natural, technological and social origin (see Refs. ). It is therefore natural to investigate the possibility of embedding scale-free networks in space. In particular in Refs. , the general problem of embedding a scale-free network of given connectivity distribution in a Euclidean lattice was studied.

In this paper we take a somewhat reversed point of view: Our starting point is a spatial distribution of points on a continuous manifold $`M`$. Such points are the nodes of a network built by joining two nodes whenever their distance is less than a certain cutoff. In the following we will call ‘spatial networks’ those obtained by this procedure (sometimes also termed as ‘ad hoc networks’ ). We study how the connectivity distribution depends on the distribution of the nodes on $`M`$. The fact that the nodes live in a continuous manifold rather than a lattice is important in generating scale-free networks since the number of neighbors that a node can have within a certain distance is not limited a priori by the lattice structure.

Spatial networks originating from a uniform distribution of the nodes were studied in Ref. , where it was shown that while the connectivity distribution takes the same Poisson form as in the Erdos-Renyi random networks (see eg. ), other important features, most notably the clustering coefficient, are radically different from the Erdos-Renyi case. The formation of giant clusters in such networks was studied in Ref. .

A natural application of this class of networks is the study of epidemics propagation. While scale-free networks in which the connectivity does not depend on a pre-existing metric structure do not display an epidemic threshold , the situation changes when geographical closeness of two nodes influences their probability to be connected . In our networks geographical closeness completely determines whether two nodes will be connected, so that we expect to find an epidemic threshold. Spatial networks appear to be a rather realistic model for epidemic propagation in animal or plant populations, where we do not expect any individual to have an interaction range very different from the average (as is the case for highly mobile human populations), while we do expect the population density to be non-uniform.

In addition, spatial networks of biological origin have also recently been constructed, especially from gene expression data obtained from micro-array experiments . The connectivity distribution of such networks turns out to be a truncated power-law with the exponent of the power-law decay often close to unity. Interestingly, also networks constructed from gene expression data by a different rule, which does not satisfy our definition of a spatial network, display a similar connectivity distribution.

These facts motivated us to study how the connectivity distribution of a spatial network depends on the distribution of the nodes in space and under which conditions a scale-free network can be obtained, and finally whether an exponent close to unity plays any special role in this context.

## II The connectivity distribution of a spatial network

### A General Expression

The $`N`$ nodes of the network are supposed to be in a $`D`$-dimensional space and we will assume that they are distributed randomly in space with density $`p(x)`$. Given a node chosen at random, the probability that it is placed within a given domain of space $`m`$ is

$$_m𝑑xp(x)$$
(1)
(we denote the integration measure by $`dx`$ independently of the dimension $`D`$). Once the nodes are distributed in this space, we have to construct the edges. We will use a simple cut-off rule: Given two nodes $`i`$ and $`j`$, located at $`x_i`$ and $`x_j`$ respectively, an edge will join them if

$$d(x_i,x_j)R,$$
(2)
where $`d(x,y)`$ is the distance between $`x`$ and $`y`$. Therefore, once the nodes have been distributed the network is completely determined by the choice of the cutoff $`R`$. This model follows strictly the rule used to construct networks based on gene expression data in Refs. , but more generally, it can be used to model the case where the interaction has a typical scale given by $`R`$.

Denoting by $`B_R(x)`$ the ball of radius $`R`$ centered in $`x`$

$$B_R(x)\{yM:d(x,y)<R\},$$
(3)
the probability that a given node is placed within $`B_R(x)`$ is

$$q_R(x)=_{B_R(x)}𝑑x^{}p(x^{}).$$
(4)
If we consider a node located at $`x`$, the probability that it will have $`k`$ neighbors is then just the probability that $`k`$ additional nodes are located in the $`B_R(x)`$. The connectivity distribution for a node placed in $`x`$ is thus given by the binomial distribution

$$P(k;x,R)=\left(\genfrac{}{}{0pt}{}{N1}{k}\right)q_R(x)^k\left[1q_R(x)\right]^{N1k}.$$
(5)
In the following, we will be concerned with the limit $`N\mathrm{}`$: In order to obtain a well defined connectivity distribution in this limit, one has to take the limit $`R0`$, so as to ensure that the product

$$NV(R)$$
(6)
where $`V(R)`$ is the volume of the ball

$$V(R)=_{B_R(x)}𝑑x,$$
(7)
tends to a finite constant $`\alpha `$. This means that $`R`$ must scale as $`N^D`$ as $`N\mathrm{}`$, and implies that the expected number of nodes found within the ball $`B_R(x)`$ remains finite in this limit:

$$Nq_R(x)\alpha p(x).$$
(8)
The constant $`\alpha `$ fixes the scale of the average connectivity $`k`$ of the network. Indeed, from Eq.(11) derived below, it is easy to obtain the following relation between $`\alpha `$ and $`k`$

$$k=\alpha 𝑑xp^2(x)$$
(9)
Let us note that although the connectivity distribution is well-defined for any value $`\alpha `$, it has been shown in the case of a uniform density that the existence of a giant connected component implies a minimum value of $`\alpha `$ (or equivalently $`k`$). We will not address this problem here but we can expect that for other densities there will be some similar conditions on $`\alpha `$ for the existence of a giant component.

In the limit $`N\mathrm{}`$, $`R0`$, the connectivity distribution for a node located at $`x`$ is Poissonian and Eq.(5) becomes

$$P(k;x,\alpha )=\frac{1}{k!}\alpha ^kp^k(x)e^{\alpha p(x)}.$$
(10)
For the whole space, the connectivity distribution of the network is then obtained as the spatial average of the former expression and is

$$P(k;\alpha )=\frac{\alpha ^k}{k!}𝑑xp^{k+1}(x)e^{\alpha p(x)}.$$
(11)
This formula solves the general problem of determining the connectivity distribution of a spatial network from the spatial distribution of the nodes. For a uniform distribution $`p(x)`$, we recover a Poissonian connectivity distribution as in Ref. . However, Eq. (11) shows that other connectivity distributions can be obtained. In the following section we will determine under which condition a spatial node distribution generates a scale-free network instead.

### B Scale-free spatial networks

Equation (11) allows us to determine the condition that must be satisfied by the node density $`p(x)`$ for the corresponding network to be scale-free. A network is scale-free if the moments $`k^\nu `$ of its connectivity distribution diverge for $`\nu `$ larger than a certain $`\nu _{\mathrm{max}}`$. It is easy to compute the moments of the connectivity distribution of a spatial network from Eq. (11): For integer $`\nu `$ we have

$`k^\nu _\alpha `$ $`=`$ $`{\displaystyle \underset{k}{}}\alpha ^k{\displaystyle \frac{k^\nu }{k!}}{\displaystyle _M}𝑑xp^{k+1}(x)e^{\alpha p(x)}`$ (13)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{\nu }{}}}{\displaystyle \underset{km}{}}{\displaystyle \frac{\alpha ^k}{k!}}𝐒_\nu ^{(m)}(k)_m{\displaystyle 𝑑xp^{k+1}(x)e^{\alpha p(x)}}`$ (14)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{\nu }{}}}𝐒_\nu ^{(m)}\alpha ^m{\displaystyle 𝑑xp^{m+1}(x)},`$ (16)
where $`(k)_m=k(k1)\mathrm{}(km+1)`$ and $`𝐒_\nu ^{(m)}`$ are the Stirling numbers of the second kind (see e.g. Ref. , p. 824), defined by

$$x^\nu =\underset{m=0}{\overset{\nu }{}}𝐒_\nu ^{(m)}x(x1)\mathrm{}(xm+1).$$
(17)
It follows that $`k^\nu _\alpha `$ exists if and only if the integrals

$$H_a𝑑xp^{a+1}(x)$$
(18)
exist for all $`a\nu `$. Conversely, the network is scale-free if there exists a $`\nu _{\mathrm{max}}`$ such that $`H_a`$ diverges for all $`a\nu _{\mathrm{max}}`$. The integral Eq. (18) is a measure of the information contained in the probability distribution $`p(x)`$, and is simply related to the Renyi entropy

$$R_q=\frac{1}{1q}\mathrm{log}H_{q1}.$$
(19)
### C Classes of spatial networks

More generally, Eq. (11) allows us to determine the type of spatial network obtained for a given spatial node distribution. The networks can essentially be distinguished by the decay of the connectivity distribution for large connectivities. For spatial networks, large connectivities are obtained in high density regions, namely maxima of $`p(x)`$. For the sake of simplicity, we will limit the discussion to the case of an isotropic distribution $`p(x)=p(r)`$ where $`r`$ is the modulus of $`x`$. We will also suppose that we have one density maximum located at $`r=0`$. We will then distinguish the following two cases:

* $`p(0)=p_0`$ is finite and decays sufficiently rapidly so that all the quantities $`H_a`$ are finite. This could be for instance the case for population density which is decreasing exponentially from the city center. In this case, an asymptotic evaluation of the integral in Eq. (11) shows that for large $`k`$ the connectivity distribution decays as

$$P(k;\alpha )\frac{(\alpha p_0)^k}{k!k^D}.$$
(20)
As expected, the low density fluctuations are reflected in the fast decay of the connectivity distribution and the corresponding spatial network will be of the ‘exponential’ type (ie. the connectivity distribution decreases at least as fast as an exponential).
* $`p(r0)r^\beta `$ with $`\beta <D`$ (if $`\beta >D`$ a cut-off is needed to normalize $`p(r)`$ and we are in the first situation where the maximum of $`p(r)`$ is finite). In this case, the information measure $`H_a`$ will diverge for $`a\frac{D}{\beta }1`$ and the connectivity distribution is a power-law: Its large-$`k`$ behavior is given by (see section III)

$$P(k;\alpha )k^{D/\beta }.$$
(21)
The large density fluctuations allow here for the existence of nodes with very large connectivities and the corresponding network is scale-free.

Even if it might appear unlikely that spatial densities behave pathologically around some points, this is actually the case in many instances where the nodes live in an abstract space, with a distance defined by correlations. Such examples are obtained in the case of gene expression networks (or the stock market) for which the distance is defined in terms of Pearson correlation coefficient between nodes. The spatial network has then a simple meaning: Two nodes are connected if their correlation is high enough. It has been observed that in this example, the connectivity distribution is a truncated power-law with exponent of order $`1`$.

### D The dependence on the average connectivity

In this section we consider the dependence of the connectivity distribution (11) of a spatial network on the parameter $`\alpha `$ that, as discussed in Sec. 2, is a measure of the average connectivity. From Eq.(11) one immediately obtains the following equation governing the dependence of $`P(k;\alpha )`$ on $`k`$:

$$\frac{dP(k;\alpha )}{d\alpha }=\frac{1}{\alpha }\left[kP(k;\alpha )(k+1)P(k+1;\alpha )\right].$$
(22)
Note that while all spatial networks obey this equation, the converse is not true: One can easily construct connectivity distributions $`P(k;\alpha )`$ satisfying Eq. (22) that cannot be obtained from a spatial network.

Inspection of Eq. (22) shows that the limiting case $`P(k)k^1`$ corresponds to a fixed point, where the connectivity distribution does not depend on $`\alpha `$. This observation might be relevant in explaining the common appearance of scale-free networks with connectivity exponent $`1`$ constructed as spatial networks from gene expression data .

## III A Simple Example of a Scale-free Spatial Network

In this section we discuss an explicit example of a one-dimensional distribution of nodes that generates a scale-free network. We calculate exactly the connectivity distribution in the limit of infinite number of nodes using Eq. (11). We also address the issue of finite-size effects and we compare the exact result for infinite size to numerical simulation for finite networks.

### A Connectivity distribution for an infinite network

We consider the space to be the open interval $`(0,1)`$ and the node distribution to be

$$p(x)=(1\beta )x^\beta ,\beta <1.$$
(23)
The information measure $`H_a`$ is given by

$$H_a=(1\beta )^{a+1}_0^1𝑑xx^{\beta (a+1)},$$
(24)
and diverges for

$$a\frac{1}{\beta }1.$$
(25)
Therefore we expect a scale-free network with $`k^\nu _\alpha `$ divergent for $`\nu \frac{1}{\beta }1`$.

The connectivity distribution can be computed explicitly from Eq. (11) and is

$`P(k;\alpha )`$ $`=`$ $`{\displaystyle \frac{\alpha ^k}{k!}}(1\beta )^{k+1}{\displaystyle _0^1}𝑑xx^{\beta (k+1)}e^{\alpha (1\beta )x^\beta }`$ (27)
$`=`$ $`{\displaystyle \frac{1}{\alpha \beta }}\left[\alpha (1\beta )\right]^{1/\beta }{\displaystyle \frac{\mathrm{\Gamma }[k+1\frac{1}{\beta },\alpha (1\beta )]}{\mathrm{\Gamma }(k+1)}},`$ (28)
where $`\mathrm{\Gamma }(x,a)`$ is the incomplete Gamma function

$$\mathrm{\Gamma }(x,a)=_a^{\mathrm{}}𝑑tt^{x1}e^t.$$
(29)
Since for $`z\mathrm{}`$

$$\frac{\mathrm{\Gamma }(z,a)}{\mathrm{\Gamma }(z)}1,$$
(30)
we have for $`k\mathrm{}`$

$`P(k;\alpha )`$ $``$ $`{\displaystyle \frac{1}{\alpha \beta }}\left[\alpha (1\beta )\right]^{1/\beta }{\displaystyle \frac{\mathrm{\Gamma }\left(k+1\frac{1}{\beta }\right)}{\mathrm{\Gamma }(k+1)}}`$ (31)
$``$ $`{\displaystyle \frac{1}{\alpha \beta }}\left[\alpha (1\beta )\right]^{1/\beta }k^{1/\beta }.`$ (32)
This result shows explicitly that scale-free networks with any value of the connectivity exponent down to 1 can indeed be obtained as a spatial networks, as it is observed in gene expression networks. Finally, we remark that the transition to a power-law behavior of the connectivity distribution (31) depends on $`\alpha `$ and $`\beta `$. Indeed, when the second argument $`\alpha (1\beta )`$ of the incomplete gamma function tends to zero, we have a power-law behavior for $`P(k;\alpha )`$ even at small $`k`$, as can be seen from Fig. (1).

### B Finite-size effects: Numerical results

Since real spatial networks contain a finite number of nodes and our analytical results were obtained in the limit of infinite networks, it is important to address the issue of finite size effects. In this section, we approach the problem from a numerical point of view by constructing finite spatial networks by Monte Carlo methods and comparing their connectivity distribution with the theoretical predictions.

The figure $`1`$ shows the result of such comparison for the one-dimensional spatial networks studied in the previous subsection, at $`\beta =0.5`$ and $`\beta =0.9`$, and $`\alpha =5`$. For a finite network, $`\alpha `$ is naturally defined as $`2NR`$, where $`N`$ is the number of nodes and $`R`$ is the distance cutoff used to define links. In each part of the figure, the connectivity distributions for $`N=1000`$ and $`N=20000`$ points are superimposed to the theoretical distribution given by Eq. (28). We can conclude that already for moderately sized networks the connectivity distribution is very close to the $`N\mathrm{}`$ behavior Eq. (28), thus indicating that our results provide a good approximation to the connectivity distribution of finite spatial networks.

## IV Discussion

We have presented a systematic analysis of the connectivity distribution of spatial networks constructed by joining nodes closer to each other than a cutoff distance, in the limit in which the number of nodes tends to infinity. The main results of our analysis can be summarized as follows. First, the connectivity distribution can be expressed as a functional of the spatial distribution of the nodes, Eq. (11). The moments of the connectivity distribution are related to a certain measure of information $`H_a`$ of the node distribution through Eq. (16). In particular, scale-free networks arise from this construction whenever the information measure $`H_a`$ diverges for some $`a`$ and we showed that scale-free networks with any exponent $`\gamma >1`$ can be constructed as spatial networks. Our results were obtained in the limit of infinite number $`N`$ of nodes but numerical results suggest that they provide an excellent approximation of real, finite spatial networks already for moderate $`N`$. Finally, the analysis of the dependence of the connectivity distribution $`P(k;\alpha )`$ on the connectivity scale $`\alpha `$ suggests that the limiting case $`P(k)k^1`$ corresponds to a fixed point, a fact that might be related to the empirical observation that several spatial networks constructed from gene expression data show a scale-free connectivity distribution with exponent close to 1.

Intuition suggests that spatial networks will be highly clustered and highly assortative. The clustering coefficient of spatial networks was studied in Ref.  for a uniform node distribution: However since the clustering coefficient is a local quantity and does not depend on the spatial density of nodes, the results of  should hold unchanged for all spatial networks, as long as the node distribution is isotropic in space. Moreover, spatial networks can be expected to display a high degree of assortativity, since high connectivity nodes are placed in high density regions and are therefore more likely to be connected to other high connectivity nodes. Indeed empirical spatial networks constructed from gene expression data show a high degree assortativity . A closely related issue is the determination of the diameter of spatial networks. Due to the embedding in a metric space, we expect the diameter to grow as a power of the number of nodes, and therefore we do not expect spatial networks to belong to the small-world networks class.

Perhaps the most interesting open problem is to clarify the role of the limit case $`P(k)k^1`$, namely to classify the node distributions that flow to this fixed point in some limit. An example is given by Eq.(28) in the $`\beta 1`$ limit. $`P(k)k^1`$ implies that the measure of information $`H_a`$ diverges for all $`a`$, which in turn implies that the node distribution must become degenerate (that is its support must shrink to zero measure). However this is not a sufficient condition, since for example a Gaussian $`p(x)`$ in the limit of zero variance, while becoming degenerate, does not flow to $`P(k)k^1`$.

Acknowledgments: One of us (MB) thanks the department of physics-INFN in Torino for its warm hospitality during the time this work was completed.


text:
# Follow-up of the 9C survey: Initial results

## 1 Introduction

It is now widely believed, as evidenced by many contributions in these proceedings, that the onset of AGN activity in galaxies is closely linked to vigorous star-formation, possibly in the presence of mergers. However, the precise mechanism that causes radio sources to trigger remains the subject of lively debate—see, e.g.,2000AAS…19710908P  versus 1995AAS…186.1502W . To advance our understanding of this issue it is clearly necessary to assemble samples of radiosources at the earliest possible times in their evolution.

The Gigahertz Peaked Spectrum (GPS) and Compact Steep Spectrum (CSS) sources are now thought to be the progenitors of the larger, older members of the population. These young sources, often morphologically classified as Compact Symmetric Objects (CSOs), have indeed in a few cases been observed growing at sub-relativistic speeds by VLBI 2000ApJ…541..112T . Unfortunately it has been difficult to obtain flux-limited samples of GPS and CSS sources because large-area surveys at frequencies higher than 5 GHz have been unavailable. Now, with the advent of cm-wavelength Cosmic Microwave Background (CMB) experiments, such surveys are being carried out to implement foreground source subtraction. The Cambridge Ryle Telscope (RT) is undertaking a survey at 15 GHz, the 9C survey (7), of CMB fields observed by the Very Small Array (VSA). Using 9C we have selected 156 sources on the basis of their 15 GHz flux alone; these have all been followed up at frequencies from 1.4 GHz to 43 GHz. Radio frequency data reduction is now complete for two of the three sample fields and we present results for these 111 sources. This sub-sample is complete to 25 mJy and covers a total area of 140 square degrees.

## 2 Radio and optical observations and analysis

The first three regions of 9C coincide with the VSA compact array fields, at $`00^h20^m`$ $`+30^{}`$, $`09^h40^m`$ $`+32^{}`$ and $`15^h40^m`$ $`+43^{}`$. Our sample is complete to 25 mJy in the $`00^h`$, the $`15^h`$ field and for one-quarter of the $`09^h`$ field. The remaining three-quarters of the $`09^h`$ field is complete to 60 mJy. Continuum snapshot maps were made at 1.4, 4.8, 22 and 43 GHz with the VLA and at 15 GHz with the RT between 2001 January and 2002 January. All observations for each region were made simultaneously to avoid problems introduced by source variability.

We have used the radio fluxes to classify each source according to the spectral index, $`\alpha `$ between 1.4 and 4.8 GHz (we take $`S\nu ^\alpha `$). We define three spectral classes: steep spectrum sources with $`\alpha _{4.8}^{1.4}0.5`$, flat spectrum sources with $`0.1\alpha _{4.8}^{1.4}<0.5`$, and a third class containing those sources with $`\alpha _{4.8}^{1.4}<0.1`$. The last of these includes the GPS sources displaying a peak in flux at 4.8, 15 or 22 GHz, and those still rising at 43 GHz, which are assumed to be extreme versions of the GPS sources. We also classify objects as compact if they have a radio size $`2^{\prime \prime }`$, and as extended otherwise. Of the 111 sources here, 55% are compact and 45% are extended in the radio. The numbers of sources in each class are summarised in Table 1.

It is already well established for samples selected at lower radio frequencies that extended radio sources have steeper spectra than compact sources (see, e.g. Peacock and Wall P1 , PW hereafter), and the trend continues for this sample. Plotting the distribution of spectral index for the two size classes (Figure 1) clearly shows this tendency: the median values of $`\alpha _{4.8}^{1.4}`$ are 0.21 and 0.76 for the compact and extended radio sources respectively.

To examine how useful our selection process is as a means of finding GPS sources we compare our sample to PW, who studied a sample of bright ($`S_{2.7}1.5`$ Jy) radio sources selected at 2.7 GHz. By combining our flat and peaked sources into the same class (i.e $`\alpha <0.5`$) we can directly compare our sample with PW. The percentages of sources falling into each class in the two samples are shown in table 2. Our sample contains a lower fraction of steep-spectrum sources than PW, but roughly the same fraction of the compact sources (note, though, that PW calculate $`\alpha `$ between 2.7 and 5 GHz, and so will be more likely to classify sources with standard convex spectra as steep-spectrum).

PW did not make any distinction between flat-spectrum and GPS sources, but O’Dea O1  summarised the frequency of occurence of GPS and compact steep spectrum (CSS) sources for flux-limited samples taken at around 5 GHz. Using the values given by O’Dea and K1 , along with the values from PW, the numbers of extended steep spectrum, CSS, flat spectrum and GPS sources a sample selected at $`5`$ GHz should contain can be estimated. These approximate percentage values, and those of this work, are given in table 3. We find roughly twice as many GPS sources as in samples selected at $`5`$ GHz, but only about half as many CSS sources.

A programme of optical imaging and spectroscopy is underway to obtain identifications and redshifts for the compact sources. To date (2003 January) we have obtained imaging with the Palomar $`60^{\prime \prime }`$ telescope of roughly two-thirds of the sample, including all targets fainter than $`R_{\mathrm{AB}}=19.2`$. While we caution that this work is still incomplete, we note that the initial results suggest that many of the compact sources have close companions or disturbed morphology at optical wavelengths. We find that almost all the low-$`z`$ host galaxies imaged to date have faint companions within $``$ 10 arcsec. As we move to more distant host galaxies, where our imaging programme is thus far insufficiently deep to detect very faint companions, several noticeably turbulent environments have been discovered (Fig. 2).

It is clearly now essential, in the light of the debate over major versus minor interactions as the trigger for radiosource activity, that we quantify these optical observations for the complete sample and compare with control samples of both larger radiosources and quiescent ellipticals.

## 3 Discussion and summary

In the standard picture of GPS sources, the energy density of the synchrotron plasma in the lobes is sufficiently high that they remain optically thick to higher radio frequencies K2 , with the spectral peak indicating the frequency at which the lobe plasma becomes optically thin. As the source grows to kpc size, the peak in the spectrum remains, but falls to MHz frequencies and below, and a CSS source is observed.

Almost all previously known GPS sources have peak frequencies of order 1—5 GHz. In our sample, however, we have found a significant number of sources which peak near the selection frequency. Of the 21 peaked spectrum sources, 13 peak at 15 GHz or above (three are still rising at 43 GHz). Two examples of these GPS spectra are shown in Figure 3. If the known trend continues—that the peak frequency increases as sources are seen younger and younger—we speculate that these very high frequency peakers are likely to be very young sources indeed. Efforts are underway to develop models which track the dynamical evolution and synchrotron spectrum to very early times.

We have found that selecting at 15 GHz is roughly twice as efficient at picking out GPS sources than at $`5`$ GHz, which makes it a valuable method of creating flux-limited samples of such objects. The 9C survey, and other forthcoming cm-wave surveys, are thus useful not only for studying the radio source population at high radio frequencies as a whole, but also as a means of singling out what are potentially the very youngest sources for further study. Our optical imaging is now uncovering evidence suggesting that a significant fraction of the host galaxies may be interacting. In the near future we hope to quantify these processes with control samples of both larger radiosources and non-AGN ellipticals. Plans are also in place to select a sample at a fainter flux limit, so that we may obtain access to sources matched in luminosity over a significant redshift range. We are hopeful that detailed studies of this sample and similar ones selected at high frequency will allow us both to model precisely the earliest stages of radiosource evolution, and to determine the exact causal connection between galaxy-galaxy interactions and the triggering of the radiosource.

### Acknowledgements

RCB thanks PPARC for a research studentship. The Ryle Telescope is funded by PPARC. The VLA is operated by the National Radio Astronomy Observatory, which is a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc.