- Email: [email protected]
Spectral line inversion as a diagnostic tool
J. Quonf .Spec~roscRadiat Printed in the U.S.A.
Trons/k Vol. 33, No. 2, pp. 93-100,
OOZZ-4073/85 $3.00 + .Xl Pergamon Press Ltd.
I985
SPECTRAL LINE INVERSION DIAGNOSTIC TOOL
AS A
A. KAVETSKY~and B. J. O’MARA Department of Physics, University of Queensland, St. Lucia, Queensland 4067, Australia (Received 6 February 1984)
Abstract-A general technique for inverting spectral line profiles to obtain physical information about a semi-infinite plasma is described. The technique allows for distinct line and continuum source functions, the presence of both stochastic and systematic velocity fields, and permits all spectral lines in a multiplet to be inverted simultaneously. 1. INTRODUCTION
Halt’ and Hearn and Halt’ developed a technique for inverting a single spectral line profile to obtain the source function directly from the emergent intensity in the line. In the next section, the technique of Hearn and Holt (hereafter referred to as H&H) is summarized. Extensions of this technique are then described which make it more realistic and capable of determining additional physical properties of the plasma producing the radiation. As the inversion technique to be described was specifically developed as a tool for analysing the solar atmosphere, the term atmosphere will be used throughout; however, the technique can be applied in any situation where radiation emanates from a plane-parallel semi-infinite plasma. 2. THE INVERSION
TECHNIQUE
OF HEARN
AND HOLT
The emergent intensity at p = cos 0 from a plane-parallel
atmosphere
is given by
(1) where S, is the total source function and z1 is the total optical depth. The total optical depth is related to t, the line optical depth at line center, by dz, = [6 (1) + r]dt,
(2)
where r is the ratio of continuum opacity to line-center line opacity and the function 4(n) is (assuming the absorption can be characterized by a Voigt profile) given by
H is the Voigt function while v is given by
v=@-&)/A,
(4)
A being the Doppler width in wavelength units and & the wavelength of line center.
In the inversion technique of H&H, the atmosphere forming the line is divided into a number of segments in line-center line optical depth t. The total source function SL(t) is assumed to be frequency independent, linear within each segment and continuous across segment boundaries, while the functions r(t), a(t), and A(t) of Eqs. (2)-(4) are represented by constant values within each segment. This model of line formation is illustrated in Fig. 1. tPresent address: Julius Kruttschnitt
Mineral Research Centre, Isles Road, Indooroopilly,
Queensland 4068, Australia. 93
Brisbane,
A. KAVETSKY and B. J. O’MARA
94
AI
SEGMENT
Tl
A2
i
1
1
SEGMENT
T2
i
2
“3
i I
SEGMENT
T3
3 t
Fig. 1.
This method of specifying the properties of each segment allows the equation of transfer to be solved analytically by integrating with respect to t through each segment, and an analytic expression for the emergent intensity from a general N-segment atmosphere to be obtained. It is, for position p = cos 8 on the disk,
Here, T, = 0 (boundary of the most shallow segment) and TN+, + co; 4j is the value of the Voigt profile in segment j, the other symbols are defined in Fig. 1. The fit between the observed and calculated profiles is optimized by minimizing, as a function of the parameters defined in Fig. 1, the function P((parameters))
i=l
p = C
Zi(obs) - Z,((parameters)) Zi(obs)
12 2
where Zi(obs) and Zi((p arameters)) are the observed and calculated intensities respectively at the ith wavelength point. The function F((parameters)) is minimized numerically using the conjugate gradient method of Powell.3 As there is a limited amount of information contained in a line profile, a set of parameters determined by minimizing F in Eq. (6) may not be unique, i.e., other sets of parameters may exist which also minimize F. An inversion yielding a given set of parameters is regarded as being unique when it yields the same solution from several widely-differing sets of starting values of the parameters. Starting with simple one-segment models, extra segments and parameters are added until the solutions obtained become dependent on the starting values of the parameters. The last solution which is independent of the starting values is considered to contain the maximum number of physically meaningful parameters which may be specified by the line profile.
95
Spectral line inversion as a diagnostic tool 3. EXTENSIONS
OF THE INVERSION
TECHNIQUE
OF HEARN
AND HOLT
(i) Wavelength dependent source function The total source function S, is
where SC is the continuum source function, S, the line source function, $x the line absorption profile at wavelength X and r is the ratio of continuum opacity to line opacity measured at the center of the line. For example, in the solar atmosphere, the continuum source function is, to a good approximation, the Planck function, which is constant across the small wavelength interval occupied by a line profile. In general, the total source function is wavelength dependent. However, under certain circumstances this is not the case, namely: (a) if the line source function is m LTE throughout the whole region of line formation and is therefore equal to the Plank function; or (b) if the line source function departs from LTE but 41 is much greater than r; in this case the total source function is essentially equal to the line source function, which is independent of wavelength if one assumes complete redistribution. This is often the situation in the cores of absorption line profiles. However, the total source function will depend strongly on wavelength if departures from LTE extend into the atmosphere to a region forming the part of the line for which c#J~, is comparable to r. A detailed discussion of redistribution functions for resonance lines and the conditions under which complete redistribution can be assumed has been given by Mihalas.4 With complete redistribution, S, is frequency-independent but due to the fact that r # 0 in Eq. (7) we still have a frequency-dependent total source function S,. We now include the possibility of such a wavelength dependent S, by explicitly representing SCand S, by linear segments with the segment boundaries assumed to be the same for both SC and S, (this is not a limitation since the slope of either source function can be kept unchanged from one segment to the next). Let S,(i) = bi + ci(t - Ti)
(8)
S,(i) =A + hi(t - TJ
(9)
and
be the continuum and line source functions in the ith segment. Then the equation for the emergent intensity becomes
where
(11) and
C,(n)= [rjcj + dJji(J+>hjl/[4j(~) + rj].
(12)
96
A.
KAVETSKY and B. J. O’MARA
F((parameters)) is then minimized in the usual way, using the conjugate gradient method of Powell3 to determine the unique set of parameters which describe a line with a wavelength dependent source function. (ii) The inversion qf’multip1et.v In this section, the method of extending the inversion technique of Hearn and Holt to invert simultaneously all nonoverlapping lines in a multiplet is described. For convenience, the following discussion is confined to a doublet: extension to triplets and higher multiplets is trivial. Consider an upper-level doublet, formed by transitions from closely spaced upper levels (1) and (2) to a lower level. At a given physical depth, it is assumed that: (a) The continuum source function is the same for the two lines, since they are closely spaced in wavelength. In addition, it is assumed that the line source functions are also equal at the same physical depth. This is not strictly true; the extent of approach to equality is determined by the amount of coupling between levels (1) and (2) (see Ref. 5, p. 92). Athay concludes that equality of line source functions to about 10% is quite common. (b) Assuming equal collisional damping for the lines, the line opacities are in the same ratio as thef-values for the two lines, i.e.,
as stimulated emission is negligible. (c) The Doppler widths are the same for the two lines. (d) The damping constants are the same for the two lines. With these assumptions, the parameters for line 2, which must be written in terms of the line-center line optical depth of this line, are related to the parameters of line 1 by
b,(2)= W),
7x2) = Hw-~~N~u)~
a4
cd21
= 4(l),
=
42) = MlW~l)lf(2)1~
f,(2) =fi(l>, ci(1WC1)/f(2)12 4) = a>, ri(2)= ri(lW(lM2)l~
(13)
Thus, the basic parameters are those of line 1, the parameters for line 2 (and any other lines of the multiplet) being derived from these by the Eq. (13). The emergent intensity in both lines is computed by Eq. (lo), thus enabling the function F((parameters)) of Eq. (6) to be calculated. The sum in this equation now extends over both lines. F is again minimized with respect to the parameters. Simultaneous inversion of all the lines in a multiplet is superior to inversion of a single line as it improves the stability of the inversion process and under the right circumstances can yield a larger set of uniquely determined parameters and therefore lead to a better description of the physical properties of the atmosphere. (iii) Determination of systematic velocity jields The presence of stochastic velocity fields has already been implicitly included in the Doppler width. A systematic motion of the atmosphere that varies with depth leads to asymmetries in the emergent line profiles and these asymmetries can be exploited to determine the systematic velocity field. Most of the methods that have been developed for doing this have the following features in common: Doppler shifts are determined at various positions in the complete profile as in the bisector shift method described by Kulander and Jefferies6 or from an analysis of filtergrams as described by Parnell and Beckers;’ the velocities associated with these Doppler shifts are then assigned to particular depths in the atmosphere by means of the Eddington-Barbier relation as described by Kulander and Jefferies or by means of the response function technique of Beckers and MilkeyE and Caccin and Marmolino.’ It is this final step in the analysis where most of the difficulties arise with these methods. Indeed, Parnell and Beckers’ have shown that the assigned optical depths depend on the form of the velocity field itself. This feedback between the depth of formation and the velocity field inherent in these methods can be traced to an
Spectral line inversion as a diagnostic tool
97
attempt to derive the velocity field in a manner which is independent of the other factors that are important in the formation of the line. In contrast, our inversion technique yields all of the parameters that are important in the formation of the line simultaneously and, as a result, no independent determination of the depth of formation, with its associated problems, is required. The immediate effect of introducing a systematic depth-dependent velocity field into a plane-parallel atmosphere is to displace the center of the local absorption profile. This, in turn, results in a modification of both the depth dependence and frequency dependence of the total source function. In our method, it will be assumed that the total source function is frequency independent. This assumption is not strictly necessary and is made purely for convenience. Any modification in the depth dependence of the source function produced by the presence of the systematic velocity field is automatically taken care of in the inversion procedure (to the extent that a “unique” set of parameters can be determined). With this assumption, the emergent intensity is still given by Eq. (1) but the relationship between the total optical depth and the line-center line optical depth given by Eq. (2) is modified. If u(t) is the depth dependent velocity in the atmosphere, then dr, = [9(n + LX)+ r]dt,
(14)
W) = &u(t)lc
(15)
where
is the depth dependent Doppler shift in the center of the absorption profile, with u taken positive towards the observer. The concept of line-center line optical depth loses its usual significance in the presence of a systematic velocity field. The variable t is the line-center line optical depth in a stationary atmosphere and does not correspond to the line optical depth at any wavelength in the observed profile; the wavelength of the observed line center depends on the specific form of the systematic velocity field. The velocity field u(t) is represented by a number of linear segments with segment boundaries which do not necessarily coincide with the source function segment boundaries. With this representation of v(t), Eq. (14) can be integrated to yield zl(t) = G(1, t, (parameters)),
(16)
where (parameters) represent the parameters involved in the piecewise continuous representation of u(t). This relation between r1 and t is then used in Eq. (1) to obtain a relationship between the emergent intensity in the line and the complete set of parameters describing the formation of the line. The inversion procedure is then identical to that already described, except that the parameters determined now include a representation of o(t). In the inversion process it is necessary to evaluate the G-function a large number of times. Efficient methods for its computation are described in the Appendix. ‘RI test the velocity inversion procedure, synthetic line profiles were constructed using a piecewise linear source function and various forms for the synthetic velocity field. The inversion procedure was then applied to these synthetic profiles. In the first series of tests, synthetic profiles were constructed in which the synthetic velocity field was represented by two linear segments with the segment boundary differing from the source function segment boundary. Different choices of the starting values of the parameters in the inversion procedure led to derived parameters in good agreement with those used to construct the synthetic profiles; which confirms that the inversion technique is capable of recovering a synthetic velocity field which has the same mathematical form as that assumed in the inversion process. In a second series of tests,‘synthetic velocity fields in which the velocity increased in proportion to the square of the continuum optical depth z, up to a certain limiting value of z, and then remained constant thereafter, were used to construct synthetic profiles. Two such synthetic profiles were constructed, one with v = 3rC2,z, < 4; u = 48 km/s, z, > 4, and
A. KAVETSKY
98
and B. J. O’MARA
1 ........LINE
BISECTOR
-
INPUT VELOCITY
---
INVERSION
1
VELOCITY
VELOCITY
! ,’
-0
0.5
1 .o
CONTINUUM
1.5 OPTICAL
2.0
2.5
DEPTH
Fig. 2.
the other with v = 12rc2, r, < 2; v = 48 km/s, r, > 2. For this range of velocities, two velocity segements were uniquely specified by the inversion. Although the quality of the fit to the synthetic profiles was not as good as that obtained in the first series of tests since the sum of the squares of the relative deviations F was a little larger, the recovered values of the parameters not associated with the velocity field were in good agreement with the parameters used to construct the synthetic profiles. This result indicates that the recovery of these parameters is not seriously affected by representing a quadratically varying velocity by two linear segments. Figure 2 compares the synthetic input velocity with the velocity recovered by inversion. Clearly the synthetic velocity is quite well represented by two linear segments and is accurately recovered by the inversion. Finally the stability of the inversion procedure was tested by adding synthetic noise, at a level of 5% of the total intensity, to the synthetic profiles prior to inversion. Although, as expected, the quality of the inversion was degraded no evidence of numerical instability was found. To illustrate the advantages of the inversion procedure over the line-bisector method the velocity field was determined by the latter method from the synthetic profiles, using the procedure suggested by Kulander and Jefferies. 6 The results are displayed in Fig. 2. Because the Eddington-Barbier approximation is used, the line bisector method is restricted to depths r, 5 1, and when v = 127,*, where the velocity varies rapidly, it fails completely. These limitations were anticipated by Kulander and Jefferies. 4.
CONCLUSlONS
Although the inversion technique has been demonstrated by us only for synthetic data, examples of its application to the inversion of solar line profiles can be found in Hearn and Halt’ and Kavetsky and O’Mara,” where the method is used to deduce physical properties of the solar atmosphere. It is hoped that the method will be equally useful in the analysis of laboratory plasmas. Acknowledgements-We would like to thank J. N. Holt for several helpful discussions on the inversion method, J. Cooper for useful comments on the manuscript, and the staff at the Joint Institute for Laboratory Astrophysics for assistance in its preparation. One of us (A.K.) was the holder of a Commonwealth Postgraduate Research Award during the course of this work, for which he is grateful. This work was supported in part by funds from Australian Research Grants Commission.
Spectral
line inversion
as a diagnostic
tool
99
REFERENCES 1. J. N. Halt, Ph.D. Dissertation, University of Queensland, St. Lucia, Queensland, Australia, published). 2. A. G. Heam and J. N. Holt, Astron. Astrophys. 23, 347 (1973). 3. M. J. D. Powell, Cornput. J. 7, 155 (1964). 4. D. Mihalas, Stellar Atmospheres, 2nd Edn. Freeman, San Francisco (1978). 5. R. G. Athay, Radiation Transport in Spectral Lines. Reidel, Dordrecht (1972). 6. J. L. Kulander and J. T. Jefferies, Astrophys. J. 146, 194 (1966). 7. R. L. Pamell and J. M. Beckers, Solar Phys. 9, 35 (1969). 8. J. M. Beckers and R. W. Milkey, Solar Phys. 43, 189 (1975). 9. B. Caccin and C. Marmolino, Astron. Astrophys. 83, 73 (1980). 10. A. Kavetsky and B. J. O’Mara, Solar Phys. (submitted). 11. A. Curtis and M. R. Osborne, Comput. J. 9, 286 (1966). 12. C. T. Fike, Computer Evaluation of Mathematical Functions. Prentice-Hall, New York (1968). 13. A. Ralston, A First Course in Numerical Analysis. McGraw-Hill, New York (1965).
1972 (un-
APPENDIX Relationship of optical depths in the presence of a velocity field For the sake of simplicity, the following analysis is confined to an atmosphere with a single linear velocity segment and two source function segments. The results are easily extended to an atmosphere with any number of velocity and source function segments. For a single linear velocity segment, the Doppler shift in the local absorption profile defined in Eq. (15) can be written as a(t) = a,+ agt.
(Al)
the explicit form of the The G-function, defined by Eqs. (14) and (16), can then be determined by substituting absorption profile [Eq. (3) along with Eq. (Al)] into Eq. (14). In segment one where t < T2, collisional damping is assumed negligible so that the a-value in the Voigt profile can be set to zero. The expression for the G-function obtained by integrating over t for the segment is then
for all values of t < T,. When act/d, is small, numerical calculations using this equation are inaccurate due to loss of significant figures in the subtraction of the two nearly equal error functions. In this situation, it is convenient to approximate the difference by a derivative to get G(I., t; r, A,, a,, aJ=rt+texp[-r$+zy]. In segment
number
two,
(A3)
t > T2,
G(X, t; r, A,, AZ, a*, a,, ag) = r,T, + rz(t + where H is the Voigt function
A2
a,sH(a,,
(A4)
0)
of Eq. (3) and the X-function
is defined
by
x(c, a)=jymarctan(c-i)exp(-y2)dy.
(A5)
Some of the mathematical properties of ~(c, a) are given in Table 1. For similar reasons to those leading (A3), the difference of X-functions in Eq. (A4) may be approximated by a derivative to give
to Eq.
WA, t; r, A,, A,, a,, a,, q)=r,T,+r,(t
(‘46) when a&t - TJA,a, is small. The error function can be easily and rapidly calculated using the rational approximation of Curtis and Osborne.” To permit rapid computation of the X-function, ~(c, a)/c was represented (for fixed values of a) by a Chebyshev polynomial expansion in c for small and moderate values of [cl, i.e., for 0 < (cl < c0 for some moderate value of c,; c, decreases with increasing a and is N 70 for typical values of a. For values of ICI greater than co, where x(c, a) approaches its limit f7?/2, it is convenient to introduce the function fi(k, a) defined by
k
(A7)
KAVETSKY
A.
100 Table 1. Some properties
of the functions
and B. J. O’MARA
x(c, a) and Q(k, a) defined by Eqs. (A$
. _ Expression
Property
FUIlCtiOll
07
(A7), and (A8).
Value
3/2 x(c.d
alternative for c > 0
expression x ev(-y2)dy aTI erfc(a)
derivative
with
t0 C, i.e., x(c/a,a)
limit
S2Ck.a)
alternative for k > 0
as
respect
exp(a2)
a,,H(a,ac)
ax/at a +
0
(7 3/2/2)
erf Cc)
312
expression
+
erfc
($
+ ;
x
n(k,a)
limit
as
k +
n(k,a)
limit
as
k + _
0
,/”
arctan[
(i
- 3
-l I
ew-y2)dy
.1/2
0
where k = l/c.
(AS)
This function was approximated for the range 0 < Ik( < l/ c,,, i.e., c,, i ICI < co by a Chebyshev polynomial expansion as a function of k for fixed values of a. Representing x(c,a) in this way ensures an accurate approximation to this function for all values of c at a given value of a including very small and very large values. The techniques and equations used in the calculation of the Chebyshev coefficients are given in Fike.12 Chebyshev coefficients for each of these functions are available from the authors for a range of a values. With the above expressions for the G-function and the Chebyshev expansions allowing its rapid calculaton the relationship between the optical depth at any wavelength and the line-center line optical depth given by Eq. (16) is fully defined. In order to do the numerical integration for the emergent intensity (Eq. [l]) it is necessary to solve Eq. (16) for t given a value of ri. To do this, Eq. (16) is written in the form TV- G(I, and solved by using the secant
method
as described
t, parameters) by Ralston.”
= 0
(A9)
Trons/k Vol. 33, No. 2, pp. 93-100,
OOZZ-4073/85 $3.00 + .Xl Pergamon Press Ltd.
I985
SPECTRAL LINE INVERSION DIAGNOSTIC TOOL
AS A
A. KAVETSKY~and B. J. O’MARA Department of Physics, University of Queensland, St. Lucia, Queensland 4067, Australia (Received 6 February 1984)
Abstract-A general technique for inverting spectral line profiles to obtain physical information about a semi-infinite plasma is described. The technique allows for distinct line and continuum source functions, the presence of both stochastic and systematic velocity fields, and permits all spectral lines in a multiplet to be inverted simultaneously. 1. INTRODUCTION
Halt’ and Hearn and Halt’ developed a technique for inverting a single spectral line profile to obtain the source function directly from the emergent intensity in the line. In the next section, the technique of Hearn and Holt (hereafter referred to as H&H) is summarized. Extensions of this technique are then described which make it more realistic and capable of determining additional physical properties of the plasma producing the radiation. As the inversion technique to be described was specifically developed as a tool for analysing the solar atmosphere, the term atmosphere will be used throughout; however, the technique can be applied in any situation where radiation emanates from a plane-parallel semi-infinite plasma. 2. THE INVERSION
TECHNIQUE
OF HEARN
AND HOLT
The emergent intensity at p = cos 0 from a plane-parallel
atmosphere
is given by
(1) where S, is the total source function and z1 is the total optical depth. The total optical depth is related to t, the line optical depth at line center, by dz, = [6 (1) + r]dt,
(2)
where r is the ratio of continuum opacity to line-center line opacity and the function 4(n) is (assuming the absorption can be characterized by a Voigt profile) given by
H is the Voigt function while v is given by
v=@-&)/A,
(4)
A being the Doppler width in wavelength units and & the wavelength of line center.
In the inversion technique of H&H, the atmosphere forming the line is divided into a number of segments in line-center line optical depth t. The total source function SL(t) is assumed to be frequency independent, linear within each segment and continuous across segment boundaries, while the functions r(t), a(t), and A(t) of Eqs. (2)-(4) are represented by constant values within each segment. This model of line formation is illustrated in Fig. 1. tPresent address: Julius Kruttschnitt
Mineral Research Centre, Isles Road, Indooroopilly,
Queensland 4068, Australia. 93
Brisbane,
A. KAVETSKY and B. J. O’MARA
94
AI
SEGMENT
Tl
A2
i
1
1
SEGMENT
T2
i
2
“3
i I
SEGMENT
T3
3 t
Fig. 1.
This method of specifying the properties of each segment allows the equation of transfer to be solved analytically by integrating with respect to t through each segment, and an analytic expression for the emergent intensity from a general N-segment atmosphere to be obtained. It is, for position p = cos 8 on the disk,
Here, T, = 0 (boundary of the most shallow segment) and TN+, + co; 4j is the value of the Voigt profile in segment j, the other symbols are defined in Fig. 1. The fit between the observed and calculated profiles is optimized by minimizing, as a function of the parameters defined in Fig. 1, the function P((parameters))
i=l
p = C
Zi(obs) - Z,((parameters)) Zi(obs)
12 2
where Zi(obs) and Zi((p arameters)) are the observed and calculated intensities respectively at the ith wavelength point. The function F((parameters)) is minimized numerically using the conjugate gradient method of Powell.3 As there is a limited amount of information contained in a line profile, a set of parameters determined by minimizing F in Eq. (6) may not be unique, i.e., other sets of parameters may exist which also minimize F. An inversion yielding a given set of parameters is regarded as being unique when it yields the same solution from several widely-differing sets of starting values of the parameters. Starting with simple one-segment models, extra segments and parameters are added until the solutions obtained become dependent on the starting values of the parameters. The last solution which is independent of the starting values is considered to contain the maximum number of physically meaningful parameters which may be specified by the line profile.
95
Spectral line inversion as a diagnostic tool 3. EXTENSIONS
OF THE INVERSION
TECHNIQUE
OF HEARN
AND HOLT
(i) Wavelength dependent source function The total source function S, is
where SC is the continuum source function, S, the line source function, $x the line absorption profile at wavelength X and r is the ratio of continuum opacity to line opacity measured at the center of the line. For example, in the solar atmosphere, the continuum source function is, to a good approximation, the Planck function, which is constant across the small wavelength interval occupied by a line profile. In general, the total source function is wavelength dependent. However, under certain circumstances this is not the case, namely: (a) if the line source function is m LTE throughout the whole region of line formation and is therefore equal to the Plank function; or (b) if the line source function departs from LTE but 41 is much greater than r; in this case the total source function is essentially equal to the line source function, which is independent of wavelength if one assumes complete redistribution. This is often the situation in the cores of absorption line profiles. However, the total source function will depend strongly on wavelength if departures from LTE extend into the atmosphere to a region forming the part of the line for which c#J~, is comparable to r. A detailed discussion of redistribution functions for resonance lines and the conditions under which complete redistribution can be assumed has been given by Mihalas.4 With complete redistribution, S, is frequency-independent but due to the fact that r # 0 in Eq. (7) we still have a frequency-dependent total source function S,. We now include the possibility of such a wavelength dependent S, by explicitly representing SCand S, by linear segments with the segment boundaries assumed to be the same for both SC and S, (this is not a limitation since the slope of either source function can be kept unchanged from one segment to the next). Let S,(i) = bi + ci(t - Ti)
(8)
S,(i) =A + hi(t - TJ
(9)
and
be the continuum and line source functions in the ith segment. Then the equation for the emergent intensity becomes
where
(11) and
C,(n)= [rjcj + dJji(J+>hjl/[4j(~) + rj].
(12)
96
A.
KAVETSKY and B. J. O’MARA
F((parameters)) is then minimized in the usual way, using the conjugate gradient method of Powell3 to determine the unique set of parameters which describe a line with a wavelength dependent source function. (ii) The inversion qf’multip1et.v In this section, the method of extending the inversion technique of Hearn and Holt to invert simultaneously all nonoverlapping lines in a multiplet is described. For convenience, the following discussion is confined to a doublet: extension to triplets and higher multiplets is trivial. Consider an upper-level doublet, formed by transitions from closely spaced upper levels (1) and (2) to a lower level. At a given physical depth, it is assumed that: (a) The continuum source function is the same for the two lines, since they are closely spaced in wavelength. In addition, it is assumed that the line source functions are also equal at the same physical depth. This is not strictly true; the extent of approach to equality is determined by the amount of coupling between levels (1) and (2) (see Ref. 5, p. 92). Athay concludes that equality of line source functions to about 10% is quite common. (b) Assuming equal collisional damping for the lines, the line opacities are in the same ratio as thef-values for the two lines, i.e.,
as stimulated emission is negligible. (c) The Doppler widths are the same for the two lines. (d) The damping constants are the same for the two lines. With these assumptions, the parameters for line 2, which must be written in terms of the line-center line optical depth of this line, are related to the parameters of line 1 by
b,(2)= W),
7x2) = Hw-~~N~u)~
a4
cd21
= 4(l),
=
42) = MlW~l)lf(2)1~
f,(2) =fi(l>, ci(1WC1)/f(2)12 4) = a>, ri(2)= ri(lW(lM2)l~
(13)
Thus, the basic parameters are those of line 1, the parameters for line 2 (and any other lines of the multiplet) being derived from these by the Eq. (13). The emergent intensity in both lines is computed by Eq. (lo), thus enabling the function F((parameters)) of Eq. (6) to be calculated. The sum in this equation now extends over both lines. F is again minimized with respect to the parameters. Simultaneous inversion of all the lines in a multiplet is superior to inversion of a single line as it improves the stability of the inversion process and under the right circumstances can yield a larger set of uniquely determined parameters and therefore lead to a better description of the physical properties of the atmosphere. (iii) Determination of systematic velocity jields The presence of stochastic velocity fields has already been implicitly included in the Doppler width. A systematic motion of the atmosphere that varies with depth leads to asymmetries in the emergent line profiles and these asymmetries can be exploited to determine the systematic velocity field. Most of the methods that have been developed for doing this have the following features in common: Doppler shifts are determined at various positions in the complete profile as in the bisector shift method described by Kulander and Jefferies6 or from an analysis of filtergrams as described by Parnell and Beckers;’ the velocities associated with these Doppler shifts are then assigned to particular depths in the atmosphere by means of the Eddington-Barbier relation as described by Kulander and Jefferies or by means of the response function technique of Beckers and MilkeyE and Caccin and Marmolino.’ It is this final step in the analysis where most of the difficulties arise with these methods. Indeed, Parnell and Beckers’ have shown that the assigned optical depths depend on the form of the velocity field itself. This feedback between the depth of formation and the velocity field inherent in these methods can be traced to an
Spectral line inversion as a diagnostic tool
97
attempt to derive the velocity field in a manner which is independent of the other factors that are important in the formation of the line. In contrast, our inversion technique yields all of the parameters that are important in the formation of the line simultaneously and, as a result, no independent determination of the depth of formation, with its associated problems, is required. The immediate effect of introducing a systematic depth-dependent velocity field into a plane-parallel atmosphere is to displace the center of the local absorption profile. This, in turn, results in a modification of both the depth dependence and frequency dependence of the total source function. In our method, it will be assumed that the total source function is frequency independent. This assumption is not strictly necessary and is made purely for convenience. Any modification in the depth dependence of the source function produced by the presence of the systematic velocity field is automatically taken care of in the inversion procedure (to the extent that a “unique” set of parameters can be determined). With this assumption, the emergent intensity is still given by Eq. (1) but the relationship between the total optical depth and the line-center line optical depth given by Eq. (2) is modified. If u(t) is the depth dependent velocity in the atmosphere, then dr, = [9(n + LX)+ r]dt,
(14)
W) = &u(t)lc
(15)
where
is the depth dependent Doppler shift in the center of the absorption profile, with u taken positive towards the observer. The concept of line-center line optical depth loses its usual significance in the presence of a systematic velocity field. The variable t is the line-center line optical depth in a stationary atmosphere and does not correspond to the line optical depth at any wavelength in the observed profile; the wavelength of the observed line center depends on the specific form of the systematic velocity field. The velocity field u(t) is represented by a number of linear segments with segment boundaries which do not necessarily coincide with the source function segment boundaries. With this representation of v(t), Eq. (14) can be integrated to yield zl(t) = G(1, t, (parameters)),
(16)
where (parameters) represent the parameters involved in the piecewise continuous representation of u(t). This relation between r1 and t is then used in Eq. (1) to obtain a relationship between the emergent intensity in the line and the complete set of parameters describing the formation of the line. The inversion procedure is then identical to that already described, except that the parameters determined now include a representation of o(t). In the inversion process it is necessary to evaluate the G-function a large number of times. Efficient methods for its computation are described in the Appendix. ‘RI test the velocity inversion procedure, synthetic line profiles were constructed using a piecewise linear source function and various forms for the synthetic velocity field. The inversion procedure was then applied to these synthetic profiles. In the first series of tests, synthetic profiles were constructed in which the synthetic velocity field was represented by two linear segments with the segment boundary differing from the source function segment boundary. Different choices of the starting values of the parameters in the inversion procedure led to derived parameters in good agreement with those used to construct the synthetic profiles; which confirms that the inversion technique is capable of recovering a synthetic velocity field which has the same mathematical form as that assumed in the inversion process. In a second series of tests,‘synthetic velocity fields in which the velocity increased in proportion to the square of the continuum optical depth z, up to a certain limiting value of z, and then remained constant thereafter, were used to construct synthetic profiles. Two such synthetic profiles were constructed, one with v = 3rC2,z, < 4; u = 48 km/s, z, > 4, and
A. KAVETSKY
98
and B. J. O’MARA
1 ........LINE
BISECTOR
-
INPUT VELOCITY
---
INVERSION
1
VELOCITY
VELOCITY
! ,’
-0
0.5
1 .o
CONTINUUM
1.5 OPTICAL
2.0
2.5
DEPTH
Fig. 2.
the other with v = 12rc2, r, < 2; v = 48 km/s, r, > 2. For this range of velocities, two velocity segements were uniquely specified by the inversion. Although the quality of the fit to the synthetic profiles was not as good as that obtained in the first series of tests since the sum of the squares of the relative deviations F was a little larger, the recovered values of the parameters not associated with the velocity field were in good agreement with the parameters used to construct the synthetic profiles. This result indicates that the recovery of these parameters is not seriously affected by representing a quadratically varying velocity by two linear segments. Figure 2 compares the synthetic input velocity with the velocity recovered by inversion. Clearly the synthetic velocity is quite well represented by two linear segments and is accurately recovered by the inversion. Finally the stability of the inversion procedure was tested by adding synthetic noise, at a level of 5% of the total intensity, to the synthetic profiles prior to inversion. Although, as expected, the quality of the inversion was degraded no evidence of numerical instability was found. To illustrate the advantages of the inversion procedure over the line-bisector method the velocity field was determined by the latter method from the synthetic profiles, using the procedure suggested by Kulander and Jefferies. 6 The results are displayed in Fig. 2. Because the Eddington-Barbier approximation is used, the line bisector method is restricted to depths r, 5 1, and when v = 127,*, where the velocity varies rapidly, it fails completely. These limitations were anticipated by Kulander and Jefferies. 4.
CONCLUSlONS
Although the inversion technique has been demonstrated by us only for synthetic data, examples of its application to the inversion of solar line profiles can be found in Hearn and Halt’ and Kavetsky and O’Mara,” where the method is used to deduce physical properties of the solar atmosphere. It is hoped that the method will be equally useful in the analysis of laboratory plasmas. Acknowledgements-We would like to thank J. N. Holt for several helpful discussions on the inversion method, J. Cooper for useful comments on the manuscript, and the staff at the Joint Institute for Laboratory Astrophysics for assistance in its preparation. One of us (A.K.) was the holder of a Commonwealth Postgraduate Research Award during the course of this work, for which he is grateful. This work was supported in part by funds from Australian Research Grants Commission.
Spectral
line inversion
as a diagnostic
tool
99
REFERENCES 1. J. N. Halt, Ph.D. Dissertation, University of Queensland, St. Lucia, Queensland, Australia, published). 2. A. G. Heam and J. N. Holt, Astron. Astrophys. 23, 347 (1973). 3. M. J. D. Powell, Cornput. J. 7, 155 (1964). 4. D. Mihalas, Stellar Atmospheres, 2nd Edn. Freeman, San Francisco (1978). 5. R. G. Athay, Radiation Transport in Spectral Lines. Reidel, Dordrecht (1972). 6. J. L. Kulander and J. T. Jefferies, Astrophys. J. 146, 194 (1966). 7. R. L. Pamell and J. M. Beckers, Solar Phys. 9, 35 (1969). 8. J. M. Beckers and R. W. Milkey, Solar Phys. 43, 189 (1975). 9. B. Caccin and C. Marmolino, Astron. Astrophys. 83, 73 (1980). 10. A. Kavetsky and B. J. O’Mara, Solar Phys. (submitted). 11. A. Curtis and M. R. Osborne, Comput. J. 9, 286 (1966). 12. C. T. Fike, Computer Evaluation of Mathematical Functions. Prentice-Hall, New York (1968). 13. A. Ralston, A First Course in Numerical Analysis. McGraw-Hill, New York (1965).
1972 (un-
APPENDIX Relationship of optical depths in the presence of a velocity field For the sake of simplicity, the following analysis is confined to an atmosphere with a single linear velocity segment and two source function segments. The results are easily extended to an atmosphere with any number of velocity and source function segments. For a single linear velocity segment, the Doppler shift in the local absorption profile defined in Eq. (15) can be written as a(t) = a,+ agt.
(Al)
the explicit form of the The G-function, defined by Eqs. (14) and (16), can then be determined by substituting absorption profile [Eq. (3) along with Eq. (Al)] into Eq. (14). In segment one where t < T2, collisional damping is assumed negligible so that the a-value in the Voigt profile can be set to zero. The expression for the G-function obtained by integrating over t for the segment is then
for all values of t < T,. When act/d, is small, numerical calculations using this equation are inaccurate due to loss of significant figures in the subtraction of the two nearly equal error functions. In this situation, it is convenient to approximate the difference by a derivative to get G(I., t; r, A,, a,, aJ=rt+texp[-r$+zy]. In segment
number
two,
(A3)
t > T2,
G(X, t; r, A,, AZ, a*, a,, ag) = r,T, + rz(t + where H is the Voigt function
A2
a,sH(a,,
(A4)
0)
of Eq. (3) and the X-function
is defined
by
x(c, a)=jymarctan(c-i)exp(-y2)dy.
(A5)
Some of the mathematical properties of ~(c, a) are given in Table 1. For similar reasons to those leading (A3), the difference of X-functions in Eq. (A4) may be approximated by a derivative to give
to Eq.
WA, t; r, A,, A,, a,, a,, q)=r,T,+r,(t
(‘46) when a&t - TJA,a, is small. The error function can be easily and rapidly calculated using the rational approximation of Curtis and Osborne.” To permit rapid computation of the X-function, ~(c, a)/c was represented (for fixed values of a) by a Chebyshev polynomial expansion in c for small and moderate values of [cl, i.e., for 0 < (cl < c0 for some moderate value of c,; c, decreases with increasing a and is N 70 for typical values of a. For values of ICI greater than co, where x(c, a) approaches its limit f7?/2, it is convenient to introduce the function fi(k, a) defined by
k
(A7)
KAVETSKY
A.
100 Table 1. Some properties
of the functions
and B. J. O’MARA
x(c, a) and Q(k, a) defined by Eqs. (A$
. _ Expression
Property
FUIlCtiOll
07
(A7), and (A8).
Value
3/2 x(c.d
alternative for c > 0
expression x ev(-y2)dy aTI erfc(a)
derivative
with
t0 C, i.e., x(c/a,a)
limit
S2Ck.a)
alternative for k > 0
as
respect
exp(a2)
a,,H(a,ac)
ax/at a +
0
(7 3/2/2)
erf Cc)
312
expression
+
erfc
($
+ ;
x
n(k,a)
limit
as
k +
n(k,a)
limit
as
k + _
0
,/”
arctan[
(i
- 3
-l I
ew-y2)dy
.1/2
0
where k = l/c.
(AS)
This function was approximated for the range 0 < Ik( < l/ c,,, i.e., c,, i ICI < co by a Chebyshev polynomial expansion as a function of k for fixed values of a. Representing x(c,a) in this way ensures an accurate approximation to this function for all values of c at a given value of a including very small and very large values. The techniques and equations used in the calculation of the Chebyshev coefficients are given in Fike.12 Chebyshev coefficients for each of these functions are available from the authors for a range of a values. With the above expressions for the G-function and the Chebyshev expansions allowing its rapid calculaton the relationship between the optical depth at any wavelength and the line-center line optical depth given by Eq. (16) is fully defined. In order to do the numerical integration for the emergent intensity (Eq. [l]) it is necessary to solve Eq. (16) for t given a value of ri. To do this, Eq. (16) is written in the form TV- G(I, and solved by using the secant
method
as described
t, parameters) by Ralston.”
= 0
(A9)