Volume 76, number 1
OPTICS COMMUNICATIONS
1 April 1990
SPECTRAL CHARACTERISTICS OF STEADY-STATE A M P L I F I E R S
Mark D. ROTTER and Roger A. HAAS Department of Applied Science, University of Cal(/brnia, Davis/Ltvermore L- 794, P.O. Box 808, L~vermore, ( ;.t 94550, ( "A'.4 and Lawrence Livermore National Laboratoo,, Livermore, USA Received 25 September 1989
The theory of steady-state optical amplification for homogeneously broadened media is extended to include broad spectral bandwidth input signals. The output power gain, amplifier efficiency and spectral narrowing of the amplified radiation is studied. In the small and large signal regimes, approximate analytic expressions are obtained and compared with numerical results.
In a recent paper [1 ] the spectral characteristics of short-pulse amplifiers were explored in the context of generating high-power, partially coherent light for fusion [2,3] with tunable solid-state laser media such as Ti:A1203. Such light has been shown to increase the laser light/target coupling in inertial-confinement fusion (ICF) experiments [4]. However, other media with broad fluorescence bands, notably organic dyes [ 5 ] and excimer lasers [ 6 ] such as KrF [2,7-9] may prove useful in generating such radiation. In these media, the fluorescence lifetime may be less than the pulse duration of interest with the result that optical amplifiers would nominally operate in the steady-state regime. The steady-state amplification characteristics [ 10-12] of homogeneously-broadened media, including excimer [8,13 ] and dye laser amplifiers [14-17] have been extensively studied. These studies have generally centered on amplification of narrow-band radiation or the spectral characteristics of noise generated within a laser amplifier. In this paper, the steady-state amplifier theory [8,10-17] is extended to include nonmonochromatic incident radiation. The amplified radiation power gain, amplifier efficiency and spectral narrowing are studied by numerical solution of the governing equations. In addition, approximate analytic expressions are obtained for these quantities in the small and large-signal limits and compared to the numerical results. Consider a light amplifier containing an opticallypumped, homogeneously-broadened laser medium 56
whose fluorescence lifetime is short compared to the laser pulse to be amplified. In this steady-state mode of operation the pump pulse is coincident with the signal to be amplified. This light signal is injected into one side of the medium and is of sufficient strength so that the forward and backward-traveling amplified spontaneous emission signals may be neglected. This assumption ultimately limits [ 8,13, 17 ] the minimum input intensity and maximum amplifier aperture for which the results obtained here apply. The limitations imposed by noise buildup are conventionally overcome by using bidirectional amplifiers [8,13]. The medium is presumed to have a broad fluorescence (gain) bandwidth. It is also assumed that the relaxation rate out of the lower lasing level is much greater than any rate into it so that no build-up of population occurs in this level. In this model, excited-state absorption as well as self-absorption of the radiation is neglected. The host medium is non-dispersive with refractive index n. The amplification of broad-spectral-bandwidth, plane wave, light pulses is governed by [ 1 ]
(n/c)Ori:(z, t, ),) +O=i:(z, t, 2) =a,(),)Nz(z,
l ) i : ( z , l, 2 ) - o ' t : ( : ,
t. 2) .
(l)
and O,N2( Z, t ) =c~a( Up ) No( z, t ) l . ( t ) / hvp t
[ ) / t 2 - N 2 ( z , [) ~ 2a,(~.)i/(z, t. 2 ) d 2 / h ( . (2)
0030-4018/90/$03.50 © Elsevier Science Publishers B.V. (North-Holland)
Volume 76, number 1
OPTICS COMMUNICATIONS
where it(z, t, 2) and N2(z, t) are the spectral light intensity (intensity per unit wavelength interval) and upper lasing level population density, respectively. The physical constants c and h are the speed of light in vacuum and Planck's constant. The stimulated emission cross-section ors(2) is assumed to be a gaussian function of wavelength 2 with center wavelength 2o and a full-width at half-maximum (fwhm) of A2F. The integral in eq. (2) ranges over the high frequency gain band of the medium. The monochromatic optical pump intensity is Ip(t) and ~a(Vp) is the absorption cross-section at the pump frequency vp. Non-saturable photoabsorption losses in the host medium are characterized by a and t2 is the lifetime of the upper lasing level. The quantity No is the population density of the ground state and for dye media, No (z, t) + N2 (z, t) = N, where N is the total concentration of lasing species (a constant). This form of pumping in eq. (2) and ground-state depletion is applicable to dye laser media. However the non-dimensional results presented in this paper also apply to other steady-state laser media, such as excimer lasers (e.g. KrF) which are not optically pumped. Specific changes in the analysis required for excimer lasers are discussed later. In steady-state, eqs. (1) and (2) reduce to the single equation
Ozit(z, 2) = [g(z, 2) --ollil(z, 2) ,
(3)
where
g(z, 2 )=go(2) / ( l + ( l + t2/ta)-'t2
× f 2~,(2)i,(z,2)da/hc)
(4)
is the space and wavelength-dependent gain coefficient. In eq. (4), ta represents the effective pump absorption time given by t~= hvp/a~(Vp)Ip. The smallsignal gain coefficient is go(2)=es(2)N(tz/ta)/ (1 + t2/ta). This steady-state amplification problem is fully specified by eq. (3) together with the boundary condition: iA0, 2 ) = i o ( 2 ) . In addition to the spectral intensity, the light is also characterized by the total intensity
I(z) = f it(z, 2)d2.
(5)
The performance of an amplifier of length L may be characterized by the radiation power gain,
1 April 1990
G=I(L)/I(0), and the spectral narrowing fraction, q= (A2o-AAL)/A2o, where A2o and AAL are the initial and final radiation fwhm spectral widths. Since the laser radiation experiences more gain at line center, its spectrum narrows during amplification: A2o>A2L and ~/> 0. The parameter q gives the fractional decrease in the radiation spectrum due to amplification. Eqs. ( 1 ) and (2) were numerically integrated using the method of lines [ 18]. The pump intensity had a hyperbolic tangent temporal dependence with the final time sufficiently large to ensure steady-state behavior. All pertinent quantities were evaluated at the final time point. This same temporal dependence was used for the input signal. In the steady-state regime, the input spectral intensity was characterized by/0(2) =Iofo(2) where fo(2) is a normalized gaussian distribution with center wavelength 2o and a fwhm of A2o. Consequently, from eq. (5), Io is the total input intensity. In general, it is found that dye laser amplifier performance depends on the dimensionless parameters A2o/A2F, A2F/2O, Io/Is, go~a, and goL, where Is- hvo ( 1 + tEl ta )/¢7s(20) t2 is the saturation intensity and go - as (2o) N( t2/ ta ) / ( 1 + t2/ ta ) is the small-signal gain coefficient evaluated at linecenter. The quantity Vo is the line-center radiation frequency. For these calculations, 22F/2O = 0.14, as is appropriate for dye media such as LDS-698 and DCM [ 19], and t2/ta = 0.01. In general, these results are relatively insensitive to 22v/2o and are approximately correct for excimer lasers provided go = o's(2o)Rpt 2 and Is-hvo/as(2o)t2, where Rp is the excimer molecule production rate per unit volume [6 ] and t2 is its lifetime. It is important to note that for dye laser amplifiers the saturation intensity Is depends on the pump intensity e.g. when t2>>ta or Ip>>hvp/• a( up ) t2 then Is~ [hvo/as(Ao)t2] (ta/ta)
(])O/Vp) [ O'a(2p)/O's (20) ]Ip. In principle this novel circumstance permits dye laser amplifiers to operate at the intensities of interest for fusion. However, it is likely that nonlinear processes, such as Raman scattering from solvent molecules, will complicate this simple picture. Figs. 1 and 2 show the dependence of power gain and spectral narrowing fraction on the normalized 57
Volume 76, n u m b e r 1
OPTICS COMMUNICATIONS
. . . .
L
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
i<,ii~=lo ~
. . . .
. . . .
I April 1990
. . . .
. . . .
. . . .
. . . .
. . . .
1
7 ~"~"~"~'" J
lai
["-"a~Q---~-"q"Qf/)(]s~+
-1
~ I-
hi
-i
I u/1,.02-I( ~-"N'-'~
F
r,,,,I .... I,,,tl,,,,I,,ttl,,,,i,,,,i,,,,i,,,,] .1 .2 .3 .4 .5 .6 .7 .8
•I
#
~lJ
.7
.8
(1 .9
I.(/ .t
.2
.3
.4
.5
.6
.t]
1.0
AL o , AX.I.
A~ o i A L F
Fig. 1. ( a ) Steady-state a m p l i f i e r power gain G for go/+ = 2, a = 0, and A2v/2o= 0.14, as a function of n o r m a l i z e d i n p u t r a d i a t i o n spectral b a n d w i d t h Lkko/ZMv. The curves are p a r a m e t r i c in n o r m a l i z e d i n p u t intensity 1o/1~.Solid lines indicate n u m e r i c a l solutions. D a s h e d lines indicate a p p r o x i m a t e analytic solutions for small-signal (lo/1~= 10-4: --- eq. (7) ) a n d large-signal (e.g. 1o/I~= 1; . . . . . eq. ( 19 ) ) regimes. ( b ) . Steady-state a m p l i f i e r spectral n a r r o w i n g fraction ~, for goL = 2, a = 0, and A2v/2o= 0.14, as a function of n o r m a l i z e d input r a d i a t i o n spectral b a n d w i d t h A2o/A2v. The curves are p a r a m e t r i c in n o r m a l i z e d i n p u t intensity lo/1~. Solid lines indicate n u m e r i c a l solutions. D a s h e d lines i n d i c a t e a p p r o x i m a t e analytic solutions for small-signal (Io/I~ = 10- 4: ___ eq. ( 10 ) or • eq. ( 9 ) ) and large-signal (e.g. Io/1~= 1; . . . . . eq. ( 2 4 ) or O eq. ( 2 3 ) ) regimes.
6(1~.... i .... m.... r .... m.... m.... m.... T.... m. . . . .
a
,,j, :,0+
t~
'°F ~ .E
t .... m.... m.... J .... ~.... m.... m.... r .... m''"t-('
,J,
tl
.Ot
+~/~
"-.
~'t
I
.2
: ,
""'''''
~ 3 0 . ~ - ' - - i i ] ~ 2 0
5
+~
~'
2 1.? i
t I. /2. i5. .~
.4
.5
.6
.7
.8
.9
"~)~0/ \'AI
1.0 .1
,2
.3
.4
.5
.(~
7
.N
.9
I+(I
\)'o / \P'I
Fig. 2. (a) Steady-state a m p l i f i e r power gain G and ( b ) spectral n a r r o w i n g fraction r/as a function of n o r m a l i z e d input r a d i a t i o n spectral b a n d w i d t h AXo/A2 F. P a r a m e t e r s are as in fig. 1 except g o L = 4 .
input radiation spectral width for low (goL= 2 ) and high (goL = 4) gain amplifiers respectively under no photoabsorption-loss conditions, o~ = O. Part (a) of
58
each figure gives the power gain as a function of A2o/ A2v and part (b) shows the spectral narrowing fraction as a function of Lk,l.o/A)+v. In each figure, the
Volume 76, number 1
OPTICS C O M M U N I C A T I O N S
curves are parametric in normalized input intensity Io/Is. For clarity, the spectral narrowing fraction curves for lo/ls=O.O05 and 0.01 have been deleted. As the input radiation bandwidth A2o/A2v is increased, the gain is reduced significantly for small input signals (Io/Is<
~ t ) the amplifier gain is relatively insensitive to the magnitude of the input radiation bandwidth. In this regime, spectral narrowing of the amplified radiation is also much reduced compared to the small-signal regime. However, for high gain amplifiers, fig. 2b, spectral narrowing in the large signal regime is significant when the input bandwidth approaches the fluorescence bandwidth. In figs. 3 and 4 the power gain and spectral narrowing fraction dependence on input radiation spectral width are plotted under the condition go/a= 10 for goL = 2 and 4, respectively. As above, part (a) of each figure gives the power gain as a function of A2o/ A2F and part(b) shows the spectral narrowing fraction as a function of Aito/A2v. As in the case for a = 0 , for clarity the spectral narrowing fraction curves for Io/Is = 0.005 and 0.01 are not shown. The power gain
6[~,o/'is 5 2..05
is reduced when loss is introduced. In the small-signal regime spectral narrowing is relatively insensitive to loss, however, the amount of spectral narrowing increases as loss is introduced in the large signal regime. This occurs because when non-saturable loss is present the small-signal regime, where the amount of spectral narrowing is greatest, is lengthened. Both the degradation of power gain and increased spectral narrowing for large signals are greater for the high gain case. The results presented in figs. 1-4 required detailed numerical calculations. Fortunately, approximate analytic solutions can be construtted that cover a wide range of interesting physical conditions. In the small signal regime eq. (3) can be solved for the spectral intensity
il(L, 2) =io(2)exp[ (go(2) - o r ) L ] . Using eqs. (5) approximately
and
(6)
(6), the power gain is
G = e x p [ (go-a)L] [1 + (A2o/A2v)ZgoL]-'/2.
(7)
Furthermore, if AAL is the spectral fwhm of the radiation at the end of the amplifier, then
it(L, )~o+_AAL/2)/it(L, 2o) = 1 / 2 .
(8)
Substituting eq. (6) into eq. (8) produces a tran-
'104' ......... ' .... t
.5
J
.4
"
-
1 April 1990
4
i
~ 3
'~
22---1.- . . . . . . . . . .
o
ii.........i?i .2
.3
.4
.5 .6 .7 A~.0 / A~.F
iI .8
.9 1.0.1
.2
.3
.4
.5 .6 .7 A~.0 / A~.F
.8
.9 1.0
Fig. 3. (a) Steady-state amplifier power gain G and ( b ) spectral narrowing fraction ~/as a function of normalized input radiation spectral bandwidth AAo/AAF.Parameters are as in fig. 1 except g o / a = 10.
59
Volume 76, n u m b e r 1
40
OPTICS COMMUNICATIONS
....
I ....
I ....
~i~/]
I ....
I ....
i ....
I ....
I ....
I ....
= ]H-;
....
1 April 1990
['TT'Ir'''I .... I .... E. . . .
t
I
~ ....
ln/I
.... I .... / .6 = l( .~
.1 ,"
~ 2 0 ''-.02
"''.
,3 ,£
~
'"
.3 '~
.
lO - J
o .l
~
.-~
q
.2
.3
.4
.5
I
.6
.7
.8
.9 1.0 .
.2
.3
A~ 0 AL F
.4
.5 ,,'t_k o i
.6
.7
.8
.9 1.0
A~. F
Fig. 4. (a) Steady-state amplifier power gain G and (b) spectral narrowing fraction q as a function of normalized input radiation spectral bandwidth A2o/A2v.Parameters are as in fig. 1 except goL= 4 and go/a = IO. scendental equation for x = (/~L//~k~0) 2, Expanding this equation to second order in x and solving the resulting quadratic equation yields (Zk2e/Zk,~o) = { [ b - ( b 2 - 4 a ) , / 2 ] / 2 a }
,/2
(9)
where b = 1 + ( ~ o / ~ v ) 2 g o L and a=ln2(A2o/A2F) 2 ( b - 1 )/2. In the limit that A2o<
(A2o/ZK;W.)22goL ] -l/2
(10)
In the small signal regime, the power gain and spectral width o f the radiation are independent of the pulse input intensity, and the output spectral width is independent of loss. To obtain an approximate analytic solution for the large signal regime, suppose the radiation spectral intensity is taken to have the form i/(z, 2 ) = S z , 2 ) l ( z ) .
(11)
In order to identify I(z) in eq. ( 11 ) as the intensity of the light, the function f ( z , 2) must be normalized to unity. This will be shown in the following. Substituting eq. (11) into eqs. (3) and (4) and integrating over wavelength for (zk,~o/2o) 2<< 1 produces
60
l(z)3zF(z)+F(z)Ozl(z) =go(z)F(z)l(z)/[l + F ( z ) I ( z ) / ~ ( z ) ] -~F(z)l(z), where F(z)=Jf(z, 2)d2, ~ ( z ) =hvo( 1 +l")/~s(Z)t2,
#s(z)=
(12)
go(z)=ds(z)Nr/(l+z). z=t2/t., and
f G(A)f(z,A)d2/' f f(z,,~)d2
(13)
is the local radiation-spectrum-averaged value of the stimulated emission cross-section. The quantity l ( z ) is now required to satisfy
3zlnI(z)=go(Z)/[l+F(z)I(z)/I-~(z)]-a.
(14)
Then from eq. ( t 2 ) , 0 : F ( z ) = 0 , or F(z) =(', a constant. Taking this constant equal to unity identifies I ( z ) in eqs. ( 11 ) and (14) as the intensity. With this, eq. (14) can be written in the form
3zl(z) =go(z)l(z)/[ 1 +l(z)/[~(z)] - a l ( z ) , ( 15 ) where go(Z) a n d / d z ) are the local radiation-spectrum-averaged small-signal gain coefficient and saturation intensity, respectively. For the case of discharge of e-beam pumped excimer lasers such as KrF, eq. (15) remains valid provided go(Z)=G(z)Rd2
Volume 76, number 1
and [~=hvo/#s(z)t2. Through #s(z), the equation for the intensity I(z) depends o n f ( z , 2). To obtain an equation f o r f ( z , 2), eqs. ( 3 ), (4), and ( 11 ) may be combined to give Ozlnf(z, 2) = [as(2)/a~(z) - 1 ] [Ozln/(z) + c e ] , (16) which is also consistent with F(z) = C. Eqs. (15) and (16) form a complete set within the ansatz made by eq. ( 11 ). They are valid in both the small and large signal regimes provided (A,,],o/2o)2<<1. Their solution is fully specified by the boundary conditions: I(0) =Io and f ( 0 , 2) =fo(2). Provided goL and A2o/A2F are not too large, in the large-signal regime an approximate solution to this coupled set of equations can be readily obtained by takingf(z, 2) =fo(2). Then eq. (13) becomes
8s(z) ,,~os-f as(Z)fo(2) d2
(17)
or
o~=a~(2o)/[ 1 + (A2O/A2F) 2] 1/2
(18)
Since #~ is independent of z, eq. (15) can be integrated. The power gain, G, is then obtained from the transcendental equation
~ln G - l n [ ( l - ~ ( l + uoG) ) / ( 1 - ~ ( l + uo) ) ] =goL~( 1 - ~ ) ,
(19)
where ~=ot/go, and Uo=Io/[~. For dye laser media ~o=6sN(t2/ta)/(l+t2/ta), and ~=hvo(l+t2/ta) /~Yst2; whereas for excimer lasers ~o=OsRpt2 and
[s=hvo/@st2. In the limit o~--,0, eq. (19) reduces to In G+ uo(G- 1 ) =goL.
(20)
The approximate spectral distribution of the radiation may be obtained from eq. (16),
f(L, 2) =fo(2) {G}~(a)exp [7(2)c~L] ,
(21)
where the exponent 7(2) = [as(2) -~]/~rs. If A2L is the spectral fwhm of the radiation at the end of the amplifier then eq. (8) becomes
f(L, 2O+_A2L/2)/f(L, 2o) = 1/2.
(22)
Following a procedure analogous to that used to derive eq. (9) it can be shown using eqs. (21 ) and (22) that
(A2L/A2O)={[B-(B2-4A)~/E]/2A} w2 ,
1 April 1990
OPTICS COMMUNICATIONS
(23)
where
B = 1 + ( ~ 0 / / ~ F ) 2 [ 1 + ( / ~ 0 / / ~ F ) 2 ] 1/2
ln(Ge '~L) and A=In2(A2o/A2F)2(B - 1 )/2. In the limit that A2O<
(24)
X (A2O/AgF)Zln(Ge'~L)} -1/2 ,
where, in general, G depends on the normalized input intensity Io/~; eq. (19). Thus, in the small-signal regime, eqs. (7) and (9) may be used to calculate the power gain, and spectral width of the radiation respectively. Eqs. (19) and (23) give the power gain, and spectral width in the large-signal regime. Comparison of these approximate formulas for laser power gain and spectral width with the numerical results presented in figs. 1-4 is good. In the small-signal regime, the disagreement between numerical integration ofeq. (6) and the numerical solutions may be traced to the 2% accuracy of the latter. In general, as expected, the approximate solutions weaken as the input radiation spectrum width approaches the fluorescence width and at high gain. The results presented in figs. 1-4 cover the input radiation spectral range 0.1
--1
(N
(25)
o
In the large-signal regime, the integral in eq. (25) can be transformed to an integral over intensity using eq. ( 15 ). Using eqs. (2), ( 15 ), and ( 19 ) in eq. (25) yields for dye laser media r/amp= (J-v/go) Uo( 1 +Z) ( G - 1 ) × [goL(1 + z ( 1 - ~ ) ) - z l n G ] - '
(26)
The quantity G is the large-signal power gain, as given 61
Volume 76, number 1
OPTICS COMMUNICATIONS
by eq. (19). In the limit a ~ 0 , eq. (26) takes on the form
?]amp = (/]-p/~'O)(1 "1-r) b/o(G-- 1 ) × [goL+ruo(G- 1)]-I
1 April 1990
This research was performed in part under the auspices of the U.S. D e p a r t m e n t of Energy under contract No. W-7405-ENG-48.
(27)
References where eq. (20) has been used. For excimer lasers the amplifier efficiency is given by
?]amp= [ I ( L ) - I o ] / ~ R p L = (hvo/~) u o ( G - 1 ) / g o L ,
(28)
where ~ is the energy required to make the excimer molecule. The amplifier efficiency values given by eqs. ( 2 5 ) - ( 2 8 ) are in good agreement with those obtained numerically. For experimentally-obtainable values of i n p u t width ( A ~ 0 / ~ F < 0 . 5 ) , the extraction efficiency values as given by eqs. ( 25 ) - ( 28 ) are within 5% of the numerical values for both goL = 2 and 4. In summary, the steady-state theory of laser amplifiers [ 8,10-17 ] has been extended to describe the amplification of broad spectral b a n d w i d t h light. The laser power gain, amplifier efficiency, and spectral width were investigated. It was found that these characteristics are sensitive to i n p u t radiation spectral width in the small-signal regime. In the large signal regime, they are relatively insensitive to the spectral width except in the case of large gain and i n p u t spectral width. Approximate analytic expressions were obtained in both the large and small-signal regimes. These expressions provide good approximations to the numerically obtained solutions. The analytic and numerical results should prove useful in the study of broad spectral b a n d w i d t h light amplifiers that operate in the steady-state limit.
62
[ 1] M.D. Rotter and R.A. Haas, Optics Comm. 71 ( 1989) 311. [2] R.H. Lehmberg and J. Goldhar~ Fusion Technology 11 (1987) 532. [3] R.H. Lehmberg and S.P. Obenschain, Optics Comm. 46 (1983) 27. [4] A.N. Mostovych, S.P. Obenschain, J.H. Gardner, J. Grun, K.J. Kearney, C.K. Manka, E.A. McLean and C.J. Pawley, Phys. Rev. Len. 59 (1987) 1193. [5] F.P. Sch~ifer, ed., Dye Lasers, 2nd Ed., Topics in applied physics, Vol. 1 (Springer, Berlin, 1978). [6 ] J.J. Ewing,Excimerlasers, in: Laser Handbook, Vol. 3, M.L. Stitch, ed. (North Holland, Amsterdam, 1979) p. 135. [7] J.J. Ewing, R.A. Haas, J.C. Swingle,E.V. George and W.F. Krupke, IEEE J. Quantum Electron., QE- 15 (1979) 368. [ 8 ] R.A. Haas, W.L. Morgan, D. Eimerl, A.E. Siegman,A. Szoke and C. Turner, Laser Program Annual Report 1979, LLNL, UCRL-50021-79, 3, 7-2 (1979). [9] D.J. Dudziak, ed., KrF Lasers For Inertial Confinement Fusion, in: Fusion Technology 1l ( 1987). [ 10] L.W. Casperson and A. Yariv, IEEE J. Quantum Electron. 80 (1972) QE-8. [ 11 ] L.W. Casperson, J. Appl. Phys. 48 (1977) 256. [ 12] L.W. Casperson, Optics Comm. 8 (1973) 85. [13]A.M. Hunter and R.O. Hunter, Jr., IEEE J. Quantum Electron. QE-17 (1981) 1879. [ 14] R.S. Hargrove and T. Kan, IEEE J. Quantum Electron. QE16 (1980) 1108. [ 15 ] T. Urisu, K. Kajiyama and Y. Mizushima, J. Appl. Phys. 48 (1977) 4168. [ 16] G. Haag, M. Munz and G. Marowsky, IEEE J. Quantum Electron. QE-19 (1983) 1149. [ 17 ] U. Ganiel, A. Hardy, G. Neumann and D. Treves, IEEE J. Quantum Electron. QE-I 1 (1975) 881. [ 18] P.N. Brown and A.C. Hindmarsh, SIAM J. Numer. Anal. 23 (1986) 610. [ 19] P.R. Hammond, Optics Comm. 29 (1979) 331.