# Quasioptical modeling of wave beams with and without mode conversion:

II. Numerical simulations of single-mode beams

###### Abstract

This work continues a series of papers where we propose an algorithm for quasioptical modeling of electromagnetic beams with and without mode conversion. The general theory was reported in the first paper of this series, where a parabolic partial differential equation was derived for the field envelope that may contain one or multiple modes with close group velocities. Here, we present a corresponding code PARADE (PAraxial RAy DEscription) and its test applications to single-mode beams in vacuum and also in inhomogeneous magnetized plasma. The numerical results are compared, respectively, with analytic formulas from Gaussian-beam optics and also with cold-plasma ray tracing. Quasioptical simulations of mode-converting beams are reported in the next, third paper of this series.

## I Introduction

Electron-cyclotron heating and current drive in fusion plasmas involve electron-cyclotron waves (ECWs) and require deposition of the wave power with high precision; hence, accurate modeling of ECWs is required. Full-wave simulations of ECWs foot:vdovin06 are computationally-demanding and largely impractical due to the fact that the ECW wavelengths (which are in the mm range) are extremely small compared to the plasma size. Geometrical-optics (GO) ray-tracing simulations ref:tsujimura15 ; ref:marushchenko14 are much faster but cannot adequately describe the structure of the wave beam in the focal region. Many alternatives have been proposed to improve reduced modeling of ECWs, including complex-eikonal methods ref:mazzucato89 ; ref:nowak93 ; ref:peeters96 ; ref:farina07 , beam tracing ref:pereverzev98 ; ref:poli01b ; ref:poli01 ; ref:poli18 , and the quasioptical approach ref:balakin08b ; ref:balakin07a ; ref:balakin07b . However, none of them account for the possible interaction of the two cold-plasma modes (linear mode conversion), which can occur in strongly sheared magnetic field at low density, for example, at the plasma edge my:xo ; ref:kubo15 . Hence, a need remains for a more systematic approach to modeling ECWs.

This work continues a series of papers where we propose a new algorithm for quasioptical modeling of electromagnetic (EM) beams propagating in inhomogeneous plasma with and without mode conversion. The general theory was reported in Paper I of our series foot:paper1 . There, a parabolic partial differential equation was derived for a certain projection of the field envelope, whose dimension equals the number of the resonantly-coupled modes. Based on that theory, which is a generalization of extended geometrical optics (XGO) phd:ruiz17 ; my:covar ; my:qdirac ; my:qdiel , we have also developed a code PARADE (PAraxial RAy DEscription) and applied it to several test problems. Paper III of our series foot:paper3 reports PARADE’s applications to quasioptical simulations of mode-converting beams with . Here, we present a simplified version of the same code for , in which case mode conversion is not considered. This “scalar” version of PARADE serves two goals: (i) it is a stepping stone towards a more general code reported in Paper III; and (ii) it is useful as an independent tool for quasioptical beam tracing without mode conversion in sufficiently dense plasma, especially when the absorption is strongly inhomogeneous and no simple ansatz can be assumed for the beam structure. In this second capacity, PARADE can also be considered as an alternative to the algorithm proposed in Refs. ref:balakin07a ; ref:balakin07b ; ref:balakin08b . The latter is similar in spirit foot:balakincomp but is limited to single-mode beams by design.

This paper is organized as follows. In Sec. II, we introduce the key equations derived in Paper I and also adjust them to numerical modeling. In Sec. III, we outline the numerical algorithm used in PARADE. In Sec. IV, we report simulation results for test problems. In Sec. V, we summarize our main conclusions.

## Ii Theoretical model

### ii.1 Basic equations

Here, we outline how the general theory developed in Paper I can be applied, with some adjustments, to describing single-mode ECWs. We start by assuming a general linear equation for the electric field of a wave,

(1) |

where is a linear dispersion operator discussed below. We also assume that the field can be represented in the eikonal form,

(2) |

where is a slow complex vector envelope and is a fast real “reference phase” to be prescribed (Sec. II.1). The wave is considered stationary, so it has a constant frequency ; then and are functions of the spatial coordinate . Correspondingly, the envelope satisfies

(3) |

where serves as the “envelope dispersion operator” ( denotes definitions). We introduce for the local wave vector, for the wavelength, for the characteristic scale of the beam field along the group velocity at the beam center, and for the minimum scale of the field in the plane transverse to this group velocity. The medium-inhomogeneity scale is assumed to be larger than or comparable with .

We adopt the quasioptical ordering, specifically,

(4) |

Using this ordering, the “envelope dispersion operator” can be expressed as

(5) |

The operator is specified in Paper I (also see below). The matrix serves as the “effective dispersion tensor” found from , and satisfies the ordering

(6) |

where the indices and denote the Hermitian and anti-Hermitian parts, correspondingly. Assuming that spatial dispersion is weak enough, can be replaced foot:paper1 with the homogeneous-plasma dispersion tensor,

(7) |

where is a unit matrix and is the homogeneous-plasma dielectric tensor book:stix ; its dependence on is assumed but not emphasized, since is constant. (Here, denotes any given wave vector, as opposed to , which is the specific wave vector determined by ; see above.) Note that Eq. (7) assumes the Euclidean metric. Although we shall also use curvilinear coordinates below, expressing in those coordinates will not be necessary as explained in Sec. VI D of Paper I.

With and given by Eq. (7), Eq. (6) implies . Hence, is the dominant part of , and Eqs. (5) yields

(8) |

As a Hermitian matrix, has enough (three) eigenvectors to form a complete orthonormal basis. Then, it is convenient to represent the envelope in this basis,

(9) |

where are the corresponding complex amplitudes. [Summation over repeating indices is assumed. For all functions derived from , such as here, the notation convention will also be assumed by default.] Then,

(10) |

where are the corresponding eigenvalues. Due to Eq. (8) and the mutual orthogonality of all , this means that for every given individually; hence, either is small or is small. In the latter case, approximately satisfies the eigenmode equation, , so it can be viewed as the local polarization vector of a GO mode. Then, the corresponding can be understood as the local scalar amplitude of an actual GO mode, so that is allowed. In this paper, we assume conditions such that only one mode can be excited at given near the fixed frequency of interest. Then,

(11) |

where the polarization vector is the eigenvector corresponding to the excited mode and is the corresponding eigenvalue. (The notation is chosen here such that our formulas extend to mode-converting waves with generalized definitions of , , and as in Refs. foot:paper1 ; foot:paper3 .) The term in Eq. (11) can be easily calculated from and is included in our theory as a perturbation foot:paper1 but will not be considered below explicitly. Hence, we shall describe the wave in terms of , which can also be expressed as follows:

(12) |

Here, the row vector is the dual polarization vector, so . Since we consider the beam dynamics in coordinates that are close to Euclidean, it will be sufficient, within the accuracy of our model, to consider the row vector simply as the conjugate transpose of the column vector .

### ii.2 Reference ray

We shall consider the beam evolution in the coordinate system that is determined by a family of rays. These rays are governed by Hamilton’s equations foot:paper1

(13) |

where serves as an effective time variable, denote the ray location in the laboratory coordinates, and denote the ray wave-vector components in the associated momentum space. The initial conditions are assumed such that is of order or smaller foot:paper1 . (A more specific definition will not be needed.) Then, it also remains small indefinitely, because is conserved by Eqs. (13).

Let us introduce the path along each ray, , where and is given by

(14) |

At any given , the coordinates of rays with different initial form a “transverse space”. Let us define a two-dimensional transverse coordinate on this space such that its origin corresponds to the location of some “reference ray” (RR). This ray can be chosen arbitrarily as long as it is representative of the rays within the beam. The coordinate and the wave vector of the RR are denoted as and , respectively. Like the coordinate and momentum of other rays [Eqs. (13)], they can be found by solving the ray equations

(15) |

where the index denotes evaluation at . Note that we require to be exactly zero initially; then, it also remains zero at all ,

(16) |

### ii.3 Ray-based coordinates

Since there is a one-to-one correspondence between the laboratory coordinates and , one can choose the latter as the new coordinates. The two sets of coordinates are connected via

(17) |

where are the basis vectors of the “tilded” coordinate system. Accordingly, the transformation matrices and

(18) |

where , are given by

(19) |

We choose the transverse coordinates such that

(20) |

(Here and further, the indices and span from 1 to 2; other Greek indices span from 1 to 3. The specific choice of used in PARADE is discussed in Sec. III.) Then,

(21) |

and

(22) |

where the prime denotes the derivative with respect to . This leads to

(23) |

where is a matrix given by

(24) |

Using that , we also obtain ,

(25) |

For future references, note that we shall use tilde also to mark the components of vectors and tensors measured in these new coordinates. As usual, these components are connected to those in the laboratory coordinates via book:landau2

We also introduce first-order partial derivatives

(26) |

and the second-order derivatives are denoted as follows:

Accordingly, , and so on.

### ii.4 Quasioptical equation

#### ii.4.1 Basic formulas

Using the fact that the metric in the assumed ray-based coordinates is Euclidean on the RR and overall smooth, the envelope equation can be represented as follows foot:paper1 :

(27) |

Here, we introduced

(28) | |||

(29) | |||

(30) |

and also the following notation for the derivatives:

(31) |

Note that all coefficients in Eq. (27) are considered as functions of , namely, for any . Accordingly, the derivatives (31) act on both arguments. Also note that and are calculated in the laboratory coordinates; otherwise, additional terms would emerge foot:paper1 .

Within the assumed accuracy, and can be replaced with those evaluated on the RR,

(32) | |||

(33) |

In contrast, , , and generally must be calculated with a higher accuracy. (An example of a numerical simulation that ignores this fact is presented in Sec. IV.3 for a reference.) In order to do that, we need to find , so we start by specifying , through which is defined (Sec. II.1). Let us adopt

(34) |

(This choice of is different from the one assumed in Paper I but similar to that in LABEL:ref:balakin07a.) Hence,

(35) | |||

(36) |

To simplify this, note that

(37) |

by definition. Hence, the transformation is symplectic (canonical), so the “tilded” coordinates satisfy Hamilton’s equations similar to Eq. (15),

(38) |

By combining Eqs. (35), (36), and (38), one arrives at

(39) |

where we used , and is a Kronecker symbol.

#### ii.4.2 Calculating and

Expanding around the RR, one obtains

(40) |

Using Eqs. (16) and (39), one can reduce this down to

which is of the order of , as needed. Since the metric on the RR is Euclidean, the second-order derivatives entering this equation are transformed as true tensors. Hence, they can be mapped to the laboratory coordinates with and . This leads to with

(41) |

where we substituted [cf. Eq. (39)]

(42) |

(Note that, to the leading order, can be considered as a true vector and transforms accordingly.)

Using Eqs. (42), one can also adopt the following approximation of Eq. (28) for :

(43) |

Here, the polarization vectors and are evaluated at the RR (hence the stars), because they are determined by , which is a slow function. In contrast, can vary substantially in the beam cross-section, so it can be important to calculate separately at each point.

#### ii.4.3 Calculating

In order to calculate the terms involving , let us represent this function as ,

(44) |

where denotes deviations from the corresponding values of the RR. Using Eq. (25), we obtain

(45) |

Also note that

(46) |

Hence, , where

(47) |

Then, since and , one obtains

(48) |

where can be calculated from Eq. (47) by taking and summing over accordingly. Also,

(49) |

#### ii.4.4 Final result

Using the above results, one can express Eq. (27) as

(50) |

which also has the following corollary:

(51) |

Equation (51) shows that is conserved when is zero. As explained in Paper I, this represents the conservation of the wave action, and equals the energy flux up to a constant coefficient. Assuming the notation

(52) |

so , one can also rewrite Eq. (50) as

(53) |

Equation (53) is the main quasioptical equation used in PARADE. As a reminder, we assume Eq. (49) for , Eq. (43) for , Eq. (41) for , Eq. (33) for , and Eq. (32) for . We also assume Eq. (47) for , with Eq. (24) for . The quantity can be calculated from Eq. (47) by taking and summing over accordingly. The matrix is given by Eq. (22), and the matrix is given by Eq. (18). The index denotes that the corresponding quantities are evaluated at the RR (i.e., at ) and thus are independent of the transverse coordinates.

## Iii Numerical algorithm

The code PARADE is based on the model described in Sec. II. The RR starting point , the orientation of the RR initial wave vector, the frequency , the initial profile of the complex amplitude , the magnetic field profile , and the plasma density profile are considered pre-determined. Currently, they are entered as analytic functions, but it is also possible to use experimental data or profiles generated numerically. The value of is determined by the local dispersion relation at . The collisionless cold-electron-plasma model is assumed for the dielectric tensor book:stix here, but the code can also accept a general dispersion tensor . In particular, warm-plasma simulations will be reported separately.

Figure 1 outlines the code operation. At each point on the RR, starting from the initial point , the dispersion tensor is calculated. Then, the eigenvalues of are calculated, and the eigenvalue is found that has the least absolute value. This is then numerically differentiated and used as the Hamiltonian in Eqs. (15) to propagate the RR using the fourth-order Runge–Kutta method. Based on the knowledge of the RR trajectory , the ray-based coordinates are constructed by choosing the basis vectors , which is done as follows.

We adopt to be the unit vector along the RR group velocity. The other two vectors (where ) can, in principle, be chosen arbitrarily as long as they satisfy Eqs. (20) and the resulting metric is smooth. The traditional Frenet–Serret frame can be used but sometimes fails the latter requirement. Instead, we choose

(54) | |||

(55) | |||

(56) |

at each given time step, where the vectors marked with are those from the previous time step. (The prime denotes the derivative with respect to , so .) For the very first time step, and can be chosen arbitrarily as long as they are orthogonal to each other and to the initial group velocity on the RR. This ensures that the basis vectors transform continuously along the RR and result in a smooth ray-based metric.

As the next step, the polarization eigenvector is calculated as the eigenvector of corresponding to the eigenvalue ; also, its complex conjugate is determined. After these vectors are numerically differentiated, the coefficients in Eq. (53) are calculated as outlined in Sec. II.4.4. The field is mapped to using Eq. (52), and is differentiated with respect to at fixed using the finite-difference approach. Then, is propagated by solving Eq. (53) using the fourth-order Runge–Kutta method. Finally, is mapped back to , and the three-dimensional profile is yielded as a data array.

## Iv Simulation results

In this section, test simulation results are reported. These results were obtained by applying PARADE to sample problems that fit the assumptions of the model described above. Simulations that involve mode conversion are presented in Paper III.

### iv.1 Beams in vacuum

As the first example, we simulated the propagation of a wave beam in vacuum. For the initial beam profile, we assumed a Gaussian profile

(57) |

where

(58) | |||

(59) |

and , with the focal lengths m and m, the waist sizes cm and cm, and the wave frequency GHz, which corresponds to the vacuum wavelength mm. The results are shown in Fig. 2. At , the beam width is determined numerically as the distance on which drops one -fold from its maximum along the -axis. As seen from Fig. 2(d), PARADE demonstrates good agreement with the analytical formula from Gaussian-beam optics book:yariv .

### iv.2 Straight beams in cold plasma

As the second example, we simulated the propagation of waves in collisionless magnetized cold electron plasma. From now on, we introduce the standard notation for the laboratory coordinates and for the unit vectors along the corresponding axes. In particular, the axis is chosen to be along the initial direction of the wave vector. We assume a slab geometry with

(60) |

Here, m, m, and T. The initial field has the Gaussian profile (57) with the focal lengths m, the waist sizes cm, and the wave frequency GHz. The simulation results are presented in Fig. 3. Since the waves propagate along the gradient of , the RR trajectory is straight and does not depend on the polarization, just like in the vacuum case. However, the intensity profiles of the waves with different polarizations, which are the O and X waves here, are different from each other and from those in vacuum. The X wave is affected by the plasma more significantly. This is explained by the fact that the wave frequency is closer to the X-wave cutoff frequency than to the O-wave cutoff frequency . In addition, although and are assumed to be equal initially [dashed line in Fig. 3(b) and (c)], they start to differ within the plasma [solid lines in Fig. 3(b) and (c)] due to the magnetic field. Also, while the focal lengths are different from those in vacuum, the waist sizes are, in fact, the same. This is because is uniform in the transverse plane for the lack of the density gradient in that plane. Simulations with transverse gradients are presented in Sec. IV.3. Then, we have also tested PARADE for the energy conservation [Eq. (51)], which is expected here because waves in collisionless cold plasma have . While the wave intensity and group velocity change considerably along the RR path, the wave-energy flux is conserved by the code with high accuracy, as seen in Figs. 3(d) and (e).

### iv.3 Bended beams in cold plasma

As another example, we applied PARADE to simulate the propagation of bended wave beams launched at generic angles relative to the density gradient and the magnetic field. The assumed density and magnetic-field profiles are as follows:

(61) | |||

(62) |

where m, m, m, T, and m. The RR starting point is m, and the target point used to fix the orientation of the RR initial wave vector is m. For the initial profile, we adopted the Gaussian profile (57) with the focal lengths m, the waist sizes cm, and the wave frequency GHz (or mm). A simulation of the X-mode propagation is shown in Fig. 4(a). Since PARADE accounts for diffraction, the beam retains a nonzero width on the whole trajectory. This is different from conventional ray tracing, where diffraction is neglected and each ray propagates independently. A comparison is presented in Fig. 4(b). The difference is essential for resonant-plasma heating and current drive, where the ECW energy deposition must be calculated accurately, as also discussed in Refs. ref:mazzucato89 ; ref:nowak93 ; ref:peeters96 ; ref:farina07 ; ref:pereverzev98 ; ref:poli01b ; ref:poli01 ; ref:poli18 ; ref:balakin08b ; ref:balakin07a ; ref:balakin07b .

The corresponding evolution of the key frequencies along the RR is presented in Fig. 5(a), and the corresponding beam widths are presented in Fig. 5(b). For a reference, Fig. 5(c) shows the same widths calculated with neglected -dependence of . (In the latter case, the waist sizes are almost the same as for vacuum waves.) These results show that capturing the transverse inhomogeneity of in practical simulations can be important.

We also note that in our cold-plasma PARADE simulations, the beams remain Gaussian with high accuracy (Fig. 6). However, it is generally not the case when resonant dissipation is included, for such dissipation can be strongly inhomogeneous in the beam cross-section. PARADE simulations of dissipative beams in hot plasma will be reported in a separate publication.

## V Conclusions

We report the first quasioptical simulations of ECWs in cold magnetized plasma using a new code PARADE. This code does not assume any particular beam structure but instead solves a parabolic equation for the wave envelope; hence, it is particularly well-suited to modeling strongly inhomogeneous resonant dissipation in the future. The general theoretical model underlying PARADE was largely derived in the previous paper foot:paper1 and describes multi-mode beams that can experience mode conversion. Here, we consider a simplified version of this theory in which mode conversion is not considered and a beam can be described as a scalar field. We adjust the theory to numerical modeling and simulate the propagation of wave beams in vacuum and inhomogeneous magnetized plasma. Our numerical results show good agreement with Gaussian-beam optics and surpass conventional ray tracing. Quasioptical simulations of mode-converting beams are reported in the next, third paper of this series.

The work was supported by the U.S. DOE through Contract No. DE-AC02-09CH11466.

## References

- (1) V. Vdovin, Role of upper hybrid resonance and diffraction effects at electron cyclotron heating in tokamaks, Proc. 14th Joint Workshop on Electron Cyclotron Emission and Electron Cyclotron Heating (Santorini, Greece, May 9-12, 2006).
- (2) T. I. Tsujimura, S. Kubo, H. Takahashi, R. Makino, R. Seki, Y. Yoshimura, H. Igami, T. Shimozuma, K. Ida, C. Suzuki, M. Emoto, M. Yokoyama, T. Kobayashi, C. Moon, K. Nagaoka, M. Osakabe, S. Kobayashi, S. Ito, Y. Mizuno, K. Okada, A. Ejiri, T. Mutoh, and the LHD Experiment Group, Development and application of a ray-tracing code integrating with 3D equilibrium mapping in LHD ECH experiments, Nucl. Fusion 55, 123019 (2015).
- (3) N. B. Marushchenko, Y. Turkin, and H. Maassberg, Ray-tracing code TRAVIS for ECR heating, EC current drive and ECE diagnostic, Comput. Phys. Commun. 185, 165 (2014).
- (4) E. Mazzucato, Propagation of a Gaussian beam in a nonhomogeneous plasma, Phys. Fluids B 1, 1855 (1989).
- (5) S. Nowak and A. Orefice, Quasioptical treatment of electromagnetic Gaussian beams in inhomogeneous and anisotropic plasmas, Phys. Fluids B 5, 1945 (1993).
- (6) A. G. Peeters, Extension of the ray equations of geometric optics to include diffraction effects, Phys. Plasmas 3, 4386 (1996).
- (7) D. Farina, A quasi-optical beam-tracing code for electron cyclotron absorption and current drive: GRAY, Fusion Sci. Tech. 52, 154 (2007).
- (8) G. V. Pereverzev, Beam tracing in inhomogeneous anisotropic plasmas, Phys. Plasmas 5, 3529 (1998).
- (9) E. Poli, G. V. Pereverzev, A. G. Peeters, and M. Bornatici, EC beam tracing in fusion plasmas, Fusion Eng. Des. 53, 9 (2001).
- (10) E. Poli, A. G. Peeters, and G. V. Pereverzev, TORBEAM, a beam tracing code for electron-cyclotron waves in tokamak plasmas, Comput. Phys. Commun. 136, 90 (2001).
- (11) E. Poli, A. Bock, M. Lochbrunner, O. Maj, M. Reich, A. Snicker, A. Stegmeir, F. Volpe, N. Bertelli, R. Bilato, G. D. Conway, D. Farina, F. Felici, L. Figini, R. Fischer, C. Galperti, T. Happel, Y. R. Lin-Liu, N. B. Marushchenko, U. Mszanowski, F. M. Poli, J. Stober, E. Westerhof, R. Zille, A. G. Peeters, and G. V. Pereverzev, TORBEAM 2.0, a paraxial beam tracing code for electron-cyclotron beams in fusion plasmas for extended physics applications, Comput. Phys. Commun. 225, 36 (2018).
- (12) A. A. Balakin, M. A. Balakina, and E. Westerhof, ECRH power deposition from a quasi-optical point of view, Nucl. Fusion 48, 065003 (2008).
- (13) A. A. Balakin, M. A. Balakina, G. V. Permitin, and A. I. Smirnov, Quasi-optical description of wave beams in smoothly inhomogeneous anisotropic media, J. Phys. D: Appl. Phys. 40, 4285 (2007).
- (14) A. A. Balakin, M. A. Balakina, G. V. Permitin, and A. I. Smirnov, Scalar equation for wave beams in a magnetized plasma, Plasma Phys. Rep. 33, 302 (2007).
- (15) I. Y. Dodin, D. E. Ruiz, and S. Kubo, Mode conversion in cold low-density plasma with a sheared magnetic field, Phys. Plasmas 24, 122116 (2017).
- (16) S. Kubo, H. Igami, T. I. Tsujimura, T. Shimozuma, H. Takahashi, Y. Yoshimura, M. Nishiura, R. Makino, and T. Mutoh, Plasma interface of the EC waves to the LHD peripheral region, AIP Conf. Proc. 1689, 090006 (2015).
- (17) I. Y. Dodin, D. E. Ruiz, K. Yanagihara, Y. Zhou, and S. Kubo, Quasioptical modeling of wave beams with and without mode conversion: I. Basic theory, arXiv:1901.00268v2.
- (18) D. E. Ruiz, Geometric theory of waves and its applications to plasma physics, Ph.D. Thesis, Princeton Univ. (2017), arXiv:1708.05423.
- (19) D. E. Ruiz and I. Y. Dodin, Extending geometrical optics: A Lagrangian theory for vector waves, Phys. Plasmas 24, 055704 (2017).
- (20) D. E. Ruiz and I. Y. Dodin, Lagrangian geometrical optics of nonadiabatic vector waves and spin particles, Phys. Lett. A 379, 2337 (2015).
- (21) D. E. Ruiz and I. Y. Dodin, First-principles variational formulation of polarization effects in geometrical optics, Phys. Rev. A 92, 043805 (2015).
- (22) K. Yanagihara, I. Y. Dodin, and S. Kubo, Quasioptical modeling of wave beams with and without mode conversion: III. Numerical simulations of mode-converting beams, submitted in parallel.
- (23) A comprehensive comparison of our model with those presented in Refs. ref:balakin07a ; ref:balakin07b ; ref:balakin08b would be cumbersome and is not attempted here also for the lack of detailed derivations in the mentioned papers. The verification of our model relies instead on the numerical results presented here and in Paper III, and also on the XGO applications to several test problems reported earlier phd:ruiz17 ; my:covar ; my:xo ; my:qdirac ; my:qdiel .
- (24) T. H. Stix, Waves in Plasmas (AIP, New York, 1992).
- (25) L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields (Pergamon Press, New York, 1971).
- (26) A. Yariv, Quantum electronics (Wiley, New Jersey, 1967).