text:
# 1 Introduction

## 1 Introduction

The study of $`g`$-factors of subatomic particles can trace its roots back to the 1921 paper by Stern which was a proposal to study space quantization with an apparatus which is now called a “Stern-Gerlach apparatus”. By 1924 the famous experiments had been done and a review paper was written summarizing their results. Their final conclusion, that “to within 10% the magnetic moment of the electron was one Bohr magneton”, meant in modern language that the $`g`$-value of the electron was 2, where the gyromagnetic ratio $`g`$ is the proportionality constant between the magnetic moment and the spin,

$$\stackrel{}{\mu }=g(\frac{e}{2m})\stackrel{}{s}.$$
(1)
The discovery that $`g_e2`$, and the calculation by Schwinger predicting that (to first order) the radiative correction to $`g_e`$ was $`\alpha /\pi `$, were important early steps in the development of Quantum Electrodynamics (see Fig. 1).

The long lifetime of the muon permits a precision measurement of its anomalous moment at the ppm level. The muon magnetic moment is given by

$$\mu _\mu =(1+a_\mu )\frac{e\mathrm{}}{2m_\mu }\mathrm{where}a_\mu =\frac{(g2)}{2}.$$
(2)
The electron anomalous magnetic moment has been measured to a few parts per billion, which can be completely described by QED of electrons and photons to eighth order, $`(\frac{\alpha }{\pi })^4`$. The contributions of virtual muons, tauons, etc. enter at the few ppb level. The calculation of the electron anomalous moment is limited by the knowledge of the fine-structure constant. With the reliability of modern QED calculations, Kinoshita has turned things around and has used the electron $`g`$ value measurement to give the best value for $`\alpha `$.

The relative contribution of heavier particles to the muon anomaly scales as $`(m_\mu /m_e)^2`$ and the famous CERN experiment, which obtained a relative error on $`a_\mu `$ of $`\pm 7.3`$ parts per million (ppm), easily observed the predicted $`60`$ ppm contribution of virtual hadrons.

In 1984 efforts began to make a new measurement of the muon anomalous moment to a precision of $`\pm 0.35`$ ppm, which would represent a 5 standard deviation observation of the electroweak contribution, and would also be sensitive to contributions from “new physics” such as muon substructure or supersymmetry.

## 2 Theoretical Contributions to $`(g2)`$

The standard model value of $`a_\mu `$ is given by $`a_\mu (\mathrm{SM})=a_\mu (\mathrm{QED})+a_\mu (\mathrm{hadronic})+a_\mu (\mathrm{weak})`$ and any contribution from new physics would be reflected in a measured value which did not agree with the standard model.

Comparison of the measurements and calculations of the electron $`g`$ value gives one great confidence in our understanding of QED to the level needed for muon $`(g2)`$. Taking the value of $`\alpha `$ from the electron $`(g2)`$, yields the total QED contribution $`a_\mu (\mathrm{QED})=\mathrm{116\; 584\; 705.7}(1.8)(0.5)\times 10^{11}`$.

The hadronic contribution to $`(g2)`$ cannot be calculated directly, but must be determined from data. The first-order hadronic vacuum polarization dominates the uncertainty in the theoretical value of $`a_\mu `$, since it is calculated using dispersion theory and data from $`e^+e^{}\mathrm{hadrons}`$ and hadronic $`\tau `$ decay as input. The various order hadronic contributions are shown in Fig. 2. Diagrams for hadroproduction and hadronic $`\tau `$ decay are shown in Fig. 3. The most precise determination of the first-order hadronic contribution is $`a_\mu (\mathrm{had};1)=6924(62)\times 10^{11}`$ which is $`59.39\pm 0.53`$ ppm of $`a_\mu `$, but there is continuing discussion of the use of CVC and the $`\tau `$-decay data. The higher-order contribution is $`a_\mu (\mathrm{had};2)=101(6)\times 10^{11}`$ The hadronic light-by-light scattering shown in Fig. 4 has now been calculated by two groups, using essentially the same model, and agreement is found: $`a_\mu (\mathrm{had};\mathrm{lbl})=85(32)10^{11}`$. However the two groups disagree on the uncertainty on the calculation, and I have taken the larger error from Ref. . The uncertainty in this contribution could be reduced substantially by the appropriate calculation on the lattice, and perhaps by other additional calculations as well.

The total hadronic contribution is given by $`a_\mu (\mathrm{had};1+2+\mathrm{lbl})=6738(70)\times 10^{11}`$ which is $`57.79\pm 0.60`$ ppm of $`a_\mu `$, with an uncertainty dominated by the uncertainty on the first-order hadronic vacuum polarization. It is precisely this contribution which is being addressed by the programs to measure $`R(s)`$ at BES and the Budker Institute. We look forward to additional high quality data from these experiments, from DAPHNE, as well as $`\tau `$-decay data from CLEO to further reduce the uncertainty on the hadronic contribution.

The standard model electroweak contribution arises from the diagrams shown in Fig. 5 (the standard model Higgs contribution is negligible). The single-loop $`W`$ and $`Z`$ contributions were calculated by a number of authors shortly after the standard model was developed. The result is $`a_\mu (\mathrm{weak};1)=195\times 10^{11}`$ or 1.7 ppm of $`a_\mu `$. Partial calculations of the two-loop electroweak contributions indicated that they might not be small. The full calculation which was later confirmed independently showed that the total first and second-order weak contribution was 20% less than the first order result. The result is $`a_\mu (\mathrm{weak};1+2)=151(4)\times 10^{11}`$ which is $`1.30\pm 0.03`$ ppm of $`a_\mu `$.

The standard model prediction for $`a_\mu `$ is $`a_\mu (\mathrm{SM})=(\mathrm{116\; 591\; 594.7}\pm 70)\times 10^{11}`$ ($`\pm 0.60`$ ppm).

A great deal has been written about the possible contribution to the muon $`(g2)`$ value from non-standard model physics. Just as the proton substructure produces a $`g`$-value which is not equal to two, muon substructure would also contribute to the anomalous moment, the critical issue being the scale of the substructure. A standard model value for $`(g2)`$ at the 0.35 ppm level would restrict the substructure scale to around 5 TeV.

In Fig. 5(a) the triple gauge vertex $`WW\gamma `$ appears, and it is through this diagram that the muon $`(g2)`$ obtains its sensitivity to $`W`$ substructure and anomalous gauge couplings. The combined sensitivity of LEP1, LEP2 and $`(g2)`$, and the unique contribution which $`(g2)`$ makes in constraining the existence of such couplings, is described by Renard et al.

Supersymmetry has become a serious candidate for physics beyond the standard model. The SUSY contribution is shown in Fig. 6. In the case of large $`\mathrm{tan}\beta `$, the chargino diagram dominates and the contribution to $`(g2)`$ from SUSY is given by

$$a_\mu (\mathrm{SUSY})\frac{\alpha }{8\pi \mathrm{sin}^2\theta _W}\frac{m_\mu ^2}{\stackrel{~}{m}^2}\mathrm{tan}\beta 140\times 10^{11}(\frac{100\mathrm{GeV}}{\stackrel{~}{m}})^2\mathrm{tan}\beta .$$
(3)
The goal of E821 is to reach a precision of $`\pm 40\times 10^{11}`$ ($`\pm 0.35`$ ppm), so the factor of 140 above corresponds to 1.2 ppm. For $`\stackrel{~}{m}=750`$ GeV and $`\mathrm{tan}\beta =40`$, $`a_\mu (\mathrm{SUSY})=100\times 10^{11}`$.

## 3 The New $`(g2)`$ Experiment

For polarized muons moving in a uniform magnetic field $`\stackrel{}{B}`$ which is perpendicular to the muon spin direction and to the plane of the orbit, and with an electric quadrupole field $`\stackrel{}{E}`$ for vertical focusing, the difference angular frequency, $`\omega _a`$, between the spin precession frequency $`\omega _s`$ and the cyclotron frequency $`\omega _c`$, is given by

$$\stackrel{}{\omega }_a=\frac{e}{m}\left[a_\mu \stackrel{}{B}\left(a_\mu \frac{1}{\gamma ^21}\right)\stackrel{}{\beta }\times \stackrel{}{E}\right].$$
(4)
The dependence of $`\omega _a`$ on the electric field is eliminated by storing muons with the “magic” $`\gamma _\mu `$=29.3, which corresponds to a muon momentum $`p_\mu `$ = 3.09 GeV/$`c`$. Hence measurement of $`\omega _a`$ and of $`B`$ determines $`a_\mu `$. At the magic gamma, the muon lifetime is $`\gamma \tau =64.4`$ $`\mu `$s, the $`(g2)`$ precession period is 4.37 $`\mu `$s, and for the central orbit radius of 7.11 m the cyclotron period is 149 ns.

The storage ring magnet is a superferric 700 ton, 14 m diameter circular “C”-magnet, with the opening facing inward towards the ring center. The field is excited by three 14 m diameter superconducting coils which carry $`5.2`$ kA from a low voltage power supply to produce the $`1.45`$ T magnetic field. The short term field stability over several AGS cycles is better than 0.1 ppm.

The magnetic field which enters in Eq. 4 is the average field seen by the muon distribution. Since direct injection of muons does not uniformly fill the phase space, we used a tracking code to calculate the distribution of muons in the storage ring. The radial distribution obtained from this tracking code was compared with the distribution obtained from observing the beam debunching in the ring at early times. The two distributions agreed quite well. In the 1999 run, a straw-tube array was operational at one detector location which provided information on the decay positron trajectories coming out of the storage region. These data will permit us to reconstruct directly the muon spatial distribution in one section of the ring.

In 1998, direct muon injection into the storage ring was employed for the first time. The AGS performance, the beamline and the inflector magnet were as described in Carey et al., and except for the muon injection many of the experimental details are the same as described there. The positive muon beam with the magic momentum is formed by collecting the highest energy muons from pion decay in a 72 m long decay section of our beamline, which results in a muon polarization of 96%. The flux incident on the inflector magnet was $`2\times 10^6`$ per fill of the ring.

The 10 mrad kick needed to put the muon beam onto a stable orbit was achieved with a currents, since usual magnetic kicker techniques would spoil the precision magnetic field. Three pulse-forming networks powered three identical 1.7 m long kicker sections consisting of parallel plates on either side of the beam. Current flowed down one side crossed over and flowed back up the other side. The kicker plate geometry and composition was chosen to minimize eddy currents, and the eddy current effect on the total field seen by the muons was less than 0.1 ppm 20 $`\mu `$s after injection. The current pulse, which was formed by an under-damped LCR circuit, had a peak current of 4100 A and a pulse base width of 400 ns. Since the cyclotron period of the muon beam from the AGS was 149 ns, the beam was kicked several times before the kicker pulse died out. With muon injection, the number of detected positrons per hour was increased by an order of magnitude over pion injection. Thus the use of a muon beam an order of magnitude less intense than the pion beam resulted in a substantial increase in stored muons per fill of the ring, with the injection related background reduced by about a factor of 50.

Positrons from the in-flight decay $`\mu ^+e^+\nu _e\overline{\nu }_\mu `$ are detected with Pb-scintillating fiber calorimeters placed symmetrically at 24 positions around the inside of the storage ring. The decay positron time spectrum is

$$N_0e^{t/\gamma \tau }\left[1+A(E)\mathrm{cos}\left(\omega _at+\varphi (E)\right)\right].$$
(5)
The normalization constant $`N_0`$ and the parity violating asymmetry parameter $`A(E)`$ depend on the energy threshold placed on the positrons. The fractional statistical error on $`\omega _a`$ is proportional to $`A^1N_e^{1/2}`$, where $`N_e`$ is the number of decay positrons detected above some energy threshold. For an energy threshold of 1.8 GeV, we measure $`A`$ to be $`0.34`$ consistent with its theoretical value, which we attribute to the good calorimeter energy resolution ($`\sigma /E=10\%`$ at 1 GeV) and a scalloped vacuum chamber which minimizes pre-showering before the positrons reach the calorimeters.

The photomultiplier tubes of the calorimeter were gated off before injection and when gated on, they recovered to 90% pulse height in $`400`$ ns, and reached full operating gain in several $`\mu `$s. With the reduced flash following injection it was possible to begin counting as soon as 5 $`\mu `$s after injection in the region of the ring 180 around from the injection point.

The calorimeter pulses were continuously sampled by custom 400 MHz waveform digitizers (WFDs), which provided both timing and energy information for the positrons. Both the NMR and WFD clocks were phase-locked to the same LORAN-C frequency signal. The waveforms were zero-suppressed, and stored in memory in the WFD until the end of the AGS cycle. Between AGS acceleration cycles the WFD data were written to tape for off-line analysis, as were the calorimeter calibration data and the magnetic field data.

A laser/LED calibration system was used to monitor calorimeter time and gain shifts during the data-collection period. Early-to-late timing shifts over the first 200 $`\mu `$s were on-average less than 20 ps, which is needed to keep systematic timing errors smaller than 0.1 ppm.

For the offline analysis, the detector response (waveform shape) to positrons was determined from our data for each calorimeter. These shapes were then fit to all pulses in the data to determine a time, an amplitude and a width parameter for each pulse. Time histograms were formed for each detector. These independent data sets were analyzed separately and were in agreement ($`\chi ^2/\nu =17.2/20`$).

A completely blind analysis was performed on the data. Arbitrary offsets were put on the muon frequency and the proton frequency from the NMR probes. Each offset was known by one person, making it impossible to determine the actual value of $`a_\mu `$. Only when the analysis of both the magnetic field and $`\omega _a`$ were completed were the offsets removed and the new value of $`a_\mu `$ determined.

After the offsets were removed, it was necessary to make two corrections to the frequency obtained from the fitting. For muons with the “magic” momentum, $`\omega _a`$ is not affected by the electric field. For the ensemble of muons in our storage ring there is a small electric field correction to $`\omega _a`$ since not all muons are at the magic momentum. There is also a pitch correction because of the vertical betatron oscillations. The sum of these two corrections for these data is ($`0.9\pm 0.2`$) ppm.

The dominant systematic errors which were reported in our first measurement have been completely eliminated. The remaining systematic errors are under study. With many of them approximated by upper limits, one obtains a total systematic error of $`1`$ ppm, with the systematic error assigned to the magnetic field of 0.5 ppm. Since the study of systematic errors is a source of much continuing work, we have chosen not to present a detailed list at this time.

We do wish to note the substantial improvement in the magnetic field quality which has been obtained by additional shimming. In Fig. 7 we show the average magnetic field from the 1997 and 1998 runs. The field uniformity over the storage aperture in 1998 was almost an order of magnitude better than was obtained in the CERN experiment.

## 4 New Results

One month after the Symposium, we finished our analysis of the 1998 data and obtained a new result at the precision of $`\pm 5`$ ppm.

Our experiment measures the frequency ratio $`R=\omega _a/\omega _p`$, where $`\omega _p`$ is the free proton NMR frequency in our magnetic field. Including the pitch and electric field corrections, we obtain $`R=\mathrm{3.707\; 201}(19)\times 10^3`$, where the 5 ppm error includes a 1 ppm systematic error estimate. We obtain $`a_{\mu ^+}`$ from $`a_{\mu ^+}=R/(\lambda R)`$ $`=\mathrm{116\; 591\; 91}(59)\times 10^{10}`$ in which $`\lambda =\mu _\mu /\mu _p=\mathrm{3.183\; 345\; 39}(10)`$. This new result is in good agreement with the mean of the CERN measurements for $`a_{\mu ^+}`$ and $`a_\mu ^{}`$, and our previous measurement of $`a_{\mu ^+}`$, which are tabulated below.

Assuming CPT symmetry, the weighted mean of the four measurements gives a new world average of $`a_\mu =\mathrm{116\; 592\; 10}(46)\times 10^{10}`$ $`(\pm 3.9\mathrm{ppm})`$, which agrees with the standard model to within one standard deviation. These results are displayed graphically in Fig. 8. Also shown is the projected error from the 1999 data, and the $`\pm 0.35`$ ppm goal of E821.

## 5 Outlook and Conclusions

The standard model has been remarkably successful in describing a wide range of phenomena including the entire body of results from LEP. The new result from E821, which has a precision of 5 parts per million, lowers the uncertainty on our knowledge of $`(g2)_\mu `$ to 3.9 ppm. At that level there is good agreement with the standard model, but there is still room for the observation of a value at the level of one ppm or better which could agree with the previous measurements, and disagree with the standard model. We expect our one ppm result to be available before the end of 2000.

This experiment is very much a work in progress. The initial design of E821 was made with the systematic error budget of 0.12 ppm. Many of the improvements which were made to the CERN technique have indeed worked, making it straightforward to obtain a systematic error less than one ppm. While we have been learning about our systematic errors from the beginning of the first pion injection run, with the sub-ppm statistics available with muon injection we are now challenged to push the analysis of systematic errors to the limit. Our ultimate goal is a statistical answer at the 0.3 ppm level, with systematic errors at perhaps half this level. Our work thus far seems to indicate that this level of systematic error will be possible.

A new inflector magnet has been installed in the ring, which has a fringe field in the storage region a factor of 5 less than the old inflector had. This will improve our ability to map the field everywhere and will help to reduce the uncertainty on our knowledge of the field. The principal magnetic field issue remaining is our ability to track the field with time, and to understand the full calibration of the NMR probes to a few tenths of a ppm.

The other principal challenge is how to handle pile-up in the detectors without discarding an unreasonable portion of the data set. The dominant analysis effort now is to understand this effect.

It has been 20 years since the CERN experiment presented its final report which verified the hadronic contribution at the eight standard deviation level, but which left it to a future experiment to verify the electroweak contribution. The new Brookhaven experiment is now approaching that goal. Within the next few years either the standard model value will be confirmed, or evidence of a new contribution to the muon $`(g2)`$ will be discovered.

I wish to acknowledge the efforts of the many collaborators who have worked on $`(g2)`$ over the past 15 years. The steady support of the funding agencies (US and abroad) and the Brookhaven Laboratory management was essential to our reaching this point. I wish to thank R. Carey, D. Hertzog, K. Jungmann, V. Hughes, J. Miller, Y. Semertzidis and E. Sichtermann, for their comments on this manuscript.


format: text

text:
# High-resolution simulations and visualization of protoplanetary disks

## 1. The Method

Extremely small temporal and spatial scales involved in the problem of accretion onto a protoplanet necessitate the use of nonuniform discretization in the vicinity of the accretor. In our study we used adaptive mesh refinement (AMR) method combined with a high-resolution Godunov-type advection scheme (amra, Plewa & Müller 2000). The AMR discretization scheme follows the approach of Berger and Colella (1989). The computational domain is covered by a set of completely nested patches occupying levels. The levels create a refinement hierarchy. As one moves toward higher levels, the numerical resolution increases by a prescribed integer factor (separate for every direction). The net flow of material between patches at different levels is carefully accounted for in order to preserve conservation properties of hydrodynamical equations. Boundary data for child patches are either obtained by parabolic two-dimensional conservative interpolation of parental data or set according to prescribed boundary conditions.

Hydrodynamical equations are solved with the help of the Direct Eulerian Piecewise-Parabolic Method (PPMDE) of Colella & Woodward (1984), as implemented in herakles solver (Plewa & Müller 2000). Simulations have been done in spherical polar coordinates in a frame of reference corotating with the protoplanet. herakles guarantees exact conservation of angular momentum which is particularly important in numerical modeling of disk accretion problems. The use of its multifluid option with tracer materials distributed within disk (not presented here) allows to identify the origin of the material accreted onto protoplanet. The amra code is written purely in FORTRAN 77 and has been successfully used on both vector supercomputers and superscalar cache-oriented workstations. Its parallelization on shared memory machines exploits microtasking (through the use of vendor-specific directives) or the OpenMP standard.

## 2. Simulation setup

The computational domain extends from 0.25 to 2.5 radii of the planet’s orbit. We employ 7 levels with refinement ratios ranging from (2,4) to (4,4). The base level contains the protoplanetary (circumstellar) disk while the 7th level contains the planet and its immediate vicinity. The base grid consist of $`128\times 128`$ cells uniformly distributed in $`r`$ and $`\theta `$. The effective resolution at the 7th level is $`131072\times 524288`$ in $`r`$ and $`\theta `$, respectively. The topmost five levels are schematically shown in Figure 1.

White lines are boundaries of the patches. There are 1, 1, 1, 1, 12, 4 and 49 patches at levels 1-7, respectively. The structure of the grid at level 7 is shown in Figure 2f

with individual cell boundaries drawn with white lines (the dark blue circle shows size of the planet).

## 3. Physical model

The simulation is initialized with a Keplerian disk. Originally the disk has a mass of 0.01 M$`_{}`$, constant $`h/r`$ ratio of 0.05 and surface density proportional to $`r^{1/2}`$. The temperature is a fixed function of $`r`$ throughout the simulation. There is no explicit viscosity in the disk. At the outer and inner boundary of the base grid the gas is allowed to flow freely from the computational domain. No inflow is allowed for. The accretion onto the planet is accounted for in a very simplified way. At every time step the mean value of the density within two planetary radii is calculated, and whenever it is higher then a preset value, the excess gas is removed. At $`t=0`$ a planet of one Jupiter mass in inserted into the disk on a circular orbit. The radius of the orbit and the mass of the planet remain constant throughout the simulation. The disk is allowed to evolve for 100 planetary orbits. A gap is cleared in it, and a secondary, circumplanetary disk is formed.

The sequence of surface plots in Figure 2 shows the final structure of both disks (the surface density distribution is displayed). The red peak in Figure 2a is the unresolved image of the very dense circumplanetary disk. We have been able for the first time to see the details of the latter (Figures 2(c)-d). The streams of gas flowing across the gap from left and right edge of the frame (light blue) collide with the outer part of the circumplanetary disk. The collision regions (green wedges) bear strong resemblance to hot spots in cataclysmic binaries. In every region two strong shock waves are excited, one of them propagating into the stream, and the other into the disk. The shocked gas flows from the collision region along a loosely wound spiral towards the planet ((Figure 2e). This picture is significantly more detailed than the one recently published by Lubow, Seibert, & Artymowicz (2000). Streamlines of the flow around the planet are shown in Figure 3,

and they are in good accordance with those of Lubow et al.

Our simulation is of preliminary nature, and its sole purpose is to demonstrate the capabilities of amra. Currently, we are improving the physics of the model. One of the problems we are going to attack is the calculation of the accurate value of the gravitational torque from the disk onto the planet in the phase preceding gap formation.

## 4. Visualization

To visualize the complicated amra output, we have chosen the AVS/Express environment for visual programming. It allows the user to quickly built simple applications employing standard library modules. Advanced users can develop their own, highly specialized modules and applications. Our amra-visualization application (visa) is partly based on modules written by Favre, Walder, & Follini (1999), which have been substantially modified, and partly on our own modules. A screenshot of visa is shown in Figure 4.

The panel and the viewer are contained in the two topmost windows, while the bottom window contains the AVS/Express programming platform. Currently we are able to read the AMR data, extract components, perform mathematical operations on data sets and coordinates, extract any subset of levels or patches, and apply to them various visualization technique (e.g. 2-D plot, surface plot, isolines, slice). Streamlines can also be calculated. The application is still under development, and new options are being added.

### Acknowledgments.

This research is supported by the Polish Committee for Scientific Research through the grant 2.P03D.004.13.

## References

Berger, M. J., & Colella, P. 1989, J. Comput. Phys., 82, 64

Colella, P., & Woodward, P.R. 1984, J. Comput. Phys., 59, 264

Plewa, T., & Müller, E. 2000, Comp. Phys. Commun. (in preparation)

Lubow, S.H, Seibert, M., & Artymowicz, P. 2000, ApJ (astro-ph/9910404)

Favre, J. M., Walder, R., & Follini, D. 1998, in Proceedings, 40th Cray User Group Conference, Stuttgart, Germany (June 1998)


text:
# An 11.6 Micron Keck Search For Exo-Zodiacal Dust

## 1 INTRODUCTION

Our sun is surrounded by a disk of warm ($`>`$150 K) “zodiacal” dust that radiates most of its thermal energy at 10–30 microns. This zodiacal dust is produced largely in the inner part of the solar system by collisions in the asteroid belt (Dermott et al. 1992) and cometary outgassing (Liou and Zook 1996). Zodiacal dust is interesting as a general feature of planetary systems, and as an indicator of the presence of larger bodies which supply it; dust orbiting a few AU from a star is quickly removed as it loses angular momentum to Poynting-Robertson drag (Robertson 1937). Understanding the extra-solar analogs of zodiacal dust may also be crucial in the search for extra-solar planets (Beichman et al. 1996) since exo-zodiacal dust in a planetary system could easily outshine the planets and make them much harder to detect.

The best current upper limits for the existence of exo-zodiacal dust disks come from IRAS measurements of 12 and 25 micron excesses above photospheric emission. Seen from a nearby star, solar system zodiacal dust would create only a $`10^4`$ excess over the sun’s photospheric emission at 20 microns. IRAS measurements, however, have typical measurement errors of 5 percent (Moshir et al. 1992) and display systematic offsets of a similar magnitude when they are compared to other photometry (Cohen et al. 1996). If there were a solar-type zodiacal disk with 1000 times the density of the disk around the sun around Tau Ceti, the nearest G star, the excess infrared emisison would barely exceed the formal 68% confidence intervals of the IRAS photometry. Moreover, all photometric detection schemes of this sort are limited by how accurately the star’s mid-infrared photospheric emission is known. For farther, fainter stars than Tau Ceti, inferring the presence of dust from the IRAS data becomes still harder.

The detection of faint exo-zodiacal-dust emission is more feasible if one can resolve the dust emitting region. The high resolution and dynamic range needed for these observations will generally require large interferometers like the Keck Interferometer, the Large Binocular Telescope, and the Very Large Telescope Interferometer. But it is already possible to resolve the zodiacal dust mid-infrared emitting regions of the nearest stars. A 10-meter telescope operating at 12 microns has a diffraction-limited resolution of 0.25 arc seconds, corresponding, for example, to a transverse distance of 2 AU at 8 parsecs.

We have begun a search for zodiacal dust around the nearest stars using the mid-infrared imaging capabilities of the Long Wavelength Spectrometer (LWS) (Jones & Peutter 1993) on the W. M. Keck telescope. The large aperture of the telescope allows us to make spatially resolved images of the zodiacal dust 11.6 micron emitting region around the stars so that we can look for dust emission above the wings of the point-spread function (PSF) rather than as a tiny photometric excess against the photosphere. We present here the results of two nights of observations, and compare them with a simple model of exo-zodiacal thermal emission to place upper limits on the amount of dust present in the systems we observed.

## 2 OBSERVATIONS

We observed six nearby stars with LWS on the W. M. Keck telescope on August 3rd and 4th, 1996 using standard mid-infrared imaging techniques. The target stars were the nearest A–K main-sequence stars observable from Mauna Kea on those dates. With the object on-axis, we took a series of frames lasting 0.8 ms each, chopping the secondary mirror between the object and blank sky 8 arcseconds to the north at a frequency of 10 Hz. Then we nodded the primary mirror for the next series of frames so that the sky was on-axis and the object off-axis. We repeated this process for 3 nods over a period of 5 minutes, for an on-source integration time of 1.1 minutes, and a typical noise of 2 mJy in one 0.11 by 0.11 arcsecond pixel due to the thermal background. The seeing was poor both nights, up to 2 arc seconds in the visible. To measure the atmosphere-telescope transfer function, we made similar observations of seven distant, luminous calibrator stars near our targets on the sky, alternating between target and calibrator every 5–10 minutes.

We increased our frame-rate for the second night of observations so that we could compensate for the seeing using speckle analysis. Figure 1 shows a cut through a single 84 ms exposure of Altair on August 4, compared to an Airy function representing the diffraction-limited PSF of a filled 10-meter aperture at 11.6 microns. The cores of the images are diffraction-limited, but the wings are sensitive to the instantaneous seeing, making speckle analysis necessary. Table 1 provides a summary of our observations.

We flat-fielded the images by comparing the response of each pixel to the response of a reference pixel near the center of the detector. First we plotted the data number (DN) recorded by a given pixel against the DN in the reference pixel for all the frames in each run. Since the response of each pixel is approximately linear over the dynamic range of our observations and most of the signal is sky background, which varies with time but is uniform across the chip, the plotted points for each pixel describe a straight line; if all the pixels had the same response, the slope of each line would equal 1. We divided each pixel’s DN by the actual slope of its response curve relative to the reference pixel, effectively matching all pixels to the reference pixel. We then interpolated over bad pixels, frame by frame.

To compensate for the differences in the thermal background between the two nod positions, we averaged together all the on-axis sky frames to measure the on-axis thermal background and subtracted this average from all of the on-axis frames—both object and sky. We used the same procedure to correct the off-axis frames.

Next, we chose subframes of 32 by 32 pixels on each image, centered on the star (or for sky frames, the location of the star in an adjacent object frame), and processed these according to classical speckle analysis (Labeyrie 1970). We Fourier transformed them, and summed the power spectra, yielding a sky power spectrum and an object power spectrum for each series. Then we azimuthally averaged the power spectra in the u-v plane—that is, we averaged over all the frequency vectors of a given magnitude, $`\sqrt{u^2+v^2}`$. This azimuthal averaging corrects for the rotation of the focal plane of the alt-az-mounted Keck telescope with respect to the sky. We then subtracted from every object power spectrum the corresponding sky power spectrum and divided each corrected target power spectrum by the corrected power spectrum of a calibrator star observed in the same manner as the target star immediately before or after the target star. Figure 2 shows an azimuthally-averaged power spectrum of Altair and the corresponding sky power spectrum, compared with a power spectrum of calibrator Gamma Aquila and its corresponding sky power.

We then averaged all the calibrated power spectra for a given target. If the object and calibrator are both unresolved, the average calibrated power spectrum should be the power spectrum of the delta function: a constant. We found that the pixels along the u and v axes of the power spectra were often contaminated by noise artifacts from the detector amplifiers, so we masked them out.

Figures 3 and 4 show the calibrated azimuthally-averaged power spectra for our target stars. To compare different power spectra from the same target, we normalized each azimuthally-averaged power spectrum so that the geometric mean of the first 10 data points in each spectrum equals 1. For Altair and 61 Cygni A and B we had more than three pairs of target and calibrator observations, i.e. calibrated power spectra, so we show the average of all the spectra and error bars representing the 68% confidence interval for each datum, estimated from the variation among the individual power spectra. The error is primarily due to differences in the atmosphere-telescope transfer function between object and calibrator. None of the calibrated power spectra deviate from a straight line by more than a typical error; all the targets are unresolved to the accuracy of a our measurements.

## 3 DISCUSSION

To interpret our observations we compared them to models of the IR emission from the solar zodiacal cloud. We constructed a model for exo-zodiacal emission based on the smooth component of the Kelsall et al. (1998) model of the solar system zodiacal cloud as seen by COBE/DIRBE, with emissivity $`ϵr^{0.34}`$ and a temperature $`T=286\mathrm{K}r^{0.467}L^{0.234}`$, where $`r`$ is the distance from the star in AU, and $`L`$ is the luminosity of the star in terms of $`L_{}`$. For a dust cloud consisting entirely of a single kind of dust particle of a given size and albedo, the $`L`$ exponent in the expression for the temperature is simply $`1/2`$ times the $`r`$ exponent (Backman & Paresce 1993).

The physics of the innermost part of the solar zodiacal dust is complicated (see Mann & MacQueen 1993), but our results are not sensitive to the details, because the hottest dust is too close to the star for us to resolve. We assume that the dust sublimates at a temperature of 1500 K, and allow this assumption to define the inner radius of the disk. We set the outer radius of the model to 3 AU, the heliocentric distance of the inner edge of our own main asteroid belt. Our conclusions are not sensitive to this assumption; decreasing the outer radius to 2 AU or increasing it to infinity makes a negligible difference in the visibility of the model, even for A stars.

The assumed surface density profile, however, does make a difference. A collisionless cloud of dust in approximately circular orbits spiraling into a star due to Poynting-Robertson drag that is steadily replenished at its outer edge attains an equilibrium surface density that is independent of radius (Wyatt and Whipple 1950, Briggs 1962). Models that fit data from the Helios space probes (Leinert et al 1981), the fit by Kelsall et al. (1998) to the COBE/DIRBE measurements and Good’s (1997) revised fit to the IRAS data all have surface densities that go roughly as $`r^{0.4}`$. This distribution appears to continue all the way in to the solar corona (MacQueen & Greely 1995). We find that in general, if we assume an $`r^\alpha `$ surface density profile, our upper limit for the 1 AU density of a given disk scales roughly as $`10^{\alpha /2}`$; disks with more dust towards the outer edge of the 11.6 micron emitting region are easier to resolve.

Likewise, the assumed temperature profile strongly affects our upper limits. Unfortunately, we know little about the temperature profile of the solar zodiacal cloud. COBE/DIRBE and IRAS only probed the dust thermal emission near 1 AU, and Helios measured the solar system cloud in scattered light, which does not indicate the dust temperature. We found that a dust cloud model with the IRAS temperature profile ($`T=266\mathrm{K}r^{0.359}L^{0.180}`$) was much easier to resolve than the model based on DIRBE measurements that we present here, especially for G and K stars.

To compare the models with the observations, we synthesized high resolution images of the model disks at an inclination of 30 degrees. We calculated the IR flux of the stars from the blackbody function, and obtained the parallaxes of the stars from the Hipparcos Catalog (ESA 1997). We inferred stellar radii and effective temperatures for each star from the literature and checked them by comparing the blackbody fluxes to spectral energy distributions based on photometry from the SIMBAD database (Egret et al. 1991). For Altair and Vega, we use the interferometrically measured angular diameters (1974) (they are 2.98 +/-0.14 mas and 3.24 mas). Stellar fluxes typically disagree with fitted blackbody curves by $`10\%`$ in the mid-infrared (Engelke 1990), but our method does not require precise photometry, and the blackbody numbers suffice for determining conservative upper limits. We computed the power spectra of the images, and normalized them just like the observed power spectra. In figures 3 and 4, the azimuthally-averaged power spectra for our target stars are compared to the extrapolated COBE/DIRBE model at a range of model surface densities. Disks with masses as high as $`10^3`$ times the mass of the solar disk will suffer collisional depletion in their inner regions, so they are unlikely to have the same structure as the solar disk. By neglecting this effect we are being conservative in our mass limits. The density of the densest model disk consistent with the data in each case is listed in table 1.

Altair

Our best upper limit is for Altair (spectral type A7, distance 5.1 pc); with 11 pairs of object and calibrator observations we were able to rule out a solar-type disk a few times $`10^3`$ as dense as our zodiacal cloud. Such a disk would have been marginally detectable by IRAS as a photometric excess.

Vega

IRAS detected no infrared excess in Vega’s spectral energy distribution at 12 microns, with an uncertainty of 0.8 Jy. This may be due to a central void in the disk interior to about 26 AU (Backman & Paresce 1993). Aumann et al (1984) suggested that Vega (A0, 7.8 pc) could have a hot grain component (500 K) with up to $`10^3`$ of the grain area of the observed component and not violate this limit. The apparent upward trend in the visibility data may be a symptom of resolved flux in the calibrator stars. We have only 3 object/calibrator pairs for Vega, not enough to test this hypothesis. Our upper limit is a solar-type disk with approximately $`3\times 10^3`$ times the density of the solar disk. This disk would have a $``$500 K emitting area of $`10^{24}\mathrm{cm}^2`$, about $`10^3`$ of the grain area of the observed component.

61 Cygni A and B

Though 61 Cygni is close to the galactic plane and surrounded by cool cirrus emission, Backman, Gillett and Low (1986) identified an IRAS point source with this binary system and deduced a far-infrared excess not unlike Vega’s. The color temperature of the excess suggests the presence of dust at distances $`>15`$ AU from either star. However these stars are dim (spectral types K5 and K7) and the region of the disk hot enough to emit strongly at 11.6 microns is close to the star and difficult to resolve; we could not detect a solar-type dust disk around either of these objects at any density, assuming the COBE/DIRBE model, or unless it had $`10^5`$ times the density of the solar disk, assuming the IRAS model.

70 Oph B

70 Oph is a binary (types K0 and K4) with a separation of 24 pixels (2.6 arcsec). We were able to assemble a power spectrum for B from 9 object/calibrator pairs, but the image of A fell on a part of the LWS chip that suffered from many bad pixels and was unusable. The image of A may also have been distorted by off-axis effects. 70 Oph B, like 61 Cygni A and B, is dim, making any dust around it cool and hard to detect at 11.6 microns.

$`\tau `$ Ceti

IRAS could have barely detected a disk with $`1000`$ times the emitting area of the solar disk around Tau Ceti (G8, 3.6 pc), the nearest G star. We have only three object/calibrator pairs for this object, not enough data to improve on this limit.

We are grateful to Dana Backman, Alycia Weinberger, Keith Matthews and Eric Gaidos for helpful discussions, and to Keith Matthews and Shri Kulkarni for assistance with the observations. This research has made use of the Simbad database, operated at CDS, Strasbourg, France. The observations reported here were obtained at the W. M. Keck Observatory, which is operated by the California Association for Research in Astronomy, a scientific partnership among California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. It was made possible by the generous financial support of the W. M. Keck Foundation.


text:
# Ionizing radiation in Smoothed Particle Hydrodynamics

## 1 Introduction

Smoothed Particle Hydrodynamics (sph) has become a numerical method widely used for addressing problems related to fluid flows in astrophysics. Due to its Lagrangian nature it is especially well suited for applications involving variations by many orders of magnitude in density. Examples for this type of applications are simulations of the collapse of molecular clouds and the formation of a stellar cluster, as performed by Klessen, Burkert & Bate \[Klessen et al. 1998\]. A comparison between grid based methods and sph was performed by Burkert, Bate & Bodenheimer \[Burkert et al. 1996\] and Bate & Burkert \[Bate & Burkert 1997\]. They applied both methods to the numerically demanding problem of gravitational collapse and fragmentation of a slightly perturbed rotating cloud core with an $`r^2`$ density profile. Both methods yielded the same qualitative results. Bate \[Bate 1998\] performed the first calculation which followed the collapse of a molecular cloud core in 3 dimensions down to a protostellar object in hydrodynamical equilibrium, thus spanning 17 (!) orders of magnitude in density. Other applications include accretion processes in massive circumbinary disks \[Bonnell & Bate 1994, Bate & Bonnell 1997\], the collapse of cloud cores induced by shock waves \[Vanhala & Cameron 1998\] or colliding clumps \[Bhattal et al. 1998\], the precession of accretion disks in binary systems \[Larwood & Papaloizou 1997\], the dynamical behaviour of massive protostellar disks \[Nelson et al. 1998\] or the formation of large scale structure and galaxies in the early universe \[Steinmetz 1996\].

A variety of physical processes are at work in the interstellar medium, like magnetic fields, radiation or thermal conductivity, necessitating their inclusion into numerical codes. This has already been achieved to a large extent in grid based methods like the magneto-hydrodynamics codes zeus \[Stone & Norman 1992\] or nirvana \[Ziegler, Yorke & Kaisig 1996\], or codes including effects of IR and UV radiation \[Yorke & Kaisig 1995, Sonnhalter, Preibisch & Yorke 1995, Richling & Yorke 1998\].

In contrast, the addition of physical processes to sph codes is just at its beginnings. Extensions achieved so far are sophisticated equations of state (e.g. Vanhala et al. 1998) and self-gravity. Some efforts were made to make sph faster and more accurate. The introduction of tree algorithms \[Barnes & Hut 1989, Press 1986, Benz et al. 1990\] and the use of GRAvity PipE (grape), a hardware device for fast computation of the gravitational N-body forces \[Umemura et al. 1993, Steinmetz 1996\], helped reducing the numerical effort for the gravitational force calculation and the determination of the nearest neighbours for each particle. Inutsuka \[Inutsuka 1995\] presented a Godunov-like solver for the Eulerian equations in sph thus enhancing the numerical treatment of shocks. The introduction of gravitational periodic boundaries \[Hernquist, Bouchet & Suto 1991, Klessen 1997\] allows the treatment of fragmentation and turbulence in molecular clouds without global collapse. The timestep problem which arises during isothermal collapse calculations at high densities is circumvented by the formation of sink particles, which substitute the innermost parts of the collapsing clump by one particle and accumulate the infalling mass and momenta \[Bate, Bonnell & Price 1995\].

The strength of sph lies in its Lagrangian nature, which makes it especially attractive for problems involving gravitational collapse and star formation. Applications like e.g. by Klessen et al. \[Klessen 1997\], which deal with the collapse and fragmentation of molecular clouds, neglect the feedback processes of newly born stars which act on their parental cloud through stellar winds, outflows and ionization. This simplification may be justified as long as the simulations deal with collapse on timescales smaller than $``$ 1 Myr, on which single and binary stars or T Tau-like clusters are formed \[Efremov & Elmegreen 1998\]. The case is different for larger timescales, on which OB subgroups and associations are formed. Neglecting feedback in these cases can lead to unphysical results, like a star formation efficiency of 100 per cent, since in the purely isothermal case all material will sooner or later be accreted onto the evolving protostellar cores. This is in strong contradiction to observations, which estimate a global star formation efficiency for ordinary molecular clouds of order 10 per cent \[Wilking & Lada 1985\]. Another possible effect of feedback is the induction of star formation due to the compression of cloud material by shock waves and ionization fronts.

In this paper we discuss the implementation of the effects of ionizing UV radiation by massive stars into sph calculations as a first step in order to perform collapse calculations on scales where OB-stars are formed in a more realistic way. This will in future applications allow us to assess questions like: How does the process of ionization by massive stars change the stellar initial mass function? What are the implications for the star formation efficiency? Can star formation be induced by ionization, and if yes, what are the time scales and the parameter space, for which induced star formation can be expected? These questions will be discussed in subsequent papers.

## 2 Physical problem

We incorporate the effects of ionizing radiation from hot stellar photospheres into sph by dividing the problem into three major substeps:

1. calculation of the UV radiation field by solving the time-independent, non-relativistic equation of radiative transfer,
2. determination of the ionization and recombination rates from the local radiation field, density and ionization fraction,
3. advancing the ionization state of the particles in time by solving the time-dependent ionization rate equation.

### 2.1 Calculation of the UV radiation field

Given a planar infall of ionizing photons from a distant source onto the border of the volume of interest with a flux $`J_0`$ Lyman continuum photons per time and square area, the resulting photon flux inside this volume is given by

$$J(s)=J_0\mathrm{exp}\left(\tau \left(s\right)\right),$$

where $`\tau (s)`$ is the optical depth for the ionizing photons along the line of sight parallel to the infall direction of the photons, and $`s`$ is the distance from the border of the integration volume along the line of sight:

$$\tau (s)=_0^s\left[\overline{\kappa }\left(s^{}\right)+\kappa _\mathrm{d}\left(s\right)\right]ds^{}.$$
(1)
We neglect the effect of ‘photon hardening’, i.e. the stronger absorption of weaker photons, and use an ‘effective’ absorption coefficient $`\overline{\kappa }`$, the mean of $`\kappa _\nu `$ over frequency, weighted by the spectrum of the source $`S_\nu `$:

$$\overline{\kappa }=n_\mathrm{H}\overline{\sigma }=n_\mathrm{H}\frac{S_\nu ^{(\mathrm{i})}\sigma _\nu d\nu }{S_{\mathrm{tot}}^{(\mathrm{i})}},$$
(2)
where $`\sigma _\nu `$ denotes the ionization cross section of hydrogen in the ground state and $`n_\mathrm{H}`$ the particle density of the H atoms.

The role of dust in Hii regions and its effect on ionizing radiation is still very uncertain \[Feldt et al. 1998\]. If dust is present, it will partially absorb UV photons, heat up and reemit the energy in the IR regime. Its first order effect can be included easily under the assumption of a homogeneous distribution of the dust in the Hii region. The corresponding contribution to the optical depth can be incorporated by adding the dust absorption coefficient at the Lyman border $`\kappa _\mathrm{d}`$ to the absorption coefficient in Eq. 1. $`\kappa _\mathrm{d}`$ depends on the dust model used and is regularily determined using Mie theory for grains with given distributions in size and shape. In this paper, we set $`\kappa _\mathrm{d}`$ to zero throughout.

We also neglect the diffuse field of Lyman continuum photons, which are being produced by recombinations of electrons into the ground level and which themselves possess sufficient energy for ionizing other H atoms. A thorough treatment of this radiation can only be achieved by detailed radiation transfer calculations as proposed e.g. by Yorke & Kaisig \[Yorke & Kaisig 1995\]. Instead we use the assumption of the validity of the ‘on the spot’ approximation as follows: due to the fact that the spectrum of the Lyman recombination photons as well as the the ionization cross section is strongly peaked at the Lyman border, a small amount of H atoms in the ionized region is sufficient to make the medium optically thick for the Lyman recombination photons. This leads to the absorption of these photons in the ultimate vicinity of their creation sites. As the creation of one photon is related to the creation of one H atom, its absorption leads to the destruction of one H atom. Thus the net effect of these photons on the local ionization structure is zero.

This assumption breaks down in regions next to OB stars, where due to the high UV flux the density of H atoms is not sufficient to make the medium optically thick to Lyman continuum photons. Next to ionization fronts, where the density of H atoms is much higher, the ‘on the spot’ assumption is nevertheless a good approximation. On further details refer to Yorke \[Yorke 1988\].

### 2.2 Ionization and recombination rates

The ionization rate in the medium is given by the sinks of the UV radiation field, since every ionization leads to the absorption of one UV photon:

$$=n_\mathrm{H}\overline{\sigma }J=𝐉,$$
(3)
where $`𝐉=J\widehat{𝐞}_s`$ is the flux vector in the direction $`\widehat{𝐞}_s`$ of the line of sight.

The recombination rate can be estimated as :

$$=n_\mathrm{e}^2\alpha _\mathrm{B}=n^2x^2\alpha _\mathrm{B},$$
(4)
with $`n`$ being the particle density of H atoms and protons together, $`n_\mathrm{e}`$ the particle density of free electrons, $`x=n_\mathrm{e}/n`$ the ionization fraction and $`\alpha _\mathrm{B}`$ the effective recombination coefficient under assumption of validity of the ‘on the spot’ approximation. The recombination coefficient $`\alpha `$ is given as the sum over the individual recombination coefficients $`\alpha _n`$, where the electron ends up in the atomic level $`n`$:

$$\alpha =\underset{n}{}\alpha _n.$$
(5)
Under the assumption of the ‘on the spot’ approximation recombinations into the ground level do not lead to any net effect and thus $`\alpha _1`$ can be neglected in Eq. 5. The resulting net recombination rate which is used in Eq. 4 is commonly called $`\alpha _\mathrm{B}`$ after the nomenclature introduced by Baker & Menzel \[Baker & Menzel 1962\]:

$$\alpha _\mathrm{B}=\underset{n=2}{\overset{\mathrm{}}{}}\alpha _n.$$

### 2.3 Ionization rate equation

Knowing the ionization and recombination rates, $``$ and $``$, the ionization fraction can be calculated from the ionization rate equation. The time dependency of the ionization fraction in the frame comoving with the corresponding particle, i.e. its Lagrangian formulation, is given by:

$$\frac{\mathrm{d}n_\mathrm{e}}{\mathrm{d}t}=.$$
(6)
### 2.4 Modeling the source

Since the spectral distribution of the UV radiation emitted by the photospheres of intermediate to high mass stars is very uncertain, we assume a black radiator with an effective temperature $`T_{}`$.

## 3 Numerical treatment

We developed two different methods for the numerical treatment of time dependent ionization in the sph calculations. Both have in common the method of finding paths from the ionizing source to the particles, along which the optical depth for the Lyman continuum photons can be calculated. They differ in the way the ionization rate is determined given the radiation field. Method A uses the sph formalism to calculate the divergence of the radiation field in Eq. 3. In method B we adopt a different approach also used in grid methods, where we derive the ionization rate from the difference in the numbers of photons entering and leaving a particle.

### 3.1 Finding the evaluation points on the path towards the source

First, we specify the position, the rate of ionizing photons $`S_{\mathrm{tot}}`$ and $`\overline{\sigma }`$ (from Eq. 2) of the source.

For each particle i we now proceed in the following way (see Fig. 1): Given the list of nearest neighbours of particle i, which has to be determined anyway for the sph formalism, we look for the particle j in the list, closest to the line of sight defined by the smallest angle $`\mathrm{\Theta }`$ between the line connecting the particles i and j and the line of sight. We choose the angle between, not the distance from, the line of sight, since we are interested in controlling the error in the direction towards the source. This is not garanteed by the latter criterion.

We store this particle in a list and determine the evaluation point Sj as the projected particle position on the line of sight. To determine the next evaluation point Sk even closer to the source we now repeat this method using the neighbor list of particle j and so forth until we reach the source.

### 3.2 Calculating the optical depth and ionization rate for the particles

#### 3.2.1 Method A: sph formalism method

Now the path from the source to particle i is known, and the integration of Eq. (1) can be discretized by using the evaluation points $`S_\mathrm{i}`$. The value for $`n_\mathrm{H}`$ can be estimated by using the sph smoothing formalism:

$$n_\mathrm{H}\left(𝐫\right)=n_{\mathrm{H},\mathrm{i}}W\left(𝐫𝐫_\mathrm{i}\right),$$
(7)
where the sum runs over the particle corresponding to the evaluation point and its nearest neighbours. $`W`$ is the weight factor for each neighboring particle provided by the smoothing kernel. We calculate the optical depth along the line of sight by applying the Trapezian Formula, until we reach particle i:

$$\tau _{\mathrm{k}+1}=\tau _\mathrm{k}+\frac{1}{2}\overline{\sigma }\left(s_{\mathrm{k}+1}s_\mathrm{k}\right)\left(n_{\mathrm{H},\mathrm{k}+1}+n_{\mathrm{H},\mathrm{k}}\right),$$

with $`s_\mathrm{k}`$ being the position of the evaluation point on the line of sight. Note that this treatment neglects the effects of scattering of the ionizing photons by recombination or dust.

The distance between two successive evaluation points is smaller or equal to the local smoothing length, which determines the largest distance of the particles included in the nearest neighbour list as well as the spatial resolution. This guarantees that the line of sight integration of Eq. (1) is discretized into a reasonable amount of substeps, consistent with the resolution given by the underlying particle distribution.

The flux of ionizing photons at the position of particle i into the direction of photon propagation $`\widehat{𝐞}_s`$ is then given by:

$$𝐉_\mathrm{i}=J_0\widehat{𝐞}_s\mathrm{exp}\left(\tau \left(s\right)\right).$$

With the ionizing flux known at the particle positions, the nabla operator in Eq. (3) can be calculated by the sph formalism. It is given for each particle i as the sum over its neighbours:

$$_\mathrm{i}=\frac{m_\mathrm{j}}{\rho _\mathrm{j}}𝐉_\mathrm{j}_\mathrm{i}W_{\mathrm{j},\mathrm{i}}.$$
(8)
Now we are able to solve Eq. 6, which we write as:

$$\frac{\mathrm{d}x_\mathrm{i}}{\mathrm{d}t}=_\mathrm{i}n_\mathrm{i}x_\mathrm{i}^2\alpha _\mathrm{B}.$$
(9)
The time scale for the establishment of ionization equilibrium is given by $`1/(n\alpha _\mathrm{B})`$, which is regularly much shorter than the dynamical and gravitational timescales we are interested in. In order to avoid small timesteps arising from the usage of explicit methods, we use an implicit scheme. The first order discretization of Eq. 9 over a time interval $`\mathrm{\Delta }t`$ is given by:

$$x_\mathrm{i}^{\mathrm{n}+1}=x_\mathrm{i}^\mathrm{n}+\mathrm{\Delta }t\left(_\mathrm{i}^{\mathrm{n}+1}n_\mathrm{i}^{\mathrm{n}+1}x_\mathrm{i}^{\mathrm{n}+1}\alpha _\mathrm{B}\right),$$
(10)
where the indices n and n$`+1`$ denote the values at the beginning and the end of the actual timestep $`\mathrm{\Delta }t`$, respectively. We already know all the values on the right hand side from advancing the particles by the sph formalism, except the value for $`_\mathrm{i}^{\mathrm{n}+1}`$. Therefore a fully consistent implicit treatment is not feasible. We use the following guess for this value:

$$_\mathrm{i}^{\mathrm{n}+1}=_\mathrm{i}^\mathrm{n}\frac{1\mathrm{exp}\left(n_\mathrm{i}^{\mathrm{n}+1}\overline{\sigma }a_\mathrm{i}^{\mathrm{n}+1}\left(1x_\mathrm{i}^{\mathrm{n}+1}\right)\right)}{1\mathrm{exp}\left(n_\mathrm{i}^{\mathrm{n}+1}\overline{\sigma }a_\mathrm{i}^{\mathrm{n}+1}\left(1x_\mathrm{i}^\mathrm{n}\right)\right)}.$$
(11)
In this equation, we assign an effective radius $`a_\mathrm{i}`$ to each particle i proportional to the mean particle separation, given by $`a_\mathrm{i}=(M_\mathrm{i}/\rho _\mathrm{i})^{1/3}`$. This is the estimate of the size of a region with the particle mass $`M_\mathrm{i}`$ and density $`\rho _\mathrm{i}`$. The factor with the exponentials on the right hand side accounts for the effect of higher absorption and hence ionization rate with decreasing ionization fraction.

We must use the effective radius $`a_\mathrm{i}`$ in Eq. 11 instead of the smoothing length $`h`$, since the method works analogous to implementations in grid codes. In contrast to the sph formalism, each particle now represents a volume of total mass $`M_\mathrm{i}`$ and density $`\rho _\mathrm{i}`$, in which ionizing radiation enters on one side and leaves on the opposite side. The size of this volume is given by $`a_\mathrm{i}`$ as defined above. It is proportional to the particle spacing.

In contrast, $`h_i`$ differs from the mean particle separation as it is defined by the condition that there is a fixed number of neighbors $`N_{\mathrm{neigh}}`$ of mass $`M`$ in the sphere with radius $`2h_i`$ and is thus given as

$$h_i=\left(\frac{3N_{\mathrm{neigh}}M}{32\pi \rho }\right)^{1/3}.$$

It depends on $`N_{\mathrm{neigh}}`$ and can therefore not be used instead of $`a_\mathrm{i}`$ in Eq. 11.

One consequence of the discretization of the ionization rate equation is that the solution in ionized regions tends to oscillate around the equilibrium value. In order to avoid small timesteps arising from this, we set the ionization fraction $`x`$ of particles with an $`x>0.95`$ to the equilibrium value $`x_\mathrm{E}`$, which is defined by $`\mathrm{d}n_\mathrm{e}/\mathrm{d}t=0`$ in Eq. 6:

$$\frac{\mathrm{d}x}{\mathrm{d}t}=\frac{1}{n}\frac{\mathrm{d}n_\mathrm{e}}{\mathrm{d}t}=\overline{\sigma }(1x_\mathrm{E})Jnx_\mathrm{E}^2\alpha _\mathrm{B}=0.$$

With $`k=\overline{\sigma }J/(n\alpha _\mathrm{B})`$ follows that

$$x_\mathrm{E}=\frac{1}{2}\left[\left(k^2+4k\right)^{1/2}k\right].$$

This method works well in absolutely smooth, noise free particle distributions. However, if one wishes to initially distribute the particles randomly in space, one runs into problems. The sum in Eq. 8 is very sensitive to noisy particle distributions. Eventually the noise can be so high, that the error of the sum introduced by noise reaches the order of the sum itself. The ionization rate then locally drops below zero for some particles, which can only be avoided by smoothing the ionization rate spatially over several smoothing lengths. The result is poor resolution. We circumvent this problem in method B.

#### 3.2.2 Method B: grid based method

In this case, a different method is used to discretize the calculation of the optical depth. We determine the positions of the evaluation points i along the line of sight as described in Sect. 3.1 and calculate the hydrogen density $`n_{\mathrm{H},\mathrm{i}}`$ at these positions using Eq. 7. The path is then divided into pieces with length $`\mathrm{\Delta }s_\mathrm{i}=(s_{\mathrm{i}+1}s_{\mathrm{i}1})/2`$, assuming a constant hydrogen density $`n_{\mathrm{H},\mathrm{i}}`$ along each interval. The optical depth for one piece can then be approximated by

$$\mathrm{\Delta }\tau _\mathrm{i}=\overline{\sigma }n_{\mathrm{H},\mathrm{i}}\mathrm{\Delta }s_\mathrm{i}.$$

These contributions to the optical depth are summed up until we reach the position located one effective radius $`a_\mathrm{i}`$ before the position of particle k. A first order approximation for the ionization rate is now given by

$$_\mathrm{k}=\frac{J_0}{2a_\mathrm{k}n_\mathrm{k}}\mathrm{exp}\left(\tau _{\mathrm{k}\mathrm{a}}\right)\left(1\mathrm{exp}\left(\mathrm{\Delta }\tau _\mathrm{k}\right)\right),$$

where $`\tau _{\mathrm{k}\mathrm{a}}=\mathrm{\Delta }\tau _\mathrm{i}`$ denotes the optical depth one effective radius before the particle’s position and $`\mathrm{\Delta }\tau _\mathrm{k}=2a_\mathrm{k}n_{\mathrm{H},\mathrm{k}}\overline{\sigma }`$ the optical depth across the particle.

With the ionization rate derived above we solve the ionization rate equation as described for case A. One can easily show that Eqs. (10) and (11) now give the exact implicit first order discretization for Eq. (9). The solution now approaches the equilibrium value $`x_\mathrm{E}`$ in the ionized regions without the instabilities mentioned in method A. It is not necessary to set $`x`$ artificially to $`x_\mathrm{E}`$.

Method A seems to be the more consistent method since it uses the sph formalism for the calculation of $``$. This is the reason why it is also discussed in this paper. Nevertheless we prefer method B due to its robustness against noisy particle distributions and higher consistency concerning the integration scheme and have applied it to a couple of test cases.

### 3.3 Computational effort

If the procedure explained above is used, the computational effort for the line of sight integration scales approximately as $`N^{4/3}`$, since the integration has to be done for each of the $`N`$ particles, and the average number of evaluation points on each line of sight scales as $`N^{1/3}`$.

We can reduce the exponent from $`4/3`$ to 1 by introducing a ‘tolerance angle’ $`\mathrm{\Theta }_{\mathrm{tol}}`$. Suppose we determine the particles along the line of sight as expalined . As soon as $`\mathrm{\Theta }`$ for a particle j along the line of sight towards the source is smaller than $`\mathrm{\Theta }_{\mathrm{tol}}`$ we stop our search here. The optical depth of this particle $`\tau _j`$ is then used as an estimate of the optical depth along the remaining part of the line of sight from the source to Sj. Thus no integration is needed for this part of the path. One only has to make sure that $`\tau _j`$ is already known, i.e. that the line of sight integration for particle j has been performed earlier. In this case, the average number of evaluation points $`I`$ per line of sight only depends on $`\mathrm{\Theta }_{\mathrm{tol}}`$ for large $`N`$. As shown in Fig. 2, $`I`$ becomes constant for large $`N`$ and decreases with increasing $`\mathrm{\Theta }_{\mathrm{tol}}`$. As soon as $`I`$ becomes independent of $`N`$ the total computational effort for all lines of sight together scales $`N`$.

We demonstrate the effects of using the tolerance angle on the accuracy of the ionization rate calculation in Fig. 3. Histograms are plotted for the errors in $``$ and $`\tau `$ for calculations with $`\mathrm{\Theta }_{\mathrm{tol}}=0.5^{}`$, $`1^{}`$, $`2^{}`$ and $`90^{}`$ compared to $`\mathrm{\Theta }_{\mathrm{tol}}=0^{}`$. As the particle distribution we chose the evolved state of a numerical simulation which studies the compression and collapse of a dense clump within the UV field of an OB association using 200 000 particles. The results of this calculation will be presented elsewhere \[Kessel & Burkert 1999\]. Note that $`\mathrm{\Theta }_{\mathrm{tol}}=90^{}`$ represents the worst case, since the tolerance angle criterion now is fulfilled for every particle with minimal $`\mathrm{\Theta }`$ per search through the nearest neighbour list.

The particles which are most affected by the tolerance angle criterion lie next to the borders of shadows cast by optically thick regions, since here the path for the integration along the line of sight may be bent through the optically thick region, thus decreasing the ionizing flux artificially. In the opposite case, the path may lead around the opaque region, increasing the ionizing flux at the position of a particle in the shadow. These extreme cases lead to the tail in the error histograms in Fig. 3. Applying the tolerance angle criterion thus numerically blurs shadows.

The mean errors in $`\tau `$ are $`1.3`$ per cent for $`\mathrm{\Theta }_{\mathrm{tol}}=0.5^{}`$, $`2.2`$ per cent for $`\mathrm{\Theta }_{\mathrm{tol}}=1.0^{}`$, $`3.4`$ per cent for $`\mathrm{\Theta }_{\mathrm{tol}}=2.0^{}`$ and $`11.2`$ per cent for $`\mathrm{\Theta }_{\mathrm{tol}}=90^{}`$. The correspnding mean errors in $``$ are $`2.8`$, $`4.1`$, $`5.7`$ and $`13.3`$ per cent, respectively. For the remaining test cases presented in this paper the choice of $`\mathrm{\Theta }_{\mathrm{tol}}`$ has no effect, since they deal with one-dimensional problems, in which the optical depth is only a function of distance from the source. Applying the tolerance angle criterion only shifts the evaluation points from the lines of sight in directions perpendicular to these, along which there is no change in the optical depth. Indeed even the choice $`\mathrm{\Theta }_{\mathrm{tol}}=90^{}`$ gives the same results in the one-dimensional test cases as for $`\mathrm{\Theta }_{\mathrm{tol}}=0^{}`$. Thus the errors introduced by the angle criterion must be checked with problems in which this symmetry is broken and shadows are present, as the one mentioned above.

### 3.4 Smoothing the ionization front

For reasons of noise reduction we smooth the ionization front, which is not resolvable by the sph representation, over a distance of the order of one local smoothing length. Nature provides a simple way for doing this. The width of the ionization zone is of the order of one photon mean free path length,

$$d=(\overline{\sigma }n_\mathrm{H})^1,$$
(12)
where $`\overline{\sigma }`$ is the net absorption cross section for ionizing photons as defined in Eq. 2.

Since we cannot resolve the ionization region anyway, we are free to adjust $`\sigma `$ in a way that the width of the ionization region given by Eq. 12 is equal to a constant factor $`C1`$ times the local smoothing length $`h`$, but never larger than the value $`\overline{\sigma }`$ given by Eq. 2:

$$\sigma =\mathrm{min}[\overline{\sigma },\left(n_\mathrm{H}Ch\right)^1].$$

Test calculations have shown that a good value is $`C=0.1`$. It has proven to sufficiently reduce numerical noise introduced into the ionization structure by noise in the particle distribution and at the same time to keep the resolution of ionization fronts better than the resolution of the sph formalism in order not to worsen the overall resolution. Note that, when “smoothing” the ionization front over $`0.1`$ times the smoothing length, the noise reducing effect is not caused by the spatial smoothing, since it is ten times smaller than the sph smoothing. It rather results from a larger number of time steps needed to ionize a particle in the front from an ionization fraction of $`x=0`$ to $`x1`$. This gives the neighbouring particles the opportunity to react to the changed state in a smoother way.

### 3.5 Heating effect

We assume that heating and cooling effects lead to an equilibrium temperature of 10 000 K in the ionized gas penetrated by ionizing radiation. The cross sections for elastic electron–electron and electron–proton scattering are of the order $`10^{13}\mathrm{cm}^2`$. Together with a mean velocity of the electrons of the order of $`600\mathrm{km}\mathrm{s}^1`$ the thermalization timescale for the energies of the ejected electrons is far less than a year for densities of 1 particle cm-3, which is many orders of magnitude smaller than the dynamical timescale. Thermalization thus occurs quasi instantaneously. This process runs even more rapidly for higher densities. Thus we are allowed to treat the gas behind the ionization front as thermalized. We set the internal energy to:

$$e=xe_{10000}+(1x)e_{\mathrm{cold}},$$

with $`e_{10000}`$ being the internal energy corresponding to a temperature of 10 000 K for ionized hydrogen, and $`e_{\mathrm{cold}}`$ to the internal energy for the 10 K cold, neutral gas. Note that this method does not properly treat recombination zones, since in this case one needs the correct inclusion of the heating and cooling processes in order to achieve the correct gas temperatures, sound velocities and pressures. Also, the equilibrium temperature in Hii-regions can vary by 20 per cent from this value. These deviations can also only be taken into account by proper treatment of heating and cooling.

## 4 Tests of the numerical treatment

Although being of one-dimensional nature, the following test problems were performed fully in three dimensions.

### 4.1 Test 1: Ionization of a slab with constant density

With this problem we test the implementation of the time-dependent ionization rate equation by ionizing a slab of Hi gas of constant density $`n`$ with ionizing radiation falling perpendicular onto one of the boundary surfaces. With hydrodynamics switched off, we let the ionization front traverse the slab with a constant velocity $`v_\mathrm{f}`$. To achieve this, we have to vary the infalling photon flux with time. It is given by

$$J(t)=J_\mathrm{f}+J_\mathrm{t}=nv_\mathrm{f}+n^2\alpha _\mathrm{B}v_\mathrm{f}t,$$

where the first term on the right hand side is the flux which provides the photons being absorbed in the ionization front. The second, time-dependent term equals the loss of photons on their way through the slab until they reach the front.

For the initial setup we place a number $`N`$ of particles randomly into a slab with length-to-height and length-to-width ratios of 10. Subsequently we let the particle distribution relax by evolving it isothermally within the slab, adding a damping term to the force law. This is necessary to diminish the numerical noise which was introduced by the random distribution. We now have an ensemble of the particles which does not possess any privileged directions and which represents a gas of constant density and temperature. We use this distribution as our starting configuration. From now on we keep the particles fixed in space and switch off hydrodynamics.

The test was performed for a total number of $`N=\mathrm{2\hspace{0.17em}000},\mathrm{16\hspace{0.17em}000}`$ and $`\mathrm{128\hspace{0.17em}000}`$ particles. Since the spatial resolution for sph calculations scales as $`N^{1/3}`$ (with number of neighbours $`N_{\mathrm{neigh}}`$ per particle fixed), this yields an increase of linear resolution of a factor of two from one simulation to the simulation with next higher resolution. The results of these tests are shown in Fig. 4.

The mean relative errors between the theoretical result and the calculations decrease linearly with increasing resolution, consistent with our first order discretization of both the line of sight integration and the time dependent ionization equation. The error also decreases with time as the representation of the ionization front gets thinner and thinner compared to the already ionized region. The spread in the errors for $`N=16000`$ and $`N=128000`$ results from the fact that in these cases the numerical solution oscillates around the theoretical solution, sometimes being larger than the latter, sometimes smaller.

### 4.2 Test 2: Ionization of a slab with density gradient

We proceed as in test 1, with the difference that we choose a slab with a constant density gradient in the direction of photon propagation. We choose the time dependence of $`J`$ such that the ionization front should travel through the gas with constant $`v_\mathrm{f}`$. $`J`$ is given by:

$`J(t)=n_0v_\mathrm{f}+\left(\alpha _\mathrm{B}n_0^2+{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}x}}v_\mathrm{f}\right)v_\mathrm{f}t+\alpha _\mathrm{B}n_0{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}x}}v_\mathrm{f}^2t^2+`$
$`\alpha _\mathrm{B}\left({\displaystyle \frac{\mathrm{d}n}{\mathrm{d}x}}\right)^2v_\mathrm{f}^3t^3,`$
where $`n_0`$ denotes the density at the surface where the radiation penetrates the slab and $`\frac{\mathrm{d}n}{\mathrm{d}x}`$ is the density gradient.

In Fig. 5 we plot the ionized mass for the theoretical solution and the numerical simulations against time. The numerical results converge against the theoretical solution with increasing resolution. The deviations at $`t>0.9`$ are caused by the ionization front reaching the rear boundary of the slab.

Note that the version of sphi used in this paper is not able to follow ionization fronts exactly which travel faster than one local smoothing length per time step. This must be taken into account during the timestep determination. In applications with fast ionization fronts (typically R-type fronts in the early phases of the evolution of Hii-regions) this criterion can lead to very small timesteps and thus to a high amount of CPU time needed. A version which circumvents this problem is being developed.

### 4.3 Test 3: Coupling of ionization and hydrodynamics

For this test we adopt the problem mentioned by Lefloch & Lazareff \[Lefloch & Lazareff 1994\]. A box filled with atomic hydrogen of particle density $`n_0=10`$ cm-3 and temperature $`T_{\mathrm{cold}}=100`$ K is exposed to ionizing radiation, with the photon flux increasing from zero linearly with time with a rate $`\mathrm{d}\mathrm{\Phi }/\mathrm{d}t=5.0710^2`$ cm-2 s-2. There exists an analytical solution to this problem, which is self similar in the sense that physical values at position $`x`$ measured in the direction of the photon flow at time $`t`$ are only functions of $`x/t`$. This means: the structure is stretched with time. The convergence of the code towards the correct solution with increasing resolution can be tested in one calculation, since for all appearing structures the ratio between structure sizes and smoothing lengths increases linearly with time.

The resulting structure is the following: an isothermal shock is driven into the neutral medium, sweeping up a dense layer of material. This is followed by an ionization front which leaves the ionized material in quasi-static equilibrium (see Figs. 6,7). Using the parameter $`\mathrm{\Lambda }=\alpha ^1(\mathrm{d}\mathrm{\Phi }/\mathrm{d}t)`$, Lefloch & Lazareff \[Lefloch & Lazareff 1994\] find the following analytical solution:

$`\mathrm{\Lambda }=n_\mathrm{i}^2V_\mathrm{i}`$
$`n_\mathrm{i}=\left({\displaystyle \frac{n_0\mathrm{\Lambda }^2}{c_\mathrm{i}^2}}\right)^{\frac{1}{5}}`$ $`V_\mathrm{i}=\left({\displaystyle \frac{\mathrm{\Lambda }c_\mathrm{i}^4}{n_0^2}}\right)^{\frac{1}{5}}`$
$`n_\mathrm{c}=n_0\left({\displaystyle \frac{\mathrm{\Lambda }}{n_\mathrm{i}^2c_\mathrm{n}}}\right)^2`$ $`V_\mathrm{s}=c_\mathrm{n}\left({\displaystyle \frac{n_\mathrm{c}}{n_0}}\right)^{\frac{1}{2}},`$
where $`n_\mathrm{i}`$, $`n_0`$ and $`n_\mathrm{c}`$ denote the particle densities of the ionized gas, the undisturbed neutral gas and the gas in the compressed layer, respectively, and $`V_\mathrm{i}`$ and $`V_\mathrm{s}`$ the velocities of the ionization front and the shock front, respectively.

We adopt $`\alpha _\mathrm{B}=2.710^{13}`$ cm3 s-1 from Lefloch & Lazareff \[Lefloch & Lazareff 1994\] in order to directly compare the results of the sphi code to those of their grid-based method using a piecewise linear scheme for the advection terms proposed by Van Leer \[Van Leer 1979\]. The resolution of 192 grid cells along the slab of their calculations, from which they derived their results, is comparable with the one used in our high resolution case. We use the same method as described in Sect. 4.3 to produce the initial conditions. No gas is allowed to enter or leave the surface.

Table 1 lists the result of this comparison. The sphi calculation slightly underestimates $`V_\mathrm{s}`$ and $`V_\mathrm{i}`$, as is also observed for the grid code. The errors of order 5 per cent are comparable to those achieved by Lefloch & Lazareff \[Lefloch & Lazareff 1994\].

In the early phases, i.e. low resolution, the poor treatment of the ionization front leads to irregularities in the ionized region and thus produces sound waves travelling back and forth between the boundary to the left and the ionization front (Figs. 6, 7), which decrease in power as time increases, i.e. at higher resolution. With increasing resolution, i.e. increasing ratio of layer thickness to smoothing length, the representation of the dense layer and the shock front improves (Figs. 8, 9).

## 5 Summary

The method presented in this paper allows the treatment of the dynamical effects of ionizing radiation in sph calculations. Thus the study of astrophysical problems arising from ionization, like the impact of ionizing radiation from newly born stars onto the evolution of their parental molecular clouds or the more consistent treatment of heating by OB associations in galaxy dynamics calculations are now feasible for the first time with sphi in 3 dimensions. We demonstrate that the code is able to treat time-dependent ionization, the related heating effects and hydrodynamics correctly. Our first applications, detailed calculations of photoionization induced collapse in molecular clouds and results obtained from them, will be presented in a subsequent paper.

To allow the correct treatment of recombination zones, one has to include the effects of time dependent heating and cooling processes by ionization and recombination, emission of forbidden lines and thermal radiation from dust. Another important aspect which was neglected here is the effect of the diffuse Lyman continuum recombination field. It can lead to the penetration of regions shielded from the direct ionizing radiation by the ionization front, which is e.g. seen in calculations of photoevaporating protostellar disks \[Yorke & Welz 1996, Richling & Yorke 1998\]. An implementation of these processes into our sphi code is planned in the future.

## Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (DFG), grant Bu 842/4. We’d also like to kindly thank Matthew Bate and Ralf Klessen for the useful discussions concerning the capabilities and implementation of the sph method.


text:
# Localized structures in coupled Ginzburg–Landau equations

## 1 Introduction

When an extended system is close to a Hopf bifurcation leading to uniform oscillations, the amplitude of the oscillations can be generically described in terms of the complex Ginzburg-Landau (CGL) equation . When there are two fields becoming unstable at the same bifurcation, coupled complex Ginzburg-Landau equations (CCGL) should be used instead. This model set of equations appears in a number of contexts including convection in binary mixtures and transverse instabilities in unpolarized lasers .

Coherent structures such as fronts, shocks, pulses, and other localized objects play an important role in the dynamics of extended systems . In particular, for the complex Ginzburg-Landau equation, they provide the building blocks from which some kinds of spatiotemporally chaotic behavior are built-up . A systematic study of localized structures in CCGL equations in one spatial dimension was initiated in .

Here we present results on one dimensional CCGL equations in parameter ranges such that they can be written as

$$_tA_\pm =\mu A_\pm +(1+i\alpha )_x^2A_\pm (1+i\beta )\left(|A_\pm |^2+\gamma |A_{}|^2\right)A_\pm .$$
(1)
Group velocity terms of the form $`\pm v_g_xA_\pm `$ are explicitly excluded, and $`\gamma `$ is restricted to take real values (without additional loss of generality, $`\alpha `$ and $`\beta `$ are also real parameters). In addition we just consider $`1+\alpha \beta >0`$ (Benjamin-Feir stable range). These restrictions are the appropriate ones for the description of transverse laser instabilities . In that case $`A_\pm `$ are related to the two orthogonal circularly polarized light components. We further restrict our study to the case $`0<\gamma <1`$ which is the range obtained when atomic properties in the laser medium favor linearly polarized emission. In terms of the wave amplitudes $`A_\pm `$, wave coexistence is preferred.

## 2 Numerical studies

Many experiments on traveling wave systems or numerical simulations of Ginzburg–Landau–type equations  exhibit local structures that have a shape essentially time–independent and propagate with a constant velocity, at least during an interval of time where they appear to be coherent structures . In order to analyze these structures it is common to reduce the initial partial differential equation to a set of ordinary differential equations by restricting the class of solutions to uniformly traveling ones. Localized structures are homoclinic or heteroclinic orbits in this reduced dynamical system, that is they approach simple solutions (typically plane waves) in opposite parts of the system, whereas they exhibit a distinct shape in between.

Instead of looking for solutions of the reduced dynamical system, we prefer here to resort to direct numerical solution of (1) under different initial conditions. A pseudo–spectral code  with periodic boundary conditions and a second–order accuracy in time is used. Spatial resolution was typically 512 modes. Time step was typically $`0.05`$. The system size was always taken to be $`L=512`$. Several kinds of localized objects which maintain coherence for a time appear and travel around the system. Different initial conditions give birth to different kinds of structures. Some of them decay shortly, and the qualitative dynamics at long times becomes determined by the remaining ones, and essentially independent of the initial conditions.

The upper part of Fig. 1 shows the spatiotemporal evolution of $`|A_+(x,t)|`$ and $`|A_{}(x,t)|`$ at parameter values $`\alpha =0.35,\beta =2.0\text{ and }\gamma =0.2`$. Time runs upwards and $`x`$ is represented in the horizontal direction. Lighter grey corresponds to the maximum values of $`|A_\pm (x,t)|`$ and darker to the minima.

This particular evolution was obtained starting from $`A_+(x,0)`$ equal to the Nozaki-Bekki hole, a known analytical solution of the single Ginzburg-Landau equation , and for $`A_{}(x,0)`$ a Nozaki-Bekki pulse . These are not exact solutions of the set of equations (1) so that this initial condition decays and gives rise to complex spatiotemporal structures. After a transient that will be described below, the configuration of the system consists in portions with a modulus nearly constant (corresponding to plane wave states) interrupted by localized objects with particle-like behavior. Dark features in $`|A_+|`$ appear where $`|A_{}|`$ has bright features, thus indicating that the localized object carries a kind of anticorrelation between the fields. The lower panels of Fig. 1 show the modulus of the two fields at $`t=399`$ and $`x300`$, where one of such objects is present. One of the components shows a maximum in the modulus, whereas the other displays a deep minimum. We can call this object a “hole–maximum pair”. It seems to be a dissipative analog of the ‘out-gap’ solitons appearing in Kerr media with a grating , and here it is the characteristic object building-up the disordered intermittent dynamics seen at long times. It is clear that these objects connect the plane wave states (that is the constant modulus regions) filling most of the system. Before reaching the asymptotic state just described, the system evolves through configurations where additional kinds of localized objects are seen. The presence of the Nozaki-Bekki hole-pulse pair as initial condition in the central part of Fig. 1 gives birth to a pair of fronts which replace the initial lateral plane-waves by new ones. Interestingly, a different kind of localized object is seen to form just where the initial hole-pulse pair was placed. A close-up of it at $`t=90`$ is displayed in Fig. 2. It is a kind of coupled maximum-maximum pair. The moduli of the two fields are superposed in the central panel showing the full object. The lateral small bumps are propagating waves that travel towards the central maxima. Thus the center of the coherent structure acts as a wave sink .

In Figure 3 the spatiotemporal evolution of $`|A_+(x,t)|`$ and $`|A_{}(x,t)|`$ was obtained using as initial conditions a sharp phase jump at the center of the system, with small random white noise added. The parameter values are $`\alpha =0.6,\beta =1.4\text{ and }\gamma =0.7`$. After a short time, the system reaches a state dominated by branching hole–hole pair structures. Lighter grey correspond to the maximum values of $`|A_\pm (x,t)|`$ and darker to the minima. The two big triangles correspond to regions of constant modulus, that is, plane waves. The bottom panels show $`|A_+|`$ and $`|A_{}|`$ in a portion of the system at these early times. Both are superposed in the central panel to show the complete matching of the two fields.

At longer times, all the hole-hole pairs disappear from the system, thus indicating that they are not stable objects at this value of the parameters. The system decays to the same state as at the end of Fig. 1: the dominant coherent structures are the maximum-hole pairs.

## 3 Exact solutions

The different spatiotemporal evolutions shown in the previous figures (1)–(3) are themselves interesting enough for a detailed study. The localized objects appearing in the simulations are clearly responsible for most of the complex dynamics in the system. We can interpret some of the observed structures from a simple ansatz:

$$A_+(x,t)=e^{i\phi }A_{}(x,t)$$
(2)
where $`\phi `$ is constant, and $`A_{}(x,t)`$ is any solution of the single CGL equation:

$$_tA_{}=A_{}+b_x^2A_{}cA_{}^2A_{},$$
(3)
where $`b=1+i\alpha `$ and $`c=(1+\gamma )+i(1+\gamma )\beta `$ .

This simple ansatz gives us a rather rich set of exact solutions: for each known analytical solution of the single CGL equation (3), there is a corresponding solution of the CCGL equation set, in which $`A_{}`$ and $`A_+`$ have essentially the same shape except for a constant global phase. In particular, hole, pulse, shock, and front solutions are localized solutions analytically known for the single equation , so that hole-hole, pulse-pulse, shock-shock and front-front pairs are immediately found as analytical solutions of the CCGL set. In particular pulse-pulse and hole-hole structures are present in Figs. 1 to 3, and turn out to be well described by the ansatz (2).

It is worthwhile to note that the studies of instability for these objects in the complex Ginzburg-Landau equation are immediately translated into instability results for the paired structures in CCGL equations.

## 4 Conclusion

In summary, we have shown numerically the existence of different kinds of localized objects, responsible for the complex behavior or solutions of the CCGL equations. Some of these objects can be understood in terms of exact solutions arising from a simple ansatz. A more detailed analysis is still needed, however. In particular, the hole-maximum structure, which appears as the dominant coherent structure at long times, can not be described by our ansatz. In addition, much more work is needed in order to establish the stability properties of the different objects, and the nature of their interactions. In a recent work new exact solutions of equation 1 were obtained by using the Painlevé expansion method. The authors describe these solutions as analogues of the Nozaki-Bekki solutions . Comparison of these solutions, different from the ansatz (2), with our numerical results is under progress.

Financial support from DGICYT Projects PB94-1167 and PB97-141-C02-01 is acknowledged. R.M. Acknowledge financial support from CONICYT-Fondo Clemente Estable (Uruguay)


text:
# References

Time invariance violation in photon-atom
and photon-molecule interactions.

V.A.Kuz’menko

Troitsk Institute of Innovation and Fusion Research,

Troitsk, Moscow region, 142092, Russian Federation.

## Abstract

A direct experimental proof of strong T-invariance violation in interactions of the photons with atoms and molecules exists in the molecular physics.

PACS number: 42.50.-p

”Even for the physicist \[a\] description in plain language will be a criterion of the degree of understanding that has been reached”. Werner Heisenberg, ”Physics and Philosophy” (Harper Bros., New York, p.168, 1958).

The searches for the P- and T- invariance violation are being conducted in the atomic physics for many years . As concerns the T-invariance violation, these searchers are based on the indirect indications, namely, they assume the neutron, electron and atoms to possess an electric dipole momentum . Time reversibility invariance demands, that the probability amplitudes of the direct and reverse processes should be equal . The problem of the direct test of T-invariance is not discussed in the elementary particle physics. Apparently, this is due to the obvious facts, that reverse in time, for example, of the $`K^0`$ meson decay is practically unreal.

On the contrary, the reverse in time of the process of interaction of the photon with atoms and molecules, in principle, can be easily implemented experimentally. The time invariance in this case, apparently, implies the absorption cross-section of the photons to be equal to that of their stimulating emission. We can effect on the atoms or molecules by laser radiation to excite them. Then we can influence by the laser radiation once more to deexcite them. After that, it is possible to measure the cross-section of the direct process (the absorption) and that of the inverse one (the stimulating emission) and to compare them. There is no principal difficulty.

The main problem of the experimental technique is that the absorption lines are usually very narrow, so that the Doppler effect hinders to measure their natural widths and cross-sections. In other cases, when the natural linewidth happens to be wider then the Doppler width, the lifetime of the excited states proves to be very small. It does not allow separating in time the processes of excitation and monitoring the excited states.

However, in the molecular physics there is a unique object, which is characterized by an unusual combination of properties. On the one hand, it has very large homogeneous spectral width of optical transition (several orders of magnitude greater than the Doppler width). On the other hand, it has a long lifetime of the excited states compared to that of the spontaneous emission ($`>1ms`$), that allows to separate in time the processes of the laser excitation and monitoring. Experiments with this object in late $`1980^{th}`$ had yielded surprising results, which were not given any theoretical explanation, the works with this object were stopped and for more than ten years nobody dared to return to it. Those experiments directly demonstrated T-invariance violation in interactions of the photons with molecules.

This unique object is the so-called wing of the spectral line or the wide component of the line . It appears as a certain wide continuum in the absorption spectrum of the medium size molecules (having usually 4–10 atoms). Its nature, apparently, is connected with some features of the vibrational motion of atoms in the molecules . It is important, that we have several reliable experimental proofs of the existence of the wide component. For example, from the results of works , where an excitation of the $`SF_6`$ molecules by radiation of a pulsed $`CO_2`$-laser in the conditions of a molecular jet was studied, it is possible to obtain precisely the absorption cross-section of the wide component $`\sigma =610^{20}cm^2`$ and to evaluate it’s lorentzian width ($`150GHz`$). Direct observation of the wide component as a continuum was conducted in the work , where absorption of radiation of the continuous wave $`CO_2`$-laser by the $`SF_6`$ molecules was studied in the conditions of the molecular beam. From the results of this work it is possible to evaluate the absorption cross-section of the continuum ($`10^{19}cm^2`$). Finally, the lorentzian shape of a continuum was confirmed by detection of the far wings of molecular absorption bands, which have a lorentzian shape and have intensity, corresponding to the intensity of the continuum near the center of the absorption band .

Thus, the continuum characterizes mainly the homogeneous absorption of molecules. And probing the excited states by radiation of the second $`CO_2`$-laser demonstrates presence in the absorption spectrum of a sharp dip with a width $`450kHz`$ . It is a typical indication of the heterogeneity of a spectrum. In late $`1980^{th}`$ this effect was not given any explanation. The only explanation is that the spectral width of the reverse optical transition should be $`310^5`$ times less, than that of the forward transition. Moreover, in fact, this ratio can be even much greater .

In the same conditions, amplification of the probe laser radiation was also observed. If the fluence of the first $`CO_2`$-laser radiation is less then $`20\mu j/cm^2`$, about 0,01% of the total $`SF_6`$ molecules are excited in the beam . Amplification of the probe laser radiation in this case is possible only if the cross-section of the reverse transition exceeds that of the direct one by a factor of $`>10^4`$. Thus, these experiments clearly show that the direct and the reverse processes differ very much in their cross-sections and spectral widths. However, integrated cross-sections of the direct and reverse processes (the Einstein coefficients) can be identical. Observations of the Rabi oscillations vote for the benefit of such an assumption.

The experiments, mentioned above, can be considered as a direct proof of the T-invariance violation. At the same time, numerous facts, accumulated in the atomic and molecular physics, can be regarded as indirect indications of this phenomenon. An existing semiclassical theory of the optical transitions  allows satisfactory description of the physical effects from the apparent, formal side. But on this basis it is practically impossible to get any concrete information on the physical nature of those processes, which underlie the observable effects.

Perhaps, most brightly this feature manifests itself in explanation of the population transfer effect at sweeping the resonance conditions. If the laser is turned precisely in a resonance with the optical transition in a two-level system, than, according to the theory, the population of levels should oscillate with Rabi frequency. If the wavelength of the laser radiation changes and passes through a resonance (sweeping of the resonance conditions), it can be regarded as decrease in the effective energy of the laser pulse. In this case, reduction in the number of the Rabi oscillations can be expected. However, only complete and irreversible transfer of population from one level on another in one sweeping cycle is observed in the experiments. Manipulations with vectors and series in the model of a rotating wave allow the theorists to declare that they can explain this effect . However, they cannot reveal the physical reason of the incipient asymmetry and irreversibility of the population transfer process. The difference in the spectral width of the forward and backward optical transitions quite suits as such a physical reason. In this case sweeping of the resonance conditions must lead in a natural and inevitable way to the complete population transfer.

A similar situation takes place in many other cases also. T-invariance violation, described above, is a good physical basis for explanation of the nature of such effects, as amplification without inversion , a phase conjugation , a photon echo , an electromagnetically induced transparency , a coherent population trapping , the Autler-Townes effect  and others.

An existing semiclassical theory of the optical transitions is based on the assumption about the complete T-invariance of the process of interaction of photons with atoms and molecules. Taking the fact of strong T-invariance violation demands rather radical revision of the theory. But more important now is to reanimate and to prolong the experiments with the line wings .