text:
# Space-time nonlinear compression and three-dimensional complex trapping in normal dispersion

## Abstract

In positive phase-mismatched SHG and normal dispersion, a gaussian spatio-temporal pulse transforms spontaneously into a X-pulse, underlies spatio-temporal compression and eventually leads to stationary 3-D propagation. Experimental and numerical data are provided.

X-pulses (the radial analogous of tilted pulses) are known from linear optics since they sustain propagation with constant spatio-temporal profile in normal dispersive media Sonajalg . Here we show that X pulses are formed spontaneously from a gaussian beam/pulse under SHG with positive mismatch ($`2k_\omega >k_{2\omega }`$) and normal dispersion.

This transformation leads to the relevant dynamics presented by the numerical results in Figure 1. The achieved spatio-temporal compression is a genuine effect of the space-time interplay. In fact, self-focusing and time broadening are expected for independent spatial and temporal dynamics. During compression the energy radiates out of the beam center, as shown by the drop in the energy content ($`\mathrm{E}=\tau \mathrm{d}^2\mathrm{I}`$, where $`\tau `$ and d are FWHM). After such transient the pulse keeps an almost constant 25 fs duration (GVM splitting lenght=0.7mm, GVD dispersion length=6mm, for a gaussian profile). Spatially, a very small residual diffraction takes place, which let the state to relax slowly toward the linear regime. Note here the constancy of the energy content (right), which proves that no radiation occurs neither in time (due to GVM) nor in space (due to off-axis components). We verified that a small perturbation of the beam shape introduces radiation losses. If we switch off the nonlinearity, than relevant beam diffraction and GVM pulse splitting occurs. Outside the crystal, the pulse broadens due to effective anomalous GVD. The scenario indicates the occurrence of a novel process, which we might call ”nonlinear diffusion”, that keeps, at any finite intensity, exact balance among the dispersions (at all orders) introduced by the material, the self-phase-modulation and the angular dispersion. Asymptotic analysis is in progress to see if the field eventually relax to the exact linear spatio-temporal eigenmode of the double-frequency field.

A preliminary experiment done in the same conditions verified the validity of the model. Figure 2 shows the measured FH and SH spatial profiles at the output of a 22mm LBO crystal, and the autocorrelation after 35 mm of further propagation in air. The agreement with the calculated profiles (see Fig. 3), without free parameters, is very satisfactory.

Preliminary calculations indicate that, in the absence of GVM, a periodical behavior occurs, consistently with the existence of a X-type pulse with finite envelope (e. g. with finite energy), which propagate without any relaxation. The intensity evolution of such a state is plotted as continuos line (no symbols) in Fig. 1. Note .


format: text

text:
# Grain Alignment in Molecular Clouds

## 1 Why do we care?

The fact that the grain got aligned has been known for more than half a century. Nevertheless, it has been always a puzzle why grains get aligned in interstellar medium. Very soon after the discovery of grain alignment by Hall (1949) and Hiltner (1949) it became clear that the alignment happens in respect to magnetic field. Since that time grain alignment stopped to be the issue of pure scientific curiosity, but became an important missing link of connecting polarimetry observations with the all-important interstellar magnetic fields111 Additional interest to grain alignment arises from recent attempts to separate the polarized CMB radiation from the polarized foregrounds (see Lazarian & Prunet 2002 for a review)..

The history of grain alignment ideas is excited (see review by Lazarian 2003) and here we do not have space here to dwell upon it. Within this very short review we will discuss the modern understanding of grain alignment processes applicable to molecular clouds. Last decade has been marked by a substantial progress in understanding new physics associated with grain alignment. The theory has become predictive, which enables researchers to interpret observational data with more confidence.

Recent theory reviews of the grain alignment theory include Roberge (1996), Lazarian (2000, 2003). Progress in testing theory is covered in Hildebrand (2000), while particular aspects of grain dynamics are discussed in Lazarian & Yan (2003). The presentation in Lazarian (2003) goes beyond molecular cloud environment and deals with the possibility of alignment in circumstellar regions, interplanetary medium, coma of comets etc. The interested reader may use these reviews to guide her in the vast and exciting original literature on grain alignment.

## 2 How do grains rotate?

Dynamics of grains in molecular clouds is pretty involved (see Fig. 1). First of all, grains rotate. The rotation can arise from chaotic gaseous bombardment of grain surface and be Brownian, or it can arise from systematic torques discovered by Purcell (1975, 1979). The most efficient among those are torques arising from H2 formation over grain surface. One can visualize those torques imagining a grain with tiny rocket nozzles ejecting nascent high velocity hydrogen molecules (see Fig. 1). Grains are known to enable atoms of hydrogen to form molecules. The reactions are believed to take place over particular catalytic sites on grain surface. Those catalytic sites act as ”Purcell rockets”. Even when the surroundings of dust grains is mostly molecular, grains can rotate suprathermally, i.e. with kinetic energies much larger that $`kT_{gas}`$, due to the variation of the accommodation coefficient. Indeed, if the temperatures of gas and dust are different, those variations allow parts of the grain to bounce back impinging gaseous atoms with different efficiencies. It is easy to understand that this also results in systematic torques. In addition, Purcell (1979) identified electron ejection as yet another process that can drive grain to very large angular velocities.

A very different process of grain spin-up can be found in a very important, but not timely appreciated work by Dolginov & Mytrophanov (1976). These authors considered differential scattering of photons of right and left circular polarization by an irregular dust grain. As the size of the irregularities gets comparable with the wavelength, it is natural that interaction of a grain with photons will depend on the photon polarization. Unpolarized light can be presented as a superposition of equal number of left and right circularly polarized photons. Therefore it is clear that the interaction with photons of a particular polarization would deposit angular momentum to the grain. The authors concluded that for typical diffuse ISM conditions this process should induce grain rotation at suprathermal velocities. However, while Purcell’s torques became a textbook stuff, radiative torques had to wait 20 years before they were reintroduced to the field (Draine 1996, Draine & Weingartner 1996, 1997).

The minimal rotational velocity of grain is the velocity of their Brownian motion. This can be characterized by a temperature that is somewhere between that of the gas and the dust. In the case of PAHs or very small grains emitting copious microwave radiation as they rotate (Draine & Lazarian 1998) the effective rotational temperature may be subthermal (see discussion in Lazarian & Yan 2003).

It was realized by Martin (1971) that rotating charged grains will develop magnetic moment and the interaction of this moment with the interstellar magnetic field will result in grain precession. The characteristic time for the precession was found to be comparable with $`t_{gas}`$. However, soon a process that renders much larger magnetic moment was discovered (Dolginov & Mytrophanov 1976). This process is the Barnett effect, which is converse of the Einstein-de Haas effect. If in Einstein-de Haas effect a paramagnetic body starts rotating during remagnetizations as its flipping electrons transfer the angular momentum (associated with their spins) to the lattice, in the Barnett effect the rotating body shares its angular momentum with the electron subsystem causing magnetization. The magnetization is directed along the grain angular velocity and the value of the Barnett-induced magnetic moment is $`\mu 10^{19}\omega _{(5)}`$ erg gauss-1 (where $`\omega _{(5)}\omega /10^5\mathrm{s}^1`$). Therefore the Larmor precession has a period $`t_{Lar}3\times 10^6B_{(5)}^1`$ s.

Nevertheless, suprathermal, thermal and subthermal grain rotation are just components of complex grain dynamics. As any solid body, interstellar grains can rotate about 3 different body axes. As the result they tumble while rotating. This effect was attracting attention of the early researchers (see Jones & Spitzer 1967) till Purcell (1979) identified internal relaxation within grains as the process that can suppress grain rotation about all axes, but the axis corresponding to the grain maximal moment of inertial (henceforth axis of maximal inertia). Indeed, consider a spheroidal grain, which kinetic energy can be presented as (see Lazarian & Roberge 1997) $`E(\theta )=\frac{J^2}{I_{max}}\left(1+\mathrm{sin}^2\theta (h1)\right),`$ where $`\theta `$ is the angle between the axis of major inertia and grain angular momentum. In the absence of external torques grain angular momentum is preserved. The minimum of grain energy corresponds therefore to $`\theta =0`$, or grain rotating exactly about the axis of maximal inertia. As internal dissipation decreases kinetic energy, it sounds natural that $`\theta =0`$ is the expected state of grain subjected to fast internal dissipation.

Purcell (1979) introduced a new process of internal dissipation which he termed ”Barnett relaxation”. This process may be easily understood. We know that a freely rotating grain preserves the direction of $`𝐉`$, while angular velocity precesses about $`𝐉`$ and in grain body axes. We learned earlier that the Barnett effect results in the magnetization vector parallel to $`\stackrel{}{\mathrm{\Omega }}`$. As a result, the Barnett magnetization will precess in body axes and cause paramagnetic relaxation. The “Barnett equivalent magnetic field”, i.e. the equivalent external magnetic field that would cause the same magnetization of the grain material, is $`H_{BE}=5.6\times 10^3\omega _{(5)}`$ G, which is much larger than the interstellar magnetic field. Therefore the Barnett relaxation happens on the scale $`t_{Bar}4\times 10^7\omega _{(5)}^2`$ sec, i.e. essentially instantly compared to the time that it takes to damp grain rotation for typical molecular cloud conditions.

Even stronger relaxation process has been identified recently by Lazarian & Draine (1999a). They termed it “nuclear relaxation”. This is an analog of Barnett relaxation effect that deals with nuclei. Similarly to unpaired electrons nuclei tend to get oriented in a rotating body. However the nuclear analog of “Barnett equivalent” magnetic field is much larger and Lazarian & Draine (1999a) concluded that the nuclear relaxation can be a million times faster than the Barnett relaxation.

Why would the actual relaxation rate matter? The rate of internal relaxation couples grain rotational and vibrational degrees of freedom. LD99b showed that this will result in grain “thermal flipping”. Such a flipping would average out Purcell’s torques and result in grain being “thermally trapped” in spite of the presence of uncompensated torques. Whether grain gets “thermally trapped” depends on its size (with the grains less than a critical size $`a_c`$ rotating thermally). While Barnett and inelastic relaxation (see also Lazarian & Efroimsky 1999) results in $`a_c`$ equal or less than $`10^5`$ cm, the nuclear internal relaxation provides $`a_c10^4`$ cm. This means that most grains rotate thermally in the presence of Purcell’s torques. The exception to this thermallization are radiative torques that are not fixed in grain coordinates. Such torques can spin-up dust in spite of thermal flipping.

## 3 What does align grains?

Paramagnetic Alignment

Davis-Greenstein (1951) mechanism (henceforth D-G mechanism) is based on the paramagnetic dissipation that is experienced by a rotating grain. Paramagnetic materials contain unpaired electrons which get oriented by the interstellar magnetic field $`𝐁`$. The orientation of spins causes grain magnetization and the latter varies as the vector of magnetization rotates in grain body coordinates. This causes paramagnetic loses at the expense of grain rotation energy. Note, that if the grain rotational velocity $`\stackrel{}{\mathrm{\Omega }}`$ is parallel to $`𝐁`$, the grain magnetization does not change with time and therefore no dissipation takes place. Thus the paramagnetic dissipation acts to decrease the component of $`\stackrel{}{\mathrm{\Omega }}`$ perpendicular to $`𝐁`$ and one may expect that eventually grains will tend to rotate with $`\stackrel{}{\mathrm{\Omega }}𝐁`$ provided that the time of relaxation $`t_{DG}`$ is much shorter than $`t_{gas}`$, the time of randomization through chaotic gaseous bombardment. In practice, the last condition is difficult to satisfy. For $`10^5`$ cm grains in the diffuse interstellar medium $`t_{DG}`$ is of the order of $`7\times 10^{13}a_{(5)}^2B_{(5)}^2`$s , while $`t_{gas}`$ is $`3\times 10^{12}n_{(20)}T_{(2)}^{1/2}a_{(5)}`$ s ( see table 2 in Lazarian & Draine 1997) if magnetic field is $`5\times 10^6`$ G and temperature and density of gas are $`100`$ K and $`20`$ cm-3, respectively. However, at the time when it was introduced ,in view of uncertainties in interstellar parameters, the D-G mechanism looked plausible.

The first detailed analytical treatment of the problem of D-G alignment was given by Jones & Spitzer (1967) who described the alignment of $`𝐉`$ using a Fokker-Planck equation. This approach allowed them to account for magnetization fluctuations within grain material and thus provided a more accurate picture of $`𝐉`$ alignment. The first numerical treatment of D-G alignment was presented by Purcell (1969). By that time it became clear that the D-G mechanism is too weak to explain the observed grain alignment. However, Jones & Spitzer (1967) noticed that if interstellar grains contain superparamagnetic, ferro- or ferrimagnetic (henceforth SFM) inclusions222The evidence for such inclusions was found much later through the study of interstellar dust particles captured in the atmosphere (Bradley 1994)., the $`t_{DG}`$ may be reduced by orders of magnitude. Since $`10\%`$ of atoms in interstellar dust are iron the formation of magnetic clusters in grains was not far fetched (see Martin 1995). However, detailed calculations in Lazarian (1997), Roberge & Lazarian (1999) showed that the alignment achievable cannot account for observed polarization coming from molecular clouds provided that dust grains rotate thermally. This is the consequence of thermal fluctuations within grain material. These internal magnetic fluctuations randomize grains orientation in respect to magnetic field if grain body temperature is close to the rotational temperature.

Purcell (1979) pointed out that fast rotating grains are immune to both gaseous and internal magnetic randomization. Thermal trapping limits the range of grain sizes for which Purcell’s torques can be efficient (Lazarian & Draine 1999ab). Grains with the sizes larger than the wavelength size can be spun up by the incoming starlight, however (see Draine & Weingartner 1996).

Mechanical Alignment

Gold (1951) mechanism is a process of mechanical alignment of grains. Consider a needle-like grain interacting with a stream of atoms. Assuming that collisions are inelastic, it is easy to see that every bombarding atom deposits angular momentum $`\delta 𝐉=m_{atom}𝐫\times 𝐯_{atom}`$ with the grain, which is directed perpendicular to both the needle axis $`𝐫`$ and the velocity of atoms $`𝐯_{atom}`$. It is obvious that the resulting grain angular momenta will be in the plane perpendicular to the direction of the stream. It is also easy to see that this type of alignment will be efficient only if the flow is supersonic333Otherwise grains will see atoms coming not from one direction, but from a wide cone of directions (see Lazarian 1997a) and the efficiency of alignment will decrease.. Thus the main issue with the Gold mechanism is to provide supersonic drift of gas and grains. Gold originally proposed collisions between clouds as the means of enabling this drift, but later papers (Davis 1955) showed that the process could only align grains over limited patches of interstellar space, and thus the process cannot account for the ubiquitous grain alignment in diffuse medium.

Suprathermal rotation introduced in Purcell (1979) persuaded researchers that mechanical alignment is marginal. Indeed, fast rotation makes it difficult for gaseous bombardment to align grains. However, two new developments must be kept in mind. First of all, a number of papers proved that mechanical alignment of suprathermally rotating grains is possible (Lazarian 1995, Lazarian & Efroimsky 1996, Efroimsky 2002). Moreover, recent work on grain dynamics (Lazarian & Yan 2002, Yan & Lazarian 2003) proved that MHD turbulence can render grains with supersonic velocities. While we do not believe that mechanical alignment is the dominant process, it should be kept in mind in analyzing observations (see Rao et al. 1998).

Alignment via Radiative Torques

Anisotropic starlight radiation can both spin the grains and align them. This was first realized by Dolginov & Mytrophanov (1976), but this work came before its time. The researchers did not have reliable means to study dynamics of grains and the impact of their work was marginal. Before Bruce Draine realized that the torques can be treated with the versatile discrete dipole approximation (DDA) code ( Draine & Flatau 1994) the radiative torque alignment was very speculative. For instance, earlier on difficulties associated with the analytical approach to the problem were discussed in Lazarian (1995a). However, very soon after that Draine (1996) modified the DDA code to calculate the torques acting on grains of arbitrary shape. His work revolutionized the field! The magnitude of torques were found to be substantial and present for grains of various irregular shape (Draine 1996, Draine & Weingartner 1996). After that it became impossible to ignore these torques.

One of the problem of the earlier treatment was that in the presence of anisotropic radiation the torques will change as the grain aligns and this may result in a spin-down. Moreover, anisotropic flux of radiation will deposit angular momentum which is likely to overwhelm rather weak paramagnetic torques. These sort of questions were addressed by Draine & Weingartner (1997) and it was found that for most of the tried grain shapes the torques tend to align $`𝐉`$ along magnetic field. The reason for that is yet unclear and some caution is needed as the existing treatment ignores the dynamics of crossovers which is very important for the alignment of suprathermally rotating grains. A recent work by Weingartner & Draine (2003) treats flipping of grains in the presence of monochromatic radiation.

### 3.1 What is Future Work?

Observational testing of alignment is extremely important. Both the dependences of the polarization degree versus wavelength that follow Serkowski law (Serkowski 1973) and studies of changes of polarization degree with the wavelength done in Far Infrared (see Hildebrand 2000) are consistent with theoretical predictions (see discussion in Lazarian, Goodman & Myers 1997). According to Lazarian (2003) the study of grain alignment at the diffuse/dense cloud interface by Whittet et al. (2001) is suggestive that grains there are being aligned by radiative torques.

Radiative torques look as the most attractive mechanism to align grains in molecular clouds. However more theoretical work is required to understand why grains subjected to anisotropic radiation get preferentially aligned with their long axes perpendicular to magnetic field.

Acknowledgments. I thank H. Yan for suggestions and help. The research is supported by the NASA grant 0830 300 N665 736.

## Index


text:
# The Proximity Effect Around High Redshift Galaxies

## 1 Introduction

In a recent study, Adelberger et al (2003, hereafter A03) have analyzed a sample composed of 8 high-resolution QSOs spectra at $`3.1<z<4.1`$ and spectroscopic redshifts for 431 Lyman-break galaxies (LBGs) at slightly lower $`z`$, to study the intergalactic medium (IGM) close to high-$`z`$ galaxies and detect imprints of galaxy-IGM interplay. Comparing the location of galaxies to the absorption lines in QSOs spectra, the authors found that the mean Ly$`\alpha `$ transmitted flux generally increases with the line-of-sight impact parameter between the foreground galaxies and the background QSOs ($`\mathrm{\Delta }r`$), reaching the mean transmissivity of 0.67 at $`\mathrm{\Delta }r\stackrel{>}{}6`$ $`h^1`$ comoving Mpc; hovever, the trend is reversed within about $``$ $`1h^1`$ comoving Mpc of the galaxy. This result indicates that the IGM contains less neutral hydrogen close to LBGs, and hence is more transparent to Ly$`\alpha `$ photons, than at large distance from them.
Similar observations at lower redshift show an opposite behavior. Lanzetta et al (1995), Chen et al (1998) and Pascarelle et al (2001) find that galaxies, at $`z<1`$, with an impact parameter $`\mathrm{\Delta }r\stackrel{<}{}180h^1`$ comoving kpc are much more often associated with Ly$`\alpha `$ absorption lines than those at larger impact parameters. Furthermore the result of VLT/UVES study of the Ly$`\alpha `$ forest in the vicinity of the LBG MS1521-cB58 at $`z=2.724`$ (Savaglio, Panagia & Padovani 2002) shows an absorption excess close to the galaxy.

A03 suggest that the relative lack of neutral hydrogen at small distances from LBGs is possibly due to the dynamical feedback from supernova-driven winds, which may be able to accelerate the surrounding gas up to velocities exceeding the escape velocity of the galaxy potential well, driving plausibly the intergalactic material out to a distance $`>0.5h^1`$ comoving Mpc.

This scenario has been investigated with numerical simulations by different groups: Croft et al (2002), Bruscoli et al (2003), and Kollmeier et al (2003). Despite the completely independent prescriptions adopted to model supernova-driven winds in these studies, none of them were able to reproduce A03 results on the entire distance range. Unless feedback schemes are still far too crude in SPH simulations, the discrepancy indicates that the energy injection from supernova explosions is not able to displace the IGM around galaxies to distances $`0.5h^1`$ comoving Mpc. The most likely explanation is that the ram-pressure of the accreting/infalling gas slows down and effectively confines winds to relatively modest distances from the explosion sites (Bruscoli et al. 2003).

Local photoionization by galaxies themselves provide the only known alternative mechanism able to decrease the neutral hydrogen fraction in the LBG environment. The first attempts to account for local photo-ionization in numerical simulations have been performed using simplified analytical techniques and/or ad-hoc prescriptions (A03, Kollmeier et al. 2002, Croft et al. 2002); these studies have lead to conclude that this process cannot account for the observed IGM transparency. However, fully self-consistent and detailed radiative transfer calculations are required before a final conclusion is drawn.

In this Letter we investigate the role of local photoionization using high resolution numerical simulations performed with the radiative transfer code CRASH ; the full description of the code and of the tests made are given in Maselli, Ferrara & Ciardi (2003).

## 2 Simulations

We have performed radiative transfer simulations in order to evaluate self-consistently the proximity effect due to the ionizing radiation emitted by high redshift galaxies on the surrounding IGM. We use the numerical code CRASH (Maselli et al. 2003) which allows to calculate the effects of ionizing radiation on a precomputed cosmological density field, specifically obtained from SPH simulations fully described elsewhere (Marri et al. 2003) and already used by Bruscoli et al. (2003). In the following we give the details of the two simulations.

### 2.1 Hydrodynamics

The simulation has been carried out with an improved version of the SPH code GADGET2 (Springel et al. 2001), designed for the detailed study of galaxy formation (Marri & White 2002). Such multi-phase (MSPH) scheme is based on an explicit separation of protogalactic gas into diffuse and dense (star forming) components and includes new prescriptions for a realistic treatment of stellar feedbacks. With respect to the standard GADGET2 version, this new implementation is able to suppress star formation, to reheat cold gas and to drive outflows from galactic disks. The presence of a UV background (UVB) produced by QSOs and filtered through the IGM is also included. The shape and amplitude of the UVB are taken from Haardt & Madau (1996) and its impact on the IGM temperature and ionization structure is calculated as in Katz et al. (1996).

The simulation uses $`128^3`$ particles in a 10.5$`h^1`$ comoving Mpc box and assumes a $`\mathrm{\Lambda }`$CDM cosmological model with $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_b=0.04`$ and $`h=0.7`$; the initial density perturbations power spectrum is cluster-normalized ($`\sigma _8=0.9`$) and periodic boundary conditions are assumed. Galaxies are identified by applying the HOP group finding algorithm (Eisenstein & Hut 1998) to the baryonic components of the simulation. Each galaxy is assigned a star formation rate (SFR) calculated by considering the properties of all particles within the group. For each of them we first check if the conditions under which a gas particle can form stars are satisfied and we calculate the rate at which its gas is converted into stars. For more details we defer the reader to Marri & White (2002).

For consistency sake, we use the same output as in Bruscoli et al. (2003). There we showed that, although on large scales the simulation reproduces a number of statistical properties of the IGM, it does not match the transmitted Ly$`\alpha `$ flux inside 1 $`h^1`$ comoving Mpc as measured by A03. Such an output has now been post-processed using CRASH to self-consistently include the effects of local photoionization produced by the galaxies identified in the computational volume.

### 2.2 Radiative Transfer

The radiative tranfer (RT) simulation has been performed on the outputs at $`z=3.27`$ of the MSPH simulation above, to allow a direct comparison to the observational data available around this redshift (A03).

The RT is simulated running the 3-D RT code CRASH , whose description is fully reported in Maselli et al. (2003). The code works on precomputed H/He density fields assigned in a 3-D cartesian grid, and calculates the time evolution of the temperature and of the ionization structure.

To include the effect of local photoionization on the cosmological configuration obtained from the MSPH simulation, it is necessary to map the particle data on a 3-D cartesian grid. We choose a grid composed of 1283 cells, which corresponds to a spatial resolution of $`\mathrm{\Delta }x82`$ kpc $`h^1`$ comoving. At each grid point we assign a value for the mass density ($`\rho `$), temperature ($`T`$), neutral hydrogen fraction ($`x_{\mathrm{H}^0}=n_{\mathrm{H}^0}/n_\mathrm{H}`$) and free electrons fraction ($`x_e=n_\mathrm{e}/(n_\mathrm{H}+n_{\mathrm{He}})`$), performing a standard SPH smoothing using the 32 neighbour SPH particles closer to the grid point. We next locate the 398 galaxies identified inside the simulation box on the cartesian grid, converting the lagrangian spatial coordinates into cartesian coordinates.
The emission properties of the galaxies have been derived using the Starburst99 code (Leitherer et al. 1999), a web based software and data package designed to model spectrophotometric and related properties of star-forming galaxies. We have assumed a continuous star formation law and a Salpeter IMF with $`\alpha =2.35`$, normalized in a mass range $`M_{}`$; we also take Z/5 as a characteristic metallicity for all the sample. These values are suggested by the analysis of the LBGs UV continuum and the Hα and Hβ emission lines by Pettini et al. (2001). As the physical connection of LBGs with galaxy type/mass is unclear (Primack et al. 2003), we assume for all galaxies in the box similar scalings. Running Starburst99 under the above conditions, we obtain the galaxy SED corresponding to the derived SFR of each galaxy. We finally assume a value for the ionizing photon escape fraction of 50%. This value represents the upper limit of the allowed range (see discussion in Ciardi, Bianchi & Ferrara 2002) from a set of observational data.

The computed grid data are assigned as initial conditions for the CRASH run. The number densities of hydrogen and helium nuclei are calculated assuming mass fractional abundances $`f_H`$=0.76 and $`f_{He}`$=0.24, respectively. Temperature and hydrogen ionization fractions are initialized as in the MSPH simulation output; for helium we assume ionization fractions $`x_{\mathrm{He}^+}=x_{\mathrm{H}^+}`$ and $`x_{\mathrm{He}^{++}}=(x_ex_{\mathrm{He}^+}x_{\mathrm{H}^+})/2`$. The physical duration of a simulation run has been set to $`t_s=6\times 10^6`$ yrs, i.e. several times the estimated box-crossing-time for the ionization front produced by a typical source in the volume.

Nearby galaxies the mean gas temperature is above 105 K, due to the energy injection into the IGM produced by outflows. The cooling time of such gas in the simulation is longer than $`t_s`$, with a minimum value equal to 1.5$`\times `$108 yrs. As the temperature does not change significantly during the RT simulations, its evolution has not been solved for to reduce the computational time.

Two different RT simulations have been run. In the first one only the UVB radiation is included (run A), while in the second one, the radiation from all galaxies in the simulation volume is added to the UVB (run B). In this way it is possible to compare consistenly the two simulations to infer the role of galaxies and distinguish the contribution to the proximity effect due to local photoionization. As considerable uncertainty is present in the determination of SFRs of LBGs, we also studied a case similar to run B in which however SFRs are boosted ad hoc to match the observed values.

The total ionizing energy, injected in the computational volume, is discretized in $`N_p=2\times 10^8`$ photon packets for run A and in $`N_p=6\times 10^8`$ packets for run B: roughly, this corresponds to $`10^6`$ photon packets per galaxy and $`2\times 10^8`$ photon packets for the UVB.

## 3 Results

As an example, Fig.1 shows a comparison among the final outputs of run A (left panel) and run B (central and right panels). The slices represent the logarithm of the neutral hydrogen fraction on a slice of a zoomed region of about 4$`h^1`$ comoving Mpc by side, centered on the most massive galaxy in the computational volume. The two maps from run B have been derived by assigning to such a galaxy different SFR: in the central panel the SFR is about 30 M yr-1 as derived from the MSPH simulation, while in the right panel it has been boosted by a factor of 10. The effect of the local photoionization produced by the galaxy is clearly visible in both cases and depends strongly on its luminosity, as expected. In the central panel the imprint of the galaxy on the surrounding gas is feeble and the content of neutral hydrogen in a surrounding region of 0.5 h-1 comoving Mpc by side, is reduced on average by factor of 2. The luminosity associated with a SFR of 30 M yr-1 is too low to affect significantly the high density environment of such a massive galaxy (M =$`1.5\times 10^{11}`$M$`{}_{}{}^{}h_{}^{1}`$). When the SFR is boosted by a factor of 10 (right panel), the region affected by the galaxy UV flux extends up to a distance of about 1$`h^1`$ comoving Mpc. In the inner region the neutral hydrogen fraction is depleted at least by an order of magnitude. In the outer regions the effect of the local photoionization is progressively reduced and vanishes at distances $`1h^1`$ Mpc, where the UVB dominates the local UV flux. The decrease in the neutral hydrogen fraction visible in other regions far from the center is produced by other galaxies, included in run B as ionizing sources.

To estimate the radius of the “sphere of influence” produced by the local emission from a typical galaxy we study the distribution of the ratio between the UV flux produced by the galaxies and the UVB flux, $`=F_{gal}/F_{bkg}`$, in the box. Fig.2 shows the histogram derived from the final output of run B. The distribution has a mean value $`\overline{}=2.64\times 10^2`$ and a median value $`_M=1.83\times 10^2`$. These values are to be considered as lower limits to the contribution of galaxies to the UVB, as larger boxes would be required to properly estimate such quantity. The volume in which $`>1`$, i.e. where the local emission by galaxies dominates the foreground UVB, is about $`4.8\%`$ of the total computational volume. This fraction corresponds to a mean radius of influence for a typical galaxy $`R_i`$, given by:

$`R_i=\left({\displaystyle \frac{4}{3\pi }}{\displaystyle \frac{N_c(>1)\times V_c}{N_g}}\right)^{1/3}0.4h^1\mathrm{Mpc}`$ (1)
where $`N_c`$ is the number of cells where $`>1`$, $`V_c`$ is the volume of a single cell and $`N_g=398`$ is the number of galaxies in the simulation box.

To allow a direct comparison of the simulation results with the observational data we construct synthetic Ly$`\alpha `$ spectra as described in Bruscoli et al. , by tracing random LOS through the box. We then estimate the mean Ly$`\alpha `$ transmitted flux as a function of the impact parameter $`\mathrm{\Delta }r`$, computing the average on all the pixels at a distance from a galaxy in a given interval centered on $`\mathrm{\Delta }r`$. The results are plotted in Fig.3 and compared with the data, represented by the black points.

The three panels show the mean Ly$`\alpha `$ flux computed on different samples of galaxies: all the 398 galaxies (left panel), galaxies with mass above $`2\times 10^{10}`$ M (central panel) and those with mass below $`9.3\times 10^8`$ M (right panel).

Solid, dotted and dotted-dashed lines are derived from the outputs of run A, run B and from a simulation analogous to run B, but with SFRs boosted ad hoc. The difference between run A and run B in the left panel is marginal: we conclude that the mean effect of the local photoionization on the Ly$`\alpha `$ transmitted flux is negligible when all galaxies are included in the analysis.

Massive galaxies are the best candidates for LBGs due to their high luminosity and clustering properties. The SFR derived for the most massive galaxies in the MSPH simulation are in the range $`10÷30`$ M$`{}_{}{}^{}\mathrm{yr}_{}^{1}`$, the highest values in the simulation. Due to their high luminosity, one would expect them to produce a strong impact on the ionization of the surrounding gas. However this is not the case. Looking at the solid line in the central panel, obtained by neglecting the local emission, we note that the high density characterizing the environment of massive galaxies suppresses the mean transmitted flux with respect to the mean trend obtained for all the galaxies in the simulation (solid line in the left panel). The UV radiation emitted by these galaxies, using the SFR derived from the MSPH simulation (dotted line in central panel), is not strong enough to enhance significantly the transparency of the surrounding IGM, because of the high recombination rate in the denser environment. Higher values of SFR are necessary to produce a significant increase of the transmitted flux nearby massive galaxies. This can be seen from the boosted case, i.e. with SFR in the range $`100÷300`$ M yr-1. Although on the high side, the values are still plausible according to observations (Shapley et al. , 2001). Despite the sharp increase in the mean Ly$`\alpha `$ transmitted flux, yet this simulation does not match the A03 data. This result suggests that the mean transmissivity strongly depends on the galactic enviroment. The right panel shows the trend for the Ly$`\alpha `$ transmitted flux obtained selecting galaxies with mass $`9.3\times 10^8`$ M. In this case the transmissivity of the gas at $`\mathrm{\Delta }r<0.7h^1`$ Mpc has an opposite trend with respect to the mean (left panel), and seems to follow the one observed by A03. The increasing flux at smaller distance from the galaxy reflects the lower density in these regions. Here the UVB flux can counteract the recombinations which occur on longer time scales. Nevertheless the SFRs associated with these galaxies are too small (of order $`0.1M_{}`$ yr-1) to further increase the transmitted flux. Boosting the SFR up to values around 50 M yr-1 (dashed-dotted line), the increase in $`F`$ is significant and closely approches the observed value in the innermost region. By no means though the flux around 1 Mpc can be matched.

## 4 Conclusions

A recent study on the statistics of the Ly$`\alpha `$ forest in the vicinity of foreground galaxies has revealed an unexpected relative lack of neutral hydrogen in the inner Mpc $`h^1`$ comoving from LBGs at $`z3`$. High energy supernova-driven winds able to displace the surrounding gas out to distances greater than $`0.5h^1`$ have been proposed as a viable explanation. However, several numerical studies failed in reproducing the observed effect of such superwinds via (M)SPH cosmological simulations. According to the simulations, the velocity of the gas closest to galaxies varies around a characterisitc value of 40 km s-1 for small objects, and around 150 km s-1 for the biggest ones. These values seem to underestimate the high velocity gas, at 600 km s-1, measured by A03 at a scale of a few kpc from the galaxy center. Nevertheless, the possibility for the wind to break through the galactic halo, and into the IGM up to distances greater than $`1h^1`$ comoving Mpc (even at such high velocity) seems quite unphysical if one takes into account the kinetic energy used to counteract the pressure of the accretion flow.

Motivated by these reasons, we have studied the role of local photoionization in determining the Ly$`\alpha `$ absorbers statistics in the vicinity of galaxies, via numerical radiative transfer simulations. We post-processed an output at $`z=3.27`$ of a multiphase SPH simulation, in which a consistent treatment of supernova feedback is implemented, running the code CRASH to account for the impact of the UV flux emitted by galaxies on the simulated IGM. We have derived the synthetic mean Ly$`\alpha `$ transmitted flux, as a function of the distance from a typical galaxy, for different samples of galaxies. We found that on average local photoionization has a negligible impact, i.e. when all the galaxies present in the box are considered in the flux distribution derivation.

Local phoionization can affect significantly the transparency of the surrounding IGM only for particular conditions of the galactic environment. We find that in order to reproduce the trend inferred by A03, the gas around LBGs should be less dense than that surrounding massive galaxies seen in the simulation. The transmissivity as a function of distance is roughly reproduced only if LBGs are identified with dwarf galaxies (with masses $`\stackrel{<}{}10^9M_{}`$) which are undergoing a vigorous ($`50M_{}`$ yr-1) burst of star formation. Similar and higher star formation rates in larger galaxies are not sufficient to overwhelm the large recombination rates that are occurring in their denser environment. This conclusion is somewhat at odd with the common wisdom that LBGs are predominantly massive systems, although different and equally valid arguments tend to identify them with the bursting population of dwarfs that we find necessary to explain the observed gas transparency. LBGs can be interpreted as outbursting dwarf galaxies (Bullock et al. 1999) if their abundance is due to an increased collision rate at high redshift; although many of the simulated collisions have relatively small masses ($`\stackrel{<}{}10^{10}M_{}`$), they tend to cluster about large-mass halos. They therefore exhibit strong clustering, similar to that observed.

The solve the puzzle, more data are needed. The poor statistics of A03 sample at small impact parameters, errors in measured redshifts and galaxy positions make the derived transmitted flux quite uncertain.


text:
# REFERENCES

## REFERENCES

* Alcock, C., et al.  1999, Astrophys. J., 521, 602.
* Evans, D.W., and Irwin, M. 1995, MNRAS, 277, 820.
* Eyer, l., and Woźniak, P.R. 2001, MNRAS, 327, 601.
* Kuijken, K., et al.  2002, The Messenger, 110, 15.
* Melchior, A.-L., et al.  1999, Astron. Astrophys. Suppl. Ser., 134, 377.
* Tomaney, A.B., and Crotts, A.P.S. 1996, Astron. J., 112, 2872.
* Udalski, A., Paczyński, B., Żebruń, K., Szymański, M., Kubiak, M., Soszyński, I., Szewczyk, O., Wyrzykowski, Ł., and Pietrzyński, G. 2002, Acta Astron., 52, 1.


text:
# Highly nonlinear dynamics of third-harmonic generation by focused beams

## I Introduction

We study the dynamics of two-color beams in a nonlinear isotropic medium. We take the initial conditions to be a one-color beam with Gaussian profile. With a cubic ($`\chi ^{(3)}`$) nonlinearity, which is the lowest order possible in an isotropic medium, the possible nonlinearities are THG, Kerr, and Raman. Any other nonlinearity requires either a different susceptibility (second-harmonics, high-harmonics), or more than two slowly-varying envelopes (sum and difference frequency generation). This is one of the simplest nonlinear optics problems, and is for that reason important; it has an application to third-harmonic generation (THG) microscopy thg\_microscopy .

The usual model of THG takes the fundamental beam to be Gaussian, and (inconsistently, or as an approximation) has THG as the sole nonlinearity Ward.1969 ; Boyd.1992 . This model has exact analytic solution in which the third-harmonic (TH) beam takes a Gaussian profile. For a phase-mismatch that is zero or negative, the energy in the TH peaks at the focus of the fundamental beam, and, after the focus, destructive interference causes all the TH to be re-absorbed by the fundamental. The latter is often said to be due to the Guoy shift Boyd.1992 , the phase-shift of $`\pi `$ that a Gaussian beam experiences in going from a far-field, through a focus, to a far-field. But, in fact, the details show it to depend on the beam’s phase and shape everywhere along its path. However, the fast (electronic) nonlinear response of optical materials that yields THG also generates Kerr effects (self- and cross-phase modulation). Therefore, one should include Kerr and perhaps also Raman effects when modeling THG. These can have quite drastic effects on a beam Berge.1998 . Although it was recognized as early as 1973 that Kerr effects could influence THG Miles.1973 , and although Kerr, THG, and dispersion have been examined in studies without transverse spatial dynamics Band.1990 ; SammutBuryakKivshar.1998 , the results of transverse nonlinear dynamics on THG have not heretofore been studied in quantitative detail.

We have derived a set of coupled nonlinear Schrödinger (NLS) equations for two slowly-varying envelopes, with the relevant nonlinearities treated rigorously and consistently TasgalBand.nls\_derivation . Because of the microscopy application, we are interested in very tightly focused beams and pulses, with large momentum and frequency spreads. The NLS equations are thus given to all orders in dispersion and diffraction; numerical simulations are carried out with a method that is accurate to all orders relevant for the grid Trippenbach.2002 ; Rothenberg.1992 ; FleckMorrisFeit.1976 . We express the field as a fundamental, $`A_\omega (𝐱,t)`$, centered about a carrier wave at frequency $`\omega _0`$, and a TH, $`A_{3\omega }(𝐱,t)`$, centered about a carrier wave at frequency $`3\omega _0`$,

$`i_zA_\omega (𝐱,t)`$ $`=`$ $`\left[{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{i^p}{p!}}\left({\displaystyle \frac{^p}{\omega ^p}}\beta (\omega )\right)_{\omega _0}_t^p{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(2m3)!!}{(2)^mm!}}{\displaystyle \frac{\mathrm{}_{}^m}{\beta (\omega _0)^{2m1}}}{\displaystyle \underset{p,m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{i^p(2m3)!!}{(2)^mm!p!}}\left({\displaystyle \frac{^p}{\omega ^p}}{\displaystyle \frac{1}{\beta (\omega )^{2m1}}}\right)_{\omega _0}\mathrm{}_{}^m_t^p\right]A_\omega `$ (1a)
$`+{\displaystyle \underset{p,m=0}{\overset{\mathrm{}}{}}}2\pi {\displaystyle \frac{i^p(2m1)!!}{(2)^mm!p!}}\left({\displaystyle \frac{^p}{\omega ^p}}{\displaystyle \frac{(\omega /c)^2}{\beta (\omega )^{2m+1}}}\right)_{\omega _0}\mathrm{}_{}^m_t^pP_\omega (𝐱,t),`$
$`i_zA_{3\omega }(𝐱,t)`$ $`=`$ $`\left[{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{i^p}{p!}}\left({\displaystyle \frac{^p}{\omega ^p}}\beta (\omega )\right)_{3\omega _0}_t^p{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(2m3)!!}{(2)^mm!}}{\displaystyle \frac{\mathrm{}_{}^m}{\beta (3\omega _0)^{2m1}}}{\displaystyle \underset{p,m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{i^p(2m3)!!}{(2)^mm!p!}}\left({\displaystyle \frac{^p}{\omega ^p}}{\displaystyle \frac{1}{\beta (\omega )^{2m1}}}\right)_{3\omega _0}\mathrm{}_{}^m_t^p\right]A_{3\omega }`$ (1b)
$`+{\displaystyle \underset{p,m=0}{\overset{\mathrm{}}{}}}2\pi {\displaystyle \frac{i^p(2m1)!!}{(2)^mm!p!}}\left({\displaystyle \frac{^p}{\omega ^p}}{\displaystyle \frac{(\omega /c)^2}{\beta (\omega )^{2m+1}}}\right)_{3\omega _0}\mathrm{}_{}^m_t^pP_{3\omega }(𝐱,t).`$
In the NLS equations (1), $`\mathrm{}_{}_x^2+_y^2`$ is the transverse Laplacian, $`n(\omega )`$ is the index of refraction at frequency $`\omega `$, and $`\beta (\omega )n(\omega )\omega /c`$ is the wavenumber. The first summation terms on the right-hand sides (RHS) of Eqs. (1) are dispersion, the second summation is diffraction, the third summation is cross-dispersion/diffraction. The nonlinear polarization includes THG, Kerr, and Raman effects SammutBuryakKivshar.1998 ; Boyd.1992 ; Agrawal.1995 . Self-steepening terms steepening\_refs ; Trippenbach.1997 , of the form $`(|A|^2A)/t`$ and self-frequency shifting terms of the form $`\frac{|A|^2}{t}A`$ are contained in the first-order time derivatives of the nonlinear polarizations $`P_\omega (𝐱,t)`$ and $`P_{3\omega }(𝐱,t)`$ in Eqs. (1). The NLS equations (1) show that self-steepening terms and self-frequency shifting terms are just the first terms of a family of higher order nonlinear terms.

The nonlinear polarization is taken to be of the form

$`P_\omega (𝐱,t)`$ $`=`$ $`3\chi ^{\mathrm{elec}}(\omega _0;\omega _0,\omega _0,3\omega _0)\mathrm{exp}\{i[3\beta (\omega _0)\beta (3\omega _0)]z\}A_\omega (𝐱,t)^2A_{3\omega }(𝐱,t)`$ (2a)
$`+`$ $`3\chi ^{\mathrm{elec}}(\omega _0;\omega _0,\omega _0,\omega _0)|A_\omega (𝐱,t)|^2A_\omega (𝐱,t)+6\chi ^{\mathrm{elec}}(\omega _0;3\omega _0,3\omega _0,\omega _0)|A_{3\omega }(𝐱,t)|^2A_\omega (𝐱,t)`$
$`+`$ $`A_\omega (𝐱,t){\displaystyle _0^{\mathrm{}}}[3\chi ^{\mathrm{nucl}}(\omega _0;\omega _0,\omega _0,\omega _0;s)|A_\omega (𝐱,ts)|^2+3\chi ^{\mathrm{nucl}}(\omega _0;\omega _0,3\omega _0,3\omega _0;s)|A_{3\omega }(𝐱,ts)|^2]𝑑s`$
$`+`$ $`A_{3\omega }(𝐱,t){\displaystyle _0^{\mathrm{}}}\mathrm{exp}(2i\omega _0s)3\chi ^{\mathrm{nucl}}(\omega _0;3\omega _0,\omega _0,3\omega _0;s)A_\omega (𝐱,ts)A_{3\omega }(𝐱,ts)^{}𝑑s`$
$`P_{3\omega }(𝐱,t)`$ $`=`$ $`\chi ^{\mathrm{elec}}(3\omega _0;\omega _0,\omega _0,\omega _0)\mathrm{exp}\{i[3\beta (\omega _0)\beta (3\omega _0)]z\}A_\omega (𝐱,t)^3`$ (2b)
$`+`$ $`6\chi ^{\mathrm{elec}}(3\omega _0;\omega _0,\omega _0,3\omega _0)|A_\omega (𝐱,t)|^2A_{3\omega }(𝐱,t)+3\chi ^{\mathrm{elec}}(3\omega _0;3\omega _0,3\omega _0,3\omega _0)|A_{3\omega }(𝐱,t)|^2A_{3\omega }(𝐱,t)`$
$`+`$ $`A_{3\omega }(𝐱,t){\displaystyle _0^{\mathrm{}}}[3\chi ^{\mathrm{nucl}}(3\omega _0;\omega _0,\omega _0,3\omega _0;s)|A_\omega (𝐱,ts)|^2+3\chi ^{\mathrm{nucl}}(3\omega _0;3\omega _0,3\omega _0,3\omega _0;s)|A_{3\omega }(𝐱,ts)|^2]𝑑s`$
$`+`$ $`A_\omega (𝐱,t){\displaystyle _0^{\mathrm{}}}\mathrm{exp}(2i\omega _0s)\mathrm{\hspace{0.17em}3}\chi ^{\mathrm{nucl}}(3\omega _0;\omega _0,\omega _0,3\omega _0;s)A_\omega (𝐱,ts)^{}A_{3\omega }(𝐱,ts)𝑑s.`$
This breaks up the nonlinear response into an electronic (fast) part and a nuclear (slow) part. It also assumes that the electronic part of the response may be considered frequency-independent on the scale of the pulse bandwidths; it does not make this assumption for the nuclear part of the response. For the calculations, we take a more specific nonlinear polarization, essentially the standard model of the third-order susceptibility of fused silica Hellwarth.1977 ; Stolen.1989 ; Blow.1989 ; Agrawal.1995  plus a generalization,

$`P_\omega (𝐱,t)`$ $`=`$ $`\mathrm{exp}\{i[3\beta (\omega _0)\beta (3\omega _0)]z\}\chi ^{\mathrm{THG}}A_\omega (𝐱,t)^2A_{3\omega }(𝐱,t)+\chi ^{\mathrm{elec}}(|A_\omega (𝐱,t)|^2+2|A_{3\omega }(𝐱,t)|^2)A_\omega (𝐱,t)`$ (3a)
$`+A_\omega (𝐱,t){\displaystyle _0^{\mathrm{}}}\chi ^{\mathrm{nucl}}(s)(|A_\omega (𝐱,ts)|^2+|A_{3\omega }(𝐱,ts)|^2)𝑑s`$
$`+A_{3\omega }(𝐱,t){\displaystyle _0^{\mathrm{}}}\mathrm{exp}(2i\omega _0s)\chi ^{\mathrm{nucl}}(s)A_\omega (𝐱,ts)A_{3\omega }(𝐱,ts)^{}𝑑s,`$
$`P_{3\omega }(𝐱,t)`$ $`=`$ $`\mathrm{exp}\{i[3\beta (\omega _0)\beta (3\omega _0)]z\}{\displaystyle \frac{1}{3}}\chi ^{\mathrm{THG}}A_\omega (𝐱,t)^3+\chi ^{\mathrm{elec}}(2|A_\omega (𝐱,t)|^2+|A_{3\omega }(𝐱,t)|^2)A_{3\omega }(𝐱,t)`$ (3b)
$`+A_{3\omega }(𝐱,t){\displaystyle _0^{\mathrm{}}}\chi ^{\mathrm{nucl}}(s)(|A_\omega (𝐱,ts)|^2+|A_{3\omega }(𝐱,ts)|^2)𝑑s`$
$`+A_\omega (𝐱,t){\displaystyle _0^{\mathrm{}}}\mathrm{exp}(2i\omega _0s)\chi ^{\mathrm{nucl}}(s)A_\omega (𝐱,ts)^{}A_{3\omega }(𝐱,ts)𝑑s.`$
Here the THG coefficient is decoupled from the other electronic susceptibilities; this is outside the usual model for pure fused silica in which the electronic contribution is considered instantaneous and the nuclear contrubution takes the form arising from a single damped harmonic oscillator, $`\chi (t;t_1,t_2,t_3)=\chi ^{\mathrm{elec}}\delta (tt_1)\delta (t_1t_2)\delta (t_2t_3)+\chi ^{\mathrm{nucl}}(t_1t_2)\delta (tt_1)\delta (t_2t_3)`$, with $`\chi ^{\mathrm{elec}}=n_2(1f_{\mathrm{Raman}})`$, $`\chi ^{\mathrm{nucl}}(t)=n_2f_{\mathrm{Raman}}(\tau _1^2+\tau _2^2)\tau _1^1\tau _2^2\mathrm{exp}(t/\tau _2)\mathrm{sin}(t/\tau _1)`$, $`n_2=n(\omega )cn_2^I/(2\pi )`$, $`n_2^I=2.8\times 10^{20}`$ m2/W, $`f_{\mathrm{Raman}}=0.18`$, $`\tau _1=12.2`$ fs, and $`\tau _2=32.0`$ fs. If the electronic THG susceptibility differs from the other electronic susceptibilities, the electronic contribution to the susceptibility is not instantaneous compared to all scales. Direct experimental measurements of the THG susceptibilities are available Milam.1998 ; GublerBosshard.2000 ; THG.older\_measurements . The more recent measurement of the THG coefficient $`\chi ^{\mathrm{THG}}`$ GublerBosshard.2000  is smaller than the electronic contribution to self-phase modulation Milam.1998  $`\chi ^{\mathrm{elec}}`$ by a factor of almost four; older measurements THG.older\_measurements  give a THG coefficient smaller than the electronic part of the self-phase modulation by a factor of about $`1.5`$. In the absense of direct experimental measurements of all the nonlinear polarization coefficients for doped silica, we use the simplest case $`\chi ^{\mathrm{THG}}=\chi ^{\mathrm{elec}}`$ in numerical simulations, but discuss how the results scale for different values of the THG coefficient. Vector effects are neglected. Inter-band Raman scattering is negligible at the carrier frequencies because of the fast relative phase oscillation; but in the pulse simulations, inter-band Raman scattering is possible between the lower frequencies within the higher frequency band, and the higher frequencies within the lower frequency band.

## II Numerical Simulations

A numerical NLS propagation scheme may be said to be accurate to all orders of dispersion and diffraction if it is accurate to as many orders as there are grid points. Accuracy up to all available orders requires the index of refraction over the entire numerically represented frequency range. In the split-step fast Fourier transform scheme, linear propagation is carried out in momentum space; the algorithm may be made accurate to all orders by putting the frequency dependent index of refraction directly into the formulas for the propagators TasgalBand.nls\_derivation , i.e.,

$`i_zA_\omega (𝐱,t)`$ $`=`$ $`^1\left\{\left[\sqrt{\beta (\omega _0+\omega )^2k_{}^2}\beta (\omega _0)\right]\{A_\omega (𝐱,t)\}\right\}+^1\left\{{\displaystyle \frac{2\pi (\omega _0+\omega )^2/c^2}{\sqrt{\beta (\omega _0+\omega )^2k_{}^2}}}\{P_\omega (𝐱,t)\}\right\},`$ (4a)
$`i_zA_{3\omega }(𝐱,t)`$ $`=`$ $`^1\left\{\left[\sqrt{\beta (3\omega _0+\omega )^2k_{}^2}\beta (3\omega _0)\right]\{A_{3\omega }(𝐱,t)\}\right\}+^1\left\{{\displaystyle \frac{2\pi (3\omega _0+\omega )^2/c^2}{\sqrt{\beta (3\omega _0+\omega )^2k_{}^2}}}\{P_{3\omega }(𝐱,t)\}\right\}.`$ (4b)
Here $``$ and $`^1`$ are Fourier and inverse Fourier transforms in $`x,y,t`$, $`k_{}=\sqrt{k_x^2+k_y^2}`$ is the transverse momentum, and $`\omega `$ is frequency. The linear dispersion effects are contained in $`\beta (\omega )=n(\omega )\omega /c`$; diffraction is contained in $`k_{}`$. The offsets by $`\omega _0`$ and $`3\omega _0`$ are due to the fact that the fields are slowly-varying envelopes about carrier waves at those frequencies; when these appear, $`\omega `$ is the frequency relative to the offset. There is no need to compute the coefficients for dispersion, diffraction, self-steepening, etc., explicitly because they are contained implicitly in the linear dispersion $`\beta (\omega )`$. Since we simulate focusing and collapse, during which diffractive, dispersive, and nonlinear length scales can easily change by factors of a hundred or more Trippenbach.1997 ; Agrawal.1995 , we allow the propagation step to vary such that it remains an order of magnitude less that the smallest relevent scale.

The effects that we are interested in are larger and clearer when the fundamental and TH are not too far from being phase-matched, $`\mathrm{\Delta }k3\beta (\omega _0)\beta (3\omega _0)0`$, or $`n(\omega _0)n(3\omega _0)`$. Optical materials may be doped to obtain desired properties Agrawal.1995 ; Nicacio.1993 . We consider silica doped with neodynium to obtain approximate phase-matching between the fundamental $`\lambda =1.5\mu `$m and TH $`\lambda =0.5\mu `$m. To model the frequency dependent index of refraction, we use the Sellmeier relation for fused silica Malitson.1965 , but add one additional resonance at $`\lambda _{\mathrm{Nd}}=0.59\mu `$m, which is neodynium’s largest resonance in the vicinity of our TH. For this, phase-matching is achieved at the Sellmeier coefficient $`B_{\mathrm{Nd}}0.0138`$, which corresponds to a few percent doping of the material. Phase-mismatch is varied by changing the dopant concentration (in this model, the Sellmeier coefficient $`B_{\mathrm{Nd}}`$). We take the conventional nonlinear coefficients for pure silica Agrawal.1995 , as given above.

We simulated the propagation of both pulses and continuous beams, over a range of light intensities, phase-mismatch values, and focusing strengths. The effects of the nonlinearity are clearly visible in Fig. 1, which shows the evolution of beams with varying intensity. For all simulations, we took initial conditions with zero power in the TH, and the fundamental in the form of a Gaussian with some radial phase factor $`A_1(x,t;z=0)=A_1\mathrm{exp}[(ic_1+2\mathrm{ln}2/W_1^2)(x^2+y^2)]`$. The radial phase factor $`c_1`$ can be related to the wavefront radius of curvature $`R_1`$ via the magnitude of the wavevector $`\beta (\omega )=n(\omega )\omega /c`$ according to the formula $`R_1=0.5\beta (\omega )/c_1`$. Because the fundamental focuses very intensely, while only a relatively small part of it is converted to TH, we illustrate the peak intensity of the fundamental beam and the power of the TH beam. In Fig. 1, the input beam has initial width FWHM = $`50\mu `$m, a radial phase factor $`c_1=5.0\times 10^9/\mathrm{m}^2`$ that brings the beam to a focus in 0.6 mm, and the material is phase-matched (the reciprocal of the phase-mismatch is more than an order of magnitude greater than the simulation distance). The fact that the normalized TH power varies with fundamental input intensity shows that the beams experience nonlinear dynamics, changing shape as the intensity changes. In Fig. 1, the initial wave front radii of curvature are small and their effect on the position of the focus overshadow the intensity dependence. One can see that the beams approach the linear limit at low intensities, as these curves start to almost overlap. Fig. 2 shows the normalized far-field TH power as a function of input power, for the data in Fig. 1 and for another set of runs with weaker focusing.

For a qualitative explanation of these figures, we first contrast our results with the model in which THG is the only nonlinearity Ward.1969 ; Boyd.1992 . Here, in the case with phase-matching, the TH power $`P_{\mathrm{TH}}(z)=(3/2)(2\pi /nc)^2(\omega /c)^4(\pi W_1^2)^3\chi _{}^{(3)}{}_{}{}^{2}I_\omega ^3/[1+(zz_{\mathrm{focus}})^2/z_1^2]`$ reaches a maximum at the focus of the fundamental and then drops off as the inverse square of $`z/z_1`$, where $`z_1=\pi W_1^2/\lambda _0`$ is the Rayleigh range of the fundamental and the TH beam, and $`W_1`$ is the width of the fundamental at its focus. Our simulations show that at high intensities, the Kerr effect causes the beams to lose reflection symmetry about the focal plane. With this asymmetry, THG from the incoming and outgoing beams does not fully interfere destructively, and allows some TH to propagate to the far-field. The lower-intensity curves in Fig. 1 do not quite drop off as the inverse square of distance because the numerical simulations did not start out from minus infinity, but began with merely a large (finite) beam width. The reflection asymmetry in Fig. 1 is due partly to the nonlinearities and partly to starting with a finite initial beam width. As long as the amount of energy in the TH is relatively small, the nonlinear dynamic effects will remain in the fundamental beam, and the TH peak intensities and beam powers may scale up or down by a uniform factor, but will be otherwise unaffected. Clearly, there is significant variation with the fundamental input power. Moreover, this variation is quite different for different focusing conditions.

In another series of simulations, phase-mismatch is varied. Fig. 3 shows the peak intensity of the fundamental and the TH beam power with initial FWHM = $`50\mu `$m, radial phase-factor $`c_1=5.0\times 10^9/\mathrm{m}^2`$, and phase-mismatch from -12.9 to 51.8 mm-1. Because the TH power first increases and then decreases with phase-mismatch, we show the former range on one plot and the latter on another. Clearly, there is a residual far-field TH and its power depends on phase-mismatch. Fig. 4 shows the far-field TH conversion efficiency as a function of phase mismatch, for the runs in Fig. 3 and another set of runs with weaker focusing, $`c_1=2.5\times 10^9/\mathrm{m}^2`$. As in the linear model Ward.1969 , THG by a focusing beam is maximized around a certain phase-mismatch; but the nonlinear dynamics complicate the results considerably.

In a further series of simulations, the intensity and phase-mismatch were held constant, and the radial phase factor varied. This mixes up a few physical effects, since initial conditions with small radial phase factors cannot be considered as starting from the far-field. Fig. 5 has extremely small phase-mismatch ($`\mathrm{\Delta }k=0.067`$ mm-1), and radial phase factors which vary from zero (i.e., starting at a focus) up to $`c_1=5.0\times 10^9/\mathrm{m}^2`$. The horizontal axis is on a log scale to help visually distinguish the superimposed sumulations. The peak fundamental intensity varies, but over the range studied, the peak TH power is rather insensitive to the wavefront radius of curvature of the fundamental, but the far-field TH power does vary considerably. Fig. 6 shows the far-field TH power as a function of radial phase factor, for the runs in Fig. 5, and also for a series of runs with significant positive and negative phase-mismatch. In increasing (decreasing) the tightness of the focusing, the nonlinear and diffraction lengths both decrease (increase); the complexity of the results for far-field THG reflect the complexity of the nonlinear dynamics of the beam.

Fig. 7 shows the dependence of the third-harmonic power on the when the THG suseptibility $`\chi ^{\mathrm{THG}}`$ is varied. The curve labeled $`\chi ^{\mathrm{THG}}=\chi ^{\mathrm{elec}}`$ corresponds to the same conditions used in Fig. 3 but with $`\mathrm{\Delta }k=0.067`$ mm-1 (i.e., almost phase-matched). The third harmonic power scales with $`|\chi ^{\mathrm{THG}}|^2`$ when the third-harmonic intensity is small since then the nonlinear dynamics of the fundamental is unaffected by TH and the generation of TH field is proportional to $`\chi ^{\mathrm{THG}}`$. The values of $`\chi ^{\mathrm{THG}}`$ used in the calculations shown in Fig. 7 are a factor of $`1.0`$, $`1/3.9`$, and $`1/1.5`$ times the value used in the previous figures, where the latter two factors correspond to the measured values of $`\chi ^{\mathrm{THG}}`$ reported in Refs. Milam.1998 ; GublerBosshard.2000 ; THG.older\_measurements

Numerical simulations showed the dynamics of pulses to differ from those of continuous waves in essentially two ways. First, the group velocities of the fundamental and TH will, except for special cases, not be the same. A TH pulse will thus generally walk off from a fundamental pulse. This tends to reduce destructive interference; it also limits the effective distance over which the pulses interact. Compared to beams, TH pulses tend to carry off a larger part of the fundamental pulse energy. Secondly, the dynamics of pulses (3+1–dimensional) are qualitatively different to those of continuous beams (2+1–dimensional). For strongly focused but not very short pulses, these differences tend to be minor. A thorough analysis of focused pulses in this system is quite involved, and will be pursued elsewhere. Thus, the detailed results for continuous-wave beams apply to pulses for which group-velocity differences are relatively small or for pulses that are relatively long.

## III Conclusions

In a nonlinear medium, intense beams or pulses of finite diameter which converge to a focus and then diverge may exhibit nonlinear dynamics that significanly affect propagation dynamics. These nonlinear effects break the beam’s reflection symmetry about the focal plane. The greater the intensity, the bigger the difference between incoming and outgoing beams. THG with such an input beam or pulse produces TH in the far-field when the phase-mismatch between the fundamental and its third-harmonic is zero or negative (and small). We have quantitatively demonstrated this for several cases. When phase-mismatch is positive, where some far-field TH power is possible in the essentially linear case, the nonlinear beam dynamics complicate the accumulation of TH power. The non-zero far-field THG for zero or negative phase-mismatch – a qualitatively new effect for homogeneous media – is affected by, and thus contains information about the medium in the region of the beam focus. For an inhomogeneous medium, and for THG microscopy, these effects should be understood, either to be utilized or better avoided.

###### Acknowledgements.

Richard Tasgal gratefully acknowledges a Kreitman Foundation Fellowship. This work was supported in part by a grant from the Israel Science Foundation for a Center of Excellence (grant No. 8006/03) and by KBN as a reasearch grant 2003-2006 (2P03B04325).

## List of Figure Captions

Fig. 1. Peak intensity of the fundamental ($`\lambda _1=1.5\mu `$m) beam, and power of the third-harmonic ($`\lambda _3=0.5\mu `$m) normalized by the cube of the input power. The initial conditions have a range of intensities (1.01, 4.05, 16.2, 36.4, 64.7, 101, 122 and 145 KW, corresponding to circled points in the appropriate curve in Fig. 2), but identical beam width, FWHM = $`50\mu `$m, radial phase factor, $`c_1=5.0\times 10^9/\mathrm{m}^2`$ (i.e., transverse phase $`\mathrm{exp}[ic_1(x^2+y^2)]`$), and all are phase-matched.

Fig. 2. far-field third-harmonic beam power, normalized by the cube of the input power, vs. input power. The medium has TH phase-matching. Two curves are shown, one the far-field results from Fig. 1, and another with smaller radial phase factor, that brings the beam to a focus in about 1.2 mm.

Fig. 3 Peak intensity of the fundamental beam, and power of the third-harmonic. The curves represent beams with identical power, width, and radial phase factor, while the phase-mismatch goes from -12.9 to 51.8 mm-1 (corresponding to circled points in the appropriate curve in Fig. 4). The middle plot shows the TH increasing with phase-mismatch, and the bottom plot shows TH decreasing with phase-mismatch.

Fig. 4. Far-field conversion efficiency vs. phase-mismatch. Two curves are shown, one the far-field results from Fig. 3, and another with smaller radial phase factor, that brings the beam to a focus in about 1.2 mm. Initial intensities are held constant.

Fig. 5. Peak intensity of the fundamental beam, and power of the third-harmonic versus propagation distance for different fundamental radial curvatures. The curves show beams with phase-matching and initially identical power and width; the initial radial phase factors go from zero (starting at a focus) up to $`5.0\times 10^9/\mathrm{m}^2`$, by increments of $`0.5\times 10^9/\mathrm{m}^2`$ (the curves here correspond to circled points in the appropriate curve in Fig. 6).

Fig. 6. Far-field conversion efficiency vs. radial phase factor. Two curves are shown, one the far-field results from Fig. 5 (phase-matching), and another in a medium with positive phase-mismatching.

Fig. 7. Peak fundamental intensity and third-harmonic power versus position $`z`$, for various values of $`\chi ^{\mathrm{THG}}`$. The third harmonic power scales with $`|\chi ^{\mathrm{THG}}|^2`$ when the third-harmonic intensity is small.


text:
# 1 Temperature profile (left) and best simulated and observed SEDs (right, short dashed and dotted lines, respectively). Details regarding the summarized target list at several wavelengths from the literature are given in Rengel et al. (2004b).

Baltic Astronomy, vol. 12, XXX--XXX, 2003.

Modelling of the continuum emission from Class 0 sources

M. Rengel1, D. Froebrich2,1, S. Wolf3 and J. Eislöffel1

1   Thüringer Landessternwarte Tautenburg, Sternwarte 5, 07778 Tautenburg, Germany

2   Dublin Institute for Advanced Studies, 5 Merrion Square, Dublin 2, Ireland

3   California Institute of Technology, 1201 East California Blvd., 105-24, Pasadena, CA91125, USA

Received October 15, 2003

## Abstract

Class 0 sources are objects representing the earliest phase of the protostellar evolution. Since they are highly obscured by an extended dusty envelope, these objects emit mainly in the far-infrared to millimetre wavelength range. The analysis of their spectral energy distributions with wide wavelength coverage allows to determine the bolometric temperature and luminosity. However, a more detailed physical interpretation of the internal physical structure of these objects requires radiative transfer modelling. We present modelling results of spectral energy distributions of a sample of nine Class 0 sources in the Perseus and Orion molecular clouds. The SEDs have been simulated using a radiative transfer code based on the Monte Carlo method. We find that a spherically symmetric model for the youngest Class 0 sources allows to reproduce the observed SEDs reasonably well. From our modelling we derive physical parameters of our sources, such as their mass, density distribution, size, etc. We find a density structure of $`\rho r^2`$ for the collapsing cores at young ages, evolving to $`\rho r^{3/2}`$ at later times.

Key words:  Radiative transfer–Stars: formation–dust

1.   Introduction

A Class 0 source is an object representing the earliest phase of the protostellar evolution (e.g. André et al. 1993). It consists of a central protostellar object surrounded by an infalling dusty envelope and a flattened accretion disk. Many (presumably all) Class 0 sources are associated with a bipolar molecular outflow. Because of their difficult detection (their photospheres are highly obscured and the objects spend only a short time ($``$ $`10^4`$ yr) in this evolutionary phase), and the lack of knowledge in model protostellar parameters (distribution, composition and size of dust, density and temperature distribution, etc.), constraining physical properties to Class 0 objects is quite difficult. Continuum sub-mm observations of Class 0 sources, however, let us detect dust emission of the massive circumstellar envelope of these sources. Information from the dust emission, a physical model of the extended envelope as a pre-requisite, and the implementation of techniques like the blackbody and the envelope fitting procedures, open a way to interpret the structure of Class 0 sources.

2.   An envelope model

How can we describe the protostellar emission from the envelope? We adopt the standard envelope model (Adams 1991) to keep the problem as simple as possible. The circular symmetry of the observed emission (Rengel et al. 2004a), and the lack of significant internal structure justifies the simplicity of a spherical model case. If the emission is optically thin, the observed intensity for a spherical symmetric protostellar envelope at an impact parameter b is given by Eq. 1 assuming single power-laws of the opacity, the radial temperature and density distribution ($`\beta `$, q and p, respectively).

$$I_\nu \left(b\right)=2\kappa _\nu _b^{r_0}B_\nu \left[T_d\left(r\right)\right]\rho \left(r\right)\frac{r}{\sqrt{r^2b^2}}𝑑r$$
(1)
$`r_0`$ is the outer radius, $`\rho `$ is the density, $`T_d`$ the dust temperature, $`\kappa _\nu `$ the opacity of the dust grains and $`B_\nu \left[T_d\left(r\right)\right]`$ the Planck function at dust temperature Td. If the emission is in the Rayleigh-Jeans limit and if $`r_0`$$``$ b, Eq. 1 can be approximated to $`I_\nu \left(b\right)/I_\nu \left(0\right)=\left(b/b_0\right)^m`$ (where m is the power-law index of the observed intensity profile, and b0 is the normalization factor to the peak emission).

Results from the standard theory provide the first direct insights into observable estimations. Nevertheless, the temperature profile will diverge from a single power-law q as the envelope becomes optically thick at the primary wavelengths of energy transport (inner portion of the envelope) (e.g. Shirley, Evans & Rawlings 2002). Here it becomes necessary to calculate the temperature distribution self-consistently by implementation of a radiative transfer code. To calculate the observables (temperature and density distributions, and intensity map) we implemented a 1D spherical symmetric radiative transfer code (MC3D, Wolf et al. 1999). The model consists of a central heating source, with stellar luminosity $`L_{}`$ and an effective temperature $`T_{}`$, embedded within a spherical symmetric dust envelope of mass $`M_{env}`$. The envelope is specified by an outer radius $`R_{out}`$, a density power-law index p, and inner radius $`r_{in}`$, which is set by the dust destruction radius $`r_{sub}`$, adopted to a temperature of 1500 K.

3.   Modeling continuum submillimetre emission

First constrains to physical properties (spectral indices, masses, radial profiles and sizes) of 15 objects are derived from continuum SCUBA imaging of six star formation regions in the Perseus and Orion molecular cloud complexes (Rengel et al. 2004a). Accurate SEDs are computed for nine sources combining existing multi-wavelength surveys with new sub-mm data (Rengel et al. 2004b), and fitting a modified blackbody curve to the fluxes (Froebrich et al. 2003). We derived the bolometric temperature and luminosity, size of the envelope and sub-mm slope for each source. Furthermore, dust temperature distribution, SED, and an intensity map are derived self-consistently by the MC3D code. Physical parameters of the sources (e.g. envelope masses, density distributions, sizes and sublimation radius) were derived by finding the consistency between observed and modelled SEDs.

We created a bolometric luminosity-temperature diagram for the nine Class 0 sources and compare the positions of the objects with an evolutionary model in order to derive the age and mass of the envelope of each source.

4.   Results

In Fig. 1, the calculated temperature profile (left), and the observed and best simulated SEDs for L1448 C are shown as an example. Table 1 lists the characteristics of the best-fitting model parameters. The envelope emission is calculated as an image from the MC3D code, and convolved by the SCUBA beam. The simulated intensity profile is calculated at 850 $`\mu `$m and compared with the observed profile to test the symmetric model.

Is there a correlation of the power-law index p with time? We plot the estimated values of p as a function of t. Fig. 2 suggests a density structure of $`\rho r^2`$ at younger ages, evolving to $`\rho r^{3/2}`$ at later times. In order to investigate this with a larger spread, a Class 1 sample is mandatory.

5.   Conclusions

Detailed physical interpretation of the internal structure of nine objects was carried out with intensive computer modeling. A simple spherically symmetric model envelope, and assumptions about density and dust distributions following the standard envelope model reproduce reasonably well the observed SED and the radial profiles of the sources. $`Tr^{0.4}`$ is a good approximation for the sample. The radial temperature distribution as function of distance, however, departs significantly from the optically-thin assumption, an observational derivation of a single-power law q of 0.4 for radii $`<`$ 10 AU. Modelling results indicate a density profile well described by a power-law between p=1.5 to 2, which is expected by all of the collapse models and numerical studies. A density structure of $`\rho r^2`$ at younger ages, evolving to $`\rho r^{3/2}`$ at later times was found.

SED and radial profile fits can constrain physical parameters of the sources such as envelope masses, density distributions, sizes and sublimation radii. Nevertheless, observations at 10-300 $`\mu `$m, and the inclusion of outflow and disk, and other geometries will decrease the differences.

| Table 1: Best envelope fit results for the embedded sources. The temperature of the central object T is set to 3500 K. $`\chi ^2`$ quantifies the agreement between model and data. The age is estimated according to the model of Smith (2000) and given in 103 yrs. | | | | | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Object | Tbol \[K\] | L \[L\] | Rsub \[AU\] | p | Rout \[AU\] | age | Menv \[M\] | $`\chi ^2`$ |
| L 1448 NW | 27 | 6 | 5 | 2.0 | 4000 | 13.8 | 2.8 | 6.4 |
| L 1448 C | 40 | 11 | 5 | 1.6 | 3000 | 19.5 | 1.7 | 2.6 |
| RNO 15 FIR | 45 | 8 | 3 | 1.6 | 3500 | 20.7 | 1.0 | 1.8 |
| NGC 1333 iras 1 | 52 | 13 | 3 | 1.5 | 4500 | 21.9 | 1.5 | 8.6 |
| NGC 1333 iras 2 | 48 | 68 | 3 | 1.5 | 7000 | 26.1 | 3.0 | 0.1 |
| HH 211-mm | 30 | 5 | 3 | 1.9 | 6000 | 13.8 | 4.1 | 4.5 |
| L 1634 | 42 | 19 | 3 | 1.6 | 6000 | 21.9 | 2.8 | 1.1 |
| L 1641 N | 45 | 85 | 3 | 1.5 | 8500 | 25.2 | 7.0 | 1.0 |
| L 1641 sms III | 50 | 80 | 3 | 1.5 | 10000 | 27.6 | 6.3 | 0.6 |

REFERENCES André, P., Ward-Thompson, D., Barsony, M. 1993, ApJ, 406, 122
Adams, F. 1991, ApJ, 382, 544
Froebrich, D. et al. 2003, MNRAS, in press
Rengel, M., Eislöffel, J., Hodapp, K.W. 2004a, in prep
Rengel, M., Froebrich, D., Wolf, S., Eislöffel, J. 2004b, in prep
Shirley, Y., Evans, N.J., Rawlings, J. 2002, ApJ, 575, 337
Smith, M.D. 2000, IAJ, 27, 1
Wolf, S., Henning, T., Stecklum, B. 1999, A&A, 349, 839