text:
# Short-range type critical behavior in spite of long-range interactions: the phase transition of a Coulomb system on a lattice

## Abstract

One- to three-dimensional hypercubic lattices half-filled with localized particles interacting via the long-range Coulomb potential are investigated numerically. The temperature dependences of specific heat, mean staggered occupation, and of a generalized susceptibility indicate order-disorder phase transitions in two- and three-dimensional systems. The critical properties, clarified by finite-size scaling analysis, are consistent with those of the Ising model with short-range interaction.

The problem under which conditions Coulomb glasses exhibit genuine phase transitions has been under controversial debate for two decades DLR82 ; V93 ; GY93 ; VS94 ; Detal00 . The significance of the kind of static disorder is unclear yet GY93 ; VS94 . In this context, the relation to the Ising model with short-range interaction is of great interest. It would be very useful to know how replacing its nearest-neighbor coupling by an “antiferromagnetic” long-range Coulomb interaction modifies the critical behavior of a system without static disorder. Due to the interplay of long-range interaction and frustration, one might alternatively expect this model to belong to different universality classes: In case of a sufficiently slowly decaying ferromagnetic power-law interaction, including the $`1/r`$-interaction, mean field behavior is obtained LB02 . But the condensation in a three-dimensional model of an electrolyte exhibits Ising-like critical behavior, possibly due to efficient screening Letal02 .

To decide the question, we numerically study lattices half-filled with localized particles interacting via the long-range Coulomb potential,

$$H=\frac{1}{2}\underset{ij}{}\frac{(n_i1/2)(n_j1/2)}{\left|𝐫_i𝐫_j\right|}.$$
(1)
Here $`n_\alpha \{0,1\}`$ denote the occupation numbers of states localized at sites $`𝐫_i`$ within a $`d`$-dimensional hypercube of size $`L^d`$. Elementary charge, lattice spacing, dielectric and Boltzmann constants are all taken to be 1. Neutrality is achieved by background charges -1/2 at each site.

For reducing finite-size effects, we impose periodic boundary conditions for $`d=1`$ and 2, using the minimum image convention Metro53 . For $`d=3`$, the same approach would give rise to an unphysical feature: The groundstate would be a layered arrangement of charges instead of the expected NaCl structure in case $`L`$ is a multiple of 4 Metal01 . Hence, for $`d=3`$, we consider the sample to be surrounded by eight equally occupied cubes.

Our numerical investigations have been performed by means of algorithms based on the Metropolis procedure Metro53 ; Metal01 ; MT97 . Such simulations are very expensive not only for the long-range interaction, but also for correlation times diverging close to phase transitions and as temperature $`T`$ vanishes. Therefore it is necessary to adapt the dynamics to the situation that the simultaneous change of $`n_i`$ for certain clusters is necessary for overcoming passes in the energy landscape. To our knowledge, a procedure similar to the Swendsen-Wang algorithm Swen.Wang  is not available for the interaction type considered here. Hence we modified the dynamics “by hand” to take into account several kinds of processes: one-electron exchange with the surroundings, one-electron hops over ($`T`$ dependent) restricted distance, and two-electron hops simultaneously changing the occupation of four neighboring sites.

At high $`T`$, we use the original Metropolis method Metro53 . But at low $`T`$, we take advantage of a hybrid procedure much accelerating the computations. It connects the direct evaluation of weighted sums over states within a low-energy subset of the configuration space with Metropolis sampling of the complementary high-energy subset MT97 .

We now turn to the qualitative behavior of specific heat, order parameter, and susceptibility. The specific heat $`c`$ was obtained from energy fluctuations utilizing $`c=(H^2H^2)/(T^2L^d)`$. The $`T`$ and $`L`$ dependences of $`c`$ are presented in Fig. 1 for dimensions $`d=1`$ to 3. This graph shows that, for $`d=2`$ and 3, sharp peaks of increasing height evolve within a small $`T`$ region as $`L`$ grows. Away from the peaks, within the $`T`$ intervals presented, $`c`$ is almost independent of $`L`$. But for $`d=1`$, there are only broad, rounded peaks with $`L`$ independent height – a logarithmic $`T`$ scale is used for $`d=1`$, in contrast to the linear scales for $`d=2`$ and 3 which display far smaller $`T`$ intervals. For $`d=1`$, finite size effects are restricted to low $`T`$ where the reliability bound decreases with $`L`$.

Hence, according to the behavior of $`c(T,L)`$ a phase transition likely occurs for $`d=2`$ and 3, in agreement with results of lattice gas simulations for $`d=3`$ Letal02 ; DS99 . However, for $`d=1`$, in spite of the long-range interaction, there seems to be no phase transition at finite $`T`$.

Analogously to an antiferromagnet, the order inherent in a charge arrangement $`n_i`$ can be characterized by means of the staggered occupation $`\sigma _i`$ relating to a NaCl structure. For example, if $`d=3`$, it is given by

$$\sigma _i=(2n_i1)(1)^{x_i+y_i+z_i}$$
(2)
where $`x_i`$, $`y_i`$, and $`z_i`$ denote the (integer) components of $`𝐫_i`$. Thus we measure the mean absolute value of the staggered occupation $`|\sigma |`$ as order parameter.

$`T`$ and $`L`$ dependences of $`|\sigma |`$ are shown in Fig. 2. For $`d=1`$, a rapid decrease of $`|\sigma |`$ with increasing $`T`$ occurs already clearly below the temperature of maximum $`c`$ (same $`T`$ scales in Figs. 1 and 2). This marked decrease shifts to lower $`T`$ with increasing $`L`$. For $`d=2`$, a qualitatively different behavior is found: $`|\sigma |`$ decreases rapidly just in that $`T`$ region where the peak of $`c(T)`$ evolves. The $`T`$ interval of rapidly diminishing $`|\sigma |`$ shrinks as $`L`$ rises, another indication of the phase transition. This interpretation is confirmed by an extrapolation $`L\mathrm{}`$: Assuming $`|\sigma |(T,L)=|\sigma |(T,\mathrm{})+A(T)/L^{d/2}`$ with $`A`$ independent of $`L`$, an almost sharp transition is obtained from data for $`L=34`$ and $`58`$, although this extrapolation is of limited accuracy in the immediate vicinity of the transition. For $`d=3`$, the results (not shown here) qualitatively resemble our findings for $`d=2`$.

The generalized susceptibility $`\chi `$ related to the order parameter $`|\sigma |`$ is given by

$$\chi =L^d(\sigma ^2|\sigma |^2)/T.$$
(3)
Fig. 3 shows $`T`$ and $`L`$ dependences of $`\chi `$. It confirms the conclusions drawn from $`|\sigma |(T,L)`$: On the one hand, for $`d=1`$, a broad peak of $`\chi (T)`$ evolves with increasing $`L`$ where $`T_{\mathrm{max}}`$, the temperature of maximum $`\chi `$, decreases. As the inset demonstrates, $`T_{\mathrm{max}}`$ can be approximated by $`a/\mathrm{ln}(bL)`$ with constants $`a`$ and $`b`$. Hence, $`T_{\mathrm{max}}`$ likely vanishes as $`L\mathrm{}`$ so that there seems to be no phase transition for $`d=1`$ at finite $`T`$. On the other hand, for $`d=3`$, as $`L`$ rises, a narrow peak grows in just that $`T`$ region where $`c(T,L)`$ has such a feature, a further indication of the phase transition. For $`d=2`$, $`\chi (T,L)`$ behaves qualitatively similar to the results for $`d=3`$.

The quantitative evaluation of our simulation data consists in a finite-size scaling analysis Betal00 ; G92 . For this aim, we first consider $`q_2=\mathrm{ln}(1\sigma ^2^2/\sigma ^4)`$ and $`q_3=\mathrm{tan}(\pi (11.5\sigma ^2^2/\sigma ^4))`$ for $`d=2`$ and 3, respectively. These quantities are directly derived from the Binder cumulant $`1\sigma ^4/(3\sigma ^2^2)`$ and exhibit similar features. The $`q_d(T)`$ have small curvature in the vicinity of the transition what alleviates a precise interpolation.

Fig. 4 shows $`q_d(T)`$. For $`d=2`$, there clearly is a common intersection point of the curves for different $`L`$ at the critical temperature $`T_{\mathrm{c},2}`$. But for $`d=3`$, only a tendency towards such a behavior is seen. However, the systematic corrections to scaling can to a large extent be taken into account in a simple way: We define a size-dependent critical temperature by $`q_d(T_{\mathrm{c},d}(L),L)=q_{d,0}`$ where the $`q_{d,0}`$ are appropriate constants, fixed below. Scaling of $`q_d(T)`$ curves for different $`L`$ with respect to $`TT_{\mathrm{c},d}(L)`$ yields very good data collapse. (Referring to $`t=0`$ instead of $`T=T_{\mathrm{c},d}(\mathrm{})`$ considerably reduces the influence of deviations from scaling on the values for the critical exponents Hase .) Thus we assume $`q_d`$ to depend only on $`t=a_d(L)(TT_{\mathrm{c},d}(L))`$ in the vicinity of the transition.

We approximate $`q_d(t)`$ by polynomials of third degree, $`q_{d,0}+t+b_dt^2+c_dt^3`$. By regression studies we first adjusted the $`L`$ independent parameters $`b_d`$ and $`c_d`$ of the ansatz, and then determined the $`a_d(L)`$ values. The latter were analyzed by means of power law fits including various $`L`$ intervals, where according to finite-size scaling $`a_d(L)L^{1/\nu }`$ was presumed. For consistency, the mean-square deviation of these fits must be understandable as resulting from random errors alone. Tab. I presents the most precise results for $`\nu `$, which were obtained from the fits safely fulfilling this requirement.

In obtaining $`T_{\mathrm{c},d}(L)`$ from $`q_d(T_{\mathrm{c},d}(L),L)=q_{d,0}`$, a deviation $`\delta `$ of $`q_{d,0}`$ from the $`L\mathrm{}`$ limit of $`q_d(T,2L)=q_d(T,L)`$ gives rise to a contribution $`\delta /a_d(L)`$ to $`T_{\mathrm{c},d}(L)`$. We chose $`q_{d,0}`$ so that this term vanishes: $`q_{2,0}=1.8933`$ and $`q_{3,0}=0.174`$. The remaining higher order corrections in $`T_{\mathrm{c},d}(L)`$ originate from imperfection of finite-size scaling. Comparing several empirical approximations, we observed that, over a wide $`L`$ range, they are almost proportional to $`a_d(L)^2`$, see Fig. 5. Corresponding extrapolations yield the following values of $`T_{\mathrm{c},d}(\mathrm{})`$: $`0.10308\pm 0.00002`$ and $`0.12881\pm 0.00003`$ for $`d=2`$ and 3, respectively. The confidence intervals include the $`3\sigma `$-random errors and cautious estimates for the systematic uncertainty of the extrapolation $`L\mathrm{}`$, see Fig. 5.

The analysis of $`c(T,L)`$, $`|\sigma |(T,L)`$, and $`\chi (T,L)`$ was performed similarly to the evaluation of $`q_d(T,L)`$: We considered $`\mathrm{ln}c`$, $`\mathrm{ln}|\sigma |`$, and $`\mathrm{ln}\chi `$ as functions of $`t`$ and $`L`$. For not too large $`|t|`$, as $`L\mathrm{}`$, scaling implies that each of these quantities is decomposable into a sum of two functions depending only on $`t`$ and $`L`$, respectively. However, for the $`L`$ regions considered here, this hypothesis proved to be well fulfilled only for $`\mathrm{ln}\chi `$. In the cases of $`\mathrm{ln}c`$ and $`\mathrm{ln}\sigma `$, there is a clear tendency towards such a behavior, but small deviations cannot be neglected. Thus we approximated $`\mathrm{ln}c`$, $`\mathrm{ln}|\sigma |`$, and $`\mathrm{ln}\chi `$ by polynomials in $`t`$ of third order, taking advantage of universalities in the coefficients as far as possible. This regression provides precise values for the observables at $`t=0`$. Simultaneously, we obtained the confidence intervals taking into account the uncertainties in the individual measurements of the observables and in the $`T_{\mathrm{c},d}(L)`$ values.

The interpolation results for $`|\sigma |`$ and $`\chi `$ were analyzed by means of power law fits, where proportionality to $`L^{\beta /\nu }`$ and $`L^{\gamma /\nu }`$, respectively, was presumed. These studies were performed analogously to the determination of $`\nu `$. Tab. I presents the values of the exponent ratios $`\beta /\nu `$ and $`\gamma /\nu `$. The high quality of these power law fits is illustrated by Fig. 6 presenting $`|\sigma |(T_\mathrm{c}(L),L)`$.

Compared to the study of $`|\sigma |`$ and $`\chi `$, the analysis of $`c`$ is more difficult: Exponent values obtained from power law fits converge only slowly with increasing $`L`$, and the mean-square deviations remain too large. Therefore we took into account a background contribution presuming $`c(T_\mathrm{c}(L),L)=f(L)=a+bL^{\alpha /\nu }`$. However, we treated $`f(L_1)`$ and $`f(L_2)`$ at fixed finite $`L_1`$ and $`L_2`$ as adjustable parameters instead of $`a`$ and $`b`$ to avoid numerical problems for almost logarithmic behavior of $`c(L)`$ (small $`|\alpha |`$). Results for $`\alpha /\nu `$ obtained this way are included in Tab. I.

Our values in Tab. I have to be regarded as effective exponents. Due to the finiteness of $`L`$, tiny systematic errors are certainly present, presumably the more relevant the smaller the exponent value. Unfortunately, our data set is not sufficient for a convincing estimate of them. However, even if only random errors are considered, our values comply with the Widom relation, $`2=\alpha +2\beta +\gamma `$, and the hyperscaling relation, $`2\alpha =d\nu `$.

For comparison, Tab. I includes also values for the Ising model with nearest-neighbor interaction G92 ; Hase . It is surprising that, in spite of the differing characters of the interactions, the critical exponent data obtained here are very close to these values: In particular, for $`\gamma /\nu `$ and $`\nu `$. the agreement is perfect within numerical accuracy. The slight deviations in $`\alpha /\nu `$ and $`\beta /\nu `$ for $`d=3`$ presumably arise from a too small sample size. Note: obtaining these values indirectly, via Widom and hyperscaling relation from $`\gamma /\nu `$ and $`\nu `$, yields again perfect agreement.

Concluding, in spite of the long-range interaction, the Coulomb system described by Eq. (1) seems to belong to the same universality class as the Ising model with short-range interaction. This suggests that the lattice Coulomb-glass model might have the same critical properties as the random-field short-range Ising model.

###### Acknowledgements.

We thank H. Eschrig, M.E. Fisher, B. Kramer, T. Nattermann, M. Richter, M. Schreiber, and T. Vojta for helpful discussions and literature hints.


format: text

text:
# Lévy Model for Interstellar Scintillations

## Abstract

Observations of radio signals from distant pulsars provide a valuable tool for investigation of interstellar turbulence. The time-shapes of the signals are the result of pulse broadening by the fluctuating electron density in the interstellar medium. While the scaling of the shapes with the signal frequency is well understood, the observed anomalous scaling with respect to the pulsar distance has remained a puzzle for more than 30 years. We propose a new model for interstellar electron density fluctuations, which explains the observed scaling relations. We suggest that these fluctuations obey Lévy statistics rather than Gaussian statistics, as assumed in previous treatments of interstellar scintillations.

PACS numbers: 95.30.Qd, 98.38.Am, 98.38.Dq, 95.85.-e

1. Introduction. Electron density fluctuations in the interstellar medium (ISM) cause scintillations of the intensity of signals arriving from distant pulsars. If the medium were completely transparent, the shape of the arriving signal would coincide with the shape of the signal emitted by the pulsar. However, the observed pulse is much broader, and this effect is attributed to the random refraction the waves experience while they travel through the medium . To investigate pulse broadening one can assume that the pulsar intrinsic signal is narrow in time, $`I_0(t)\delta (tt_0)`$, where $`I_0(t)`$ is the signal intensity. The observed signal is broad and asymmetric, with a sharp rise and a slow decay; see Fig. 1. Observations show that broadened shapes of the pulses are similar for different pulsars (after proper rescaling), suggesting that the density fluctuations statistics along different lines of sight are to some extent universal.

For estimates assume that the pulsar distance is $`d10\text{kpc}`$, the typical electron density is $`n0.03\text{cm}^3`$, and the observational wave frequency is $`\nu 500\text{MHz}`$. Then the plasma electron frequency $`\omega _{pe}=(4\pi ne^2/m_e)^{1/2}`$ is much smaller than $`\nu `$, and density fluctuations change the wave phase only slightly. To estimate the time delay one can use the approach of geometric optics, where the propagating ray is refracted (scattered) by small prisms of density inhomogeneities . At each scatter event occurring at a mean-free path $`l`$, the propagation angle changes by a small amount, $`\mathrm{\Delta }\theta \lambda ^2r_0\mathrm{\Delta }n`$ \[see below\], where $`\lambda `$ is the wavelength, $`r_0=e^2/m_ec^2`$ is the classical radius of the electron, and $`\mathrm{\Delta }n`$ is the density difference at characteristic separation $`l`$. Using the standard assumption that $`\mathrm{\Delta }\theta `$ is random and Gaussian, one finds that the path direction deviates from a straight line by $`\theta \lambda ^2r_0\mathrm{\Delta }n_0(d/l)^{1/2}`$, where $`\mathrm{\Delta }n_0`$ is the characteristic amplitude of density-difference fluctuations, and the path length deviates from the distance $`d`$ by $`\mathrm{\Delta }dd\theta ^2\lambda ^4d^2`$. The broadening time can be estimated as $`\tau _d\mathrm{\Delta }d/c`$, which gives the standard scaling $`\tau _d\lambda ^4d^2`$.

Observations show that the signal width, $`\tau _d`$, estimated at the half-amplitude level, scales with the wave length according to the obtained formula, $`\tau _d\lambda ^4`$, while the scaling with distance is close to $`\tau _dd^4`$, contradicting the analytical prediction, as is seen in Fig. 2 . This paradox was first discussed by Sutton , and although the theory of scintillations has been developed for more than 30 years, the contradiction has resisted analytical understanding .

In this paper we propose that the anomalous scaling with the distance is an evidence of non-Gaussian density fluctuations in the ISM. We suggest that the probability distribution of density gradients has a power-law decay, and its second moment is divergent. Such probability distributions are common in theories of turbulence, as is consistent with the argument that the density statistics are governed by turbulent motions in the ISM . The sum of many angular deviations caused by such fluctuations does not have a Gaussian distribution; instead, the limiting distribution is of the Lévy type, and the ray angle performs a Lévy flight instead of a conventional random walk. We present a solvable model of scintillations that allows us to unify and extend to a non-Gaussian case the standard analytical approaches, see, e.g., . We then apply this model to Lévy density statistics, compare it to the observational data, and demonstrate that the model naturally produces correct scalings of the signals. We report main results here, the detailed discussion is presented in .

2.Wave equation in a random medium. The Fourier amplitude of electric (or magnetic) field, $`E_\omega (𝐫)`$, in the isotropic ISM with dielectric permittivity $`ϵ_\omega `$ obeys the wave equation:

$`\left[\mathrm{\Delta }{\displaystyle \frac{\omega ^2}{c^2}}ϵ_\omega (𝐫)\right]E_\omega (𝐫)=0,`$ (1)
where $`ϵ_\omega (𝐫)=1\omega _{pe}^2(𝐫)/\omega ^2`$, and the electron plasma frequency $`\omega _{pe}(𝐫)`$ changes slowly on the wave scale $`\lambda `$. Assuming that the wave propagates in the line-of-sight direction, $`z`$, we separate the quickly-changing phase of the wave from the slowly-changing amplitude, $`E_\omega (𝐫)=\mathrm{exp}(iz\omega /c)\mathrm{\Phi }_\omega (z,𝐱)`$, where $`𝐱`$ is a coordinate perpendicular to $`z`$. Substituting this into the wave equation (1), we derive the equation for the wave amplitude,

$`\left[2i{\displaystyle \frac{\omega }{c}}{\displaystyle \frac{}{z}}+\mathrm{\Delta }_{}4\pi r_0n(𝐱,z)\right]\mathrm{\Phi }_\omega (𝐱,z)=0,`$ (2)
where $`\mathrm{\Delta }_{}`$ is a two-dimensional Laplacian in the $`𝐱`$ plane. Following , introduce the function $`I(𝐫_1,𝐫_2,t)=\mathrm{\Phi }(𝐫_1,t)\mathrm{\Phi }^{}(𝐫_2,t)`$, whose Fourier transform with respect to time is $`I_\mathrm{\Omega }(𝐫_1,𝐫_2)=\sqrt{2/\pi }\text{d}\omega \mathrm{\Phi }_{\omega +\mathrm{\Omega }/2}(𝐫_1)\mathrm{\Phi }_{\omega \mathrm{\Omega }/2}^{}(𝐫_2)`$. For coinciding coordinates this function is the intensity of the radiation whose variation in time we seek. To find this function, we may first solve the equation for $`V_{\omega ,\mathrm{\Omega }}(𝐫_1,𝐫_2)\mathrm{\Phi }_{\omega +\mathrm{\Omega }/2}(𝐫_1)\mathrm{\Phi }_{\omega \mathrm{\Omega }/2}^{}(𝐫_2)`$, which can be derived from Eq. (2). Assuming that $`\mathrm{\Omega }\omega `$, we obtain

$`i{\displaystyle \frac{V}{z}}={\displaystyle \frac{2k+\mathrm{\Delta }k}{4k^2}}{\displaystyle \frac{^2V}{𝐱_{\mathrm{𝟐}}^{}{}_{}{}^{2}}}{\displaystyle \frac{2k\mathrm{\Delta }k}{4k^2}}{\displaystyle \frac{^2V}{𝐱_{\mathrm{𝟏}}^{}{}_{}{}^{2}}}+{\displaystyle \frac{2\pi r_0}{k}}\mathrm{\Delta }nV,`$ (3)
where we denoted $`k=\omega /c`$, $`\mathrm{\Delta }k=\mathrm{\Omega }/c`$, and $`\mathrm{\Delta }n=n(𝐱_1,z)n(𝐱_2,z)`$.

Equation (3) is hard to solve without further simplification since $`n(𝐱,z)`$ is an unknown random function. The standard procedure is to assume that the density fluctuations are Gaussian with a specified correlator in $`x`$ and only short-scale correlations in $`z`$, see . Eq. (3) can then be averaged over the Gaussian ensemble of density fluctuations and over different positions in $`𝐱`$-space. However, the resulting solution yields a scaling of $`\tau _d\lambda ^4d^2`$ that contradicts observations, as noted above.

We propose that the turbulent gas motions in the ISM give rise to strongly intermittent and non-Gaussian density fluctuations. If the distribution function of $`\mathrm{\Delta }n`$ has a power-law decay as $`|\mathrm{\Delta }n|\mathrm{}`$ and has no second moment, then the sum of many independent ray angle deviations does not behave as a Gaussian variable (the central limit theorem does not hold). Instead, the limiting distribution, if it exists, is the Lévy distribution. A random walk whose increments are Lévy distributed is called a Lévy flight. Such processes are common in various random systems, and often replace Brownian motion in turbulent systems .

The Fourier transform (the characteristic function) of a symmetric Lévy distribution $`P_\beta (\mathrm{\Delta }n)`$ has the simple form,

$`F(\mu )={\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}\mathrm{\Delta }nP_\beta (\mathrm{\Delta }n)\mathrm{exp}(i\mu \mathrm{\Delta }n)=\mathrm{exp}\left(C|\mu |^\beta \right),`$ (4)
where $`0<\beta <2`$, and $`C`$ is some positive constant. Eq. (4) can be taken as the definition of a symmetric Lévy distribution. The sum of $`N`$ Lévy distributed variables scales as $`^N\mathrm{\Delta }nN^{1/\beta }`$, which becomes diffusion in the Gaussian limit $`\beta =2`$. For $`\beta <2`$, the probability distribution function has algebraic tails, $`P_\beta (\mathrm{\Delta }n)|\mathrm{\Delta }n|^{1\beta }`$ for $`|\mathrm{\Delta }n|\mathrm{}`$, and its second moment is divergent. We thus assume that the random density-gradient fluctuations are Lévy distributed, and are short-scale correlated in $`z`$. In the next section we first show how Eq. (3) can be solved for a general case of non-Gaussian random density field. In Sec. 4 we apply our method to Lévy distribution.

3. Batchelor approximation. Propagation as described by Eq. (3) cannot be simplified in general, when $`\mathrm{\Delta }n`$ is not a short-scale correlated Gaussian random variable. However, analytical investigation is possible in the important case of smooth turbulent fluctuations. This case is analogous to the Batchelor limit, in the problem of turbulent random advection . For that approximation in this paper, we neglect all effects other than those of density gradients: $`n(𝐱_1)n(𝐱_2)\sigma \sigma \sigma (z)(𝐱_1𝐱_2)`$, where the density gradient $`\sigma \sigma \sigma (z)`$ is a random variable with correlation length $`ld`$ along $`z`$. In this approximation, the variables separate in Equation (3) and it can be solved exactly. We leave analysis of more complicated cases for further communcation, and present here results for this simple case, which captures the essential physics.

As a further simplification, consider one-dimensional variables $`𝐱_1`$ and $`𝐱_2`$. Since the variables separate, we can look for the solution in the factorized form $`V(x_1,x_2,z)=U_1(x_1,z)U_2(x_2,z)`$. Then the equation for $`U_1`$ reads:

$`i{\displaystyle \frac{U_1}{z}}={\displaystyle \frac{2k\mathrm{\Delta }k}{4k^2}}{\displaystyle \frac{^2U_1}{x_1^2}}+{\displaystyle \frac{2\pi r_0}{k}}\sigma (z)x_1U_1.`$ (5)
The analogous equation for $`U_2`$ is obtained by changing $`kk`$. The solution of equation (5) is sought in the form $`U_1(x_1,z)=A(z)\mathrm{exp}\left[iB(z)x_1+iC(z)x_1^2\right]`$, with the initial condition $`U_1(z=0)=\delta (x_1)`$, if the refracting medium extends all the way up to the pulsar. Substituting this ansatz into (5), we find

$`A(z)`$ $`=`$ $`{\displaystyle \frac{A_0}{\sqrt{z}}}\mathrm{exp}\left({\displaystyle \frac{i(2k\mathrm{\Delta }k)}{4k^2}}{\displaystyle \underset{0}{\overset{z}{}}}B^2(z^{})\text{d}z^{}\right),`$ (6)
$`B(z)`$ $`=`$ $`{\displaystyle \frac{2\pi r_0}{kz}}{\displaystyle \underset{0}{\overset{z}{}}}\sigma (z^{})z^{}\text{d}z^{},C(z)={\displaystyle \frac{k^2}{(2k\mathrm{\Delta }k)z}}.`$ (7)
Note that this solution describes the path of a single ray through a sequence of density gradients $`\sigma (z)`$. Effects of multiple rays can be found from superposition. The intensity of received radiation can be calculated from the Fourier transform $`I_\omega (z,t)=_{\mathrm{}}^{\mathrm{}}V_{\omega ,\mathrm{\Omega }}(x=0,z)\mathrm{exp}(i\mathrm{\Omega }t)\text{d}\mathrm{\Omega }/\sqrt{2\pi }.`$ In this Fourier transform, individual ray paths will yield contributions with phase proportional to $`\mathrm{\Omega }=\mathrm{\Delta }kc`$, with coefficient equal to the travel time for that path. Cross-terms describing interference of paths yield contributions that oscillate rapidly with frequency and average to zero, see e.g. . The intensity, averaged over an ensemble of statistically independent rays, is then given by the average over individual travel times or equivalently over different realizations of $`\sigma (z)`$.

This leads to

$`I_\omega (z,t)\delta \left(t{\displaystyle \frac{d}{c}}{\displaystyle \frac{1}{2k^2c}}{\displaystyle \underset{0}{\overset{z}{}}}B^2(z^{})\text{d}z^{}\right),`$ (8)
where the angular brackets denote the statistical average. Note that $`B(z)`$ is proportional to the deflection angle $`\theta `$ of the ray. Formula (8) gives the shape of the signal observed at the Earth; if the scattering medium were absent, the signal would be undistorted, $`I(t)\delta (td/c)`$.

Lee & Jokipii investigated Eq. (3) for short-scale correlated Gaussian density fluctuations . For averaging over Gaussian $`\sigma (z)`$, our solution (8) reproduces those obtained by Williamson  with a phenomenological approach. Thus, Williamson’s solutions are applicable under the assumption of smooth Gaussian density fluctuations, and when one keeps only the linear term in the expansion of $`\mathrm{\Delta }n`$. In the next section we apply our approach to the Lévy distributed density fluctuations.

4. Scintillations as Lévy flights through the interstellar medium. The averaging in formula (8) can be performed for the Lévy distributed short-scale correlated density gradients $`\sigma (z)`$. To do this, we represent integrals is (7), and (8) in the discretized forms, i.e., we assume that $`d=nl`$, $`z^{}=ml`$, and $`z^{\prime \prime }=sl`$, where $`l`$ is the correlation length of density fluctuations, and change $`_0^zf(z^{})\text{d}z^{}_{m=1}^nf(lm)l`$ for an arbitrary function $`f(z)`$. The right hand side of (8) is the probability distribution of the propagation time delay, $`\tau =td/c`$. For a continuous medium, this time delay is given by

$`\tau ={\displaystyle \frac{r_0^2l^3\lambda ^4}{8\pi ^2c}}{\displaystyle \underset{m=1}{\overset{n}{}}}\left[{\displaystyle \frac{1}{m}}{\displaystyle \underset{s=1}{\overset{m}{}}}s\sigma _s\right]^2.`$ (9)
The solution (8) can now be calculated numerically as the probability distributions of this variable, under the assumption that $`\sigma _s`$ are distributed independently, identically, and according to the Lévy law (4). We however need to specify the parameter $`\beta `$ in the Lévy formula. We do this by comparing our model with observations. The observed scaling of pulse broadening is close to $`\tau _d\lambda ^4d^4`$, while our model gives $`\tau \lambda ^4d^{(2+\beta )/\beta }`$, as is seen from the scaling for sums of Lévy-distributed variables, $`^m\sigma _sm^{1/\beta }`$, following Eq. (4). Thus, we obtain, $`\beta 2/3`$.

Note that standard Gaussian models of density distribution were not able to satisfy both observational scalings, $`\lambda ^4`$ and $`d^4`$, simultaneously. Various studies of non-smooth density fluctuations within these models have not reproduced this scaling either . Our model reproduces the anomalous $`d`$-scaling naturally. Moreover, it predicts that the probability distribution function of electron density gradients in the interstellar turbulence decays as $`P(\sigma )|\sigma |^{1\beta }|\sigma |^{5/3}`$. Power-law distributions $`P(\sigma )`$ with $`\beta <2`$ are indeed observed in numerical simulations of compressible turbulence , however, no one has yet derived them from first principles. To date theories of scintillations have exploited only second-order correlators of the density fluctuations, while in our approach these correlators do not exist (or do not matter) and one must work with the whole probability distribution function.

Although our goal was to explain the scalings of the signals, it is interesting to see to what extent we can predict their shapes. The delay time is proportional to the square of the typical deflection angle of the ray trajectory, $`\tau \theta ^2`$, where $`\theta `$ has the Lévy distribution $`P_\beta (\theta )`$. Therefore, the distribution of arrival times is $`I(\tau )P_\beta (\tau ^{1/2})\tau ^{1/2}`$, with the asymptotic form $`I(\tau )\tau ^{1\beta /2}`$ as $`\tau \mathrm{}`$. Fig. (3) shows the distribution of $`\tau `$ from numerical calculation, with a power-law decay at long times, as expected. Because the observed shape of the scattered pulse is directly related to the probability distribution of gradients in electron density in the ISM, observational data offer the possibility of characterizing interstellar turbulence .

The curve in Fig. (3) closely resembles the observed signals, although the presented analytical shapes are the result of averaging over an ensemble of non-interfering rays, corresponding to an observational average over an infinite amount of time. In practice, the averaging time is finite, and the long tail of the distribution, dominated by rare events, may not have converged. We also ignore instrumental response . Moreover, a non-analytic density field is more natural for a turbulent cascade . Also, the small-scale density fluctuations that produce scintillations should be collisionless, and elongated along the local magnetic field ; the scattering is likely to be highly anisotropic, and locally nearly 1-dimensional. We will consider these effects in future work. Interestingly, some scattered pulsars show power-law declines at long times, such as that in Fig. 1. Scintillation of nearby pulsars also shows evidence for weak large-angle scattering . Some interferometric studies suggest a “halo” surrounding the source at large scattering angles and excess scattering at small angles relative to a Gaussian , as might be expected for a Lévy distribution of scattering angles. Intrinsic source structure, and the relatively short observatonal averages, may complicate this interpretation.

Finally, we wish to comment on the original explanation of the anomalous $`d`$-scaling by Sutton . Sutton suggested that encounters with much more-strongly-scattering HII regions become more probable on longer lines of sight. This however requires a perhaps surprisingly close coordination of DM (over 1.7 orders of magnitude) with $`\tau _d`$ (over 8 orders of magnitude). Sutton’s proposal assumes essentially non-stationary statistics for the density distribution along $`z`$. Our proposal also invokes rare, large events, but in a statistically stationary way.

To summarize, we propose that the observed anomalously strong time-broadening of pulsar signals is evidence for non-Gaussian distribution of electron density gradients in the ISM. We argue that this distribution is of the Lévy type, in accord with the turbulent origin of density fluctuations, and we present a simple model that explains the observational scalings of pulsar signals.

The work of SB was supported by the ASCI Flash Center at the University of Chicago, under DoE subcontract B523820.


text:
# The Surroundings of Gamma-Ray Bursts: Constraints on Progenitors

## 1 Introduction

Some of the best evidence for nature of gamma-ray burst (GRB) progenitors has come from the identification of the Type Ic supernovae SN 1998bw with GRB 980425 Ge98  and of the recent SN 2003dh with GRB 030329 Se03 . The finding of these events supports models of long-duration GRBs originating from stripped massive stars MWH01 . The surroundings of massive stars are expected to be shaped by the winds emanating from the progenitor stars. Clear evidence for the signature of a wind has been difficult to establish, and the possibility remains that there is more than one type of progenitor for the long-duration bursts. Direct interaction with the interstellar medium might be expected if the progenitor involves a compact binary system.

Another possibility is that the GRB occurs in a massive star weeks to years after it has become a supernova VS98 . The progenitor of the GRB may be a rapidly spinning neutron star that spins down and eventually collapses, leading to a burst. This can create a pulsar wind nebula immediately surrounding the burst progenitor, which has possible advantages for producing a GRB KG02 .

## 2 Afterglows

### 2.1 Afterglows and the circumburst medium

The afterglows of GRBs provide a probe of the immediate surroundings of GRBs. The evidence that we have from bursts related to supernovae is that the progenitors are Type Ic supernovae. In addition to SN 1998bw and SN 2003dh, the burst GRB 021211 has possibly been identified with a Type Ic supernova like SN 1994I DVe03 . These supernovae are thought to have massive star progenitors that have been stripped of their hydrogen envelopes, i.e. Wolf-Rayet stars. This type of massive star progenitor is also suggested by the argument that the collimated flow from a burst be able to pass through the star in a time that is not significantly longer that the duration of the GRB. Even for the long duration bursts, this implies a relatively compact stellar progenitor like the Wolf-Rayet stars MWH01 ; M03 .

Wolf-Rayet stars are found to have winds with typical mass loss rate $`\dot{M}=10^5M_{}\mathrm{yr}^1`$ and wind velocity $`v_w=1000\mathrm{km}\mathrm{s}^1`$ ending in a termination shock where the wind runs into the surrounding medium CL99 . A steady wind produces a density distribution $`\rho =Ar^2`$; the value of the density can be scaled to the corresponding value for the standard wind parameters, $`A_{}=A/(5\times 10^{11}\mathrm{gm}\mathrm{cm}^1)`$. For a surrounding medium of pressure $`p`$, the wind termination shock occurs at a radius

$$R_t=5.7\times 10^{19}\left(\frac{v_w}{1000\mathrm{km}\mathrm{s}^1}\right)\left(\frac{p/k}{10^4\mathrm{cm}^3\mathrm{K}}\right)^{1/2}A_{}^{1/2}\mathrm{cm}$$
(1)
where $`k`$ is Boltzmann’s constant. The pressure is normalized to a typical value in the interstellar medium of our Galaxy. A higher pressure can occur as a result of the wind bubble evolution, or as the result of an especially high pressure interstellar medium.

For a spherical explosion with energy $`E`$, the blast wave in the free wind reaches the radius $`R_t`$ after a time

$$t_t=1.35\times 10^4\left(\frac{v_w}{1000\mathrm{km}\mathrm{s}^1}\right)^2\left(\frac{p/k}{10^4\mathrm{cm}^3\mathrm{K}}\right)^1\left(\frac{E}{10^{53}\mathrm{ergs}}\right)^1A_{}^2\mathrm{days}.$$
(2)
This shows that over typical times of observation, the blast wave may be expanding into the free wind for standard parameters. In view of this, interaction with a wind has become one of the models that is investigated in modeling the afterglows of GRBs CL00 . This model is compared to results for a constant density medium, which was initially taken in afterglow modeling as the simplest assumption. Taking the surrounding density of the form $`\rho r^s`$, these cases can be designated $`s=0`$ (uniform) and $`s=2`$ (wind). In afterglow models, both wind and constant density models provide adequate fits to the data in some cases, although the constant density model is usually favored PK01 ; PK02 .

At first sight, this ambiguity is surprising. In the context of models with constant efficiencies evolving before a jet break sets in, the two cases have distinctive behavior: the synchrotron self-absorption frequency, $`\nu _a`$, evolves to higher frequency with $`s=2`$, but remains constant for $`s=0`$; the peak flux drops with time for $`s=2`$, but remains constant for $`s=0`$; the synchrotron cooling frequency, $`\nu _c`$, evolves to higher frequency for $`s=2`$, but evolves to lower frequency for $`s=0`$. There are various reasons why these differences have not provided clear tests of the models. The evolution of $`\nu _c`$ requires good light curve information at optical/IR and X-ray wavelengths, which is usually lacking. The evolution of $`\nu _a`$ requires early radio data; these data are usually sparse and are affected by interstellar scintillation. Another problem is that jet breaks are observed in light curves and the evolution in the post-break regime can mimic some of the features of evolution in an $`s=2`$ medium: the peak flux drops with time and $`\nu _a`$ evolves to higher frequency SPH99 .

Another possibility for distinguishing between the models is to go to very early times, within about a minute of the GRB burst. In addition to the older case of GRB 990123, this has recently been achieved for GRB 021004 Fox03  and GRB 021211 Le03 ; We02 ; Pe02 . The advantage of these early times is that jet effects do not play a role in the evolution. Li & Chevalier LC03  suggested that the early flat optical light curve of GRB 021004 could be interpreted in terms of wind interaction, in which the critical frequency $`\nu _m`$ had not yet moved down through optical wavelengths. Although this model has some promise, the case remains ambiguous in that the relatively flat evolution might also be produced by a combination of emission from the reverse shock wave and the later forward shock emission in the case of interaction with a $`s=0`$ medium KZ03 . To distinguish between these possibilities, both good light curve data and color information are needed.

An additional problem with GRB 021004 is that the optical light curve showed variability superposed on the overall trend He03 . This makes it difficult to clearly specify characteristic times in the evolution of the afterglow. The variability, which is not seen in all afterglows, may be due to density inhomogeneities in the circumburst medium and may thus provide additional diagnostics for that medium HP03 . The winds of Wolf-Rayet stars are known to be inhomogeneous, with clumps filling $`1/4`$ of the volume HK98 . However, the observed degree of inhomogeneity refers to a region close to the stellar surface and the inhomogeneity may decrease in the outer parts of the wind where the afterglow occurs.

### 2.2 Jets

Jet breaks themselves provide a possible diagnostic for the medium. The usual assumption has been that an afterglow light curve steepens due to geometric and spreading effects when its Lorentz factor is about 1/(opening angle) to the form $`t^p`$ SPH99 . Kumar & Panaitescu KP00  find that the transition to the asymptotic jet evolution requires a factor $`>10`$ in time for expansion in a uniform medium, but expansion by $`10^4`$ in time for expansion in a wind medium. However, they used a simplified treatment of jet evolution. Numerical simulations indicate that sharp jet breaks do occur in a uniform medium, but that most of the emission remains within the initial opening angle of the jet Ge01 . Granot & Kumar GK03  have recently considered structured jets in uniform and declining density media, and again found that jets in an $`s=2`$ medium cannot give sharp jet breaks.

The recent burst GRB 030329, which was clearly associated with a supernova, showed a sharp break in the light curve at $`t0.5`$ day, which has been interpreted as a jet break Be03 . However, it has become clear that there is considerable structure in the light curve of GRB 030329 over the first 10 days Ue03 ; Ma03  and the identification of an early jet break cannot be made with certainty.

### 2.3 Afterglow parameters

There are uncertainties in the basic assumptions involved in standard afterglow modeling, which include constant values of electron energy efficiency, $`ϵ_e`$, magnetic energy efficiency, $`ϵ_B`$, and particle spectral index, $`p`$, in the evolution of one source (e.g., SPN98 ). If these parameters remain constant during the evolution of one source, the expectation would be that they tend toward “universal” values that apply to various sources. Standard models developed for observed afterglows do not show this. As an example, I take the results of Panaitescu & Kumar PK02 , who treat a set of the 10 best observed afterglows with the standard assumptions. The values of spectral index, $`p`$, cover the range $`1.362.78`$. The presence of values $`p<2`$ is noteworthy because most of the particle energy is at high energy for this case, although the number of particles is dominated by the low energy particles. The values of $`ϵ_B`$ and $`ϵ_e`$ cover the ranges of $`4\times 10^50.07`$ and $`0.010.4`$, respectively.

The theoretical values of these parameters are poorly known. The production of the magnetic field requires some mechanism at the forward shock front to build up the magnetic field. The mechanism remains uncertain, although the Weibel instability has been suggested ML99 , and recent simulations show some promise for this mechanism Ne03 ; Fe03 . Studies of Fermi-type particle acceleration in ultarelativistic shocks have yielded a preferred value of $`p=2.22.3`$ in the test particle limit Ae01 . The way electrons are injected into the acceleration process remains uncertain, and there are indications that the acceleration may differ from the Fermi process Ne03 ; Fe03 .

The fact that a range of parameters is needed to explain the various afterglows suggests that the parameters depend on the physical conditions. It might be expected that the shock velocity and preshock density are important determinants of the physical conditions, so that the model parameters should vary during the evolution of a burst. Yost et al. Ye03  have recently considered models in which $`ϵ_B`$ is allowed to vary as a power law of the shock Lorentz factor and in which the value of the preshock density parameter $`s`$ is allowed to cover a large range. Their modeling of 4 sources shows that a wide variety of models are potentially possible, including ones in which the density increases steeply with radius. It appears that more extensive observations, including spectral information over a wide time range, are needed to further constrain the models. One possibility is to follow the evolution of the characteristic frequency, $`\nu _m`$, from optical to radio wavelengths. This frequency is typically observed at radio wavelengths on a timescale $`10`$ days. Its passage at optical wavelengths requires very early observations; as noted above, this may have been observed in GRB 021004. The typical frequency is not sensitive to the density, but it is sensitive to the efficiency factors, so that constraints on their evolution may be obtainable.

### 2.4 A shocked wind environment

In view of the evidence for the association of GRBs with massive stars and the evidence from afterglows for interaction with a constant density medium, consideration must be given as to how a massive star might produce a uniform density surrounding. The most plausible way to do this is the approximately constant density region expected downstream from the wind termination shock W01 . Wijers W01  suggested two ways for creating a smaller value of $`t_t`$ (eq. 2): reducing the mass loss rate from the progenitor to $`10^6M_{}\mathrm{yr}^1`$ because of the low metallicity of the progenitor star and increasing the pressure by interacting with dense molecular gas, especially if the progenitor star is moving. However, the metallicity dependence of the mass loss from Wolf-Rayet stars is uncertain: WN type stars in the lower metallicity Large Magellanic Cloud HK00  and the Small Magellanic Cloud C00  have similar mass loss rates to those in the Galaxy, although WC stars do seem to show a metallicity effect Ce02 . Also, photoionizing radiation during the life of a massive star tends to clear a region around the star to a moderately low density. Another way of increasing $`p`$ is by having the burst occur in a high pressure starburst region F03 , where the pressure can reach values of $`p/k>10^8`$ cm-3 K CF01 . In this case, there should be a relation between the properties of the afterglow (relatively dense surroundings) and the position of the burst relative to a region of very active star formation.

One expectation of the models with a termination shock is that some bursts should be observed to make a transition from an $`s=2`$ to an $`s=0`$ medium, with a density jump between them. There has been little evidence for such a transition. The expectation for such a transition is that the light curve should evolve to a flatter asymptotic decline after a jump in flux. Wijers W01  mentioned GRB 970508 because it had a bump in the optical light curve at an age of 1 day. However, it did not show the expected flattening of the light curve. A burst that showed a steepening with a possible bump is GRB 030226. Dai & Wu DW03  suggested that the transition was due to the interaction with a large density jump, which might occur at the contact discontinuity between the shocked progenitor wind and the dense red supergiant wind from a previous evolutionary phase. In this picture, the steepening of the light curve is due to the sideways expansion of the jet in the dense medium. However, the data on GRB 030226 are not of sufficient quality to clearly show the expected features at the time of transition with the density jump. In the case of interaction with the termination shock of the stellar wind, the density jump is such that the effects of the reverse shock are not expected to be important, as opposed to the high density jump case.

## 3 Optical/ultraviolet absorption lines

A recent development relevant to the surroundings of GRBs is the observation of strong optical/ultraviolet absorption lines in a few cases. The best case is GRB 021004, which has a redshift $`z=2.32`$ so that strong ultraviolet lines in the rest frame are redshifted to optical wavelengths Me03 ; Sce03 . Strong lines of Ly$`\alpha `$, Ly$`\beta `$, C IV, and Si IV are found blueshifted relative to the host velocity by $`450`$, $`990`$, and $`3155\mathrm{km}\mathrm{s}^1`$ Me03 . The lines have not been observed to vary, so they cannot be directly tied to the immediate circumburst environment, but both Mirabal et al. Me03  and Schaefer et al. Sce03  argue that they are likely to be formed in the nearby environment. One argument for this is that the strong lines are unusual for intervening systems observed in the spectra of quasars. The lines may be formed close to the host galaxy if the burst occurred in a starburst region with a strong galactic superwind. However, the maximum velocity shift is higher than has been observed in galactic superwinds. In addition, the high velocity would require a high initial temperature for the gas if the wind is thermally driven. The gas would be completely ionized and it is unlikely that it would be able to cool to allow the observed ions.

For a circumstellar origin, there are two possibilities: the high velocities are related to the wind velocities of the progenitor system, or the velocities are due to radiative acceleration by the GRB light; combinations of these models are also possible. Schaefer et al. Sce03  argue that the high blueshifted velocity can be naturally produced by the Wolf-Rayet star wind velocity and the lower velocity components can be identified with denser shells swept-up by the fast wind. Mirabal et al. Me03  argue that the abundances deduced from the lines disfavor the wind model. Hydrogen is present in the observed lines, but it is also observed in the spectra of some WN stars SSM96 ; however, the absence of N V lines in the observed spectra indicates that N is not overabundant. In the radiatively accelerated model, the accelerated clumps must be at an initial distance of several 0.1’s pc from the progenitor star. The acceleration, primarily by bound-free transitions, must occur early in order to avoid the observation of time variability.

The problem with both of these scenarios is that the strong radiation field from the GRB and its afterglow is able to completely photoionize the gas out to a distance $`>10^{18}`$ cm Le02 . At a sufficiently high density ($`>10^7`$ cm-3), the recombination time becomes shorter than the age of the burst. This density might be present in clumps of the swept-up red supergiant wind, especially if the progenitor is in a high pressure region, but further exploration of this topic is needed.

GRB 021004 is not alone in showing these line features. High excitation, high velocity absorption features have been found in GRB 020813 Ba03  and GRB 030226 Ge03 ; Pr03 ; CF03 . The absorption lines of CIV in GRB 020813 are at $`0\mathrm{km}\mathrm{s}^1`$ and $`4320\mathrm{km}\mathrm{s}^1`$ relative to the host. In this case, the blueshifted absorption is also present in a number of lower ionization species (Si II, Al II, Fe II, Mg II, and Mg I); there is no coverage of Ly$`\alpha `$. In the case of GRB 030226, strong absorption line systems are present at a velocity separation of $`2300\mathrm{km}\mathrm{s}^1`$, with C IV and Si IV present, as well as numerous lower ionization species and Ly$`\alpha `$. The velocity separation seen in these sources is consistent with expectations for the velocity of a Wolf-Rayet star wind. However, the presence of H does not support this origin for the lines.

## 4 Pulsar wind bubble environment

An interesting possibility for a GRB environment is that created by a pulsar wind nebula. This possibility was proposed in the context of the supranova model in which the supernova precedes the GRB VS98 ; KG02 . In this scenario, the supernova core contracts to a massive, rapidly rotating neutron star which spins down and collapses to a black hole after a period of weeks, months, or years. If the neutron star has a magnetic field similar to that of radio pulsars, it can create a shocked bubble of relativistic electron/positron fluid and magnetic field before it collapses. The bubble accelerates the supernova ejecta, so that the ejecta can play a role in producing the X-ray lines that have possibly been observed in some bursts.

One issue is how effectively the pulsar nebula can accelerate the supernova gas, because the situation is subject to Rayleigh-Taylor instabilities KG02 ; IGP03 . The X-ray line features are typically blueshifted by $`0.1c`$, which is higher than the velocities of the heavy element ejecta that would be expected from the supernova itself. A supernova energy of $`10^{51}`$ ergs in an ejecta mass of $`10M_{}`$ leads to a typical velocity of $`0.01c`$. Königl & Granot KG02  suggest that the rotational energy of a rapidly rotating neutron star, $`10^{53}`$ ergs, is transferred to the ejecta, giving the observed velocity. Although a pulsar nebula can certainly shock and compress the ejecta gas, the ability to further accelerate the ejecta is less certain. However, if the supernova explosion energy is high ($`>10^{52}`$ ergs as inferred for some supernovae) and the ejecta mass is low, a typical velocity of $`0.1c`$ can be attained.

Another issue is the fact that the observation of X-ray line features at an age of $`1`$ day requires dense gas at $`r10^{16}`$ cm, but observations of afterglow emission at an age of a week or more imply a radial scale $`>3\times 10^{17}`$ cm. Königl & Granot KG02  suggest that pulsar nebula and supernova may be highly elongated along the axis along which the GRB flow propagates.

Although the observational evidence for X-ray lines remains controversial, there are other reasons for considering a pulsar wind nebula environment, as articulated by Königl & Granot KG02 . One is that the pulsar nebula is composed of just the ingredients that are necessary for the synchrotron emission from a GRB afterglow. There is no problem with the efficiencies for production of the synchrotron emission. In addition, the bubble density can be constant with radius, or drop with radius. The first case can occur in the same way that a constant density occurs downstream from an ordinary stellar wind. The decreasing density occurs in regions where the magnetic field pressure becomes important, although the structure of such regions in pulsar nebulae remains uncertain.

Although the pulsar bubble model has some appeal, it does not apply to cases where the GRB occurs close in time to the supernova, as apparently was the case with GRB 030329 and SN 2003dh (e.g., Ma03 ).

## 5 Discussion and conclusions

Despite several years of effort, the study of GRB environments from their interaction has not clearly pointed to the progenitor objects. Reasons for this include the uncertainties in the basic model parameters, the fact that the GRB ejecta appear to be collimated and the similarity in the surrounding densities expected in different scenarios. The clearest progenitor information we have comes from the association of SN 2003dh with GRB 030329. The similarity of the supernova to SN 1998bw suggests that the progenitor object is a similar massive star. Analysis of the radio emission from SN 1998bw showed compatibility with expectations of a wind-like surrounding medium LC99 . The radio observations of GRB 030329 are slightly better fit by a uniform medium than a wind-like medium, although the difference between the fits is not large Be03 . The host galaxy of GRB 030329 appears to be a starburst dwarf galaxy Ma03  and the burst is positioned near the edge of the star forming region Fue03 , so it is not clear whether a high pressure surroundings is expected. Detailed modeling will be needed to determine whether the afterglow features can be explained by a complex explosion, or whether the surroundings need special properties.

An outstanding question is whether there is any need for a progenitor of long-duration bursts other than massive stars. Since massive stars are expected to modify their surroundings through winds, this evidence would be an incompatibility with the wind effects expected around a massive star. Frail et al. Fr03  and Yost et al. Ye03  have noted that a number of afterglows that can best be modeled as expanding into a uniform density medium with $`n1030`$ cm-3, and that this density is typical of Galactic interstellar clouds or the interclump medium of molecular clouds. However, such a medium would be modified by the winds from a massive star progenitor and a significant fraction of GRBs probably occur in starburst regions where clouds are denser than in Galactic case. One possibility is that the uniform medium is created by a shocked stellar wind in a high pressure medium.

A possible problem for massive star models is the low density inferred around some GRBs. Even if a shocked wind is present, it has higher density at a given radius than the corresponding free wind, so strong limits on the wind density can be set. Two afterglows with low densities are GRB 990123 PK02  and GRB 021211 KP03 , which require $`A_{}<10^310^2`$. Such low densities have not observed around Wolf-Rayet stars. GRB 021211 shows evidence for a supernova DVe03 , but the evidence is not conclusive. The low densities are not a problem for a burst that interacts directly with the hot interstellar medium in a galaxy. Another problem with the massive star models is the lack of evidence for bursts crossing the interface between a free wind and a shocked wind.

If direct interaction with the interstellar medium is required, a plausible progenitor object is a binary of compact objects. This requires that two different progenitor objects can give GRBs that appear similar, presumably from the formation of black holes. The application of the pulsar wind nebula model to some bursts requires both that different progenitor objects give rise to similar bursts, but also that the expansion of the burst into a different kind of medium (pair plasma and magnetic field) can give rise to similar afterglows. The finding of a supernova (SN 2003dh) occurring at approximately the same time as a GRB (GRB 030329) mitigates against the pulsar wind nebula picture for this case.

I am grateful to Jon Marcaide and the other conference organizers for their hospitality in Valencia, and to Zhi-Yun Li and Claes Fransson for continuing collaboration on these topics. This work was supported in part by NASA grant NAG5-13272 and NSF grant AST-0307366.

## Index


text:
# Galaxy Formation and Baryonic Dark Matter

## 1. Introduction: main problems of the $`\mathrm{\Lambda }`$CDM paradigm

One of the main problems encountered by the CDM-model has already been debated actively in this conference: the observed dark halo profiles of DM-dominated galaxies (dwarfs and LSB) are much better fitted with constant-density cores than with the steep power-law slopes predicted by CDM simulations (e.g. Navarro, Frenk & White 1997). Rotation curves have now sufficiently large spatial resolution (their dynamical range is even larger than in simulations), and it appears that it is no longer possible to explain the discrepancy by observational artefacts (de Blok et al. 2003). While the observed slope is in average $``$ -0.2, recent simulations show that the predicted slopes are not flatter than -1.5 (Fukushige & Makino 2003). The amount of discrepancy is however disputed (e.g. Swaters et al 2003).

A second problem is the low angular momentum of condensed baryons in CDM models, and consequent small predicted radius of galaxy disks, much too small compared to observed ones (Navarro & Steinmetz 2000). CDM simulations reproduce the Tully-Fisher slope and scatter, but not the zero point. They end up with dark halos too concentrated, and therefore a dark-to-visible mass ratio much too large within the observable disk.

A third problem is the high predicted number of small haloes, a power spectrum for structures advantaging the small-scales, leading to a large number of substructures that are not observed (Moore et al. 1999).

Can the hypothesis that a large fraction of dark baryons are in the form of cold condensed gas help to solve the problems? If we believe the primordial nucleosynthesis constraints, 90% of baryons are dark, and several hypothesis have been explored for the nature of these dark baryons; baryons in compact objects (brown dwarfs, white dwarfs, black holes) are either not favored by micro-lensing experiments or suffer major problems (Lasserre et al. 2000, Alcock et al. 2001, Afonso et al. 2003). The best hypothesis today is gas, either hot gas in the intergalactic and inter-cluster medium, or cold gas in the vicinity of galaxies (Pfenniger et al. 1994, Pfenniger & Combes 1994). The presence of large amounts of warm and diffuse gas in space cannot solve the problems, as shown by most $`\mathrm{\Lambda }`$CDM cosmological simulations based on this hypothesis, so we will concentrate on the second option.

## 2. Cusps in Galaxy centers

This central problem has been the subject of many efforts to solve it, either redistributing the dark matter or changing its nature, but up to now no one is perfectly satisfying: – black hole binaries can heat and flatten cusps (Milosavljevic & Merritt 2001), but they do not exist in dwarf galaxies, – stellar feedback (Gnedin & Zhao 2002), but even generous amounts of feedback are insufficient to destroy the central halo cusp, while the inner density is lowered only by a modest factor of 2 to 6; – bars (Weinberg & Katz 2002) have only weak effects and are not likely in dwarf Irr again. Among the various alternatives to the $`\mathrm{\Lambda }`$CDM model, warm dark matter WDM appears not promising (Knebe et al 2002), nor self-interacting dark matter SIDM (Spergel & Steinhardt 2000). Changing the cosmological parameters, matter density or power spectrum, allowing little power on galaxy scales (McGaugh et al 2003), bring only partial solutions.

Dwarf Irr galaxies are dominated by dark matter, but also their gas mass is dominating their stellar mass. It is remarkable that they obey the relation $`\sigma _{DM}`$/$`\sigma _{HI}`$ = cste, revealing that the dark-matter and HI gas surface densities have the same radial distribution (e.g. Hoekstra et al. 2001). All rotation curves can be explained, when the observed surface density of gas is multiplied by a constant factor, around 7-10. This $`\sigma _{DM}`$/$`\sigma _{HI}`$ ratio is almost constant from galaxy to galaxy, being slightly larger for early-types (Combes 2002). A solution to the problem of non-observed CDM cusps is that CDM would not be dominating in LSB galaxy centers, as is already the case in more evolved early-type galaxies, dominated by the stars. On the contrary, the mass in the center of LSB HI-dominated galaxies could be predominantly baryonic, in the form of cold condensed molecular gas. The column density of this H2 gas would be about one order of magnitude larger than the HI column density, which is already the case in early-type galaxies (where the larger metallicity of the gas allows to observe the H2 tracer, i.e. the CO molecule). That the gas mass has an essential dynamical role is supported by the success of the baryonic Tully-Fisher relation (McGaugh et al. 2000). It is interesting to note that dark haloes at cluster scale are also observed with a core (and not a cusp, Sand et al. 2002), and this could be due to the mass domination of baryons in the center (El-Zant et al 2003).

The influence of the baryon fraction $`\mathrm{\Omega }_b`$ on galaxy formation has been investigated through numerical simulations by Gardner et al. (2003): they show that increasing $`\mathrm{\Omega }_b`$ does not change the amount of diffuse gas in the intergalactic medium, but on the contrary the condensed phase associated to galaxies increases faster then $`\mathrm{\Omega }_b`$, due to the increased cooling rate. This leads to a considerably larger baryon representation in virialized systems with respect to their universal fraction. Another factor is the limited spatial resolution of the simulations, which could reduce the concentration of the dissipative material.

## 3. Angular momentum and disk formation

Most CDM cosmological simulations encounter a fundamental problem to form galaxy disks: due to low specific angular momentum for baryons, disks end up with too small radii (e.g. Abadi et al. 2003). Similarly elliptical systems end up too concentrated (Meza et al. 2003). In the commonly adopted paradigm, baryons at the start have the same specific angular momentum than dark matter. But then, in the course of evolution, baryons lose their angular momentum to the CDM, through dynamical friction. This is due to the main formation mechanism of galaxies, through hierarchical merging. The problem is the too early concentration of baryons, and the subsequent formation of galaxy disks as the final outcome of a sequence of merger events. Instead disks could be formed essentially through gas accretion from large-scale filaments, and this accretion could be prolonged until late times.

### 3.1. Feedback and cold gas accretion

The main solution to the angular momentum problem has been to increase the efficiency of feedback processes due to star formation. Supernovae should provide enough energy in the interstellar medium to prevent further collapse and star formation, and maintain disks with angular momentum close to the value required to fit observations (Thacker & Couchman 2001). In the same vein, Binney et al (2001) propose that, after the formation of protogalactic disks, massive winds expell considerable amounts of baryons. This contributes to the expansion of the inner dark matter (and the formation of cores), while the presently observed disks form late, when the dark halo contribution in the inner parts has been reduced.

Weil et al (1998) varied the gas cooling epoch, and found that discs can form by the present day with correct angular momenta if radiative cooling is suppressed until z=1. This means that feedback processes must be efficient enough to prevent gas from collapsing until late epochs. Eke et al. (2000) delay the gas cooling until z=1, and succeed to form disks with angular momentum approximately conserved during collapse. They show however that keeping enough specific angular momentum depends on the cosmology adopted, on the inhomogeneity of the haloes after the gas begins to collapse, and of course on the value of $`\mathrm{\Omega }_m`$, larger values of $`\mathrm{\Lambda }`$ for flat universes being much more favourable.

How does the gas collapse? In semi-analytical models (e.g. Mo et al. 1998), the gas is assumed to be hot and shock heated to the Virial temperature of the halo. But another way to accrete mass is cold gas mass accretion (Katz et al. 2002, Binney 2003). Gas is channeled through filaments, and only moderately heated by weak shocks, and radiating quickly. It never heats to the Virial temperature of its halo. Accretion is not spherical, gas keeps angular momentum and rotates near the galaxies, making it more easy to form disks. Cold mode of gas accretion is the most efficient at z larger than 1.

### 3.2. Gas accretion and bar strength

The bar frequency in nowadays galaxy disks can provide a way of quantifying gas accretion. It is well known that bars provoke their own destruction, by driving gas towards the center (e.g. Hasan & Norman 1990). After 5 Gyr of galaxy disk evolution, there should not be any bar left. Why do we observe so many bars today? more than 2/3 of galaxies are barred. There are even more bars detected in NIR images, when dust extinction does not blur them.

The only way to renew a bar in a galaxy is to accrete external cold gas. The gas will increase the disk mass, and reduce the bulge-to-disk mass ratio, destabilizing the disk again. Several bar episodes can process in a Hubble time (cf Fig 1; Bournaud & Combes 2002).

We have recently quantified the frequency of bars, in a sample of 163 galaxies imaged in the NIR by Eskridge et al. (2002). The gravity torques have been estimated by computing first the potential of the galaxy, through a Fourier transform of its NIR image. The intensity $`Q_2`$ of the $`m=2`$ torques is defined by the maximum over the disk of the ratio $`F_t(m=2)/F_r`$, with $`F_t`$ the tangential force, and $`F_r`$ the radial force. Similarly, the total torques are estimated by the maximum of $`F_t/F_r`$, with all $`m`$ included (Block et al. 2002). The surprising result is the dearth of axisymmetric galaxies (cf Fig. 2). Simulations of external gas accretion have been made in a sample of objects very similar in types to the observed sample. The frequency of bars was estimated in the simulated sample in the same way as for the observed sample. Comparison of the two histograms of torque intensity allows a quantification of the accretion rate (Fig. 2). It reveals that a galaxy doubles its mass in about 10 Gyr, through cold gas accretion.

### 3.3. Avoidance of dynamical friction

If the gas is condensed very early in galaxies, and is tightly bound to massive entities, it will experience strong dynamical friction at each merging stage, against the dark matter particles. A way to avoid this is to eject hot gas in supernovae feedback processes, but this might be insufficient (e.g. Gnedin & Zhao 2002, Maller & Dekel 2002).

Another scenario is that the gas accretes slowly in a cold phase onto galaxies, the hierarchical merging will lose less angular momentum through dynamical friction (cf. Fig 3). Gas in this scenario is only slightly bound to galaxies, and does not form stars too early. Effective accretion occurs later. In a merging, the gas is rapidly stripped, and does not experience friction.

## 4. Disruption of small structures

CDM models predict for a large-scale halo thousands of substructure clumps, which at cluster scale correspond to galaxies and are compatible with observations. However, this substructure feature is also expected at galactic scale, and raises many problems, when applied for example to our Galaxy. These small objects are not observed, and could be only dark halos (Chiu et al 2001), but then would yield to an excessive heating if present (Moore et al 1999). The assumption of some warm dark matter (WDM) helps to reduce the problem (Colin et al 2000). Sub-structures do form in WDM models, but with an average concentration parameter which is approximately twice as small as that of the corresponding CDM sub-halos. This difference in concentration leads to a higher satellite destruction rate in WDM models, and a better agreement with observations. Models of self-interacting dark matter (SIDM), with primordial velocity dispersion also reduce substructure and eliminate central cusps (Hogan & Dalcanton 2000), but at the expense of new problems and fine tuning.

In the scenario where baryons dominate the mass in the center of galaxies, there is more cold gas in dwarf haloes, and much less mass concentration. Baryonic clumps heat dark matter through dynamical friction and smooth any cusp in dwarf galaxies (El-Zant et al. 2003, in prep).

The material inside galaxies is more dissipative, more resonant, and more prone to disruption and merging. The high merging rate may change the mass function for low-mass galaxies (e.g. Mayer et al. 2001). The disruption of dwarf satellites explain the formation of the stellar halo of the Milky Way (Bullock et al. 2001). The solution of the cusp problem, leading to more fragile galaxies, easier to disrupt, can lead also to the solution of the substructure problem.

## 5. Dark Matter in galaxy clusters

The fact that the baryonic gas could be highly condensed in galaxies is supported by the high baryon fraction in clusters. In rich clusters, the mass of the hot gas in the intra-cluster medium (ICM) can be much larger than the stellar mass in galaxies. Most of the baryons are in the hot ICM, which represents nearly the baryon fraction of the matter in the universe $`f_b=\mathrm{\Omega }_b/\mathrm{\Omega }_m0.16`$. The relative radial distribution of dark and visible matter is reversed, the mass becomes more and more visible with radius, contrary to what occurs in galaxies (David et al. 1995, Ettori & Fabian 1999, Sadat & Blanchard 2001).

Moreover, the hot gas has a relatively high metallicity (1/3 solar), and most of the metals in a cluster comes from the ICM: there is about twice more Fe in the ICM than in galaxies (Renzini 1997, 2003).

Since the metals are synthetized in galaxies, this means that either part of the hot gas comes directly from galaxies (by stripping, or disruption), or that stellar winds and supernovae have enriched the ICM. In fact, both sources of metals must be there, since metals expelled by SNe would not be sufficient. There is the same MFe/LB in clusters and galaxies, implying the same processing in clusters than in the field. Clusters have not lost iron, nor accreted pristine material since the ratio between $`\alpha `$ elements (made in SNII) and iron is about solar in the ISM of all clusters (with no or little variation from cluster to cluster). This means that the metals come from normal stellar nucleosynthesis, and observations are compatible with the origin of the gas from the outskirts of galaxies.

## 6. Conclusion

The physics of the baryonic gas plays a crucial role in the formation of galaxies. The usual assumption that gas is shock heated to the virial temperature of the dark haloes might not be true. Instead, there could be cold gas accretion onto galaxies, with the consequence of more baryons accreted at a later time, and generally, more baryons condensed in galaxies, including dark baryons (Valageas et al. 2002). The main consequences are:

* baryons dominate in the center of galaxies, masking the dark matter cusps,
* existence of large gas extent around galaxies, and less angular momentum lost by dynamical friction,
* small galaxies are more prone to disruption and merging, which reduces substructures within massive galaxies.

## References

Abadi, M. G., Navarro, J. F., Steinmetz, M., Eke, V.R.: 2003, ApJ 591, 499

Afonso, C., Albert J.N., Andersen J. et al.: 2003, A&A 400, 951

Alcock, C., Allsman, R.A., Alves, D.R. et al.: 2001, ApJ 550, L169

Binney, J.: 2003, MNRAS, in press (astro-ph/0308172)

Binney, J., Gerhard, O., Silk, J.: 2001, MNRAS 321, 471

Block, D., Bournaud, F., Combes, F.. Puerari, I., Buta, R.: 2002, A&A 394, L35

Bournaud, F., Combes, F.: 2002, A&A 392, 83

Bullock, J. S., Kravtsov, A. V., Weinberg, D. H.: 2001, ApJ 548, 33

Chiu, W. A., Gnedin, N. Y., Ostriker, J. P.: 2001, ApJ 563, 21

Col n, P., Avila-Reese, V., Valenzuela, O.: 2000, ApJ 542, 622

Combes, F.: 2002 NewAR 46, 755

David, L. P., Jones, C., Forman, W.: 1995, ApJ 445, 578

de Blok, W. J. G., Bosma, A., McGaugh, S.: 2003, MNRAS 340, 657

Eke, V., Efstathiou, G., Wright, L.: 2000, MNRAS 315, L18

El-Zant A., Hoffman Y., Primack J., Combes F., Shlosman I.: 2003, ApJL, submitted (astro-ph/0309412)

Eskridge, P. B., Frogel, J. A., Pogge, R. W. et al. 2002, ApJS 143, 73-112

Ettori, S., Fabian, A. C.: 1999 MNRAS 305, 834

Fukushige, T., Makino, J.: 2003, ApJ 588, 674

Gardner, J. P., Katz, N., Hernquist, L., Weinberg, D. H.: 2003, ApJ 587, 1

Katz, N., Keres, D., Dave, R., Weinberg, D.H.: 2002, in The IGM/Galaxy Connection: Kluwer, p. 185 (astro-ph/0209279)

Knebe, A., Devriendt, J. E. G., Mahmood, A., Silk, J.: 2002, MNRAS 329, 813

Gnedin O.Y., Zhao H.: 2002, MNRAS 333, 299

Hasan, H., Norman, C.: 1990 ApJ 361, 69

Hoekstra, H., van Albada, T. S., Sancisi, R.: 2001, MNRAS 323, 453

Hogan, C. J., Dalcanton, J. J.: 2000, Phys Rev D 62, 3511

Lassere, T., Afonso, C., Albert J.N. et al. : 2000, A&A 355, L39

Maller A.H., Dekel A.: 2002, MNRAS 335, 487

Mayer, L., Governato, F., Colpi, M. et al. 2001 ApJ 547, L123

McGaugh, S. S., Barker, M. K., de Blok, W. J. G.: 2003, ApJ 584, 566

McGaugh, S. S., Schombert, J. M., Bothun, G. D., de Blok, W. J. G.: 2000, ApJ 533, L99

Meza, A., Navarro, J. F., Steinmetz, M., Eke, V. R.: 2003, ApJ 590, 619

Milosavljevic, M., Merritt, D.: 2001, ApJ 563, 34

Mo, H. J., Mao, S., White, S. D. M.: 1998 MNRAS 295, 319

Moore, B., Ghigna, S., Governato, F., et al.: 1999, ApJ 524, L19

Navarro, J.F., Frenk, C.S., White, S.D.M.: 1997, ApJ 490, 493

Navarro, J. F., Steinmetz, M.: 2000, ApJ 538, 477

Pfenniger, D., Combes, F., Martinet, L.: 1994, A&A 285, 79

Pfenniger, D., Combes, F.: 1994, A&A 285, 94

Renzini, A.: 1997, ApJ 488, 35

Renzini, A.: 2003, in Clusters of Galaxies: Probes of Cosmological Structure and Galaxy Evolution, Carnegie Observatories Symposium III. (astro-ph/0307146)

Sadat, R., Blanchard, A.: 2001, A&A 371, 19

Sand, D.J., Treu, T., Ellis, R.S.: 2002, ApJ 574, L129

Spergel, D. N., Steinhardt, P. J.: 2000, Phys Rev Let 84, 3760

Swaters, R. A., Madore, B. F., van den Bosch, F. C., Balcells, M.: 2003, ApJ 583, 732

Thacker, R. J., Couchman, H. M. P.: 2001, ApJ 555, L17

Valageas, P., Schaeffer, R., Silk, J.: 2002, A&A 388, 741

Weil M.L., Eke, V. R., Efstathiou, G.: 1998, MNRAS 300, 773

Weinberg, M. D., Katz, N.: 2002, ApJ 580, 627


text:
# Lattice-QCD-based Schwinger-Dyson Approach for Chiral Symmetry Restoration at Finite Temperature

## 1 Introduction

Dynamical chiral-symmetry breaking (DCSB) is one of the most important nonperturbative features in QCD. For the study of DCSB, the Schwinger-Dyson (SD) formalism  has been used as an interesting and powerful method. For QCD, the SD formalism consists of an infinite series of nonlinear integral equations which determine the $`n`$-point Green’s function of quarks and gluons, and therefore it includes the infinite-order effect of the QCD gauge coupling constant $`g`$.

For instance, the SD equation for the quark propagator $`S(p)`$ is described with the nonperturbative gluon propagator $`D_{\mu \nu }(p)`$ and the nonperturbative quark-gluon vertex $`g\mathrm{\Gamma }_\nu (p,q)`$ as

$`S^1(p)=S_0^1+g^2{\displaystyle _q}\gamma _\mu S(q)D_{\mu \nu }(pq)\mathrm{\Gamma }_\nu (p,q),`$ (1)
where $`S_0(p)`$ denotes the bare quark propagator and the simple notation $`_q\frac{d^4q}{(2\pi )^4}`$ has been used in the Euclidean metric.

In the practical calculation for QCD, however, the SD formalism is drastically truncated: the perturbative gluon propagator and the one-loop running coupling are used instead of the nonperturbative quantities in the original formalism. This simplification seems rather dangerous because some of nonperturbative-QCD effects are neglected.

In this paper, we formulate the SD equation based on the recent lattice QCD (LQCD) results, i.e., the LQCD-based SD equation, and aim to construct a useful and reliable analytic framework including the proper nonperturbative effect in QCD. Using the LQCD-based SD equation, we investigate also DCSB at finite temperatures.

## 2 The Quark Propagator in Lattice QCD

First, we briefly review the quark propagator $`S(p)`$ in lattice QCD. The inverse quark propagator in the Landau gauge is generally given by

$`S^1(p)=Z(p^2)/\{\mathit{}+M(p^2)\}`$ (2)
in the Euclidean metric. Here, $`M(p^2)`$ is called as the quark mass function, and $`Z(p^2)`$ corresponds to the wave-function renormalization of the quark field. In the quark propagator, DCSB is characterized by the mass generation as $`M(p^2)0`$.

The quark mass function $`M(p^2)`$ in the Landau gauge is recently measured in lattice QCD at the quenched level , and the lattice data in the chiral limit is well reproduced by

$`M(p^2)=M_0/\{1+(p/\overline{p})^\gamma \}`$ (3)
with $`M_0`$=260 MeV, $`\overline{p}`$=870MeV and $`\gamma `$=3.04. The infrared quark mass $`M(0)=M_0260\mathrm{M}\mathrm{e}\mathrm{V}`$ seems consistent with the constituent quark mass in the quark model. Using this lattice result of $`M(p^2)`$, the pion decay constant is calculated as $`f_\pi `$ 87 MeV with the Pagels-Stokar formula , and the quark condensate is obtained as $`\overline{q}q_{\mathrm{\Lambda }=1\mathrm{G}\mathrm{e}\mathrm{V}}(220\mathrm{M}\mathrm{e}\mathrm{V})^3`$. These quantities related to DCSB seem consistent with the standard values.

## 3 The Schwinger-Dyson Equation

In this section, we formulate the SD equation for quarks in the chiral limit in the Landau gauge. By taking the trace Eq.(1), one finds

$`{\displaystyle \frac{M(p^2)}{Z(p^2)}}={\displaystyle \frac{g^2}{4}}{\displaystyle _q}\mathrm{tr}\{\gamma _\mu {\displaystyle \frac{Z(q^2)}{\mathit{}+M(q^2)}}\mathrm{\Gamma }_\nu (p,q)\}D_{\mu \nu }(\stackrel{~}{q})`$ (4)
with $`\stackrel{~}{q}pq`$. For the quark-gluon vertex, we assume the chiral-preserving vector-type vertex,

$`\mathrm{\Gamma }_\mu (p,q)=\gamma _\mu \mathrm{\Gamma }((pq)^2),`$ (5)
which keeps the chiral symmetry properly. (In contrast, to be strict, the Higashijima-Miransky approximation  explicitly breaks the chiral symmetry in the formalism.) Then, one obtains

$`{\displaystyle \frac{M(p^2)}{Z(p^2)}}=C_Fg^2{\displaystyle _q}{\displaystyle \frac{Z(q^2)M(q^2)}{q^2+M^2(q^2)}}\mathrm{\Gamma }(\stackrel{~}{q}^2)D_{\mu \mu }(\stackrel{~}{q})`$ (6)
with $`C_F=4/3`$ being the color factor for quarks. In the Landau gauge, the Euclidean gluon propagator is generally expressed by

$`D_{\mu \nu }(p^2)={\displaystyle \frac{d(p^2)}{p^2}}\left(\delta _{\mu \nu }{\displaystyle \frac{p_\mu p_\nu }{p^2}}\right),`$ (7)
where we refer to $`d(p^2)`$ as the gluon polarization factor. Therefore, Eq.(6) is rewritten as

$`{\displaystyle \frac{M(p^2)}{Z(p^2)}}=3C_Fg^2{\displaystyle _q}{\displaystyle \frac{Z(q^2)M(q^2)}{q^2+M^2(q^2)}}\mathrm{\Gamma }(\stackrel{~}{q}^2){\displaystyle \frac{d(\stackrel{~}{q}^2)}{\stackrel{~}{q}^2}}.`$ (8)
Here, we define the kernel function

$`K(p^2)g^2\mathrm{\Gamma }(p^2)d(p^2)`$ (9)
as the product of the quark-gluon vertex $`\mathrm{\Gamma }(p^2)`$ and the gluon polarization factor $`d(p^2)`$. Then, the SD equation is expressed as

$`{\displaystyle \frac{M(p^2)}{Z(p^2)}}=3C_F{\displaystyle _q}{\displaystyle \frac{Z(q^2)M(q^2)}{q^2+M^2(q^2)}}{\displaystyle \frac{K((pq)^2)}{(pq)^2}}.`$ (10)
In the Landau gauge, the quark wave-function renormalization is not so significant and seems to be approximated as $`Z(p^2)=1`$, which reduces the SD equation to

$`M(p^2)=3C_F{\displaystyle _q}{\displaystyle \frac{M(q^2)}{q^2+M^2(q^2)}}{\displaystyle \frac{K((pq)^2)}{(pq)^2}}.`$ (11)
## 4 Extraction of the Kernel Function in the SD Equation from Lattice QCD

In this section, we extract the kernel function $`K(p^2)g^2\mathrm{\Gamma }(p^2)d(p^2)`$ in the SD equation (11) using the quark mass function $`M(p^2)`$ obtained in lattice QCD. By shifting the integral variable from $`q`$ to $`\stackrel{~}{q}pq`$, we rewrite Eq.(11) as

$`M(p^2)=3C_F{\displaystyle _{\stackrel{~}{q}}}{\displaystyle \frac{M((p\stackrel{~}{q})^2)}{(p\stackrel{~}{q})^2+M^2((p\stackrel{~}{q})^2)}}{\displaystyle \frac{K(\stackrel{~}{q}^2)}{\stackrel{~}{q}^2}}.`$ (12)
Therefore, we obtain

$`M(p^2)={\displaystyle \frac{3C_F}{8\pi ^3}}{\displaystyle _0^{\mathrm{}}}𝑑\stackrel{~}{q}^2\mathrm{\Theta }(p,\stackrel{~}{q})K(\stackrel{~}{q}^2),`$ (13)
where $`\mathrm{\Theta }(p,q)`$ is defined with $`M(p^2)`$ as

$`\mathrm{\Theta }(p,q){\displaystyle _0^\pi }𝑑\theta \mathrm{sin}^2\theta `$
$`{\displaystyle \frac{M(p^2+q^22pq\mathrm{cos}\theta )}{p^2+q^22pq\mathrm{cos}\theta +M^2(p^2+q^22pq\mathrm{cos}\theta )}}.`$ (14)
Since the quark mass function $`M(p^2)`$ is given by Eq.(3) in lattice QCD, we can calculate the kernel function $`K(p^2)`$ from Eq.(13).

As shown in Fig.1, we numerically obtained the kernel function $`K(p^2)g^2\mathrm{\Gamma }(p^2)d(p^2)`$ extracted from the lattice QCD result for the quark propagator in the Landau gauge.

As remarkable features, we find “infrared vanishing” and “intermediate enhancement” in the kernel function $`K(p^2)`$ in the SD equation. In fact, $`K(p^2)`$ seems consistent with zero in the very infrared region as

$`K(p^20)0,`$ (15)
while $`K(p^2)`$ exhibits a large enhancement in the intermediate-energy region around $`p`$ 0.5GeV.

These tendencies of infrared vanishing and intermediate enhancement in the kernel function $`K(p^2)g^2\mathrm{\Gamma }(p^2)d(p^2)`$ are observed also in the direct lattice-QCD measurement for the polarization factor $`d(p^2)`$ in the gluon propagator in the Landau gauge .

Note here that the usage of the perturbative gluon propagator and the one-loop running coupling $`g_{\mathrm{run}}(p^2)`$ corresponds to $`K(p^2)=g_{\mathrm{run}}^2(p^2)`$, which largely differs from the present result based on lattice QCD both in the infrared and in the intermediate-energy regions. In fact, the simple version of the SD equation using the perturbative gluon propagator and the one-loop running coupling would be too crude for the quantitative study of QCD.

## 5 Chiral Symmetry at Finite Temperature

Finally, we demonstrate a simple application of the LQCD-based SD equation to chiral symmetry restoration in finite-temperature QCD.

At a finite temperature $`T`$, the field variables obey the (anti-)periodic boundary condition in the imaginary-time direction, which leads to the SD equation for the thermal quark mass $`M_n(𝐩^2)`$ of the Matsubara frequency $`\omega _n(2n+1)\pi T`$ as

$`M_n(𝐩^2)={\displaystyle _{m,𝐪}}{\displaystyle \frac{3C_FM_m(𝐪^2)}{\omega _m^2+𝐪^2+M_m^2(𝐪^2)}}{\displaystyle \frac{K(\omega _{nm}^2+\stackrel{~}{𝐪}^2)}{\omega _{nm}^2+\stackrel{~}{𝐪}^2}},`$ (16)
with $`_{m,𝐪}T_{m=\mathrm{}}^{\mathrm{}}\frac{d^3q}{(2\pi )^3}`$, $`\omega _{nm}\omega _n\omega _m`$ and $`\stackrel{~}{𝐪}𝐩𝐪`$.

Using the kernel function $`K(p^2)`$ obtained in the previous section, we solve Eq.(16) for the thermal quark mass $`M_n(𝐩^2)`$. Figure 2 shows the preliminary result for the thermal infrared quark mass $`M_0(𝐩^2=0)`$ plotted against the temperature $`T`$. We thus find chiral symmetry restoration at a critical temperature $`T_c100\mathrm{M}\mathrm{e}\mathrm{V}`$.


text:
# Polarization effects in diffraction of light on a planar chiral structure

## Abstract

We analyzed polarization changes of light diffracted on planar chiral array from the standpoint of the Lorentz reciprocity lemma and found the bi-orthogonality in the polarization eigenstates for the waves diffracting though the grating in the opposite direction. Both reciprocal and non-reciprocal component in the polarization azimuth rotation of the diffracted light were identified. The structural chirality of the array arrangement and the chirality of individual elements of the array give rise to polarization effects.

planar chirality, nonreciprocity, diffraction

Recently we reported that planar (2D) chiral structures affect the polarization state of light in an enantiomeric fashion, similarly to three–dimensional chiral media Papakostas et al. (2003). However polarization phenomena on diffraction from planar chiral structures have never been studied theoretically before, leaving the fundamental properties of 2D chirality not fully understood. Here we report on the results of a theoretical investigation of polarization changes for light diffracted on regular arrays of planar chiral metallic structures from the standpoint of the Lorentz reciprocity theorem. By analyzing propagation of light in two opposite directions we have identified a strong component in polarization effect on diffraction that can be induced either by the chirality of the individual elements of the array, or by arranging non-chiral elements of an array in a chiral fashion.

Let us consider a planar square periodic array of metallic elements of thickness $`t`$ with equal pitch $`d`$ along the axes $`x`$ and $`y`$ placed between planes $`z=0`$ and $`z=t`$ (see Fig. 1 and Fig. 2). If a plane electromagnetic wave

$$𝐄_i=𝐀_ie^{i𝐤_i𝐫}$$
(1)
of unit amplitude and polarization vector $`𝐀_i`$ is incident on the array from the region corresponding to $`z>0`$, the transmitted field may be written as a summation over all diffracted waves, numbered by integer indices $`q`$ and $`p`$

$$𝐄_t=\underset{q,p=\mathrm{}}{\overset{\mathrm{}}{}}𝐚_{qp}e^{i𝐤_{qp}(𝐫+𝐞_zt)},z<t$$
(2)
where $`𝐚_{qp}`$ and $`𝐤_{qp}`$ are amplitudes and wave-vectors partial of diffracted waves and

$`𝐤_{qp}=𝐠+𝐡_{qp}𝐞_z\sqrt{k^2|𝐠+𝐡_{qp}|^2},`$
$`𝐡_{qp}=2\pi (q𝐞_x+p𝐞_y)/d.`$
Here $`𝐞_x`$, $`𝐞_y`$ and $`𝐞_z`$ are unit vectors along the axes $`x`$, $`y`$ and $`z`$, $`𝐠`$ is the component of $`𝐤_i`$ transverse to the axis $`z`$ and $`k=|𝐤_i|`$. Let us now consider a ”reversed” wave with polarization vector $`𝐀_r`$ approaching the array from the opposite side of the structure ($`z<t`$) along the direction of one of the partial diffracted waves of the ”direct” scenario, with indices $`s`$ and $`l`$ and the wave vector $`𝐤_r=𝐤_{sl}`$ ,

$$𝐄_r=𝐀_re^{i𝐤_r(𝐫+𝐞_zt)}.$$
(3)
In the region $`z>0`$ this wave will produce diffracted waves with amplitudes $`𝐛_{qp}`$ . This corresponds to the ”reversed” scenario of diffraction. The Lorentz reciprocity lemma Kong (1986) applied to the field superposition in the volume bounded by surface $`S`$ that consists of planes $`x=\pm d/2`$, $`y=\pm d/2`$, $`z=z_1>0`$, and $`z=z_2<t`$ may be written in the following form:

$$_S\left\{[\stackrel{~}{𝐄}_i\times \stackrel{~}{𝐇}_r][\stackrel{~}{𝐄}_r\times \stackrel{~}{𝐇}_i]\right\}𝑑𝝈=0$$
(4)
where $`\stackrel{~}{𝐄}_i`$, $`\stackrel{~}{𝐇}_i`$ and $`\stackrel{~}{𝐄}_r`$, $`\stackrel{~}{𝐇}_r`$ are electric and magnetic fields created by waves incident from opposite directions. By using corresponding field expressions it may be shown from formula (4) that

$$\sqrt{k^2|𝐠|^2}(𝐀_i𝐛_{sl})=\sqrt{k^2|𝐠+𝐡_{sl}|^2}(𝐀_r𝐚_{sl}).$$
(5)
The equality (5) constitutes the universal relation between the amplitudes of partial waves in the ”direct” and ”reversed” diffraction scenarios. Scattering processes are often described in terms of $`2\times 2`$ transformation matrices, relating Cartesian components of electric fields in coordinate frames of incident and scattered waves. For the ”direct” ($`\widehat{D}`$) and ”reversed” ($`\widehat{R}`$) scenarios these matrices for the incident and partial diffracted waves can be introduced as follows: $`𝐚_{sl}=\widehat{D}𝐀_i,`$ $`𝐛_{sl}=\widehat{R}𝐀_r.`$ It may be shown from equation (5) that these matrices are linearly related and mutually transposed:

$$R_{nm}=c(2\delta _{mn}1)D_{mn},$$
(6)
where $`\delta _{mn}`$ is Kroneker index and $`c=\sqrt{k^2|𝐠+𝐡_{sl}|^2}/\sqrt{k^2|𝐠|^2}`$. For the purpose of analyzing the polarization eigenstates of the diffraction process it is instructive to present both scattering matrices in the coordinate frame of the direct scenario (see Fig. 1,a) where the operator of the reversed scattering process acts on the complex-conjugated field amplitudes. Here the polarization eigenstates are simply two linearly independent eigenvectors of matrices $`\widehat{R}^{}`$ and $`\widehat{D}`$. The relation between them may be derived from equation (6), which, when converted to the coordinate frame of the direct scenario, gives $`R_{nm}^{}=cD_{mn}^{}`$. It follows from the theory of matrix operators that eigenvectors of the Hermitian-conjugated matrices with elements $`D_{nm}`$ and $`D_{mn}^{}`$, and therefore of matrices $`\widehat{R}^{}`$ and $`\widehat{D}`$, are biorthogonal or, in terms of polarization eigenstates, are represented by antipode points on the Poincare Sphere, as shown on Fig. 3. In general, the point representing the first eigenstate in the direct scenario $`1d`$ is an antipode to one of the points, which represents the eigenstates of the reversed scenario (this point is designed as $`2r`$ in Fig. 3). However eigenstate $`1d`$ does not necessary coincide with eigensate $`1r`$ of the reversed scenario, nor eigenstate $`2d`$ coincide with $`2r`$. Therefore, the polarization eigenstates in the direct and reversed scenarios presented in the coordinate frame of the direct scenario could be different. Such a situation takes place if the complex matrix $`\widehat{D}`$ is an asymmetric or even a non-diagonal matrix.

We found that scattering matrices of non-zero order diffraction on periodic plano-chiral arrays, where chirality is due to either structural chirality, or chirality of individual elements of the array, are either asymmetric or non-diagonal. Below we will illustrate these properties by numerical modelling of the diffraction process for various planar chiral gratings. We calculated the fields and polarization characteristics of light diffracted on gratings numerically using the method described in reference Prosvirnin (1999). It is based on a vector integral equation for the surface current induced by the light wave on the array particles. The equation is derived with boundary conditions for ideal metallic structures that assume a zero value for the tangential component of the electric filed on the metal. The integral equation is then reduced to an algebraic equation set by use of the Galerkin technique.

In our modelling we concentrated on planar chiral arrays of 442 symmetry wallpaper group and calculated the polarization eigenstates of the diffraction process and polarization changes occurring in the diffracted wave for different incident polarizations.

We studied diffraction for two different incident angles $`\beta `$, at $`\beta _1=0`$ and at $`\beta _2=\mathrm{arcsin}(\pi /kd)`$. In the first case, the diffracted wave ($`q=1`$, $`p=0`$) propagates at angle $`\xi =\xi _1=\mathrm{arcsin}(2\pi /kd)`$ to the array. In the second case the diffracted wave ($`q=1`$, $`p=0`$) the same angle makes to the array as the incident wave $`\xi =\xi _2=\beta _2`$ (for definitions of angles see Fig.1). The wave’s polarization azimuth $`\theta `$ and degree of ellipticity $`\eta `$ were calculated from the Cartesian field amplitudes using the standard definitions: $`\mathrm{tan}2\theta =s_2/s_1`$, $`\mathrm{sin}2\eta =s_3/s_0`$, where $`s_i`$ are the Stokes parameters. The results of our analysis for $`d=4\mu \text{m}`$, $`\lambda =2\pi /k=1520`$ nm, $`\beta _1=0`$, $`\xi _1=22.3^{}`$, $`\beta _2=11.0^{}`$, $`\xi _2=11.0^{}`$ are summarized in Table 1 and Fig. 4. We considered an array without a substrate. The width of the metal strips was equal to 0.05 $`\mu `$m.

In optics polarization elements are often classified as reciprocal or non-reciprocal depending on whether their effect on the polarization state of the transmitted light is the same or different for light propagating in the opposite directions. This understanding of optical reciprocity which we will use below is somewhat different from the general, more tolerant definition of reciprocity based on the Lorentz lemma. For the purpose of comparison of the polarization transformations for opposite directions of light propagation, in the table and figures the polarization parameter of the waves are converted into the coordinate frame of the direct scenario. In such a presentation if the values of polarization azimuth rotation in the direct and reversed scenarios are the same, the rotation is truly non-reciprocal. On contrary, a difference between the values of polarization azimuth rotation in the direct and reversed scenarios would represent a reciprocal component of the polarization change, that is analogous to the optical activity effect in a chiral liquid.

The calculations revealed that:

i) For all diffraction processes involving twisted or non-twisted arrays, equalities (6) are held to within the numerical accuracy of the method. They are thus compatible with the Lorentz lemma.

ii) For the arrays of straight crosses polarization azimuth rotation in opposite directions have opposite signs due to the difference in the efficiency of diffraction for perpendicular polarization components (line C on Fig. 4(B)). This non-reciprocity of polarization azimuth rotation is analogous to the polarization rotation non-reciprocity in dichroic media due to anizotropic dissipation.

iii) No polarization rotation is seen in the non-diffracted part of the beam at the normal incidence. It’s polarization eigenstates are the same in both directions and for any type of array.

iv) From Table I, one can see that for non-zero order asymmetrical diffraction ($`\xi \beta `$) when individual structural elements of the array are twisted rosettes polarization azimuths of eigenstates for the ”direct” and ”reversed” scenarios are resolutely different. The difference between the polarization azimuths of the eigensatates depends on the rosette curvature angle $`\phi `$ and reaches a maximum of about $`57^{}`$ at $`\phi =120^{}`$. The difference in the eigenstates vanishes at rosette bending angle $`\phi =95^{}`$.

v) When individual structural elements of the array are twisted rosettes, polarization azimuth rotation on diffraction has both reciprocal and nonreciprocal components. The non-reciprocal component of the polarization azimuth rotation is due to a difference in the efficiency of diffraction for perpendicular polarization components. The corresponding oscillating dependence of the nonreciprocal rotation on the incident angle is shifted in respect of line C corresponding to straight crosses. It is shifted along the incident polarization azimuth axis towards left for left rosettes and towards right for right rosettes. The split between corresponding solid and dashed lines in Fig. 4 indicates the reciprocal component of the polarization azimuth rotation analogous to optical activity.

vi) Nonreciprocity of polarization rotation in the diffraction process is evident when a diffracted light wave is reflected straight back towards the twisted planar structure by a mirror, and then diffracts again. The polarization state of the returning light after the second diffraction is different from that of the incident light, even if the incident light was an eigenstate in the forward direction. For an array of rosettes with $`\phi =120^{}`$, the two incident eigenstates and corresponding returning polarizations have azimuths different by $`27^{}`$ and $`21^{}`$.

Therefore, polarization effects on diffraction from planar chiral grating can be induced by either structural chirality or the chirality of individual elements of the array. However, no polarization rotation compatible with the Lorentz lemma is possible for a wave transmitting through or reflected from a planar chiral structure at normal incidence, as the scattering matrices are diagonal in this case.

Finally, we shall note that our analysis is underpinned by the Lorentz lemma while our computational method is compatible with it. It shall be noted, however, that the recent observation of broken parity and time reversal evident in polarized optical images of planar chiral structures Schwanecke et al. (PRL, in press) calls for re-examining the validity of the Lorentz lemma for planar chiral structure.

###### Acknowledgements.

The authors thank A.Papakostas, A.Potts, D.Bagnall and K.MacDonald for fruitful discussions and acknowledge the support of the Science and Engineering Research Council (UK).