text:
# Metallicity, planetary formation and migration

## 1 Introduction

Planet searches, mostly via radial velocity surveys of nearby stars, have so far yielded some 110 planets (e.g. Schneider 2003). The stars that are harboring planets have been known for some time to be relatively metal-rich (e.g. Gonzalez 1998, 1999; Gonzalez, Wallerstein & Saar 1999; Reid 2002). While some of the early studies suffered from clear selection effects, since the planet hunters tended to concentrate on metal-rich stars, more recent studies tend to be merely brightness limited. In particular, Fischer & Valenti (2003) examined 754 stars in the solar neighborhood and found that the probability of hosting a planet rises from about 5 per cent when the iron abundance is $`1/3`$ that of the sun, to about 20 per cent, when the iron abundance is $``$3 times that of the sun. We will generally assume that the metal content of the stellar host reflects the metallicity of the planets/protoplanetary disc (although see Saumon et al. 1996). Suggestions have been made that the observed metallicity ‘excess’ may reflect ‘pollution’ of the star’s convective envelope by the accretion of planets (e.g. Laughlin & Adams 1997; Conzalez 1998). If that were true, however, one would expect stars with shallow convection zones (e.g. F stars) to show correspondingly higher metallicity enhancements. Such an effect is not observed (Fisher & Valenti 2003).

Broadly speaking, there are two scenarios for planet formation. One scenario involves a multi-stage process, in which dust grains coagulate to form rocks, then planetesimals, and the latter coalesce to form planetary embryos that accrete the gaseous envelopes (e.g. Goldreich & Ward 1973; Pollack et al. 1996; Lissauer 2001). In the other scenario, giant planets form from the gravitationally collapsing regions of unstable protoplanetary discs (e.g. Boss 2001; Mayer et al. 2002; Rice et al. 2003). On the face of it, the observed correlation between the metallicity of a star and its probability of hosting planets can be taken as supporting the multi-stage formation process. However, it should be realized that all the observed extrasolar planets (around main sequence stars) are believed to have migrated inward from larger radii (e.g. Lin, Bodenheimer & Richardson 1996; Armitage et al. 2002; Trilling, Lunine & Benz 2002), because planet formation at the observed radii appears to be difficult (e.g. Bodenheimer, Hubickyj & Lissauer 2000; Mayer et al. 2002). This raises the possibility, in principle, that the higher metallicity is required for the migration, rather than for the formation process (Sigurdsson et al. 2003). In the present Letter we examine this possibility.

## 2 Migration and metallicity

In order for the observed metallicity-planet frequency correlation to be an effect of migration alone, we need to explore the possibility that a low metallicity somehow inhibits migration. Mechanisms of migration that involve planet-planet scattering (e.g. Rasio & Ford 1996; Papaloizou & Terquem 2001) do not depend on the metallicity. The scattering of planetesimals by massive planets (e.g. Murray et al. 1998) does depend on the metallicity, but only if one requires the planetesimals to form via the multi-stage process. The only mechanism we need to study, therefore, is migration due to gravitational interaction with a gaseous, viscous, protoplanetary disc. The process we are interested in is “Type II” migration, in which the planet opens a gap in the disc, and subsequently follows its viscous evolution. We take the view that an observed Jupiter-sized planet is essentially the *last one* to have formed as the migration came to a halt (Armitage et al. 2002; Trilling, Lunine & Benz, 2003). In this case, assuming that Jupiter-sized planets have to form close to, or beyond, a radius of $``$5 AU, (e.g. Mayer et al. 2002; Boss 1995) we are interested in the properties of remnant discs at a time when the remaining mass of the disc was a few Jupiter masses at a radius of $``$5 AU. If the disc is much more massive than this, then the subsequent evolution of the disc will sweep the planet into the star.

There are two obvious ways, in principle, in which migration can be significantly affected: (i) If the viscosity in the disc is reduced to extremely low values due to the MHD turbulence dying away, and (ii) if migration becomes significantly slower due to global changes in the disc structure. In what follows, we examine each of these possibilities in turn.

### 2.1 Dead zones

The viscosity in accretion discs most probably originates from the magneto-rotational instability (MRI; Balbus & Hawley 1991). Numerical simulations have shown that the MRI, and the concomitant MHD turbulence are suppressed when the magnetic Reynolds number is smaller than some critical value (Gammie 1996; Gammie & Menou 1998; Fleming, Stone & Hawley 2000). This leads in some parts of the disc to a layered disc structure, in which the gas near the disc mid-plane is cold, shielded from ionizing radiation and non viscous. Accretion then occurs only in a very thin surface layer that is ionized by cosmic rays.

In a shear flow, the magnetic Reynolds number is defined as

$$Rm=LV/\eta ,$$

where $`L`$ and $`V`$ are typical length and velocity scales, respectively, and $`\eta `$ is the resistivity. Adopting $`LH`$ (the disc half-thickness), $`Vc_s`$ (the speed of sound), as is appropriate for the simulations that make use of the shearing box approximation (e.g. Fleming et al. 2000), the requirement $`RmRm_{\mathrm{crit}}`$ corresponds to balancing the growth rate of the MRI by Ohmic diffusion. When no externally imposed vertical magnetic field is present, the critical value is of order $`Rm_{\mathrm{crit}}2000`$ (Hawley, Gammie & Balbus, 1996). When one allows for an external vertical uniform field, $`Rm_{\mathrm{crit}}100`$ (Fleming et al. 2000). The decay of MHD angular momentum transport results in a “dead zone,” in which the viscosity is very low, and migration might be expected either to stop or to be slowed down significantly. Of course, if the dead zone does not significantly change the rate of migration, then the effect of metallicity on the dead zone is not relevant to our current considerations. Since the Reynolds number depends on the resistivity, which, in turn, depends on the electron fraction, a change in the metallicity could lead to changes in the formation of a dead zone that could lead to changes in the properties of migration. The resistivity is given by (e.g. Matsumura & Pudritz 2003)

$$\eta =234T^{1/2}/x_e\mathrm{cm}^2\mathrm{s}^1,$$

where $`T`$ is the mid-plane disc temperature, and $`x_e=n_e/n`$ is the electron fraction ($`n_e`$ and $`n`$ are the number densities of electrons and neutral atoms respectively). Thus, at fixed $`x_e`$, since $`Hc_s/\mathrm{\Omega }`$, where $`\mathrm{\Omega }R^{3/2}`$ is the angular velocity of the disc, we have that

$$RmT^{1/2}R^{3/2}.$$

Matsumara & Pudritz (2003) consider the radiative, hydrostatic disc model developed by Chiang & Goldreich (1997), and Chiang et al. (2001). This model has an assumed surface density of the form

$$\mathrm{\Sigma }=\mathrm{\Sigma }_o(R/\mathrm{AU})^{3/2},$$

with a fiducial value of $`\mathrm{\Sigma }_o=10^3`$ g cm-2 which gives a disc mass of a few Jupiter masses within about 5 AU. The flared, radiative equilibrium disc of Chiang & Goldreich (1997) has $`TR^{3/7}`$ for $`R<84`$ AU. Assuming a state of ionization balance (where the ionization can be caused by either the central star or cosmic rays), and using the equation of the electron fraction given by Oppenheimer & Dalgarno (1974), Matsumura & Pudritz (2003) solved for the radius of the dead zone in the two limiting cases of no metals at all, $`x_m=0`$, and essentially all metals ($`x_mx_e`$), where $`x_m`$ is metal fraction. Assuming a critical Reynolds number appropriate for a disc with a superimposed poloidal field (which they take to correspond to a critical Reynolds number $`Rm_{\mathrm{crit}}=1`$–10), they found that even at these extremes, the radius of the dead zone differed by at most a factor $`2`$ ($`1`$ AU in the metal dominant case, and $``$2 AU in the no-metals case).

For the fiducial disc we see that $`RmR^{9/7}`$. Thus, if one increases the critical Reynolds number by a factor of $`20`$–100, as is inferred from simulations with no imposed poloidal field, and which is more likely to be the case for a protoplanetary disc at $`R5`$ AU, then the radius of the dead zone becomes independent of the metallicity, and is determined solely by the surface density ($`\mathrm{\Sigma }100`$ g cm-2) that allows for ionization by cosmic rays (Gammie 1996). In the disc model used by Matsumura & Pudritz (taken from Chiang et al. 2001), this corresponds to a dead zone radius of $`5`$ AU. We therefore find that, if one just considers the effect of dead zones on the disc structure, and hence on the migration rate, then the migration rate is unaffected, even if the mass fraction in metals is changed from $`Z=1`$ to $`Z=0`$.

### 2.2 Disk structure and metallicity

We now consider whether a change in metallicity leads to a sufficient change in disc structure that the disc inflow timescale, $`\tau _\nu `$, and hence the planet migration rate, changes significantly. Since the disc models of Chiang & Goldreich (1997) are static, they are not suitable for this purpose. Ideally, one would like to consider a time-dependent model of a proto-planetary disc when its mass has reached the relevant value (of a few Jupiter masses at $``$5 AU). However, for simplicity and to get some indication of the sensitivity of the disc structure to metallicity, we use the disc models of Bell et al. (1997), who consider the detailed disc structure for a variety of values of steady accretion rates. In these models no account is taken of the possibility of “dead zones,” and the viscous process is parameterized using the usual $`\alpha `$-prescription (Shakura & Sunyaev, 1973). Recent calculations by Fleming & Stone (2003) indicate that this might in fact be a reasonable approach. Their model with an accretion rate of $`\dot{M}=10^8`$ M y-1 has a surface density of $`\mathrm{\Sigma }140`$ g cm-2 at a radius of 5 AU. Thus the mass of the disc at this radius is of order a few Jupiter masses, as required. What is relevant for our purposes is the viscous inflow timescale given by (Pringle 1981)

$$\tau _\nu R^2/\nu ,$$

where $`\nu \alpha c_s^2/\mathrm{\Omega }`$ is the viscosity. The Bell et al. (1997) models assume that the dimensionless viscosity parameter is $`\alpha =0.01`$. Thus at a given radius, $`\tau _\nu T^1`$, the mid-plane disc temperature. In the same model the (optically thick) disc temperature at $`R5`$ AU is $`T30`$ K. In this temperature range, Bell et al. (1997) find, using Henning & Stognienko (1996), that the opacity is approximately of the form $`\kappa =\kappa _oT^{2.1}`$, and that $`T\kappa _o^{0.34}`$. Since at these temperatures the opacity is predominantly due to dust (Pollack et al. 1994; Henning & Stognienko 1996), and although there may be complications due to the detailed structure of the dust and the associated chemistry, it seems reasonable as a first approximation to assume that $`\kappa _o`$, is proportional to the metallicity, i.e. $`\kappa _oZ`$. In this case a change in metallicity by a factor of $`10`$ leads to a corresponding change in the inflow timescale, and hence migration timescale by a factor of $`2`$.

## 3 Discussion and conclusions

The observed correlation between the metallicity of a star and its probability for hosting a planet seems to indicate that metallicity is an important factor in planet formation. However, the sample of detected extrasolar planets suffers from obvious selection effects (e.g. Zucker & Mazeh 2001; Tabachnik & Tremaine 2002; Lineweaver & Grether 2002; Armitage et al. 2002). In particular, radial velocity techniques have almost no sensitivity to planets with orbital periods longer than the duration of the survey. Consequently, if the migration process can be inhibited in the low metallicity systems, this would clearly produce a correlation of the type that has been observed.

In the present work, we have therefore examined potential effects of low-metallicity on migration. By considering the basic physical processes involved, we have shown that in low-metallicity systems, migration could be slowed down, but at most by a factor $`2`$ (in terms of timescale). Such a modest sensitivity is not expected to produce the observed rise in probability for hosting a planet (Armitage et al. 2002 have considered models in which the viscous timescale was changed by similar factors, with no significant consequences for the resultant orbital distribution). Therefore, our tentative conclusion is that the metallicity observations do support the idea that a higher metallicity is associated with a higher rate of formation of planets rather than with migration. Our conclusion is further supported by the fact that there is no distinguishable difference in metallicity between the stars hosting planets with semi-major axes longer than the observed median, and those with shorter ones (Fischer & Valenti 2003). To make this conclusion ironclad would require a detailed, time-dependent, migration computation, that would include a full treatment of the thermal structure of the disc.

It is interesting to note that the non-detection of transiting planets in the globular cluster 47 Tuc (Gilliland et al. 2000) is entirely consistent with the dependence of hosting a planet on metallicity. Originally, Gilliland et al. estimated an expected number of 17 detections (from monitoring 34,000 stars). If one takes account, however, of the fact (Marcy et al. 2003) that only 0.6% of the F, G, K, M stars in the solar neighborhood have “hot Jupiters” (giant planets with orbital periods $`P<5`$ days; Gilliland et al. assumed 1%), and the dependence of probability on metallicity, the number of expected detections is reduced to $`2`$.

If the formation of planets indeed requires metals, one needs to explain the recent observation of a planet in the globular cluster M4 (Sigurdsson et al. 2003). That cluster has a metallicity that is only 5% that of the sun, and therefore the formation of planets in it would be expected to be suppressed. It is important to note, however, that even in such a low metallicity environment, there are circumstances where dust and metals are abundant. One example is supernovae (note that the planet in M4 is a companion to the pulsar B1620$``$26), which can produce copious amounts of dust (2–4 solar masses in the case of Cas A; Dunne et al. 2003). Another dust-rich environment is around asymptotic giant branch stars. In fact, the white dwarf companion to the pulsar B1620$``$26 had to evolve through such a phase. It is not impossible that the circumbinary planet was in fact formed during that phase (and not around a main sequence star, as suggested by Sigurdsson et al.), and it was later pushed via interactions into a non-coplanar orbit. Finally, planets around pulsars (or white dwarfs) can form from the disruption of white dwarf companions (Livio, Pringle & Saffer 1992; Podsiadlowski et al. 1991).

## Acknowledgements

We are grateful to Jeff Valenti and Ron Gilliland for helpful discussions. JEP thanks STScI for hospitality and support from the Visitors’ Program. ML acknowledges support from John Templeton Foundation Grant 938-COS191. We acknowledge useful comments from an anonymous referee.


format: text

text:
# Parity-violating Electron Deuteron Scattering and the Proton’s Neutral Weak Axial Vector Form Factor

The SAMPLE Collaboration

## Abstract

We report on a new measurement of the parity-violating asymmetry in quasielastic electron scattering from the deuteron at backward angles at $`Q^2=0.038`$ (GeV/c)2. This quantity provides a determination of the neutral weak axial vector form factor of the nucleon, which can potentially receive large electroweak corrections. The measured asymmetry $`A=3.51\pm 0.57(\mathrm{stat})\pm 0.58(\mathrm{sys})`$ ppm is consistent with theoretical predictions. We also report on updated results of the previous experiment at $`Q^2=0.091`$ (GeV/c)2, which are also consistent with theoretical predictions.

preprint: Draft v4.5

Parity-violating electron scattering provides a unique probe of the electroweak structure of the nucleon. It has been well established that elastic scattering studies yield new and interesting information on the strange vector matrix elements KAP88 ; MCK89-BEC89 . This is the basis for a substantial program of experiments at modern electron accelerator facilities, beginning with the SAMPLE experiment MUE97-SPA00 ; HAS00  at MIT-Bates.

The primary goal of SAMPLE is to determine the proton’s strange magnetic form factor $`G_M^s`$ through parity-violating electron scattering from the proton at backward angles. However, the parity-violating asymmetry – the asymmetry in the scattering cross section with respect to the helicity of the incident electron – is not only sensitive to $`G_M^s`$, but is also sensitive to the proton’s neutral weak axial form factor. As pointed out in Ref. MUS90 , the neutral weak axial form factor as measured in electron scattering, $`G_A^e`$, can potentially receive large electroweak corrections that are absent in neutrino scattering. These corrections include the anapole moment, which is identified as the effective parity-violating coupling of a photon to the nucleon. Determining $`G_A^e`$ is important not only for a reliable extraction of $`G_M^s`$, but also because of its sensitivity to the hadronic effects on the electroweak radiative corrections. The adequate understanding of such effects is essential to proper interpretation of other precision electroweak measurements such as neutron and nuclear $`\beta `$ decay MCK03-RAM03 . Parity-violating quasielastic electron-deuteron scattering at backward angles is predominantly sensitive to $`G_A^e`$ and thus can be used to determine $`G_A^e`$ BEI91 .

The SAMPLE collaboration previously performed an experiment on a deuterium target (SAMPLE II) as well as on a hydrogen target (SAMPLE I) at 200 MeV \[$`Q^20.1`$ (GeV/c)2\]. Combining the results from these two experiments allows separate determination of $`G_M^s`$ and $`G_A^e`$. Our data HAS00  indicated that, while the overall contribution from strange quarks to the proton’s magnetic form factor is small, the size of the electroweak radiative corrections to the axial form factor is significantly larger than anticipated from theory ZHU00 .

These results stimulated considerable interest among theorists. Many different processes and effects were studied for their potential contributions to the axial form factor or to the parity-violating asymmetry in electron-deuteron scattering. These include the anapole moment ZHU00 ; MAE00a-MAE00b ; RIS00 , nuclear effects including two-body currents DIA01 , and the parity-violating hadronic interaction SCH03a ; LIU03 . None of the effects studied here were significant enough to explain the discrepancy.

In order to experimentally confirm these results, we performed a third SAMPLE experiment (SAMPLE III), with a deuterium target at a lower beam energy of 125 MeV \[$`Q^2=0.038`$ (GeV/c)2\]. As was the case for SAMPLE II, the dominant scattering process is quasielastic scattering, and the asymmetry is predominantly sensitive to $`G_A^e`$. Since the parity-violating asymmetry in the cross section is proportional to $`Q^2`$ to first order, the expected asymmetry was roughly 3 times smaller than that for 200 MeV. The cross section, however, is larger by a factor of 2 with the same background level, resulting in an experiment sensitive to the same physics with roughly the same sensitivity but with very different systematics.

The experiment was carried out at the MIT Bates Linear Accelerator Center. The experimental method and apparatus were identical to SAMPLE II, except for the incident beam energy. A 125 MeV longitudinally polarized electron beam was incident on a 40 cm long liquid deuterium target, and electrons scattered at backward angles were detected by an air Čerenkov detector covering angles between 130 and 170 (solid angle $`1.5`$ sr). The detector consists of the radiator air volume and 10 detector elements, each with an ellipsoidal mirror to focus Čerenkov light onto a corresponding 8-inch photomultiplier tube (PMT).

A remotely controlled light shutter was used to cover each PMT for the background measurements. About 10% of the data were taken with the shutters closed. In addition, in order to further study the background, an additional measurement of the shutter closed asymmetry was made with a plate of plastic scintillator placed in front of each PMT to enhance the statistics. Also, the non-Čerenkov sources of light in the detector signal, mostly due to scintillation light in the air, were studied by covering the mirrors. The background level was the same as SAMPLE II, consistent with bremsstrahlung radiation in the target being the dominant origin of all the different background components.

The incident electron beam was pulsed at 600 Hz, and the average beam current was 40 $`\mu `$A. The polarized electron beam was generated by directing a circularly polarized laser beam onto a GaAs crystal. The helicity of the beam was pseudorandomly chosen for each pulse. The helicity of the beam with respect to the electronic signal was manually reversed every 2-3 days by inserting and removing a halfwave plate in the laser beam path to check for and reduce possible systematic effects. (These two configurations are called “in” and “out”.) The polarization of the electron beam was measured daily with a transmission polarimeter, and occasionally with the Møller polarimeter on the beam line.

Various beam parameters, including the intensity, energy, position, and angle, were monitored continuously for helicity correlated differences. Four Lucite Čerenkov counters (luminosity monitors) downstream of the target at the forward angles ($`12^{}`$) detected low $`Q^2`$ scattering which has negligible parity-violating asymmetry, thus serving as monitors for false asymmetries.

As in the past, the yield for each detector for each beam pulse (integrated over the duration of the pulse) was normalized to the beam charge measured in front of the target, and then was corrected for the beam position, angle, and energy. A linear regression technique was used to determine the dependence of the detector yields on each parameter. Then, the asymmetry was computed for the appropriate pulse pairs. In addition, the normalized detector yields were also corrected for transmission of the beam (defined as the ratio of the beam intensities at the target and at the end of the accelerator) to account for an observed beam intensity asymmetry caused by differential scraping at an energy-defining slit. The results of this analysis are shown in Fig. 1, where the detector asymmetry is plotted as a function of time.

The asymmetry was further corrected for the beam polarization, the background dilution, and electromagnetic radiative effects (effects due to the bremsstrahlung radiation of the incident and scattered electrons) to obtain the physics asymmetry. The average beam polarization during the experiment was $`P_e=38.9\pm 1.6`$%. The background dilution factor, determined for each detector from the ratio between the shutter open and closed detector yields and the mirror covered studies, was typically $`1.41.7`$ with a relative uncertainty of 4.5%. Electromagnetic radiative effects were evaluated using a spin-dependent modification to Ref. MO69  within the context of a geant GEANT  simulation of the detector geometry. In the simulation, scattered electron events were generated uniformly in energy, angle, and along the length of the 40 cm target. The scattered electron kinematics were selected after accounting for energy loss in the target. Each event was weighted according to the scattering cross section and the detector efficiency, and was assigned an asymmetry according to its kinematics. The correction factor for electromagnetic radiative effects was evaluated for each detector by comparing the (weighted) asymmetry with and without the radiative effects included in the simulation, and was typically 1.09 with a relative uncertainty of 3%.

The systematic error in the corrections procedure was estimated by comparing results from two different methods that are mathematically equivalent for normally distributed infinite data: one computes the dependence of the detector signal on the beam parameters for normalized yields, and the other for asymmetries. We assign a relative systematic error of 11.2% for out and 2.1% for in, the larger error for out naturally reflecting the larger correction due to the larger beam intensity asymmetry.

Additional uncertainties were assigned to the resulting physics asymmetry to account for two systematic effects observed during the experiment. The first is the residual asymmetry in the luminosity monitors. Some of the luminosity monitors showed non-zero asymmetries even after the corrections procedure was applied, potentially indicating the existence of a helicity correlated difference in some unmeasured beam parameter(s) that caused false asymmetries in the luminosity monitor signal. The size of the false asymmetry that this effect could cause in the Čerenkov detector signal was estimated from the observed luminosity monitor asymmetries and the correlation between the Čerenkov detector asymmetry and the luminosity monitor asymmetry, and was assigned as the systematic error. Relative systematic errors of 20.0% and 19.2% were assigned for the out and in data, respectively, and the errors were treated as uncorrelated when combining the two data sets.

The second is that, although the measured shutter closed asymmetry for all 10 detectors combined was consistent with zero, the individual detectors showed a non-zero shutter closed asymmetry. The detector-by-detector distribution showed a definite pattern dependent on the azimuthal angle, indicating that this asymmetry is of parity-conserving nature, and hence cancels out when averaged over all 10 detectors that are symmetrically arranged azimuthally. The shutter closed asymmetry was estimated from the “high-statistics” shutter closed data taken with plastic scintillator and was subtracted from the shutter open asymmetry for each detector. The value of the final asymmetry is very insensitive to this procedure because of the symmetry of the detector arrangement, and the associated systematic error was estimated to be 5%.

The resulting physics asymmetry is

$$A(Q^2=0.038)=3.51\pm 0.57\pm 0.58\mathrm{ppm},$$
(1)
where the first uncertainty is statistical and the second is the estimated systematic error as summarized in Table 1.

Since the SAMPLE detector does not have energy resolution for the scattered electrons, the measured asymmetry contains contributions not only from quasielastic scattering but also from elastic scattering and threshold breakup. In order to construct the theoretical expression of the asymmetry as a function of the quantities of interest, i.e., $`G_M^s`$ and $`G_A^e`$ of the nucleon, we did the following. First we performed a full nuclear calculation according to Ref. SCH03a  to obtain the parity-conserving and parity-violating response functions for the total inelastic processes (quasielastic scattering and threshold breakup) for selected kinematics. The dependence on $`G_M^s`$ and $`G_A^e`$ was explicitly kept track of in the calculation. Electroweak radiative corrections were included. In particular, the isoscalar axial radiative correction was taken to be $`R_A^0=0.03\pm 0.05`$ from Ref. ZHU00 . We use $`\mathrm{sin}^2\theta _W=0.23113(15)`$ PDG .

The parity-violating asymmetry was computed on an event-by-event basis in the geant simulation, and separately for the elastic (from Ref. POL90 ) and inelastic (using the above obtained response functions) processes. The resulting asymmetry distributions represented an average over the detector acceptance and incident electron energies. The physics asymmetry was then computed as a combined average of the elastic and inelastic distributions weighted by the appropriate cross sections. The resulting theoretical asymmetry is

$$A(Q^2=0.038)=2.14+0.27G_M^s+0.76G_A^{e(T=1)},$$
(2)
where the asymmetry is in parts per million and the form factor is in nuclear magnetons (n.m.). In this expression, we retain explicitly the isovector ($`T=1`$) component of $`G_A^e`$. The small isoscalar component is absorbed into the first term. The dependence on the nuclear model is small.

The radiative corrections and theoretical asymmetry for the SAMPLE II data were also re-evaluated with the geant simulation. In addition, background dilution factors coming from pion photoproduction were re-examined in light of recently published data BER98 . Such processes contribute to the detector yield through their decay products, but have negligible parity-violating asymmetry CHE01 . The largest contribution is from coherent $`\pi ^0`$ photoproduction on the deuteron, which had been neglected in Ref. HAS00 , but was found in Ref. BER98  to be significantly enhanced relative to the corresponding incoherent process. Including this effect increased the background dilution factor by 9%. The re-evaluated electromagnetic radiative corrections resulted in an additional 2% increase in the background dilution factor. Finally, improved determination of the scintillation component of the detector signal resulted in another 2% increase. Thus, the final physics asymmetry increased by 13% in magnitude compared with our previously published results HAS00 , giving

$$A(Q^2=0.091)=7.77\pm 0.73\pm 0.62\mathrm{ppm},$$
(3)
where the first error is statistical and the second is the estimated systematic error. (Note that this asymmetry value contains the contribution from the non-quasielastic processes, which was estimated to be $`1.5`$% and removed in Ref. HAS00 .) The re-evaluated theoretical value for the asymmetry, using the nuclear calculation as described above, has resulted in an expected value that is 2% smaller:

$$A(Q^2=0.091)=7.06+0.77G_M^s+1.66G_A^{e(T=1)}.$$
(4)
In Fig. 2, the physics asymmetries measured in SAMPLE II (updated results) and SAMPLE III are plotted as a function of $`Q^2`$. Also plotted are the theoretical predictions with the value of $`G_A^e`$ taken from Ref. ZHU00  \[$`G_A^e(Q^2=0.038)=0.91\pm 0.28`$ and $`G_A^e(Q^2=0.091)=0.84\pm 0.26`$\], and $`G_M^s=0.15`$ n.m. The dependence of the theoretical values on $`G_M^s`$ is small.

The results from SAMPLE III (125 MeV deuterium run) and the updated results from SAMPLE II (200 MeV deuterium run) both agree with the theoretical prediction on the electroweak radiative correction on the neutral weak axial form factor of the nucleon by Zhu et al. ZHU00 . In addition to these two experimental results, various theoretical efforts also support the theoretical prediction by Zhu et al. The confirmation on the theoretical value of $`G_A^e`$ not only allows us to extract $`G_M^s`$ reliably from the data from SAMPLE I (200 MeV hydrogen run), but also is important for interpreting results from future parity-violating electron scattering experiments at JLab and Mainz.

We gratefully acknowledge the skillful efforts of the staff of the MIT-Bates laboratory to support the experiment and provide a high quality beam. We thank M. J. Ramsey-Musolf for insightful discussions. This work was supported by NSF, DOE, Jeffress Memorial Trust, and CNRS of France.


text:
# A WARM HOT INTERGALACTIC MEDIUM TOWARDS 3C 120?

## 1. Introduction

Around half of the baryons in the universe (by mass) should reside in intergalactic filaments of a warm/hot intergalactic medium (WHIGM) (Cen et al. 1995, Davé et al. 2001). These filaments are believed to be shock-heated to $`10^5\mathrm{K}<\mathrm{T}<10^7\mathrm{K}`$ (see e.g. Cen & Ostriker, 1999). Most WHIGM may lie at the hotter end of this temperature range (Cen et al. 2001), so high spectral resolution X-ray detectors such as those aboard Chandra and XMM-Newton are best placed for investigating the ‘hot’ component of the missing baryons. ‘Hot’ WHIGM may recently have been discovered for the first time at X-ray energies (Nicastro et al. 2002; Fang et al. 2002a,b; Rasmussen et al. 2002; McKernan et al. 2003).

We observed the X-ray source 3C 120 with the High Energy Transmission Grating Spectrometer (or HETG–Markert, et al. 1995) and ACIS aboard Chandra . 3C 120 (z=0.033, Michel & Huchra 1988) is classified both as a Seyfert 1 active galactic nucleus (AGN) and a broad line radio galaxy (BLRG). Here we discuss the serendipitous discovery of an O viii Ly$`\alpha `$ absorption signature in the 3C 120 spectrum, likely due to WHIGM along the line of sight.

## 2. Observations and Data

We observed 3C 120 with the HETGS on board Chandra on 2001 December 21 for $`58`$ ks, beginning at UT 11:13:52. The HETGS consists of two grating assemblies, a High-Energy Grating (HEG) and a Medium-Energy Grating (MEG). Only the summed, $`\pm 1`$ order Chandra grating spectra were used in our analysis. We found the counts distribution as a function of cross-dispersion angle could be approximated by a Gaussian model with widths in the range $`0.460.49\mathrm{}`$ for the $`\pm 1`$ orders of the MEG. This Gaussian model is consistent with a point source, although there is emission of $`6\%`$ peak emission at $`+3\mathrm{}`$, possibly due to the extended jet. We accumulated grating spectra along the dispersion direction and within $`\pm 3.6\mathrm{}`$ of the peak in the cross-dispersion direction. The mean MEG and HEG total count rates were $`0.8134\pm 0.0038`$ and $`0.3770\pm 0.0031`$ cts/s respectively. We extracted spectra and spectral responses exactly as described in Yaqoob et al. (2003) and obtained a net exposure time of 57,224 s (including a deadtime factor of 0.01618). During the observation 3C 120 did not vary much; the light-curve yielded an excess variance of $`0.0043\pm 0.0007`$ (Turner et al. 1999). Since the HEG bandpass only extends down to $`0.8`$ keV we use the MEG as the primary instrument. We used $`C`$-statistic for finding best-fit model parameters, and quote 90$`\%`$ confidence, one-parameter statistical errors. The harder spectrum, including the Fe-K emission line will be discussed elsewhere.

The ACIS CCDs have been undergoing a low-energy QE degradation, due to absorption by contaminants 111http://cxc.harvard.edu/cal/Links/Acis/acis/Cal\_prods/qeDeg/index.html. There is a model of the contamination pertinent for pure ACIS data (without the gratings). We can use this model as a ‘worst case’ to estimate the effects on the continuum. Thus, we will give results with and without corrections using the ACIS QE degradation model in XSPEC v11.2 (acisabs), using default absorption by the contaminants C, H, O, and N in the ratio 10:20:2:1 by number of atoms.

All spectral fits were done in the 0.5–5 keV energy band, excluding the 2.0–2.5 keV region, which suffers from systematics due to limitations in the calibration of the X-ray telescope 222http://asc.harvard.edu/udocs/docs/POG/MPOG/node13.html. First we fitted spectra binned coarsely at $`0.32`$A, in order to compare with previous CCD spectra. A simple model consisting only of a single power–law ($`\mathrm{\Gamma }=1.69\pm 0.02`$) and Galactic absorption of $`1.23\times 10^{21}\mathrm{cm}^2`$ (Elvis, Lockman & Wilkes 1989) fit the data adequately, with residuals $`<`$ 30 $`\%`$ over the fitted energy band. The inverted photon spectrum and residuals are shown in Fig. 1. The residuals at $``$ 0.8 keV coincide with sharp changes in effective area and are likely due to the limitations of the calibration. Including the ACIS degradation model requires a soft excess in the data, which we modeled as a broken power–law. We found a best-fit soft X-ray index of $`2.58\pm 0.02`$, a hard X-ray index of $`1.76\pm 0.02`$ and a break energy of $`1.06_{0.06}^{+0.04}`$ keV.

We compared our data with previous CCD observations by modeling the MEG data with two absorption edges and a single power law (and Galactic absorption). We find threshold optical depths of $`\tau <0.06`$ and $`\tau =0.06_{0.05}^{+0.05}`$ at $`90\%`$ confidence for O vii and O viii respectively. These small optical depths do not alter the power-law index. These optical depths agree reasonably well with $`\tau <0.01`$ (O vii ) and $`\tau <0.05`$ (O viii ) found by ASCA (Reynolds, 1997). Including ACIS degradation yielded $`\tau <0.11`$ and $`\tau <0.07`$ for O vii and O viii respectively. Finer spectral binning ($`0.02`$A) does not change these results.

## 3. Discrete Spectral Features

The MEG spectrum of 3C 120 is shown in detail in Fig. 2(a) (0.6–0.8 keV) and Fig. 2(b) (0.8–1.4 keV). Strong absorption appears at $`0.645`$ keV and $`0.71`$ keV in the observed frame in Fig. 2(a). We identify the feature at $`0.645`$ keV with O viii Ly$`\alpha `$ (0.653 keV rest–frame) blueshifted by $``$ 5500 km $`\mathrm{s}^1`$ from 3C 120 systemic velocity. The feature at $`0.71`$ keV remains unidentified since all absorption transitions considered for this feature should yield strong absorption features elsewhere in the spectrum. For all spectral fitting hereafter we use a binsize of $`0.02`$A ( $``$ MEG FWHM resolution). Since there are bins with zero counts in the O viii Ly$`\alpha `$ feature, the data are of limited statistical quality, modeling the absorption profile with a Voigt function is not warranted. Therefore we modeled the profile with a simple window function, whereby over an energy interval $`W`$, centered on some energy, $`E`$, a fraction $`f`$, of the continuum is absorbed.

Table 1 shows the results of spectral fitting the O viii Ly$`\alpha `$ profile. Also in Table 1 are the corresponding results when ACIS degradation is included, and the best–fit to the unidentified feature at $`0.71`$ keV (also fit with a window function plus ACIS degradation). The O viii Ly$`\alpha `$ profile parameters are not strongly influenced by ACIS degradation and the C–statistic improves by $`>18`$ upon the addition of the absorption-line model (which corresponds to $`>3\sigma `$ significance for three additional parameters). The O viii Ly$`\alpha `$ feature has FWHM$`=1010_{265}^{+295}`$ km $`\mathrm{s}^1`$, EW$`=2.17_{0.63}^{+0.57}`$ eV and is strongly blueshifted ($`\mathrm{\Delta }v=5500\pm 140`$ km $`\mathrm{s}^1`$) from 3C 120 systemic velocity. The absorption occurs at $`z=0.0147\pm 0.0005`$ if it is due to intervening WHIGM. O viii Ly$`\beta `$ absorption is not statistically significant (EW$`=0.70_{0.70}^{+1.19}`$ eV) when the data are fit with a Gaussian model with the same FWHM and $`z`$ as O viii Ly$`\alpha `$ in Table 1. We tested the data for O vii (r) , Ne ix (r) and Ne x Ly$`\alpha `$ absorption in the same way and found EW$`<`$4.12 eV and $`<`$1.78 eV for O vii (r) and Ne x Ly$`\alpha `$ respectively. There is weak evidence of Ne ix (r) emission so we could not find a meaningful upper limit on the absorption EW for this feature.

An independent confirmation of the significance of the O viii Ly$`\alpha `$ feature was obtained as follows. From the best-fitting continuum, the instrument effective area, and the exposure time, we estimated if there were no absorption feature we should expect a mean of 13.8 photons in the fitted energy interval containing the feature (Table 1). We observed 2 photons in this energy interval (13 bins out of 126 over the range $`z=0`$ to $`z=0.033`$). The Poisson probability of obtaining this result is $`0.00975\times (126/13)\%`$ so the confidence level of the detection of the feature is $`>99.9\%`$ (i.e. $`>3\sigma `$). We verified this result by performing $`10^7`$ simulations of the continuum, with no absorption feature. We note that a 9% decrease in the level of the continuum is allowed before the detection significance of the absorption feature drops below $`3\sigma `$. Both the O viii Ly$`\alpha `$ and 0.71 keV absorption features are detected separately in the -1 and +1 arms of the MEG. Absorption models improved the C–statistic by $`17.3,16.0`$ for O viii Ly$`\alpha `$ and by $`24.8,17.4`$ for the $``$ 0.71keV feature in the $`1,+1`$ arms respectively. Thus it is extremely unlikely that either feature is due to a statistical fluctuation since the probabilities above do not take into account detections in *both* arms of the MEG.

Fig. 2 (c) shows a velocity spectrum centered on the O viii Ly$`\alpha `$ feature. Negative velocity indicates blueshift relative to the observed energy of the O viii Ly$`\alpha `$ feature. Superimposed is the best-fitting continuum model and absorption line profile as described below (see also Table 1). The dashed lines correspond to rest–frame velocities at $`z=0`$ and $`z=0.033`$ respectively. The O viii Ly$`\alpha `$ is blueshifted by $`5500`$ km $`\mathrm{s}^1`$ from 3C 120 systemic velocity ($`z=0.033`$) and there may be absorption at $`4500`$ km $`\mathrm{s}^1`$ in Fig. 2 (c), close to $`z=0`$, presumably due to a hot medium in the Local Group, as observed elsewhere (eg Nicastro et al. 2002, Fang et al. 2002b, Rasmussen et al. 2002, McKernan et al. 2003). However, the relativistic jet in 3C 120 could also account for the large outflow velocity (see §4). There is no absorption near the 3C 120 systemic velocity ($`z=0.033`$).

## 4. Discussion

We identify the most prominent absorption feature in the soft X-ray spectrum of 3C 120 with blueshifted O viii Ly$`\alpha `$ ($`5500`$ km $`\mathrm{s}^1`$ relative to systemic, see Table 1). There may also be absorption at $`z=0`$ due to O viii Ly$`\alpha `$ in the Local Group of galaxies and there is unidentified absorption at $`0.71`$ keV. The O viii Ly$`\alpha `$ outflow velocity is considerably larger than in ‘warm absorbers’ in Seyfert 1 AGN (typically a few hundred km $`\mathrm{s}^1`$). O viii Ly$`\alpha `$ absorption possibly due to intervening, non-local WHIGM has been observed in the spectrum of the BL Lac PKS-2155-304 by Fang et al. (2002a) who report a similar width ($`<1450`$ km $`\mathrm{s}^1`$) to our feature ($`1010_{265}^{+295}`$ km $`\mathrm{s}^1`$), but weaker (EW$`=0.48_{0.19}^{+0.25}`$ eV versus EW$`=2.17_{0.78}^{+0.57}`$ eV). Absorption by a hot medium in the Local Group ($`z0`$) has also been detected towards PKS-2155-304 (Fang et al. 2002a, Nicastro et al. 2002), 3C 273 (Rasmussen et al. 2002, Fang et al. 2002b) and NGC 4593 (McKernan et al. 2003), but these features are generally weaker than the O viii Ly$`\alpha `$ feature in 3C 120 .

From Table 1, O viii Ly$`\alpha `$ FWHM$`<1305`$ km $`\mathrm{s}^1`$, so the path length of a putative WHIGM filament has an upper limit $`19h_{70}^1`$ Mpc, otherwise differential Hubble flow would broaden the line. Also, using EW$`>1.39`$eV and a curve–of–growth analysis, we obtain a column density of $`N_{O\mathrm{viii}}>3.4\times 10^{16}\mathrm{cm}^2`$ (at 90$`\%`$ confidence). The EW upper limits on O vii (r), O viii Ly$`\beta `$ and Ne x Ly$`\alpha `$ discussed in §3 above, are consistent with this curve–of–growth analysis. Assuming all O is in O viii , and an O abundance of $`0.1`$ solar ( or $`8.5\times 10^5\mathrm{A}_\mathrm{H}`$), the corresponding neutral Hydrogen column is $`4.0\times 10^{20}\mathrm{cm}^2`$. For a filament depth $`<19h_{70}^1`$ Mpc, we obtain an electron density $`n_e>6.7\times 10^6\mathrm{cm}^3`$. If we assume the mean baryon density in the universe ($`\mathrm{\Omega }_b`$) relative to the critical density ($`\mathrm{\Omega }_{\mathrm{crit}}=9.2\times 10^{30}h_{70}^2`$ g $`\mathrm{cm}^3`$) is $`\mathrm{\Omega }_bh_{70}^2=0.0224\pm 0.0009`$ (Spergel et al. 2003), then using $`m_\mathrm{H}=1.7\times 10^{24}`$ g, the mean number density of baryons ($`\overline{N}_b`$) in the universe is $`\overline{N}_b=1.2\times 10^7\mathrm{cm}^3`$. Equating $`n_e`$ with $`\overline{N}_b`$, we find a baryon overdensity ($`O_b>`$56$`\mathrm{\Omega }_bh_{70}^2`$) in the WHIGM. From WHIGM simulations for a variety of $`\mathrm{\Lambda }`$-dominated cold dark matter universes, $`O_b`$ in WHIGM is expected to peak in the range $`1030`$, for distributions where $`7080\%`$ of WHIGM baryons lie in the range $`O_b5200`$ (Davé et al. 2001). Our results indicate that WHIGM may indeed lie in diffuse, hot, intergalactic filaments, slightly denser than predicted by simulations. The high ionization state of the WHIGM filament is consistent with previous observations of 3C 120 in the UV with IUE , which exhibit no obvious Ly$`\alpha `$ absorption (Kinney et al. 1991). There is also no evidence for absorption in the optical band (Baldwin et al. 1980) in 3C 120. The absence of strong absorption features from Ne ix (r) and Ne x Ly$`\alpha `$ may indicate that metal abundances in WHIGM are considerably different from solar abundances.

Is it possible that the absorption feature in 3C 120 is actually due to intrinsic absorption in the jet? The O viii Ly$`\alpha `$ feature could originate in an outflow of $`5500`$ km $`\mathrm{s}^1`$ relative to systemic velocity and many QSOs can exhibit such high velocity outflows. Optically bright QSOs exhibit C iv absorption signatures at 5,000–65,000 km $`\mathrm{s}^1`$ (Richards et al. 1991). Very broad ($`30,000`$ km $`\mathrm{s}^1`$) absorption features have been observed in BL-Lacs with ROSAT (Madejski et al. 1991), however these are much broader than the features in 3C 120. XMM-Newton observations have found high-velocity outflows in PG1211+143 ($`0.08c`$) and PG0844+349 ($`0.2c`$) (Pounds et al. 2003a,2003b). If absorption signatures in 3C 120 arise in the jet, they might extend nearly to systemic velocity due to jet acceleration or deceleration. A hint of absorption around $`2500`$ km $`\mathrm{s}^1`$ in Fig. 2 (c), suggests we cannot rule this out. Interestingly, all the other X-ray sources where WHIGM detection has been claimed (Fang et al. 2002a,b; Nicastro et al. 2002; Rasmussen et al. 2002) possess a jet closely aligned to the line-of-sight. However X-ray sources with jets exhibiting WHIGM absorption may simply represent a selection effect. These are more luminous, distant X-ray sources, so are more likely to have detectable WHIGM filaments along their line-of-sight.

We acknowledge support from NSF grant AST0205990 (BM), NASA grants NCC-5447 (TY), NAG5-7385 (TJT) and CXO grant GO2-3133X (TY, BM). We used HEASARC online data archive services, supported by NASA/GSFC and the NASA/IPAC Extragalactic Database (NED) supported by JPL. Thanks to the referee for valuable suggestions to improve this letter.


text:
# The Simulation of the Inelastic Impact

## 1 Introduction

The coefficient of restitution (COR) $`e`$ is introduced to determine the post-collisional velocity in the normal collision of two materials and defined by $`𝐮^{^{}}𝐧=e𝐮𝐧`$, where $`𝐮`$ and $`𝐮^{^{}}`$ are respectively the velocity of the contact point of two colliding materials before and after the impact, and $`𝐧`$ is the normal unit vector of the tangential plane of them. In the oblique impact, it has become clear that COR depends on the incident angle as well as on the impact velocity. Louge and AdamsLouge and Adams (2002) recently reported that COR increases as a linear function of the tangent of the incident angle in the oblique impact of a hard aluminum oxide sphere on a thick elastoplastic plate. They also suggested that COR can exceed $`1`$ in most grazing impacts. In this proceeding, we carry out the two-dimensional simulation of the oblique impact and investigate the dependency of COR on the incident angle based upon the theory of elasticity.

## 2 Model

Our numerical model is a two-dimensional model which is composed of an elastic disk and an elastic wallKuninaka and Hayakawa (2003). Each of them is composed of randomly placed mass particles connected by nonlinear springs each other. Numbers of mass particles are $`400`$ for the disk and $`2000`$ for the wall. The width and the height of the wall are $`8R`$ and $`2R`$ respectively where $`R`$ is the radius of the disk.

Both sides and bottom of the wall are fixed. Spring potential is described as $`V(x)=\frac{1}{2}k_ax^2+\frac{1}{4}k_bx^4`$, where $`x`$ is a stretch from natural length. The values of $`k_a`$ are $`1.0\times mc^2/R^2`$ for the disk and $`k_a=1.0\times 10^2mc^2/R^2`$ for the wall, where $`c`$ is the one-dimensional sound of velocity. And we adopt $`k_b=k_a\times 10^3/R^2`$ for each of them. Poisson’s ratio $`\nu `$ can be evaluated from the strains of the band of random lattice in vertical and horizontal directions to the applied force. We obtain Poisson’s ratio $`\nu `$ of the model as $`\nu =(7.50\pm 0.11)\times 10^2`$.

In our simulation, we define the incident angle $`\gamma `$ by the angle between the normal vector of the wall and the initial velocity vector of the disk(Fig.(1)). We fix the initial colliding velocity of the disk as $`|𝐯(0)|=0.1c`$ to control the normal and tangential components of the initial colliding velocity as $`v_t(0)=|𝐯(0)|\mathrm{sin}\gamma `$ and $`v_n(0)=|𝐯(0)|\mathrm{cos}\gamma `$, respectively. From the normal components of the contact point velocities before and after the collision, we calculate COR for each $`\gamma `$. We use the fourth order symplectic numerical method for the numerical scheme of integration with the time step $`\mathrm{\Delta }t=10^3R/c`$.

## 3 Results and Discussion

Figure 2 is the COR against the tangent of the incident angle $`\mathrm{\Psi }_1=\mathrm{tan}\gamma `$ in our simulation. The cross points are the average and the error bars are the standard deviation of 100 samples for each incident angle. This result shows that the COR increases as $`\mathrm{\Psi }_1`$ increases to exceed $`1`$, and has a peak around $`\mathrm{\Psi }_1=5.0`$. This behavior is contrast to that in the experiment by Louge and AdamsLouge and Adams (2002).

To explain this result, we consider the correction of COR by the local deformation of the wall. We assume that the normal unit vector $`𝐧`$ to the surface of the wall rotates toward the incoming disk by an angle $`\alpha `$ to become $`𝐧^\alpha `$. If we define $`e=(𝐮^{^{}}𝐧)/(𝐮𝐧)`$ and $`e^\alpha =(𝐮^{^{}}𝐧^\alpha )/(𝐮𝐧^\alpha )`$, the relation between $`e`$ and $`e^\alpha `$ becomes

$$e=(e^\alpha +\mathrm{\Psi }_2^\alpha \mathrm{tan}\alpha )/(1\mathrm{\Psi }_1^\alpha \mathrm{tan}\alpha ),$$
(1)
where $`\mathrm{\Psi }_1^\alpha =(\mathrm{\Psi }_1\mathrm{tan}\alpha )/(1+\mathrm{\Psi }_1\mathrm{tan}\alpha )`$ and $`\mathrm{\Psi }_2^\alpha =(\mathrm{\Psi }_1\mathrm{tan}\alpha )/(1+\mathrm{\Psi }_1\mathrm{tan}\alpha )3(1+e^\alpha )(\mu +\mathrm{tan}\alpha )/(1\mu \mathrm{tan}\alpha )`$. As for $`\mathrm{\Psi }_2^\alpha `$, we use the phenomenological theory for the oblique impact by Walton and BraunWalton and Braun (1986). The correction angle $`\alpha `$ can be estimated by the theory of elasticity. If we express the contact area by a parabora, $`\mathrm{tan}\alpha `$ equals to $`|x_cx_a|(12\theta )/R(22\theta )`$ with $`\theta =(1/\pi )\mathrm{arctan}(12\nu )/(\mu (22\nu ))`$, where $`\mu `$ is the coefficient of friction and $`x_c`$ and $`x_a`$ are the x coordinates of both ends of the contact area. $`\mu `$ can be calculated from the simulation data through the definition $`\mu =|J_t|/|J_n|`$, where $`J_n`$ and $`J_t`$ are the normal and tangential components of the impulse. The relation between $`\mu `$ and $`\mathrm{\Psi }_1`$ are cross points in Fig. 3. The solid curve in Fig. 2 is Eq. (1) with $`e^\alpha =0.95`$ which is COR in the normal impact. Our numerical results can be reproduced by our phenomenological theory.

The relation between $`\mu `$ and $`\mathrm{\Psi }_1`$ can be explained as follows. We assume that jags are uniformly placed on the surface of the wall with the density $`\rho `$ per unit length and the tangential velocity of the disk is decreased by $`\eta `$ when the disk interacts with the one jag. The tangential and normal impulses can be calculated by calculating the number of jags the disk interacts during collision timeHayakawa and Kuninaka (2002) as $`J_t=m\eta \rho |𝐯(0)|\mathrm{sin}\gamma \pi (R/c)\sqrt{\mathrm{ln}(4c/|𝐯(0)|\mathrm{cos}\gamma )}`$ and $`J_n=m(e+1)|𝐯(0)|\mathrm{cos}\gamma `$. We also assume that the tangential impulse decreases by $`J_t^{^{}}=m\zeta |𝐯(0)|\mathrm{sin}\gamma `$ which is proportional to the initial tangential velocity with the proportionality constant $`\zeta `$. Thus, $`\mu `$ can be calculated by the ratio of $`|J_tJ_t^{^{}}|`$ to $`J_n`$ as

$$\mu =\left|\zeta \mathrm{tan}\gamma \eta \rho \mathrm{tan}\gamma \pi \frac{R}{c}\sqrt{\mathrm{ln}\left(\frac{4c}{|𝐯(0)|\mathrm{cos}\gamma }\right)}\right|/(e+1).$$
(2)
The solid curve in Fig.3 is Eq.(2) with $`\zeta =0.317`$ and $`\eta \rho =0.0416c/R`$. Our numerical results can be well reproduced by our phenomenological theory of the coefficient of friction.

## 4 Summary

In the present study, we have shown that COR can exceed $`1`$ also in our two-dimensional simulation and depends on $`\mu `$ in the oblique impact. Our results can be explained by our simple phenomenological theory.

This study is partially supported by the Grant-in-Aid of Ministry of Education, Science and Culture, Japan (Grant No. 15540393).


text:
# The trilinear neutral Higgs self-couplings in the MSSM. Complete one-loop analysis

## Introduction

Two basic elements of the gauge boson and fermion mass generation in the Standard Model (SM) and its minimal supersymmetric extension (MSSM) are (1) the spontaneous symmetry breaking (SSB)  and (2) the Higgs mechanism . Gauge bosons and fermions gain masses via gauge-invariant interaction with scalar fields, which have non-zero vacuum expectation values. The SSB and Higgs mechanism respect gauge invariance and renormalizability of the model. In order to confirm experimentally the mechanism of mass generation, it is needed (1) to detect the Higgs bosons and measure their masses (2) to measure the Higgs boson couplings with gauge bosons and fermions (3) to determine the Higgs boson self-couplings. The last step is important especially for the case of MSSM, where the self-couplings are determined by the soft supersymmetry breaking principle. Tree-level analysis shows that some neutral Higgs boson self-couplings can be measured  with relatively high accuracy on the future high-luminosity colliders. The leading one-loop corrections to the neutral Higgs self-couplings may be expected to be sizable in some regions of MSSM parameter space, with the potential to change substantially the tree-level results.

We calculate two self-couplings of the neutral MSSM Higgs bosons $`h^0`$, $`H^0`$ taking into account the complete set of one-loop diagrams. Our main interest is the analysis of modifications in the Higgs potential by the corrections coming from the scalar and sparticle sectors of the model.

In section 1 the Higgs boson interaction Lagrangian and the self-couplings are presented. In sections 2 and 3 we discuss technical details of the approach. Some numerical results and their comparative analysis are contained in section 4.

## 1 Higgs boson interaction Lagrangian

The Higgs sector of MSSM includes five physical fields:

$``$ two neutral $`CP`$ even Higgs fields $`\{h^0,H^0\}`$;

$``$ neutral $`CP`$ odd Higgs field $`A^0`$;

$``$ charged fields $`H^\pm `$.

The Higgs boson interaction Lagrangian has the form

$$_{Int}^{Higgs}=_{Int}^{(3)}+_{Int}^{(4)}$$
(1)
$`_{Int}^{(3)}`$ Lagrangian of the triple Higgs boson interactions,
$`_{Int}^{(4)}`$ Lagrangian of the quartic Higgs boson interactions.

$`_{Int}^{(3)}={\displaystyle \frac{\lambda _{hhh}}{3!}}hhh+{\displaystyle \frac{\lambda _{hhH}}{2!}}hhH+{\displaystyle \frac{\lambda _{hHH}}{2!}}hHH+`$
$`{\displaystyle \frac{\lambda _{HHH}}{3!}}HHH+{\displaystyle \frac{\lambda _{hAA}}{2!}}hAA+{\displaystyle \frac{\lambda _{HAA}}{2!}}HAA+`$
$`\lambda _{hH^+H^{}}hH^+H^{}+\lambda _{HH^+H^{}}HH^+H^{}`$ (2)
The tree level Higgs boson self-couplings are represented in (2) as functions of the two free parameters: mass of Higgs field $`A^0`$, $`M_A`$ and mixing angle $`\beta `$.

$`\lambda _{hhh}`$ $`=`$ $`3ia_0\mathrm{cos}(2\alpha )\mathrm{sin}(\alpha +\beta )`$ (3)
$`\lambda _{hhH}`$ $`=`$ $`ia_0[2\mathrm{sin}(2\alpha )\mathrm{sin}(\alpha +\beta )\mathrm{cos}(2\alpha )\mathrm{cos}(\alpha +\beta )]`$ (4)
$`\lambda _{hHH}`$ $`=`$ $`ia_0[2\mathrm{sin}(2\alpha )\mathrm{cos}(\alpha +\beta )+\mathrm{cos}(2\alpha )\mathrm{sin}(\alpha +\beta )]`$ (5)
$`\lambda _{HHH}`$ $`=`$ $`3ia_0\mathrm{cos}(2\alpha )\mathrm{cos}(\alpha +\beta )`$ (6)
$`\lambda _{hAA}`$ $`=`$ $`ia_0\mathrm{cos}(2\beta )\mathrm{sin}(\alpha +\beta )`$ (7)
$`\lambda _{HAA}`$ $`=`$ $`ia_0\mathrm{cos}(2\beta )\mathrm{cos}(\alpha +\beta )`$ (8)
$`\lambda _{hH^+H^{}}`$ $`=`$ $`igM_W\mathrm{sin}(\beta \alpha )ia_0\mathrm{cos}(2\beta )\mathrm{sin}(\alpha +\beta )`$ (9)
$`\lambda _{HH^+H^{}}`$ $`=`$ $`igM_W\mathrm{cos}(\beta \alpha )ia_0\mathrm{cos}(2\beta )\mathrm{cos}(\alpha +\beta )`$ (10)
where $`M_W`$, $`M_Z`$ \- gauge boson masses, $`\theta `$ \- Weinberg angle, $`g`$ \- $`SU(2)_L`$ gauge constant. Connections between the mixing angles $`\alpha `$ and $`\beta `$ and also the parameter $`a_0`$ are determined by the conditions

$`\mathrm{tan}2\alpha =\mathrm{tan}2\beta {\displaystyle \frac{M_A^2+M_Z^2}{M_A^2M_Z^2}},a_0={\displaystyle \frac{igM_z}{2\mathrm{cos}(\theta _W)}}`$ (11)
## 2 Perturbation theory. Vertex function

Main object of our study is the vertex function (VF) of the triple Higgs interaction.

Following the standard normalization conventions , the amplitude of the given diagram can be represented as follows (see Fig. 1):

$$iA_j=iG_1(p_1,\{a_1\})(i\mathrm{\Gamma }_{[3]}(p_1,p_2,p_3,\{c\}))G_2(p_2,\{a_2\})G_3(p_3,\{a_3\})$$
(12)
where $`G_i(p_i,\{a_i\})`$ – the Green’s functions beyond the given diagram, $`p_i`$ – the four-vector of $`ith`$ Higgs boson. $`\{a\},\{c\}`$ – sets of model parameters, characterizing given diagram. $`\mathrm{\Gamma }_{[3]}(p_1,p_2,p_3,\{c\})`$ – vertex function of the triple Higgs interaction. It can be represented as

$$\mathrm{\Gamma }_{[3]}(p_1,p_2,p_3,\{c\})=\underset{l=0}{\overset{\mathrm{}}{}}\mathrm{\Gamma }_{[3]}^{(l)}(p_1,p_2,p_3,\{c\})$$
(13)
In the last expression the decomposition in $`\alpha _e`$ is performed. Each term coincides with the sum of one-particle irreducible diagrams of the order of $`\alpha _e`$. Our analysis is carried out for the decomposition

$$\mathrm{\Gamma }_{[3]}(p_1,p_2,p_3,\{c\})=\mathrm{\Gamma }_{[3]}^{(0)}(p_1,p_2,p_3,\{c\})+\mathrm{\Gamma }_{[3]}^{(1)}(p_1,p_2,p_3,\{c\})$$
(14)
Graphical interpretation of the last formula is shown in Fig. 2.

The condition (13) in terms of self-couplings has the form

$$\lambda _{Hg_1Hg_2Hg_3}=\lambda _{Hg_1Hg_2Hg_3}^{(0)}+\mathrm{\Delta }\lambda _{Hg_1Hg_2Hg_3}^{(1)}=\lambda _{Hg_1Hg_2Hg_3}^{(0)}(1+\alpha _e\mathrm{\Delta }\stackrel{~}{\lambda }_{Hg_1Hg_2Hg_3}^{(1)})$$
(15)
The complete one-loop contribution (Feynman gauge) is determined by contributions of physical fields (standard fermions, sfermions, gauge bosons, chargino, neutralino), as well as Goldstone modes and gost fields.

$`\mathrm{\Delta }\lambda _{Hg_1Hg_2Hg_3}^{(1)}=\mathrm{\Delta }\lambda ^{(1)}(f)+\mathrm{\Delta }\lambda ^{(1)}(\stackrel{~}{f})+\mathrm{\Delta }\lambda ^{(1)}(Gb)+\mathrm{\Delta }\lambda ^{(1)}(Ch)+`$
$`\mathrm{\Delta }\lambda ^{(1)}(Neu)+\mathrm{\Delta }\lambda ^{(1)}(Hg)+\mathrm{\Delta }\lambda ^{(1)}(Gs)+\mathrm{\Delta }\lambda ^{(1)}(Gh)`$ (16)
The following analysis is based on the general formulas for the one-loop contributions to two- and three-point vertex functions. Our algorithms are implemented in the computer algebra programs that allow to operate with bulky intermediate expansions.

## 3 Scheme of calculations

Main features of our calculation are

1. The center-of-mass frame is used for explicit phase space formulas, when one Higgs boson is a virtual particle.

2. Feynman gauge was used.

3. The solutions of renormalization group equations for gauge constants and third generation quarks masses were employed to take into account the energy dependence more precisely .

4. Tensor reduction for the scalar one-loop integrals was used .

5. On-shell renormalization procedure was employed .

The experimental data were taken from .

## 4 Results and their analysis

Two types of dependences are represented in this paper: Higgs self-couplings dependences on $`\mathrm{tan}\beta `$, $`\sqrt{s}`$.

In Fig. 3 three kinds of dependences on $`\mathrm{tan}\beta `$ are represented: tree-level self-couplings ($`\lambda ^0`$), self-couplings with leading $`t\stackrel{~}{t}`$ one-loop correction ($`\lambda ^0+\mathrm{\Delta }\lambda ^{Lead}`$) and self-couplings with complete one-loop contribution ($`\lambda ^0+\mathrm{\Delta }\lambda ^{Comp}`$). Analytic expressions for ($`\mathrm{\Delta }\lambda ^{Lead}`$) can be found in  in the case of large mass approximation for virtual particles.

Obviously, differences between ($`\lambda _{hhh}^0+\mathrm{\Delta }\lambda _{hhh}^{Lead}`$) and ($`\lambda _{hhh}^0+\mathrm{\Delta }\lambda _{hhh}^{Comp}`$) are insignificant (9.5 % from curve ($`\lambda _{hhh}^0+\mathrm{\Delta }\lambda _{hhh}^{Lead}`$) for $`\mathrm{tan}\beta =45`$). It is important to note that the sizable value of correction does not contradict the perturbation theory applicability.

In  it was shown that the self-coupling ($`\lambda _{hhh}^0+\mathrm{\Delta }\lambda _{hhh}^{Comp}`$) in the decoupling limit (our situation for $`M_A=500`$ GeV satisfy this limit) can be represented in terms of Higgs boson $`h^0`$ mass. By redefinition of ($`\lambda _{hhh}^0+\mathrm{\Delta }\lambda _{hhh}^{Lead}`$) the correction value is being insignificant and do not exceed $`6\%`$. The analogous conditions for ($`\lambda _{hhh}^0+\mathrm{\Delta }\lambda _{hhh}^{Comp}`$) give the correction of $`11\%`$. Another situation is observed in $`\lambda _{HHH}`$ case. The curve coinciding to ($`\lambda _{HHH}^0+\mathrm{\Delta }\lambda _{HHH}^{Lead}`$) demonstrates huge one-loop leading contribution in $`\lambda _{HHH}`$, however the complete one-loop contribution is very small. The reason of this discrepancy is in the usage of the large masse approximation for virtual particles. We have not used any approximations for one-loop scalar integrals. The applicability of the large masse approximation is from our point of view questionable in the situation under consideration (maybe except the $`hhh`$ \- case) because of too large masses of virtual Higgs bosons.

In Fig. 4 the self-couplings dependences on the process energy passing by virtual Higgs boson are represented. One can observe that $`\lambda _{HHH}`$ is very weekly dependent on $`\sqrt{s}`$. For a given situation a future high-luminosity colliders can not detect these dependences because we can certainly suppose that $`\lambda _{HHH}`$ is not a running parameter. Another situation is observed in $`\lambda _{hhh}`$ case. For instance for $`\sqrt{s}=186`$ GeV the self-coupling $`\lambda _{hhh}=205`$ GeV and for $`\sqrt{s}=1500`$ GeV the self-coupling $`\lambda _{hhh}=175`$ GeV. If the accuracy of the self-coupling measurement on the future colliders will be better than 30 GeV, that given dependence could be detected.

The summing error of results is determined by determination error of $`t`$ quark mass ($`\delta m=5.1`$ GeV) and by generality of soft supersymmetry breaking parameters values. Errors for $`\lambda _{hhh}`$ case are equal 5 GeV and 8 GeV, and for $`\lambda _{HHH}`$ case 1 GeV and 2 GeV accordingly (for $`M_A=500`$ GeV).

## Conclusion

In this work the neutral Higgs boson self-couplings with complete set of the one-loop diagrams are represented and analyzed in the framework of the MSSM. The effective dependences of the Higgs boson self-couplings $`\lambda _{hhh}`$, $`\lambda _{HHH}`$ on $`\mathrm{tan}\beta `$ and energy $`\sqrt{s}`$ are evaluated with one of the Higgs bosons is off-shell. The authors have avoided any approximations in one-loop scalar integrals calculations. It has been shown, that the tree-level self-couplings can acquire sizable one-loop corrections. These must be taken into account for detailed comparative analysis of theory and experimental data, which could be phenomenologically important.

## Acknowledgements

The authors thank E.E. Boos, M.N. Dubinin and A.A. Biryukov for the valuable discussions and useful critical remarks.


text:
# Gap soliton formation by nonlinear supratransmission in Bragg media

## Introduction.

Light propagation in dielectric media with periodically varying index acquires intriguing properties when Kerr nonlinearities come into play. From a linear theory Brill  we know that such a medium is a photonic band gap structure, or Bragg medium, which is totally reflective when the frequency of the incoming wave falls in one of the gaps.

However, for sufficiently energetic irradiation, the Kerr effect becomes sensible and the Bragg medium may become transparent in the gap, switching from total reflection to high transmissivity, which has been predicted in 1979 winful  by using the stationary coupled mode approach of kog , and which has been experimentally realized for the first time in 1992 sankey .

In 1987 the existence of effective wave propagation in a gap has been associated to the existence of gap solitons derived from slowly varying envelope limits of the Maxwell equation in a periodic structure chen ; mills . For unidirectionnal propagation, fast Kerr nonlinearity and slowly varying evelope approximation, one of these limit models is the coupled mode system desterke  which describes the evolution of the envelopes of the contrapropagating electric fields. A rigorous derivation of from the anharmonic Maxwell-Lorentz equations by the method of multiple scales is given in goodman . The explicit soliton-like solution of the coupled mode equation has been discovered in 1989 christo , and at the same time generalized to a two-parameter family in aceves . This two-parameter gap soliton solution, that we use here, has motivated many experimental searches of the Bragg soliton.

After the pioneering experiments of sankey  in steady-state regime, the first experimental observations of nonlinear pulse propagation in a fiber Bragg grating have been performed in 1996 eggleton  under laser pulse irradiation at a frequency near (but not inside) the photonic band gap. In that case the nonlinearity acts as a source of modulationnal instability generating trains of pulses eggleton2 . Particularly interesting is the observation of pulse propagation at a velocity less than the light velocity. Later in 1998 taverner  a convincing experiment with a fiber Bragg grating demonstrated the repeated formation of gap soliton under quasi-constant wave irradiation inside a band gap, this is the problem we are interested in.

There has been a great deal of discussions concerning the exact role of the gap soliton in the switching to high transmissivity, see e.g. souk . A rather natural intuition is that the nonlinearity shifts the gap and hence allows the medium to become transparent and to create a gap soliton. We shall discover that this intuition is wrong as the process at the origin of high transmissivity in Bragg media is the nonlinear supratransmission recently discovered in the sine-Gordon chain jg-alex ; nst-prl  and experimentally realized (together with application to Josephson junctions arrays) in nst-jpc . As a consequence, the switching to high transmissivity is indeed accomplished by means of gap soliton generation as a result of a fundamental instability instab  which are effectively the objects experimentally detected in the output in taverner .

This result allows us to predict analytically by formula (8) the threshold of incident energy of a cw-radiation above which high transmissivity is reached, as an explicit simple function of the departure (denoted by the dimensionless angular frequency $`\mathrm{\Omega }`$) from the Bragg frequency (gap center). In particular we obain that the process requires less energy for frequencies close to the band edge $`\mathrm{\Omega }=1`$, a result which was previously attributed to pulse shape properties. Interestingly enough, in the limit $`\mathrm{\Omega }1`$ the required energy vanishes, opening the way to experimental realizations.

## The model and its solution.

Our starting point is the basic coupled mode system desterke ; goodman  governing the forward $`E_f`$ and backward $`E_b`$ slowly varying envelopes of the electric field

$$E(Z,T)=[E_fe^{ik_0Z}+E_be^{ik_0Z}]e^{i\omega _0T}.$$
(1)
$`E(Z,T)`$ is the transverse polarized component propagating in the direction $`Z`$ at frequency close to the Bragg frequency $`\omega _0`$. Following goodman  we write the coupled mode system in reduced units as

$`i[{\displaystyle \frac{e}{t}}+{\displaystyle \frac{e}{z}}]+f+({\displaystyle \frac{1}{2}}|e|^2+|f|^2)e=0,`$
$`i[{\displaystyle \frac{f}{t}}{\displaystyle \frac{f}{z}}]+e+({\displaystyle \frac{1}{2}}|f|^2+|e|^2)f=0.`$ (2)
The reduced units are $`z=\kappa Z`$ and $`t=\kappa cT`$ where $`c`$ is the light velocity in the medium and $`\kappa `$ the coupling constant. The reduced field variables are $`e=E_f\sqrt{2\mathrm{\Gamma }/\kappa }`$ and $`f=E_b\sqrt{2\mathrm{\Gamma }/\kappa }`$ where $`\mathrm{\Gamma }`$ is the nonlinear factor. The (linear) dispersion relation of (The model and its solution.), namely $`\omega ^2=1+k^2`$, possess the gap $`[1,+1]`$.

The solitary wave solution of (The model and its solution.) given in aceves  has a simple stationnary expression

$`e=\sqrt{2/3}\lambda e^{i(\mathrm{\Omega }t\varphi )}\mathrm{cosh}^1[\lambda (zz_0){\displaystyle \frac{i}{2}}q],`$ (3)
$`f=\sqrt{2/3}\lambda e^{i(\mathrm{\Omega }t\varphi )}\mathrm{cosh}^1[\lambda (zz_0)+{\displaystyle \frac{i}{2}}q].`$ (4)
The real valued parameters determining this ‘gap soliton’ are $`\mathrm{\Omega }`$ (frequency), $`\varphi `$ (initial phase) and $`z_0`$ (center), and we have

$$\mathrm{\Omega }^2=1\lambda ^2,\mathrm{\Omega }=\mathrm{cos}q.$$
(5)
Then $`\lambda `$ plays the role of the ‘wave number’ of the evanescent wave.

## Nonlinear supratransmission threshold.

For a Bragg medium, extending in the region $`z[0,L]`$, and initially in the dark, we set the initial data

$$e(z,0)=0,f(z,0)=0.$$
(6)
The boundary value problem that mimics the scattering of an incident radiation at frequency $`\omega _0+\kappa c\mathrm{\Omega }`$ on the medium in $`z=0`$ is

$$e(0,t)=Ae^{i\mathrm{\Omega }t},f(L,t)=0,$$
(7)
where the second requirement means no backward wave incident from the right in $`z=L`$. The constant $`A`$ (in general complex valued) is the dimensionless amplitude of the incoming radiation. Equations (6) and (7) constitute a well posed initial-boundary value problem for the PDE (The model and its solution.).

In the linear case such a boundary forcing would simply generate the evanescent wave $`e(z,t)=Ae^{i\mathrm{\Omega }t\lambda z}`$. In the nonlinear case however this boundary forcing generates the solution (3) for the value of $`\varphi `$ and $`z_0<0`$ such that the amplitude in $`z=0`$ be precisely $`A`$. Then, for each fixed forcing frequency $`\mathrm{\Omega }`$, there exists a maximum value $`A_s`$ of $`A`$ beyond which there is no solution $`\{\varphi ,z_0\}`$. This threshold $`A_s`$ is given by $`|e(0,t)|`$ from (3) evaluated at $`z_0=0`$, namely by

$$A_s=2\sqrt{2/3}\mathrm{sin}(\frac{1}{2}\mathrm{arccos}\mathrm{\Omega }).$$
(8)
This is the threshold amplitude of the incident envelope above which the system develops an instability and generates a propagating nonlinear mode.

## Nonlinear dispersion relation prediction.

By seeking a solution of (The model and its solution.) $`e=a\mathrm{exp}[i(kz\omega t)]`$, $`f=b\mathrm{exp}[i(kz\omega t)]`$ one gets a nonlinear algebraic system for the unknowns $`a`$ (incident amplitude) and $`\omega `$ (incident frequency) expressed in terms of the wave number $`k`$ once the amplitude $`b`$ has been eliminated. This gives the nonlinear dispersion relation $`\omega (a,k)`$ as a solution of a third order algebraic equation. One of the 3 solutions must be discarded as it is singular in the linear limit $`a0`$, the other two giving the deformations of the linear branches $`\omega (k)=\pm \sqrt{1+k^2}`$. It can be shown that they correspond to the relations $`b=a`$ (upper branch) and $`b=a`$ (lower branch).

Being interested in the shift of the gap, it is sufficient to study the solutions at $`k=0`$ where the system simplifies, and to stick with the upper branch (as we look at right-going incident waves) for which $`b=a`$. Then (The model and its solution.) readily provides the value $`\omega (0,a)=1\frac{3}{2}a^2`$ of the gap opening.

For our purpose it is more convenient to express the value $`A_m`$ of the incident amplitude $`a`$ for which the incident frequency $`\mathrm{\Omega }`$ touches the nonlinear gap edge $`\omega (0,a)`$, namely

$$A_m=\sqrt{\frac{2}{3}(1\mathrm{\Omega })}.$$
(9)
This formula provides a prediction of transparency by nonlinear shift of the gap. It is now compared to our prediction (8) by means of numerical simulations.

## Numerical simulations.

The system (The model and its solution.) is solved by means of an semi-implicit third order finite difference scheme and boundary values are taken into account at first order. Explicitely we first rewrite (The model and its solution.) for two functions $`u(z,t)`$ and $`v(z,t)`$ and replace the operators $`_t+_z`$ by the set of differences for $`u_n(t)=u(nh,t)`$

$`\dot{u}_1+{\displaystyle \frac{1}{h}}(u_1u_0),\dot{u}_2+{\displaystyle \frac{1}{2h}}(u_2u_0),\mathrm{}`$
$`\dot{u}_n+{\displaystyle \frac{1}{h}}\left[{\displaystyle \frac{2}{3}}(u_{n+1}u_{n1}){\displaystyle \frac{1}{12}}(u_{n+2}u_{n2})\right],\mathrm{}`$
$`\dot{u}_{N1}+{\displaystyle \frac{1}{2h}}(u_{N+1}u_{N1}),\dot{u}_N+{\displaystyle \frac{1}{h}}(u_Nu_{N1}),`$
so as for $`v_n(t)=v(nh,t)`$. The length is $`L=hN`$ and overdot means time differentiation. Equation (The model and its solution.) results as a system of $`2N`$ coupled ODE then solved through the subroutine dsolve of the MAPLE8 software package that uses a Fehlberg fourth/fifth order Runge-Kutta method. Finally the solution, e.g. $`e(z,t)`$, is obtained from $`u_n(t)`$ by $`e(z,t)=[u_{n+1}(t)+u_n(t)]/2`$. Although such code sends some numerical noise in the solution, it is uneffective for sufficiently small $`h`$ depending on the required time of integration.

Next, in order to avoid initial shock, the boundary condition $`e(0,t)`$ is smoothly turned on and smoothly turned off by assuming instead of (7)

$$e(0,t)=\frac{A}{2}\left[\mathrm{tanh}(p(tt_0))\mathrm{tanh}(p(tt_1))\right].$$
(10)
A practical interest of such an incident wave is that it reproduces the quasi-constant wave irradiation of the experiments of taverner .

Most of the results presented here are obtained with $`N=120`$ spatial mesh points, $`h=0.05`$ grid spacing over (hence a length $`L=6`$ normalize units) a time of integration $`t_m=200`$. The parameters of the boundary field are a real-valued amplitude $`A`$, a slope $`p=0.2`$, an ignition time $`t_0=20`$ and an extiction time $`t_1=180`$.

We display in figure 1 an instance of two different numerical solutions (the modulus of the right-going envelope) for an incident frequency $`\mathrm{\Omega }=0.8`$ and amplitudes $`A=0.51`$ (no transmission) and $`A=0.52`$ (gap soliton generation) when the theoretical threshold predicted by (8) is $`A_s=0.5164`$.

Using this simple diagnostic, the bifurcation predicted by expression (8) is numerically checked for a series of frequency values (in the range $`[0.1,0.995]`$) and we obtain the figure 2 where the dots represent the smallest value (with absolute precision of $`10^2`$) of the amplitude $`A`$ for which nonlinear supratransmission is seen to occur. The expression (9) is also plotted (dashed line) which shows that the nonlinear shift provides a wrong answer.

## Bifurcation of transmitted energy.

The system (The model and its solution.) possess the conservation law

$$_t(|e|^2+|f|^2)+_z(|e|^2|f|^2)=0.$$
(11)
For a given boundary condition (7) (i.e. for fixed $`A`$ and $`\mathrm{\Omega }`$), we define the incident ($`I`$), reflected ($`R`$) and transmitted ($`T`$) energies at given arbitrary time $`t_m`$ by

$`I(A,\mathrm{\Omega })={\displaystyle _0^{t_m}}𝑑t|e(0,t)|^2,`$ (12)
$`R(A,\mathrm{\Omega })={\displaystyle _0^{t_m}}𝑑t|f(0,t)|^2,`$ (13)
$`T(A,\mathrm{\Omega })={\displaystyle _0^{t_m}}𝑑t|e(L,t)|^2.`$ (14)
Then, upon integration of (11) on the length $`[0,L]`$ and time $`[0,t_m]`$ we obtain

$$R+TI+_0^L𝑑z\left(|e(z,t_m)|^2+|f(z,t_m)|^2\right)=0.$$
(15)
If $`t_m`$ is larger than the irriadiation duration, the energy injected eventually radiates out completely, namely $`e(z,t_m)`$ and $`f(z,t_m)`$ vanish, and we are left with the conservation relation $`R+T=I`$ which can be written

$$\rho (A,\mathrm{\Omega })+\tau (A,\mathrm{\Omega })=1,\rho =R/I,\tau =T/I.$$
(16)
The figure 3 shows a typical numerical simulation where, for $`\mathrm{\Omega }=0.95`$, we have computed the reflection and transmission factors, together with their sum, for 50 different values of the amplitude $`A`$ in the range $`[0.20,\mathrm{\hspace{0.17em}0.45}]`$.

These simulations showt the sudden energy flow through the medium, as a result of nonlinear supratransmission. Note from the sum $`\rho +\tau `$ before the bifurcation on figure 3 that the numerical code conserves the energy with a precision of $`2\%`$.

The transmissivity is due to the generation and propagation of light pulses shown figure 4 representing the energy density $`|e(L,t)|^2`$ measured at the output as a function of time. We have also plotted the input energy density given in (10) with parameters $`t_m=100`$, $`t_0=20`$, $`t_1=90`$ and $`p=2`$.

We show now that these are gap solitons travelling at a fraction of the light velocity (slow light pulses).

## Travelling gap solitons.

Although the system (The model and its solution.) is not Lorentz invariant, the propagating solution of aceves  can still be written in terms of the boosted variables barash . As we are intersted here only in the energy flux, we write the soliton solution $`|e(z,t)|^2`$ moving at velocity $`v<1`$ with frequency $`\mathrm{cos}q`$ as

$`|e(z,t)|^2={\displaystyle \frac{2\mathrm{sin}^2q\left(\frac{1v^2}{3v^2}\right)\left(\frac{1+v}{1v}\right)^{1/2}}{\left|\mathrm{cosh}\left(\frac{\mathrm{sin}q}{\sqrt{1v^2}}(zz_0vt)\frac{i}{2}q\right)\right|^2}}.`$ (17)
Such an expression allows to fit a given simulation by seeking the two parameters $`v`$ and $`q`$, and the initial position $`z_0`$, that reproduce with the explicit solution the results of the numerical simulation.

An instance of such a fit is displayed in figure 5 that shows the analytic soliton $`|e(z,t)|^2`$ compared with the numerical simulation for two fixed values of time. This furnishes the following velocities

| $`\mathrm{\Omega }`$ | 0.95 | 0.9 | 0.8 | 0.7 | 0.6 | 0.5 |
| --- | --- | --- | --- | --- | --- | --- |
| $`A`$ | 0.28 | 0.38 | 0.53 | 0.65 | 0.75 | 0.83 |
| $`v`$ | 0.25 | 0.34 | 0.42 | 0.50 | 0.58 | 0.66 |

(18)
These are slow light pulses and naturally, small velocities are obtained close to the gap edge where the required energy is small.

## Conclusion.

The property of a Bragg medium in the nonlinear Kerr regime to sustain nonlinear supratransmission has allowed us to obtain the analytic expression (8) that fixes the threshold amplitude of an incident cw-beam in terms of its departure from the Bragg resonance.

This constitutes a practical tool to investigate switching properties of a Bragg medium and allows at the same time to understand the process at the origin of sudden transmissivity of the Bragg mirror. It is the nonlinear supratansmission that results from a nonlinear instability intrinsic to boundary value problems instab .

We expect our result to be usefull for experiments by the ability of the Bragg medium to become partly transparent by gap soliton generation for an incident beam of low energy if its frequency is chosen close to the gap edge.

## Acknowledgments.

This study was initiated during a stay of one of us (J.L.) at the department of applied mathematics, Boulder University. It is a pleasure to aknowledge invitation and enlighting discussions with M.J. Ablowitz.

This work received support from Fundação para à Ciência e à Tecnologia, grant BPD/5569/2001, and from Açoes Integradas Luso-Francesas, grant F-4/03.