We derive analytic expressions for the Brillouin thresholds of square pulses in optical fibers. The equations are valid for pulse durations in the transient Brillouin scattering regime (less than 100 nsec), as well for longer pulses, and have been confirmed experimentally. Our analysis also gives a firm theoretical prediction that the Brillouin gain width increases dramatically for intense pulses, from tens of MHz to one GHz or more.
© 2014 Optical Society of America
Stimulated Brillouin scattering is one of the principal nonlinear processes that affect optical fiber systems [1, 2]. It occurs when narrow bandwidth light guided in an optical fiber interacts with the fiber itself to produce a density grating that backscatters the light. The stimulated Brillouin wave grows exponentially and, if the intensity of the initial beam is greater than a certain threshold value, the Brillouin scattering depletes most of the light from the beam.
The Brillouin threshold presents an enormous obstacle to transmitting high-power, single-frequency pulses in a fiber. Pulsed fiber lasers and amplifiers such as those used for micro-machining  achieve peak powers of tens of kilowatts or more, and when the pulse duration is greater than about one nanosecond, the Brillouin threshold puts a firm upper limit to the power that can be delivered.
Until now, however, there has been no easy way to calculate the Brillouin threshold intensity for a pulse in a single mode fiber. By contrast, the Brillouin threshold intensity IP for a continuous (c.w.) pump beam is simply related to the Brillouin gain gB, the fiber length L, and a threshold parameter Θ by the formula :
This equation is accurate when the pulse is longer than about 100 nsec, but it severely underestimates the Brillouin threshold for shorter pulses. The reason is that there is a time scale relevant to Brillouin scattering, the phonon lifetime TB, that causes substantial deviations from Eq. (2) when τ ≤ ΘTB.
Relatively few authors have investigated Brillouin scattering at such short time scales (known as transient Brillouin scattering). Some studies pertain to small interaction cells that are much shorter than the pulse length [2,4,5], the limit opposite to the one relevant for fiber optics. Other papers look at the Brillouin scattering of pulses in optical fibers, but use numerical integration [6, 7] or experiments  to obtain results.
None of these authors give an elementary equation analogous to Eq. (2) for the Brillouin threshold for a pulse in fiber. One of the goals of this paper is to provide such an equation, relevant to the high-powered nanosecond pulses produced in fiber laser systems today. Such pulses typically have peak intensities greater than 10 W/μm2, and the fibers used to deliver these pulses are typically less than 10 meters in length. We consider only Brillouin scattering in passive fiber; our results would be modified somewhat for fibers with gain, such as fiber lasers and amplifiers.
We have found that under three assumptions, it is possible to derive an algebraic expression for the Brillouin threshold of a pulse in passive fiber. The assumptions are, first, that the pump pulse is square (that is, it has negligible rise and fall times, and the body of the pulse has constant intensity); second, that a negligible amount of power is depleted from the pump pulse; and finally, that the round trip time in the fiber (2L/υ) is greater than the pulse duration τ. Given these three assumptions, we have derived a simple equation, Eq. (14) below, that generalizes Eq. (2) and is valid for all pulse lengths.
The second assumption, the undepleted pump approximation, is valid as long as (1) fiber losses are negligible, which is the case for the short lengths of fiber used for high power pulse delivery, and (2) the equations are not used to model the pulse after it reaches the Brillouin threshold. This second limitation does not affect our ability to calculate the Brillouin threshold itself, however, and this threshold is the most important quantity to predict, since pulses beyond the threshold are generally too noisy for their intended application (such as material processing).
Extending our analysis, we have found that the Brillouin gain width for high intensity pulses is significantly greater than the narrow width usually associated with Brillouin scattering. This result is presented in Section 3. We have also found that when the third assumption above is relaxed—that is, when relatively short fibers are used—we can still find an analytic solution to the differential equations governing Brillouin scattering, though the resulting expression for the Brillouin threshold, Eq. (38), requires numerical root-finding techniques to evaluate.
Finally, we have conducted experiments to test these predictions. To compare the experiment and the theory, accurate values of the Brillouin gain gB and the phonon lifetime TB were needed. These values depend on the composition of the fiber, and were not known in advance; we therefore fit our equations for the Brillouin threshold to the data, treating gB and TB as free parameters. The best fit values for gB and TB are similar to the values found for other fibers, and the results are shown in Fig. 1.
2. Brillouin threshold in long fibers
We begin with the fundamental equations for the Brillouin scattering of a square pump pulse with duration τ. As the pulse travels a distance z for a time t in the fiber, the envelope of the pulse is described by an amplitude function AP(z, t). In the undepleted pump approximation, the pulse propagates with group velocity υ without changing its shape:1]:
The amplitudes AP and AB are normalized so that |AP|2 and |AB|2 are the instantaneous powers of the pump and Brillouin pulses. In Eq. (5), a term describing the propagation of phonons (proportional to ∂Q/∂z) has been omitted, because the phonons decay before they can travel a significant distance .
The above equations are difficult to solve because the pump and Brillouin pulses travel in opposite directions, and the phonons are essentially stationary. Therefore three different frames of reference are important, and there is apparently no way to describe all three waves simply in a single frame of reference.
However, when the fiber is longer than the pulse and when the pump is undepleted, a simplification is possible. The basic insight is that to an observer that is co-moving with the pump pulse, the Brillouin amplitude at the observer’s position does not change. To see this, consider an observer positioned at the trailing edge of the pump pulse: at any given instant he sees a Brillouin wave that has built up over the time τ that has elapsed since the leading edge of the pulse passed that point in the fiber. This is true no matter how far along the fiber he has traveled; therefore the Brillouin amplitude at the back of the pulse is independent of time (as long as the front of the pulse has not yet reached the exit end of the fiber). Another observer placed in the middle of the pulse would get a similar result. This invariance of the Brillouin amplitude for any co-moving observer can be expressed as an equation similar to Eq. (3):
Equation (7) is the key simplification that makes it possible to derive an algebraic expression for the Brillouin threshold of the pulse. Adding Eq. (7) to Eq. (4), the spatial derivative ∂AB/∂z is eliminated:Eq. (5) to eliminate ∂Q*/∂t from the above equation, Eq. (8) to eliminate Q*: Eq. (6) to eliminate κ1κ2, defining the pump intensity IP = |AP|2/Aeff, and using ΓB = 1/TB:
Equation (9) is an ordinary differential equation that can easily be solved. The solution at a given coordinate z is an amplitude that grows exponentially in time:
To find the complete behavior of AB(z, t) as a function of both z and t, note that Eq. (7) requires that AB(z, t) be a function of t − z/υ. Therefore the complete solution is:
The Brillouin radiation grows exponentially from thermal noise of power |AB0|2 at the beginning of the pulse . (The noise power |AB0|2 is calculated below in section 4.) The Brillouin threshold occurs when the Brillouin power at the back of the pulse is equal to the pump power ; that is, whenFig. 1. It gives the Brillouin threshold for a pulse whose duration is less than the relevant transit time in the fiber. Since the front of the pulse exits a fiber of length L at time t = L/υ, this would seem to require that τ ≤ L/υ. However, Eq. (14) is still valid for greater values of τ, because of the time it takes for the information to propagate that that the front of the pulse has reached the end of the fiber. As long as the back of the pulse enters the fiber before it receives the information that the front of the pulse has left, Eq. (14) can still be used. That is, Eq. (14) is valid for pulse durations less than the round trip time in the fiber, or τ ≤ 2L/υ.
A possible concern about the above approach is that the group velocity of the Brillouin wave is substantially smaller than the group velocity of the pump, due to the narrowness of the Brillouin gain bandwidth . This leads to the objection that perhaps the velocity υ appearing in Eq. (4) should not be just the group velocity of the material, but should include the effects of the Brillouin gain as well. However, it turns out that the correct group velocity in Eq. (4) is indeed the group velocity of the material only; all gain-dependent effects are automatically included in the solution to the equations.Eq. (2), the c.w. Brillouin threshold for an interaction length of υτ/2, as expected.
3. Brillouin gain width
Brillouin scattering is well known to have a very narrow gain bandwidth, on the order of tens of MHz [1, 2]. However, we have found that at high intensitites, the Brillouin bandwidth becomes much larger, on the order of 1 GHz or more.
To see why, we return to Eq. (5) and generalize it slightly. Equation (5) is valid only when the Brillouin wave is exactly on resonance; that is, when the Brillouin optical frequency ωB is equal to the difference between the pump and acoustic frequencies: ωB = ωP − ΩA. To investigate the Brillouin gain bandwidth, we allow ωB to vary, and look at the Brillouin gain at different frequencies. Defining Ω = ωP − ωB, Eq. (5) is replaced by :Eqs. (4) and (7), and following the same procedure as before, we find, in place of Eq. (9),
The solution to Eq. (16) is an exponentially growing wave similar to before:
The Brillouin power at the back of the pump pulse is:
On the other hand, for high intensity pump pulses (γ ≫ 1), Eq. (18) implies that the gain β + β* has a half-width of , orEq. (22), with υ = 0.2 m/nsec, gB = 31 μm2/(W-m), and TB = 4 nsec (see Section 7 for a discussion of these values), we find a half-width of:
This is an important result, because it means that Brillouin radiation is more difficult to suppress than previously thought: Some methods of reducing Brillouin scattering rely on changing the frequency of the Brillouin gain peak by more than the gain bandwidth before the Brillouin radiation can reach threshold. (These methods include introducing a temperature  or strain  gradient.) Given the increased gain bandwidth at high powers, this change will have to be much larger than expected.
Other authors have noted that the Brillouin gain width increases under certain circumstances. For example, the width can increase due to saturation effects  or due to the Fourier transform limit of short pulses [7,13]. However, the width given by Eq. (22) is caused by neither of these, since it occurs in the unsaturated regime (the pump is undepleted), and does not depend on the pulse duration.
What, then, is the intuitive explanation for this dramatic increase of gain width at high powers? The right hand side of Eq. (22) is approximately equal to the peak gain β + β* on resonance (when δ = 0). We define the inverse of this quantity as the gain time Δtg, the time it takes for the Brillouin radiation to increase e-fold; Eq. (22) can then be rewritten:Eq. (21) can be written:
For the next section, we need to calculate the width of the actual Brillouin radiation, ΔωB′, as opposed to the gain width ΔωB. The gain width ΔωB is larger than the width of the Brillouin pulse due to the spectral narrowing that is typical of amplified light. To calculate the final width ΔωB′ of the Brillouin pulse, we approximate the power gain for high power pulses (γ ≫ 1) as:Eq. (11). Substituting Eq. (23) into Eq. (20) and using the definition of δ, we find that:
4. The threshold parameter Θ
The value of the threshold parameter Θ has been calculated to be about 21 by Smith  for powers and fiber lengths appropriate to telecom applications. However, in view of the substantial changes introduced at high powers, it is worth revisiting this calculation.
The Brillouin power PB at the back of the pump pulse can be modeled as if it is built up from thermal noise at all frequencies at the front of the pump, amplified by the frequency-dependent gain discussed above :Eq. (23) and performing the integral, we find: 1].
The threshold parameter Θ is defined as the logarithmic gain when the final Brillouin power is equal to the pump power PP; that is,
To calculate the value of Θ to be used for the experiments described below, we use the following values: υA = 6 km/sec, υ = 0.2 m/nsec, kT = 4 × 10−21 J, PP = 30 W, TB = 4 nsec, τ = 20 nsec, and γ = 50. The result is: Θ ≈ 22. Considering how different the present powers and Brillouin gain bandwidths are from those relevant to telecommunications, it is remarkable how close this result is to Smith’s result of 21.
5. Brillouin threshold in short fibers
Once the front of the pump pulse exits the fiber, the invariance condition, Eq. (7), no longer holds, and so the key assumption behind Eq. (14) is no longer valid. Therefore, to find the Brillouin threshold in fibers whose round trip time is shorter than the pulse duration, we must then return to Eqs. (4) and (5) and begin again. To solve these equations, we first change to coordinates ẑ and t̂ that are co-moving with the Brillouin radiation:Eqs. (4) and (5) become:
We define new phonon and optical Brillouin variables,Eqs. (29) and (30) become: Eq. (33), and use Eq. (34) to eliminate ∂q*/∂t̂. We end up with the equation:
We have therefore reduced the two equations describing the Brillouin and phonon waves to the single equation above. This equation can be solved using Riemann’s method, the details of which are given in the Appendix. The solution is most easily written in terms of dimensionless variables: the dimensionless gain γ introduced in Eq. (19), γ = 2gBIPυTB; a dimensionless fiber length Λ,
Using these dimensionless variables, the Brillouin threshold condition for pulses longer than the roundtrip time in the fiber can be written in terms of the modified Bessel function I0 (see Appendix):Equation (38) is the third main result of this paper, the Brillouin threshold in short fibers. Unfortunately, this equation is not as simple as the threshold in long fibers, Eq. (14). However, for a given pulse duration (related to τ̄) and fiber length (proportional to Λ), Eq. (38) can be solved numerically for γ to obtain the intensity IP at the Brillouin threshold.Eq. (1).
Equation (38) is valid whenever the pulse is longer than the roundtrip time in the fiber (τ̄ > 0, or τ > 2L/υ). The equation interpolates between the short pulse threshold given by Eq. (14) when τ = 2L/υ and the c.w. threshold when τ → ∞. Examples are shown in Fig. 1. The solid curve is the threshold for pulses shorter than the fiber round trip time, Eq. (14); the dashed curves are the thresholds for longer pulses (Eq. (38)) in fibers of length 1, 5, and 25 m. The circles and squares are data points from the experiments described below.
To test this analysis, we measured the Brillouin thresholds for different pulse durations and fiber lengths. The experimental setup is shown in Fig. 2. The c.w. output of a seed diode was modulated to create pulses; the pulses were amplified by a series of fiber amplifiers, then sent through a length of passive fiber to induce Brillouin scattering.
The seed diode was a DFB laser operating at 1060 nm, such as is available from QPC Lasers or Eagleyard Photonics. The output was formed into a train of pulses by a Mach-Zehnder modulator from Photline Technologies (model NIR-MX-LN10) driven by a pulse generator. For long pulses (τ > 20 nsec) we used a Hewlett Packard 8082A pulse generator; for short pulses, we used an Avtech AVMP-2-C-EPIA because of its shorter rise and fall times. A bias voltage was also provided to the Mach-Zehnder modulator to optimize the on/off extinction ratio at each repetition rate.
The pulses were amplified, first in a two-stage fiber pre-amplifier, then in a two-stage main fiber amplifier, the last stage of which used a cladding-pumped photonic crystal fiber available from NKT Photonics (DC-200/40-PZ-Yb). The output was put through an isolator, then coupled into polarization maintaining passive fiber (PM 980 XP from Nufern) that was angle-cleaved at the far end. A half-wave plate was used to align the polarization of the beam with one of the principal axes of the fiber. (In practice, this was done by rotating the half-wave plate until the Brillouin scattering was maximized.)
The output of the fiber was then collimated and its average power measured by a power meter. To examine the shapes of the pulses, a glass pickoff was placed in the output beam, reflecting in p-polarization, so that less than 2% of the power was deflected toward a fiber-coupled fast detector (ThorLabs SIR5) whose output was viewed on a 1 GHz oscilloscope.
To measure the Brillouin threshold, the average power of the pump beam was increased by slowly increasing the power of the main amplifier. The pulses were monitored on the oscilloscope for any sign of Brillouin scattering. When Brillouin scattering started to occur, the trailing edge of the pulse became noisy and generally reduced in intensity. The power was then further increased until the pulses reached Brillouin threshold.
First, however, an experimental criterion for the Brillouin threshold had to be defined. Mathematically, the Brillouin threshold occurs when the Brillouin intensity equals the pump intensity in the undepleted pump approximation. However, this is an idealized calculation; in reality, the pump is depleted as the Brillouin wave grows, so the Brillouin intensity never actually reaches the pump intensity, but asymptotically approaches it.
Therefore, the Brillouin threshold must be defined in practice as when the Brillouin intensity equals some fraction of the pump intensity. The exact value of this fraction is somewhat arbitrary; but fortunately, it is not too important since the threshold value only weakly depends on the fraction chosen. For example, in the c.w. case, if the Brillouin intensity at threshold is selected to be 20%, as opposed to 50%, of the pump intensity, the threshold pump power changes by only 6%.
For experimental convenience, the Brillouin threshold was defined near the visible onset of Brillouin scattering, when the back of the pump pulse was reduced by 20%. At all pulse widths, this introduced only a small loss in the pump pulse energy (estimated to be < 3%), so the launched pump power at threshold was then approximated as being equal to the power exiting the fiber. Although the threshold measurement had some uncertainty due to the noise inherent in Brillouin scattering, it was reproducible to within a few percent. An example is shown in Fig. 3.
Brillouin thresholds were measured for 50 m and 5 m lengths of fiber, with measurement errors estimated to be on the order of ±10%. The pulse repetition rate was 500 kHz or 1 MHz, which left plenty of time for the last of the Brillouin radiation to escape the fiber and for the phonons to decay back down to their thermal level before the next pulse arrived.
The results of the experiments are plotted in Fig. 1. To compare the experiment to the theory, the values of the Brillouin gain gB and the phonon lifetime TB were needed. In bulk silica, gB = 50 μm2/(W-m) and TB = 5 nsec [1,2]. However, these values are different in optical fibers: The gain gB is decreased due to the imperfect overlap of the acoustic and optical modes , so that typical values of gB range from 10 to 30 μm2/(W-m) [14–17]. The lifetime TB is also reduced in the fiber [14, 18], with typical values in the range of 2–3 nsec when the optical wavelength is 1060 nm.
Since the values of gB and TB for our fiber (Nufern PM-980 XP) were not known in advance, we fit them to the data. To obtain the gain gB, we used the two longest pulses (τ = 300 and 400 nsec) measured in the 5 m fiber, since for these pulses, the steady state threshold equation (gBIPL = Θ) holds to an excellent approximation. We used the average threshold intensity for these two pulse lengths, and calculated the Brillouin gain to be: gB = 31 μm2/(W-m), close to the values reported above.
Equation (14) was then fit to the short pulse data. This included all of the pulses measured in the 50 m fiber, and the shorter pulses (τ ≤ 50 nsec) from the 5 m fiber. For this fit, gB was fixed at the above value of 31 μm2/(W-m) and TB was determined by the least squares method. The result was: TB = 4.1 nsec, a little longer than the 2 to 3 nsec reported in other fibers, but less than the bulk silica value of 5 nsec.
Finally, the fitted values of gB and TB were used to plot Eq. (14) and Eq. (38) with L = 5 m on the same graph as the data, as shown in Fig. 1. The theory and experiment closely agree. For reference, the thresholds for 1 m and 25 m length fibers are also shown.
In conclusion, the equation for the Brillouin threshold for short pulses in long fibers, Eq. (14), has proved accurate as well as simple. The more complicated expression for the Brillouin threshold for long pulses in short fibers, Eq. (38), also agrees with experiment. Furthermore, the surprising finding that the Brillouin gain bandwidth increases dramatically for high powered pulses, Eq. (22), will be important when considering methods to prevent Brillouin scattering. All of these results will be useful when designing pulsed fiber systems.
Appendix: Riemann’s method
This appendix describes the solution of Eq. (35), subject to the appropriate boundary conditions. To find the boundary conditions, we note that for a fiber of length L, the Brillouin wave is described by Eq. (12) until the front of the pulse reaches the exit end of the fiber at z = L and t = L/υ, or, in terms of the hatted coordinates, ẑ = L and t̂ = 2L/υ. As mentioned before, Eq. (12) continues to hold at later times t for those parts of the pulse that have not yet received the news that the front of the pulse has reached the end of the fiber; that is, for all ẑ provided that t̂ ≤ 2L/υ. We therefore have:
The second boundary condition is that the Brillouin power at the far end of the fiber is due to thermal noise. That is, AB(ẑ = L, t̂) = AB0, or
We now introduce one more change of variables: we define new distance and time coordinates z̄ and t̄ that are scaled so they are dimensionless, and are chosen so that the boundary conditions occur more conveniently at z̄ = 0 and t̄ = 0. Since the first boundary condition, Eq. (39), occurs at t̂ = 2L/υ, we setEq. (40)) states that the Brillouin wave starts with thermal noise at ẑ = L; the Brillouin amplitude then grows in the negative ẑ direction to reach its maximum at z = 0. We invert the z-coordinate, and define
We also define a dimensionless Brillouin amplitudeEq. (19): γ = 2gBIPvTB. In terms of these new variables, Eq. (35) becomes Eqs. (39) and (40), can then be written: Eq. (11), and
The problem defined by the three equations above can be solved by Riemann’s method . The method consists of first finding a so-called Riemann function w(z̄, t̄) that has the following properties: first, it satisfies Eq. (44),
Assuming that such a function can be found, then according to Riemann’s method, the solution for B(z̄, t̄) is:Eq. (44) and the necessary boundary conditions, due to the properties of the Riemann function given in Eqs. (47) and (48).
The only remaining step is to find the Riemann function w(z̄, t̄). This can be done by expanding the function in a power series:Eq. (47) gives a condition on the functions cn(t̄): Eq. (48), and reasoning by induction, we find that 20], so the Riemann function can be expressed as: Eq. (43)) and 𝒜 (Eq. (32)), the threshold condition becomes: Eq. (52) into this equation yields the threshold condition, Eq. (38).
We would like to thank Jim Morehead for helpful conversations.
References and links
1. G. P. Agrawal, Nonlinear Fiber Optics, 5thedition (Academic,Oxford, 2013).
2. R. W. Boyd, Nonlinear Optics, 2nd edition (Academic, San Diego, 2003).
3. M. J. Leonardo, M. W. Byer, G. L. Keaton, D. J. Richard, F. J. Adams, K. Monro, J. L. Nightingale, S. Guzsella, and L. Smoliar, “Versatile, nanosecond laser source for precision material processing,” presented at the 28th International Congress on Applications of Lasers and Electro-Optics (ICALEO), Orlando, Florida, 2–5 Nov.2009, paper #M103.
4. N. M. Kroll, “Excitation of hypersonic vibrations by means of photoelastic coupling of high-intensity light waves to elastic waves,” J. Appl. Phys. 36, 34–43 (1965). [CrossRef]
5. D. Pohl and W. Kaiser, “Time-resolved investigations of stimulated Brillouin scattering in transparent and absorbing media: Determination of phonon lifetimes,” Phys. Rev. B. 1, 31–43 (1970). [CrossRef]
6. H. Li and K. Ogusu, “Dynamic behavior of stimulated Brillouin scattering in a single-mode optical fiber,” Jpn. J. Appl. Phys. 38, 6309–6315 (1999). [CrossRef]
7. V. P. Kalosha, E. A. Ponomarev, L. Chen, and X. Bao, “How to obtain high spectral resolution of SBS-based distributed sensing by using nanosecond pulses,” Opt. Express 14, 2071–2078 (2006). [CrossRef] [PubMed]
8. M. S. Bigelow, S. G. Lukishova, R. W. Boyd, and M. D. Skeldon, “Transient stimulated Brillouin scattering dynamics in polarization maintaining optical fiber,” CLEO2001, paper CTuZ3.
12. D. Williams, X. Bao, and L. Chen, “Characterization of high nonlinearity in Brillouin amplification in optical fibers with applications in fiber sensing and photonic logic,” Photon. Res. 2, 1–9 (2014). [CrossRef]
13. V. Lecoeuche, D. J. Webb, C. N. Pannell, and D. A. Jackson, “Transient response in high-resolution Brillouin-based distributed sensing using probe pulses shorter than the acoustic relaxation time,” Opt. Lett. 25, 156–158 (2000). [CrossRef]
14. M. Niklès, L. Thévenaz, and P. A. Robert, “Brillouin gain spectrum characterization in single-mode optical fibers,” J. Lightwave Technol. 15, 1842–1851 (1997). [CrossRef]
15. V. I. Kovalev and R. G. Harrison, “Means for easy and accurate measurement of the stimulated Brillouin scattering gain coefficient in optical fiber,” Opt. Lett. 33, 2434–2436 (2008). [CrossRef] [PubMed]
16. M. D. Mermelstein, “SBS threshold measurements and acoustic beam propagation modeling in guiding and anti-guiding single mode optical fibers,” Opt. Express 17, 16225–16237 (2009). [CrossRef] [PubMed]
17. V. Lanticq, S. Jiang, R. Gabet, Y. Jaouën, F. Taillade, G. Moreau, and G. P. Agrawal, “Self-referenced and single-ended method to measure Brillouin gain in monomode optical fibers,” Opt. Lett. 34, 1018–1020 (2009). [CrossRef] [PubMed]
18. V. I. Kovalev and R. G. Harrison, “Waveguide-induced inhomogeneous spectral broadening of stimulated Brillouin scattering in optical fiber,” Opt. Lett. 27, 2022–2024 (2002). [CrossRef]
19. S. L. Sobolev, Partial Differential Equations of Mathematical Physics (Dover, New York, 1989).
20. G. Arfken, Mathematical Methods for Physicists, 3rd edition (Academic, Orlando, 1985).