format: text

text:
# The high energy tail of 𝛾-ray burst 941017: Comptonization of synchrotron self absorbed photons

## 1 Introduction

In fireball models of GRBs (Piran, 2000; Mészáros, 2002; Waxman, 2003) the energy released by an explosion is converted to kinetic energy of a thin baryonic shell expanding at an ultra-relativistic speed. The GRB is most likely produced by internal shocks within the expanding shell. At a later stage, the shell impacts on surrounding gas, driving an ultra-relativistic shock into the ambient medium. This shock continuously heats fresh gas and accelerates relativistic electrons, which produce synchrotron emission that account for the X-ray, optical and radio emission (the ”afterglow”) following the GRB. The initial interaction of fireball ejecta with surrounding gas produces a reverse shock which propagates into and decelerates the fireball ejecta (Mészáros, Rees & Papathanassiou, 1994). As the reverse shock crosses the ejecta, it erases the memory of the initial conditions, and the expansion then approaches self-similarity (Blandford & McKee, 1976), where the expansion Lorentz factor decreases with radius, $`\mathrm{\Gamma }_{\mathrm{BM}}=(17E/16\pi nm_pc^2)^{1/2}r^{3/2}`$. Here $`E`$ is the total explosion energy (assuming spherical symmetry), $`n`$ is the number density of the ambient medium.

The duration $`T`$ of the stage of transition to self-similar expansion, during which the reverse shock ”lives”, is comparable to the longer of the two time scales set by the initial conditions (Waxman, 2003): The (observed) GRB duration $`T_{\mathrm{GRB}}`$ and the (observed) time $`T_\mathrm{\Gamma }`$ at which the self-similar Lorentz factor $`\mathrm{\Gamma }_{\mathrm{BM}}`$ equals the original ejecta Lorentz factor $`\mathrm{\Gamma }_i`$, $`\mathrm{\Gamma }_{\mathrm{BM}}(T_\mathrm{\Gamma })=\mathrm{\Gamma }_i`$. Since the characteristic time over which radiation emitted by the fireball at radius $`r`$ is observed by a distant observer is $`r/4\mathrm{\Gamma }^2c`$ (Waxman, 1997), $`T`$ is determined by (Waxman, 2003)

$$T\mathrm{max}[T_{\mathrm{GRB}},10\left(\frac{E_{53}}{n_1}\right)^{1/3}\left(\frac{\mathrm{\Gamma }_i}{300}\right)^{8/3}\mathrm{s}],$$
(1)
where $`E=10^{53}E_{53}`$ erg and $`n=0.1n_1\mathrm{cm}^3`$. Note, that the duration is increased by a factor $`1+z`$ for a burst at redshift $`z`$.

Observations of GRB941017 show two distinct spectral components (González *et al.*, 2003): A low energy component, with photon energies $`\epsilon _\gamma <\mathrm{\hspace{0.33em}3}`$ MeV, and a high energy component, $`\epsilon _\gamma >\mathrm{\hspace{0.33em}3}`$ MeV. The low energy component shows rapid variability, has a characteristic GRB spectrum peaking at $`\epsilon _\gamma 0.5`$ MeV, and decays over $`100`$ s. The high energy component has a very hard spectrum, number of photons per unit photon energy $`dn_\gamma /d\epsilon _\gamma \epsilon _\gamma ^1`$, and persists over $`200`$ s. The different temporal behavior suggests that the two components are produced in different regions of the expanding fireball. The characteristics of the low energy component suggest that it is produced by internal shocks, similar to other GRBs. The temporal behavior of the high energy component suggests that it is produced during the transition to self-similar expansion (Granot & Guetta, 2003).

The hard, $`dn_\gamma /d\epsilon _\gamma \epsilon _\gamma ^1`$, spectrum is difficult to account for in models where emission is dominated by shock accelerated electrons (González *et al.*, 2003). This has lead González et al. to suggest that the high energy tail is due to electromagnetic cascades initiated by the interaction of photons with ultra-high energy shock-accelerated protons (Waxman, 1995; Vietri, 1995; Bötcher & Dermer, 1998; Totani, 1998; Waxman & Bahcall, 2000). We present here an alternative explanation: Electrons accelerated in the forward shock inverse-Compton scatter optical photons, emitted by the reverse shock electrons, to create the observed spectra. A key point, which allows to reproduce the observed hard spectrum, is the modification of the synchrotron spectrum by self-absorption in the reverse shock. We show that this effect allows to reproduce the observed hard, high energy tail also in the internal shock phase (see § 4). However, we consider the latter explanation less likely, due to the weak time dependence of the high energy component.

Granot & Guetta (2003) have recently considered inverse-Compton emission from reverse shock electrons during the transition to self-similarity, as an explanation to the high energy tail of GRB941017. They have found that in order for such an explanation to be viable, the Lorentz factor associated with fireball expansion should be higher, $`\mathrm{\Gamma }10^4`$, and the magnetic field in the fireball plasma should be much lower (well below equipartition), compared to values typically inferred from early afterglow observations (see, e.g., Zhang, Kobayashi & Mészáros, 2003). We show here that the high energy tail may be explained as emission from the forward shock electrons, with fireball plasma parameters which are typical to those inferred from GRB observations: $`\mathrm{\Gamma }10^2`$ and magnetic field close to equipartition (within the fireball plasma). In this scenario, the inferred density of plasma surrounding the fireball is lower than inferred for most other GRBs (the magnetic field strength in the shock driven by the fireball into the surrounding gas is poorly constrained by current observations). We consider the latter scenario more likely, since it requires a modification of the parameters of the environment external to the fireball, rather than modifications of the fireball physics. Our analysis further improves on that of Granot & Guetta (2003) in including the effects of self-absorption, and in carrying out detailed (numerical) calculations of the spectra, which are necessary given the inferred parameter range (see below).

This paper is organized as follows. In § 2 we briefly discuss the dynamics of transition to self-similar expansion and the plasma conditions during the transition. In § 3 we analytically derive the constraints that should be satisfied by model parameters in order to allow an explanation of the observed high energy tail as inverse-Compton emission from forward shock electrons during this phase. In § 4 we present the results of detailed numerical calculations of the spectrum, which demonstrate that the observed spectrum may be reproduced when the constraints derived in § 3 are satisfied. Such calculations are necessary since the spectral shape near the inverse-Compton up-scattered self-absorption frequency is not well described by simple power-law approximations. Our conclusions are summarized in § 5.

## 2 The transition to self-similar expansion

During the transition, plasma shocked by the reverse shock expands with Lorentz factor close to that given by the self-similar solution,

$$\mathrm{\Gamma }140\left(\frac{E_{53}}{n_1}\right)^{1/8}T_2^{3/8},$$
(2)
where $`T=10^2T_2\mathrm{s}`$. The unshocked fireball ejecta propagate with the original expansion Lorentz factor, $`\mathrm{\Gamma }_i>\mathrm{\Gamma }`$. The forward shock which propagates into the ambient medium is highly relativistic, propagating with Lorentz factor $`\mathrm{\Gamma }_{\mathrm{BM}}`$, while the reverse shock is mildly relativistic, propagating into the ejecta with Lorentz factor $`1+\mathrm{\Gamma }_i/\mathrm{\Gamma }_{\mathrm{BM}}`$.

The radiation luminosity and spectrum are determined, for given $`E`$, $`n`$ and $`\mathrm{\Gamma }_i`$, by the fraction $`ϵ_B`$ ($`ϵ_e`$) of shock thermal energy carried by magnetic field (relativistic electrons), and by the shape of the electron distribution function, which is commonly assumed to be a power law of index $`p\mathrm{d}\mathrm{ln}n_e/\mathrm{d}\mathrm{ln}\gamma _e`$ for electron Lorentz factors $`\gamma _e`$ exceeding $`\gamma _{em}`$. The values of the parameters $`ϵ_B`$ and $`ϵ_e`$ can not be determined from first principles, since the micro-physics of collisionless shocks which determines them is not fully understood. We therefore leave them as free parameters in the following analysis, and derive the constraints imposed on them by observations. Since the forward shock is highly relativistic, while the reverse shock is only mildly relativistic, we allow different values for $`ϵ_B`$ and $`ϵ_e`$ in the two shocks. Thus, for example, $`\gamma _{em,r}=ϵ_{e,r}(m_p/m_e)(\mathrm{\Gamma }_i/\mathrm{\Gamma }_{\mathrm{BM}})(p2)/(p1)`$ for the reverse shock, and $`\gamma _{em,f}=ϵ_{e,f}(m_p/m_e)\mathrm{\Gamma }_{\mathrm{BM}}(p2)/(p1)`$ for the forward shock (for $`p>2`$).

Synchrotron emission from electrons accelerated in the reverse shock typically peaks in the optical/IR band (Mészáros & Rees, 1997; Sari & Piran, 1999),

$$\epsilon _{\gamma ,\mathrm{RS}}^{ob.}\mathrm{}\mathrm{\Gamma }\gamma _{em,r}^2\frac{3}{2}\frac{qB}{m_ec}=0.2n_1^{1/2}\mathrm{\Gamma }_{i,2.5}^2ϵ_{e,r,1}^2ϵ_{B,r,2}^{1/2}\mathrm{eV}.$$
(3)
Here, $`\mathrm{\Gamma }_i=10^{2.5}\mathrm{\Gamma }_{i,2.5}`$, $`ϵ_{e,r}=10^1ϵ_{e,r,1}`$, and $`ϵ_{B,r}=10^2ϵ_{B,r,2}`$. This energy is similar to the photon energy $`\epsilon _{\mathrm{ssa}}`$ at which the optical depth for synchrotron self-absorption in the reverse shock equals unity,

$$\epsilon _{\mathrm{ssa}}^{ob.}0.8E_{53}^{1/3}n_1^{1/3}T_2^{2/3}\mathrm{\Gamma }_{i,2.5}^{1/3}ϵ_{e,r,1}^{2/3}ϵ_{B,r,2}^{1/3}\mathrm{eV}.$$
(4)
Denoting by $`\gamma _{c,r}`$ the Lorentz factor of reverse shock electrons for which the synchrotron cooling time is comparable to the ejecta expansion time, and assuming that electrons with $`\gamma >\gamma _{c,r}`$ emit nearly 100% of their energy during the dynamical time, the photon energy density is (for $`p=3`$)

$$u_{ph,syn,r}2\frac{\gamma _{em,r}}{\gamma _{c,r}}ϵ_{e,r}u=2.5\times 10^2E_{53}^{1/2}n_1^{3/2}T_2^{1/2}\mathrm{\Gamma }_{i,2.5}ϵ_{e,r,1}^2ϵ_{B,r,2}\mathrm{erg}\mathrm{cm}^3,$$
(5)
where the thermal energy density is $`u=4\mathrm{\Gamma }_{\mathrm{BM}}^2nm_pc^2`$.

## 3 Constraints on model parameters

Electrons accelerated in the forward shock, with Lorentz factors $`\gamma _{em,f}ϵ_{e,f}(m_p/m_e)\mathrm{\Gamma }_{\mathrm{BM}}/210^{4.5}`$, inverse-Compton scatter the optical photons, boosting their energy by $`\gamma _{em,f}^2`$ to $`1`$ GeV. The spectral shape of inverse-Compton scattered photons at observed energy range, $`\epsilon _\gamma <1`$ GeV, is similar to the spectral shape of the reverse shock synchrotron spectrum at $`\epsilon _\gamma <\mathrm{\hspace{0.33em}1}`$ eV, as the energy of all photons is boosted by roughly the same factor $`\gamma _{em,f}^2`$. However, the spectral shape is expected to smoothen somewhat, because Compton scattering is not exactly monochromatic. The low energy synchrotron spectrum is affected by self-absorption (Eq. 4). The optical depth increases as photon energy decreases, leading to a thermal spectrum $`dn_\gamma /d\epsilon _\gamma \epsilon _\gamma ^1`$ at $`\epsilon _\gamma \epsilon _{\mathrm{ssa}}`$. This spectrum is harder than observed. However, it becomes softer as $`\epsilon _\gamma `$ approaches $`\epsilon _{\mathrm{ssa}}\epsilon _{\gamma ,\mathrm{RS}}`$, which may allow to produce a spectrum consistent with observations.

Thus, in order to explain the observed spectrum as due to inverse-Compton scattering of self-absorbed synchrotron photons, the following conditions must be met:

1. The energy of synchrotron photons for which the reverse shock optical depth is unity, $`\epsilon _\gamma =\epsilon _{\mathrm{ssa}}`$, should be increased via inverse-Compton scattering by the lowest energy electrons in the forward shock, $`\gamma _e=\gamma _{em,f}`$, to $`200\mathrm{MeV}`$, i.e.

$$\gamma _{em,f}^2\epsilon _{\mathrm{ssa}}200\mathrm{MeV}.$$
(6)
2. The lowest energy electrons in the forward shock, $`\gamma _e=\gamma _{em,f}`$, should lose energy (via synchrotron and inverse-Compton scattering) on a time scale longer then $`T`$. If this condition is not met, cooling of the electrons would lead to a spectrum $`\mathrm{d}\mathrm{ln}n_e/\mathrm{d}\mathrm{ln}\gamma _e2`$ at $`\gamma _e<\gamma _{em,f}`$, which will produce an inverse-Compton spectrum $`dn_\gamma /d\epsilon _\gamma \epsilon _\gamma ^{3/2}`$ below 1 GeV, softer than observed. This condition can be written in the form

$$\begin{array}{c}\gamma _{c,f}\gamma _{em,f},\gamma _{IC,c,f}\gamma _{em,f},\hfill \end{array}$$
(7)
where $`\gamma _{IC,c,f}`$ ($`\gamma _{c,f}`$)is the Lorentz factor of forward shocked electrons for which the cooling time due to inverse-Compton (synchrotron) emission is comparable to the expansion time. The inverse-Compton cooling time is estimated via $`t_{cool,IC}=t_{cool,syn}(u_{B,f}/u_{ph,syn,r})`$, where $`u_{B,f}=ϵ_{B,f}u`$ is the energy density in the magnetic field at the forward shock.
3. Synchrotron emission from the forward shock electrons, as well as inverse-Compton emission from reverse shock electrons, should not exceed the observed flux at the 10 KeV to 1 MeV range. These requirements can be met by requiring either (a) The flux to be lower than observed, or (b) The photons to be emitted below 10 KeV. For synchrotron emission from the forward shock this implies either

$$\frac{L_{IC,f}}{L_{syn,f}}=\frac{u_{ph,syn,r}}{u_{B,f}}>10^3,\mathrm{or}\epsilon _{\gamma ,\mathrm{FS}}^{ob.}1\mathrm{KeV}.$$
(8)
For inverse-Compton emission from the reverse shock it implies either

$$\frac{L_{IC,f}}{L_{IC,r}}=\frac{ϵ_{B,r}}{ϵ_{e,r}}\frac{\gamma _{c,r}}{\gamma _{em,r}}\gamma _{em,f}^2\tau _{T,f}>10^3,\mathrm{or}\epsilon _{\gamma ,IC,\mathrm{RS}}^{ob.}=\epsilon _{\gamma ,\mathrm{RS}}^{ob.}\gamma _{em,r}^21\mathrm{KeV}.$$
(9)

A lengthy, but straight forward, manipulation of the above constraints lead to the following constraints on model parameters111 A detailed description of the derivation may be found at http://www.weizmann.ac.il/$``$asaf/ :

$$\begin{array}{cccc}a.\hfill & \hfill E_{53}& & 4T_2^{3/2},\hfill \\ b.\hfill & \hfill \left(\frac{E_{53}}{n_1}\right)& & 50T_2^4,\hfill \\ c.\hfill & \hfill ϵ_{e,r,1}& & 1.5E_{53}^1T_2^{5/2},\hfill \\ d.\hfill & \hfill ϵ_{B,r,2}& & 15T_2^{37/8},\hfill \\ e.\hfill & \hfill ϵ_{e,f}& & 0.33,\hfill \\ f.\hfill & \hfill ϵ_{B,f,2}& & (2.5\times 10^4,5\times 10^4T_2^{3/4},5\times 10^3T_2^{1/2},1.2\times 10^2T_2^{5/4}),\hfill \\ g.\hfill & \hfill \mathrm{\Gamma }_{i,2.5}& & 0.76T_2^{1/8}.\hfill \end{array}$$
(10)
It is sufficient to meet 1 of the 4 conditions (f). It should be pointed out here that the redshift of the source of GRB941017 is considered a free parameter in the above analysis. A particular choice of model parameters implies a particular choice of redshift given the measured flux, as illustrated in few examples in figure 1.

## 4 Numerical results

In order to obtain an accurate description of the spectrum, we have carried detailed numerical calculations of the emission of radiation from a fireball in transition to self-similar behavior. Our numerical model contains full description of cyclo-synchrotron emission, synchrotron self absorption, inverse- as well as direct Compton scattering, pair production and annihilation and cascade processes occurring at high energies. Full description of the model appears in Pe’er & Waxman (2004). In the present version, photons emitted after acceleration of electrons at the reverse shock, serve as seed photons for physical processes occurring after acceleration of electrons in the forward shock.

Figure 1 presents three examples of the results of our calculations for three different sets of parameter values, consistent with the constraints derived in § 3. It demonstrates that the process described above may indeed account for the observed hard spectrum of the high energy component.

An alternative explanation of the spectrum, emission from internal collisions within the expanding fireball, is also presented. For this scenario, we assumed $`ϵ_B`$ and $`ϵ_e`$ to be similar in the two shock waves created as two plasma shells collide. The observed spectrum may be reproduced adopting rapid variability, $`\mathrm{\Delta }t_{var}10^5`$ s, and very low value of the equipartition fraction, $`ϵ_B10^7`$. Here too, the hard high energy tail is due to inverse-Compton scattering of self-absorbed synchrotron spectrum. Accounting for $`1`$ MeV GRB emission, like the variable low energy component observed at the first $`100`$ s in GRB941017, typically requires in the internal shock scenario equipartition fraction close to unity. This, and the temporal characteristics of the high energy component, suggest that the high energy tail is produced during the transition to self-similar expansion, rather than in internal shocks.

The following point should be noted here. González *et al.* (2003) present spectra averaged over 5 distinct time intervals. The flux carried by the high energy tail shows variations between different time intervals, while large uncertainties in the determination of the spectral index in individual time intervals make it difficult to determine whether the spectral index is time dependent. In principle, time dependence of the high energy flux and spectrum provides additional constraints on the model. Such variations would arise in the model due to the expansion of the plasma during the transition phase, and due to deviations from homogeneity in the fireball shell as well as in the surrounding medium. Thus, the observed time dependence may allow to constrain the distributions of density and Lorentz factor within the fireball ejecta, and the distribution of ambient medium density. Given the relatively large uncertainties in the determination of the time dependence of the flux and of the spectral index, we find that such analysis is not warranted by the current data.

## 5 Discussion

We have shown that the high energy spectral component of GRB941017 may be explained as inverse-Compton emission from electrons accelerated by the shock driven by the fireball into its surrounding medium during the transition to self-similar expansion (Figure 1). The high energy tail may also be explained as inverse-Compton emission from electrons accelerated in internal shocks (Figure 1), and in the reverse shock during transition to self-similar expansion (Granot & Guetta, 2003). As explained in the introduction, we consider the latter explanations less likely, since they require modifications of the fireball physics (in particular, magnetic field well below equipartition in the fireball plasma), and since the temporal behavior is hard to account for in the internal shock scenario.

We have shown that the high energy data, when interpreted as emission from forward shock electrons, allow to put stringent constraints (Eq. 10) on fireball parameters, such as the total (isotropic equivalent) energy $`E`$, initial Lorenz factor $`\mathrm{\Gamma }_i`$ and ambient medium density $`n`$. Such constraints can not be inferred from keV—MeV observations alone. The inferred values of $`E`$ and $`\mathrm{\Gamma }_i`$ are similar to those typical to cosmological GRBs, while the inferred constraint $`E/n10^{56}\mathrm{erg}\mathrm{cm}^3`$ implies a ratio $`E/n`$ which is higher than typically obtained, $`E/n10^{54}\mathrm{erg}\mathrm{cm}^3`$ (Bloom, Frail & Kulkarni, 2003). The inferred values of $`ϵ_{e,f}`$, $`ϵ_{B,r}`$ and $`ϵ_{e,r}`$ are consistent with those inferred from other GRB and afterglow observations. The value of $`ϵ_{B,f}`$ is usually less well constrained by observations, and is required to be well below equipartition in our case. This, and the large ratio of $`E/n`$ may account for the rareness of GRB941017-type high energy tails.

The unknown redshift of the source of GRB941017 is treated in our analysis as a free parameter: A choice of model parameters implies a choice of redshift (see figure 1). Values of $`E`$ close to the lower limit, $`10^{54}`$ erg, imply low redshift, $`z0.1`$. This is a direct consequence of the fact that the gamma-ray fluence of GRB941017 is unusually high, rather than of the fact that its high energy spectrum is unusually hard. The fluence, $`6.5\times 10^4\mathrm{erg}/\mathrm{cm}^2`$, implies an (isotropic equivalent) gamma-ray energy release of $`10^{55}`$ erg for $`z=1`$.

Figure 1 demonstrates that different scenarios accounting for the high energy component of GRB941017, as well as different model parameters in a given scenario, lead to different model predictions for the fluxes at optical, X-ray and sub-TeV energy bands. The predicted fluxes are well within the detection capabilities of the SWIFT (in optical) and GLAST (at 10–100GeV) satellites, and of sub-TeV ground based Cerenkov telescopes (e.g. HESS222http://www.mpi-hd.mpg.de/hfm/HESS/HESS.html , MAGIC 333http://hegra1.mppmu.mpg.de/MAGICWeb , MILAGRO (McEnery, 2003), VERITAS (Weekes et al., 2002)). Therefore, optical and sub-TeV observations on minute time scale will allow to put more stringent constraints on explosion parameters than those given in Eq. 10.

We thank J. Granot & D. Guetta for helpful discussions. This work was supported in part by a Minerva grant and by a grant from the Rosa and Emilio Segré fund.


text:
# Stability of vortex solitons in a photorefractive optical lattice

## Abstract

Stability of off-site vortex solitons in a photorefractive optical lattice is analyzed. It is shown that such solitons are linearly unstable in both the high and low intensity limits. In the high-intensity limit, the vortex looks like a familiar ring vortex, and it suffers oscillatory instabilities. In the low-intensity limit, the vortex suffers both oscillatory and Vakhitov-Kolokolov instabilities. However, in the moderate-intensity regime, the vortex becomes stable if the lattice intensity or the applied voltage is above a certain threshold value. Stability regions of vortices are also determined at typical experimental parameters. OCIS codes: 190.0190, 190.5330.

Vortex solitons are ubiquitous in many branches of physics such as optics  and Bose-Einstein condensation . In a homogeneous medium, bright vortex rings are unstable , and only dark vortex solitons are possible with defocusing nonlinearity . However, in the presence of a periodic optical lattice, stable lattice vortices become possible due to the guiding properties of the lattice. Indeed, recent theoretical work  has shown that in an optical lattice with Kerr nonlinearity, both on-site vortices (vortices whose singularity is located on a lattice site)  and off-site vortices (vortices whose singularity is located between sites)  are stable within certain ranges of parameters. These theoretical studies are quickly followed by experiments in photorefractive crystals, where vortex lattice solitons were observed very recently . For a review of other nonlinear localized states in one- and two-dimensional periodic optical waveguides, see .

Stability of vortex lattice solitons in photorefractive crystals is clearly an important issue. This question was considered in , where the evolution of a particular on-site lattice vortex under random-noise perturbations was simulated. It was found that the on-site vortex was stable to very long distances. However, we know that lattice vortices in photorefractive crystals can not be all stable. For instance, when the peak intensity (or power) of the vortex is high, the lattice is effectively weak, thus the lattice vortex would become the familiar ring vortex, which is known to be unstable \[see Fig. 1(b) below\] . The natural questions to ask then are: what lattice vortices are stable? If lattice vortices are unstable, what are the sources of their instability? So far, these questions have not been addressed comprehensively for either of the on-site and off-site lattice vortices.

In this paper, we study the off-site vortex solitons in a photorefractive optical lattice. Off-site lattice vortices are more closely packed — the diagonal distance between their four main lobes is $`\sqrt{2}`$ times shorter than that of on-site vortices. Thus their dynamics is stronger and more interesting. We show that these vortices are not only unstable in the high-intensity limit, but also in the low-intensity limit. However, they do become stable in the moderate-intensity regime if the lattice intensity or the applied voltage reaches over a certain threshold. We also determine the stability regions of vortices for a wide range of experimental parameters, and show that the stability region expands when the applied voltage increases.

The mathematical model for light propagation in a photorefractive crystal has been known for some time . Here we make the usual paraxial assumption, and the assumption that the photorefractive screening nonlinearity acts isotropically along the two transverse directions, both of which are justified in many experiments. If the probe beam is extra-ordinarily polarized, while the lattice is ordinarily polarized, then the probe beam does not affect the linear lattice. In this case, the governing equation for the probe beam is

$$iU_z+\frac{1}{2k_1}\left(U_{xx}+U_{yy}\right)\frac{1}{2}k_0n_e^3r_{33}E_{sc}U=0,$$
(1)
where $`U`$ is the slowly-varying amplitude of the probe beam, $`z`$ is the distance along the direction of the crystal, $`(x,y)`$ are distances along the transverse directions, $`k_0=\frac{2\pi }{\lambda _0}`$ is the wavenumber of the laser in the vacuum ($`\lambda _0`$ is the wavelength), $`n_e`$ is the refractive index along the extraordinary axis, $`k_1=k_0n_e`$, $`r_{33}`$ is the electro-optic coefficient for the extraordinary polarization, $`E_{sc}`$ is the space-charge field, $`E_{sc}=E_0/[1+I_l(x,y)+|U|^2]`$, $`E_0`$ is the applied DC field, and $`I_l`$ is the field intensity of the optical lattice. Here the intensities of the probe beam and the lattice have been normalized with respect to the dark irradiance of the crystal $`I_d`$. Material damping of the probe beam is very weak in typical experiments since the crystals are fairly short (up to 2 cm). Hence it is neglected in Eq. (1). If the lattice is periodic along the $`x`$ and $`y`$ directions (rectangular lattice), then $`I_l(x,y)=I_0\mathrm{sin}^2\left(\pi x/D\right)\mathrm{sin}^2\left(\pi y/D\right),`$ where $`I_0`$ is its peak intensity, and $`D`$ is its spacing.

Eq. (1) can be non-dimensionalized. If we measure the transverse directions $`(x,y)`$ in units of $`D/\pi `$, the $`z`$ direction in units of $`2k_1D^2/\pi ^2`$, and the applied voltage $`E_0`$ in units of $`\pi ^2/(k_0^2n_e^4D^2r_{33})`$, then Eq. (1) becomes

$$iU_z+U_{xx}+U_{yy}\frac{E_0}{1+I_0\mathrm{sin}^2x\mathrm{sin}^2y+|U|^2}U=0.$$
(2)
Consistent with the experiments , we choose physical parameters as $`D=20\mu m,\lambda _0=0.5\mu m,n_e=2.3,r_{33}=280pm/V.`$ Thus, in this paper, one $`x`$ or $`y`$ unit corresponds to $`6.4\mu m`$, one $`z`$ unit corresponds to 2.3 mm, and one $`E_0`$ unit corresponds to 20 V/mm in physical units.

Lattice vortices of Eq. (2) are sought of in the form $`U=u(x,y)e^{\mu z}`$, where $`\mu `$ is the propagation constant. We determined these vortices by a Fourier iteration method . At lattice intensity $`I_0=I_d`$ and applied voltage $`E_0=8`$, these vortices are shown in Fig. 1. We see that when the vortex’ peak intensity $`I_p`$ is high, the vortex becomes a familiar ring vortex \[see Fig. 1(b)\] since the optical lattice is relatively negligible in this case. As $`I_p`$ decreases, the vortex develops four major lobes at four adjacent lattice sites in a square configuration \[see Fig. 1(c,e)\]. When $`I_p`$ is low, the vortex spreads over to more lattice sites and becomes less localized \[see Fig. 1(f)\]. The phase fields of all these lattice vortices, however, remain qualitatively the same as in a regular ring vortex \[see Fig. 1(d)\]. Note that the singularities (centers) of these vortices are not on a lattice site, thus these vortices are off-site vortices. An interesting fact we found is that, for given lattice intensity and applied voltage values, lattice vortices with $`I_p`$ below a certain threshold $`I_{p,c}`$ do not exist. In the present case where $`I_0=I_d`$ and $`E_0=8`$, this threshold value is $`I_{p,c}0.28I_d`$. This fact indicates that, unlike fundamental lattice solitons, lattice vortices do not bifurcate from infinitesimal Block waves.

We can further determine the power and peak-intensity diagrams versus the propagation constant $`\mu `$. Here the power is defined as $`P_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}|u|^2𝑑x𝑑y`$. When $`I_0=I_d`$ and $`E_0=8`$, the results are shown in Fig. 2(a). We see that the peak-intensity is a monotone-decreasing function of $`\mu `$, but the power is monotone-decreasing only when the peak intensity is above $`0.34I_d`$. Below this intensity value, the power starts to increase with $`\mu `$. This behavior qualitatively holds also for other $`I_0`$ and $`E_0`$ values, and it is similar to that in the Kerr medium .

Now we address the critical question of linear stability of these vortices in a photorefractive lattice. High-intensity lattice vortices clearly should be linearly unstable because they approach the regular ring vortex \[see Fig. 1(b)\] . The instability is oscillatory (i.e., the unstable eigenvalues are complex). At low intensities, $`dP/d\mu >0`$, hence the lattice vortices are expected to be linearly unstable as well according to the Vakhitov-Kolokolov (VK) criterion . The VK instability is purely exponential (i.e., the unstable eigenvalues are purely real). How about the stability behaviors of vortices at moderate peak intensities? To answer this question, we have simulated the linearized equation of (2) around lattice vortices $`u(x,y)e^{\mu z}`$ for very long distances, and obtained the unstable eigenvalues $`\sigma `$ of small disturbances (the real part of $`\sigma `$ is the growth rate). The results for $`I_0=I_d`$ and $`E_0=8`$ are shown in Fig. 2(b). We find that lattice vortices are linearly unstable when $`I_p>2.1I_d`$ and $`I_p<0.70I_d`$, consistent with our expectations. In addition, the instability for $`I_p>2.1I_d`$ is oscillatory \[Im$`(\sigma )0`$\], and the VK instability for $`I_p<0.34I_d`$ is purely exponential \[Im$`(\sigma )=0`$\], as we would expect. However, Fig. 2(b) reveals another oscillatory instability for $`0.34I_d<I_p<0.70I_d`$, which was not anticipated. This additional oscillatory instability has been seen in the Kerr medium before .

A more important result revealed by Fig. 2(b) is that for $`0.70I_d<I_p<2.1I_d`$, lattice vortices are linearly stable. This is an important result. It implies that lattice vortices with such moderate intensities could be observable in experiments.

When lattice vortices are linearly unstable, what is the outcome of the instability? To address this question, we select the linearly-unstable vortex soliton with $`I_0=1,E_0=8`$ and peak intensity $`I_p=3I_d`$, and perturb it by random noise. The noise has Gaussian intensity distribution in the spectral $`k`$-space with FWHM 2 times larger than the soliton FWHM spectrum. The noise power is 1% of the soliton’s. The simulation result on the evolution of this vortex under noise perturbations is shown in Fig. 3. We see that this vortex breaks up into a fundamental lattice soliton plus some radiation. This breakup scenario is typical of unstable lattice vortices under noise perturbations.

When lattice vortices are linearly stable, how do they evolve nonlinearly? To answer this question, we select the linearly-stable vortex soliton with $`I_0=1,E_0=8`$ and $`I_p=1.5I_d`$, and perturb it by the same random noise as described above. The simulation result on the evolution of this perturbed vortex is shown in Fig. 4. We see that this vortex does propagate stably. In addition, its phase structure is maintained throughout the evolution. Evolution of other linearly-stable lattice vortices under weak perturbations is similar. This means that linearly-stable vortex solitons could be observed in experiments, as the work  has shown.

Above at specific lattice intensity and applied voltage values $`I_0=I_d`$ and $`E_0=8`$, we have revealed the sources of instability of lattice vortices, and obtained stable lattice vortices. The next question quickly follows: if the lattice intensity and voltage values are varied, how would they affect the stability properties of lattice vortices? This question is important for experiments. To find the answer to this question, we have systematically determined the linear stabilities of lattice vortices at various lattice intensity, applied voltage and vortex peak-intensity values. The results are summarized in Fig. 5. Here at two applied voltage values $`E_0=8`$ and 10, the stability boundaries are presented in the $`(I_p,I_0)`$ plane. This figure reveals several important facts. First, high-intensity and low-intensity vortex solitons are always linearly unstable, as we have observed in Fig. 2 above. Second, when the applied-voltage value $`E_0`$ is fixed, there is a threshold lattice intensity $`I_{0,c}`$, below which all lattice vortices (including moderate-intensity ones) are linearly unstable. When $`E_0=8`$, this threshold value is $`I_{0,c}0.7I_d`$; while when $`E_0`$ is increased to 10, $`I_{0,c}`$ decreases to $`0.44I_d`$. Similarly, when the lattice intensity $`I_0`$ is fixed, there is also a threshold applied-voltage value below which all lattice vortices are linearly unstable. Thirdly, when the applied voltage increases, the region of stable lattice vortices expands. In other words, higher-applied voltage stabilizes lattice vortices. Fig. 5 should be helpful to experimentalists on their choices of physical parameters for the observation of lattice vortices.

In summary, we have carried out a stability analysis on off-site lattice vortices in photorefractive optical lattices. We showed that high- and low-intensity lattice vortices suffer oscillatory and VK instabilities, but moderate-intensity vortices can be stable when the applied voltage or lattice intensity is above a certain threshold. Higher applied voltage stabilizes lattice vortices.


text:
# SNUTP-03-025 Description and comparison of (Fat7)̄ and HYP fat links

## 1 INTRODUCTION

Despite its prominent role as a useful tool in the non-perturbative study of QCD, lattice field theory has long struggled with the difficulties of accurately representing a continuum field theory in a discrete spacetime. The occurrence of discretization errors has lead to the development of a great many improvement schemes, some more useful than others. One extremely popular method of improvement has been the use of “fat links”, in which each ordinary or “thin” link of the lattice, $`U_\mu (x)`$, is replaced by a linear combinaton of the thin link and the links adjacent to it. Two particular fattening algorithms will be of interest to us in this report, HYP or hypercubic blocking , and $`\overline{\mathrm{Fat7}}`$ blocking (that is, SU(3)-projected Fat7 blocking ). These two algorithms have been shown to be perturbatively equivalent at one-loop level . We will present the results of a preliminary numerical investigation to determine whether there are any signifiant non-perturbative differences between the two algorithms. These results will point the way to more detailed investigations.

## 2 HYPERCUBIC BLOCKING

The HYP algorithm was proposed as a modification of standard APE smearing . APE smearing involves constructing a fattened link $`V_\mu (x)`$ by adding the sum of staples, weighted by some factor $`\alpha `$, to the thin link,

$$V_\mu (x)=(1\alpha )U_\mu (x)+\frac{\alpha }{6}\underset{\nu \mu }{}W_\nu (x)$$
(1)
followed by projection back to SU(3) (where $`W_\nu (x)`$ refers to the sum of three-link staples in the positive and negative $`\nu `$ directions, with central link parallel to $`U_\mu (x)`$). Clearly $`n`$ repeated applications of such an algorithm will access the gauge fields at a distance $`n`$ lattice spacings from the original thin link (although smoothing algorithms such as cooling and smearing can be viewed as diffusive “random-walk” processes which primarily affect physics only within a range $`a\sqrt{n}`$, as investigated in ). To localize APE smearing within the smallest possible volume, HYP blocking was defined as being equivalent to three iterations of APE smearing with the caveat that the staples at each stage could not be constructed from links which are fattened in the same plane as the staple itself. The original thin link is therefore fattened by pairs of positively-and-negatively-oriented staples in three directions. Each of these staples is constructed from links fattened by staples in only two (positive and negative) directions, and each of these staples is constructed from links fattened in a single plane (for further details please refer to ). HYP blocking therefore uses only the staples defined by the boundaries of hypercubes attached to the original thin link (Fig 1). Projection back to SU(3) is performed after each level of fattening.

We shall henceforth refer to the links fattened by the outermost set of staples as the initial-stage fat links, the links fattened by staples formed from the initial-stage links as the middle-stage fat links, and the completely fattened links (i.e. those fattened by staples formed from middle-stage fat links) as the final-stage links.

Since its introduction, HYP blocking has shown itself to be extremely effective at reducing taste-symmetry breaking effects, and has gained widespread attention and use.

## 3 FAT7 AND $`\overline{\mathrm{FAT7}}`$ BLOCKING

Fat7 blocking  is a form of fattening which incorporates not only standard three-link staples but also five-link and seven-link staples (Fig. 2). Clearly, with staples up to length seven it is possible to traverse all three directions orthogonal to the original thin link, as is the case with HYP blocking. However Fat7 blocking as it was originally envisaged did not incorporate SU(3) projection. $`\overline{\mathrm{Fat7}}`$ is a modification of this algorithm which incorporates SU(3) projection of the final sum of link and staples .

If we construct $`\overline{\mathrm{Fat7}}`$ fat-links in an analogous manner to the iterative construction of HYP links we have greater freedom in our application of SU(3)-projection. An iterative approach allows us to include SU(3) projection at each level of the construction of the link, rather than just at the final level. The only difference between HYP and $`\overline{\mathrm{Fat7}}`$ is then that in HYP blocking each link in each staple is fattened at the subsequent level, whereas in $`\overline{\mathrm{Fat7}}`$ it is only the links parallel to the original thin link that are fattened (see Fig. 3).

We can see that the second covariant derivative operator  may be used to recursively define a fattened link

$$L_\nu (\alpha )U_\mu (x)=(12\alpha )U_\mu (x)+\alpha W_\nu (x)$$
(2)
($`W_\nu (x)`$ defined as in Eq. (1)). With the parameter $`\alpha `$ taking the value one-quarter this operator can be interpreted as suppressing flavour-changing gluon interactions by vanishing in the limit as gluon momentum approaches the lattice cut-off. The Fat7 link can easily be constructed by repeated application of this operator as

$$V_\mu =\frac{1}{6}\underset{\mathrm{perm}(\nu ,\rho ,\lambda )}{}L_\nu (\alpha )\left(L_\rho (\alpha )\left(L_\lambda (\alpha )U_\mu \right)\right)$$
(3)
where we note explicitly that the sum over permutations of the directions is taken. This repeated use of the same operator accounts for the use of a single weighting parameter $`\alpha `$ at each stage, whereas HYP blocking uses three parameters $`\alpha _1`$, $`\alpha _2`$, $`\alpha _3`$.

As noted above, to construct the $`\overline{\mathrm{Fat7}}`$ links from the Fat7 links, SU(3) projection must be performed after the final summation is taken, but it can also optionally be applied after each or either of the first two fattening stages, leaving us with four possible combinations of projection schemes to define our $`\overline{\mathrm{Fat7}}`$ link (i.e. projection at the final, initial+final, middle+final, and initial+middle+final stages).

## 4 COMPARISON OF FATTENING APPROACHES

Fat links have shown themselves to be extremely useful at overcoming problems associated with the lattice formulation. Specifically they have shown a reduction in the severity of taste-symmetry breaking associated with staggered fermions, and a reduction in the severity of the exceptional configuration problem associated with the Wilson fermion formulation. Both of these problems can be directly related to short-scale fluctuations of the link values. In , the values of the minimum plaquettes calculated on an ensemble of configurations were used as a probe of the most severe link fluctuations. The maximum plaquette values are bounded above, and hence we expect an increase in the minimum plaquette value of any configuration to indicate a smaller range of link fluctuations across that configuration, and hence a reduction in the scale of the largest short-range link fluctuations.

The process of fattening tends to increase the minimum plaquette values of different configurations to differing degrees, spreading the minimum values over quite a wide range (as can be seen by comparing Figs. 4 and Figs. 7 and 7). We expect that certain measurements will achieve a lower level of statistical uncertainty if performed on an ensemble of configurations having a narrow spread of minimum plaquette values. One of the aims of this preliminary study will be to identify whether the $`\overline{\mathrm{Fat7}}`$ blocking algorithm is capable of producing a narrower spread of minimum plaquette values than HYP, and is hence worthy of a more detailed investigation.

## 5 RESULTS

We present the minimum plaquette values for an ensemble of 194 quenched configurations generated on an $`8^3\times 32`$ lattice at $`\beta =5.7`$. We also present the minimum plaquette values for the ensembles obtained after blocking each configuration with the HYP algorithm (Figs. 7 and 7), and with the $`\overline{\mathrm{Fat7}}`$ algorithm (Figs. 7 to 10). As noted above the natural choice for the $`\overline{\mathrm{Fat7}}`$ weighting parameter is $`\alpha =0.25`$. For HYP we are interested in two choices of parameters, those found non-perturbatively in , namely $`\alpha _1=0.75`$, $`\alpha _2=0.6`$, $`\alpha _3=0.3`$ (HYP-I), and those found perturbatively in , namely $`\alpha _1=7/8`$, $`\alpha _2=4/7`$,

$`\alpha _3=1/4`$ (HYP-II), although we note the caveat that choices of $`\alpha _1`$ larger than 0.75 may tend to destabilize the smearing algorithm, making the blocked configuration rougher rather than smoother . For the $`\overline{\mathrm{Fat7}}`$-blocked configurations we have analysed the four different SU(3) projection options described in section 3. The histograms of minimum plaquette values are presented in Figs. 4 to 10. All figures are to the same scale.

We can see clearly that the initial unblocked data produce a very narrow peak. After HYP blocking the distribution of minimum plaquette values has been shifted towards (but not into) the positive region, and taken on a Poisson distribution-like outline. The data obtained from $`\overline{\mathrm{Fat7}}`$ with SU(3) projection at the final stage, and at all three stages are especially interesting from an analytical basis . It is clear that all four $`\overline{\mathrm{Fat7}}`$ fattening schemes lead to minimum plaquette distributions that have a substantially narrower spread of values than the HYP-I and HYP-II data. Although the number of configurations used in this study is small, the results seem significant enough to warrant further investigations with $`\overline{\mathrm{Fat7}}`$.

## 6 CONCLUSIONS

Although this work is only preliminary, we find encouraging signs that $`\overline{\mathrm{Fat7}}`$ may be capable of producing results which are as good as HYP, and possibly better for some calculations. Both HYP and $`\overline{\mathrm{Fat7}}`$ algorithms can be made substantially faster and more efficient by pre-calculating and storing the staples across the entire lattice, and using these to construct the links at the next level of fattening. This approach avoids calculating the same staple more than once, as one would be forced to do if one fattened in a naive site-by-site manner. Using this approach $`\overline{\mathrm{Fat7}}`$ can be made less memory-intensive than HYP, as it only requires the staples in a limited number of directions to be pre-calculated at each stage of fattening. An efficiently-designed algorithm for $`\overline{\mathrm{Fat7}}`$ does not sacrifice speed, and enables us to use only one-quarter of the memory that is needed to store the pre-calculated staples for HYP. In future work we hope to determine whether this cheaper form of fattening merely equals, or in fact can produce results which are superior to HYP for certain calculations.

## 7 ACKNOWLEDGEMENTS

We wish to thank Anna Hasenfratz for discussions concerning the relation of this work to her own, and Derek Leinweber for bringing refence  to our attention.


text:
# Structure of densified amorphous germanium dioxide

## Abstract

Classical molecular dynamics simulations are used to study the structure of densified germanium dioxide ($`GeO_2`$). It is found that the coordination number of germanium changes with increasing density (pressure) while pressure released systems do not exhibit a marked change in local structure as compared to the virgin system. The structural modification with pressure appears to be stepwise and gradually affects long range (through the reduction of the long-range correlations as seen from the shift of the first sharp diffraction peak), intermediate range (by angular reduction) and finally short range structure (by tetrahedron distorsion).

Pacs: 61.43.Fs-46.30.Cn

Structural transitions in minerals are known to take place under various geological conditions. In the Earth’s interior, silicate and alumino-silicate melts change their local structure causing strong density and viscosity modifications in magmas and silicon exhibits at high pressure a change of its coordination number. These results have been mostly obtained with high pressure-temperature experiments but also with computer simulations reproducing extreme conditions. Both have inferred the nature of these structural transitions, the structure and the phase portrait of the liquid state. However, while a majority of studies have been devoted to the silica and silicate chemistry, little has been done to elucidate the corresponding behavior in germanium dioxide ($`GeO_2`$) even though this material is a structural analog of silica in many respects: both materials exhibit at ordinary conditions a tetrahedral local structure, they can also exist in $`\alpha `$ as well as $`\beta `$ quartz phases and the change of germanium coordination from four to six also occurs at high pressure. The structure of in situ densified or permanently densified $`GeO_2`$ remains however controversial. While the global increase with pressure of the germanium-oxygen distance in $`GeO_2`$ has been related from x-ray diffraction  with the conversion of tetrahedral Ge(4) into octahedral Ge(6), it seems that this structural change is reversible as no Ge(6) is found in decompressed samples, a situation which does not occur in the decompressed rutilelike $`cGeO_2`$ The interest in germanium coordination change has been also motivated by the observation of the so-called “germanate anomaly” which corresponds to a maximum in density and refractive index when $`1516mol\%`$ $`Na_2O`$ are added into the basic $`GeO_2`$ network former. These binary systems have been investigated by various spectroscopic tools - and it is suggested that the increasing presence of $`GeO_6`$ octahedra (Ge(6)) within the network is responsable of the anomaly. On the other hand, micro-Raman applied on the same systems suggests that the anomaly is due to a massive conversion of 4-membered rings into more close-packed ones such as 3-membered rings, with no, or at least very few, Ge(6) present.

A preliminary task if one wishes to describe density induced structural changes in germanates, is first to understand how the basic network former $`GeO_2`$ changes with densification. While crystalline phases of $`GeO_2`$ have been studied from numerical calculations, we are not aware of any published result on simulated liquid and amorphous $`GeO_2`$. It is therefore of striking interest to see what Molecular Dynamics (MD) can tell about a certain number of experimental open questions which remain at this stage. How does densification affect the local structure in the amorphous material ? Which thermodynamical quantity (temperature, Pressure) controls mostly the existence of a 6-fold coordinated germanium ? Does octahedral germanium exists in pressure released (permanently densified) $`GeO_2`$ ? This Letter attempts to address some of these basic issues by providing the first MD study of amorphous germania. The results show several main features. They demonstrate the reversible nature of pressure induced changes and no major difference with the virgin structure is found upon pressure release, while the local coordination number of the germanium increases smoothly from 4 to 5.5 with applied pressure, a result that would be accessible from in situ experiments as already realized for silica. Furthermore, several régimes of densification can be clearly identified: a first régime \[$`P1.8GPa`$\] during which no global change in the local structure is found to occur whereas the increase in the position of the first sharp diffraction peak (FSDP) suggests a global reduction of the longer range correlations. In the window $`1.8GPaP2.8GPa`$, a sharp decrease of the intertetrahedral bond angle permits to densify the structure with no distorsion of the basic $`GeO_{4/2}`$ tetrahedron. This leads to the buckling of the network connected tetrahedra. Finally, for larger pressures, distorsion of the tetrahadron sets in. As a consequence, additional constraints appear for $`P>3GPa`$ which produce a global stiffening of the network

The system consists of 256 germanium and 512 oxygen atoms interacting via a Born-Huggins-Mayer type potential which has been fitted in the case of $`GeO_2`$ by Oeffner and Elliot to recover the crystalline phases of $`GeO_2`$ and its vibrational spectra. The atoms have been first confined in a cubic box of length $`L=23.044\AA `$ in order to recover the experimental value of the density ($`\rho _g=3.66g.cm^3`$). After having thermalized the system at $`3000K`$ for 20000 time steps (20 ps) the system has been cooled to $`300K`$ with a linear cooling schedule at a quench rate of $`2.5\times 10^{12}K.s^1`$. Integration has been done using a leap-frog Verlet algorithm. Various configurations (positions and velocities) have been saved at different temperatures which have been used as starting configurations for production runs of $`10^5`$ steps.

Densification has been realized reproducing experimental conditions, i.e. starting from an initial configuration at $`T=300K`$ and $`\rho _g`$ and increasing the density during $`10^4`$ time steps. At the density $`\rho _g`$, the glass transition temperature was about $`1660K`$, and as expected, the $`T_g`$ shifts to the higher temperatures with increasing density ($`1894K`$ at $`\rho =1.1\rho _g`$). The sample data of Price and co-workers has been almost recovered, i.e. after MD pressure release from 15.16GPa determined from the virial. The final density of the decompressed system at zero pressure was $`\rho =4.5g.cm^3=1.25\rho _g`$. The corresponding configuration (termed in the following as “permanently densified system”) was taken for comparison with the experimental study of permanently densified $`GeO_2`$

The results for the partial structure factors $`S_{ij}(Q)`$ are displayed in Fig. 1 at $`\rho _g=3.66g.cm^3`$ which show a fair agreement with the experimental results obtained from a combination of x-ray data.

One has however to note that these measurements have been refined but without any complete resolution of all the partials. The principal features in the structure factors of $`GeO_2`$ are the peaks occuring at $`Q_P=1.52\AA ^1`$ (the first sharp diffraction peak, FSDP, corresponding to a correlation length $`L_1=2\pi /Q_P4.1\AA `$). Peaks at $`2.6\AA ^1`$ and $`4\AA ^1`$ appearing in the partial structure factors of the Ge-O and O-O can be associated with chemical short range order (CRSO) and topological short-range order as derived from scaling considerations. For the former, the similar height in $`g_{GeGe}`$ and $`g_{OO}`$ but with opposite sign suggests the CRSO nature of the network. The evolution of the total scattering function and the position of the FSDP with respect to pressure is discussed below.

The calculation of the pair correlation functions permits to extract the simulated bond lengths which also agree with experimental findings: $`d_{GeGe}=3.32\AA `$, $`d_{GeO}=1.72\AA `$ and $`d_{OO}=2.81\AA `$ (to be respectively compared to the experimental values of $`3.16\pm 0.03\AA `$, $`1.73\pm 0.03\AA `$ and $`2.83\pm 0.05\AA `$). Weak changes are found when comparing the virgin and the permanently densified systems. The most striking difference is found in the Ge-Ge correlator which differs even at low distances whereas the loss of peak heights in $`g_{OO}`$ produces a slight shift in the corresponding running coordination number $`n_{OO}`$. The number of oxygen neighbours around an germanium atom remains however unaffected by pressure release. More dramatic is the change in structure in the pressurized system at $`16.6GPa`$ (dots in the inserts of Fig. 2) which show substantial differences in the bond distances ($`d_{GeGe}=3.25\AA `$, $`d_{GeO}=1.75\AA `$ and $`d_{OO}=2.56\AA `$) and an increase of the number of oxygen neighbours in the vicinity of a Ge atom (about 5 in the insert representing $`n_{GeO}`$, see also Fig. 3). Evidence of supplementary atoms in the first shell surrounding a central $`GeO_{4/2}`$ unit is also provided by the increase of $`n_{GeGe}`$.

Using the present simulation, it is now possible to analyze in more depth the low pressure behaviour. Focus is made on the distorsion of the tetrahedron parameter defined by $`\delta =d_{XO}/d_{OO}`$ \[X=Si,Ge\] which provides a direct measure of the effect of the pressure on the local structure of the network. For an ideal tetrahedron, $`\delta =\sqrt{3/8}`$. The latter quantity is of central interest when studying the flexibility of the glass under pressure and the rigid unit modes in the context of pressure induced rigidity. In Fig. 3 (left axis) are represented the variation of the tetrahedron parameter $`\delta `$ with pressure for both silica and germania. At low pressure, $`\delta `$ remains almost constant suggesting that the tetrahedral environment is preserved, slightly higher however than the value of a perfect tetrahedron. For $`P2.8GPa`$, there is increasing distorsion of the tetrahedron in germania. The way of distorsion appears however to be radically different as compared to the silica system. In the former it is found a stepwise increase (most noticable from the jump of $`\delta `$ at around $`3GPa`$), in contrast with the more or less smooth increase of $`\delta `$ for the latter. In the insert of Fig. 3 are shown the variation of the mean bond angles O-Ge-O and Ge-O-Ge with respect to pressure. For both the virgin and the permanently densified system, the O-Ge-O bond angle peaks at around $`109^o`$. As pressure is increased, the intratetrahedral bond angle decreases since the Ge-O bond distances increases. The constant value of the O-Ge-O bond angle at low pressure correlates of course with the absence of distorsion of the basic tetrahedron. The main feature provided by the angular analysis comes from the variation of the intertetrahedral mean bond angle Ge-O-Ge with pressure, which exhibits a sharp drop at around $`1.8GPa`$, from $`158^o`$ to $`135^o`$ followed by a stabilization at around $`140^o`$. Densification first applies on the angles connecting the tetrahedra and preserves the latter. If pressure keeps increasing, the distorsion of the tetrahedra sets in, which in turn stabilizes more or less the Ge-O-Ge bond angle.

Further interpretation is provided from the variation of the position of the FSDP with respect to pressure. In Fig. 4 are represented the scattering functions for different applied pressures which show a global broadening of the peak at $`4\AA ^1`$ whereas the position of the FSDP at $`1.5\AA ^1`$ is shifted to the higher values in $`Q`$, already at very low pressure and even before the angular drop at $`P=1.8GPa`$. With increasing pressure, the calculated FSDP broadens and becomes less intense, as currently observed from X-ray or neutron diffraction. It is worth to note that the peak at $`2.5\AA ^1`$ obtained in the present simulation is only weakly observable in the experiments from neutron diffraction displayed in Fig. 4. However, the simulated double peak distribution between $`1.5`$ and $`2.5\AA ^1`$ has been observed by different authors (see also discussion in ). The evolution with pressure of the position $`Q_P`$ of the FSDP is represented in the insert of Fig. 4. Both experiments and simulation show that $`Q_P`$ already increases at low pressures and then stabilizes at around $`Q=1.7\AA ^1`$ ($`2\AA ^1`$ in the simulation). This suggests that intermediate range order is immediately affected by the densification and then remain unaffected with further densification.

In conclusion, we have shown that simulated $`GeO_2`$ under pressure shows several main features with applied pressure and pressure release: i) a global increase of the number of oxygen neighbours in the vicinity of a germanium atom, ii) a stepwise change in the local structure with applied pressure, made of a reduction of long range correlation (seen from the position of the FSDP), a sharp reduction of the intertetrahedral bond angle and then a progressive distorsion of the $`GeO_{4/2}`$ tetrahedron and iii) no noticable change in local structure between a virgin and a permanently densified system. This clearly draws the following picture: Pressure applies on different length scales. With increasing magnitude, densification is realized by a successive deformation of long range structure, intermediate (angular) and finally short range structure (tetrahedral).

It is a pleasure to acknowledge instructive discussions with Yves Guissani and Bertrand Guillot. The Laboratoire de Physique Théorique des Liquides is Unité Mixte de Recherche n. 7600 du CNRS.


text:
# Electronic transport and vibrational modes in the smallest molecular bridge: H2 in Pt nanocontacts

## Abstract

We present a state-of-the-art first-principles analysis of electronic transport in a Pt nanocontact in the presence of H2 which has been recently reported by Smit et al. in Nature 419, 906 (2002). Our results indicate that at the last stages of the breaking of the Pt nanocontact two basic forms of bridge involving H can appear. Our claim is, in contrast to Smit et al.’s, that the main conductance histogram peak at $`G2e^2/h`$ is not due to molecular H2, but to a complex Pt2H2 where the H2 molecule dissociates. A first-principles vibrational analysis that compares favorably with the experimental one also supports our claim.

Prediction ability is, possibly, the most important attribute associated with first-principles numerical implementations of quantum transport theory in molecular and atomic-scale systems. This desirable property relies on the correct handling of two closely-related factors: The atomic structure and the electronic structure of the nanoconstriction or molecular junction. The calculation of the electronic structure in current-carrying situations is a topic of increasing interest and has received due attention in recent yearsLang (1995); Yaliraki et al. (1999); Taylor et al. (2001); Palacios et al. (2001, 2002); Damle et al. (2001); Brandbyge et al. (2002). The computation of the atomic structure, which is well-established on standard geometry relaxation methods, is, however, much more demanding from a first-principles computational point of view and it is not subject to scrutinity on a regular basis in molecular electronics. Besides, there are inherent uncertanties associated with the various experimental procedures that can only be addressed by statistical analysis.

In this paper we show how a combination of state-of-the-art current-carrying first-principles electronic structure calculations and a careful structural analysis can shed light on a recent experimentSmit et al. (2002); Smit  of electronic transport in Pt nanocontacts where, apparently, H2 molecules anchor themselves to the metal and modify the standard conductance histogram of Pt. A H2 molecule is the simplest molecule and is ideal to carry out a detailed structural, vibrational, and transport analysis. On the basis of the extensive ab initio calculations we have performed, our main conclusion is the following: Contrary to the authors’ explanation for the conductance histograms where the main peak around $`𝒢𝒢_0`$ (for $`𝒢_0=2e^2/h`$) is attributed to a configuration where a H2 molecule bridges the gap between Pt electrodes (see also Ref. Heurich et al., ), we attribute this peak to a situation where the H2 molecule dissociates at the nanocontact and forms a stable complex Pt2H2 where the two H atoms chemisorb between electrode tip atoms \[see left inset in Fig. 2\]. Our proposed structure presents vibrational modes in accordance with the reported inelastic effectsSmit et al. (2002) and with a more recent vibrational analysisSmit . Finally, we speculate that the bridge with H in molecular form could account for the peak around $`𝒢0.1𝒢_0`$ reported in Ref. Smit et al., 2002, but not commented on there.

For the structural and vibrational analysis of the Pt nanocontact with H we have used the GAUSSIAN98 codeFrisch et al.  and for the transport calculations we have employed the Gaussian Embedded Cluster Method (GECM), recently developed by the authorsPalacios et al. (2001, 2002); Louis et al.  which, in turn, makes use of GAUSSIAN98. The GECM consists of a first-principles evaluation at the density functional theory (DFT) level of the density matrix and the Fock matrix $`F`$ (or selfconsistent hamiltonian) of the region determining the transport (here the H atoms and the Pt electrode atoms in proximity to them) in a current-carrying situation. As regards the DFT calculations for the core clusters (see insets in Figs. 1 and 3) we used the B3LYP functionalFrisch et al.  , the basis sets and core pseudotentials described in Ref. Pacios and Christiansen, 1985; Hurley et al., 1986; Ross et al., 1990 for Pt and the cc-pVTZ basis for HFrisch et al.  . Here we are concerned with the linear response conductance which is calculated in the standard mannerAgrait et al. :

$$𝒢=\frac{2e^2}{h}\mathrm{Tr}[t^{}t].$$
(1)
In this expression, Tr denotes the trace over all the orbitals in the cluster. The matrix $`t`$ is given by

$$t=\mathrm{\Gamma }_\mathrm{R}^{1/2}G^{(+)}\mathrm{\Gamma }_\mathrm{L}^{1/2}=\left[\mathrm{\Gamma }_\mathrm{L}^{1/2}G^{()}\mathrm{\Gamma }_\mathrm{R}^{1/2}\right]^{},$$
(2)
where the Green function matrix (in orthogonal basis) is

$$G^{(\pm )}=\left[(E\pm i\delta )1F\mathrm{\Sigma }^{(\pm )}\right]^1,$$
(3)
while $`1`$ is the unity matrix. The matrices $`\mathrm{\Gamma }_\mathrm{R}`$ and $`\mathrm{\Gamma }_\mathrm{L}`$ are given by $`i(\mathrm{\Sigma }_\mathrm{R}^{()}\mathrm{\Sigma }_\mathrm{R}^{(+)})`$ and $`i(\mathrm{\Sigma }_\mathrm{L}^{()}\mathrm{\Sigma }_\mathrm{L}^{(+)})`$, respectively, where the self-energy $`\mathrm{\Sigma }^{(\pm )}=\mathrm{\Sigma }_\mathrm{R}^{(\pm )}+\mathrm{\Sigma }_\mathrm{L}^{(\pm )}`$ represents the contribution of the semiinfinite left and right leads not included in the self-consistent evaluationPalacios et al. (2001, 2002). In order to single out the contribution of individual channels to the current we diagonalize the matrix $`t^{}t`$.

As a starting point we consider the conductance of a clean Pt nanocontact without H2. It would be desirable to simulate the whole breaking(formation) process of the nanocontact on stretching(indentation). This is, unfortunately, computationally not feasible from first principles for Pt. Out of the many possible atomic arrangements, we have chosen to study the conductance for the following nanocontact atomic configurations: 9-4-1-4-9 (single-atom contact) and 9-4-1-1-4-9 (chain contact) for the (001) direction and 6-3-1-3-6 (single-atom contact) and 6-3-1-1-3-6 (chain contact) for the (111) direction (the numbers indicate the number of atoms per plane). We have considered two possible plane orderings for the (111) direction: 6B-3A-$`\mathrm{}`$-3A-6B and 6C-3A-$`\mathrm{}`$-3C-6A (remember that the bulk one is $`\mathrm{}`$ A-B-C-A-B-C $`\mathrm{}`$). These nanocontacts correspond to cases with and without mirror symmetry with respect to the middle plane perpendicular to the stretching axis $`z`$ \[the (001) nanocontact has also been chosen with mirror symmetry\]. Specifically, we have taken bulk atomic distances for the atoms in the outer planes, but we have performed full relaxation in all the other coordinates, including the distance between the two outer planes. The resulting structures are thus representative of the ones responsible for the last(first) conductance plateau on stretching(indentation) cycles.

The conductance reaches $`2𝒢_0`$ at the Fermi energy (set to zero) for the two short-chain cases with mirror symmetry and drops slightly for the other one \[see Fig. 1(b)\]. Two conduction channels mainly contribute to the conductance in the three clusters considered with the two-atom chain \[see Fig. 1(d)\]. Given the polycrystalline nature of the electrodes, one should not expect mirror symmetry to be the dominant situation in the experiment so we can expect an average conductance around $`𝒢1.75𝒢_0`$. This value nicely coincides with the position of the lowest conductance histogram peak reported in Ref. Smit et al., 2002. The conductance of the single-atom contacts presents similar values for the (111) direction (although, as Fig. 1(c) shows, three channels now contribute to transport), but is appreciable higher in the (001) direction \[see Fig. 1(a)\]. Although we cannot afford a statistical analysis, the fact that the position of the lowest conductance histogram peak shifts to higher values when single-atom contacts are favoured (on only considering indentation cycles or higher bias voltages as reported in Ref. Nielsen et al., ) is consistent with our results. The conductance is also here noticeably dependent on the detailed atomic arrangement of the electrodes (note that the conductance increases appreciably on the scale of 1 eV below and above the Fermi energy). This is due to the fact that the Fermi level lies on the edge of the $`d`$ band for bulk and the contact atom has almost bulk coordination. This is also fairly consistent with the fact that Pt histogram peaks are broad in comparison to, e.g., those of AuAgrait et al. .

When the same conductance measurements are performed in a H2 atmosphere the histograms change noticeably. Inelastic vibrational spectroscopy along with theoretical support lead the authors in Ref. Smit et al., 2002 to conclude that a H2 molecule can position itself along the stretching axis of the Pt nanocontact right before it breaks apart. The authors claim that this configuration accounts for the main histogram peak for conductance where $`𝒢𝒢_0`$. With this claim in mind we consider nanocontacts like that shown in the inset of Fig. 1(b), but with a H2 molecule anchored in between electrode tip atoms. The results presented below are for the smallest cluster in the (001) direction, 4-1-H2-1-4, but we have performed similar calculations for (111) nanocontacts for reassurement. We maintain bulk interatomic distances for the four atoms in the two outer planes whose variable respective distance $`d`$ mimics the effect of the piezo voltage. For every $`d`$ we let the positioning of the H atoms and that of the two contacting Pt atoms relax, starting from a configuration where the H2 molecule sits along the stretching axis. Full blue circles in Fig. 2 represent the total energy vs. $`d`$ for all the stable configurations we have found in which the molecule stays on the $`z`$-axis (see right inset in Fig. 2). Beyond $`d9.5`$Å the molecule sticks to one of the electrodes or dissociates. Below $`d8.25`$Å the molecule tends to loose its alignment with the $`z`$-axis. (As a guide to the eye we also present with empty circles the energy for configurations where the molecule is forced to stay along the $`z`$-axis). In all these points the H2 retains its molecular form. At the middle point ($`d9`$ Å) the H-H bond distance and the H-Pt distance turn out to be similar to those proposed in Ref. Smit et al., 2002. Figure 3(a) shows $`𝒢(E)`$ for $`d=9.25`$Å (we chose the distance that more closely reproduced the vibrational mode frequency reported in Ref. Smit et al., 2002, see Table I). The conductance at the Fermi energy turns out to be remarkably smaller than $`𝒢_0`$. A similar calculation for the whole range of stable configurations shown in Fig. 2 gave a conductance in the range 0.2-0.5$`𝒢_0`$. To our ease this result agrees with the intuitive picture: Closed-shell molecules cannot conduct if the molecular character is maintained after contact with the electrodes and the charge transfer is much smaller than one. Both conditions apply here. The charge transfer from the H2 molecule to the electrodes is typically $`0.15`$ or smaller and the bonding and antibonding molecular levels maintain their character since they are clearly visible in the density of states (not shown here) around -6 eV and 17 eV, respectively. A value of conductance close to $`𝒢_0`$ is still possible through hybridization of the Pt states with those of the moleculeHeurich et al. . In our calculations, however, this does not occur at the Fermi energy, but at $`E1.5`$ eV. We believe that this H2 bridge actually accounts for a secondary peak around $`𝒢0.1𝒢_0`$ not discussed, but clearly present in the reported dataSmit et al. (2002).

So, what accounts for the main conductance peak at $`𝒢_0`$?. There is little doubt that this peak is H-related, but here we propose an alternative to the H2 molecular bridge discussed above. We now perform calculations similar to the previous ones, but initially placing the molecule axis perpendicular to the stretching $`z`$-axis. We let the Pt tip atoms and the H atoms relax towards equilibrium for several values of $`d`$ without assuming the H2 to be in molecular form. The total energy of the stable configurations we have found are represented by full black circles in Fig. 2. (Again, for completeness, we also present results when the H atoms are not allowed to leave the plane perpendicular to the $`z`$-axis). The two H atoms are only in molecular state at the upper limit of this curve ($`d8`$ Å). For smaller values of $`d`$ the H2 molecule breaks apart and the H atoms position themselves up to 2.0 Å apart from each other forming a double atomic H bridge as the one shown in the left inset of Fig. 2. The results of conductance for $`d=6.5`$Å (a distance for which vibrational modes are closer to the experimental ones, see below and Table I) are depicted in Fig. 3(b). The conductance at the Fermi energy approaches closely $`𝒢_0`$. The eigenchannel transmission results are also consistent with the noise analysis in Ref. Smit et al., 2002 since, the fact that fluctuations do not strictly vanish at $`𝒢_0`$ may indicate that more than one channel is contributing to the conductance. At this point the distance between tip Pt atoms is approximately that of bulk Pt and that between tip Pt atoms in Fig. 1(b). When comparing with the results presented in Fig. 1(b) we see that the H seems to block out one of the channels. It is important to finally stress here that the exact value of the conductance is not very dependent on the basis set and the functional. However, the atomic configuration of the electrodes in the cluster and the model used to represent the bulk of these electrodes introduces a source of uncertainty in the results. We believe that this is the source of the discrepancy between our results and those presented in Refs. Smit et al., 2002 and Heurich et al., . After considering various other possibilities in terms of cluster sizes and orientations, we can safely conclude that the stark difference between the conductance of the two basic bridge structures studied in this work is robust.

We have also performed a full analysis of the H-H vibrational modes in a configuration in which the H2 axis is perpendicular to the $`z`$-axis and for the distance that gives a conductance close to one quantum ($`d=6.5`$Å see Figs. 2 and 3). The results for the mode frequencies are reported in Table I. The first mode shown in the table (H atoms moving in opposite directions along the stretching axis) is the one that has a frequency closer to the experimental data (a frequency of 516$`\pm `$21 cm-1 was reported in Smit et al. (2002)). The isotopic effect shows up clearly with frequencies 439 and 386 cm-1 for H-D and D-D, to be compared with those obtained by multiplying the H-H frequency by the mass ratii (445 and 392 cm-1, respectively). It is worth mentioning that, although for D-D the isotopic effects shows up for all modes, this is no longer true for H-D. This mode induces a large variation of the conductance (see Table I), again, an effect that is not the same for all modes. Note also that the mode with the two H atoms moving in the same direction along the stretching axis, which also induces a significant change in the conductance (see second row in Table I), has a frequency close to that of a second mode observed in the experiments (around 1200 cm-1, see Ref. Smit, ). The latter is not reproduced by configurations with a molecular H2 bridge (see Table I).

In conclusion we have analyzed the conductance and vibrational modes of a Pt nanocontact in the presence and absence of H. Our results for clean Pt agree with most data reported in the literature. On the contrary, our results on Pt with H question the current interpretation of the experiments reported in Ref. Smit et al., 2002. While H2 seems to block the current, a novel complex Pt2H2 formed at the neck could be responsible for the reported histogram peak near $`𝒢_0`$. The vibrational modes of this complex also agree with the experimental ones.

We acknowledge support from the Spanish Ministry of Science and Technology under Grants No. 1FD97-1358 and No. MAT2002-04429-C03 and from the Universidad de Alicante. We thank C. Untiedt for sharing with us his point of view on the experimental results.