Hydrogen Stark Broadening

Calculations With the Unified

Classical Path Theory

U.S.

•ARTMENT

OF

OMMERCE

National

Bureau

of

Standards
NATiONAL BUREAU OF STANDARDS

The National Bureau of Standards ' was established by an act of Congress March 3, 1901. Today,

in addition to serving as the Nation's central measurement laboratory, the Bureau is a principal

focal point in the Federal Government for assuring maximum application of the physical and

engineering sciences to the advancement of technology in industry and commerce. To this end

the Bureau conducts research and provides central national services in four broad program

areas. These are: (1) basic measurements and standards, (2) materials measurements and

standards, (3) technological measurements and standards, and (4) transfer of technology.

The Bureau comprises the Institute for Basic Standards, the Institute for Materials Research, the

Institute for Applied Technology, the Center for Radiation Research, the Center for Computer

Sciences and Technology, and the Oflfice for Information Programs.

THE INSTITUTE FOR BASIC STANDARDS provides the central basis within the United

States of a complete and consistent system of physical measurement; coordinates that system with measurement systems of other nations; and furnishes essential services leading to accurate and

uniform physical measurements throughout the Nation's scientific community, industry, and commerce. The Institute consists of an Office of Measurement Services and the following technical

divisions:Applied Mathematics—Electricity—Metrology—Mechanics—Heat—Atomic and Molecular Physics—Radio Physics -—Radio Engineering -'—Time and Frequency -—Astrophysics - —Cryogenics. -' THE INSTITUTE FOR MATERIALS RESEARCH conducts materials research leading to im- proved methods of measurement standards, and data on the properties of well-characterized

materials needed by industry, commerce, educational institutions, and Government; develops,

produces, and distributes standard reference materials; relates the physical and chemical properties of materials to their behavior and their interaction with their environments; and provides

advisory and research services to other Government agencies. The Institute consists of an Office of Standard Reference Materials and the following divisions: Analytical Chemistry—Polymers—Metallurgy—Inorganic Materials—Physical Chemistry.

THE INSTITUTE FOR APPLIED TECHNOLOGY provides technical services to promote

the use of available technology and to facilitate technological innovation in industry and Government; cooperates with public and private organizations in the development of technological

standards, and test methodologies; and provides advisory and research services for Federal, state, and local government agencies. The Institute consists of the following technical divisions and

offices: Engineering Standards—Weights and Measures — Invention and Innovation — Vehicle Systems Research—Product Evaluation—Building Research—Instrument Shops—Measurement Engineering—Electronic Technology—Technical Analysis.

THE CENTER FOR RADIATION RESEARCH engages in research, measurement, and application of radiation to the solution of Bureau mission problems and the problems of other agencies and institutions. The Center consists of the following divisions: Reactor Radiation—Linac Radiation—Nuclear Radiation—Applied Radiation. THE CENTER FOR COMPUTER SCIENCES AND TECHNOLOGY conducts research and

provides technical services designed to aid Government agencies in the selection, acquisition,

and effective use of automatic data processing equipment; and serves as the principal focus

for the development of Federal standards for automatic data processing equipment, techniques, and computer languages. The Center consists of the following offices and divisions: Information Processing Standards—Computer Information — Computer Services — Sys- tems Development—Information Processing Technology.

THE OFFICE FOR INFORMATION PROGRAMS promotes optimum dissemination and

accessibility of scientific information generated within NBS and other agencies of the Federal Government; promotes the development of the National Standard Reference Data System and a system of information analysis centers dealing with the broader aspects of the National Measure- ment System, and provides appropriate services to ensure that the NBS stafT has optimum accessibility to the scientific information of the world. The Office consists of the following

organizational units:

Oflrice of Standard Reference Data—Clearinghouse for Federal Scientific and Technical Information —Office of Technical Information and Publications—Library—Office of Public Information—Office of International Relations.

' Heaciquarlers and Laboiatories at Gait hers burg, Maryland, unless otherwise noted; mailing address Washington, D.C. 20234,

= Located at Boulder, Colorado 80302. Located at 5285 Port Royal Road. Springfield, Virginia 22151,
UNITED STATES DEPARTMENT OF COMMERCE

Maurice H. Stans, Secretary

National Bureau of Standards Lewis M. Branscomb

Hydrogen Stark Broadening Calculations With

the Unified Classical Path Theory

C. R. Vidal, J. Cooper, and E. W. Smith

Institute for Basic Standards

National Bureau of Standards

Boulder, Colorado 80302

^.S National Bureau of Standards Monograph 116 .

Nat. Bur. Stand. (U.S.) Monogr. 116, 143 pages (May 1970)

CODEN: NBSMA

Issued May 1970

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402

(Order by SD Catalog No. C 13.44:116). Price S1.50.
NAL'ONAL BUREAU OF STANDARDS

^ ft

^ '

Library of Congress Catalog Card Number: 75-606028
TABLE OF CONTENTS

I. Introduction 2

II. Basic Relations 5

III. The Tetradic Notation 9

IV. The Generalized Unified Theory 13

V. The No-Quenching Approximation for Hydrogen 19

VI. The Thermal Average S^'^^t) For Hydrogen 23

VII. The Multipole Expansion of the Classical

Interaction Potential 2.6

VIII. The Spherical Average of the Time Development

Operator ^^(t) 2«

-(1

)

IX, Evaluation of the Thermal Average 3" (t) for

Hydrogen 34

X. The Fourier Transform of the Thermal Average 46

XI. The One-Electron Limit for Hydrogen and the

Asymptotic Wing Expansion 57

XII. The Unified Theory for Hydrogen 62

Appendix A 73

Appendix B 92

Appendix C 103

References 119

Figures 122

iii

HYDROGEN STARK BROADENING CALCULATIONS WITH

THE UNIFIED CLASSICAL PATH THEORY

C. R. Vidal, J. Cooper, and E. W, Smith

The unified theory has been generalized for the

case of upper and lower state interaction by introducing

a more compact tetradic notation. The general result

is then applied to the Stark broadening of hydrogen.

The thermal average of the time development operator

for upper and lower state interaction is presented. Except for the time ordering it contains the effect of finite

interaction time between the radiator and perturbers to

all orders, thus avoiding a Lewis type cutoff. A simple

technique for evaluating the Fourier transform of the

thermal average has been developed. The final calculations based on the unified theory and on the one-electron theory are compared with measurements in the high

and low electron density regime. The unified theory calculations cover the entire line profile from the line center

to the static wing and the simpler one-electron theory

calculations provide the line intensities only in the line

wings.

Key words: Classical path; hydrogen lines; line

wings; one-electron theory; Stark

broadening; unified theory.

* This research was supported in part by the Advanced

Research Projects Agency of the Department of

Defense, monitored by Army Research OfficeDurham under Contract No. DA-31-124- ARO-D-139.

1
I. INTRODUCTION

For the first few Balmer lines of hydrogen, recent papers

(Cerardoand Hill, 1966; Bacon and Edwards, 1968; Kepple and Griem,

1968; Birkeland, Oss and Braun, 1969) have demonstrated fairly good

agreement between measurements in high electron density plasmas

16-3

(n > 10 cm ) and improved calculations of the so called "modified

e impact theory". The experimental and theoretical half-widths differ

less than about 10%. However, measurements of the Lyman-a wings

(Boldt and W. Cooper, 1964; Elton and Griem, 1964) and low electron

13 -3

density measurements (n =10 cm ) of the higher Balmer and

e Paschen lines (Ferguson and Schliiter, 1963; Vidal, 1964; Vidal, 1965)

have revealed parts of the hydrogen line profile, for which the modified impact theory appears to break down. For the higher series

members better agreement has been obtained with quasi- static calculations (Vidal, 1965). The reason the current impact theories

break down is that these theories correct the completed collision

assumption by means of the Lewis cutoff (Lewis, 1961) which is only

correct to second order. With this cutoff it was possible to extend

the range of validity for the impact theory beyond the plasma frequency.

However, in the distant wings, where the electron broadening becomes

quasistatic, the second order perturbation treatment with the Lewis

cutoff breaks down because the time development operator must then

be evaluated to all orders. Attempts to correct the second order theory

have been made already (Griem, 1965; Shen and J. Cooper, 1969), but

these theories still make the completed collision assumption by replacing the time development operator by the corresponding S-matrix,

and so it has to be emphasized that in conjunction with the Lewis cutoff

these theories would only be correct to second order. The impact

theory in its present form is intrinsicly not able to describe the static

wing and the transition region to the line center where dynamic effects

2
' cannot be neglected. To overcome this problem, several semiempirical procedures (Griem, 1962; Griem, 1967a; F. Edmonds, Schluter

and Wells, 1967) have been suggested to generate a smooth transition

from the modified impact theory to the static wing.

been reinvestigated in two review papers ( E. Smith, Vidal and J. Cooper, 1969a, 1969b), which are from now on referred to as papers

I and II. The purpose of I and II was to state clearly the different

approximations which are required to obtain the classical path theories

of line broadening and to find out where these theories are susceptible

to improvements. In a manner similar to the Mozer-Baranger treatment of electric microfield distribution functions (Baranger and Mozer,

1959, I960), it was shown that the general thermal average can be

expanded in two ways, one of which leads to the familiar impact tieory

describing the line center (Baranger, 1958, 1962; Griem, Kolb and

Shen, 1959, 1962). The other expansion represents a generalized

version of the one electron theory (J. Cooper, 1966), which holds in

the line wings. It is also shown that there is generally a considerable

domain of overlap between the modified impact theory and the one

electron theory. Based on these results, a "unified theory" was then

developed (E. Smith, J. Cooper and Vidal, 1969), henceforth referred

to as paper III, which presents the first line shape expression which is

valid from the line center out to the static line wing including the

problematic transition region. The line shape obtained by the unified

theory has the form

Recently the classical path methods in line broadening have

(I. 1)

3
where d, Auu and jb(Aaj) are operators. In paper III it was shown that

the familiar impact theories, which hold in the line center, may be

obtained by making a Markoff approximation in the unified theory,

while the one electron theory describing the line wings is just a wing

expansion of the unified theory. Consequently the crucial problem for

any line broadening calculation is to evaluate the matrix elements of

jL(AtD), which is essentially the Fourier transform of the thermal average

(see Eq. (46) and (47) of paper III). This will be done in detail in this

paper for the general case of upper and lower state interactions.

In the following Sec. II we start with a brief summary of the

basic relations which are required for the classical path approach

pursued here. We then generalize the results of the unified theory to

include lower state interaction (Sec. IV) after introducing a more

compact tetradic notation (Sec. III). From this general result we turn

to the specific problem of hydrogen by discussing briefly the no

quenching assumption (Sec, V) and deriving the thermal average 5^^^(t)

(see Eq. (47) of paper III) for the general case of upper and lower state

interaction (Sec. VI). We next investigate the multipole expansion of

the classical interaction potential in the time development operator

(1

)

(Sec. VII). The thermal average 5 (t) is then evaluated in two steps

by first performing a spherical average (Sec. VIII) and then an average

over the collision parameters: some reference time t , impact parao meter p and velocity v (Sec. IX). Appendix A gives the computer

program which we used in calculating the thermal average for dipole

interactions including lower state interactions. The large time limit

of the thermal average, which leads to the familiar impact theories in

the line center, is investigated in detail in Appendix B for different

cutoff procedures and compared with the results in the literature. In

4
Sec. X, a method is developed for performing the Fourier transformation of the thermal average and it leads us to the crucial function for

any classical path theory of Stark broadening. This function is finally

applied in Sec. XI to the one electron theory, which forms the basis

for the asymptotic wing expansion, and in Sec. XII to the unified theory,

which describes the whole line profile from the line center to the static

wing. Numerical results are given for the hydrogen line profiles as measured by Boldt and W. Cooper, 1964; Elton and Griem, 1964, and

j Vidal, 1964, 1965. The computer program for the unified theory calculations and the asymptatic wing expansion is given and explained in

Appendix C.

11. BASIC RELATIONS

In this section we will briefly outline the basic relations which

^ are used in our classical path treatment of line broadening.

As discussed in Sec, 2 of paper I, we are considering a system

containing a single radiator and a gas of electrons and ions. We will

make the usual quasi- static approximation for the ions by regarding

their electric field 6. as being constant during the times of interest

1/Au) . This approximation is usually very good because the region

where ion dynamics are important is normally well inside the half

width of the line except for a few cases such as the n-a lines of

hydrogen (Griem, 1967b). The complete line profile I(uu ) is then given

by the microfield average (see Eq. (3) of paper II)

00 i((w)= /'p(e.) I ((ju,e.) de. (II. 1)

where the normalized distribution function P{F..) is the low frequency

component of the fluctuating electric microfields. Due to shielding effects

Pit.) depends on the shielding parameter r /D where r and D are the

1 o o mean particle distance and the Debye length (for electrons CTnl^jrespectively.

5
With the static ion approximation we have reduced the problem

to a calculation of the electron broadening of a radiator in a static

electric fields . The resulting line profile I(aj,e.) is then simply

1 1 averaged over all possible ion fields to give the complete line profile

I((jd). The static ion field will be used to define the z-axis for the

radiator and the ion- radiator interaction will be taken to be the dipole

interaction eZP, . where -eZ denotes the Z-component of the radiators

1 dipole moment.

If the unperturbed radiator is described by a Hamiltonian H , a we may then define a Hamiltonian for a radiator in the static field

P.. by

H = H + eZp,. {II. 2) o a 1 The complete Hamiltonian for the system is then given by

H = H + V (R, X, V, t) (II. 3) o e

where denotes the electron radiator interaction. In this equation,

X and V are 3N vectors x = (x^ , 5c^, • . .

, x^), v = (v^,^^' ' * *

' ^n^'

which denote the positions and velocities of the N electrons and R

denotes some internal radiator coordinates. For one-electron atoms,

R is the position of the "orbital" electron relative to the nucleus.

The interaction V will be regarded as a sum of binary interactions,

V (R, X, V, t)= y V.(R, X., v., t) (II. 4)

J where denotes the interaction between the radiator and a single

electron. As is well known the line shape I(aJ,£.) may be given by

the Fourier transform of an autocorrelation function C(t) (Baranger, 1962)

fe'^^C(t)dt (II. 5) o
In the classical path approximation, the correlation function for electric

dipole radiation is given by

C(t)=Tr I d (t"'" (t) d T (t) > p \ , (II. 6) a I a a av a /

where and p denote the dipole moment and the density matrix for

the radiator. The thermal average denoted by the subscript av

represents the average over electron states (see Eq. (47) of paper I):

<T"'"(t)dT (t)) =/dxdvP (x)W(v) t'^(R,x, V, t) d T (R, x, v, t) (II. 7) a a av I a a

where P(x) and W(v) are the position and velocity distribution functions

for the electron perturbers (defined by Eqs.(37) to (40) in paper II).

The time development operator for the system T (R, x, v, t) is the

solution of the differential equation

ih^T (t)= [h +V (t)l T (t) (11. 8) 9t a

L ° e I a

and it may be written in an interaction representation defined by

T (R, X, V, t) = exp (-itH /h) U (R, x, v, t) (II. 9) a o. a

\ where

fti-^ u (t) = V (t) u (t) (11.10) 3t a e a

7
and

V (t) = exp (itH /h)V (t) exp (-itH /h) . (II. 11)

e o e o

It should be noted that V (t) is identical with V (t) in paper II except

e e that we have not yet made the no quenching assumption which removes

the unperturbed part H in the Hamiltonian H in Eq. (II. 11), Using

a o the time ordering operator O, U (t) may be written in the form

U (R, X, 'v, t) = exp f

- ^ /'v (R, X, V, t') dt' !. (11.12) a

To evaluate the trace over atomic states in Eq. (II. 6), it is convenient

to use the H eigenstates I a o Hence, using U (t) we haveto use the H eigenstates la), lb), • • • with the eigenvalues E , E , o ° a b

-iju , t dc

C(t)= ^<a!d|b) (crd|d>e (11.13)

abed

(b|u"^(t)lc) <dl U (t)|a>

2i Si <a|p la)

av a

where

0) , = (E,-E )/h (II. 14) dc d c

In paper II and III, the correlation function C(t) was evaluated

for the case of no lower state interactions in order to keep the mathematics as simple as possible because one of the U (t) operators in

3. Eq. (II. 13) may then be replaced by a unit operator. In this paper we

will give a more general evaluation of C(t) which includes lower state

interactions. For this purpose we introduce in the next section a more

compact tetradic notation. Furthermore, it should be noted already

at this stage that we will interchange the sequence of approximations

8
with respect to paper II by deriving the generalized unified theory

before making the no quenching approximation. This makes the results

of the unified theory also useful for situations where the no quenching

approximation cannot be made like, for example, microwave lines.

III. THE TETRADIC NOTATION

The purpose of the tetradic notation which we shall use is to

write the product of the U (t) operators in Eq. (II. 13) in terms of a

single operator. To do this we first consider the product of the matrix

elements (a I aI a' ) and ( B I B! g ' ) where A and B may be any arbitrary

operator. This product may be written in terms of the direct product

A|^B according to

(alAla') <e!B|e'> = (a^\ABB\a'B') , (IILl)

where the product states laB) = la)lB) are essentially the same as the

states of Barangers "doubled atom" (Baranger , 1962). This direct

product, Ag)B, is a simple form of tetradic operator. K one of the

operators A or B happens to be a unit operator I, we may conveniently

denote this fact by means of superscripts t and r according to

<aSl A2;i|a'B')= <aB ! A'^ la 'B ') = (a\A\a') (III. 2)

Bp

<aBll^-Bla'p') = (aSlB^'la'B') = <b!bIb') 5 , . (III. 3)

act

That is, a superscript ^ denotes a "left" operator which operates only

on the "left" subspace (in this case the la), la') subspace) and a superscript r denotes a "right" operator which operates on the "right"

9
subspace. It is thus clear that any "left" operator will commute with

any "right" operator:

[A^,B^] = 0. (III. 4)

With this notation, the thermal average in Eq. (II. 13) can now be

written in the more compact form

<blu"^{t)lc> (dIU (t)|a) a a

(cd| TJ^ "\t) U'^jba)

av

(III. 5)

(clu''(t)|b> (dIU (t)la>

av

av

= (cd| U (t) (t) a a av I ba> .

con- We have chosen to write (blu'^(t)Ic) as < c I U ='T(t)

j b) simply for

venience in the derivation given in later sections. Noting the definition

of U^(t) given in Eq. (II. 12), we define operators (R,x, v, t) and

V (R, X, v", t) so that

(R, ^, V, t) = exp

I

-

t

. , (R,5^, ^,t') dt'

(III. 6)

U (R, 5c, V, t) = O exp ) / V "^(R, ^, ^, t

•t

') dt'j

Since any "left" operator commutes with any "right" operator,

have

we

U^;^"~(R, t V, t) U''(R, X, V, t) = <^ exp

a a

=^(R,x,i;,t) (in. 7)

10
where

r{R,t,v,t)= V"" (R,x, V, t) - V^''(R,x, V, t). (III. 8) e e e

We have now succeeded in replacing the two U (t) operators by a more

general tetradic operator "^(t) which operates in both "left" and "right"

subspaces. Eq. (III. 5) thus becomes

r<blU^t)lc) <d|U (t)la>l = <cdim(t)l Iba).

(III. 9)

It is important to realize that the tetradic operator ^(t) is formally

the same as the operator U (t); that is, it satisfies the same type of

3. differential equation

ih^^(R,5,v, t) = r^ (R,x, v,t)^(R,x,v,t) . (III. 10)

This means that all of the line broadening formalism which has been

developed for U (t), will be directly applicable to ^^(t).

To make the formal correspondence more complete we use the

operators , h"", V (R, X, V, t) and V (R, x, v, t) to define the tetradics

o o e e K and y (R, x, v, t) according to

o e

-K: = - H^^ (III. 11)

o o o

r (R,x, V, t) = v'"(R,x, V, t) - V^^'^ (R,x, V, t) , (III.12)

e e e

11
Since any left operator commutes with any right operator we have

(t) = exp ; itHVh ^ V^(t) exp |-itH^/h> (III. 13)

^

o

I e I o

I

I

it^cyh

I

v^{t) exp

I

= exp -i«ryh|

Hence

(III. 14)

rjR, X, V, t) = exp

I

i^^qA

j "^J^' ^) |"^^o^^|

which is formally the same as Eq. (II. 11), It is also obvious that both

V and V will be given by a sum over binary interactions \ or 7/" e e 11

just as in Eq. (II. 4).

1r (R, x, V, t) = (R, X ., v., t) III. 15

e 4-^ 1 J J

J

r.(R, x.,v.,t)= vJ;(R,x., v.,t) - v'[^''(R,x., v.,t) (in. 16)

The formal similarity between the operators H , V (t), V (t), U (t), o e e a

etc. and the tetradic operators K^, ^g^*)' ^\(^)> will greatly

simplify the treatment of the thermal average foi* the general case of

upper and lower state interactions.

12
IV. THE GENERALIZED UNIFIED THEORY

Using the tetradic operators as defined in the previous section

jj we have for the correlation function

C(t)= 2^ (aldlb) <cld!d> e """"(alp^la)

abed

<cd|"F(t) I ba) (IV. 1)

where 5 (t) denotes the thermal average of^*|^(R, x, v, t):

= yd^dv P(x)W(v)^(R, X, V, t) . (IV. 2)

This tetradic operator 3 (t) is formally identical to the operator F(t)

defined in Sec. (2. A) of paper III. It would also be formally identical

y to the F(t) defined by Eq. (19) of paper II if we would make the no

I quenching approximation at this point. To preserve generality, however,

the no quenching approximation will be deferred until a later section

when we specify the la), lb), * * * eigenstates to be H^ eigenstates for

hydrogen.

13
Following the formalism developed in Sec. 2 of paper III, we

define an operator ^{H, x, v, t) by

g;(R,x,v:t)= P(x) W(v)^(R,x, ^?,t) (IV. 3)

so that

:3{t) = / d5"dv 3;(R, x, v, t) (IV. 4)

(cf. Eqs. (11 ) and (12) of paper III). From Eq. (III. 1 0) we see that

-

„

. ih

^

5 (R, % v,t) = r X, V, t) 3?(R, X, ^, t) (IV. 5)

o t e

which is formally the same as Eq. (13) in paper III. We next introduce

a projection operator P which is identically the same as the operator P

defined by Eq. (14) of paper III (the fact that P now operates on tetradics

does not change its definition). That is, for any function of electron

variables f(x, v) we have

Pf(x, v) = P(x) W(^/)ydx'dv' f (x', v') , (IV. 6)

This relation holds whether f is a matrix, tetradic or any other type

of operator. With this operator we can follow the derivation in Sec. (2. B)

of paper III replacing H , V , V etc. hyK ,">/ ,'if , etc. As a result

o e e o e e

(cf. Eq. (27) in paper III) we have

14
r^r 3(t) = - h"^/ exp (itV A) y (t-t') Q(t-t')r (0) ht X o L

e e J av

exp {-it'K /h) 3(t') dt' (IV. 7) o where

•t-t' Q(R,x,^, t-t') = O exp |-^^(l-P)r^{R,x, V, t") dt"| . {IV. 8)

Returning to Eqs. (II. 5) and (IV. 1) we see that the quantity of interest

is not 3(t) but rather its Fourier transform.

/•« . -iju t (cdlJ(aj)lba) =/ e"^ e ( cdl W(t)l ba> dt

Jo

e ^^(cdl exp (-ittr /h)'^(t)Iba) dt

=J

e"^^ <cdl 3^t)Iba) dt (IV. 9)

where

g;(t) = exp (-iuc /h) 3(t) . (IV. 10)

o

From Eq. (IV. 7) we see that

-|-?(t) = -(i^- /h)'?(t)-h"7 exp (-i(t-t')ar /h)

ht o' o' I exp (-i(

•'o

(t-f) Q(t-t')r (0)1 ^(t') dt'. (IV. 11)

j e Q I av

Solving this equation by Fourier transforms gives

-9(uu) = i

I Auu -£(Auu )\ (IV. 12)

op op J

15
where

X{AUJ ) = - ih

op

exp (itAuu ) r (t)Q(t)r (0)

op av

dt (IV. 13)

and Auu is an operator defined byop

tm = ud-k: /h = aj-(H^-H^-' )/h . op o o o

(IV. 14)

With these results, the line shape given in Eq. (II. 5) becomes

(cfEq. (I. D)

I((JU,e.) = — Im (a!dlb) <c]dld) (ajp la)

(cdl A«3 -Z^m )1 ^

I ba) . [op op

J

' ' (IV. 15)

We next simplify jb(Auu^p) by means of the impact approximation

(see Sec. (3.2) of paper II). Basically this approximation assumes

that the average collision is weak, that strong collisions do not overlap

in time and that a weak collision overlapping a strong one is negligible

in comparison (weak collisions are those interactions for which a low

order perturbation expansion in T/" provides a good approximation to or Q ; for strong collisions the full exponential must be retained).

It should be emphasized again that we make a distinction between the

impact approximation and the impact the'ory . The latter contains the

impact approximation as well as other approximations like the completed collision assumption which will not be made here. We also

assume that the electron perturbers may be replaced by statistically

independent quasi particles (e. g. , shielded electrons). In Sec. (3)

and Appendix B of paper III, it is shown that these approximations

16
reduce ^(Atu ) to

op

•t itAuD e °^ 3 (t) dt AO) {IV. 16)

op

u where

^\t)= d5?^dx^^W(v^) j^<^^(R,5rJ,^?^,t) - ij (IV. 17)

and

^^(R, x^, v^,t) = (yexp Z^j(R, x^, v^,t')dt'|(IV. 18)

and n denotes the electron density,

e Equations {IV. 15) through {IV. 18) give the line profile of the

generalized unified theory. To obtain the impact theory we simply

replace i!{Auu ) by Z{0) and as discussed in Sec. 4 of paper III, we have

op

the familiar result (cf. Eq. (44) of Baranger, 196Z).

£(0) = iy[S^'" S^^ - 1) dv (IV. 19)

where S^^ denotes an S-matrix for a binary (completed) collision and

^

dv denotes the integral over collision variables, as defined in the

Appendix of paper II.

00 ydv - n^ ^ dv vfCv) dp 2rrpJ'dfl (IV. 20)

In comparing Eq. (IV. 19) with Barangers result it is important to note

that Barangers operators S. and operate only on "initial" and "final"

states respectively, whereas our operators and S^^ operate on all

17
possible H eigenstates. This difference occurs because we have not

o made the no quenching assumption yet.

The other limit of the one electron theory is obtained by making

a wing expansion of the unified theory; that is, the operator

[Am -I (Acu )]'

op op

order this gives

[A(« (Acu )1

^ is expanded in powers of [jb(Atw )/A(A) 1. To lowest

op op op op"

[m -ziisw )'\~^ = + zim )-r- + op op A.UJ Auu op A^

op op op

I rt ^^^^ n ) = ^ - i/ e V)dt+ ••• . (IV.21)

op •/o

The first term, ^/A^u^p, gives a delta function when one takes the

imaginary part required by Eq. (IV. 15). To get this delta function we

approximate radiation damping effects by using (a^d + ie) in place of

op

AOJ (see Sec. (3. A) of Smith and Hooper, 1967); the imaginary part of

op

1 /Muu + ie ) is just -nbiAiU ) when e->0. When this delta function term

op op

is averaged over ion fields according to Eq. (II. 1) it will produce the

line broadening due to the static ions alone (see Sec. 5 of paper II).

The influence of the electrons as well as electron-ion coupling is

contained in the second term of Eq. (IV.21). Hence one is interested

in the matrix elements of the Fourier transform of 3^^^(t), which is

also the quantity of interest in the unified theory (see Eq. (IV. 16)). The

primary difference between calculations made by the unified and oneelectron theories is therefore the matrix inversion of f Auu -X(Auu )1

op op

which is required by the unified theory but not by the one electron

-(1

)

theory. Since the matrix elements of the Fourier transform of 3" (t)

play such a central role in any classical path theory (including the impact theory), the evaluation of these matrix elements for hydrogen will

be discussed in detail in the following sections.

18
V. THE NO-QUENCHING APPROXIMATION FOR HYDROGEN

In the preceding section we have derived the thermal average

3f{t) and its Fourier transform .9 (uo ) for the general case of upper and

lower state interaction. In order to evaluate I((jo,S.) in Eq. (IV. 15)

we have to consider the complete trace over all H eigenstates

la), lb), ••• . However, in looking at the Eqs. (IV. 9), (IV. 18) and

(III. 14) one realizes that due to the exponential factors only a few of all

the possible matrix elements will contribute significantly to the final

line profile at a particular frequency uu. That is, we can neglect those

matrix elements for which the argument of the exponential factor is

so large that it gives rise to rapid oscillations within the range of the

time integral. Hence, if one treats well isolated lines, only those

matrix elements of Uj(t) between either "initial" or "final" states

have to be considered. We may therefore state the no-quenching

approximation as

^j(t)= V^^\t)\j\{t) (V.l)

where now no longer operates on the complete "left" or "right"

subspace, but only on "initial" or "final" states ( see also Sec. 2.2

and 7. 2 of paper II).

Further approximations cannot easily be generalized and depend on

the particular problem investigated. We now apply our general results

to the problem of hydrogen. In this case the no-quenching assumption

states that we need to consider only those matrix elements of U^^t)

and '\^^(R, ^jj'^j' ^) which are diagonal in the principal quantum number

n. As shown already in paper II this is a good approximation as long

as the lines investigated are well separated. For calculating the line

wings it is furthermore required that there is no appreciable overlap

19
with wings of adjacent lines in the region of interest. The same is

true also in any reliable measurement of line wings.

To show this we can proceed as in Sees. 2. 2 and 7. 2 of paper

II with the difference that now we are dealing with the operator

H = H + eZ£. rather than just H . Since H does not commute with

o a 1 a a Z we introduce a projection operator P (see Sec. 2. 2 of paper II)

n

which picks out the part of an operator which is diagonal in n. Using

this operator we split H into a part which is diagonal in n

H = H + eP Ze. (V. 2)

on a n 1

and a part which is not diagonal in n

H . = e{l-P ) Ze.. (V. 3)

off n 1

H now commutes with P Z because both operators are diagonal in a n

parabolic Coordinates. We therefore specify their eigenstates completely by the principal quantum number n, the magnetic quantum

number m and the quantum number q which is defined to be

q^n^-n^; (V. 4)

n and n are the usual parabolic quantum numbers which obey the

relation

n = n + n + I m| + 1 . (V. 5)

20
Knowing the solution of the eigenvalue problem

H nqm) = E I nqm) (V. 6) on' nqnn'

with

E = E + eZ p,. (V. 7) nqm n nq i

we see from a second order perturbation approach (cf. Chap. 1 6 of

Merzbacher, 1961) that the energy correction

2

(2)

^ j<nqm!H^^^ln'q'm')|

I ^^nqm=2.,^-^ -E . . . (^-^^

rFn* nqm nqm

can always be neglected as long as the ion fields do not become too

large. This is again equivalent to stating that the lines have to be well

separated.

As a result one is left with the eigenvalues E , E , a b

of the Hamiltonian H whose eigenstates la), lb), • • • I d) are the

on

parabolic states I nqm) . This allows us to rewrite the autocorrelation

function C(t) in Eq. (IV. 1) for hydrogen in the form:

C(t,r..) (nq^m^l dl n'q^m^ ) < n'q^m^l dl nq^mj^ )

exp f- — [e -E /+e{Z -2 , , )e.lt|(nqm !p Inqm) ^

y h L n n' nq^ ^^b ) ^ a a a a'

^n'q'bbm' ; nq^m^!?{t)!n'q'm' ; nq m ) (V. 9) bb aa aa

where quantum numbers which refer to the lower state are distinguished

from the upper state quantum numbers by a prime.

21
The matrix elements of P Z are given by (see Bethe and

n

Salpeter, 1957)

(nqmlzlnqm) = Z = — nqa (V. 10)

nq 2 o

2 2 with a^ = h /(me ) being the Bohr radius. As a further definition the

ion field 6. will be normalized to the Holtsmark field strength 6

1 ^ o

P.. = B-e (V. 11)

1 o

where

( 4tjY

2 l-T") e n^ . (V. 12)

This yields

-1-{Z -Z , = AUD.(n,q ,n',q') • R {V.13)

h

^ nq^ n'q^ lib b

with

/ - -

Aa).(n, q^, n', q^) = . -\ (nq^-n^q^) ^ n^ .

(V. 14)

Am. is now the frequency shift of a particular Stark component characterized by the quantum numbers ^> "^^j ^' ^^'^ to the Holtsmark

field strengths . Introducing the frequency shift t\\Xi = , where

nn

the frequency of the unperturbed line uu , is given by

nn

(D , = (E -E ,)/h (V. 15)

nn n n

22
the line profile I(Auu, p) can be written in the form

I(AU3,B) = — y <nq m I dl n'q'm') <n'q'mM d| nq m )

TT ^—^ a a a a bbbb

(V. 16)

(nq m Ip I nq m ) ('n'q'm' ; nq^m, | j(tju )| n'q'm' ; nq m ) a a a aa bb bb' aa aa

where

(n'q'm''; bbbbnq m L9(uu) In'q'm'; aaaanq m ) = ^^^^^^ X

dt exp < i(Acu -Auu .B )t > (n'q/m/ ; nq, I 5(t)l n'q'm' ; nq m ) . ^ i| bbbb aa aa

Performing the ion field average according to Eq. (II. 1) will then give

us the desired line profile once we know the thermal average ^(t).

VI. THE THERMAL AVERAGE ^^^V) FOR HYDROGEN

In Sec. ,IV. we saw that the crucial problem in any classical

path theory of line broadening is the evaluation of the matrix elements

of 3^^^(t). With the no-quenching approximation for hydrogen a typical

matrix element in parabolic states is given by

(n'q'm'; nq, m_ I 3?^^ ^(t)] n'q'm' ; nq m > = bb bb aa a a

^ ^

^/^^l^^l^^^l^ (n'q^m; ; nq^m^l^2/|^(t).l

] n'q^m^; nq^m^>

To simplify the evaluation we transform to the natural collision variables p, V and t which denote the impact parameter electron velocity

o ' and some reference time of the collision (see the appendix of paper II).

The orientation of the collision axes with respect to the radius vector

R of the orbital electron is specified by the three Euler angles

23
represented by Q. Furthermore we assume a spherically symmetric

distribution of perturbing electrons; this is a good approximation as long as the impact parameters are not too small. The velocity distribution function W(v) is related to the Maxwell distribution function

f{v) by

f{v) = 4ttv^W(v). (VI. 2)

With the preceding definitions Eq. {VI. 1) can be rewritten as

(n'q'm' ; nq m, j'S^^ V)l n'q'm'; nq m ) = b b b b a a a a

{VI. 3)

^Jdnjdv vf(vydp p^dt^ ^^'^^{^^^J\}^ithi\n'<^^m^^

Next we have to know the matrix elements of the time development

operator Ulj^{t) defined by Eq. (IV. 18). This requires the matrix

elements of the interaction potential ^^(t). In order to save some

writing we consider for the moment only Uj^(t) and Vj^(t) which after

making the no-quenching assumption may be the "initial" or "final"

part of the corresponding tetradic operators (see Eq. (V. 1)). A

typical matrix element of Vj^(t) is given by

<nq m I V, (t)l nq ,m ,) = expJi— (Z -Z )e.tWnq m j V, (t)j nq ,m , ^ ^c c' r ' ^d d'^ h nq nq / i I ^ ^c c' 1 ' ^d d

^ c d ) ^^j^^j

With the no-quenching assumption the unperturbed energy eigenvalues

E have cancelled. At this stage we now make another simplification

n

by dropping the exponential in the latter equation; this has been done

24
jin all previous Stark broadening calculations but it is rarely stated

explicitly. This will be a good approximation in the line wings where

the times of interest 1/auu are small and Auu is much larger than the

average ion field splitting. In the line center, however, the argument

of the exponential can easily be on the order of unity or larger in

which case V^(t) effectively vanishes due to rapid oscillations of the

exponentials. This effect was first noted by Van Regemorter, 1 964 , who

shows that this effectively introduces another cut off which may easily

be smaller than the usual Debye or Lewis cutoffs . This additional

cutoff has been included in recent calculations( Kepple and Griem, 1968).

However, as discussed in Sec. XII it turns out that its influence on the

final line profile is in most cases negligible.

Neglecting the ion field exponentials in Eq, (VI. 4). the time

development operator U^^ is now given by

where the time ordering is still required because P Vj^(t) need not

commute with P V^(t'). In paper II it was shown that this time ordering

is negligible for weak collisions (to second order) as well as quasistatic collisions (i. e. , in the distant line wings). Time ordering is not

negligible for strong collisions; however, when the thermal average is

performed, the errors due to neglecting time ordering are expected to

be small. The reason for this is that the time development operator

still retains its unitarity (cf. Sec. 8 of paper 11).

25
VII. THE MULTIPOLE EXPANSION OF THE CLASSICAL

INTERACTION POTENTIAL

Before evaluating the thermal average 5^^\t) we briefly consider

the classical interaction potential V^(t) due to a single electron. If the

perturber does not "penetrate" the radiator, Vj^{t) is given by the well

known multipole expansion

^[r(t)1

where I R| is the distance of the orbital electron fronn the nucleus, r(t)

00 ^ ]^ V,(t) = T; — • Pk [cos e(t)l (VII. 1)

is the instantaneous distance of the perturbing electron, the P are

Legendre polynomials and 0(t) is the instantaneous angle between R

and r(t).

In most cases it is sufficient to consider only the dipole (k=l)

term. However, to account for some asymmetries of a line, it may be

necessary to keep some of the higher multipole terms as well. In any

case, one can show that this multipole expansion is terminated after

some finite number of terms due to symmetries of the radiator.

-> i> To show this we specify the angular positions of R and r(t) by

9.,, cp., and 0 , cp respectively and we apply the spherical harmonic

addition theorem (Eq. (4.6.7) of Edmonds, I960)

+k

P^(cos0)=^ (-1)^ (0^, cpj) • 0^02, cp^) (vn.2)

p=-k ^

^

where

cos 0 = cos 0 COS 0 + sin 0 sin 0 cos (cp -cp_). (VII. 3)

X Ld X ui L iL

26
We may simplify the mathematics without loosing generality by choosing

a coordinate system in which cp = 0. Using the relation

(VIL 4)

one then obtains

p fcos e)= (e,,cp,)c^ (9,) + k o 1 1 o Z

p=i

which gives for the interaction potential

\tl^

k=l [r(t)]

Cok. P, (cos e-(t)) +

d

The dipole case {k=l) gives the well known result

r (t)

cos 9,(t)+^jc^j-C;}si

(VII. 5)

(VII. 6)

sin Q^{t}

Z • cos Q^{t) + X sin B^{t)

]

(VII. 7)

The y- component vanished because cp = 0. Similarily one can write

down the higher order multipole terms. The necessary matrix elements

of are given by

P

/ /-t'k^\/f,'k/f\

27
From the last 3j- symbol we see that these matrix elements exist only

if k, I,' satisfy the triangle condition and their sum is an even integer.

Therefore it turns out that, within the no-quenching assumption where

one needs the matrix elements of P V, (t), only a finite number of

n 1 multipole terms exist. The summation index k in equation (VII. 6) has

to obey the condition

1 < k < 2 (n-1) . (VIL 9)

As an example we see that a calculation of the upper state interaction

of Lyman a requires only the dipole and quadrupole terms. This

condition also illustrates the well known fact that there is no ground

state interaction for the Lyman series.

VIIL THE SPHERICAL AVERAGE OF THE TIME

DEVELOPMENT OPERATOR Ul (t)

_(1)

In our evaluation of the thermal average 5 {t)jdefined in

Eq. (VI. 3)jWe first perform the spherical average represented by the

integral over the Euler angles Q, because it greatly simplifies the

remaining integrals over t , p and v. This is due to the spherical

o symmetry of the time development operator Uj(t) defined in Eq. (VI. 6).

It should be noted that this symmetry was achieved by dropping the ion

field exponentials in Eq. (VI. 4), thus replacing Vj^(t) by V^(t). We

will perform this average by means of a rotation technique used by

J. Cooper, 1967, andBarangei, 1958, for S-matrices. Although we

are working with the more general time development operators

Uj^(t) or ']j^j(t), the rotation technique is the same.

In terms of the collision variables p , v, t and Q, the dipole

o interaction between the radiator and a perturber is given by

3^ (VIIL 1

)

V^(t) = e R p + V (t+t ) o .2 2, ,2

/ p +v t+t )

' o

28
(see the appendix of paper II). The three Euler angles denoted by

describe the orientation of the collision frame relative to the atomic

frame. It is therefore convenient to perform a rotation of the atomic

axis, through the angles Q in such a way that R points in the same direction as ^ and the x axis of the rotated atomic frame points in the

same direction as "v. In this rotated frame, the interaction potential

takes the form

V (t) = e^

c Zp + Xv (t+t^)j / j^p^+v^{t+t^)^ j

^ (VIII. 2)

I This rotation transforms the time development operator into a new

I operator U, , where U and U, are related to one another by

^ Ic 1 Ic

\J = & ~^(Q) m) (VIII. 3)

1 1 c

I where & (Q) is a rotation operator (see Chap. 4 of Edmonds, I960). The

j

i time development operator in the rotated frame, U^^, is given by

U, = exp

1 c !- - IP V (f) dt' > . (VIII. 4)

To make the form of U, more explicit, we perform the integral over

1 c

t' and we obtain

J h pv 1^ n

U, = exp <- — — P ZA(t, t ,p, v)-P XB(t, t ,p, v)J (VIII. 5)

Ic JhDv n o n o *

(v/p)(t+t) (vt/p)

A(t,t^,p,v)= ° - . (VIII. 6)

V

l+(v/p )^(t^+t)^ ijlHvtJpy

29
and

B(t, t^, p,v)= - —. (VIII. 7)

^fl7(v/p )^(t^+t)^ ^1 + (vtyp)^

Substituting Eq. (VIII. 3) into Eq. (VI. 3) we see that the integral over

Q in Eq. (VI. 3) involves only the matrix elements of four rotation

operators. Since it is convenient to use spherical states I nJf.m) when

taking matrix elements of j9 (Q) , we make use of the unitary transformation from parabolic to spherical states discussed by Hughes, 1967.

I nqm) = I nf.m ' ) ( n^,m ' j nqm)

/ , /I \ c / , .^(1+m-q-n) / . , (rLf.m' I nqm) = 6 ,(-1) J ^1+^ mm %

2 2

using 3j- symbols and the definitions in the Eqs. (V. 4) and (V. 5). (An

error in the phase factor has been pointed out by H. Pfennig, private

communication). Noti-ng that S5(0) is diagonal in the angular momentum Jf,,

the Q integral in Eq. (VI. 3) may now be written

y(n'q'm'; nq m | Ul (t)| n'q 'm ' ; nq m > dQ = bb bb l aa aa

(^')-l

E/dd(n'q'mMn'^'m') J5 " ^' (n'^, 'm' I U^"'" I n'^ 'm')

J"

\ aa aa mm aclc bd

^ ^ (VIII. 9)

no^ , , (n'^'m'ln'q'm') m'm' ^ ^b b b b

d b

m

(nq^m^ln^^n:^)

1 m, bd lcac aaaai

h d mm J

30
where the summation 7^ denotes sums over ^ , ^

' ^ 9^ rn m' m a a b ' b ' c' c' d

and m'. The Q dependence of the integrand in Eq. (VIII. 9) is contained

d

entirely in the four rotation operators ^ (Q). Using Eqs. (4. 3.2) and

(4, 6. 1 )of Edmonds, 1 960, we obtain the identity

/

a b b a _ & dQ = mm m ' m mm m m ac db bd ca

9 --^ m'-m' + m -m, + M-M' S (-1) ^ " ^ ^ (2L+1)

L, M, M'

(VIII. 10)

aa Uaa i/bb Vbb

• m ' m My \m' m M/ Vm' m , My\m' m, M

,

aa dd ^bb

Hence Eq. (VIII. 9) becomes

J^'^'^^h^W ''%'^bl'^l^^^'''''^I'^a' '''^a'^a^ "

E(n'q'm' 1 n'-f/m') (n'^/m/ I n'q/m/') (nq^ m, ! nP,, m, ) <frLf, m | nq m ) a a a a b b b b b b b b a a a a

^ , m'-m' + m.-m, +M-M'

8^2 (.1) = a d b

(VIILU)

•m'm My Vm' m M/ \-m', m^MyVm' m, M, cc aa dd bb

(n' l' m' ; nf , m , !"7n , (t)l n';/m'; n; m ) ^b d b d UVlc a c a c

31
This result is spherically symmetric; that is, any further rotations of

the atomic coordinate system leave this expression unchanged. One

may verify this rotational invariance by rotating U through some

-1

arbitrary angle Q ' so that U, = ^ (Q')U' j&(Q'). Taking matrix

1 c 1 c elements of the new rotation operators and making use of the orthogonality properties of the 3j- symbols one sees that the right hand side

of Eq. (VIII. 11) did not change. Since we are free to perform further

rotations on U, without altering Eq. (VIII. 11), it is convenient to

1 c rotate the X-Y plane through an angle e = arctg (B/A) where A and B

are given by Eqs. (VIII. 6) and (VIII. 7). This rotation transforms U

into an operator U^^ given by

= exp — e P Zg(t, t ,p,v) h n o

Ic

(VIII. 12)

where

_L

^

pv

g(t, t^, p,v) = yA + Bpv

1 + (v/p) t (t + t) o o

(VIII. 13)

1 - —

y 1 + (v/p)^(t^ + t)^ ^1 + (vtyp)^

The operator U^^ has the important property that it is diagonal in

parabolic states (because it contains only P Z). Hence a typical

n

matrix element of U^^ is given by

(nqml U^(t)| nqm) = exp |-i "|" ^q ^

g (t, t^, p, v)| (VIII. 14)

We also realize that one and the same rotation through the angle e= arctg

(B/A) diagonalizes simultaneously both time development operators

32
acting on initial and final states respectively. As a result a typical

matrix element of the corresponding tetradic operator iJf^(t) is given by

<n'q'm';nqna|'Ul^(t)| n'q'm';nqm) = exp

|

-i -|- (nq-n'q ' )^g't, t^p , v)|

(VIII. 15)

Substituting this identity into Eq. (VIIL 11) the spherical average of the

time development operator "U[j^(t) finally becomes

f(n'^'m'; nq mJ lA (t)| n'q 'm' ; nq m ) dO =

/ bb bb'l aa a a

X^<n'q'm'l nV.'m <n'^ 'm' I n'q' m' > (n'q'm' I n'^/m') (n'^/m/ | n'q/m/) X y aa aa ac cc cc' be bb bb

(nq^mj n^^m^Xn^^mJ nq^m^Xnq^mJ n^^m^> (n^^m^I nq^m^)

Srr^-l) ^ Vl+1)( ^ ^ Y

^ ^ irb 'b^Y ^

^ Um^ m^ Myy-m^ m^ MyVm^ m^ M7Vm^ M

exp { -i-|-(nq -n'qO — g (t, t ,p, v)[ (VIII. 16)

' Z c cm o I

where the unitary transformations are given by Eq. (VIII. 8). This

result greatly simplifies if there is no lower state interaction (e. g. , Lyman lines), in which case one obtains

/

(nq m IU(t)lnq m > dO =\ (na m I n^, m >(n^, m ! nq m > 5 D b 1 a a / J ^b b a a a a a a^ m ma b

pn^^mjnq^m^)j

f^-^^

exp j-i^ nq^ ^ g(t.t^,p,v)j

(VIII. 17)

33
This simplified relation may also be used for the higher series members of the Balmer^Paschen etc. series where lower state interactions

contribute only a negligible amount of broadening to the final line

profile.

IX. EVALUATION OF THE THERMAL AVERAGE

"i^^^(t) FOR HYDROGEN

Having performed the spherical average over the Euler angles

Q we can rewrite Eq. (VI. 3) in the form

^'''^b"^b' "^S^h^ ^ (t)|n'q^m^; m^rn^) = (-1) ^ ^

(ZL+l)

<n'q'm'|n'^'m'><n';/m'|n'q'm')<n'q'm'|n'^'m'><n'^'m/In'q/m

a a aa ac cc cc be bb b

(nq m \nl m )(rLf, m ! nq m )(nq m I n^, m ) (nt, m, I nq, m, ) a a a a a c c c c c b c b b b b

-^c ^c k "^a^/V^c ^c^/V^b^b m/^^*'''''^^'''''^^^

(IX. 1)

where

F(t, n, q , n', q') = Znn / dv vf(v) /dp p /dt $(t, t,p,v) (IX. 2)

c c ej J J

°

and

$ (t,t ,p,v)= exp {-i-^ j^(nq^-n'q^)^g(t,t^,p,v)j-: (nq^-n'qM — g (t, t ,p,v)\-l. (IX. 3)

Thus, the problem is now reduced to evaluating F(t), which will be

done in this section. It is interesting to note the similarity between

Eqs. (Vm. 13) and (IX. 2) and the i]/ -function of Anderson and Talman,

1955, which is the crucial function in their classical adiabatic theory.

34
We first realize that due to the symraetry of the line profile we

only have to evaluate the real part of $ (t, t^ , p , v); that is, for every

positive value of (nq -n'q') there will be the corresponding negative

^ c c value. Hence we are left with

Mt.t^.p, v)= cos

I

|-(nq^-n'qM ^ t^, p , v)

|

-1 . (IX. 4)

In performing the integrals over p and t^ in Eq. (IX. 2) we account for

shielding by setting the interaction potential V and hence also $ equal

to zero whenever the distance of the perturbing electron is larger than

the Debye length D. We also introduce a strong collision cutoff p ^^^^

In principle we can let the impact parameter go to zero because the

functions $ and F(t) do not diverge for small impact parameters as they

do in some second order theories. However, for numerical purposes

this would result in very large computer times due to the growing

fluctuations in the integral. For this reason we will choose p . to be mm

small enough so that when we are interested in large frequency perturbations Auo where perturbers at small impact parameters are quasistatic, the rest of the integral from 0 to p . may then be replaced by mm

the static limit. In the dipole approximation this gives rise to the

-5/2

well known Holtsmark AX -wing (see also Sec. X). According

to the validity conditions of the classical path theories (see Paper I) the

minimum impact parameter p . will be of the order of mm

2

p^ = + n a^ (IX. 5)

where ^ is the De Broglie wavelength.

35
We now concentrate our attention on the integral

/

G(t, p, v) = / dt $ (t, t , p, v). o o

(IX. 6)

For convenience we consider the collision sphere as shown in Fig. 1.

The perturbing electron moves along the classical straight line

trajectory L and we are interested in the interaction from some time

t to some time t + t. Due to the Debye cutoff the t - integral extends

o o ^ o

^ from -T to + T where

(IX. 7)

and the interaction potential vanishes if the electron is outside the

sphere of radius D. The corresponding time integration limit t due to

the strong collision cutoff p . is given by ^ ^mm ^

'

1 rr

(IX. 8)

Based on this model of the collision sphere we split the integral G

into two parts

G(t, p,v)= U(p>p ) -G (t, p,v)+ U(p >p) -Gjt, p,v)

o a o b

(IX. 9)

where the step function U is defined to be

r U(a>b)

1 if a^b

0 if a<b

(IX. 10)

36
In order to evaluate G (t, p, v) we have to distinguish the following

3, four cases depending on whether the initial and final times of interaction are inside or outside the sphere.

Case I:-T<t;t +t<T o o

p. V

1 +—2: t_ (t + t) o o 14

(IX. 11a)

This is the same general expression as given in Eq. (VIII. 13).

Case II: -T < t ; T < t + t o o D p

P

(IX. lib)

Case III: t <-T; t +t<T

o Op.

vr D p

1 —

^

(IX. 11c)

Case IV: t < - T; T < t + t o o

g4(t, t^, p,v) = pv 4 (IX. lid)

D

37
After defining

(t, t , p, v) = cos I ^ (nq -n'q') g ko (t, t p, v) S Id c cm. °k o 1

the integral G is given by

3.

-1

(IX. 12)

T-t T -T

G (t, p, v) = U(2T>t)| /$,dt +/$_dt +/ $_ dt }

T-t

-T

(IX. 13)

+ U(t>2T)^

-T -T-t T-t

where we have separated the cases where the time of interaction is

longer or shorter than the time 2T required to cross the collision

sphere.

In a similar manner we evaluate G^ distinguishing between the

following cases:

Case I: -T<t ;t + t <-t

o o

or T<t;t +t<T o o

g (t, t^, p,v)= g^(t, t^, p,v) (IX. 14a)

Case II: t < t < T; T < t + t o o

g (t, t ,p,v)= g (t, t ,p,v) (IX, 14b)

O w o

38
Case III: t <-T;-T<t + t<-T o o

g (t, t , p, v) = g (t, t , p, v)

O JO (IX. 14c)

Collisions which enter the strong collision sphere are neglected because

o£ the strong oscillations. This yields

Gjt, p,v)= U(T-T>t)

D

r-T-t rT-t fT ' f-T

flollol U.dt +|f dt +1 f 2ol3o dt +1 $ dt

J-rp JX-t •'TT-tt

, . {IX. 15) T-t

+ U (t>T-T)

where again interaction times longer or shorter than (T-t) have been

separated. In the expressions for G and G we realize after a change

a b

of variables that the corresponding integrals over $ and $ are

identical. From the Eqs. (IX. 11a) and (IX. 12) it is also clear that

^ ^ is a symmetric functi(

integral one finally obtains

i -G (t,p, v) = U(2T>t)

is a symmetric function in z = t^ + . Performing the 1^-

^T-t \ (IX. 16a)

[ r$,dt +r $^ dt \

>l 1 o I 2 o(

^^;*'-t/2 •'T-t ^

+ U (t>2T)|Jl^^^o"^ ^"T"

'^^ ' ^ 4I

and

h' G^(t, p, v) = U(T-T>t)| r$,dt +A-dt \

If

+ U(t>T-T)<; l^^dt^ S . (IX. 16b)

39
We now introduce the following dimensionless variables

D

_o

"o D J1^

kT

4un e

e

s = t with uj = */2'

D P V

UU = P

Srrn e

em

y = t_ ,/T

u = — with V

V av

av

1 m

(IX. 17)

and the following abbreviations

R = u • s

P = -r: = T

X - X

(IX. 18)

With these definitions the preceding relations can be rewritten as

gj(s,y,x,u)=

2 2

X + (1-x )y (y+R)

2 2, 2

X + (1-x )y

2

x + (1-x )(y+R)^

g2(s, y, X, u)

g^(s, y, X, u)

xu

1-

X + y (1-x ) V-2

r- 2"

X + y (1 -X )

(IX. 19)

2 r—2^

= — • -t/ 1 -X

xu f 40
and

^(s, y, X, u) = cos

I

C- g^{s, y, x, u)

|

-1 (JX. 20)

where

C = C • C (IX. 21a)

and

C = 4" -n'q') (IX. 21b)

1 Z c c

C = ^ ^= 0. 03043t/-^i^ . J^t^ 2kT Y 2 m- D- V 2kT _ , ^.18 T

av '10

Similarly we have for the integrals over t^

G^(s, X, u) = U(2>R)| r$ ^dy +[ ^^dy

| 'i^ ^1-R (IX. 22a)

+ U(R>2)|J'^$^dy + ( ^ - 1)

and

G^(s, X, u) = U(1>R+P)

1/

^^"^^ "^^1

+ U(R+P > 1)<{a ay}

^ $^ dy> (IX. 22b)

which leads to the thermal average

2rl

F(s, n^, T) = 2TTn^D / du =,u • e " ^ / dxx-^l-x 2-G(s,x. u) (IX. 23)

•'o ^

41
with

, G(s, X, u) = U(x>x )

. G (s, X, u) + U (x >x)- G, (s, x, u). (IX. 24)

o a o b

These integrals have been evaluated numerically using the program

(FORTRAN IV) discussed in Appendix A. This program calculates the

thermal average F as a function of the normalized time s for the

parameters (nq -n'q'), n and T. The upper cutoff is given in units of

c c e the Debye length and the lower cutoff in units of the strong collision

cutoff p of Eq. (IX. 5). oBefore we discuss the methods for obtaining the Fourier transform

of F(t) and the actual intensity profile, it is useful to derive the small

and large time limit of F(t). The small time limit is determined by the

integrals over f ^ and gives the asymptote of the thermal average for

the static wing. The large time limit depends only on the $ ^

integrals

and yields the thermal average as required by the iinpact theory.

In the small time limit $ , reduces to the form

-1

(IX. 25)

where

^ + v^t^ . (IX. 26)

o

This expression depends only on the instantaneous distance r as expected in the static limit and the thermal average is therefore

obtained immediately by the integral over r

f^t), ^=4nn pr^$ (t. r) dr . (IX. 27)

t-»0 el 1 t-»0

•'o

= cos — (nq -n'q') j 2 c c m 2

r

42
3 h t In the small time limit where —(nq -n'q') — ' —- 0 we can

2 c c m Z

then perform the integral with the result c

2 r 1—

n

= " pMrn -n'q' ) — t 2 (IX. 28)

t-*(J 3 el c c m J

For the limit of large times of interaction we have to solve the

integral

p(t) = 2^n I dv vf(v)| dp p t . . (IX. 29)

o For simplicity we set p^ equal to zero (for p 0 see Appendix B).

After a change of variables and a partial integration the integral can

be rewritten as r 4 2 -u^ T "° j F(t) = - 2nn tD^v /I du == u"ee " CIl — ^ —' dz

t-^ e avf

(IX. 30)

The z -integral is known as Raabe's integral (see p. 144 of Bateman, 1953)

and can be expressed in terms of exponential integrals. Furthermore,

from Eq. (IX. 21) we realize that for most practical situations C <<1.

Keeping only the leading term in C we have

(IX. 31)

43
where

2

2nn^ / 3(nq^-n'q'he \2

4C = -7^^ Tt~"^ J and B= 0.27 (IX. 32)

The large time limit of the thermal average in Eq. (IX. 31) is required

for the calculation of the line center and all modern impact theories give

the same result except for the additive constant B whose value depends on the particular cutoff procedure applied. Appendix B gives

a summary of the different constants obtained in the literature which

vary considerably. To what extent this uncertainty shows up in the

final line profile depends on the value of the constant C. The influence

2

will be small if - In (4C ) is considerably larger than the uncertainty

in the additive constant B. Furthermore, the large time limit of the

thermal average affects primarily the center of the line profile and

its contribution vanishes when moving into the line wings.

Finally we show numerical results for F(t) as obtained by means

of the program in Appendix A. Most of the calculations shown in this paper

have been performed for the following electron density and temperature

parameters.

case n^[cm-^] ^e [^] experiment

A 8. 4- 10^^ 12200 Boldt and W. S. Cooper, 1964 (cascade arc)

B 3. 6.

10^"^ 20400 R. C. Elton and H. R. Griem, 1964(T-shock tube)

C 1.3.10^^ 1850 C. R. Vidal, 1964, 1965 (RF-discharge)

44
These parameters correspond to experiments which, as stated already

in the introduction, have revealed the largest discrepancies between

experiment and the modified impact theory. We will concentrate our

attention on the high density case A and the low density case C, since

case B is regarded as being less accurate because of lacking absolute

intensity calibrations.

Figures 2 and 3 show the normalized thermal average F/ F as

a function of the dimensionless variable s = of . t for the cases A and C.

P

Figure 3 shows the results for three different Stark components

specified by the quantum numbers n = nq -n'q' • F is the small

k ceo

time limit according to Eq. {IX. 28) whose Fourier transform leads to

the static wing. The dashed lines are obtained with a lower cutoff

2

n . = 0 ='^+na. It can be seen that for case C the dashed curves mm o o get closer to the static limit F^ than for case A. In order to obtain

the thermal average F for the limit n . 0 the numerical calculations mm

were finally performed with typically n . =^ 0. 01 n so that F , and ^ ^ y ^mm ^o calc

F differed less than about 0. 1% over at least one order of magnitude

o

in s. Far smaller values of s, where F , and F start to differ

calc o again, F is then replaced by F . In this maiuier we obtain the solid

calc o curves in Eqs. 2 and 3 which are used in the following.

It should be noted that these curves are calculated on the basis

of the dipole approximation. It is clear that for impact parameters ^ 2

p < n a higher multipole terms have to be considered. Since the

values s of interest are approximately given by s < of /ao), one expects

P

higher multipole terms to be less important the closer F gets

_ calc —

2

to F for p = n a . This is consistent with the experimental fact

o min o

that in case A an asymmetry of the line has been observed which cannot

be explained within the dipole approximation, while in case C no

asymmetry has been observed.

45
For large s Figs. 2 and 3 show the transition to F as given

00 in Eq. (IX. 31), which forms the basis for the familiar impact theories,

X. THE FOURIER TRANSFORM OF THE THERMAL AVERAGE

Having calculated the thermal average F(t) we now focus our

attention on the evaluation of its Fourier transform

i(Aa)j^) =— / exp (iAojj^s) F{s) ds (X.l)

as required by Eq. (V. 17) (see also Eq. (IV. I6))where the dimensionless variable

A'D^Rip = (AUD-AU).- e)/(ju (X. 2)

is the frequency separation from a particular Stark component

(cf. Eq. (V. 14)) for an ion field strength 6 in units of the plasma

frequency uu . . The thermal average F(s) does not immediately allow a straightforward Fourier transformation because for large s F(s) is proportional to s according to Eq. (IX. 31), hence i{/\UJ-^) diverges. This diver- R

gence is due to the fact that we neglected the finite lifetimes of the unperturbed states involved which naturally terminate the maximum time

of interaction s. This may be taken care of by introducing a convergence

factor exp (-gs) which can be obtained by replacing the delta function

in the power spectrum of Eq. (3) in paper I by a narrow Lorentzian

line with a natural width e (E. Smith and Hooper, 1967). In the final

line profile, however, natural line broadening is always negligible with

respect to Stark broadening which allows us to set e to zero without

46
affecting the shape of the profile. For this reason we will evaluate

On/ i(AuUj^) =^lim

^ / e "^^

e F(s) ds (X. 3)

F(s) is known numerically and there are many ways to perform the

Fourier transform. In order to find the most convenient method we

notice that according to Eqs. (IX. Z8) and (IX. 31), 'F(s) has the following

asymptotes _ 3/2

for s-0: F (s) = p s

' o 1

and for s-»oo: F (s) = p s , 00 2

{X.4)

where

and

2 3 3/2

p = - — n D (2nC) ' (X. 5)

i o e

= - 4^n^ [b - In (4C^] ,

The transition from F to F is very smooth because the power in s O 00 changes only by l/2 over the entire range. It has been found that F(s)

may be approximated by a function G(s) whose Fourier transform can

be given analytically and whose parameters may be determined by a

least square fit. The function G(s) can be given in terms of the series

G(s)=^Gj^(s) , (X. 6)

k

where the number of terms in the series depends on the required

47
accuracy. As a first approximation Eq. (X. 4) suggests

2

a^s

G (s)= . ^ , (X. 7)

2b^s

with

aj_ = P2 (X. 8)

and

bj = i(p^/Pj)2

Gj^(s) has the small and large s behavior of F(s). It then turns out

that

for s-0 F^(s)-G^(s) = p^s^^^

and for s^oo F(s) - Gj(s) = p^ , (X. 9)

where p^ and p^ now have to be determined numerically. Consequently

we take G (s) to be

5

3. S °2'«'=

, 2 ' ,5/2 • '^-l"*

(s ^-Zb^s) '

It then becomes apparent that G (s) is given by

3k-l

=

, 2^^!; ,21c- 3/2

'^•">

(s +2b, s) ' k

48
4k- 3

such that one obtains

for s-»0: G{s) = y p^^ ^ s

^

00 2-k— K

k=l

In this manner the Fourier transform of any G (s) can be expressed

in terms of modified Bessel functions K and K, . For all situations

o 1 calculated it was found that G (s) and G (s) were sufficient to keep

the deviation F(s)-G{s) smaller than 1% for all values of s. In some

situations a fit better than 2% was obtained with Gj^(s) alone. As a further advantage it should be noted that this method tends to suppress

"noise" introduced by the numerical evaluation of F(s).

In the following we evaluate the Fourier transform i(k, AuJj^) of

any G (s) as defined by

1 r -es ^^^R^

i(k,A«;j^)= Hm — / e e G^(s) ds (x. 14)

•'O Their sum will then give us the desired Fourier transform i (Aa)-o)' In

particular we are interested in i(k=l,Au)„) and i{k = 2, Auu^^). We have

2

xAuu ^ s as

-es " R 1 ws +2b s o If 1 J

iAti) „ s

1 r «

49
Introducing

(X. 16)

one finally obtains

2

= a^-b^ (cos Z^- i sin Z^) (X. 17)

where H^^^and H^}^ are Hankel functions and J , J , Y and Y Bessel

o 1 o 1 o 1 functions. These functions like all the other functions used in this

report are consistent with the definitions as given, for example, by

the NBS Handbook of Mathematical Functions (Abramowitz, 1969).

For large arguments is also useful to have the asymptotic

expansion

2

• / X 3 ^1*1 1 (i L .i. 5 ., 3. 5- 7 .,

9- 25- 7 9- 25 - 49* 11

-d-i) +. — (1+i) + - - 2-8 Zl 8^Z^ ) (X.18)

. LJT, 8t , Pi , R

L _L (1^) . WM W I 4 ^pj A.^ 32 \p^; 2

2/ AUJj^

50
Using Eq. (X. 5) for p the latter relation gives us exactly the Holtsmark

-5/2 ^ AX wing for all Stark components

yAuOj^) = nn^D^C^/^- Aa)"j^^^ . (X. 19)

In a similar way one derives

1 -es ^^"^R^

^2^"

^(AuUt,) = i(k=2, Aa)„) = lim

4 e, e —5

E~7^ ds

^ ^

R

rrjf^ (s^2b2s)^/2

(X. 20)

With

^2 " V^R (X.21)

one finally obtains

^2^2 '^h \ n) 2

}

(1 ) 2 + H^^ '{Z^) {IGZ^ + 1282^-3)

a^b^ (X.22)

= —

^ (cos Z^-i sin Z^)

if- 36Z J (Z ) + J.(Z ) (l6Z^-3) - Y (Z ) (16Z^-15) - 28Z Y (Z ) 11 2o2 12 2 o2 2 2-12

+ i

I JJZ^) (16Z^-15) + 28Z2Jj(Z2) -36Z2 Y^(Z2)+Y^(Z2)(l6Z^-3)

51
The asymptotic expansion for large is given by

128Z^ )

(X. 23)

If one requires an even better fit of G(s) to F(s) the general transform

i(k, A«J ) defined in Eq. (X . 14 ) is given by R

4^'02k-2)! 1 . - .k+1 d

(4k-4)!

a, lirrx —( - 1

)

TT

k+1 ^2k-2 / b^(e-iAti3j^) e ^ m (e-iAt«p) k+1 ,2k- 2, d e d b

j^b^(e-iAt«j^)|

"o I "k^" " R

(X. 24)

Finally we want to show that this technique always gives the static wing

according to Eq. (X. 19) for large Atu. For this purpose one has to

perform the Fourier transform of the small time limit of G{s) as given

in Eq. (X. 13).

-^""-R' - Air-^oo ^^^-R' on U A,., = 2 P2k-1 »L n ' R

k=l

e e s "^ds

(Zk+1)!

2k+l P2k-1

2 k!

exp i-^{4+k)l 2^2'

"'J

^^R

-(k+-)

(X, 25)

5 P3

2 Auu„ 4 2

R AoJj^

5- 7- 9 7

8

(1-i) - + R

52
One recognizes that the first two terms are identical with the first

terms in the Eqs. (X. 18) and (X. 23). Hence, we always obtain the

static wing for large AuUj^.

Another important property of i(A(JUT-,) is that for small its

leading terms in the expansion are

^ 1 . <^N-^2'

M.v%^i--o«-K>=-^;;;r- ....

^ ^^^^^

j In this manner it smoothly goes over to the Lorentz profile of the unj modified impact theory.

Before discussing the numerical results of i(AUU^) we first list R

I the constants a, and b, for the cases A, B and C as specified at the end k k

I of Sec. IX. a 11and b are determined from Eq. (X. 8), where p, is

1

1 given by Eq. (X. 5) and p^ is taken from the large time limit of the

computed F(s). p as calculated by the program of Appendix A

2 comp.

may differ slightly from p^ as defined in Eq. (X. 5), if C is not very

much smaller than unity because p^ is based on Eq. (A. 18),

2 comp.

which is true for any value of C and goes over to Eq. (X. 5) for small

C. a and b are determined numerically by a least squares fit.

If —

' The maximum deviations from F(s) obtained with G^(s) alone and with

k Gj(s) + G^ls) are listed too.

53
I

COo

I—

CO

-vrI o

.X o

IT) 00

II

*\ O

00 (Mt ^ o rH

CO

CO

I—

oo

II o o o o 1—

1—

i-H sO un 1—1 tM o CO CO o O t—i CO ro o O o o o m o oI I—

CO o

CO

1—

1

1—4t O

I—

O o CO O o 1—1 o o o

00oo

-MinOI CO

1—

ooV

LD a- 1-1

COI CO

CO

1—1 o

00

CO o

CO

11 ooo

II

-IH

COI

II

IT) CO o CO

00 o o CO

CO o o

no in CO o o o o o o o o o

COo in oV I

CO'so

1—1

00*

II oo

II o in o

vD CO in

1—1 p-H CO CO CO in CO CO 1—1 o o o o CO o o o o o

vDO

in CO

in o

-Vco • r—

I—

C/} ao H

u

Bou

CO

II I—

a

CO cj

54
In presenting the numerical results of i(AtD^) we concentrate on the real

part which turns out to be the most important part. We have chosen

two different normalizations. In Figs. 4 and 5, i(Auu ) is normalized R

with respect to the large frequency limit i {/Sw^) to show the useful

CO K

range of the static theory. The short vertical lines mark the position

of the Weiskopf frequency

Ati) = V V f4-<^^ -^'q') ~l = ^ /C (X. 27)

c av ' \ 2 c c m/ p' ^

for a particular component (nq^-n'q^) which according to classical

arguments determines roughly the range of validity for the static

theory (see p. 321 of Unsold, 1955 and paper II). It should be pointed

out that Auu is usually defined in terms of an average Stark splitting. In

c both cases A and C /\uu describes the range of the static theory very

well. If one allows for a deviation of about 10% at the most from the

static asymptote, Aoj^ may be lowered effectively by more than an order

of magnitude. A more detailed discussion is given later with the final

line profile calculations.

The other normalization with respect to the small frequency

limit i (Auu^) is shown in Figs. 6 and 7 for cases A and C again. These

o R ° ° plots demonstrate the useful range of the unmodified impact theory,

which is based on i (AuUt^) and is expected to break down around the

o R

plasma frequency, as can be seen in Figs. 6 and 7. In order to extend

the range of validity, the modified impact theory makes an impact parameter cutoff at v/auu ( the Lewis cutoff) whenever this is smaller than

the Debye length D; this cutoff accounts for the finite time of interaction

to second order. More details are given in Appendix B. The corresponding function i, . (Auu^) has been included in Figs. 6 and 7. Since Lewis R

the usual derivation of i^ . (AO)-.) is based on the limit of very small C, Lewis R

one expects the best agreement between the Lewis result and our result,

which considers the finite time of interaction to all orders, for the

55
situation with the smallest C. That this is in fact true can be seen

from the low density case with nq^-n'q^ = 3. This component is plotted

again in Fig. 8, in order to demonstrate the importance of G (s) for

those cases where the deviation of G^(s) from F(s) is large (Table I gives a maximum deviation of 13%).

Figures 6 and 7 also contain the static limit i (Auu„) (dashed

00 R

lines) and the Weiskopf frequency /\w . It gives an idea how close the

Lewis results get to the static limit. One notices that with increasing

values of C the deviation of i^ . (Auu„) from the static limit becomes Lewis R

larger. In his line wing calculations (Griem, 196Z, 1967a)Griem adjusts

his "strong collision term" E ^ , in such a manner that the Lewis result

BP

is identical with the static limit at the Weiskopf frequency. In the Figs.

6 and 7 this means that the straight line representing i^ . (auu^) is Lewis R

shifted to the right until it cuts Auu . We use here Au) as defined in

c c Eq. (X. 27) for every individual component instead of the average value

A^^ - kT/(hn ) used by Griem. Since the Lewis line would then lie

appreciably above the curve i(Auo ) one realizes that this procedure R

definitely overestimates the electron broadening as already observed

experimentally (Vidal, 1965; see also Pfennig, Trefftz and Vidal, 1966).

A better method would have been to adiust E ^ , such that i^ . (Auu„)

Bp Lewis R

forms a tangent of the static limit. However, it is clear that any

adjustment of -^g^/ effectively changes the range of the unmodified

impact theory and also defeats the purpose of the Lewis cutoff, namely

to correct the completed collision assumption to second order.

Finally it ought to be emphasized again that except for the time

ordering the Fourier transform of the thermal average i(A«3T3) ^® pre- R

sented here takes into account the finite time of interaction to all orders.

Hence, for small Aw_ it goes over to the impact theory limit and for R

large Ati^-o it gives the static limit without requiring a Lewis cutoff. R

56
XI. THE ONE-ELECTRON LIMIT FOR HYDiROGEN AND THE

ASYMPTOTIC WING EXPANSION

Having obtained the Fourier transform of the thermal average

i(AuJTD) now prepared to calculate the actual line intensity by R

evaluating 1(a), according to Eq. (IV. 15) and averaging it over all

ion fields according to Eq. (II. 1). As explained in Sec. IV. this problem

is greatly simplified in the one electron limit where no matrix inversion is required and the intensity I(A(Ju) is given by

/•oo

I(A(Ju) = UAoj) +

1 P(e) I (Am, B) dQ . (XLl)

I(Auu) is the static ion contribution originating from the first term

i1/Auu^P in Eq. (IV. 21 ) and I(A(D, B) is given by

I(AU0,S)= — y^(nq m jdln'q ' m ' ) ( n 'q ' m ' I dl nq m ^

Tj ^ aa aa bb bb

(XL 2)

dt expfiAtu^tJ (n'q'm'; bbbbnq m l'^^^ ^(t)| n'q'm'; aaaanq m >

using the definitions of Eqs. (V. 14) and (X. 2). The density matrix p

is assumed to be constant over the relevant initial states. With Eq.

(IX. 1) the last expression can be rewritten as

I(Atu,P)= V(n; m Idln';/m')(n'q'mMn';/m'>(n'/, 'm'ldlrtf, m )

' aa a a bb cb cb cb

(n^ cbbb m, Ina m, ) (n';, 'm' I n'q'm') (n'q'm' I n'/ 'm')<f n'^ 'mM n'q'm'^ ac cc cc be bb bb

-m -m,

( nq, m, I nJ , m, ) < nJ^^, m Inqm)(nqm Irtf, m)(-l) (2L+1)

b b b b"^ b c c c c c a c

'I' r L \ /;/ ^ L \ /; ' i. L \ /

' i L\ aa \/aa \/bb \/bb\

li(Aa.j^,6,n,n', q^, q;>q^.q^)

„ 'm'm M/\-m'm M/\-m/m, M/

c c ^^t^ ^ (XI. 3)
where the dipole matrix elements have been transformed from parabolic

to spherical states and the summation over intermediate states I nq m ) a a and In'q'm') has been performed. We next apply the Wigner Eckart

theorem (see Eq. ( 5. 4. 1 ) of Edmonds, I960) to the dipole matrix

elements and replace the reduced matrix elements by the corresponding

radial matrix elements (see Bethe and Salp.eter 1957).

/1 1 I \/l 1 I \

(n^.mld In';/m') = ( - 1 )""^ ij(2^+l )(2^ ' + 1 1 l| )<n^, | r! n'^'>

. Vm|.i myyO 0 0

y

(XI. 4)

Inserting this relation into Eq. (XI. 3) and using the orthogonality

properties of the 3j- symbols we have

I(Auu, B) = "V^ (n'q/m/ I n'

'm/) <nf, m^ | na m N (n'^ 'mM n'q'm') / J bb cb cbbb ac cc

(n'q^mM n 'i;m^) (n't^m^l n'q^m;^ (nq^m^l m^m^Vnt^mJ

nq^m^)

(nq m Inf m > (2;, +1 )(2|/ + 1 )(2^ +1 )(2r ' + 1 ) j < n^. ! r I n ' ) < n ' ccac I r I nl, )

i a a c cl a a c c

\p 0 0/ \ 0 0 0/ \-m' m MV \-m' m M7 \ -m/ m, M/ \ -m/m Mi ^ ^cc/\cc/\bb/^bb

i(Auu^, 8, n, n', q^, q' q , q' ) (XI. 5) R b b c c

58
; J£ we finally replace the unitary transformation by the corresponding

3j-symbols according to Eq. (VIII. 8) the result is

I(AaJ,e)= V(2^ +1)(2^' + 1)(2; +1)(2; ' + 1)(2^ +l)(2f' + l) jL^ a a b b c c

a a \/ c c \l a a \/ b b \/ b b \/ c c

0 0 0/ \ 0 0 0/\-m' m M7\-m'm ]w\ -m/ m^M/Vrn/ M> cc'^ cc/^ bb/^bb

n-1

2

2

n-1

2

m -q

c c m + qc

n^l

y~- ^b

m -q m +q

c c c c ' -r m ; ^ 2 ci

/

I'

c

2 -'"ci

b ^b b ^b , M bjb b ^b , m' A — ;— -m'

2 'b/\ b / \ 2 2 b

/ \ /

/ n'-l n'-l . \

I ~1 — -^a \

' 2

m'-q' m'+q'

c c c c m

n'-l , \

m'-q' m'+q'

c c c c m

2 c/

(np I rl n'; '>(n'^ ' I rl nt ) a a c c

i(^'^R' ^' ^' ^b' %' < ^ (XL 6)

59
The preceding relations hold for the general case of upper and

lower state interactions. They simplify considerably if there is no

lower state interaction (e. g. , Lyman lines). Then one obtains

^ I ^

^ n-1 n-1

MAuo,R)= l(nllrllO)l^y^

n-1

2

m -q m +q

c c c c m

L <

i (Ao)^, B, n, q q ) u -R- be (XI. 7)

with

i uK(AuUp, B, n, q DCq ) = \Km^, B, n, n' = 1, q q ' = 0, q , q' = 0). K bb cc

(XI. 8)

Equation (XI. 7) may be further simplified by evaluating the 3j- symbols

and summing over m, and m with the result. ^ be

I (Am, 8)

u

I(nl|rll0)|

4(n^-l)^n^

n-1

E

q +n

n

2

+(-1)

^22

(n -2q ) b

n-1 q +n

2

*^ 2 2 n+(-l) (n''-2q^)

c q^=-(n-l)L

i (Auu„, B, n, q q ) u XV be

(XI. 9)

60
These simplified relations may also be used for the higher series

members of the other series, whose transitions do not end on the

ground state if lower state interactions contribute only a negligible

amount of broadening to the final line profile.

The foregoing relations for the one electron limit essentially

represent the asymptotic expression for the intensity in the line wings.

If one is interested in frequency perturbations Auu which are significantly

larger than the average ion field splitting Eq. (XI. 1) can be simplified

by replacing the ion field average of the electron contribution by the

electron contribution for the average ion field 6

av

KAud) = I. (AO)) + I (Am, B^^) (XI. 10)

with

OS

0 = / BW(e) de . (XI. 11)

av Jo K A'i) is very much larger than the average ion field splitting, then

according to Eq. (X. 2) Auu^, Am/uu and I(A(U, may be replaced by

° rl p av

I(A(JU, B= 0).

KAo)) = I.(A(u) + I (Aud, P= 0) (XL 12)

In the limit B 0 the Eqs. (XI. 5) to (XI. 95 simplify drastically because

i(A'i)j^) depends no longer on the quantum numbers q^ and q^ which

specify the Stark components shifted by the quasistatic ion fields. This

allows us to sum in Eq. (XL 7) over q^ and m^ which gives us for the

case of no lower state interaction

61
I im, R=o) = u

1 <nll rl 10)

qm

/ n-1 n-1 A

m-q m+q

i (auu, B=0, n, q) u i

{XI. 13)

(n-1)

ni l rl 10 ) V^/ 2 , ,,q+n, 2 ^ 2,\ . , r y n +(-1)^ (n -2q ) i (Am,B=0, n, q)

2(n^-l)n ^\ / ^2(n -1)

q=-(n-l)

For the general case of upper and lower state interaction we can sum in

Eq. (XI. 6) over q^, q^, m^, and M and after applying Eq. (XL 4)

we finally sum over the intermediate spherical states to obtain

I(A(JU ,3=0) = ^^ |(nqmldln'q'm')l^ i(Aa),B=0, n, n', q, q') ,

(XI. 1 4)

m, m

How far into the line center the simplified relations (XI. 10) and (XI. 12)

may be used, depends on the required accuracy. Numerical results,

which compare the asymptotic wing expansions with the more rigorous

unified theory calculations describing the entire line profile, are given

at the end of the next section.

XIL THE UNIFIED THEORY FOR HYDROGEN

In those cases, where the entire line profile including the line

center is required, the line intensities have to be calculated on the basis

of the unified theory. It has to be pointed out that in principal even in

calculating the distant line wings the unified theory has to be used whenever AdJj^ in Eqs. (X. 1) or (V. 17) is no longer large compared to unity.

This will happen in the final integration over ion fields whenever 8 is

close to

B = Aoj/Aoj. (n, q n,' q') (XII. 1)

62
However, it was shown in the last section, that for large Auu one may use one of the asymptotic expansions in Eq. (XI. 10) and (XI. 12).

In the unified theory we have to evaluate the following expression

I(AUJ,e) = — \^(nq m I d| n'q'm') (n'q'm' I dl nq, m, >

TT / J aa aa bb bb

(XII. 2)

(n'q'm'; nq^m^l AO) -£(AllJ ) In'q'm'; bb bb lop op J aa nq m ). a a The matrix elements of Auo as defined in Eq. (IV. 14) are diagonal in

op

parabolic states and are given by

(n'q'm'; nq m I AcD I n

' q ' m ' ; nq m ) = Ac« - AO) .(n, q , n ' , q^) R

^ ^a a a a

" op a a a a i a

(XII. 3)

where AOJ . is defined in Eq. (V. 14). The matrix elements of i!(Aa) ) are

given by

^n'q'm'; nq m I ^(Aa) ) I nq'm'; nq m ) =

b b b b op a a a a

-iTT.(-l)

•m -m,

a b

-.2

j^AaJ-Aa).(n,q^,n',q^) Rj ^(2L + 1

(n'q'm'In'1'm'Xn't'm'ln'q; m ' ) ( n ' q 'm '

I n '^.^m^X n 't^m M n 'q^m

aa aa ac cc ^

(nqaaaaaccc m 1 nt m ' ) ^ nq m I n^ m ) (n; m ! nq^ m^) ccu«- ^

/

I' I L

a a

\ -m' m M'

I' I L

a a

• m' m M a a

-ro.' m M' • m/ m, M

b b

i

( Atw j^, e, n, n, q^,q

(XII. 4)

63
Eqs. (IV. 16), (IX. 1) and (X. 1). This relation simplifies significantly in case of no lower state interaction in which case we need the

nnatrix elements

(nq.m li:(Auu )Inq m ) = - in AtJU -Auo .(n, q ) R b b op a a i ^ I

/ -T,—Ti— ( ^^1."^!,

^ m Inq m > (nt m Inq m > '^a b b a a a a a a a c c c J ^ , q , m a c c

i (AuJt3. B, n, q q ) u R be (XII. 5)

Due to the delta function the matrix of the operator £ is then block

diagonal in m, which reduces the size of the matrices to be inverted

to n X n or (n-1) x (n-1) depending on the quantum numbers n and m.

Furthermore, Eq. (XII, 2) simplifies in case of no lower state interaction.

After transforming the dipole matrix elements from parabolic to spherical

states, applying the Wigner-Eckart theorem (see Eq. (XI. 4)) and using

Eq. (VIII. 8) one obtains

l{m,Q)= UnllrllO)!^

^

(-1)

n+m- 1 -

1 n-1

m2 2

1

- m

n-1 n-1

\

\ m-q m+q

\ b 1 - m

/

Im (nq.ml [auu -Jl(Auu )| Inq m) . (XII. 6) b j op op I

' a

In order to keep the mathematics simpler we concentrate in the

following on the case of no lower state interaction, because it covers the

experimental situations of case A and B and is also a good approximation

to the higher Balmer lines of case C (see the list of references at the

64
jend of Sec. IX). Including lower state interactions means at this stage

jonly a more extensive summation over 3j- symbols because the crucial

jfunction i(Ati)j^, i^), the Fourier transform of the thermal average, has

jalready been evaluated for the general case of upper and lower state

{interactions.

Using the unitary transformation of Eq. (VIII, 8), Eq. (XII. 5)

imay be rewritten as

n+m- 1 - (nq,ml£(Auu )Inq m) = - iTT(-l) b op a

|^Atu-/\uu.(n, q^) ^+1

)

/ n-1 n-1 \ / n-1

2

m-q m+q

a a \

-m

/

n-1

m-q^ m+q^

7

q > 0

c

(Amj^,8,n, ^^^^

n-1

m-q m+q c c c c m

(XII. 7)

where we have used the fact that i uK 8, n, q^bcuK , q ) = i (Aa)_, S, n, q,be, -q ) and that i (AuUt^, 6, n, q, , q = 0) = 0. We also realize that u R be

<'nq -mi<C(Auu )l nq -m> = (nq,ml£(Auu )|nq m)

b op a b op ' a

(XII. 8)

A computer program (FORTRAN IV), which evaluates I(Aud, B) according

to Eqs. (XII. 6) and (XII. 7) and also performs the final ion field average

according to Eq. (II. 1) is presented in Appendix C. The ion microfield

65
distribution function employed is the one given by Hooper, 1968a, 1968b,

which differs less than about 1% from the values determined independently by Pfennig and Trefftz, 1966.

For the experimental parameters of case C, Figs„ 9 and 10

show numerical results of I(Auo, P = 10) for n = 6 and n = 10. The fat

vertical lines indicate the relative intensities and the positions of the

Stark components for the static field (ion field) 3=10 and it demonstrates the electron broadening.

Figures 11 to 13 show the final line profiles I(At«) after performing the ion field average for the experimental cases A, B and C

(see end of Sec. IX). As a first result it turns out that for numerical

accuracies of about 2% it is in all 3 cases sufficient to consider only

Gj(t) meaning that i(A(Dj^, B, n, q^, q^) may be replaced by

ij(Aujj^, 6, n, q^, q^) as given in Eq. (X. 17). Although according to

Table I, G^(t) may differ from F(t) for some components of case C

by up to 13%, it turns out that after summing over all Stark components

and averaging over ion fields this difference F(t) - Gj(t) is apparently

smeared out over the entire line profile and affects the final line

profile by not more than about 2%. This is very convenient for

practical calculations, because it no longer requires an extensive

evaluation of the thermal average anymore, but for most practical

situations it is sufficient to calculate the line intensities directly on

the basis of Gj^(t) whose specifying constants a^^ and b^^ are given

immediately by the Eqs. (X. 8) and (X. 5).

This is even more true in view of the fact that the final line

profile is partially affected by an uncertainty in the constant B as defined in Eq. (IX. 31) or (B. 19). As summarized in Table II of

66
Appendix B its actual value depends on the cutoff procedure applied,

a problem, which has not yet been solved satisfactorily. The upper

the frame work of the classical path theory (see also Chappell, J. Cooper and E. Smith, 1969). The lower cutoff parameter, however,

which essentially replaces the dynamic strong collisions not amenable

in a classical path theory, can only be determined conclusively from a quantum mechanical theory which is also able to handle strong collisions

and which does not yet exist. The constant B adopted here is based on

3 2

a lower cutoff parameter p . =^ + -r- n a , which specifies approxi- mm Z o mately the region of validity for the classical path theories (see paper

I). Numerical results based on other values for the constant B as used in the literature (see summary of Appendix B) are also included

in Figs. 11 and 12 for the cases A and B. The largest value B = 1. 27

is the one adopted in the recent calculations of Kepple and Griem,

3 2

1968, while the smallest value of B is obtained for p . = + -r-n a mm 2 o and choosing an upper cutoff of p ^* ^06D as proposed by Chappell,

J. Cooper and E. Smith, 1969. For case C this variation of the constant B does not show up in Fig. 13 and amounts to an intensity change

of at the most about 4%. These variations indicate the reliability of

the classical path theories and demonstrate that for some cases the

error estimates given in the literature are too optimistic. The effect

on the final line profile due to the uncertainty of the constant B will

be small if either according to Eq. (IX. 31)

cutoff parameter a = p /D (see Appendix B) and therefore also the

'V (t') dt' can in principal be decided within

-2 In (2C) » 1 (XII. 9)

67
or if (Like for the higher series members) the number of Stark components is large which tends to smear out the influence of the constant

B. It should be pointed out again that the unified theory is intrinsicly

normalized independent of the value of the specifying constants of a particular line. Hence, any variation of the constant B does not affect

the normalization of the line profile.

In comparing the numerical results for case C, with the

experiment it has to be kept in mind that we are comparing the higher

Balmer lines with calculations for the higher Lyman lines, because

our final line profile calculations have not yet taken into account lower

state interactions. This means that in a plot of the intensity versus the

wavenumbers which is essentially an energy scale, the line profiles

cannot be expected to coincide because of the difference in the Stark

effect. This gives rise to different static wings as explained in detail

by Vidal, 1965. Hence, we have to rescale the Lyman profiles preserving normalization in order to be able to compare the measured profiles

of the Balmer lines with the calculated profiles of the Lyman lines.

This means that in a plot of log I versus log /\v we can compare the

line shape of the corresponding lines directly. The agreement is

remarkable. For the higher lines, n > 8, where Doppler broadening

was shown to be negligible and where lower state interactions no

longer affect the line shape noticeably, the agreement is better than

2% over the entire measured line profile, which for n = 8 extends over

3 orders of magnitude in intensity. In particular, the calculations show

-5/2

also the surprisingly large range of the f\iu - wing, which extends

to l/lO of the maximum intensity. This fact is not explainable by a purely static theory considering also shielding effects. For the lower

lines the calculated profiles have to be folded into a Doppler profile in

order to achieve similarly good agreement. For the lower line we

68
also expect in the line center some influence due to lower state interj actions on the line shape, which is partially removed again by Doppler

broadening.

-5/2 A more detailed study of the Aut) -wings reveals some other

interesting facts. In Fig. 13, the dashed lines indicate the asymptotic

5/2 5/2

A(D -wings; except for n = 5 and n = 6, what appears to be a -wing in the measurements and calculations is not the asymptotic

-5/2 Holtsmark /\(U -wing in the region of interest. If one extends the

calculations to even larger frequencies At)U> a.11 the wings will , eventually

coincide with their asymptotic limit. In the paper of Vidal, 1965,

Table II gives a list of the electron densities, which were evaluated

-5/2 under the erroneous assumption that the measured Atw -wing was

the asymptotic Holtsmark wing; it was stated that for to ^-^^ electron densities coincide within i 4%, A more careful analysis of

the values, which have been plotted again in Fig. 14 reveals a systematic trend. For large and very small principal quantum numbers the

electron density values rise above the average value, while the minimum value was obtained for n = 7. From Fig. 13, it now becomes

apparent that the electron densities based on the asymptotic

Holtsmark wing will go up for increasing n. For smaller n the

quantum number dependence of the electron density is masked by

I Doppler broadening which raises the wings again and explains the

increasing values of electron density for small n. Another important

result can be seen from Fig. 13. For small principal quantum numbers

the line intensities are much smaller than predicted by a quasistatic

^ theory. This was observed first by Schluter and Avila, 1966 and the

effective electron densities for a quasistatic theory as a function of AX

show the qualitative behavior measured by them after unfolding the

69
Doppler broadeningo This observation together with the measurements

of Boldt and Cooper, 1964^ suggested the semiemperical procedure

proposed by Edmonds, Schluter and Wells, 1967. A detailed quantitative

comparison requires for the first series members a consideration of

lower state interactions, which is in process.

For the parameters of case C, Kepple and Griem, 1968, have

already calculated the lines H and H . These calculations have been

o 7

extended to H^^^ by Bengtson, Kepple and Tannich, 1969, using identically

the same computer program. The results are plotted in Fig, 15 and

comparing the line shape for the higher series members, for which

lower state interaction becomes negligible, with our results in Fig. 13

one realizes a significant difference . In particular, their calculations

-5/2

do not reveal the AX decay in the near line wing for intensities

smaller than about l/lO of the maximum intensity at /\uu = 0 which is

discussed above. It should be pointed out that the ion field dependent

cutoff, which has been introduced by Kepple and Griem, 1968, to account

for the usually neglected exponential in Eq. (VI. 4) cannot be responsible

for it. This has been tested in our calculations. One can understand

this by realizing that for the higher series members the effect of

dynamic broadening due to the electrons as described roughly by the

constants p^ in Eq. (X. 5) turns out to be much smaller than the halfwidth of the total line, which is essentially determined by quasistatic

broadening. - As another interesting result, Fig. 16 shows a plot of a calculated

Lyman-p profile for two different values of the constant B(B = 1, 27 and

3 2 , B for p . = ^ + —— n a ), which allows also some qualitative state- mm 2 o ments concerning H^. We realize that changing B affects the very line

center, where the profile shows the two humps and the near line wing,

but it does not change the intensity around the halfwidth significantly,

70
which may be understood as an effect of the normalization. This is in

agreement with experimental observations of Wende, 1967, which show

that the calculations of Griem, Kolb and Shen, 1962 overestimate the

near line wing. It also explains the good agreement of experimental and

theoretical halfwidths in high density plasmas (see Gerardo and Hill,

1966) because the line intensity around the half width is rather insensitive

to the exact value of B.

Finally in Fig. 17 to 19, we compare the unified theory calculations (solid curves) with calculations based on the one-electron

theory in order to see how far into the line center the asymptotic

wing expansions as given in Eq. (XL 10) or (XI. 12) may be used. In

all Figures the short vertical line indicates the position of the outermost, unperturbed Stark component for an average ion field B^^.

which is given by Ati) = 6 Auu . (n, q = n- 1 ) where Auu . is defined in

av 1 1 Eq. (V. 14). The dashed lines correspond to the one-electron theory

calculations for B = 0 according to Eq. (XL 12), while the dash-dotted

lines give the results for 6=R according to Eq. (XL 10), First of

av

all we realize that, as expected, the one-electron result for R = B av

diverges when Auu approaches 6 ^ti) .{n, q = n-1). However, in all three

3. V X cases we see that for frequencies

Auu > 5 B Auu. (n, q = n-1) (XILIO)

av 1

the one-electron theory calculations according to Eq. (XL 10) coincide

with the unified theory calculations to within 1% and better. For slightly

released accuracy requirements one may also apply the simpler oneelectron theory calculations based on Eq. (XI. 12) with 3=0. In

particular we see that for small principal quantum numbers the useful

range is very much larger than for the one-electron theory calculations

71
with B = because, for p = 0, the one-electron theory diverges only

at AtA) = 0. We also realize that for the line intensity range of practical

interest both asymptotic wing expressions with B = 0 and S = B av

become less useful with increasing principal quantum number.

72
APPENDIX A

PROGRAM FOR CALCULATING THE THERMAL AVERAGE F(t)

This section gives a complete listing of the program which was

used to calculate the thermal average F(t) as discussed in Sec. IX. In

I order to understand the program the following explanations may be of

help.

1. Calculation of $ (s, y, x, u)

Function PHI( KAY) calculates $ , $ and ? as defined by

the Eqs. (IX, 17)-(IX. 21). In order to assure sufficient numerical

accuracy for computing $ (s, y, x, u), series expansions have been

applied whenever one of the different g (s, y, x, u) becomes very small.

The following expansions have been used abbreviating

w

= Yp + V t^ /D = X + y (1-x )

= p/ "Yp + v t^ = x/r

/ 2 2

' ^ / = vt / + V t = ^jl-x y/r

(A.l)

,J"2 2 2

' / and Y

- '^^/ TP +vt =us/

(a) g^ {s, y, X, u)

a < 0. 01 and y > 0

2 3 = ^s

^ 1 + V^^+v) ^ + y(16+80y+40y +7y ) ^4

r^(l+Y) ( 8(1+y)^ 128tH-Y)^

73
(96+448y+1024y^+1520y^+880y^+252y^+33y^ ) 6

1024(1 + y)

(A. 2)

a < Oo 01 and y < - 2

2 (

^

2 y^(16-5y) 4 y ^32-192y + 112y^-21y^)

"""^ ( 8(1-y) 128(1-y) 1024 (1-y)

(A. 3)

Y < 0, 01

-2 7 3

1 - w-Y - (3-llw ) \ + (9-17w )

+ (31-350w^+447w^) + } (A. 4)

Y > 100

8

4

128

I

V2(l-wy ^ L (1+w ) r 5w-l (3+4w-15w^) xu • I

" 2y [

" 4y " 3^2

(ll-83w-59w^+195w^)

3

64y

(b) g^is.y, X, u)

2

(x/y) < 0. 0002 and y > 0

]!

(A. 5)

74
3 6

- (187-443y+297y^-33y^ • + . . . 128

(A. 6)

(x/y)^< 0. 0002 and y < 0

.2 , V 2 . .2 , V 4

2

2 xu

(A. 7)

I 1-yl < 3. lO"^

!

2 3 + (l-x^) (l--^x^)(l-y) +«..> (A. 8)

2. Calculation of G(s, x, u)

Function GN(AX, AXS) calculates G(s, x, u) as defined by the

Eqs. {IX. 22) and (IX. 24). This integral has been solved by rewriting

it in the following way.

G(s, X, u) = G^(s, X, u) + G^Cs, X, u) + G^(s, x, u) (A, 9)

The new integrals G (s, x, u) are only functions of $ (s, y, x, u) and

turn out to be

1-R 1-R

G, = U(x> X ) U(2 > R) I $, dy + U(x > x) U(l > R+P) | $,dy

1 o J1 o J1 R P

"2 (A. 10a)

75
= U{l-R > K) dy with K = P if X < X —R

if X > X

2 o

+1

G2=U(x>x^) j^U(2 >R) j^^dy + \J(R>

2) J

^^j

1-R 1 -1

1 + U{x^>x) j^U(l > R+P)J* $^dy + U{R+P >l)y $2^^^

1-R

1 I

G = I $ dy with B = MAX (1-R, K) and K=

B

P

- 1 if X >x

o P if X < x

and

G = U,x>x^).u(f>l) -(f-l) 4

The first integral is split up at most into three parts

-X +x 1-R

G

r

lj=J ?jdy+J $jdy+J l^dy

K -X

(A. 10b)

(A. 11a)

(A. lib)

(A. 12)

(A. 13)

if every upper limit is larger than the lower limit using a different

convenient change of variables in every part. The integration is performed in all cases by means of Weddle's rule (subroutine WEDDLE).

The number of points besides the initial point is firstly taken to be 6

and is doubled in every iteration by calculating only the values in between

the old values until the integral changes less than a preset relative

value called DG.

76
3. Function GUS (Au) calculates the x-integral and is defined by

GUS(AU)= =4= u^»e"^ H(s,u) (A. 14)

where 1

H(s, u) = zfx-Vl-x G(s,x,u)dx (A. 15)

IJ x^l-x^

'

(

o The integral is performed in two parts

x^ -^n(x^)

H(s,u) =

/*^1-z 'G{s,yJT^u) dz + zf e'^^-^l-e'^^ G(s, e"^, u) dz

o o (A. 16)

by means of Weddle's rule again using a technique as described before

for G(s, X, u). The new test parameter, which determines the number

of iterations, is called ACC and is given in the main program. The

old test parameter DG has been made a function of ACC, x and u in

order to calculate those values of G(s, x, u) with the highest accuracy

which give the largest contribution in the final velocity integral. For

very small values of u* s the integral has been approximated by

Z 2 H (s, u) = - C • s a

(A. 17)

This limit does not lead to the static wing, as F(t)^ ^ in Eq. (IX. 28)

does because of the lower cutoff x . As explained in Sec. IX. this cut- o

off is necessary in order to avoid rapid fluctuations in the integrals

which require more integration steps and longer computer time. For

practical calculations x is chosen small enough so that over at least

o one order of magnitude in t the calculated F(t) comes as close as required to the asymptotic limit F(t) . t -*0

77
For us > 2 the following relation is applied

1

2 f

xyi-x^'G^(s,x,u)dx= -C sf -^iiLi|l£.) dz

J J 1+z

(A. 18)

where

+ ^(l-b K (b) + b^K (b) ) 3 1 o

u

For simplicity x has been set to zero for this relation without

o affecting the final result noticably. For very large s Eq. (A. 18)

leads to Eq. (IX. 31). However, Eq. (A. 18) does not require C to be

small as Eq. ( IX. 31) does. The Raabe integral is calculated by the

function SNZ and the modified Bess el functions K and K by the

o 1 subroutine BESMOD.

4. The final thermal average F(s, n^^, n^, T) as defined in Eq.

(IX. 23) is calculated in the main program FSTEST. The best results

for the velocity integral have been obtained by Gauss's quadrature

formula (function GLQUAD). The values FS in the program are given

by

FS= F(s,n,, n , T) . (A. 19)

2nn D-^ ^

®

e The main program reads in the temperature T, the electron density

n , the quantum number n , two cutoff parameters, which specify

e k

the lower and upper cutoff of the x- integral by

D . = STRONG • p ^"^^^ ° (A. 20)

p = CUT . D

max

78
where is given by Eq. (IX. 5). Finally the main program reads in

the initial value of s for computing the thermal average and the number

of values which proceed according to s = s • 10 . The program also

x^"r J, xC gives the asymptotic thermal average leading to the static wing which

is called GS and calculates the relative thermal average in units of

this asymptotic value. Furthermore, the function KLOCK provides

a means to test the computer time for every individual value of FS.

The results are shown and discussed in Sec. IX.

79
PROGRAM FSTEST

C CALCULATION OF THE THERMAL AVERAGE FS

C INTEGRAL FS OVER U IS TESTED FOR DIFFERENT L IN GAUSS QUADRATURE

COMMON/ PDS/ S/ PDCON /C 1 CZ » CON BCON /PDR AD/ BR AD / PDACC/ ACC DG » NGN

COMMON /PDSTR/ STRONG

EXTERNAL GUS

500 READ(60»100) T EMP » DENS T Y , QNUM Q , CUT STRONG » SS K

IF ( DENSTY.EQ.O. ) CALL EXIT

WRI TE { 61 104 ) TEMP* DENSTY* QNUM, Q, CUT» STRONG

T = TEMP*l.E-4

RELDEN = DENSTY^l .E-18

CI = 1.5 * Q * QNUM

C2 = 0.03043 ^ SORT ( RELDEN )/( T * CUT)

WRITE (61,105) C1,C2

NK = QNUM * 0 + 0.01

CON = Cl * C2

BCON = 1.414213562 * CON

CT3 = CUT 3

BRAD = 7.67E-3*SQRT (RELDEN/T )*QNUM^-*2*STR0NG/CUT

NUM = 0

ACC = 3.E-4

7 DO 1 I = 1,K

NUM = NUM + 1 SS = SS * 10.0

S = SS / CUT

FSOLD = 0.0

GS = -1.671085516 ^ (CON*CUT*SS) 1.5

GSAS = -( CON->^SS ) **2* ( 1 . 1 2 8 3 79 1 6 7 1 / ( C2 *ST RONG ) - 1.) * CUT

PRINT 200, SS, GS, GSAS

DO 20 L = 3*5

LLL = KLOCK(O)

FS = GLQUAD ( GUS , 0 . 0 , 5 . 0 , L ) * CT3

FFGG = FS/GS

DACC = ABS( (FSOLD - FS)/FS)

LLL = KLOCK(O) - LLL

PRINT 300, FS, FFGG, DACC LLL

IF ( DACC.LT. ACC ) GO TO 50

20 FSOLD = FS

50 PUNCH 400, SS, FS, DACC, NK , CUT, STRONG, DENSTY, TEMP, NUM

1 CONTINUE

GO TO 500

100 FORMAT( 2E10.2, 4F10.1, ElO.2, 15)

104 FORMAT (1H1»* TEMPERATURE = *E14.5, lOX, * DENSITY = *E14.5//

1 * QUANTUMNUMBERS N = *F6.1,10X,* (Nl - N2 ) = »-F6.1//

2 * DEBYE CUTOFF FACTOR = *F8.3,10X,* STRONG COLLISION FACTOR = 3 F8.3)

105 FORMAT (//* CI = * E14.6 , lOX C2 = * E14.6)

200 FORMAT {///^ S = *E15.7,9X,* GS = *E17.9,9X,* GSAS = *E17.9/)

300 FORMAT (* FS = ^E17.9,* FS/GS = *E17.9,* DACC = *E17.9,

1 * LLL = *I8/

400 FORMAT (ElO.3, E15.7, ElO.2, 15, F7.1, F7.3, 2E11.3, 14)

END

C

c

80
FUNCTION GUS(AU)

COMMON/PDS/S/PDCON/Cl »C2 jCON » BCON/PDU/U/PDXM/XMl yXfAZ

COMMON/ PDRAD/BRAD/PDACC/ AC C »DG» NGN/ PDSTR/ STRONG

DIMENSION F( 192 ) H ( 192 ) GUS = 0.0

U = AU

IF(S.EQ.O.O) RETURN

XMl = C2*STR0NG/U

IF ( XMl .GE . 1 . ) XMl = 1. - l.E-9

XM2 = XMl + BRAD

IF(XM2.GE.l. ) XM2 = 1. - l.E-9

CONS = CON * S

U2 = U * U

FU = 2.2567583342 * U2 * EXPF(-U2)

US = U * S

IF( (US/STRONG) .GT. 2. E-6) GO TO 5 GUS = -CONS * CONS * (l./XM] - 1 . ) * FU

RETURN

5 NGN = 0

GD = 0.0

IF(US.LT.2. ) GO TO 8

NGN = 1 PA = 2. * CON/U

XMM = XM2 * XM2

SXM = SORT ( 1 . - XMM

)

CALL BESMOD (PA, F I 0 , F 1 1 , FKO , FK 1

)

GD = -CONS * SNZ(PA) + 2

. * ( 1 . +PA*

(

PA*FK0 - FKl))/3.

8 GUOLD = 0.0

GIOLD = 0.0

G20LD = 0.0

N = 3

DO 100 K = 1»6

GUS = GD

DG = ACC/FU

N = 2 * N

NN = 2 IF(K.EQ.l ) NN = 1 AN = N

ANN = NN

I F ( K. EG. 1 ) GO TO 10

DACC = ABS{ DGl /GUOLD

)

I F ( DACC.LT. ACC ) GO TO 40

GO TO 20

10 FO = GN(0.0»1.0)

DQl = XMl * XMl

20 DQO = DOl/AN

DQ = DQO * ANN

Q = DOO - DQ

DO 30 J = 1»N,NN

Q = Q + DQ

X = SORT(Q)

XS = SORT (1. - Q)

30 F(J) = GN(X,XS) * XS

81
CALL WEDDLE ( DQO , N F » GUSl , FO

40 GUS = GUS + GUSl

IF(K.GE.6 ) GO TO 50

DGl = GUSl - GIOLD

GIOLD = GUSl

IF{ DGl. EG. 0. ) GO TO 50

DO 45 J = 1»N

L = N + 1 - J

45 F(2 * L) = F(L)

50 I F ( XMl .GE .0.999999999 ) GO TO 90

IF(K.EQ.l ) GO TO 60

DACC = ABS{ DG2/GU0LD)

IF(DACC.LT.ACC) GO TO 80

GO TO 6 5 60 XO = SQRTd. - DQl)

HO = GN(XM1»X0) * XO DOl

DYl = -LOGF ( XMl

65 DYO = DYl/AN

DY = DYO ^ ANN

Y = -DYl + DYO - DY

Nl = N - 1 DO 70 J = 1»N1»NN

Y = Y + DY

X = EXPF(Y)

X2 = X * X XS = SQRTd. - X2)

DG = ACC/XS

70 H(J) = GN(X»XS) * XS * X2

H ( N ) = 0.0

CA'.L WEDDLE ( DYO » N H OUT HO ) GUS2 = 2.0 * OUT

80 GUS = GUS + GUS2

IF(K.GE.6) GO TO 90

DG2 = GUS2 - G20LD

G20LD = GUS2

IF(DG2.EQ.O. ) GO TO 90

DO 85 J = 1»N

L = N + 1 - J

85 H(2 * L) = H(L)

90 DACC = ABS((GUS - GUOLD)/GUS)

IF( DACC.LT.ACC ) GO TO 120

100 GUOLD = GUS

120 GUS = GUS * FU

RETURN

END

C

C

C

FUNCTION GN{AX»AXS)

COMMON/ PDS/S/PDU/U/PDCON/ CI C 2 CON , RCON/ PDXD/X , XS X 2 » DX

C0MM0N/PDACC/ACC»DG»NGN/PDXM/XM1 XM2

DIMENSION F(768), G{768)» H(768)» D(768)

GN = 0.0

X = AX

82
XS = AXS

X2 = X * X 0X2 = 1. - X2

GiMOLD = 0.0

GIOLD = 0.0

G20LD = 0.0

G30LD = 0.0

G40LD = 0.0

R = U * S / XS

Rl = 1. - R R2 = 0.5 * R PI = -1.

P2 = -1.

N = 3

DO 100 K = 1*8

GN = 0.0

N = 2 * N

NN = 2 IF(K.EQ.l ) NN = 1 AN = N

ANN = NN

QA = -R2

IFCX.GE.XMl) GO TO 5 IF(K.EQ.l) PI = SQRT(XM1*XM1

QA = PI

5 IF(Rl.LE.Pl) GO TO 50

QB = -X

IF(Rl.LT.-X) OB = Rl

IF(QA.GE.QB) GO TO 20

IF(K.EQ.l ) GO TO 7

DGG = ABS(DG1/GN0LD)

IF( DGG.LT.DG) GO TO 13

GO TO 8

7 DQl = (QB + R2) ( 1./3. ) 8 DQO = DQl/AN

DQ = DQO * ANN

0 = DQO - DQ

DO 10 J = 1»N»NN

Q = Q + DO

02 = Q * 0

Y = QR - 0 * 02

10 F ( J ) = PHI ( 1 »Y) * Q2

CALL WEDDLE ( DQO » N » F » OUT , 0 . 0

Gl = 3.0 OUT

13 GN = Gl

IF(K.GE.8) GO TO 18

DGl = Gl - GIOLD

GIOLD = Gl

IF(DG1.EO.O» ) GO TO 18

DO 15 J = 1»N

L = N + 1 - J 15 F(2 * L) = F(L)

18 QA = -X

20 QB = X
IF(R1.LE«X) QB = Rl

IF(QA.GE.OB) GO TO 35

IF(K.EQ.l ) GO TO 23

DGG = ABS(DG2/GN0LD)

IF(DGG.LT.DG) GO TO 31

GO TO 25

23 FO = PHI ( 1 ,QA

DY2 = QB - OA

25 DYO = DY2/AN

DY = DYO * ANN

Y = OA + DYO - DY

DO 30 J = 1»N»NN

Y = Y + DY

30 G ( J ) = PHI { 1 Y

CA'-L WEDDLE ( DYO » N , G » G2 » FO ) 31 GN = GN + G2

IFIK.GE.S) GO TO 33

DG2 = G2 - G20LD

G20LD = G2

IF(DG2.EO.O. ) GO TO 33

DO 32 J = 1»N

L = N + 1 - J

32 G(2 * L) = G( L)

33 OA = X 35 IF(QA.GE.Rl) GO TO 50

IF(K.EQ.l) GO TO 36

DGG = ABS(DG3/GN0LD)

IF(DGG.LT.DG) GO TO 42

GO TO 37

36 HO = PH I ( 1 Rl ) * Rl * Rl

QO = l./Rl

DQ3 = l./QA - QO

37 DQO = DQ3/AN

DO = DQO * ANN

Q = QO + DQO - DO

DO 40 J = 1»N»NN

0 = Q + DQ

Y = l./Q

40 H( J) = PHI ( 1 »Y) ^ Y * Y

CALL WEDDLE ( DQO N H G3 HO

42 GN = GN + G3

IF(K.GE.a) GO TO 50

DG3 = G3 - G30LD

G30LD = G3

IF(DG3.EQ.O. ) GO TO 50

DO 45 J = 1»N

L = N + 1 - J

45 H(2 * L) = H(L)

50 IF(R.LT,l.E-6) GO TO 80

OA = Rl

IF(X.GE.XM2) GO TO 55

IFCK.EQ.l) P2 = SORT(XM2*XM2 - X2 ) /XS

55 IF(R1.LT,P2) QA = P2

IF(K.EQ.l) GO TO 57

84
DGG = ABS ( DG4/GN0LD

)

IF( DGG.LT.DG) GO TO 63

GO TO 58

57 DA = PHI(2»QA)

DY4 = 1. - QA

58 DYO = DY4/AN

DY = DYO * ANN

Y = QA + DYO - DY

DO 60 J = 1»N»NN

Y = Y + DY

60 D{ J) = PHI (2»Y)

CALL WEDDLE ( DYO » N D » G4 DA

)

63 GN = GN + G4

IF(K.GE.8) GO TO 70

DG4 = G4 - G40LD

G40LD = G4

IF ( DG4.EQ.0. ) GO TO 70

DO 65 J = 1»N

L = N + 1 - J

65 D(2 * L) = D(L)

70 IF(R2.LE.l. ) GO TO 80

IF(K.EQ.l) G5 = ( R2-1. )*PHI (4»Y)

GN = GN + G5

80 DGG = ABS((GN - GNOLD)/GN)

IF( DGG.LT.DG) GO TO 120

100 GNOLD = GN

120 IF(NGN.EQ.l) GN = GN - G5

RETURN

END

FUNCT ION PHI ( K »AY

)

COMMON/ PdS/S/PdU/U/PDCON/ CI » C 2 » CON , BCON/ PDXD/ X » XS , X 2 » Y = AY

PHI = 0.0

FAC = 0.0

IF(K.GT.2) GO TO 30

RI = l./SQRT(X2 + Y * Y -i^ DX2 ) IF{K.GT.l ) GO TO 10

C = X * RI

W = XS * Y * RI

G = u * S * RI

IFCC.GT.O.Ol ) GO TO 3

IF ( Y.LT.O . ) G = -G

AG = 1 . /ABS ( 1 . + G)

IF( AG.GT.2. ) GO TO 70

GI = (C * AG) 2 IF(G.LT.-1.) GO TO 1 FAC = CON * S * RI * RI * AG

ARG = FAC * (1. + GI * G * 0.125 * (4. + G))

I F ( C.LT. 1 .E-3 ) GO TO 40

AG =( ( (0.0546875*G+0.3125)*G+0,625)*G+0.125)*G*GI*GI

GO TO 35

85
1 FAC = 2.0 * CON/{X*U)

ARG = FAC * (1. - GI * 0.125 * G * G)

IF( C.LT. 1 .E-3 ) GO TO 40

AG = -(0. 125 + 0. 0390625*G) * GI * GI * G-*^*3

GO TO 35

3 IF(G.GT.0.01 ) GO TO 5 FAC = CON * RI * RI * S

W2 = W *• W

ARG = FAC * (({1.375*W2 - 0.375)*G - W)*G + 1.)

IF ( G.LT.^ .E-4 ) GO TO 40

AG = { 3. 4921875*W2-2. 734375 ) *W2+0. 2421875

AG = (AG*G - (2.12 5*W2 - 1.125 )-^W) G ** 3

GO TO 35

5 IFCG.LT. 100. ) GO TO 8

FAC = BCON * SQRT(1. - W)/(X*U)

GI = l./G

ARG = FAC * ( ( ( 0.625-!^W-0. 125 ) ^GI-0.5 ) «-GI* ( 1 .+W ) + 1 . ) IF(GI .LT.3.E-4) GO TO 40

AG= ( (1. 52 34375* W- 0.4609375) *W-0.648 4 3 75 ) *W+0. 08593 75

AG = { AG*GI+ (( 0. 25-0. 9375*W ) *W+0. 1 875 ))"( 1 .+W ) *GI**3

GO TO 35

8 A = (1. + G*W)/SQRT(1. + G-^(2.*W + G))

GO TO 20

10 AG = ABStY)

BY = ABS( 1. - Y)

IF(AG.LT.l.E-lO) GO TO 15

GI = (X/Y) 2 IF ( GI .GT.2.E-4) GO TO 15

IF{Y.LT.O.) GO TO 13

FAC = CON * BY / (AG * U)

ARG = FAC * ( l.-0.125*( l.+Y)*(3.-Y)*GI

IF(GI.LT.l.E-6) GO TO 40

AG = ( (0.0546875*Y-0. 328125 )*Y + 0. 2421875 ) *(GI*( l.+Y) ) «-*2

GO TO 35

13 FAC = 2.0 * CON/(X*U)

ARG = FAC * ( 1 .-0. 125*GI*BY**2

IF(GI .LT.l.E-6 ) GO TO 40

AG = (( 2. -Y ) *0. 0390625*Y+0. 0859375 )*( GI*BY ) **2

GO TO 3 5 15 IF(BY.GT.3.E-3) GO TO 18

I F ( BY.LT. 1 .E-11 ) RETURN

FAC = CON * XS * BY / U

ARG = FAC * ( l. + ( ( 1 .-1 . 375^^X2 )*BY + l.)*DX2*RY)

IF(BY.LT.2.E-4) GO TO 40

AG = (1. - 2.125*X2)*BY*(BY*DX2) 2 GO TO 35

18A= (X2+Y* DX2 ) * R

20 IFCA.GE.l.) RETURN

ARG = BCON * SQRTd. - A)/(X*U)

GO TO 40

30 IF(X.FQ.O.) RETURN

ARG = CON * 2.0 «• XS / (X ^- U)

GO TO 40

35 ARG = ARG + FAC * AG

86
40 IF( ARG.LT.O. ) GO TO 60

IF(ARG.GT. 0.05) GO TO 5 0

PHI = -(0.5 - ARG*ARG/24. ) *ARG*ARG

RETURN

50 PHI = COS(ARG) - 1.0

RETURN

60 WRITE (61»1^U) ARG» FAC» K, X» U, Y, S

100 FORMAT (^^ ERROR NEGATIVE VALUE OF ARG IN PHI» ARG = *E17.9,

1 * FAC = *E17.9/* K = *I1»* X = *E17.9,* U = *E17.9,

2 * Y = *E17.9,* S = -^E17.9)

RETURN

70 WRITE (61*200) G

200 l-ORMAT (* PHI NOT ACCURATE ENOUGH G = *E17.9)

CALL EXIT

END

C

C

C

SUBROUTINE WEDDLE ( DX N, F» A, FO

)

C F IS THE FUNCTION TO BE INTEGRATED BY WEDDLES RULE

C FO IS THE VALUE OF THE FUNCTION TO BE INTEGRATED AT SOME STARTING

C POINT WHICH IS NOT INCLUDED IN THE INPUT ARRAY F DIMENSION F(N)

1 A = 0.0

2 K = N - 1 3 DO 15 I = 1»6

4 SUM = 0.0

5 DO 6 J = I » K » 6

6 SUM = SUM + F(J)

7 GO TO (8» 1^» 12» 10, 8, 14), I 8 A = A + 5.0 * SUM

9 GO TO 15

10 A = A + SUM

11 GO TO 15

12 A = A + 6.0 * SUM

13 GO TO 15 14A=A+2.0* SUM

15 CONTINUE

16 A = 0.3 * DX * ( A + FO + F(N))

17 RETURN

END

C

C

c FUNCTION SNZ(X)

C SNZ(X) CALCULATES RAABES INTEGRAL OVFR S I N ( X*Z ) / ( 1 . +Z*Z ) DZ

DIMENSION AA(6), BB(4), CC ( 4

)

DATA( ( AA( I ) , 1= 1,6) = -0.57721566, 0.99999193, -0.24991055,

1 0.05519968, -0.00976004, 0.00107857)

DATA ( ( BB( I ) I = 1,4) = 8.5733287401, 18.0590169730, 8,6347608925,

1 0.2677737343)

lATA( (CC( I ) 1= 1,4) = 9.5733223454, 2 5.6329561486,

1 21.0996530827, 3.9584969228)

SNZ = 0.

87
IF(X.FO.O.) RETURN

A = ABS(X

I F ( A.GT.40. ) GO TO 100

A2 = A * A

EMA = EXPF(-A)

IF(A,GT.0.2> GO TO 10

EEA = ( ( { (A2/72.+1. )*A2/42.+l. )*A2*0.05+1. )*A2/6.+l. )*A

GO TO 20

10 EEA = 0.5 * (EXPF(A) - EMA)

20 IF( A.GT.l. ) GO TO 30

SUM = AA( 1

Z = A

DO 25 J = 2»6

SUM = SUM + AA(J) * Z 25 Z = Z * A

FIT = SUM - LOGF(A)

SNZ = FFA * FIT

GO TO 40

30 SUM = ( ( ( A + BB ( 1 )

) *A + BB ( 2 ) ) *A + BB { 3 ) ) *A + BB ( 4

SUMM = { ( ( A+CC ( 1 )

) *A+CC ( 2 )) *A+CC ( 3 )

) *A+CC (4

FIT = SUM/(SUMM * A)

SNZ = 0.5 * FIT <- (1. - EXPF ( -2.-!^A ) ) 40 PROD = A

SUM = A

DO 45 J = 1*200

AJ = 2 * J + 1 PROD = PROD * A2/ ( AJ* ( AJ-1 . ) ) SUM = SUM + PROD/AJ

PT = PROD * l.F+10

IF(SUM.GT.PT ) GO TO 50

45 CONTINUE

50 PT = SNZ + SUM * EMA

SNZ = SIGNF(PT»X)

RETURN

100 PROD = l./X

A2 = PROD * PROD

SUM = PROD

DO 1 10 J = 1 » 100

AJ = 2 * J

PROD = PROD * A2 * AJ * (AJ - 1.)

SUM = SUM + PROD

IF{AJ.GF.A) GO TO 150

PT = ABS(PROD * l.E+10)

IF( ABS( SUM) .GT.PT) GO TO 150

110 CONTINUE

150 SNZ = SUM

RETURN

END

C

C

C

SUBROUTINE BESMOD ( X F I 0 F I 1 FKO » FK 1

C BESMOD CALCULATES THE MODIFIED BESSELFUNC T I ONS 10, II, KO AND Kl

C BY MEANS OF POLYNOMIAL APPROXIMATIONS AS GIVEN IN THE

88
NBS HANDBOOK OF MATHEMATICAL FUNCTIONS* PAGE 378

DIMENSION A(7)»B(9)»C(7),D(9) £{?), F(7)»G(7),H(7)

DATA ((A(I)»I = 1,7) = 0.0045813, 0.0360768, 0.2659732,

1 1.2067492, 3.0899424, 3.5156229, 1.0)

DATA {(B(I)»I = 1,9) - 0.00392377, -0.01647633, +0.02635537,

1 -0.02057706, 0.00916281, -0.00157565, 0.00225319, 0.01328592,

2 0.39894228)

DATA ({C(I)»I = 1,7) = 0.00032411, 0.00301532, 0.02658733,

1 0.15084934, 0.51498869, 0.87890594, 0,5)

DATA ((D(I),I = 1,9) = -0.00420059, 0.01787654* -0.02895312,

1 0.02282967, -0.01031555, 0.00163801* -0.00362018,-0.03988024,

2 0.39894228)

DATA ((E(I),I = 1,7) = 7.4E-6, 1.075E-4, 0.00262698* 0.0348859*

1 0.23069756* 0.4227842* -0.57721566)

DATA ((F{I)*I = 1*7) = 0.00053208* -0.00251540, 0.00587872,

1 -0.01062446, 0.02189568, -0.07832358, 1.25331414)

DATA ((G(I)*I = 1*7) = -4.686E-5* -0.00110404, -0.01919402*

1 -0.18156897* -0.67278579* 0.15443144* 1.0)

DATA ((H(I)*I = 1*7) = -0.00068245* 0.00325614* -0.00780353,

1 0.01504268, -0.0365562* 0.23498619* 1.25331414)

XS = 0.0

XE = 0.0

IF(X.LT.2. ) GO TO 10

XS = SORT(X)

XE = EXPF(X)

10 IF(X.GT.3.75 ) GO TO 20

Y = X/3.75

Y2 = Y * Y

FIO = ( ( ( ( ( Y2*A ( 1 )+A ( 2 )

) *Y2+A { 3 )

) *Y2 + A (4 )

) *Y2 + A( 5 )

) * Y2 + A ( 6 )) * Y2 + 1

.

FIl =((((( Y2*C ( 1 ) +C ( 2 ) )*Y2 + C(3) ) *Y2+C ( 4 )

) *Y2 + C ( 5 )

) *Y2+C ( 6 )

) *Y2 + 0. 5 FIl = FIl * X GO TO 40

20 Y = 3.75/X

FIO = 8(1)

FIl = D{ 1 ) DO 30 K = 2*9

FIO = FIO * Y + B(K)

30 FIl = FIl * Y + D(K)

XEX = XE/XS

FIO = FIO * XEX

FIl = FIl * XEX

40 IF(X.GT.2.) GO TO 60

Y = X/2.

Y2 = Y * Y

FKO = E ( 1 ) FKl = G( 1 ) DO 50 K = 2*7

FKO = FKO * Y2 + E(K)

50 FKl = FKl * Y2 + G(K)

XEX = LOGF(Y)

FKO = FKO - XEX * FIO

FKl = FKl/X + XEX * FIl

RETURN

60 Y = 2.0/X

89
FKO = F( 1 ) FKl = H{ 1 ) DO 70 K = 2»7

FKO = FKO ^ Y + F(K)

FKl = FKl * Y + H(K)

XEX = 1 ./ { XS * XE)

FKO = FKO * XEX.

FKl = FKl * XEX

RETURN

END

FUNCTION GLQUAD(F,A,B,L)

GAUSSIAN-LEGENDRE QUADRATURE

40 OR 80 NODES FOR L = 1» 2»

COMMON/GLQDAT/Xl (2) »W1(2) »X2

X5 ( 20 ) »W5 ( 20 ) XS { 40 ) we (

DATA

DATA

DATA

DATA

DATA

( ( XI ( I

( ( Wl ( I

( ( X2 ( I

( ( W2 ( I

( ( X3 ( I 1 .973906528b

DATA ( ( W3 ( I 1 ,0666713443

DATA ( ( X4{

1 .9931285992

2 .6360536807

DATA ( ( W4{ I 1 .0176140071

2 ,1181945320

DATA ( ( X5 (

1 ,9982377097

2 .9020988070

3 .6719566846

4 ,3419940908

DATA ( { W5 (

1 ,0045212771

2 ,0334601953

3 ,0574397691

4 ,0728865824

DATA ((X6(I

1 ,9995538227

2 ,9749091406

3 ,9132631026

4 .8169541387

5 .6896376443

6 .5361459209

7 .3623047535

8 .1747122918

DATA (IW6(I

1 .0011449500

2 .0086839453

I = 1^2) =

I = 1 ^2 ) =

I = 1^3) =

I = 1^3) =

I = 1^5) = .8650633667

I = 1^5) =

. 1494513491

I = 1»10)

.9639719273

.5108670020

I = 1 10

.0406014298

.1316886384

I = 1^20)

.9907262387

.8659595032

.6125538897

.2681521850

I = 1^20)

.0104982845

.0387821680

.0613062425

.0747231691

I = 1^40)

.9976498644

. 9654850890

.8966755794

.7938327175

.6608598990

. 5028041 1 19

. 3256643707

.1361640228

I = 1^40)

.0026635336

.0101617660

OF F FROM A

3^ 4^ 5^ 6

( 3 ) WP ( 3 ) »X3

•0 ) 8611363116

3478548451

9324695142

1713244924

6794095683

2 190863625

9122344283

3737060887

0626720483

1420961093

.9772599500

8246122308

5494671251

.1926975807

.0164210584

. 0438709082

.0648040135

.0761 103619

.9942275410

9545907663

.8787225677

.7695024201

6310757730

.4686966152

.2885280549

0974083984

,0041803131

,0116241141

TO R WITH 4

( 5 ) W3 ( 5 ) X4

.3399810436

.6521451549)

.6612093865

. 3607615730

.4333953941

.2692667193

.8391169718

.227785851

.0832767416

. 1491 729865

.9579168192

.7783056514

.4830758017

.1160840707

.0222458492

.0486958076

.0679120458

.0770398182

.9892913025

.9422427613

.85943 14067

. 7440002976

,6003306228

,4338753708

.2509523584

.0585044372

.0056909225

.0130687616

6»10^20»

( lO ) W4( 10

.2386191861 \ .4679139346

.1488743390

.2955242247

.7463319065.

.0 76 5265211:

.1019301198.

.1527533871

.9328128083-

,7273182552

,4137792044.

.0387724175

,0279370070<

.0532278470-

.07061 16474,

,0775059480

.9828485727,

.9284598772

.8388314736,

.7173651854.

.5686712681

.3983934059,

.2129945029,

.0195113833

,0071929048,

,0144935080,

90
3 .0158961836* . 4 .0225050902, . 5 .G28259816l» . 6 .0329419394, . 7 .0363737499, , 8 .038424993U, , TO = ( A + B) /2.

Tl = ( R - A) /2.

Y = 0.

GO TO (1,2,3,4,

1 DO 10 K = 1»2

10 Y=Y+W1 { K

) * ( F ( TO

GO TO 100

2 DO 20 K = 1»3

20 Y=Y+W2 (K)*( F( TO

GO TO 100

3 DO 30 K = 1»5

30 Y=Y+W3 ! K) * ( F ( TO

GO TO 100

4 DO 40 K = 1»10

40 Y=v+W4 ( K

) * ( F ( TO

GO TO 100

5 DO 50 K = 1»20

50 Y=Y+W5 ( K) * ( F ( TO

GO TO 100

6 DO 60 K = 1»40

60 Y = Y + W6(K)*-(F(T0

100 GLQUAD=Y*T1

RETURN

END

0172746521

0237318829,

0292883696 » 0337332150,

0368977146 » 0386617598 »

.0186268142 » .0249225358 , .0302723218 » .0344731205 » .0373654902

.0388396511 ,

5,6) L -Tl*Xl(K) )+F(TO+Tl*Xl (K)

-Tl*X2(K) )+F ( T0+T1*X2 (K)

-T1*X3 ( K ) ) +F { T0+T1*X3 ( K

)

-Tl*X4(K) )+F(T0 + Tl^^X4(K)

-Tl*X5 ( K ) ) +F ( T0 + T1^^X5 ( K )

-Tl*X6(K) ) +F ( T0+T1*X6 ( K

) .0199506109,

.0260752358

,

.0312101742

.0351605290

,

.0377763644,

.0389583960

.0212440261

,

.0271882275

.0321004987,

.0357943940,

.0381297113,

.0390178137)

i

91
APPENDIX B

THE LARGE TIME LIMIT OF THE THERMAL AVERAGE F(t)

In Eq. (IX. 31) the large time limit of the thermal average has

been given, which is of the form

This form has been obtained by most modern impact theories. The

additive constant B varies depending on what type of cutoff has been

usedo In the following we derive the different constants B for the

different cutoff procedures which have been used and compare them

with the numbers given in the literature.

F(t, n , n , T), differ essentially in three respects, namely by the

k e upper and lower limits of the p -integral and by the limits of the t'-

integral in Eq. (VIII. 4). Based on the completed collision assumption

(Baranger, 1962), the limits of the latter integral are usually extended

from -00 to + <». This approach, however, is not quite consistent with

the cutoff at the Debye length, which would rather require the integral

to go from -T to + T as done in this report (T is defined by Eq. (IX. 7)).

We therefore have to investigate the following integrals:

F(t)^ = At (B-^n (4C^)). (B.l)

The various methods to evaluate F ,the large time limit of

00

a.) (case a)

00 1

X mm (B.2)

92
and

b. ) for/rj PV^(t')dt' PV^(t') dt' ^/ PV^(t')dt' (case b)

1 du ==, u e I dx X cos I

J - "~b Z F = 2^n (a. D) V t /du =^ 00 e av

(B.3)

Xmm

where

Xmm. = p mm. / (a- D) . (B. 4)

The factor a = p /D is usually taken to be one and has in some max

papers (e. g. , Griem et al. 1962) been varied to 1. 1. or to 0. 606

as proposed by W. -R. Chappell, J. Cooper and E. Smith, 1969.

As a lower cutoff we consider in particular the three cases of

p . = 0, p . - K and p . = 3(nq-n'q')>^ 2CD/u by setting mm mm mm ^ > i 7 &

(B.5) mm a u

In the following we will set a=l. In order to recover the dependence

on the upper cutoff parameter a, we only have to replace in all the

following relations C by C/a.

First of all, one realizes that with x . < 1 the lower limit on the mm

u- integral is given by

u = b- 2C . (B. 6) o

Hence, we have to evaluate the following two integrals

2

=y du Je'^^^j dxx l^cos

^

^-x^'j-lj (B. 7)

o/ u

,a

u

o

93
and

dx X

u u /u

cos (B. 8)

o o- The second integral can be simplified after a change of variables and

a partial integration to

du (B.9)

After expanding the cosine and another change of variables we have

?. 2 00 k

^ r-n /?r.\2k / ^ o dz (B. 10) 4 Z-^ (2k)! U 1 / k=l ^ o / A

which can be expressed in terms of exponential integrals

2k

b C

I = - E ,2 ""o (-1) /2C\

k=2 ^ o /

E (u ) . (B.ll) k o

With the lower cutoff parameters stated above (b < 1) and typical

densities and temperatures of interest one usually has u < 0. 1 (see

2 2 ° Eq. (IX. 21)). Since for k> 2 E(u ) = l/(k-l ) + 0(u ) one obtains

2 ^ ^

°

to lowest order in u

o 2k-2

(B. 12)

.2 r

k-l

1 (2k)! (2k-l)! (2k-2)(2k-2)!

94
+ 2 fr)f""(f)(f)|-H(-^)-(t)|]-

which yields

3{y-1) + (4C^) + 2K ( ^ + 0 (u") | (B. 13)

where K is defined as

, \ 1 -cos z , sin z ,

. , .X K (z) = — + Ci (z) (B. 14)

L. Z

z

and Ci is the cosine integral. Eq. (B. 13) was obtained already by Shen

and Cooper, 1969. Their constant A is identical with our constant 2C.

The other integral I of Eq. (B. 7) one can obtain by evaluating

Al = I du u e I dx X cos I '%/1-x l-cos

J J L V""' /

so that

= + Al . (B. 16)

If we again expand the cosine functions, aI can be given by

00 y, k-1 . no _ 1 . \-> (2C) '^"^ k! (-l)J

/ , 3-2k -u / l-2j

\ V" (B.17) —— e

95
2k-2

S

I /2c

(2k)!

k=2 ^

2C\ E (u )-E ( k-1 o k o

I k-1

dz

2k-2

^

4

2 k! / 2C\ (-l)-^

""o (2k)! I u I Z^j!(k-j)!|

k.3 \ 0/

(2j-2) k-1 o k-j o (B.17)

which gives us to lowest order in u

C

1 + 0(u ) o (B. 18)

This means that for the same lower cutoff case a. ) and b. ) as defined

in the Eqs, (B. 2) and (B. 3) differ only by a constant 1 in their additive

constants B. As a result we have

F = - (-l-(nq-n'q')-^l n t

00 12 m I e B-;,n (4C ) (B. 19)

'1

where the constant B for the different cutoff parameters is compiled in

the following Table 2.

Table 2. The constant B for different cutoff parameters.

+ T

/ PV^(t') dt'

- T

+00

/ PV (f) dt'

J

"

-00

p . = 0 mm 0. 27 1.27

p . = ^ mm 0. 23 - 0. 27 1. 23 1.27

p . = 3n ^ mm k

-1. 66 - 0. 66

96
In order to compare our results with the numbers given in the literature we rewrite Eq. (B. 19) as

F = - A-t- ^B -Y-; n(y . )| (B. 20)

00 o mm

where y . as introduced by Griem, Kolb and Shen (Griem et al. , 1959) mm

is given by

4TTn / 2\ ^

^ I 2 \^

y . =

e

I/ ehn

, \

I = -r

2

1-77

I n

TtA12, 4C . B.21 mm 3 m V kT / ^ \ 3(^*1"^ 4 )/

Consequently, B and B are related by the following relation

B = B +Y+^ (B. 22)

Comparing Eqs. (B. 19) and (B. 20) one notices that for a particular

line the value of the square bracket as derived here depends on the

quantum number n, for that particular state. This is also true for the k

paper of Shen and Cooper, 1969, who consider our case (b) with infinite

limits on the t' -integral. Otherwise the constants given in the literature are independent of n because the lower cutoff parameter is usually

2

based on an average Stark splitting. If we set nq-n'q' = n /2, which

corresponds approximately to the average Stark splitting and also

gives the results for the Stark shifted component of Lyman-a, we have

B = B- 0. 64 . (B. 23)

o

2 This yields directly for n^ - n /Z the B values corre'sponding to the

B values in Table 2.

97
The following constants B have been given in the literature

o Griem, Kolb, Shen, 1959; (Eq. 29) : B = 0 o Griem, Kolb, Shen, 1962; (Eq. 2) : B = 1.0

Griem, 1965

Kepple, Griem, 1968

(neglecting quadrupole term) : B = 0. 58

o

Shen, Cooper, 1969; B = 0. 58

o Recently the time development operator (S-matrix) has been evaluated

for Lyman-a including time ordering by solving the differential equations

for the S-matrix elements (Bacon, 1969)., Again the square bracket

depends on (nq-n'q') and the average value B =1,1 considering only

the dipole term. It should, however, be stressed that one should not

overinterpret these numbers because within the classical path approximation there is always some uncertainty about the "correct" constant

B because of the ambiguous lower cutoff. This is due to the fact

that the classical path approximation breaks down roughly for p< ^

(for details see Paper I), For most cases this has no significant effect

for the Stark broadening of hydrogen because the dynamic broadening is

primarily due to weak collisions. More details are given with the

discussion at the end of Chap. XII. The situation is quite different for

the broadening of ionized lines where strong collisions are very important and where the uncertainty of the classical path approximation

accounts for part of the still existing discrepancies between theory

and experiment, which are large compared with the Stark broadening

of hydrogen.

So far we have considered F , which is the basis of the un00 modified impact theory. In order to extend the range of validity beyond

the plasma frequency the modified impact theory introduces the Lewis

cutoff by considering only those collisions for which the duration of a

collision, which is typically p/v, is smaller than the time of interest

being typically 1/aud. For this reason the modified impact theory

introduces an upper cutoff

98
Pmax= ^^^<^'^a>^^ 24)

or

It should be noted at this stage that in the following relations we not

only have to replace C by C/a but also AU)_^ by a Auu_, in order to obtain

R R

the dependence on the upper cutoff parameter a. Considering the usually applied case (b) (Eq. (B. 3)) we have to evaluate the following

integral for ^{u > I R

duu^-y dxx[cos(^-l] ,B.25,

u AO) T3 ^ / o R o/u

where the lower limit on the u-integral is determined by the condition

u /u< 1/auu„. After a change of variables and a partial integration

o R

one obtains similar to Eq. (B. 9)

^ZCAm^^

du (B.26)

1 r -n' (^^'">r\

:r I ue cos I 1 - 2 Am

R u auUt, o R

Expanding the cosine again and performing another change of variables

the result is

2 2

2 « 2k °° -u Aw^z

4

2 « 2k °° -u Aw^z

1 This gives us then

2 2k

(B. 28)

99
Evaluating the exponential integrals E for small arguments only we

'nally have

L 2 -^i(v^r) + ^1 ("o' + I

^2 3 = C j-.n AuOj^ + i (B.29)

''^^'ch gives us for ^w ^ > 1 the log-dependence of the $ -matrix elements R ab

in the modified impact theory, A more appropriate way for applying the

Lewis cutoff, which avoids the discontinuity at Aoj = 1 , is to take as an upper cutoff

p = MIN (D, v/auu) max

or p /D=MIN(1, u/aoj^) max R(B. 30)

which for case b leads to the following integral.

f \ 3 -uY / r /2c\ 1

I =1 du u e f dx X cos I 1 -1 Lj

J I \:.nj

J

./ u^e- Y ax . [cos M -:1

A(«o u /u R o

(B.31)

100
These integrals are identical to

duu e J

dxx

V"" (B.32)

0= 2 -V^^R

-J

duu^-^

I

dxx j^cos

^-^^

-ij .

and after a partial integration we have

00

= —— I du u e L 2

AO)

cos

*

' o R

{B. 33)

-u^r /2c\ ,1 + i| du ue cos "1

^'•"r b Expanding the cosine functions, I can be evaluated in a similar

manner as above with the result

_C_

2 E^{A^j^)-2E^(u^Ac«j,)

+ 2C

(2k)!

k

k=2

2k-2

(B. 34)

E (u Auu„) 2k-l o R

00 2k-2

101
For fSW-^ 1 with 2C«1 and b'< 1 the latter result may be simplified

to give

.2 r

+ I (B. 35)

b b which for fsa^ ~* 0 reduces down to I =1 and which has no discontinuity

at = 1

to

Furthermore, we see that for Aw R

00 Eq. (B. 34) goes over

2

u

L, Aw R

2AaJ

o R

e cos

R

/2C\ e /2C \

(B. 36)

which does not lead to the static limit.

Similar results can be obtained for case a, which are not

included because they are no longer required. The derivation for case

b has been included, in order to obtain consistent relations which allow

a comparison with the calculations dons in this paper. The results for

case b as given here differ slightly from the results in the literature

which also vary from paper to paper depending on the average matrix

elements used and on what lower cutoff and average velocity has been

applied.

102
APPENDIX C

PROGRAM FOR CALCULATING THE FINAL LINE PROFILE I(A(ju)

In the following a complete listing of the program is given which

j was used to calculate the final line profile I(Aa)) on the basis of the unified

! theory for the case of no lower state interactions.

1. The Fourier transform of the thermal average

The complex function FOUTR calculates essentially the Fourier

J transform of the thermal average as defined by

FOUTR =iAtUpi (A(JUp,B, n, q q ) (C.l) K u K- be

i

' where

DOM= Auu„ = (Auu-Aa).( n, q. ) B)/uu . (C,2) R X b p

I It uses the Eqs. (X.17) and (X. 22) for calculating i (A(i)„) and i {f\(xi^)

I 1 R 2 R

respectively. The required Bessel functions J , J, , Y and Y, are

o 1 o 1 evaluated by the subroutine BSJYOl. For large and small arguments

i ^^R ^^^^^ relations are replaced by their asymptotic expansions (X. 18),

(X. 23) and (X. 26). The specifying constants p , p , b , a and b

(PI, P2, Bl, A2, B2) are set in the function AIIM and are calculated

once for all the Stark components in the main program STBRHY. The

function FOUTR can be replaced by another short function FOUTR listed

at the end of the program, which then makes the program calculate

line profiles according to the modified impact theory.

103
2. Calculation of I(AijJ, p)

The function AIIM calculates I(AtJO, B) as defined by Eq. (XII. 6).

It establishes first of all the matrix of the jb-operator according to

Eq. (XII. 7) and calculates the array

AMATR (NBN„ NAN) = 5 Auu„-(nq^mU(A(W )Inqm) (C. 3)

q q R b op a

3, D

for the m values 0 and ±1. The required 3j- symbols are calculated

once in the main program and their values are stored in three different

arrays according to the following definitions

n-1 n-1 Y

SAR (NLA, NQC) = ( Zl^ + 1 )^

and

m

I

m -q m +q

c c c c •m

(C.4)

i

where NLA = l and NQC = q

SJQL (MCT, NBN, NLA) =

/ n-J_ n-1 \

\

m-q m+q

a a •m

(C.5)

/

SSJJ (MCT, NBN, NAN)

q_+q

n+m-1- —

(-1)

'b n-1

2

m-q

n-1

2

m+q

m

/ n-1 n-_l

2 2

m-q^ m+q^

•m

( C. 6)

104
where MCT = m and NAN and NBN are nurabering indices specifying

the q and . The matrix AMATR is then inverted by the subroutine

a b

' CGAUSSEL, which is able to solve systems of linear, complex equations

by Gaussian elimination. Multiplying the inverted matrix by the 3jsymbols according to Eq. (XII. 6) yields finally I{Aw, 6).

3. The final line profile I(Auu)

The main program STBRHY calculates the final line profile by

performing the ion field average according to Eq. (II. 1). It first of

all reads in the ion microfield distribution function for 0 < B ^ 30 in

steps of 0. 1, which has been calculated in a separate program for the

particular shielding parameter r /D. It then reads in the density n , o e the temperature T, the upper principal quantum number n, the wavelength X , the average value of the static ion fields 8 , the initial

value AUU) the logarithmic stepwidth, the total number of points, a parameter which specifies the number of ion field integration points

and finally 6 numbers, which specify the G -function and hence i (aoj^)

for all Stark components and which may in practically all cases be set

to zero. These 6 numbers are obtained from the thermal average

described in Appendix A,

The program then calculates the constants p^ , p^, b^ , a^ and b^

for all Stark components and stores them in the array FPAR. p^ is

calculated according to Eq. (X. 5). p^ is determined on the basis of

Eq. (B. 19) where the constant B is given by

B = 0.27- 2K(-^) . (C.7)

b

105
The K-£unction is defined in Eq. (B. 14) and is calculated for a lower

cutoff p . = -X + -7 n a by setting

^ mm O

The necessary cosine integral is calculated by the function COSINT.

As a next step the main program evaluates all the required

3j-symbols by means of the function S3 J, which in turn uses the function

FCTRL, to calculate all the necessary factorials. The numbers are

stored in arrays according to the definitions in (C. 4) , (C. 5) and (C. 6).

In performing the ion field average the microfield distribution

function is calculated by the function WFLD, which uses a 5 point

interpolation on the values read in initially for B< 30 and otherwise

uses the asymptotic expansion given by Hooper, 1968b. As a function

of Auu and |3, which determine the shape of the ion field integrand, the

ion field integral is subdivided in intervals, which are integrated

separately by means of Weddle's rule (subroutine WEDDLE) using a convenient change of variables in every interval. Furthermore, the

program calculates the asymptotic wing expansion according to

Eq. ( XI. 12) and the unified theory for 8=0 and 3 = B^^ performing

not the ion field integral. All three values are normalized with

-5/2 . respect to the asymptotic -wing.

106
PROGRAM STBRHY

C

C PROGRAM FOR CALCULATING THE STARKBROADEN I NG OF HYDROGEN ON THE

C BASIS OF THE UNIFIED THEORY FOR NO LOWER STATE INTERACTION

DIMENSION FF(1000)» PFAC{6)» SSJ(20), STRONG(20)

COMPLEX FOUTR

C0MM0N/FDAT/P1»P2»B1 »A2 »B2 »PPFF

COMMON/ PSTR/NNN»NM1 »BET tFPAR ( 6 »20

)

COMMON/PSJD/SSJJ( 2»20»20) . SJQL ( 2 »20 » 20 ) » SAR(20»20)

COMMON/ PFW/F I ELD ( 301

)

FIELD(l) = 0.0

READ 100» (FIELD(I)» I = 2,301)

100 FORMAT (6E12.4)

120 READ 150,DEN»TEMP,NNN,ALAM»BAV»GIN»DGG»NTOT.NFAC» (PFAC{ I ) »I = 1 »6)

150 FORMAT (2E10.2, 15, 4F10,2, 2I5/6F10.5)

IF (EOF, 60) 577, 170

170 PRINT 180, DEN, TEMP, NNN» ALAM

180 FORMAT (1H1»* DENSITY = *E12,A* TEMPERATURE = *E12.4,

1 * QUANTUMNUMBER =*I2* WAVELENGTH =*F8.2* ANGSTROM*//

2 13X,»Pl*18X,*P2*18X,*Bl*18X,*A2*18X,*B2*17X,*STROWG*/

)

IF (NNN.LE.20) GO TO 200

PRINT 190

190 FORMAT (» PROGRAM NOT EXECUTED BECAUSE N IS LARGER THAN 20*)

577 CALL EXIT

200 SDEN = SQRT{DEN)

FAC = 2064.936 * TEMP * SQRT ( TEMP/DEN

)

NMl = NNN - 1 NEVODD = M0D(NNN,2)

AN = NNN

AN2 = AN * AN

ANIM = 0.5 * (AN - 1.)

CFAC = 4.5645E-7 * AN * SDEN/TEMP

DEBROG = 2. 1027E-6/SQRT ( TEMP)

RMIN = DEBROG + AN*AN*7 • 9376E-9

DO 230 K = 1, NNN

AK = K — 1 230 SSJ(K)=(aN2+( (-1. )**M0D(NMi+K,2) ) * ( AN2-2

.

*AK*AK ) ) / ( 2

.

*AN* ( AN2-1 • )

)

BET = 5.6558E-5 * AN * DEN**(l./6.)

C

C ARRAY FOR G-FUNCTION CONSTANTS

ASY = 0.0

DO 270 < = 1»NM1

QC = K

C = CFAC * QC

PI = -1.671086 * FAC * C * SQRT(C)

BS = 3. * AN * QC * DEBROG/RMIN

STRONG (K) = 0.269-2.* ( ( ( l.-COSF(BS)

)

/BS+S INF ( BS ) ) /BS'COS iNT ( BS )

)

PPFF = -1.128379 FAC * C * C

P2 = PPFF * (STRONG(K) - 2 . *LOGF ( 2 . *C )

)

FPAR(1,K) = PI

FPAR{2»K) = P2

FPAR{3»K) = 0.5 * (P2/P1)**2

FIN = LOGF(AN*QC)

FPAR(4,K) = P2 * ( (PFAC(3)*FIN + PFAC(2)) * FIN + PFAC(l))

107
FPAR(5»K) = (PFAC(6)*FIN + PFAC(5))»FIN + PFAC(4)

FPAR(6»K) = PPFF

270 ASY = ASY + 2. * FPAR(1,KJ * SSJ(K+l)

PRINT 2A0 » ( { FPAR( K» I) = 1 5 ) »STRONG ( I ) » I = l.NMl)

240 FORMAT (6E20.4)

PRINT 280» FAC, CFAC , BET* ASY, DEBROG» NFAC

?.8U FORMAT(/* FAC =*E12.4,* CFAC =*E12.4,* BET =*E12.4,* ASY =*E12.4

1 * DEBR06 =*E12.4,* INTEGRAT lONFACTOR =*I2//

2 5X*DOM*8X*DLAMi^8X*ITOT*6X*lHOLTS*8X*ASY*10X*WING*7X*WHOLTS*7X

3 *WSTAT*7X*WWOO*8X*WWBB*8X*WWNG*/

AOLFAC = 4.23538E-15 ^ SDEN * ALAM * ALAM

C

C 3JSYMB0L--ARRAYS SS J J ( MCT , NBN ,NAN ) AND S JQL ( McT .NBN ,NLa ) DO 650 MCT =1,2

AMA = MCT - 1 NLIM = NNN + 1 - MCT

NOB = -NLIM - 1 DO 450 NBN = 1,NLIM

NQB = NQB + 2

QB = NQB

FMBl = (AMA - QB) * 0.5

FMB2 = (AMA + QB) * 0.5

DO 325 NLA = 1»NNN

ALA = NLA - 1 325 SJQL(MCT,NBN,NLA) = S3 J ( AN IM AN IM » ALA » FMB 1 , FMB2 . -AMA

450 CONTINUE

DO 375 NBN = 1,NLIM

DO 350 NAN = 1,NLIM

AABB = (-1. )^^*M0D(NAN+NBN.2 ) * S JQL ( MCT , NBN , 2 ) * SJQL ( MCT .NaN ,2 ) 350 SSJJ(MCT,NBN,NAN) = AABB

375 CONTINUE

650 CONTINUE

C

C 3JSYMB0L-ARRAY SAR(NLA»NQC)

DO 780 NQC = l-NMl

QC = NQC

DO 720 NLA = 1,NNN

ALA = NLA - 1 FBB =0,

DO 680 NMC = 1»NLA

IF(NEV0DD.NE.M0D(NMC+NQC»2) ) GO TO 680

AMC = NMC - 1 FCF = 2.

IF{NMC.EQ.l) FCF = 1,

FMCl = 0.5 * (AMC - QC)

FMC2 = 0.5 * (AMC + QC

FBB = FBB + FCF * ( S3J ( ANIM, ANlM ALA tEMCl »FMC2 ,-AMC) ) **2

680 CONTINUE

720 SAR(NLA,NQC) = FBB * (2.*ALA + !•)

780 CONTINUE

C

C CALCULATION OF THE lONFIELD INTEGRAL

NN6 = 6 * NFAC

ANN6 = NN6

108
N12 = 12 * NFAC

AN12 = N12

N30 = 30 * NFAC

AN30 = N30

G = GIN - DGG

DO 950 MM = 1»NT0T

G = G + DGG

DOM = 10. ** G

DLAM = ADLFAC * DOM

WING52 = -0.2992067103 * AS Y/ ( SORT ( DOM ) * DOM * DOM)

FHOLTS = 0.

AWING = 0.

DO 815 NQC = 1,NM1

PI = FPAR(1»NQC)

P2 = FPAR(2»NQC)

Bl = FPAR(3»NQC)

A2 = FPAR(4»NQC)

B2 = FPAR(5»NQC)

PPFF = FPAR(6»NQC)

AWING = AWING + SS J ( NQC + 1

) *A I MAG ( FOUTR ( DOM ) ) *2 . / ( DOM*DOM

)

QC = NQC

BETFAC = BET * QC

BCRIT = DOM/BETFAC

815 FHOLTS = FHOLTS + SSJ ( NQC + 1) * WFLD(BCRIT) / BETFAC

AIRES = 0.

IF(DOM.GT. (-3.*P2 ) ) GO TO SAO

BCRIT = (DOM - P2)/((AN - 1 . ) *BET

)

DB = BCRIT/ANN6

B = 0.

DO 820 J = 1»NN6

B = B + DB

820 FF(J) = AIIM(D0M,B) * WFLD(B)

CALL WEDDLE ( DB NN6 FF . A 1 1 I » 0 . ) AIRES = AI II

DY = l./(BCRIT*AN30)

Y = 0.

DO 830 J = 1»N30

Y = Y + DY

B = l./Y

830 FF(J) = B » B * AIIM(D0M»B) * WFLD(B)

CALL WEDDLE ( DY N30 FF » A I I I » 0 • ) AIRES = AIRES + AIII

GO TO 980

840 BCRCR = DOM/BET

EPSPS = -P2/BET

DO 957 NQ = 1»NM1

ANQ = NQ

BCR = BCRCR/ANQ

EPS = EPSPS/ANQ

IF(NQ,EQ.l) GO TO 907

SLl = l./(GAM - BCR)

GO TO 908

907 SLl = 0.

908 SL2 = l./EPS

109
SL3 = l./(BCR + EPS)

SL^ = 1,/(BCR - EPS)

GAM = 0.5 * (SCR - EPS + (BCRCR + EPSPS)/(ANQ + !•))

IF(NQ.EQ.NM1 ) GAM = 0.5 * (BCR - EPS^

SL5 = l./(BCR - GAM)

CRIT = SL2 - SL5

Y = SLl

IF(NQ.EQ.l) GO TO 913

B = BCR + l./Y

FA = AIIM(DQM,B) * WFLD(B)/(Y * Y)

GO TO 914

913 FA = 0.

914 DY = (SL2 - SL1)/ANN6

DO 917 J = 1»NN6

Y = Y + DY

Yl = l./Y

B = BCR + Yl

917 FF(J) = Yl * Yl * AIIM(DOM»B) * WFLD(B)

CALL WEDDLE ( DY ^NNe »FF » A 1 1 I » FA

)

AIRES = AIRES + AJII

Y = SL3

B = l./Y

FA = B * B AIIM(DOM»B) * WFLD(B)

DY = (SL4 - SL3)/ANN6

DO 927 J = 1»NN6

Y = Y + DY

B = l./Y

927 FF(J) = B * B » AIIM(D0M,B) * WFLD(B)

CALL WEDDLE ( DY .NN6 . FF » A III » FA

)

AIRES = AIRES + AI I I IF(CRIT.LE.O. ) GO TO 977

Y = SL5

B = BCR - l./Y

FA = AIIM(D0M»8) * WFLD{B)/(Y * Y)

DY = CRIT/AN12

DO 937 J = 1»N12

Y = Y + DY

Yl = l./Y

B = BCR - Yl

937 FF(J) = Yl * Yl * AIIM(DOM»B) * WFLD(B)

CALL WEDDLE ( DY »N 12 .FF » A 1 1 I » FA

)

957 AIRES = AIRES + AI I I IF(GAM.LT. 5.) GO TO 968

Y = l./GAM

FA = GAM * GAM * A I I M ( DOM » GAM ) * WFLD(GAM)

DY = (0.2 - Y)/AN12

DO 967 J = 1»N12

Y = Y + DY

B = l./Y

967 FF(J) = B * B * AIIM(D0M»B) * WFLD(B^

CALL WEDDLE ( DY N12 » FF » A I I I , FA

)

AIRES = AIRES + All I SL4 = 0.2

GO TO Q"77

110
966 SL4 = 1,/GAM

977 B = 0.

[j DB = l./(SL4 * AN30)

' DO 947 J = 1»N30

B = B + DB

947 FF(J) = AIIM{DOM»B) * WFLD(B)

CALL WEDDLE ( DB NSO FF » A I I I 0 .

)

AIRES = AIRES + AI I I 980 WING = AIRES/WING52

WISTAT = AlRES/FHOLTS

WINHOL = FHOLTS/WING52

WWOO = (AIIM(DOM» 0.) + FHOLTS ) / W I NG52

WWBB = ( AIIM( DOM»BAV) + FHOL TS ) / W I NG52

WWNG = (AWING + FHOLTS ) /WING52

950 PRINT 978»D0M,DLAM,AIRES»FH0LTS»WING52»WING»WINH0L»WISTAT,

2 WWOO* WWBB» WWNG

978 FORMAT (11E12.4)

GO TO 120

END

C

FUNCTION AIIM(DOM»B)

I C

'C CALCULATION OF I(DOM»B)

COMPLEX DFTR(2U) » AMATR ( 20 40

)

»F0UTR,aREF

COMMON/FDAT/Pl ,P2 ,31 »A2 B2 ,PPFF

C0MM0N/PSTR/NNN,NM1 ,BET»FPAR(6,20)

C0MM0N/PSJD/S5JJ( 2»20,20) S JQL ( 2 20 20 ) SAR(20»20)

AI IM = 0.

DO 800 MCT = 1,2

AMCT = MCT

NLIM = NNN + 1 - MCT

NL22 = 2 * NLIM

NQB = -NLIM - 1 DO 750 NBN = 1,NLIM

NQB = NQB + 2

QB = NQB

DOMRB = DOM - BET » QB » B

DO 220 NQC = 1,NM1

PI = FPAR(1,NQC)

P2 = FPAR(2»NQC)

Bl = FPAR(3,NQC)

A2 = FPAR(4,NQC)

B2 = FPAR(5»NQC)

PPFF = FPAR(6,NQC)

220 DFTR(NQC) = F0UTR(D0MRBJ

DO 700 NAN = 1»NLIM

AREF = (0.,0,)

DO 600 NQC = 1,NM1

FAA = 0.0

DO 500 NLA = 1,NNN

500 FAA = FAA + S JQL ( MCT ,NAN ,NLA ) * S JQL ( MCT , NBN , NLA ) * SAR(NLA,NQC)

600 AREF = AREF + FAA * DFTR(NQC)

AMATR{NBN,NAN + NLIM) = (0,,0.)

700 AMATR(NBN,NAN) = 6.2831853072 * (

( - 1 • ) **MOD ( NAN+NBN , 2 ) ) * AREF

111
AMATR(NBN»NBN+NLIM) = (1,»0.)

750 AMATR{NBN»NBN) = AMATR ( NBN ,NBN ) -i- DOMRB

C MATRIX INVERSION

CALL CGAUSSEL( AMATR » 20 NLIMjNLZZ^NRANK)

DO 795 NBN = 1»NLIM

DO 793 NAN = 1»NLIM

793 AIIM = AIIM + AMCT*SSJJ ( MCT »NBN tNAN ) *A IMAG ( AMATR ( NBN »NAN+NL IM )

)

795 CONTINUE

800 CONTINUE

AIIM =-AI IM * 0.3183099

RETURN

END

C

C THE FOLLOWING FUNCTION FOUTR MAY BE REPLACED BY THE FUNCTION FOUTf

C AT THE END OF THE LISTING TO OBTAIN THE MODIFIED IMPACTTHEORY

FUNCTION FOUTR (DOM)

C

C FOURIERTRANSFORM OF THERMAL AVERAGE FOR UNIFIED THEORY

COMPLEX FOUTR

COMMON/FDAT/Pl ,P2 ,B1 ,A2»B2 »PPFF

ARG = ABSF(DOM)

Z = Bl * ARG

IF (Z.LE. 0.001) GO TO 600

IF (Z.LE. 40.) GO TO 300

FACl = -0.2992067103 * P 1 /( SORT ( ARG ) * ARG * ARG)

CC = FACl * ((1. - 1.3125/Z)*0.625/Z + 1.)

SS = FACl * ((-1.- 1.3125/Z)*0.625/Z + 1.)

GO TO 500

300 CALL BSJYOl (Z» AJO» Y0» AJl» Yl)

FACl = Y1/{2.*Z) + AJl - YO

FAC2 = AJO + Yl - AJ1/(2.*Z)

CINE = COSF(Z)

SINE = SINF(Z)

CC = P2 * Bl * Bl * (CINE * FACl + SINE * FAC2)

SS = P2 * Bl * Bl * (CINE * FAC2 - SINE * FACl^

IF (A2.EQ.0.) GO TO 500

Z = B2 * ARG

IF (Z.GT. 10. ) GO TO 400

CALL BSJYOl (Z» AJO» Y0» AJl» Yl)

FACl = ( ( AJ1-Y0)*16.*Z-36.*AJ0-28.*Y1 )^Z+15.*YO-3.*AJl

FAC2 = ( ( AJO+Yl )*16.*Z-36.*Y0+28.*AJ1 )*Z-15.*AJ0-3.*Y1

CINE = COSF(Z)

SINE = SINF(Z)

CC = CC + A2*B2*(CINE * FACl + SINE * FAC2)/6.

SS = SS + A2*B2*(CINE * FAC2 - SINE * FACl)/6.

GO TO 500

400 FACl = 0.1322319336 * A2 * B2 * Z**{-3.5)

CC = CC + FACl * (1. - (3.9375/Z + l.)*4.375/Z)

SS = SS - FACl * (1. - (3.9375/Z - 1,)*4.375/Z)

500 IF (DOM.LT.O.) SS = -SS

FOUTR = ARG * ARG * CMPLX ( -SS»CC

RETURN

600 SS = (P2*B1 - A2) * DOM

FOUTR = 0.3183099 * CMPLX ( SS »-P2

112
RETURN

END

C

C

c

123

4

56

7

8

9

10

11

12

13

14

15

16

17

SUBROUTINE WEDDLE ( DX N» F» A» FO

)

INTEGRATION SUBROUTINE

DIMENSION F{N)

A = 0.0

K = N - 1 DO 15 I = 1» 6

SUM = 0#0

DO 6 J = I » K » 6

SUM = SUM + F(J)

GO TO (8» lU, 12» 10. 8,

5.0 * 14) I A = A + GO TO 15

A = A

SUM

+ SUM

200

70

74

75

GO TO 15

A = A + 6.0 * SUM

GO TO 15

A = A + 2.0 * SUM

CONTINUE A=0.3*DX* <A+FO+F(N))

RETURN

END

FUNCTION WFLD(B)

CALCULATION OF THE ION MICROFIELD DISTRIBUTION FUNCTION USING A

5P0INT INTERPOLATION FOR THE DATA READ INTO THE MAINPROGRAM

COMMON/PFW/FIELD( 301

)

WFLD = 0.0

IF (B.LE.30»0) GO TO 200

SBS = l./(B * SQRT(B)

)

WFLD = ((21.6 * SBS + 7.639) * SBS + 1.496) * SBS/B

RETURN

IF (B.LE.0.0) RETURN

J = (B + 0.2 ) * 10.0

L = J - 1 IF (J.GT.2) L = J - 2

IF (J.GT.3) L = J - 3

IF (J.GT.300) L = 297

LLL = L + 4

DO 75 K = L»LLL

AK = K - 1 TERM = 1.0

DO 74 M = L»LLL

IF (K.EQ.M) GO TO 74

AM = M - 1 TERM = TERM * (10. *B - AM)/(AK - AM)

CONTINUE

TERM = TERM * FIELD{K)

WFLD = WFLD + TERM

RETURN

END
c FUNCTION S3J (FJl, FJ2» FJ3» FMl » FM2 » FM3 ) C

C CALCULATION OF 3J-SYMB0L

S3J=0#0

IF(ABS(FM1 + FM2 + FM3 ) .GT. 0.001) GO TO 153

FM3=FM1+FM2

A=FJ2+FJ3+FMl+.O0 5 B=FJl-FMl+.005

C=-FjH-Fj2+FJ3+.005

D=FJ3+FM3+«005

E=FJl-FJ2-FM3+«005

F=FJl-FJ2+FJ3+.005

G=FJl+FJ2-FJ3+.005

H=FJ1+FJ2+FJ3+1 .0+.005

E2=FCTRL( B)*FCTRL(FJ1+FM1 )*FCTRL( FJ2-FM2 ) *FCTRL( FJ2+FM2)

IF (E2) 153»153»145

145 E1=(FCTRL(C)*FCTRL(F)/FCTRL(H) ) *FCTRL ( G ) *FCTRL ( D ) *FCTRL ( F J3-FM3

IF (El) 153» 153» 150

150 E1=SQRT(E2)/SQRT(E1)

I1=XMAX1F(0.0» -E+0.01)

I2=XMIN1F{A» CD)

IF (12-11) 153, 151f 151

151 DO 152 I=IliI2

FI = I E2=FCTRL(FI J *FCTRL(C-FI )*FCTRL( D-FI ) / FCTRL ( A-F I

152 S3J = S3J+( n-l,0)**XMODF( I»2) )/E2)*FCTRL(B+ FI )/FCTRL(E+FI

U=ABS(FJ1+FM2+FM3 ) +0.001

S3J=S3J*( (-1.0)**XMODF(XFIXF(U) »2) )/El

FM3=-FM3

153 RETURN

END

C

FUNCTION FCTRLCA)

C

C CALCULATION OF FACTORIALS REQUIRED BY FUNCTION S3J

DIMENSION FCTI (20)

C

DATA ( (FCTI ( I ) . I = l»20) = 1 . 0 , 2 . 0 6 . 0 24 . 0 1 20 . 0 720 . 0 » 5040 . 0

1 40320.0 » 362 880.0 » 362 8 800.0 » 399 16800 •0*47900 1600.0

2 622 7020 800 .0» 87 178291200.0 , 130 7674368000.0

,

3 2.0922789888E13, 3 . 55687428096E 14 , 6 . 402 373705 72 8E 1 5 4 1.2164510040883 E17» 2.432 9020081766 E18)

C

IF(A) 50»60»70

50 IF( A.GE. (-0.1 ) ) GO TO 60

FCTRL = 0.0

RETURN

60 FCTRL = 1.0

RETURN

70 I = A + 0.1

IF ( I .EQ. 0) GO TO 60

IF (1-20) 140*140.130

130 F=20.0

114
FCTRL=FCTI ( 20

)

DO 131 J=21»I

F=F+1.0

131 FCTRL=FCTRL*F

GO TO 150

140 FCTRL = FCTI { I ) 150 RETURN

END

C

SUBROUTINE CGAUSSEL ( C » NRD NRR ,NCC » NSF

)

C

DIMENSION C(NRD,NCC) »L( 128»2)

TYPE COMPLEX C , DET » P D Q

R

DATA (BITS = 1755 4000 OOOO 0000 B)

CALL ROLLCALL (48HCGAUSSEL 6/5/68 l-BANK BITS=2**-18 ) NR=NRR $ NC=NCC

IF(NC.LT,NR.OR.NR.GT,128.0R,NR.LE.O) CALL Q8QERR0R ( 0 9HB AD CALL.)

C INITIALIZE,

NSF = 0

NRM=NR-1 S NRP=NR+1 S D=(1.»0.) $ LSD=1

DO 1 KR=1»NR

L(KR.l)=KR

1 L(KR.2)=0

CALL Q9EXUN(EXUN)

IF(NR.EQ.l) GO TO 42

C ELIMINATION PHASE,

DO 41 KP=1»NRM

KPP=KP+1 S PM=0. $ MPN=0

C SEARCH COLUMN KP FROM DIAGONAL DOWN, FOR MAX PIVOT.

DO 2 KR=KP»NR

LKR=L (KR» 1

)

PT=CABS(C(LKR»KP)

)

IF(PT.LE.PM> GO TO 2 PM=PT S MPN=KR $ LMP=LKR

2 CONTINUE

C IF MAX PIVOT IS ZERO, MATRIX IS SINGULAR.

IFCMPN.EQ.O) GO TO 9

NSF=NSF+1

IF(MPN.EQ.KP) GO TO 3

C NEW ROW NUMBER KP HAS MAX PIVOT.

LSD=-LSD

L(KP»2 ) =L(MPN,1 )=L{KP,1

)

L(KP,l)=LMP

C ROW OPERATIONS TO ZERO COLUMN KP BELOW DIAGONAL.

3 MKP=L(KP,1)

P=C(MKP,KP) S D=D*P

DO 41 KR=KPP,NR

MKR=L(KR)

Q=C {MKR,KP ) /P

IF(REAL(Q) .EQ.O..AND.AIMAG(Q) .EQ.O. ) GO TO 41

C SUBTRACT Q * PIVOT ROW FROM ROW KR.

DO 4 LC=KPP»NC

R=Q*C{MKP,LC)

C(MKR,LC> =C(MKR»LC)-R

115
4 IF(CABS(C(MKR»LC) ) .LT,CABS(R)*BITS) C ( MKR » LC » = ( 0 . , 0 .

41 CONTINUE

C LOWER RIGHT HAND CORNER.

42 LNR=L(NR,1) S P=C{LNR»NR)

IF(REAL(P) .EQ.O..AND.AIMAG(P) .EQ.O, ) GO TO 9

NSF=NSF+1

D=D*P*LSD

IF{NR,EQ.NC) GO TO 8

C BACK SOLUTION PHASE.

DO 61 MC=NRP»NC

C(LNR.MC)=C(LNR,MC) /P

IF(NR.EQ.l) GO TO 61

DO 6 LL=1»NRM

KR=NR-LL S MR=L{KR,1) $ KRP=KR+1

DO 5 MS=KRP»NR

LMS=L{MS»1)

R=C (MR»MS )*CLMSjMC)

C(MR.MC)=C(MR,MC)-R

5 IF(CABS(C(MR»MC) ) .LT.CABS(K)*BITS) C ( MR » MC ) = ( 0 . , 0 .

6 C(MR»MC)=C(MR,MC) /C(MR»KR)

61 CONTINUE

C SHUFFLE SOLUTION ROWS BACK TO NATURAL ORDER.

DO 71 LL=1»NRM

KR=NR-LL

MKR=L(KR»2)

IF(MKR.EQ.O' GO TO 71

MKP=L(KR,1)

DO 7 LC=NRP»NC

Q=C(MKR»LC)

C(MKR,LC)=C(MKP,LC)

7 C{MKP,LC)=Q

71 CONTINUE

C NORMAL AND SINGULAR RETURNS. GOOD SOLUTION COULD HAVE D=0.

8 C( 1 ,1 )=D $ GO TO 91

9 C(1.1)=(0.»0.)

91 CALL S9FAULT{ EXUN) $ RETURN

END

C

SUBROUTINE BSJYOl (X» AJO» YO AJ1» YD

C

C CALCULATION OF THE BESSEL FUNCTIONS JO YO, Jl, AND Yl FOR AN

C ARGUMENT X

DIMENSION A(7)» B{7)» C(7), D(7), E(7), F(7). G(7)» H(7)

DATA ((All), I = 1,7) = 0.00021, -0.0039444, 0.0444479,

1 -0.3163866, 1.2656208, -2.2499997, 1.0)

DATA ((B(I)» I = 1,7) = -0.00024846, 0.00427916, -0.04261214,

1 0.25300117, -0.74350384, 0.60559366, 0.36746691)

DATA (IC(I)» I = 1,7) = 0.00014476, -0.00072805, 0'.00137237,

1 -0.00009512, -0.0055274, -0.00000077, 0.79788456)

DATA ((D(I)» I = 1,7) = 0.00013558, -0.00029333, -0.00054125,

1 0.00262573, -0.00003954, -0.04166397, -0,78539816)

DATA ((E(I)» I = 1,7) = 0.00001109, -0.00031761, 0.00443319,

1 -0.03954289, 0.21093573, -0.56249985, 0.5)

DATA ((F(I)» I = 1,7) = 0.0027873, -0.0400976, 0.3123951,

116
1 -1,3164827» 2.16827U9* 0.2212091, -0.6366198)

DATA ((G(I)» I = 1»7) = -0, 00020033* 0,00113653» -0.00249511»

1 0.0O017105» 0.01659667* 0,00000156» 0.79788456)

DATA ((H(I)» I = 1»7) = -0.00029166, 0.00079824, 0.00074348,

1 -0.00637879, 0.0000565, 0.12499612, -2.35619449)

AX = ABSF(X)

IF (AX.GT.Q.O) GO TO 10

AJO =1.0

YO = -l.E+030

AJl = 0.0

Yl = -l,E+030

RETURN

10 IF IAX.GT.3.0) GO TO 50

XX = (AX/3.0) ** 2 AJO = All)

YO = Bd)

AJl = Ed)

Yl = F( 1 ) DO 20 M = 2»7

AJO = AJO * XX + A(M)

YO = YO * XX + B(M)

AJl = AJl * XX + E(M)

20 Yl = Yl * XX + F(M)

AJl = AJl * X

ALF = 0.6366197724 * LOGF(0.5 * AX

)

YO = YO + ALF * AJO

Yl = Yl/X + ALF * AJl

RETURN

50 X3 = 3.0/AX

FO = C( 1)

THO = D( 1 ) Fl = G( 1)

THl = H(l

)

DO 60 M = 2»7

FO = FO * X3 + C(M)

THO = THO * X3 + D(M)

Fl = Fl X3 + G(M)

60 THl = THl * X3 + H(M)

THO = THO + AX

THl = THl + AX

XS = l./SQRT{AX)

AJO= XS * FO * COSF(THO)

YO = XS * FO * SINF(THO)

AJ1= XS Fl * COSF(THl)

Yl = XS * Fl * SINF(THI)

RETURN

END

C

FUNCTION COSINT(X)

C

C CALCULATION OF THE COSINE INTEGRAL

TYPE DOUBLE Y2 , PROD , SUM , PT , DK

IF(X.LE.O.) GO TO 50

X2 = X * X

117
IF(X.GT.20.^ GO TO 30

Y2 = DBLE(X2)

PROD = -Y2 * 0.5

SUM = PROD * 0.5

DO 10 K = 2»50

DK = 2 * K

PROD = -PROD * Y2/(DK*(DK - 1.))

SUM = SUM + PROD/DK

PT = ABSIPROD * 1.D+10>

IF(ABS(SUM) .GT.PT) GO TO 20

10 CONTINUE

20 SS = SNGL(SUM)

COSINT = SS + 0.5772156649 + LOGF(X>

RETURN

30 FA = 1.

FB = 1.

PO = 1.

X2 = 1./X2

DO 40 K = 1»10

AK = 2 * K

PO = -PO * AK X2

FA = FA + PO

PO = PO ^^ (AK + 1.)

FB = FB + PO

PA = ABS(PO * I.E+IO)

IF(PA.LE.FB) GO TO 45

40 CONTINUE

45 FX = FA/X

GX = FB * X2

COSINT = FX SIN(X) - GX * COS(X)

RETURN

50 WRITE (61»100) X 100 FORMAT (* X LESS OR EQUAL TO ZERO, X = *E17.9)

RETURN

END

FUNCTION FOUTR (DOM)

FOURIERTRANSFORM OF THERMAL AVERAGE FOR MODIFIED IMPACTTHEORY

COMPLEX FOUTR

COMMON/FDAT/Pl ,P2,B1 ,A2 ,B2 PPFF

ARG = ABSF(DOM)

CC = P2 ^ ^„

IF (ARG. GT. 1.22474) CC = P2 - 2

. *LOGF ( ARG/ 1 . 22474 ) *PPFF

FOUTR = 0.3183099 * CMPLX ( 0 .

» -CC

)

• RETURN

END

118
REFERENCES

j F. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions

ed. by M. Abramowitz and I. Stegun (NBS Applied Mathematics Series

55, 1969).

P. W. Anderson and J. D. Talman, Proceedings of the Conf. on the

Broadening of Spectral Lines, p. 29, Pittsburgh 1955, also Bell

Telephone System Monograph 31 17.

M. E. Bacon and D, F. Edwards, Phys. Rev. 170 , 125 (1968).

j M. Bacon, K. Y. Shen and J. Cooper, Phys. Rev. , to be published,

i M. Baranger, Phys. Rev. Ill , 494 (1958); 11?, 855 (1958).

I'

i M. Baranger and B. Mozer, Phys. Rev. 115 , 521 (1959); 118 , 626

I (I960).

M. Baranger, Atomic and Molecular Processes p. 493, ed. D. R. Bates

(New York, Academic Press 1962).

H. Bateman, Higher Transcendental Functions Vol. II (McGraw Hill,

' New York, 1953 ).

R. D. Bengtson, P. Kepple and J. D. Tannich, to be published.

A. Bethe and E. Salpeter, Quantum Mechanics of One-and TwoElectron Atoms (Springer Verlag, Berlin 1957).

J. W. Birkeland, J. P. Oss and W. G. Braun, Phys. Rev. 178 , No. 1,

368 (1969).

I G. Boldt and W. S. Cooper, Z. Naturf or schung 19a , 968(1964).

W. R. Chappell, J. Cooper and E. W. Smith, J. Q. S. R. T. 9, 149

(1969).

J. Cooper, Phys. Rev. LettersJ/7, 991 (1966).

J. Cooper, Rev. Mod. Phys. 39_, 167 (1967).

A. R. Edmonds, Angular Momentum in Quantum Mechanics (Princeton

University Press I960).

119
F. N. Edmonds, H. Schluter and D. C. Wells, Mem. Roy. Astr. Soc.

71^, 271 (1967).

R. C. Elton and H. R. Griem, Phys. Rev. US, A1550 (1964).

E. Ferguson and H. SchMter, Ann. Phys. 22^, 351 (1963).

J. B. Gerardo and R. A. Hill, Phys. Rev. Lett. _17, 623 (1966).

H. R. Griem, A. C. Kolb and K. Y. Shen, Phys. Rev. 116 , 4 (1959).

H. R. Griem, Astrophys. J. 136, 422 (1962).

H. R. Griem, A. C. Kolb and K. Y. Shen, Astrophys. J. 135 , 272

(1962).

H. R. Griem, Phys. Rev. 140 , A1140 (1965).

H. R. Griem, Astrophys. J. 147 , 1092, (1967).

H. R. Griem, Astrophys. J. 148 , 547 (1967).

C. F. Hooper, Phys. Rev. J_65, 215 (1968).

C. F. Hooper, Phys. Rev. 169 , 193 (1968).

J. W. B. Hughes, Proc. Phys. Soc. 9J^, 810 (1967).

P. Kepple and H. R. Griem, Phys. Rev. 173 , 317 (1968).

M. Lewis, Phys. Rev. 121 , 501 (1961).

Eo Merzbacher, Quantum Mechanics, Chap. 16 (J. Wiley, New York

1961).

H. Pfennig, E. Trefftz and C. R. Vidal, J. Q. S. R. T. _6^ 557 (1966).

H. Pfennig and E. Trefftz, A. Naturforschung 21a , 697 (1966).

H. Pfennig, Z. Naturforschung 21a , 1648 (1966).

H. Schluter and C. Avila, Astroph. J. 144, 785, 1966.

120
K. Y. Shen and J, Cooper, Astroph. J. 155 , 37 (1969)o

E. W. Smith and C, F. Hooper, Phys. Rev. 157 , 126 (1967).

E. W. Smith, C. R. Vidal and J. Cooper, J. Res. Natl. Bur. Std.

73A, 389 (1969).

E. W. Smith, C. R. Vidal and J. Cooper, J. Res. Natl. Bur. Std.

73A, 405 (1969)

E. W. Smith, J. Cooper and C. R. Vidal, Phys. Rev. to be published.

A. Unsold, Physik der Sternatmosph'aren (Springer Verlag, Berlin

1955).

H. Van Regemorter, Comptes Rendus 259 , 3979 (1964).

C. R. Vidal, Z. Naturfors Chung 19a, 947 (1964).

C. R. Vidal, Proceedings of the 7th. Intern. Conf. on Phenomena in

Ionized Gases, p. 168 Belgrade 1965.

B. Wende, Z. f. Angew, Physik 22, 181, 1967.

121
Figure 1. Schematic picture of the collision sphere showing the

Debye sphere, a strong collision sphere and a straight line classical

path trajectory.

122
Figure 2. The thermal average F of the time development operator

normalized with respect to the static, small interaction time asymptote F as a function of the normalized time s =*(ju • t. The two curves

are obtained with two different lower ^utoff paran?eters in the pintegral, n . = 0 and n . = ^ + n a . mm mm o

123
F(S)

Fo (S)

Log,o S

Figure 3. The thermal average F of the time development operator

normalized with respect to the static, small interaction time asymptote F as a function of the normalized time s = % t. The two sets of o p curves are obtained with two different lower cutoff parameters in the

2

p -integral, p . = 0 and p . - + n a . The three different curves mm mm o

in every set correspond to different Stark components characterized by

the quantum number n^ = nq-n'q'.

124
Log,o Aoip —

Figure 4. The Fourier transform of the thermal average i, normalized

with respect to the static, large frequency limit i as a function of the

normalized frequency Aujp = (Atu-Auu.« B)/uj' .

125
n,^ = 90 16 3

Figure 5. The Fourier transform of the thermal average i, normalized

with respect to the static, large frequency limit i as a function of the

normalized frequency Aui)-^ = (Aw-AU).* P • The three different curves correspond to different Stark corn.pone:&s characterized by the quantum number n, = nq~n'q'.

126
3

(U

o

)4 o

o -

c o

. f!

"+-1

CUD TM

>

1—1niau

Oh

cr

-l-> i-H —

1

3

<a

i)

•1-4 s >^

CO O

=1 ri 0 <uqu 0 <u •w u

Eh 4JU

<u <u ds ISl •H

re

re T—

1

0 B

u • iH •H O h

q: 3

cr 3

127
3

to

O 0) C > ° !^ ^

!->

O •->

<-* ^5 O ^ >

0)

CO CO "J ^

•H ^ O

(-> ^ fl} •H +J ^ £ (u 5 O to - cd . OS

+^ .rt \ O

cn g<

y 5 o C r u

(U TO ^ II o ^-

CO

M 1-4 "

1—1 3 O

(U

CO

U

01

u o o

CO .H

iXj

T3 Tlo

a,

CO

IhOU

cr

cr

II usds

(ti cr

tr 3or3

128
1 EoO

O

ro

If) 00

II II

<u c 1- 3

0)a

-l-> )^oS

;4o >

.—I t ^

U- -l-i

!->Oa

;-i o 1 to 3

<^

J ^ a " ^ o C.

cr W

14 ? ^ ,-1 a

nJu

14(U •HUO

4:1 (0 .s

1-1

a

CO

2i -

+j o

OJ o

>-i 3o C n] Li ^ ^ u ^ .t5 o

0) (0

u

CO

(T. 3 3

129
130
ro

1E

II oroO O

* lO O ro• 00

II N 11 c c

CD

CD

Q.

i3

^ 3

<

ro

CVJ

00

i-> nt

M

£J isi

o . •H M-( t)JO fi

<o t o

rO

. -H O O ^ h

<u .2 <i>

h :S

(OI='Jfi^'«V)I

131
- o< o s O 0) ^ uo

isl ^

•H cm

dUO

ft mO

-i-> o

(U X!o

CO

ss

(U

(U I—i —

(U 5- 2 X u c o o

PQ S ^ o

-i->C

flOO

(U

Oi CO

d ^ 3

0)

(U

CO

Tl

<U • H

r—

to (fl

h > t«

or

132
lO o m^

rj TO

U ^ <U

0) ^, O

tg

O 3

.-I 43

CO

O ^ .2 ^ o I 0) In 3

1-1 ocH

«5 pq 0) +J j3. U CO ^ O o

flJ Q,e (U 3 C 43 M ><

CO fa rt

O

or

133
A<«J/a3p

Figure 13. The final line profiles for the density and temperature

parameters of Vidal, 1965.

134
4 5 6 7 8 9 10 11 12 13 14

Figure 14. Plot of the electron density values as a function of the

principal quantum number n, which have been evaluated by Vidal,

1965, under the assumption that the AX~^/^ -wing revealed by the

experiment is identical with the asymptotic Holtsmark /^,X"^/^- wing

for electrons and ions.

135
10l! I I I I I

k 11.

12

ne = 1.3-l0'^ - T=1850 K : KEPPLE a GRIEM

BENGTSON ft KEPPLE ^ — ASYMPTOTIC Aar^/2-wiNG

10 0.1

Figure 15. The Balmer line profiles for the density and temperature

parameters of Vidal, 1965, as calculated by the indicated authors.

136
Figure 16. The Lyman-p profile for two different strong collisions

cutoffs p . . mm

137
T3

I—

• H

u o

(XJ

(1) (H

0)or

o

o

•)-> CO

T5

o;

CO

^< c!

o A

o

l-H O

(U r-H ^^ (U d I 00 OJ h o

o

0)

CO °, ^ \ U I cd CD G CD.

138
Q.

?33<

CO .H o ^ o

XI

(U •H^3

4-> +J os ^ ^

>^

TM (4 a o E5 5o

00 i

l-H

O

0)

u <u 3 I h o

flo

a.Hoo

4-> CO

CO

COOsu

+J

13o

0) 43Oc > O («

•rj ca

•H II CO OD o

(U

139
lO

OJ o

OJ o lO« o

13 o

OJc

0) ^ X!

,E) to .1-1 X ^ 5 ° S ?i S O Pi o o

(U

•4-> Xi

OJoaoo

AS . u

00 03

II

w

a; o ^tj mm ° I - O «3 O

.2 ^ S u u u o >^ ^ o ^ O •4-' O ^ el £-

<U CD

q:

140
THE NATIONAL ECONOMIC GOAL

Sustained maximum growth in a free market economy, without inflation, under conditions of full employment

and equal opportunity

THE DEPARTMENT OF COMMERCE

The historic mission of the Department

is "to foster, promote and develop the

foreign and domestic commerce" of the United States. This has evolved, as a

result of legislative and administrative

additions, to encompass broadly the re- sponsibility to foster, serve and promote

the nation's economic development and

technological advancement. The Depart- ment seeks to fulfill this mission through

these activities:

MISSION AND

FUNCTIONS

OF THE

DEPARTMENT OF

COMMERCE

"to foster, serve and

promote the nation's

economic development

and technological

advancement"

Participating with other government

agencies in the

creation of national

policy, through the

President's Cabinet and its subdivisions.

• Cabinet Committee

on Economic Policy

• Urban Affairs Council

•Environmental

Quality Council

Promoting progressive business policies and

growth.

• Business and

Defense Services

Administration

• Office of Field Services

Assisting states, communities and

individuals toward economic progress.

• Economic

Development

Administration

• Regional Planning Commissions

• Office of Minority Business Enterprise

NOTE: This schematic Is neither an organization chart nor a progrann outline for budget purposes. It is a general statement

of the Department's mission in relation to the national goal

of economic development.

JULY 1969

Strengthening

the international economic position

of the United

States.

• Bureau of International Commerce

• Office of Foreign Commercial

Services

• Office of Foreign

Direct Investments

• United States

Travel Service

• Maritime

Administration

Assuring effective use and growth of the

nation's scientific and technical resources.

• Environmental Science Services

Administration

• Patent Office

• National Bureau of Standards

• Office of Telecommunications

• Office of State

Technical Services

Acquiring, analyzing and disseminating

information concerning the nation and

the economy to help achieve increased

social and economic

benefit.

• Bureau of the Census

• Office of Business Economics

6 U. S. GOVERNMENT PRINTING OFFICE : 1970 O - 381-921

ll



NBS TECHNICAL PUBLICATIONS

PERIODICALS

JOURNAL OF RESEARCH reports National

Bureau of Standards research and development in physics, mathematics, chemistry, and engineering.

Comprehensive scientific papers give complete details of the work, including laboratory data, experimental

procedures, and theoretical and mathematical analyses. Illustrated with photographs, drawings, and

charts.

Published in three sections, available separately: • Physics and Chemistry

Papers of interest primarily to scientists working in these fields. This section covers a broad range of

physical and chemical research, with major emphasis on standards of physical measurement, fundamental

constants, and properties of matter. Issued six times

a year. Annual subscription: Domestic, $9.50; foreigii, $11.75*. • Mathematical Sciences

Studies and compilations designed mainly for the mathematician and theoretical physicist. Topics in mathematical statistics, theory of experiment design,

numerical analysis, theoretical physics and chemistry, logical design and programming of computers

and computer systems. Short numerical tables. Issued quarterly. Annual subscription: Domestic,

$5.00; foreign, $6.25*. • Engineering and Instrumentation

Reporting results of interest chiefly to the engineer

and the applied scientist. This section includes many

of the new developments in instrumentation resulting

from the Bureau's work in physical measurement,

data processing, and development of test methods.

It will also cover some of the work in acoustics,

applied mechanics, building research, and cryogenic

engineering. Issued quarterly. Annual subscription:

Domestic, $5.00; foreign, $6.25*.

TECHNICAL NEWS BULLETIN

The best single source of information concerning the Bureau's research, developmental, cooperative and

publication activities, this monthly publication is designed for the industry-oriented individual whose

daily work involves intimate contact with science and

technology

for engineers, chemists, physicists, re- search managers, product-development managers, and

company executives. Annual subscription: Domestic,

$3.00; foreign, $4.00*.

• Difference in price is due to extra cost of foreign mailinff. Order NBS publications from:

NONPERIODICALS

Applied Mathematics Series. Mathematical tables, manuals, and studies. Building Science Series. Research results, test methods, and performance criteria of building materials, components, systems, and structures. Handbooks. Recommended codes of engineering

and industrial practice (including safety codes) de- veloped in cooperation with interested industries,

professional organizations, and regulatory bodies. Special Publications. Proceedings of NBS conferences, bibliographies, annual reports, wall charts,

pamphlets, etc. Monographs. Major contributions to the technical

literature on various subjects related to the Bureau's

scientific and technical activities. National Standard Reference Data Series.

NSRDS provides quantitive data on the physical

and chemical properties of materials, compiled from

the world's literature and critically evaluated. Product Standards. Provide requirements for sizes, types, quality and methods for testing various industrial products. These standards are developed cooperatively with interested Government and industry groups

and provide the basis for common understanding of product characteristics for both buyers and sellers. Their use is voluntary.

Technical Notes. This series consists of communications and reports (covering both other agency and

NBS-sponsored work) of limited or transitory interest. Federal Information Processing Standards Publications. This series is the official publication within

the Federal Government for information on standards

adopted and promulgated under the Public Law

89-306, and Bureau of the Budget Circular A-86

entitled. Standardization of Data Elements and Codes

in Data Systems.

CLEARINGHOUSE

The Clearinghouse for Federal Scientific and

Technical Information, operated by NBS, supplies

unclassified information related to Government-generated science and technology in defense, space, atomic energy, and other national programs. For

further information on Clearinghouse services, write: Clearinghouse

U.S. Department of Commerce

Springfield, Virginia 22151

Superintendent of Documents

Government Printing Office

Washington, D.C. 20402
U.S. DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20230

OFFICIAL BUSINESS POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERCE