Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.astro.spbu.ru/staff/ilin2/INTAS2/PAP/P1_4.pdf Äàòà èçìåíåíèÿ: Fri Nov 19 16:18:49 2010 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 02:58:54 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: rho ophiuchi

Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921­ 937

www.elsevier.com/locate/jqsrt

Computation of light scattering in young stellar objects
P.W. Lucas
Department of Physical Sciences, University of Hertfordshire, College Lane, HatÚeld AL10 9AB, UK Received 1 June 2002; accepted 13 August 2002

Abstract A Monte Carlo light scattering code incorporating aligned non-spherical particles is described. The major e ects on the ux distribution, linear polarisation and circular polarisation (CP) are presented, with emphasis on the application to young stellar objects (YSOs). The need for models with non-spherical particles in order to successfully model polarisation data is reviewed. The ability of this type of model to map magnetic Úeld structure in embedded YSOs is described. The possible application to the question of the origin of biomolecular homochirality via UV CP in star forming regions is also brie y discussed. ? 2003 Elsevier Science Ltd. All rights reserved.
Keywords: Monte Carlo

1. Introduction The technique of Monte Carlo modelling of light scattering in spatially resolved nebulae has proven very useful in the Úeld of astrophysics. This stochastic method works by generating millions of photons and allowing them to propagate in random directions, undergoing single or multiple scattering into directions sampled from a theoretical or empirical probability distribution function. Those photons which eventually escape the bounds of the system are mapped on to a grid using their location and direction after the last scattering. This grid then represents a 2D image, such as that obtained with a telescope using an optical CCD camera or an infrared array camera. Monte Carlo modelling is used to investigate the physical structure and composition of both the dusty nebulae surrounding young stellar objects (YSOs), which are stars in the process of formation, and protoplanetary nebulae and planetary nebulae (PPNe and PNe), which surround stars near the end of their life cycle that have ejected their outer atmospheres. Another application is imaging and spectroscopic modelling of electron scattering for hot stars with an ionised stellar wind, e.g. [1,2] but in this paper only imaging and polarimetric imaging models of scattering by dust are considered.


Corresponding author. Fax: +44-1707-284644. E-mail address: pwl@star.herts.ac.uk (P.W. Lucas).

0022-4073/03/$ - see front matter ? 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0022-4073(02)00329-1


922

P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

Until recently, published models considered only spherical dust grains, using the Mie theory to generate the appropriate amplitude and Stokes matrices for scattering in di erent directions. Such models have been successful in explaining most of the major observed features of YSOs, e.g. [3­8]. However, several features of the linear and circular polarisation (CP) data have not been successfully reproduced by models using spherical grains (see Section 2), so a number of researchers have begun to incorporate scattering by aligned non-spherical particles (hereafter dichroic scattering) into their codes. Wolf, Voshchinnikov and Henning [9] published the Úrst demonstration of such a model, incorporating dichroic scattering but not fully treating the e ect of dichroic extinction on the polarisation state. Such a model is strictly applicable only to the optical thin case (e.g. a planetary nebula). Whitney and Wol [10] have made a detailed investigation of the effect of dichroic scattering and extinction on images and linear and CP maps, applied to spherical clouds and to YSOs for the case that the particles are uniformly aligned with a magnetic Úeld along the axis of the system. The code presented here is conceptually similar to [10], so this paper deal mainly with results and details of the algorithm which were not already covered by those authors. 2. Limitations of spherical grain treatment and applications of non-spherical grains Models with spherical grains can successfully reproduce most features of imaging data but are less successful with polarimetric data. SpeciÚc examples where previous Monte Carlo work either fails to match data or leads to a misleading match to data are: (1) The phase function (distribution of scattered light) in imaging data. Spheres and even spheroids tend to overestimate the amount of back-scattering since their symmetry leads to constructive interference at de ection angles near 180 which are unlikely to occur for realistic grains with imperfect symmetry. A consequence is that the grain size distribution inferred by Útting an observed phase function is likely to include too many large grains (since large grains back-scatter less). This error may have occurred in recent papers on the dust ring around the GG Tau A binary YSO system [11,12]. (2) The CP produced by spherical grains is very low and arises only from multiple scattering if the source of radiation is essentially unpolarised, as is the case for most stars. Observations of high-infrared CP in YSOs (up to 20%, e.g. [13­16]) require the presence of signiÚcantly non-spherical grains, which can produce high polarisation via both single and multiple scattering or by dichroic extinction of linearly polarised light. (3) Near-infrared linear polarisation (LP). The extended regions of aligned vectors seen in a minority of YSOs, e.g. L1551 IRS5, IRAS 04361 + 2547, IRAS 04239 + 2436 and YLW16a [8,17] have not been reproduced with Monte Carlo models using spherical grains. Such models have more success (at least qualitatively) in reproducing the thin region of aligned vectors or `polarisation disk' seen more frequently in YSOs by a combination of multiple scattering and instrumental blurring of the polarisation pattern due to limited spatial resolution. In addition, long standing direct evidence for emission from non-spherical grains at a range of temperatures and distances from YSOs is provided by the many published observations of LP at wavelengths 5 ¡ ¡ 1500 m (e.g. [18,19]). At these wavelengths, scattering becomes unimportant due to the small size parameter (and hence low albedo) of grains which predominate in the disks and


P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

923

envelopes around YSOs. This has been modelled in detail by Aitken et al. [20]. The size parameter is deÚned as x =2 a= . The motivation for modelling with non-spherical grains stems from the important information that can be extracted from polarisation data. In the astrophysical context it is routine to infer the direction of the tangential component of the magnetic Úeld in the interstellar medium or near a YSO from LP data. The magnetic Úeld is thought to play a major role in the star formation process but this is poorly understood. Observations of magnetic Úeld structure in the central regions of the circumstellar envelope of a YSO could potentially provide valuable tests of di erent theoretical models of star formation via the collapse of a large molecular cloud. While the mechanism for dust grain alignment is often unclear, e.g. [21­24], it is generally accepted that the alignment is such that the rotation axis giving the largest moment of inertia is on average parallel to the magnetic Úeld. Linear polarimetry in mid-infrared and sub-mm wavebands permits a direct measurement of the direction of tangential component of the magnetic Úeld along the line of sight to a star but it is generally not possible to map the magnetic Úeld in circumstellar matter, due to poor spatial resolution in the sub-mm and the small size of the emitting region in the mid-infrared (see [25] in this issue). In the near-infrared, YSOs often display bright extended nebulae which can be observed at high signal to noise. Monte Carlo modelling of the linear and CP produced by dichroic scattering and extinction in this waveband is at present the only way that the magnetic Úeld structure can be studied on scales of tens to hundreds of astronomical units (AU), so this should be pursued despite the complexity of the problem. Another important application of this type of Monte Carlo model is to the possible connection between high CP in YSOs and the origin of biomolecular homochirality [13]. The homochirality of DNA and amino acids in terrestrial organisms may well have an extraterrestrial cause (this is supported by data from meteorites, e.g. [26]) and it is thought to be a prerequisite for the origin of life. Circularly polarised light in star forming regions is one possible explanation for an extraterrestrial origin of homochirality. The relevant waveband for this hypothesis is the near-ultraviolet: photons at these wavelengths can preferentially dissociate organic molecules with left-handed or right-handed (L or D) chiral structure, depending on the CP state. It is very di cult to observe YSOs in this waveband, owing to the very high optical depth of the ambient medium, so Monte Carlo simulations are the only way to investigate the degree of CP at the relevant wavelengths. 3. The code and dust grain model The algorithm incorporates several measures to improve e ciency which have been developed by others who work with stochastic models, see e.g. [27,5] and references therein. In this work, we model a homogeneous sphere and an axisymmetric 2D disk+envelope distribution of matter which is a simple model of a YSO (see Fig. 1). Non-axisymmetric 3D systems can be modelled with only a minor modiÚcation of the code (see [8]), but such models are far more computationally intensive since output photons must then be binned in the azimuthal ordinate as well as polar viewing angle. A unique 3D solution is hard to identify with 2D imaging and polarisation data, though observations at several wavelengths provide some depth perception by penetrating to a range of optical depths. The 2 â 2 amplitude matrix used in the phase function and the 4 â 4 Stokes matrix which provides the new Stokes vector of each photon after scattering are precalculated using the T-matrix codes of Mishchenko (see http://www.giss.nasa.gov/crmim). The complexity of the code is increased com-


924

P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

Fig. 1. Schematic diagram of a young stellar object. The magnetic Úeld is vertical and parallel to the rotation axis of the system. The density of the envelope declines monotonically with distance from the protostar and distance from the disk plane. The density spans several orders of magnitude, as indicated by the wrapped greyscale shading. The cavity is taken to be evacuated.

pared to the treatment for spheres because all 4 elements of the amplitude matrix and all 16 elements of the Stokes matrix are non-zero. At present the Monte Carlo code works only with oblate spheroids, whose appearance does not change with rotation about the axis of greatest inertia. However, it can treat an arbitrary magnetic Úeld structure. The quadruple precision version of Mishchenko's code `amplq.new.f' is required for ultraviolet work with grains which are highly attened (at least in the astronomical context) with axis ratio of 3:1 and have equivalent surface area sphere radii up to 0:75 m, which may be suitable for dust in the envelopes around YSOs. For longer wavelengths the double precision version `ampld.new.f' is used. The Stokes matrix elements are then used to calculate the 3 non-zero elements of the 4 â 4 extinction matrix which modiÚes the polarisation state of each photon during ight between scatterings (see [28,10,29]). The phase function, Stokes matrix elements and extinction matrix elements are stored as look-up tables in the code, usually sampled at intervals of 2 (essential for ultraviolet work at large size parameters) or sometimes 4 for work at small size parameters (x¡ 1), appropriate in the infrared. The Mishchenko code is used to produce matrix elements averaged over a distribution of grain sizes for the full range of polar scattering angles (de ection angle, D) and the full range of azimuthal and polar grain orientation angles ( , Ú). The algorithm proceeds as follows: (i) Millions of `unpolarised photons' with Stokes vector (IQUV ) = (1; 0; 0; 0) are generated at random points on the stellar surface or a hot accretion disk and are assigned random directions of ight. Real polarised photons could be used but unpolarised photons are more e cient, representing the average of many polarised photons, like a classical light wave.


P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

925

(ii) Each photon steps along its path, gathering optical depth until it is absorbed or scattered after reaching a randomly generated optical depth ( = -ln(1 - R), for random number 0 ¡R¡ 1) or else escapes the system. To improve e ciency, photons are not allowed the albedo-dependent chance to be absorbed and discarded until several scatterings have occurred. Instead, the `weight' or importance of the photons is reduced by multiplying the Stokes vector by the albedo, until Stokes I ¡ 10-3 , at which point the weight is held constant and absorption occurs in a random manner. (iii) Treatment of the Stokes vector upon scattering. The azimuthal scattering angle of the photon, , is set at zero for purposes of calculating the new Stokes vector after scattering, since and are not independent, i.e. the azimuthal grain orientation is calculated relative to the azimuthal de ection angle. This provides an essential saving in the size of the Stokes matrix look-up tables, which is permissible since the Stokes vector is rotated into the scattering plane for the calculation, and then rotated back into the systemic frame, as described by Hovenier and van der Mee [30], and in [31,32]: I I Q Q = L( - 2 )ZL(- 1 ) ; (1) U U V
scat

V matrix Z Z Z Z
12 22 32 42

inc

where the Stokes Z11 Z21 Z= Z31 Z
41

Z Z Z Z

13 23 33 43

Z

Z24 Z34 Z
44

14



and L is the simple rotation matrix described in those papers, acting on the rotation angles 1 and 2 . These angles and their signs are as deÚned by [30]. The deÚnitions of LP position angle and the sense of CP are also as deÚned in that work. The position angle is zero along the projection of the symmetry axis of the system on to the plane perpendicular to the Poynting vector. Position angle increases in a clockwise rotation from the point of view of the distant observer, so this is a left-handed deÚnition if the thumb points along the Poynting vector. If the zero angle is regarded as due north this is an increase to the west (in the astronomical, upward looking sense of west, opposite to the terrestrial sense). This is opposite to the usual astronomical observational convention of polarisation increasing from north through east but we retain it for consistency with previous modelling work. Positive Stokes V corresponds to right-handed CP, with the thumb indicating the Poynting vector. (iv) Determination of scattering direction. Upon scattering, the polar and azimuthal de ection angles of the photon (D, ) are randomly sampled from the full scattering phase function by rejection sampling. is a left-handed rotation about the direction of the Poynting vector, with the zero of parallel to the position angle of LP. The new direction of ight in the systemic frame is then calculated by transforming to the frame deÚned by the incident photon. The phase function is dependent on the grain orientation angles ( , Ú) and the linear and CP state of the incident photon (see Appendix A).


926

P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

(v) During ight, the polarisation state of each photon is modiÚed at each step by the dichroic extinction of the medium. To do this the Stokes vector is rotated into the plane deÚned by the magnetic Úeld (with which the grains are aligned) and the direction of ight (Poynting vector), modiÚed as shown below and then rotated back into the systemic frame. Thus dichroic extinction by an arbitrary magnetic Úeld structure can be treated using only 3 independent non-zero elements of the 4 â 4 extinction matrix, after the manner described by [10], see Eqs. (2)­(5). For computational purposes this is more e cient than calculating the full extinction matrix in the manner expressed by [27], which incorporates the rotations. I
k +1

=0:5 â{(Ik + Qk ) exp[ - n(K11 + K12 ) s]+(Ik - Qk ) exp[ - n(K11 - K12 ) s]}; =0:5 â{(Ik + Qk ) exp[ - n(K11 + K12 ) s] - (Ik - Qk ) exp[ - n(K11 - K12 ) s]}; = exp[ - nK = exp[ - nK
11

(2) (3) (4) (5)

Q U V

k +1 k +1

s]{Uk cos(nK s]{Vk cos(nK

34

s) - Vk sin(nK s)+ Uk sin(nK

34

s)}; s)};

k +1

11

34

34

K11 , K12 and K34 are the 3 independent elements of the extinction matrix, evaluated at grain orientation ( ; Ú)=(0;Ú). The k subscript denotes the k th step of the photon along the path, s is the physical length of the step and n is the number density of grains. (vi) Those photons which eventually escape the system are output into Úle bins corresponding to the full range of polar viewing angles (0 ­180 ). In 2D axisymmetric models, all photons are assumed to be output into the observer's azimuthal viewing angle. This means that the azimuthal coordinate of each photon is not speciÚed when the photon is generated and the code only calculates the change in azimuthal coordinate during ight. For systems which have mirror symmetry about the equatorial plane, the equivalence of the output into the northern and southern hemispheres is a useful test of the code. A separate program is then used to map the photons on to a 2D image grid, using the location and direction of the photon at last scattering. For each grid element, the 4 Stokes components are summed for all photons mapped on to that element. This generates the raw (I; Q; U; V ) model output maps. For comparison with real data the raw (I; Q; U; V ) maps are then convolved with a point spread function appropriate for the observations. Accurate measurement and selection of the point spread function is crucial for useful modelling in cases where the optical depth to the central star 6 6, since the strong central ux peak (which is often highly polarised by dichroic extinction) tends to modify or obscure features arising from scattered ux. This algorithm is very similar to that developed in [10] and both groups chose to employ the Mishchenko code to provide the input matrices. The codes of the two groups have been bench tested against each other for the case of the uniform axial magnetic Úeld, applied to both uniform spherical clouds and axisymmetric 2D disk+envelope distributions of matter suited to YSOs. The only signiÚcant known di erences are: (a) the algorithm in [10] can handle precessing oblate or prolate spheroids, whereas this code uses only perfectly aligned oblate spheroids at present; (b) this algorithm can handle arbitrary magnetic Úeld structures, whereas [10] dealt only with the case of an axial Úeld. Any other di erences are believed to be trivial, though the source codes have not been compared. The precision of the code is determined by the length of photon steps in stage (ii) and the sampling frequency of the look-up tables. It can be altered to suit the desired spatial resolution and angular binning of the output maps. The necessary precision for a given application


P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

927

is best determined empirically. We note that bench testing is highly desirable for computations with non-spherical particles, since the many of the symmetries in photon output which can be relied upon to test computation with spheres no longer apply, particularly if the magnetic Úeld is non-uniform. The algorithm was coded in Fortran 77, and can be run in series or in parallel. Nearly 100% parallelisation e ciency is achieved by evenly distributing the photons between CPUs. The parallel code is compiled with the Silicon Graphics MIPSpro Fortran 77 compiler but it should run in series with any compiler compatible with Fortran 77. The program speed depends mainly on the complexity of the structure of circumstellar matter: for complicated systems like a YSO, the move from spherical grains to non-spherical grains does not signiÚcantly increase the run time. The large look-up tables lead to far greater RAM requirements, however, and a sampling of the matrices at intervals of 1 would be beyond the capacity of most desktop PCs at present. All simulations shown in this paper are run for a wavelength of 2:2 m and a grain mixture consisting of (1) silicate particles (n =1:7 - 0:03i) with equivalent surface area radii in the range 0.005 ­0:35 m, and (2) small amorphous carbon particles [33] (n =2:77 - 0:80i) in the size range 0.005 ­0:030 m. All particles are oblate with a 3:1 axis ratio. For both grain types the size distribution is given by n(r ) r -3:5 . The maximum grain size of 0:35 m is towards the lower end of the plausible range in YSOs, comparing with the widely adopted maximum size of 0:25 m for interstellar dust. The choice of a mixture of grains with fairly small size parameters contributes to the exotic behaviour of CP in the results shown here (see Section 4.2). The scattering is dominated by the shiny silicates but the absorptive amorphous carbon particles serve to reduce the albedo to approximately 0.16. This mixture is used as a compromise between possible dielectric and highly absorptive mixtures. The true composition of circumstellar dust near YSOs is uncertain but is known to include water ice and CO ice (e.g. [34]), believed to be deposited in mantles around a silicate core. The ices are likely to have only a small imaginary component of refractive index, like the silicate used here, but the mantles may have metallic inclusions. In the models shown here, the opacity of the grain mixture is set to unity and absorbed into the product of opacity and the mass density of the medium. 4. Model results 4.1. Homogeneous sphere with a uniform magnetic Úeld The spatially resolved (I; Q; U; V ) map of a sphere illuminated by a point source at its centre has been discussed by [10] for the case of 2:1 prolate and oblate spheroids with varying degrees of alignment. Their results are reproduced by this code (considering only perfectly aligned oblate grains, which show the strongest e ects of non-spherical symmetry). Here some additional features arising from 3:1 grains are noted which arise in the case of fairly high optical depth. 4.1.1. Stokes I : a sphere can appear as a bipolar nebula In [10], it is described how an optically thin sphere appears slightly elongated along the axis of the magnetic Úeld, due to the greater optical depth and hence larger scattered ux along the axis (Fig. 2(a)). For the case of high optical depth the e ect of enhanced scattering by grains located along the axis is outweighted by the strong extinction in that direction, which prevents light from reaching the upper and lower parts of the sphere, as shown in Fig. 2(b). Thus the sphere appears


928

P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

Fig. 2. Images of a homogeneous sphere with vertical magnetic Úeld. (a) =0:5, showing slight vertical elongation. (b) = 10, showing horizontal elongation and a bipolar structure in the central region. Optical depth is less along horizontal paths than vertical paths due to 3:1 oblate shape of the aligned grains.

(a)

(b)

Fig. (b) are the

3. Linear polarisation vector maps for a homogeneous sphere. (a) =0:5, showing LP dominated by dichroic scattering. = 5, showing LP dominated by dichroic extinction. The magnetic Úeld is vertical and the long axes of the grains horizontal. 100% polarisation corresponds to a vector length of 0.8 pixels. The direct ux from the central source at (0,0) position is not included here.

fairly round near an optical depth of unity and the image becomes elongated perpendicular to the magnetic Úeld at high optical depth. Here it is shown that for 3:1 grains with an optical depth = 10 along the axis from the source to the edge of the sphere, the greater extinction along the axis squashes the ux distribution so much that it becomes bipolar. In this model, the optical depth perpendicular to the axis is =2:9 from source to edge of sphere. The bipolar structure begins to appear at =5. 4.1.2. Stokes Q and U For attened grains the LP pattern departs from the centrosymmetric pattern commonly observed in circumstellar nebulae. In this model (Fig. 3(a) and (b)), vectors oriented along the magnetic Úeld


P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

929

axis are deÚned as having a position angle of zero, i.e. positive Q and zero U . The scattered electric vector is enhanced in the plane containing the two long axes of the grains, leading to enhanced negative Q at low optical depth (i.e. vectors perpendicular to the axis). At high optical depth, the e ect of dichroic extinction becomes more important, preferentially absorbing photons with a negative Q, and leading to a vector pattern dominated by positive Q. The Stokes U component is absorbed more than positive Q, so at high optical depth the vectors have an axial orientation, parallel to the magnetic Úeld. This is described in more detail by [10]. 4.1.3. Stokes V and U : the sense of circular polarisation can reverse several times between the centre and edge of the sphere Non-spherical grains readily produce high levels of CP. In the more familiar case of spherical grains the CP is weaker and the pattern is quadrupolar: the sense of CP alternates between adjacent quadrants of the image of a spherical cloud, with zero CP along the vertical and horizontal diameters. This has been observed in several YSO nebulae (e.g. [35]) and is simply due to alternations in the sign of the Stokes U component as perceived in the scattering plane. At low optical depth, spheroidal grains produce the same quadrupolar pattern as spheres (Fig. 4(a)). However, at high optical depth dichroic extinction plays a prominent role in modifying the Stokes vector after scattering (Eqs. (4), (5)). In the case shown (Fig. 4(b), (c)), the CP produced by scattering is quite high so V=I is large at the edges of the image, where dichroic extinction has little e ect on the polarisation state due to the short path length to the edge of the system after scattering. The CP falls as the line of sight moves inward due to conversion of Stokes V to Stokes U as the path length for dichroic extinction of the scattered light increases. If the path length is long enough (as in Fig. 4(b), (c)) then Stokes V will reverse in sign one or more times between the centre and edge of the image as Stokes V is converted to U , then to negative V , then to negative U , then back to positive V and so on. The behaviour of Stokes U is easy to understand. Dichroic extinction introduces a phase di erence between the electric Úeld components parallel and perpendicular to the magnetic Úeld axis, generating CP. CP is maximised when the phase di erence, , is 90 , declines to zero when =180 , and has the opposite sign for ¿ 180 , peaking at 270 and so on at 450 and 630 , etc. Reversals in Stokes V were also noted by [10] for the case of 2:1 grains, though their origin was not discussed. The e ect of the reversals in U on the LP vector plot (Fig. 3(b)) is fairly subtle: the vector pattern is dominated by Q at high optical depth, so the vectors merely twitch slightly to either side of the vertical as the sign of U changes. While Fig. 4(b), (c) shows the e ects of mutual interconversion of Stokes U and Stokes V initially produced by scattering, there may be cases where scattering produces little CP, e.g. if grains are not aligned, or if the Z41 Stokes matrix element happens to be small. The observed CP may then be produced almost entirely by dichroic extinction by aligned grains in the foreground. In summary, CP can be produced either by scattering of initially unpolarised light (both single and multiple scattering) or by dichroic extinction of linearly polarised light. At high optical depth the latter is likely to be the more important mechanism. 4.2. A young stellar object with a uniform axial magnetic Úeld Model results for an embedded YSO are shown in Fig. 5(a) ­ (d). The circumstellar matter has the following components, illustrated in Fig. 1.


930

P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

(a)

(b)

(c)

Fig. 4. Stokes V and U reversals for a homogeneous sphere. Bright regions correspond to positive U and V , dark regions to negative U and V . (a) at =0:5 a simple quadrupolar pattern is seen. This is similar to that produced by spheres but the CP is much higher, reaching 30%. (b) A denser sphere ( = 3) produces a radial reversal in the sense of CP due to dichroic extinction. CP oscillates between 0% and 20%, near the middle, and exceeds 40% at the edge. (c) Stokes U for the ( = 3) case, illustrating the exchange between U and V , giving 2 radial reversals in the sign of U .

(i) (their radius radius

an accretion disk of radius 400 AU in vertical hydrostatic equilibrium, with density as in [3] Eq. (5)), for the case of steady accretion. For such a disk the volume density declines with as r -15=8 and scale height increases with radius r 9=8 . The disk also has a small inner hole of 2 â 1010 m, which is a factor 20 larger than the adopted stellar radius.


P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

931

(a)

(b)

(c)

(d)

Fig. 5. A YSO with a uniform axial magnetic Úeld (see Fig. 1). The axis ratio of the grains remains 3:1. The YSO is viewed in the equatorial plane and appears bipolar in Stokes I due to the high optical depth of the envelope in that plane. The image (a) shows geometrical limb brightening at the walls of the evacuated cavity. The vector map (b) shows vertical vectors near the middle due to strong dichroic extinction of light in the dense central regions. The pattern is more centrosymmetric near the edges, being dominated by scattering there. (c) and (d) Stokes V and U both reverse sign over short distances in the dense central region. The e ect of the reversals in U can be seen as the radial vectors in (b) approximately halfway from the middle to the edge of the plot. Raw, unsmoothed model output is shown here to reveal the Úne spatial details. CP exceeds 50% at peak.

(ii) a large, equatorially condensed envelope of radius 5000 AU. This has the following density distribution:
env

= CR

-3= 2

(1= (

k

+0:05));

where C is a free parameter (3:7 â 106 Kg m-3=2 in this case), R is the 3D radius for the natural system of cylindrical polar coordinates deÚned by the disk and rotation axis (R = r 2 + z 2 ) and = |z=R|. The index k determines the degree of equatorial condensation and is set to 2.0 here.


932

P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

This density distribution was derived empirically following previous comparison of observations and models [7,8]. Many other algebraic forms could be chosen. (iii) an evacuated conical bipolar cavity with opening angle of 90 and radius of 50 AU in the disk plane. The behaviour of the four Stokes components is similar to that seen in the homogeneous sphere at high optical depth. The main di erences are: (a) that the image is less strongly a ected by the non-spherical shape of the grains, since its bipolar or cometary appearance [4,6,8] is determined more by the envelope and cavity structure of the system (Fig. 1); (b) the radial sign changes in Stokes U and V between the centre and edge of the image are conÚned to a small central region, since the density of the envelope (and hence optical depth) is greatest near the centre of the system. The azimuthal sign changes between quadrants remain prominent however. The radial sign reversals in Stokes U occur at slightly larger radii than those of Stokes V , but would be di cult to observe since the LP pattern is dominated by Stokes Q in the central region. The radial sign reversals in Stokes V are certainly observable in principle via imaging polarimetry at high spatial resolution, e.g. with Adaptive Optics. However, the sign reversals will only be seen if the grains have a fairly large axis ratio, so that the dichroic component of optical depth ( n(s)K34 d s) reaches = 2 before the total optical depth ( n(s)K11 d s) becomes high enough to make observation impossible. This is most likely to occur for dielectric or weakly absorptive substances at size parameters x¡ 1, where the ratio K34 =K11 is high, and may exceed unity. The number of radial sign reversals in U and V is reduced if the maximum size of the silicate grains is increased from 0.35 to 0:75 m (not shown) but the reversals still occur. The grains would also need to be well aligned with a magnetic Úeld which is reasonably uniform on the observed spatial scales. Dichroic extinction appears to be ubiquitous in YSOs, judging by observation of the less embedded systems such as T Tauri stars where the central source is seen directly (e.g. [36]). The degree of LP produced by dichroic extinction is usually ¡ 10% in T Tauri stars, but can be much higher in more embedded systems, exceeding 30% in extreme cases such as IRAS 04361 + 2547 (also known as TMR-1), see [7]. In such systems, it is possible that these sign reversals could be observed, but if not their absence will act as a constraint on models of the dust grain shapes and the density of the nebula. 4.3. Changing the magnetic Úeld: a YSO with a uniform toroidal Úeld The structure of the magnetic Úeld in YSO nebulae is expected to be axial on scales of thousands of AU due to preferential collapse of the parent molecular cloud core along the Úeld lines towards a attened envelope and disk structure. This is conÚrmed by observations of LP produced by emission from non-spherical grains at millimetre and sub-mm wavelengths [18,19], as noted in Section 2. At smaller distances from the protostar the magnetic Úeld is likely to be pinched in by the infall of the circumstellar envelope and may also be twisted by the rapid rotation of the inner accretion disk, e.g. [37,38]. Here we simply show (Fig. 6) that for a uniform toroidal magnetic Úeld structure the LP vector pattern near the centre is rotated through 90 relative to the axial magnetic Úeld case so that the vectors mostly lie in the plane of the accretion disk. This is the expected result for LP dominated by dichroic extinction. In reality, the magnetic Úeld may only become toroidal very close to the centre of the system. Test models (not shown) indicate that in such a case this inner magnetic Úeld


P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

933

Fig. 6. A YSO with a uniform toroidal (horizontal) magnetic Úeld. The vector pattern in the central region is rotated through 90 relative to the axial Úeld case (Fig. 5(b)).

structure can lead to a rotation of the LP vectors in the middle of the image, relative to those further out, providing a useful diagnostic of magnetic Úeld structure. However, in the case of high optical depth (i.e. the younger, more embedded YSOs) this inner structure might be obscured from view by the dichroic extinction of grains aligned with the nearly axial magnetic Úeld in the outer parts of the envelope. This topic will be investigated in future papers by Útting models to data on individual YSOs. 4.4. Implications for the connection with biomolecular homochirality The Stokes V plot for a YSO Fig. 5(c) shows that high degrees of CP can be emitted in some directions, as required for e cient asymmetric photolysis of prebiotic molecules in star forming regions. However, Fig. 5(c) is an infrared simulation. At UV wavelengths, the phase function for scattering by sub-micron sized grains becomes highly forward throwing, and the ux emerging from the axisymmetric model system is dominated by photons which have been de ected through small angles and therefore have weak CP and LP (only a few percent after single scattering). The axisymmetry of the idealised 2D system will also produce zero net CP due to the opposite contributions from adjacent quadrants but this is not a serious problem since many asymmetric nebulae are known, e.g. Cha IRN [39,40]. The hypothesis of Bailey [13] seems to require fairly high (though at present poorly quantiÚed) degrees of UV CP. Hence, a contribution to CP from the dichroic extinction mechanism discussed here may be required, which will be a topic for a future investigation. Another possibility is generation of CP by optically active carbonaceous molecules, perhaps associated with the almost ubiquitous 2175 A interstellar absorption peak. 5. Conclusions Monte Carlo models incorporating dichroic scattering and extinction are likely to prove a valuable tool for the study of dusty nebulae. The results shown here are in good agreement with recently


934

P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

published results using a similar algorithm [10], though the larger axis ratio employed here (3:1 rather than 2:1) leads to some interesting new results at high optical depth, caused by dichroic extinction. The large number of model parameters (e.g. grain size distribution, shape, degree of alignment and any spatial variation in these) means that it will often be necessary to make run a large number of simulations in order to Út observational data. The model polarisation maps shown here contain many complicated features, even though very simple magnetic Úeld structures have been used. Further, abstract investigation of the e ects of non-spherical grains therefore seems warranted, as well as comparison with observation. Acknowledgements I wish to thank Tim Gledhill, Allan McCall and Antonio Chrysostomou for helpful discussions about dichroic scattering. I also thank Barbara Whitney and Mike Wol for invaluable bench testing and discussion of the Monte Carlo codes. The simulations shown here were run on MIRACLE, the SGI Origin 2000 computer at the Hiperspace Centre at University College, London. The author acknowledges support by the UK Particle Physics and Astronomy Research Council (PPARC). Appendix A. The phase function for spheroids The phase function for spheroids may be derived from the Úrst line of the Stokes matrix: Isc = Z11 I + Z12 Q + Z13 U + Z14 V; (A.1) where the (real) Stokes matrix elements Zij ( ; Ú; D; 0) are functions of grain azimuthal orientation angle (deÚned as a right-handed rotation about the Poynting vector relative to the azimuthal de ection angle, ), the grain polar orientation angle Ú and photon de ection angle D. The Stokes matrix elements Zij are formed from non-linear combinations of the complex amplitude matrix elements Sij ( ; Ú; D; ), with set to zero in the Sij and Zij elements. (We adopt the sign conventions of [29] for the matrix elements since these are used in the Mishchenko code.): Z11 =0:5 â (|S11 |2 + |S12 |2 + |S21 |2 + |S22 |2 );
Z13 = -Re(S11 S12 + S22 S21 );

Z12 =0:5 â (|S11 |2 -|S12 |2 + |S21 |2 -|S22 |2 );

Z14 = -Im(S11 S12 - S22 S21 ):

The coordinate system of [30] is used so is a left-handed rotation about the Poynting vector. Since is set to zero in the matrix elements, the orientation of the Stokes vector must be deÚned relative to , with = 0 corresponding to de ection in the plane containing the LP vector and Poynting vector of the incident photon. Again, following [30], the Stokes vector is a left-handed rotation about the Poynting vector. This leads to the relations: Q = I â LP cos(2 ); U = -I â LP sin(2 ); V = I â CP; (CP = V=I ): (LP = Q2 + U 2 =I );


P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937
Polarisation vector

935

Azimuthal grain orientation


Azimuthal deflection direction

Poynting vector of scattered photon D

Poynting vector of incident photon

Fig. 7. Illustration of some of the angles associated with the phase function. The polarisation vector, azimuthal grain angle and azimuthal de ection angle lie in the plane perpendicular to the incident photon's Poynting vector. is a right-handed rotation but is a left-handed rotation.

These relations between the three azimuthal angles ( , frame of the incident photon are illustrated in Fig. 7. Eq. (A.1) then rearranges to

and direction of LP) in the reference

Isc d = Iin {(1 - LP )= 2 â (|S11 |2 + |S12 |2 + |S21 |2 + |S22 |2 ) +LP â [(|S22 |2 + |S12 |2 ) sin2 ( )+(|S11 |2 + |S21 |2 ) cos2 ( )
+sin(2 ) Re(S11 S12 + S22 S21 )] - CP â Im(S11 S12 - S22 S21 )} d

(A.2)

giving the scattered intensity per unit solid angle d . Eq. (A.2) may also be derived directly from the amplitude matrix by separating the electric Úeld amplitudes into components parallel and perpendicular to the scattering plane, after the method used by [5] to derive the simpler phase function appropriate for Mie scattering by spheres. Eq. (A.2) is equivalent to the following probability density function: P (D; )d D d = sin(D){(1 - LP )= 2 â (|S11 |2 + |S12 |2 + |S21 |2 + |S22 |2 ) +LP â [(|S22 |2 + |S12 |2 ) sin2 ( )+(|S11 |2 + |S21 |2 ) cos2 ( )
+sin(2 ) Re(S11 S12 + S22 S21 )] -CP â Im(S11 S12 - S22 S21 )} d D d =I max

( ; Ú; LP; CP );

(A.3)


936

P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

where Imax ( ; Ú; LP; CP ) is the intensity at the most probable scattering direction (D; ) which is di erent for each combination of ( ; Ú; LP; CP ). This term is used to normalise the probability density function relative to unity. Eq. (A.3) is sampled by rejection sampling: values of both D and are randomly generated and then accepted if P (D; ) is greater than a third random number R(0 ¡R¡ 1). This method requires calculation of Imax ( ; Ú; LP; CP ) for each possible combination of incident ( , Ú, LP and CP). Eq. (A.2) must often therefore be calculated for all values of all six variables to do this, and the peak values stored in a look-up table. For large size parameters, it may be assumed that the most probable de ection angle will always be small, which reduces the number of calculations required. Rejection sampling is ine cient for highly forward throwing phase functions. However, even at size parameters up to 20 the run time of the code is not a ected, since it is dominated by the computation of the photon location and Stokes vector as it moves in many small steps along its path, not by the scattering calculation. Eq. (A.3) is usually sampled every 2 in each of the 4 angles, every 2% in LP and every 20% in CP, which generally has only a small e ect on the phase function via conversion of some of the incident Stokes V to the direction-dependent Stokes U . This leads to 150 billion calculations of Eq. (A.2), which takes about an hour on an SGI Origin 2000 computer, running with several CPUs in parallel. This speed is achieved by simplifying Eq. (A.2) by precalculating components of the equation at the same time that the matrix elements are Úrst calculated using a slightly modiÚed version of the Mishchenko code, ampld.new.f. For small size parameters (x¡ 1) the 4 angles and LP may be sampled a factor of 2 less often, which reduces the number of iterations of A.2 by a factor of 32. References
[1] Woods JA. Winds in cataclysmic variables. D. Phil thesis, University of Oxford, 1991. [2] Knigge C, Woods JA, Drew JE. The application of Monte Carlo methods to the synthesis of spectral line proÚles arising from accretion disc winds. Mon Not R Astron Soc 1995;273:225. [3] Whitney BA, Hartmann L. Model scattering envelopes of young stellar objects. I--method and application to circumstellar disks. Astrophys J 1992;395:529. [4] Whitney BA, Hartmann L. Model scattering envelopes of young stellar objects. II--infalling envelopes. Astrophys J 1993;402:605. [5] Fischer O, Henning Th, Yorke HW. Simulation of polarization maps. 1: protostellar envelopes. Astron Astrophys 1994;284:187. [6] Fischer O, Henning Th, Yorke HW. Simulation of polarization maps. II. The circumstellar environment of pre-main sequence objects. Astron Astrophys 1996;308:863. [7] Lucas PW, Roche PF. Butter y star in Taurus: structures of young stellar objects. Mon Not R Astron Soc 1997;286:895. [8] Lucas PW, Roche PF. Imaging polarimetry of class I young stellar objects. Mon Not R Astron Soc 1998;299:699. [9] Wolf S, Voshchinnikov NV, Henning Th. Multiple scattering of polarized radiation by non-spherical grains: Úrst results. Astron Astrophys 2002;385:365. [10] Whitney BA, Wol M. Scattering and absorption by aligned grains in circumstellar environments. Astrophys J 2002;574:205. [11] Silber J, Gledhill TM, Duchne G, Menard F. Near-infrared imaging polarimetry of the GG tauri circumbinary ring. Astrophys J 2000;536:L89. [12] Krist JE, Stapelfeldt KR, Watson AM. Hubble space telescope/WFPC2 images of the GG tauri circumbinary disk. Astrophys J 2002;570:785.


P.W. Lucas / Journal of Quantitative Spectroscopy & Radiative Transfer 79­80 (2003) 921 ­ 937

937

[13] Bailey J, Chrysostomou A, Hough JH, Gledhill TM, McCall A, Clark S, Menard F, Tamura M. Circular polarization in star-formation regions: implications for biomolecular homochirality. Science 1998;281:672. [14] Chrysostomou A, Gledhill TM, Menard F, Hough JH, Tamura M, Bailey J. Polarimetry of young stellar objects--III. Circular polarimetry of OMC-1. Mon Not R Astron Soc 2000;312:103. [15] Menard F, Chrysostomou A, Gledhill TM, Hough JH, Bailey J. High circular polarization in the star forming region NGC 6334: implications. In: Lemarchand G, Meech K, editors. Bioastronomy 99: a new era in the search for life in the universe. Proceedings of the ASP Conference, vol. 213. San Francisco: ASP, 2000. p. 355. [16] Clark S, McCall A, Chrysostomou A, Gledhill T, Yates J, Hough J. Polarization models of young stellar objects--II. Linear and circular polarimetry of R Coronae Australis. Mon Not R Astron Soc 2000;319:337. [17] Aspin C, Casali MM, Walther DM. Infrared imaging and polarimetry of star forming regions. In: Reipurth B, editor. Low mass star formation and pre-main sequence objects. ESO Conference Workshop Proceedings, vol. 33. 1989. p. 349. [18] Aitken DK, Wright CM, Smith CH, Roche PF. Studies in mid-infrared spectropolarimetry. I--magnetic Úelds, discs and ows in star formation regions. Mon Not R Astron Soc 1993;262:456. [19] Tamura M, Hough JH, Hayashi SS. 1 millimeter polarimetry of young stellar objects: low-mass protostars and T tauri stars. Astrophys J 1995;448:346. [20] Aitken DK, Efstathiou A, McCall A, Hough JH. Magnetic Úelds in discs: what can be learned from infrared and mm polarimetry? Mon Not R Astron Soc 2002;329:647. [21] Purcell EM. Suprathermal rotation of interstellar grains. Astrophys J 1979;231:404. [22] Roberge WG. Grain alignment in molecular clouds. In: Roberge WG, Whittet DCB, editors. Polarimetry of the interstellar medium. ASP Conference Series, vol. 97. 1995. p. 401. [23] Lazarian A. Mechanical alignment of suprathermal grains. In: Roberge WG, Whittet DCB, editors. Polarimetry of the interstellar medium. ASP Conference Series, vol. 97. 1995. p. 425. [24] Lazarian A. Physics of grain alignment. In: Franco J, Terlevich L, Lpez-Cruz O, Aretxaga I, editors. Cosmic evolution and galaxy formation: structure, interactions, and feedback. ASP Conference Proceedings, vol. 215. 2000. p. 69. [25] Hough JH, Aitken DK. Polarimetry in the infrared: what can be learned? JQSRT 2002;this issue. [26] Cronin JR, Pizzarello S. Enantiomeric excesses in meteoritic amino acids. Science 1997;275:951. [27] Whitney BA. Thomson scattering in a magnetic Úeld. PhD thesis, Wisconsin University, Madison, 1989. [28] Martin PG. Interstellar polarization from a medium with changing grain alignment. Astrophys J 1974;187:461. [29] Mishchenko MI, Hovenier JW, Travis LD. Light scattering by nonspherical particles: theory, measurements, and applications. San Diego: Academic Press, 2000. [30] Hovenier JW, van der Mee CVM. Fundamental relationships relevant to the transfer of polarized light in a scattering atmosphere. Astron Astrophys 1983;128:1. [31] Voshchinnikov NV, Karjukin VV. Multiple scattering of polarized radiation in circumstellar dust shells. Astron Astrophys 1994;288:883. [32] Code AD, Whitney BA. Polarization from scattering in blobs. Astrophys J 1995;441:400. [33] Preibisch Th, Ossenkopf V, Yorke HW, Henning Th. The in uence of ice-coated grains on protostellar spectra. Astron Astrophys 1993;279:577. [34] Pollack JB, Hollenbach D, Beckwith S, Simonelli DP, Roush T, Fong W. Composition and radiative properties of grains in molecular clouds and accretion disks. Astrophys J 1994;421:615. [35] Chrysostomou A, Menard F, Gledhill TM, Clark S, Hough JH, McCall A, Tamura M. Polarimetry of young stellar objects--II. Circular polarization of GSS 30. Mon Not R Astron Soc 1997;285:750. [36] Sato S, Tamura M, Nagata T, Kaifu N, Hough J, McLean IS, Garden RP, Gatley I. Infrared polarimetry of dark clouds. II--magnetic Úeld structure in the Rho Ophiuchi dark cloud. Mon Not R Astron Soc 1988;230:321. [37] Galli D, Shu FH. Collapse of magnetized molecular cloud cores. I. Semianalytical solution. Astrophys J 1993;417:220. [38] Galli D, Shu FH. Collapse of magnetized molecular cloud cores. II. Numerical results. Astrophys J 1993;417:243. [39] Ageorges N, Fischer O, Stecklum B, Eckart A, Henning Th. The Chamaeleon infrared nebula: a polarization study with high angular resolution. Astrophys J 1996;463:L101. [40] Gledhill TM, Chrysostomou A, Hough JH. Linear and circular imaging polarimetry of the Chamaeleon infrared nebula. Mon Not R Astron Soc 1996;282:1418.