The evolution of the zodiacal dust cloud under plasma drag and lorentz forces in the latitudinally asymmetric solar wind

Planet. Space Sci., Vol. 43, Nos. 314, pp. 301-312, 1995 Copyright 0 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0032-0633/95 $9.50+0.00

Pergamon

0032-0633(94)00174-X

The evolution of the zodiacal dust cloud under plasma drag and Lorentz forces in the latitudinally asymmetric solar wind H. J. Fahr,’ K. Scherer’ and M. Banaszkiewicz’ ’ Institute for Astrophysics and Extraterrestrial Research of the University of Bonn, Auf dem Huge1 71, D-53 121 Bonn, Germany * Space Research Center of the ‘Polish Academy of Sciences, Bartycka 18A, 00-716 Warsaw, Poland Received 16 March 1994 ; revised 25 July 1994 ; accepted 28 July 1994

Introduction The solar system is surrounded

by a cloud of zodiacal

dust

particles. This dust configuration has been extensively investigated by observations of the zodiacal light distribution (see, e.g. Leinert, I975 ; Levasseur-Regourd and Dumont, 1980; Giese, 1980; Leinert ef al., 1978, 1981). Recently ail-sky infrared observations have also been carried out with the infrared astronomical satellite (IRAS) (Reach, 1991). The same zodiacal dust particles with radii in the range between 20 and 200 pm, which are responsible for the visible zodiacal light, also produce the infrared zodiacal emission. The faint optical light emission is generated through scattering of solar photons at the dielectric dust particles orbiting in quasi-Keplerian orbits. Due to its origin through this scattering of solar photons at orbiting dust particles, measurements of the optical zodiacal light permit the derivation of the zodiacal dust density distribution and the underlying distribution of dust particle orbits, i.e. of inclinations, eccentricities and semimajor axes (see, e.g. Haug, 1958; Leinert et al., 1976; Fahr et al., 1981; Buitrago et al., 1983 ; Fahr and Ripken, 1985 ; Kneissel and Mann, 1991; Banaszkiewicz et al., 1994). The observed dependence of the dust density n(r) on heliocentric distance according to n(r) CC(v/r)” with p > 1 was confirmed by Leinert et aZ. (1981), Toller and Weinberg (1985) and Reach (1991), and can most probably be explained only by non-uniform dust sources in combination with a Poynting-Robertson migration towards the sun (Dohnanyi, 1978; Burns et al., 1979; Leinert, 1985). The derived latitude-dependent density profiles showing a deficiency (Schwehm and Rohde, 1977 ; Mukai and Yamamoto, 1979), or an excess (Buitrago et al., 1983) of dust particles at high latitudes (e.g. see Kneissel and Mann, 1991) are generally ascribed to the source distribution. However, density distributions might also be caused to vary with latitude for the following reasons:

302

H. J. Fahr et

orbital inclinations vary due to any existing normal component of forces. Particles also might move with radial drift rates depending on inclination. The first process was discussed on the basis of a Lorentz scattering (Mukai and Giese, 1984), however, the resulting effect is found to be fairly small. No attempts up to the present can be found in the literature where latitude-dependent forces affecting the evolution of dust orbits at different inclinations. Such a force, however, results from the plasma Poynting-Robertson effect of dust grains interacting with the three-dimensional solar wind flow. This force is proportional to the momentum transferred from the plasma to the grain. Therefore, the asymmetry in the solar wind results in an asymmetric force. Observations of solar-wind plasma scintillations (Kojima and Kakinuma, 1987) or of interplanetary Lyman-alpha glow asymmetries (Lallement et al., 1985; Ajello, 1990) clearly permit the conclusion that the mass outflow from the solar poles, at least during the solar minimum periods, is smaller by about 3&50% than that from ecliptic regions. Thus it can be concluded that the plasma Poynting-Robertson force at lower latitudes is weaker decreasing the Poynting-Robertson drift rates at low inclinations compared to high inclinations. In this paper, we discuss the effect of plasma Poynting-Robertson forces which are latitude-dependent. In addition we study Lorentz scattering effects and systematic Lorentz drifts in competition to this plasmainduced effect. The motion of dust particles under the influence of the above-mentioned forces is described by orbit-averaged equations yielding the evolution of the orbital elements of the dust particles. We also solve the kinetic equation for the distribution function of particles in orbital element space following the approach by Leinert et al. (1983) and develop an expression for the spatial density of the zodiacal dust cloud as a function of latitude.

The plasma drag force The total force K acting on orbiting dust grains is composed of contributions by solar gravity KK, the electromagnetic Poynting-Robertson force K,, the Lorentz force KL, and the solar-wind-induced plasma Poynting-Robertson force KP (see, e.g. Burns et al., 1979) and thus for dust particles leads to the following equation of motion : mii: = KK+K,+KL+KP=

-t y

forces

The last force, K,, in equation (1) denotes the plasma Poynting-Robertson force ; it comprises momentum transfers by adsorption and specular or diffuse reflection of solar-wind protons and cr-particles, as well as sputtering of atoms or molecules from the grain. This force KP thus results from ions impacting on the surface (47~‘) or the cross section (7~‘) of the dust grain. One can treat this ion impact in a geometrical approximation though the dust grains are electrically charged, and their Debye-length AD is large compared to its radius s (i.e. s2 <
KPW =

W(V,,,W,r>(7ts2)U,,l(v)lmiv,,(v) d3u. (2) sIw3

To calculate the integral we choose a coordinate system which is co-moving with the solar wind plasma bulk because in this system the distribution function of the plasma ions is adequately approximated by a bi-Maxwellian. The individual relative velocity with respect to the solar wind rest frame is given by w = v,,-v

(3)

where v is the velocity of an individual ion and v,, is the solar wind bulk velocity. Then the solar wind ion distribution function is represented by : f(r,w) =n(r)

GM,m -r

&,[(I- i)r-

a/. : Evolution of zodiacal dust under asymmetric

x exp iv]+L(,+K,,

(1)

where G is the gravitational constant, M, and m are the masses of the sun and of the dust particle, respectively. S, = 1.36 x lo6 erg cmp2 s-’ denotes the solar electromagnetic energy flux. A = K? is the geometrical cross section of the particle with radius s, c is the velocity of light, and Qpr is the radiation pressure coefficient which is calculated by number at a later part of this paper (for its derivation see, e.g. Burns et al., 1979). The vectors r and v describe the position and velocity of the particle, respectively.

= n(r)%/&’

-

miw2 cos’ /I

miuj2 sin2 p

2kTi,,

2kT,,

exp (-

$$(coZ

8+Esin2B))

(4)

n(r) is the ion density. Tilland Ti, are the local ion temperatures parallel and perpendicular to the local magnetic field B. The angle between B and the velocity w is denoted by fl. The quantity d is the temperature anisotropy defined by B = Ti,,/Ti,. With U being the velocity of the dust grain with respect to the Sun, its relative velocity vector with respect to the

H. J. Fahr et al. : Evolution of zodiacal dust under asymmetric forces

solar wind rest frame is given by u = v, -U. The ion velocity relative to the dust grain is obtained by V,,I= v-u

= v,,-w-u

= u-w

(5)

with its modulus given by o,,, = J~&21(wcosB,

(6)

where 0 is the angle between the velocity vectors w and u. With polar coordinates rr 9, # and a solar wind ion distribution represented by a mono-Maxwellian (i.e. assuming T,,, = TiL ; d = 1) one obtains : 2 f(w,

r)

=

n(r)(xc:)

- 3’2

exp

-K (

, C:

(7)

)

where cf is the mean ion thermal speed. With the above representations equation (2) attains the form K, = - AcDeu

(8)

with A=

&n(r)m,s’cf.

(9)

and c,, denoting the drag coefficient and given by

+A4

(

l+&-

1 4M4 )

erf(M)

1

(10)

Here M = v,,/cj is the solar wind sonic Mach number. Its average value at solar distances of relevance here, i.e. equal or larger than 1 AU, is M 2 l&15. Thus in this regime as shown by equation (10) the drag coefficient cD is practically independent on the Mach number of the plasma flow and hence on the sound velocity c,. The above result in equation (10) asymptotically approaches that one given by Morfill and Griin (1979) for such large Mach numbers. It differs, however, from their result at low Mach numbers, i.e. in the subsonic solar wind regime. This we have explicitly demonstrated in a different paper (Banaszkiewicz et al., 1994). As is expressed in equation (8) the vector orientation of K, is co-planar to the orbital plane of the dust particle, i.e. no normal component of the force exists. This result in equation (8) is a consequence of the use of a monoMaxwellian for the solar wind ion velocity distribution. Only in this case the force KP turns out to be collinear to II, i.e. coplanar with the orbital plane of the dust particle. Normal components of KP would exist if bi-Maxwellian velocity distributions with different thermal dispersions parallel and perpendicular to the interplanetary field are used instead. For regions close to the sun (i.e. r < 0.3 AU) this representation is known to yield much better fits to the observations (Marsch et LIZ.,1982).

The Lorentz force acting on charged dust particles

Lorentz forces act on electrically charged dust particles which move through the interplanetary magnetic field.

303

Due to the frequent changes in the field direction in regions close to the interplanetary current sheet, i.e. within about + 15” from the ecliptic, these forces have the character of stochastic electromagnetic forces and here most easily can be treated as giving rise to random diffusion processes (see, e.g. Hassan and Wallis, 1983 ; Pellat et al., 1984; Morfill et al., 1986). The diffusion rates, however, are only relevant for dust particles with radii s 6 1 pm. Unfortunately, there is no clearcut answer in the literature how important these diffusive drift motions are compared to Poynting-Robertson migration rates. Following Leinert (1985) or Morfdl et al. (1986) the root-mean-square variation of the inclination of dust at 4 AU, with s = 1 pm and charged to 6 V, amounts to 4” in 10 years, whereas, according to Morfill and Griin (1979) or Barge et al. (1982) it amounts to only 0.25” or 0.07”, respectively. Based on this latter result the change in inclination during an associated Poynting-Robertson lifetime for such a particle would remain smaller than 3”. Recalling that the size-dependence of the Lorentz force roughly follows the relation KL N (s~/s)~with 1.5 < p < 2.0 (Morfill and Gri.in, 1979) it can be estimated that dust particles with s > 10 pm, which are just the ones relevant for the zodiacal light, experience inclination changes even smaller than a factor of 100 compared to those mentioned above. Thus we can conclude that migration rates (di/dt) induced by Lorentz scattering of (s 2 10 pm) particles within heliographic latitudes - 15” 6 9 < + 15”, i.e. with inclinations i < 15”, are completely negligible (see also Consolmagno, 1979; Hassan and Wallis, 1983; Pellat et al., 1984; or Morfill et al., 1986). On the other hand, particles with i 2 15” are subject to systematic Lorentz forces caused by the different polarities of the interplanetary field at northern and southern ecliptic hemispheres. The very general description of dust particle motions in such high inclination orbits under the influence of Lorentz forces induced by a dipolar field with Parker’s Archimedian spiral configuration is fairly complicated, but here we can first give approximate descriptions of these motions. Radial, azimuthal and latitudinal components of this Parker field (Parker, 1963) are given by : B = B,&(~,O,~)

= B,e,+B,e,

(11)

where c = f 1 regulates the polarity of the northern or southern ecliptic field, and where R is the angular velocity of the rotating sun. B,, is the magnetic field magnitude at r=r,=lAU. The Lorentz force of a dust particle with charge q is given by : K,_ = :(v-v,)

A B.

(12)

To take into account the above Lorentz force at an integration of the equation of motion (1) requires a somewhat laborious decomposition of this force in radial, tangential and normal components which we shall derive in a later paragraph. There is, however, a simple argument which can demonstrate clearly that the above Lorentz-forces K,

H. J. Fahr et al. : Evolution of zodiacal dust under asymmetric forces

case= Jl

i sin* (o + cp) = A.

-sin*

(15)

The projection of the Lorentz force KL given in equation (13) onto the normal vector C, of the orbital motion gives the component which is solely of interest in our case here k,,

=5(&v,,)

A B-Co

= ___ mck

Fig. 1. Schematic illustration of symmetry properties of a circular dust particle orbit and use of special coordinates to denote space and velocity vectors. For further details see text

when integrated over a complete orbit exactly cancel for circular and quasicircular orbits (a << 1). Take a circular orbit with inclination i and intersect it with two planes 9,,Z parallel to the ecliptic with distances +z and -z, respectively. The two intersection points Y,, r2 with the upper plane have space vectors r,,2 = z f p, and those r3, r, with the lower plane r3,4 = -zf p, where fp are vectors with equal moduli but opposite directions in the cutting lines of the orbital plane with the planes .Y1,?, respectively (for illustration purposes see Fig. 1). It is then evident from the given symmetries and from the fact that the field polarity cr in the northern and southern hemispheres are opposite implying that B, changes its sign (whereas B, is identical), that the Lorentz forces in points Ye and Ye, as well as those in points rf and r, exactly cancel each other (see Fig. 1). This is easily seen from the expression (12) for the Lorentz force when written in its polar coordinate components KL = (qlcM@,)

+e&J%

- %J$J

+e,(

-

k,iv =

(14)

with 8 being the heliographic latitude, i the inclination of the orbit, w the longitude of the perihelion and cp the true anomaly. From this relation one easily deduces

(16)

qB,r~Qn& mcv,,(l -E*)“~

sin cp

= (a:sincp-a$)(l+scosq)A. Furthermore, the Gaussian perturbation equation secular changes of the inclination i reads : di

(17) for the

~cos(u+~)~ nu

Zt-

l+&COScp

The average of the above quantity is given by :

LN

(18)

over one orbital period

(19) which leads to (l-&2)2

(13)

sine = sinisin(o+cp)

A k)

where p = m andp = a (1 - 8’) with the mean motion n, the semimajor axis a and the eccentricity E, respectively. After some lengthy vector algebra (see Scherer, 1995) we are finally left with

@,)I

and when keeping in mind that the proper vectors e,, e,, e+, are just opposite to each other in points Y, and r3, or points r2 and r4, respectively. This argument now can, of course, analogously be applied to all other points at the orbit generated by planes cutting at different values z, and hence one can conclude that the Lorentz forces when integrated over the complete orbit vanishes, leading to the result that normal forces are ineffective and thus (di/dt) = 0. For highly eccentric orbits with eccentricities E >> 0 this argumentation does not strictly hold because the orbitintegrated action of the Lorentz forces [TJK,(s) ds/v(s)] in these cases does not vanish exactly. In the following we derive the average drift rate in the inclination i due to the action of the normal component of the Lorentz force. For a point on the sphere the following relation is valid

((t A B) - (v,, A B))*(r

* *a cos (w + cp) sin cp

=271nau1

X

X

Jl

2n

so

-sin*

(1 +ECOS(p)* (l-E?)* xuf

isin* (o+cp)dq-

cos(w+cp)

p+‘,mdq. (1 +ECOS~)2 s0

(20)

We expand the denominator function under the integral into a Taylor series about E = 0 and, thereby, obtain 1

= l-~ECOS~~+~~*COS*(P...

(l+&COS(p)* =~~o(-l)“(n+l)E~COS~(p. Introducing

(21)

q = o + cp, and

k = sini, u,

k’=J1-k2,

A=,/-,

&l!Yu* 2nnll

”

a2

(1-E2Y 2*

=75&y

(22)

H. J. Fahr et al. : Evolution of zodiacal dust under asymmetric forces

.-a

.4

____.__._._.____....-

:

: i.--...---..--.---.--.--.-.--.----.---.-------.--...-

.A

+

i

.d

ad

---_______------___-_--l 11 50

: .i

F

30.01

..-I

t

”

0

Do--

100

Dynamical

time

200

Dynomicol

step

Fig. 2. Shown are the inclination drifts of dust particles for four different initial values (i = 0” solid, i = 30” dashed, i = 60” i = 90” dashed-dotted). initial parameters : A, = -0.2;

100

0

200

time

step

Fig. 3. Shown is the inclination drift for a dust particle with an initial value of i = 30” including the polarity changes of the magnetic field with a period of 11 years (solar cycle). The legend is the same as in Fig. 2

particle diameter

is:s=8pm

then leads to the following order in E

type of integrals

up to the first

2n+w I=

sin(v,l-o)}

cosq{a,--a,

sw x [ 1-2~~0s

(u--)]Ady

(23)

which can be solved in terms of elliptic integrals, and then leads to the results : di

(->

= -:(2saZcosw-a,

sino)

dt P

x

(

cotg2

solar wind outflow from the corona. Even though the Ulysses probe is just now passing towards high ecliptic latitudes and is at present collecting preliminary data from there, it is unclear how much latitudinal variations are entangled with solar-cycle-induced variations. Furthermore, the most recent Ulysses data are not yet published in elaborated form. Nonetheless, from interplanetary scintillation measurements (Kojima and Kakinuma, 1987) one can deduce that the solar wind bulk velocities at r, = 1 AU are typically 7&100% higher at polar regions compared to ecliptic regions. This means uSw= 400 km s-’ in the ecliptic, and between 600 and 800 km SC’ over the poles. Thus in a simple approximation we take the following representation :

l+sin2i iK(k)-

------E(4 i

sin’

(24)

u,JS) = a,,(1 +B, sin’s),

(25)

)

with the complete elliptic functions K(k) and E(k) of the first, and of the second kind, respectively. For a critical value o, elf w the first bracket on the right hand side of equation (24) vanishes. When passing through this critical value o, the sign of the i-drift is reversed and thus the migration in the inclination is turned to the opposite sense. In Fig. 2 we have shown which changes in inclination i are suffered by dust particles according to the action of systematic Lorentz forces. It is evident that due to field polarity changes within the solar magnetic cycle of about 22 years this accumulated e:ffect in the i-drift will in fact be even much smaller because drift rates will be reversed every 11 years. This we have demonstrated in Fig. 3. In view of these results we may thus simply neglect the effect of Lorentz drifts in the description of secular changes in the zodiacal dust cloud.

Modelled solar wind asymmetries

Up to now there is a lack of good and direct observational knowledge on the 3-dimensi.onal dynamic structure of the

where r&,,is a reference at r, = 1 AU, and 9 = O”, with 9 being the ecliptic latitude. The amplitude B,of the variation also may periodically change with the phase of the solar cycle. For maximum and minimum conditions, one may select the following values for the parameters B,and - . %v .

B,(Max) = 0.5 ;

a&Max) = 400 km SC’

(26)

B,(Min) = 1.0;

D,(Min) = 350 km s-r.

(27)

Furthermore, based on Lyman-alpha glow asymmetries (Kumar and Broadfoot, 1979 ; Witt et al., 1979 ; Lallement et al, 1985; and Ajello, 1990) the solar wind mass flow appears to vary with ecliptic latitude according to : aqr, 9) = &,n~

4 - (1 +A2 sin’ 9), 0r2

(28)

where n,, is a reference value at ro. In a more general form the mass flow can be represented by a Legendre polynomial series : 1+ 2 AZnsin2” 9

@(v,9) = $ (

n=l

(29) )

H. J. Fahr et al. : Evolution

306

in which at present terms describing “even asymmetries be taken account, since direct north-south solar asymmetries not yet nor do up to require any north-south Nevertheless, the detection north-south monthly averaged solar emission features (Pryor et al., 1992) may indicate the existence of such asymmetries also in the solar wind parameters. Longitudinal variations in the solar wind flow which would be associated with the solar rotation period or with coronal activity phenomena of even shorter periods are neglected here. The solar rotation period r,, = l/(27&) is much shorter than the mean orbital period of the dust particles. Hence one can safely neglect wind structures of this origin, assuming that longitudinal asymmetries associated with it are smeared out and only yield an averaged effect on dust particles. When treating the equation of motion for the dust grains we will keep only the first two even terms (AZ, A4) of the series expansion given in equation (22). The following typical values for the solar maximum and solar minimum conditions are taken : A*(Max) = -0.3;

n,(Max) = 6 cm-3

A,(Min) = -0.4;

n,(Min) = 5 cmm3.

This set of parameters nicely fits the white light coronograph observations of Munro and Jackson (1977) implying a coronal density decrease with latitude by a factor of 0.3 from the ecliptic to the poles. For the high Mach number supersonic solar wind regime which we are going to treat in this paper (I 2 1 AU) no information on the interplanetary ion temperatures is needed for the calculation of cn and of the plasma Poynting-Robertson effect.

of zodiacal dust under asymmetric forces = FlompAu2r~cD.

YSW

(33)

With the use of the vector decomposition

;=,,=++e ( ) v,,-v

u

V

where Kc=

2

(0

-2--

1-J

u

(v,;v)

V SW

%v

VSWV

-Ii2

(35) )

.

In the following calculation we keep terms of the first order in (i/vSw) and thus develop (l/rc) at v = 0 into a Taylor series up to this order arriving at the following result k: 2 1 +&w.

(36)

Insertion of equation (35) then brings equation (34) into the form u/u = e, - (rp/Qe,,

(37)

where /I is the azimuth angle in the orbital plane of the dust particle. We decompose the total force per mass into an effective Keplerian term K; and perturbing radial terms, K,, and tangential terms, K,, by : K,,, = K;+K,+K,

(38)

with K, = (~‘-np(ro)mpA~2r~~~)

K, =

Ler r2

-2y,le,

(41)

Introduction of the plasma Poynting-Robertson force KP from equation (8), together with the drag coefficient cn derived in equation (lo), into the equation of motion (1) one then with the neglect of KL is led to the following equation

_p’f-.y j-3

Y+ysw ;-2y,:;,

IC

where e, and e, are tangential and radial unit vectors, respectively. We have assumed that the solar wind flow, i.e. the vector v,,, is radial. From equations (25) and (28) one derives the following expression for n,(ro)

(30)

1 + A, sin’ 9 nJr0)

where the dots on top of the space vector r denote derivatives with respect to time and where the coefficient p’ of the effective Kepler force is represented by

Sor’A Q PI

GM,mc

(31)

&!g

The plasma Poynting-Robertson the following drag coefficient

epr.

=

(32)

force is connected with

n0

(42)

l+B,sin’S’

The ratio of the radial solar wind drag force to the effective gravity attains the order : np(ro)mpAu2r~cD P’

and where the coefficient of the electromagnetic PoyntingRobertson drag force is given by Yr =

(39)

(40)

Y2C

The effective equation of motion

)“2i:=

(34)

3

N 2x lo-‘5-

1 m’

(43)

where m is the mass of the dust particle in g. As discussed earlier the mass flow shows a latitudinal dependence of the form Q = O,,(l +A2 sin’ $+A4 sin4 S),

(44)

where Q0 is a reference value at r. = 1 AU and 9 = 0”. With this flow we find

H. J. Fahr

Y sw

do

= @AK~c~.

dt

The relative importance of the plasma Poynting-Robertson effect is then revealed by the coefficient q = (c/v,,) (Y&J which for quasicircular orbits with (r/u,,) < 0.01 can be estimated by y-

%V Yr

(>

(45)

K4 r2

yI _ 5

307

: Evolution of zodiacal dust under asymmetric forces

t 1: 0.386 - c2vvsw CT)

QP

SoQpr

’ - lTznE

{A\ sin’ i+ AI, sin4 i>_ (55)

‘p

Here the functions H, = H,(E ; i, w) and H, = H,(i, o) were introduced and the quantities H,, and H,, denote H,,=H,(~=O;i,0)=2(1+v])+A\sin~i+~A’,sin~i (56)

CZ0.5,. . ,0.7. Ha2 = A\ sin’ i(i-icos2

o)+A’, sin4 i(i--icos2 o)

(46) Here geometric optics instead of more sophisticated Mie scattering results is used yielding Qpr = 1. This reveals that the forces K, and Kp turn out to be of about the same order. Inserting equations (45) and (44) into equations (40) and (41) the radial and tangential components are gained in the form (47)

Kt = - 7

[Ab + 11: sin2 9 + A\ sin“ $1

(48)

with CI= : = 3.55 x IO-* Qpr (AU2/yr)

(49)

A’, = 2

(50)

A; ==l+q

(51)

AI=qAi

i=2,4.

(52)

The perturbing forces K, and K, are much smaller than the effective gravity force KK. Thus one can study the orbit evolution using Gauss’ equations (see Neutsch and Scherer, 1992). To follow these procedures we refer to the classical orbital elements a, E,i, cu (semimajor axis, eccentricity, inclination and argument of the perihelion, respectively). Secular changes in the inclination and in the longitude of the node 0 are proportional to the normal component of the force. In our approach this component vanishes, unless the solar wind ion distribution has to be represented by a bi-Maxwellian. In addition, since changes in the other elements integrated over one orbital revolution are very small we may average the Gauss equations over one revolution period. Thus we obtain (see Banaszkiewicz et al., 1994, for more details) da

(> dt

‘1 =

Ha(E;&W)

-_____

a(1 -,E~)~‘~

‘p

= ---

x

(&I + ~~fL2>

a(1 -.E~)~‘~ de - s&H,(i, Z @=

0

(57) He=

-a2\

{:+$+A\ 2

sin2 i(i-icos’o)

+A: sin’ i(G--acos2 CD)>. (58) In Fig. 4(a)-(c) we have shown calculations of the secular evolution of (a(t)), (E(t)) and (w(t)) carried out on the basis of the above formulae (53)-(55) for dust particles at different inclinations i. The calculations were started with the initial element conditions a, = 1 AU, c0 = 0.2, w0 = 45”. In addition we used the values A2 = -0.2, A4 = -0.02 and s = 8 pm. For four values of i (i = O”, 30”, 60”, 90”) the element evolutions are shown over a total evolution period of 19,600 years, i.e. the PoyntingRobertson lifetime for 8 pm dust particles at 1 AU. One may notice that the changes in time of the semimajor axis a(t) and of the eccentricity c(t) are less pronounced at increasing inclinations i, i.e. for particles spending more time per orbit at higher latitudes. Hereby, especially the evolutionary period in the argument of the perihelion, co(t), shows a pronounced dependence on i. One can also recognize a tendency of low inclination particles to circularize their orbits faster and to reduce their main axes more rapidly compared to high-inclination ones indicating the fact that the probability to find particles with higher inclinations increases with decreasing main axes. The kinetic equation and the density of the zodiacal dust

To describe the statistical properties of the zodiacal dust we use a distribution function f(a, E,i) depending on the three orbital element variables : semimajor axis, a ; eccentricity, E; and inclination, i as already done by Leinert et al. (1983). This is very much different from the distribution function given as a function of phase space coordinates as it was used by authors like Hassan and Wallis (1983) or Pellat et al. (1984). Following Leinert et al. (1983) we shall also keep for the sake of generality their function g(c) here distributing the dust sizes s. The kinetic equation for the steady state distribution function of dust particles, as a kind of Liouville equation in element space, can then be given in the form :

div (u&f@, e, 9) = NUT’@,E, i) g”(s)(53)

&f@, z

(n

E, 4 g(s)

Ei

s)

.

C 3 32 (59)

a)

(54)

Here dlv is the divergence operator expressed in coordinates of element space and ct is the semimajor axis nor-

308

H. J. Fahr et al. : Evolution

Dynomicol

(a)

time

of zodiacal

dust under asymmetric

forces

and on g(s). Due to the poor knowledge on the source function, the simplified form of the above kinetic equation may just be useful for an appropriate description of the distribution function, allowing to study the role of the solar wind asymmetries on the spatial distribution of the zodiacal dust. In equation (58) Nr denotes the total number of dust particles in the zodiacal cloud, whereasfand g are both normalized functions. N” is the total production rate (i.e. total number of particles per unit of time), f”(a, E, i) and g”(s) are the normalized yield functions per unit volume of element space and the initial size distribution, respectively. ~,(a. E, i, s) is the collisional lifetime of particles with corresponding orbital elements a, E, i and size s. The drift velocity vector of dust particles in element space is denoted by u and is equal to

step

where the above time derivatives of the elements are given by equations (53)-(55). Introducing these drift vector components and calculating the divergence operator applied to it, equation (59) is then obtained in the form N 0.00’

L

0

’

’

L

100

Dynamical

(b)

c

J

g

200

time

I.

g

I’

t

’

af( a, E, i) &Z,(E ;i, 0)

T

(

au

a(1 -&2)3/2

step

a&.(& w)

+ af(a, 8, i)

a&

a'( 1 -&2)“2

’

=NT

u2(1 -&2)3’2 x f(u,

‘.:

0 1

44.0

(c)

0

1

0

1

’

’

h

100

Dynamical

200

time

1

2

E

step

E,

i) - N”f”(u,

E,

g”(s) + &fh

i)~

g(s)

pi)

~,(a, 8, is) .

(61)

This partial differential equation can be solved by integration along the mathematical characteristics of this equation which are not identical with dynamical trajectories of dust particles (for more details see Banaszkiewicz et al., 1994). First from one of these characteristics the following relation is derived

Fig. 4. (a) The semimajor axis a(t). For the initial values see Fig. 2 ; (b) the eccentricity E(t). For the initial values see Fig. 2 ; (c) the longitude of perihelion w(t). For the initial values see Fig. 2

c= where the quantities

malized to 1 AU (we omit the bars in the following text). The right hand side of the above equation represents source and sink terms for dust particles of size s at the element space point {a, E, i} One should clearly remember that this above differential equation is very much different from the FokkerPlanck type equation used by Hassan and Wallis (1983) or Pellat et al. (1984) as approximation of the Boltzmann integrodifferential equation describing the phase-space distribution function. The source term in general had to be represented as a collision integral dependent on the functionf(a, E, i) itself

(fL(&;L~)-HH,(i,~))

EP’

a( 1 - &‘)P?

(62)

p, and p2 are defined by PI =

Ha, + Ha2 W

p2=$_

(63)

(64)

E

Thus, choosing now E as the ing to equation (61), one inhomogeneous linear f(a) = ~(u(E), E, i) from the teristic of equation (60)

independent variable accordthen obtains the following equation differential for other mathematical charac-

H. J. Fahr et

df -= ds

: Evolution

a’(1 -E’)“’ aH,(i, ~0)s

zodiacal dust

asymmetric forces

~[H,(s;~,~)-HH,(~,~)l a2(1-&2)3’2

+ L f r,

L?(l -&*)I’* Wg”(s) pp. - -aH,(i, co).2 &g(s)

(65)

Hereby the following functions were used A(E) =

a2(1-&*)“2 aH,(i, ~)s

a[~i,(E;i,O)-HH,(i,o)l a*(1 -&2)3’2

+ L r,

(66)

a*(l--&2) 1’2N”g”(s)

B(E) = - -__

Pp. NH, t.i,w) s NT g(s)

(67)

The general solution of equation (64) then is obtained by f(s) = exp (% A(&‘)ds’)fW

+~~B(E’)exp(~A(E^)dE”)dEl.

and is very helpful because it substantially simplifies the exponential factor in equation (67) and permits to find the solution forf(s) from equation (68) in the form

(1

_,$y

&I

-PI

II,,

I

,

,,

0.50

0.10

t.00

Eccentricity

function f(a = 1 AU, E,i = 0) in arbitrary units vs the eccentricity for four different values of ratio rPR/rco (T,,&, = 0 solid, z,,/z,~ = 0.3 dashed, T,,&, = 1 dotted, rpR/rco= 3 dashed-dotted). The source function is given by equation (75) Fig. 5. Distribution

(73) Equation (69) can now be averaged over CO.A quick inspection shows that all o-dependent functions like the terms pl, p2, p3, p4 contain the argument (aO+ a, cos2w)/(bo+b, cos* co) which does not enable an analytic integration of equation (69). Instead we have calculated the numerical average according to (74) In the following we shall thus omit the above brackets. Since the coefficients A2 and A4 in an asymmetric solar wind do not vanish, we must expect the distribution function to be dependent on i via the functions H, = H,(i), Ha2 = H,,(i) and pj = pj(Ha2, H,), Jo { 1,2,3,4). Thus it is obvious that the resulting distribution function J(a, E,i) cannot be factored into two functions, one depending on i and the other independent on i. In Figs 5 and 6 we present calculations of the distribution function f(a, E,i) on the basis of formula (68)

P4

(-;)[1-E’] dE’ ;s’

I

0.05

(68)

In this solution for-f(s), the evolution of the distribution function is described along the drift trajectory of the dust grain in element space from sOto E. It is now profitable to use the expression for the collisional lifetimes ~~developed by Leinert et al. (1983) in which collisions of dust particles with identical main axis values considered and no collision-induced changes of the inclination i are taken into account. This expression is also motivated by the work of Dohnanyi (1978) or the one by Grtin et al. (1985) and has the form

x

I

IO' 0.01

1 _p

C70)

where we have introduced Ha, --He + TPRr, ’ P3 =He

’

(71) 0.9

0

I

20

I

40

I

50

I

80

I

Inclination

(72) and

Fig. 6. Ratio of distributions

f(a = 1 AU, E = 0.5, i)if(a = 1 AU, E = 0.5, i = 0) as a function of the inclination. The source function is the same as in Fig. 5 and given by equation (75)

310

H. J. Fahr et al. : Evolution of zodiacal dust under asymmetric forces

for different ratios of Poynting-Robertson lifetimes zpR over collision times zCO covering the range from very weak to fairly strong collisional losses. In all cases a source function was assumed which for the above-mentioned Figs 5 and 6 is adopted to be f” (a, E) = f{(a/a,)-

’ 3 .zexp (- 2~~)

(75)

describing a dust production due to injection rates decreasing monotonically with semimajor axis a, with a, = 1 AU, with an explicit dependence on E peaking at eccentricities of about E,,, = 0.7. Such an injection would be associated with dust evaporation from cometary comae close to cometary perihelion passages, assuming that the occurrence frequency of comet perihelia with values of about a,, and the associated cometary dust production as a function of the perihelion distance vary in this manner. In Fig. 5 we have plotted f=f(a = a, ; E; i = 0’) for various ratios of zp&,. For vanishing collisional lossesf strongly increases towards small eccentricities, whereas for more effective collisional losses, due to their s-dependence, this increase is changed into a decrease at the value of maximum feed-in. In Fig. 6 the ratio f(a = a, ; E = 0.5 ; i)/ f(a = a, ; E = 0.5 ; i = 0’) is shown as a function of inclination i for various effectivities of collisional losses. It turns out that the f(u, E,i) has a stronger increase with i the less effective are the collisional losses. This manifests the effect of the solar wind momentum flow asymmetry causing higher lifetimes at higher values of i. In Fig. 7 we show results similar to Fig. 5, however, this time calculated for an alternative form of the source function f” given by p(E,u)

=exp(-5(:-lr)exp(-5s*)exp(-tg’i) (76)

with a, = 3.5 AU. This source function describes a dust particle production associated with the asteroidal belt, peaking at a semimajor axis of‘a = 3.5 AU with a mean

0

10

20

30

40

50

60

70

80

90

inclination

Fig. 7. The same as Fig. 6, however, with the source function given by equation (76). The full curves are for rPR/z,,, = 0 ; the two curves are caldotted curves for TV&, = 3.0. The upper culated for an i-independent dust injection, the lower two curves for an i-dependent injection according to the function f”(i) = exp(-tg*i)

spread corresponding to Au = 3.5/,/5 AU and eccentricities and inclinations strongly peaked at E = 0 ; i = 0”. Figure 7 shows that the f monotonically drops off with increasing values of s for all values of the ratio r.p&,0. Here we do not aim at a calculation of the spatial dust density which can be obtained from the distribution by an integration over the 3-dimensional element space. We only want to mention here that this density n(r,9) is a function of configuration space coordinates, i.e. of heliocentric distance Yand ecliptic latitude 9, and is related to the distribution function f (a, E, i) by an integral equation given in the following form (see Haug, 1958)

(77) Results obtained from equation (76) by integration with the distribution function obtained in equation (70) are presented elsewhere (Banaszkiewicz et al., 1994).

Outlook and conclusions

We have considered effects on orbit evolutions of zodiacal dust particles due to the plasma Poynting-Robertson drag force induced by the latitudinally asymmetric solar wind momentum flow. A 40% decrease of this momentum flow from the ecliptic to the poles is manifest in observations causing about 10% increases of corresponding PoyntingRobertson lifetimes of zodiacal dust moving at high inclinations i. Particles moving in low-inclination orbits are drifting towards the sun correspondingly faster compared to those in high-inclination orbits. In line with the higher radial drift the evolving density profile with radial distance I in the ecliptic also turns out to be steeper than at higher latitudes. One thus finds a density excess at higher latitudes compared to lower ones if no preference of particle productions with respect to inclination i existed, i.e. if a spherically symmetric production of dust particles would prevail. The dust particle production can, however, be expected to be concentrated to lower inclinations, i.e. the source functionf” will definitely depend on i. Since we have proven systematic and stochastic Lorentz drifts which would change the inclination to be ineffective, and since in our treatment of the plasma drag only monoMaxwellian ion distributions were considered not leading to normal components of plasma drag forces there hence is no coupling between dust particles which orbit in different planes ip Therefore, the Poynting-Robertson evolution of dust in an orbital plane i, is separated from that in all other planes with ik # i,. It can, therefore, be concluded that the relative contribution of particles at high inclinations i due to their enhanced Poynting-Robertson lifetimes is increased relative to what could simply be derived from the i-dependent source function. In addition this idependence in the lifetimes has the effect that the resulting

H. J. Fahr er al. : Evolution of zodiacal dust under asymmetric forces

dust density distribution cannot anymore be factorized by setting n(r, 9) = n,(r) ~(9) .as was done by many authors following Haug (1958) and Leinert et al. (1983). We have neglected longitudinal variations in the plasma drag forces caused by imprints of the solar wind stream structures (fast and slow solar wind) and of the solar wind structures corotating with the solar corona and passing over fixed space points with.in a typical period of 27 days. For particles with c1b 1 AIJ and orbital periods of z > 1 year these longitudinal variations will be smeared out, and only longitudinally averaged plasma Poynting-Robertson drag forces need to be taken into account. For dust particles, however, with a < 0.176 AU and thus z < 27 days no smearing-out of the plasma drag will take place, since these particles are about corotating with the corona and staying nearly in phase with corotating structures. Thus they may be subject to plasma drag forces either connected with slow or wit!h fast solar wind plasma. This interesting effect will be studied by us in a later paper. Moreover, these particles are close to the 3 : 1 resonance with Mercury and their density may be decreased by this resonance. We have also mentioned that the plasma drag coefficient cD developed in this paper is asymptotically independent on the Mach number of the solar wind flow if Mach numbers grow large (i.e. MP 2: 10 ; supersonic flow!). For small Mach numbers, especially for subsonic solar wind flows with IV,, < 1.0, this drag coefficient is, however, dependent on this number and increases strongly with decreasing Mp. Thus it would be interesting to study in more detail the specific plasma drag effect on dust particles within the critical sonic point of the solar wind (Y< 20 Ye). The study of the zodiacal dust cloud very close to the Sun may thus be used as an independent and highly attractive indicator for the llocation of the sonic surface of the solar wind flow and on its 3-dimensional circumsolar shape. In addition in this region also the non-Maxwellian character of the distributicln function will strongly influence the results with respect to dust drift rates in element space, especially migrations in the inclination angle here are induced close to the Sun due to the existence of normal components of plasma drag forces. For further studies of the plasma drag effect on zodiacal dust one might speculate within this context on the specific evolutionary processes occurring in the protosolar dust accretion disk during the primordial state of the solar system when a strong solar T-Tauri wind, probably with strongly increased momentum flows outside of the ecliptic plane region, can be expected. Also a speculation on conditions for circumstellar dust around other young stars with bi-polar outflows would be interesting as a followup to these afore-mentioned calculations. Acknowledgements. This work was sponsored by financial support of the Deutsche Forschungsgemeinschaft both within the frame of the project II-C-9 r(Fa 97/9-4) (DFG-Schwerpunkt: Kleine Korper im Sonnensystem) and within the bi-national cooperation project : 436-POL-113/14/O between Germany and Poland. References Ajello, J. M., Solar minimum

Lyman tl sky background observations from Pioneer Venus orbiter ultraviolet spectrometer :

311 solar wind latitude variation, 14861, 1990.

J. Geophys. Res. 95, 14855-

M., Fahr, H. J. and Scherer, K., Evolution of dust particle orbits under the influence of solar wind outflow asymmetries and the formation of the zodiacal dust cloud,

Banaszkiewicz,

Icarus 107,358-374, 1994. Barge, P., Pellat, R. and Millet, J., Diffusion of Keplerian motion by a stochastic force, Astron. Astrophys. 109, 228232, 1982. Buitrago, J., Gomez, R. and Sanchez, F., The integral equation approach to the study of interplanetary dust, Planet. Space Sci. 31, 373-376, 1983. Burns, J. A., Lamy, P. L. and Sorer, S., Radiation forces on small particles in the Solar System, Icarus 40, 148, 1979. Conscdmagno, G., Lorentz scattering of interplanetary dust. Icarus 38,398410, 1979. Dohnanyi, J. S., Particle dynamics, in Cosmic Dust (Edited by

J. A. M. McDonnell), pp. 527-587. John Wiley & Sons, New York, 1978. Fahr, H. J. and Ripken, H. W., Dust-plasma gas interactions in heliosphere, in Properties and Interactions of Interplanetary Dust (Edited by R. H. Giese and P. Lamy), pp. 304324. Astrophys. and Space Sci. Library, D. Reidel, Dordrecht, 1985. Fahr, H. J., Ripken, H. W. and Lay, G., Plasma-dust interactions in the solar vicinity and their observational consequences, Astron. Astrophys. 102, 359-370, 1981. Giese, R. H., Optical investigation of dust in the solar system, in Solid Particles in the Solar System (Edited by I. Halliday and B. A. McIntosh), pp. 1-13. D. Reidel, Dordrecht, 1980. Griin, E., Zook, H., Fechtig, A. and Giese, R. H., Collisional balance of the meteoritic complex, Icarus 62, 244255, 1985. Hassan, M. H. A. and Wallis, M. K., Stochastic dust particle motions: Planet. Space Sci. 31, l-14, 1983. Haug, U., Uber die Haufigkeitsverteilung der Bahnelemente bei den interplanetaren Staubteilchen. Z. Astrophys. 44, 71-97, 1958.

Kneissel, B. and Mann, I., Spatial distribution and orbital properties of zodiacal dust, in Origin and Evolution of Interplanetary Dust (Edited by A. C. Levasseur-Regourd and H. Hasegawa), pp. 139-146. Kluwer Academic, 1991. Kojima, M. and Kakinuma, T., Solar cycle evolution of solar wind speed structure between 1973 and 1985 observed with the interplanetary scintillation method, J. Geophys. Res. 92, 7269-7279, 1987. Kumar, S. and Broadfoot, A. L., Evidence from Mariner 10 of solar wind flux depletion at high ecliptic latitudes, Astrophys. J. 228,302-311, 1979. Lallement, R., Bertaux, J. L. and Kurt, V. G., Solar wind decrease at high heliographic latitudes detected from Prognoz interplanetary Lyman Alpha mapping, J. Geophys. Res. 90, 1413-

1423, 1985. Leinert, C., Zodiacal light measure of the interplanetary environment, Space Sci. Rev. 18, 281-339, 1975. Leinert, C., Dynamics and spatial distribution of interplanetary dust, in Properties and interactions of Interplanetary Dust (Edited by R. H. Giese and P. Lamy), pp. 369-375. D. Reidel, North-Holland, 1985. Leinert, C., Link, H., Pitz, E. and Giese, R. H., Interpretation of a rocket photometry of the inner zodiacal dust light, Astron. Astrophys. 47,221-230, 1976. Leinert, C., Hanner, M. and Pitz, E., On the spatial distribution of interplanetary dust near 1 AU, Astron. Astrophys. 63,183-

187, 1978. Leinert, C., Richter, I., PItz, E. and Planck, B., The zodiacal light from 1.0 to 0.3 A.U. as observed by the Helios space probes, Astron. Astrophys. 103, 177-183, 1981. Leinert, C., R&r, S. and Buitrago, J., How to maintain the spatial distribution of interplanetary dust, Astron. Astrophys. 118,345-357,

1983.

H. J. Fahr et al. : Evolution of zodiacal dust under asymmetric forces

312 Levasseur-Regourd, A. C. and Dumont, R., Absolute photometry of zodiacal light, Astron. Astrophys. 84,277-219, 1980. Marsch, E., Mtihlhiiuser, K. H., Schwenn, R., Rosenbauer, H., Pilipp, W. and Neubauer, F. M., Solar wind protons : three-

dimensional velocity distributions and derived plasma parameters measured between 0.3 and 1 AU, J. Geophys. Res. 87,52-72, 1982. Mortill, G. and Grtin, E., The motion of charged dust particles in the interplanetary space--I. The zodiacal cloud, Planet. Space Sci. 27, 1269-1282, 1979. Morfill, G. E., Grtin, E. and Leinert, C., The interaction of solid particles with the interplanetary medium, in The Sun and the Heliosphere in Three Dimensions (Edited by R. G. Marsden),

pp. 455474. D. Reidel, Dordrecht, 1986. of spatial distribution of interplanetary dust grains by Lorentz forces, Astron. Astro-

Mukai, T. and Giese, R. T., Modification

phys. 131,355-363, 1984. Mukai, T. and Yamamoto, T., A model of the circumsolar cloud. Astron. Sot. Japan 31, 585-595, 1979. Munro, R. H. and Jackson, B. V., Physical properties of a polar coronal hole from 2 to 5 R, Astrophys. J. 213, 874886, 1977. Neutsch, W. and Scherer, K., Celestial Mechanics. Biblio-

graphische Institut, Mannheim,

1992.

Parker, E. N., Interplanetary Dynamical Processes. Interscience

Publishers, John Wiley & Sons, New York, 1963. Pellat, R., Barge, P., Hornung, P. and Millet, J., Adv. Space Res. 4,95-104, 1984. Pryor, W. R., Ajello, J. M., Barth, C. A., Hord, C. W., Stewart, A. I. F., McClintock, W. E., Sandel, B. R. and Shemansky, D. E., The Galileo and Pioneer Venus ultraviolet spectrometer

experiments: solar Lyman-cc latitude variation at solar maximum from interplanetary Lyman-a observations, Astrophys. J. 394,363-377, 1992. Reach, W. T., Zodiacal emission. II. Dust near ecliptic, Astrophys. J. 369,529-543, 1991. Scherer, K., The influence of the interplanetary magnetic field

on the inclination

drift of orbiting zodiacal dust particles,

Zcarus (submitted). Schwehm, G. and Rhode, M., J. Geophys. 42,721-735, 1977. Toller, G. N. and Weinberg, G. J., The change in the near-

ecliptic zodiacal light brightness with heliocentric distance, in Properties and Interactions of Interplanetary Dust (Edited by R. H. Giese and P. Lamy), pp. 21-26. Astrophysics and Space Science Library, D. Reidel, Dordrecht, 1985. Witt, N., Ajello, M. and Blnm, P. W., Solar wind latitudinal variations deduced from Mariner 10 interplanetary H (1216 A) observations, Astron. Astrophys. 73, 212-28 1, 1979.