872 J. Opt. Soc. Am. A/Vol. 14, No. 4 /April 1997 M. L. Calvo and V. Lakshminarayanan Light propagation in optical waveguides: a dynamic programming approach Maria Luisa Calvo Departamento de Optica, Facultad de Ciencias Fisicas, Universidad Complutense, 28040 Madrid, Spain Vasudevan Lakshminarayanan School of Optometry and Department of Physics and Astronomy, University of Missouri–St. Louis, 8001 Natural Bridge Road, St. Louis, Missouri 63121-4499 Received February 22, 1996; revised manuscript received October 23, 1996; accepted November 19, 1996 We apply techniques of optimal-control theory, namely, the methods of dynamic programming, to the problem of light propagation in optical waveguides. This formulation is equivalent to the resolution of an eikonal equation. We illustrate this optimization technique for the case of an ideal parabolic refractive-index-profile distribution. We discuss the possibility of extending this procedure to other types of optical waveguides and optical media. © 1997 Optical Society of America [S0740-3232(97)00105-1] 1. INTRODUCTION The use of optimization techniques to describe light propagation phenomena through optical media with dif- ferent optical properties and geometry is at present a relatively new area of investigation. This methodology concerns optimized design of arbitrary optical media for which an a priori characterization of significant param- eters would be required. The problem of characterizing the output of a system from partial input data is defined as a control theoretical problem. To the best of our knowledge, this approach has not been used to study light propagation. (Other relevant applications include opti- cal system design, laser cavities, laser beam modification, minimizing of geometrical aberrations, etc.) Currently, optimization techniques applied to wave- guide design are based on variational analysis related to the Sturm–Liouville theory and the Rayleigh–Ritz method.1 Recently, for example, Sharma and Bindal ap- plied such a technique to optimize diffuse planar and channel waveguides.2 Another approach that consider- ably simplifies this analysis is the well-known ray-tracing technique.3 Ray trajectories are solutions from the eiko- nal approach to the partial differential equation for the electromagnetic field propagating inside an optical me- dium. This technique has the restrictions inherent in the integration conditions of the ray-path integral and in the required initial ray conditions. Also, an inverse pro- cedure is not always applicable, for obvious reasons. One of the main problems in optimizing the ray-trajectory so- lution is the application of optimal interpolant proce- dures. This usually requires the implementation of com- putational methods, such as the Runge–Kutta, with the cost of the time and design of complex computational techniques.4,5 The search for new simplified techniques requires a comparison between light propagation phe- nomena in inhomogeneous media and mechanical models. In this type of treatment the dynamic properties of refrac- 0740-3232/97/040872-10$10.00 © tive phenomena are given in terms of two Hamiltonian equations. The solutions give both the vector position and the tangent to the trajectory.6 It should be noted that equivalent Lagrangian formulations for studying ray paths in cylindrically symmetric media have been pro- posed and have been applied to obtain exact ray paths in bent slabs as well as in bent fibers with separable profiles.7 Ghatak and Sauter8,9 have discussed the harmonic- oscillator problem in the domains of classical and quan- tum mechanics and the parabolic-index optical-waveguide problem in the domains of ray and wave optics, using a Lagrangian formulation. They showed that the time evo- lution of a coherent state has a close relationship to the propagation of a Gaussian beam in a parabolic-index waveguide. With the aim of finding simple alternative models, we have applied the method of dynamic program- ming to light propagation in inhomogeneous media. There is a former contribution to this idea in the earlier paper of Kalaba,10 where it was demonstrated that the ei- konal equation may be derived directly from Fermat’s principle of least time by using Bellman’s principle of optimality.11,12 In the present paper we analyze the pro- cedure and extend it to inhomogeneous media with arbi- trary refractive-index distribution. The paper is organized as follows: Section 2 presents background with the fundamentals of the adaptive- control process. In Section 3 we revisit the dynamic- programming approach to optical systems and demon- strate that the eikonal solution is a Hamilton–Jacobi-type equation. Specific solutions derive from the Euler– Lagrange formulation, which is obtained from initial con- ditions imposed on the system. Section 4 is dedicated to application of these results to homogeneous optical media and Section 5 to general cases of gradient-index optical waveguides. In Section 6 we analyze the case of a para- bolic profile and present some numerical results. We 1997 Optical Society of America M. L. Calvo and V. Lakshminarayanan Vol. 14, No. 4 /April 1997 /J. Opt. Soc. Am. A 873 compare the results for the ray transit time with those predicted by the geometrical optics approximation. Fi- nally, Section 7 consists of discussion and conclusions re- sulting from the various methods. 2. HISTORICAL BACKGROUND: FUNDAMENTALS OF THE MATHEMATICAL THEORY OF THE ADAPTIVE-CONTROL PROCESS The mathematical theory of the adaptive-control process was first introduced by Bellman in 1957.11 Later, Bell- man and Kalaba13 and Dreyfus,14 as well as Bellman and Vasudevan15 and others,16 extended the applications to a wide variety of problems in physics, chemistry, and biol- ogy as well as in engineering, economics, and manage- ment science. In brief, control theory determines the be- havior of arbitrary systems through optimized output responses, Fout . Consider an arbitrary system repre- sented by the action of an operator S (linear or nonlin- ear); then Fout 5 S@F in#, where F in is any arbitrary input. More-sophisticated representations are also feasible by introduction of a higher number of inputs interacting with diverse categories of systems and networks (a mul- tistage control process). In general, the mathematical framework used to solve a control process deals with the minimization of a qua- dratic functional17 of the form J~u, v ! 5 ET dt@u2~t ! 1 v2~t !#, (1) where u(t) and v(t) have a ligature and T is an arbitrary time interval. In this case u(t) and v(t) might, for ex- ample, be related by a linear differential equation of the form du dt 5 au 1 v, (2) with the initial condition u(0) 5 C (C arbitrary constant). This class of quadratic nonlinear functional is charac- teristic of equations derived by invariant imbedding techniques.17 Let us assume a system S having a certain state at any time t. The description of this state is a first global constraint, u 5 u(t), whose solution is determined by an initial condition. u(t) is a scalar function whose physical meaning is a displacement. Its first derivative, u8 5 du/dt, has dimensions of velocity. In this way, one can determine other global constraints such as the accel- eration or maybe other higher-order derivatives, by for- mulating a nonlinear differential equation (at least of the second degree): u9 1 u8 1 u 5 g~u, u8!. (3) In Eq. (3) the initial conditions are u(0) 5 c1 and u8(0) 5 c2; c1 ,c2 constants. One may search for an explicit solution for g(u, u8) by minimizing the quadratic functional: J~u, u8! 5 ET ~u2 1 u82!dt. (4) In the simplest case, the description of the system is with one single-state variable, the displacement u(t), and the specific problem will be the minimization of the functional together with the initial condition u(0) 5 c1. The solu- tion is not always a trivial one. The most common and well-known methods are those based on linear partial dif- ferential equations and the calculus of variations.18,19 For example, a formal method would be through the Eu- ler equation, and the uniqueness of the solution as a de- sired minimum has to be proven. A different approach to solving the functional is based on the theory of dynamic programming, defined as an extension of the feedback- control process.10 The formalism of dynamic program- ming requires the definition of a minimum pathway. Di- rect application of Fermat’s principle leads to the classical problem of solving an eikonal equation. Figure 1 dis- plays the interrelationships schematically. We now introduce the idea of ‘‘bang-bang’’ control. In bang-bang control theory the initial condition assigns bounded values to the origin of the trajectory (they are discrete and can be either positive or negative). Minimi- zation of time is ensured by performing first an accelera- tion from the initial point up to an intermediate one and then a deceleration from this intermediate point up to the last position. For light propagation phenomena this pro- cedure has no specific physical meaning, and one assumes only that it ensures a defined behavior for the slope, e.g., linear concerning the trajectory. We refer the readers to the many excellent texts on dynamic programming and adaptive-control process listed in the references for fur- ther details. Fig. 1. Schematic diagram showing the connection between dy- namic programming and optimized trajectory (time minimiza- tion) for light propagating in an arbitrary optical medium. 874 J. Opt. Soc. Am. A/Vol. 14, No. 4 /April 1997 M. L. Calvo and V. Lakshminarayanan 3. DYNAMIC PROGRAMMING REVISITED: GENERAL DESCRIPTION AND FORMULATION OF THE PROBLEM IN OPTICAL SYSTEMS Let us consider Fig. 2(a). The graph z 5 z(x) represents an arbitrary path for light propagating from point P(x, y) to point P0(x0 , y0). We have to search for an optimal path minimizing the time associated with the trajectory PP0. The bang-bang control process needs a control strategy; one has to search for specific values of the slope of the trajectory dz/dx and look for the most efficient pro- cedure. In this fashion the half-trajectory PQ is recov- ered with constant time h, whereas path QP0 has associ- ated with it a variable time t and therefore an arbitrary speed. We shall formulate the problem in one dimension for simplicity. The optical medium has refractive-index distribution n 5 n(x). The speed of light associated with any point P(x, z) inside the medium is v@x~z !# 5 c n@x~z !# , (5) where c is the speed of light in vacuum. We define the minimum time required to optimize the trajectory PQP0. The first policy is to establish an initial decision on the value of the slope trajectory, tan u. Following Kalaba’s formulation,10–12 for the first policy to hold the half-trajectory, QP0 may have an associated minimum time t (optimized): t~x, z ! 5 u min$h 1 t@x 1 e~x, z !cos u, z 1 e~x, z ! 3 sin u#% 1 O~h2!. (6) In Eq. (6) the trajectory e(x, z) 5 v(x, z)h represents the variable path associated with QP0 and may be a mini- mum. O(h2) is the contribution of powers of h @ 1. The term between curly brackets is the minimum time re- quired for the trajectory PP0 to be an optimal process. Applying Taylor’s formula in Eq. (6) and after simplifi- cation, we obtain u min$h 1 e~x, z !@tx cos u 1 tz sin u#% 1 O~h2! 5 0, (7) where tx 5 ]t/]x, tz 5 ]t/]z. As a first approximation, the second derivatives ]2t/]x2, ]2t/]z2 in Eq. (7) are ne- glected if the time of the trajectory has a linear depen- dence on the ray trajectory. Dividing in Eq. (7) by h and taking the limit as h → 0, we get u min$tx cos u 1 tz sin u% 5 2 1 v~x, z ! . (8) Equation (8) represents an eikonal equation. It is re- lated to the general formulation of Maxwell’s theorem for absolute instruments.20 The solution t(x, z) associates a general path with an integral. The kernel is a function that depends on the vector position of a point (x, z) on the ray trajectory and its first derivatives, dx/ds and dz/ds, with respect to s. Here s is the distance along the ray path. This statement implies that Eq. (8) has a unique solution that is a true minimum as shown by Caratheodo- ry’s theorem.21 Differentiating in Eq. (8) with respect to u gives tan u 5 ~]t/]x ! ~]t/]z ! . (9) Equation (9) represents the slope minimizing the path- way PP0. For time t a constant in the z direction, the pathway is orthogonal to the ray trajectory coinciding with the wave front generating at P0 ; and for time t a constant in the x direction, the pathway follows the ray trajectory. The condition for the slope to be a minimum is that u also be a minimum. Rearranging Eq. (9) and according to Eq. (8), we obtain S ]t ]x D 2 1 S ]t ]z D 2 5 n2~x, z ! c2 5 1 v2~x, z ! , (10) with v(x, z) the speed of light propagating in a medium with refractive index n(x, z) (Ref. 22). Equation (10) is the key equation that relates the tran- sit time t to the index distribution n(x, z). It is a Hamilton–Jacobi equation in two (or more) variables. In general, it has no explicit known rigorous solution. The solution depends on the optical properties of the medium under study. Solutions of the form t 5 t(x, z) are known as geometrical time surfaces or geometrical time Fig. 2. (a) Introduction to the problem: choosing an arbitrary trajectory z 5 z(x) (z axis, longitudinal direction; x axis, trans- verse direction) as an arbitrary path for light propagating from an initial point P(x, z) (initial time t 5 0) to a final point P0(x0 , z0). The slope tan u is determined. (b) Contour S8 de- fines an external region with surface S2 , where there is propaga- tion of light and n2(x) . 0. Contour S defines an internal re- gion with surface S1 , where also n 2(x) . 0. ū is a unit vector normal to S. For the line contour S8, n2(x) 5 0; outside this re- gion, n2(x) , 0 and Eq. (10) loses the meaning of a real transit- time equation. M. L. Calvo and V. Lakshminarayanan Vol. 14, No. 4 /April 1997 /J. Opt. Soc. Am. A 875 trajectories. Their counterparts are the geometrical wave fronts from the eikonal equation.23 The following additional condition also holds: F S ]e ]x D 2 1 S ]e ]z D 2G 1 F e~x, z ! v~x, z !G 2F S ]v ]x D 2 1 S ]v ]z D 2G 2 2F e~x, z ! v~x, z !GF S ]e ]x D S ]v ]x D 1 S ]e ]z D S ]v ]z D G 5 1. (11) It is easy to check that Eq. (11) is analogous to Eq. (10), expressing the eikonal equation in terms of the spatial variables e and v. After some convenient operations (see Appendix A) it also reads F ] ]x S ln e e0 v0 v D G2 1 F ] ]z S ln e e0 v0 v D G2 5 1 e2 , (12) where e0 is a constant initial path having a constant ini- tial velocity v0 for certain given initial coordinates, x 5 xi , z 5 zi . Equation (12) is a nonlinear partial dif- ferential equation whose solutions are not trivial, because to be trivial it would require additional conditions for (e/v). For example, assuming (e/v) independent of one of the two variables, say, (]/]z)(ln e/v) ' 0 or (e/v) de- pending only on Ax2 1 z2, which is the case of a medium with spherical symmetry, Eq. (12) simplifies to the one- dimensional case that has analytical solutions. Notice that there is an equivalent equation in quantum mechan- ics with approximate solutions obtained by applying per- turbative methods. This is the case of a classical Hamil- tonian in two dimensions, the so-called Toda’s potential or Toda lattice,24–25 for which the equivalent of Eq. (12) has an explicit solution. The latter is expressed as a sum of exponential functions with exponents having linear de- pendence on x and z variables. At present there is no evidence, nor available data in the literature, of optical media having such a similar refractive-index distribution. Nevertheless, establishing the correspondence between quantic potential V and electromagnetic dielectric perme- ability e, interesting similarities appear between Toda’s potential and Luneburg’s lens,26 in which the source is placed at the surface of the unit sphere and the focus at infinity. This lens is well known in microwave optics ap- plications as a generator of parallel rays from a virtual point source.27 It should be noted that Luneburg’s lens is a generalization of Maxwell’s ‘‘fish-eye’’ absolute instru- ment. Summarizing the above calculations, it seems more suitable, at least in the present study, to formulate the solution of the trajectory for the time variable t(x, z) as given in Eq. (10), since application of Liouville’s theo- rem ensures integrability of the equation.28 Let us assume that the optical medium has a defined geometry and is finite (at least in one of the two dimen- sions, 2x0 < x < 1 x0) and that the refractive index has bounded values, ncl < n(x) < n0. For Eq. (10) to have a complete solution, we have to introduce the boundary con- ditions. Equation (10) has a physical meaning only if n2(x) . 0. Nevertheless, if the medium has a defined geometry and is finite, as in the case of an optical wave- guide, there can be a contour S8 limiting the region where light propagates for which n2(x) 5 0. In the classical ex- ample that we are giving, this is the condition for the presence of an evanescent wave. Outside this region, n2(x) , 0, and Eq. (10) loses the physical meaning of real transit-time solutions for geometrical ray trajectories. Figure 2(b) describes two physically interesting types of boundary S and S8. By assumption, S corresponds to the possibility that n2 has a finite discontinuity across it, so that we have two regions, S1 and S2, on both sides of S: in both S1 and S2, n 2 . 0. For points on line S8, n2(x) 5 0; outside S8, n2(x) , 0. Let us discuss the boundary conditions for contour S. At the surface of discontinuity one has to consider inci- dent and reflected waves, to be matched by the refracted wave (the law of reflection and Snell’s law have to hold). In terms of the ray transit time, we have on S t in 5 tref 5 ttrans , (13a) where t in is the ray transit time associated with the inci- dent ray and tref and ttrans are the transit times for the reflected and the transmitted rays, respectively. One should interpret that the physical meaning of the condi- tion in Eq. (13a) is the conservation of the eikonal phase at the surface of discontinuity. The boundary condition for contour S8 can be deduced by applying Eq. (10) to this discontinuity surface so that S ]t ]x D 2 1 S ]t ]z D 2 5 0. (13b) One possible solution to Eq. (13b) is that (]t/]x)2 5 (]t/]z)2 5 0, with t 5 0, so that at the external sur- face of discontinuity the transit time becomes zero. No- tice that there is another solution, (]t/]x)2 1 (]t/]z)2 5 0, for which the ray transit time becomes imaginary. At this point one has to recall that there is a ray-optical formalism for which the eikonal solution can be complex. This is called Felsen’s general complex solution and in- cludes the evanescent wave in the general solution.29 This leads to a general eikonal equation, F S ]t ~r ! ]x D2 1 S ]t ~r ! ]z D2G 2 F S ]t ~i ! ]x D2 1 S ]t ~i ! ]z D2G 5 n2~x, z ! c2 , (13c) where t (r) and t (i) are the real and the imaginary parts of the ray transit time, respectively. Equation (13c) can be interpreted in the sense of the so-called homogeneous waves, where equiphase surfaces and equiphase contours can be defined.30 Also, if x 5 x0, dx/dz 5 0, as x 5 x(z) would be a con- stant. This represents a constraint or function of the po- sition. Then ]t ]x U x 5 x0 5 n~x0! c . (13d) Note that a more general constraint could be r 5 r(z), with r 5 Ax2 1 y2, which should require the formulation of the boundary conditions in cylindrical coordinates. In the next section we shall prove that n(x0) coincides with the ray-invariant parameter, with x0 the turning point defined in ray trajectories in a gradient-index optical waveguide. 876 J. Opt. Soc. Am. A/Vol. 14, No. 4 /April 1997 M. L. Calvo and V. Lakshminarayanan 4. HOMOGENEOUS MEDIA Consider the simplest case of a homogeneous planar opti- cal waveguide with refractive index n0 defined in the x re- gion, transverse to the longitudinal axis of the waveguide z. Then Eq. (10) simplifies to S ]t ]x D 2 1 S ]t ]z D 2 5 1 v2 , (14) with v 5 c/n0 a constant value, n0 the constant refractive index of the medium, and c the speed of light. The solution of Eq. (14) is the simplest one for a partial differential equation of hyperbolic type.31 Since in the x direction the time is a constant, ]t/]x 5 0, then in Eq. (14) the characteristics are straight lines: t 2 ~z/v ! 5 j, (15a) t 1 ~z/v ! 5 h, (15b) with j and h constants. Equations (15a) and (15b) amount to global constraints for the system (defining a family of lines or linear trajectories). We emphasize this idea with the notation t → t 6 ~z/v !. (16) Also, we note that Eqs. (15a) and (15b) represent Fer- mat’s principle: tmin 5 6E zmin ~1/v !dz, (17) with the initial condition t(z, 0) 5 0. Let us recall here the transit time for the ray trajectory in a homogeneous planar, or slab, waveguide, as predicted by the geometrical optics approximation32: t 5 S n0 c D S 1 cos uz D E zmin dz, (18) with 0 < uz(z) < uC , and uC the complementary of the critical angle. An initial condition over the slope is u(z) 5 u(z0), z0 constant. The initial strategy in this case re- duces to selection of the critical angle. Two main param- eters are introduced in the process of searching for a minimum pathway: (1) angle u between the trajectory and the axial direction and (2) longitudinal variable z. The slope of the trajectory as defined in Eq. (9) is also a constant. The initial global constraint is x 5 z0 tan u. 5. GRADIENT-INDEX OPTICAL WAVEGUIDES Consider an arbitrary gradient-profile planar waveguide with refractive-index distribution following a p law: n2~x ! 5 n0 2@1 2 2D~x/r!p#, (19) with p > 1, D the height of the profile, r the half-width of the waveguide, and n0 the refractive index of the core. In this case Eq. (10) reads S ]t ]x D 2 1 S ]t ]z D 2 5 n2~x ! c2 5 1 v2~x ! , (20) with n(x) following the law defined in Eq. (19). Equation (20) is the equation of the characteristics for the optimized time of the trajectory. It depends on the variable x (cross section) and the z longitudinal direction of the waveguide. The type of solution depends on the value of p. Operating on Eq. (20) and taking into account the con- straints on the system, we arrive at dt dx 5 6S 1c D n~x ! @1 1 z82#1/2 z8, (21) with z8 5 dz/dx the slope of the trajectory, z 5 z(x). This global constraint represents a one-dimensional policy since we are defining z with a single variable. One has to evaluate z 5 z(x) for each particular trajectory, say, for each particular medium. Intuitively, function z 5 z(x) defines the ray trajectory as defined by the geo- metrical optics approximation. The correct formulation provides identical results. Let us analyze Eq. (21). The Euler–Lagrange equa- tion becomes, in this case, d2t dx2 5 S 1c D d dx F n~x ! ~1 1 z82!1/2 z8G 5 0, (22) which implies that n~x ! ~1 1 z82!1/2 z8 5 c1 , (23) with c1 a constant. Then Eq. (23) is z~x ! 5 c1E 1 @n2~x ! 2 c1 2#1/2 dx. (24) Equation (24) coincides with the ray-path parameter characterizing ray propagation in graded-profile waveguides with the condition c1 5 b̄. (25) In Eq. (25), b̄ is the ray-invariant or ray-constant parameter.33 Equation (24) defines the slope of the tra- jectory z8. Defining the first strategy as z 5 z(x) (and z8), we can reasonably start from Eq. (22) to obtain ana- lytical and numerical solutions for the minimization of the quadratic functional. The Euler equation [Eq. (22)] together with the boundary condition [Eq. (25)] assigns a unique solution that gives the absolute minimum of the quadratic functional. The procedure is a variational technique, e.g., the Rayleigh–Ritz method.34 The time minimizing the longitudinal trajectory for the light inside the waveguide is t~z ! 5 t8~zp! 1 1 c E n~z ! ~1 1 z82!1/2 dz. (26) Notice that in Eq. (26) we are integrating over the vari- able trajectory x 5 x(z) and that the constant of integra- tion t8(zp) is an initial condition (e.g., the value for dt/dz is zero) and holds for z 5 zp . For example, for symmetrical profiles, zp represents the half-period (longi- tudinal section) of the trajectory between two turning points. This is not a restrictive condition, since for asym- metrical profiles one could apply the appropriate initial condition for an entire period. The first property to be M. L. Calvo and V. Lakshminarayanan Vol. 14, No. 4 /April 1997 /J. Opt. Soc. Am. A 877 checked for the integrability of Eq. (26) is that the conver- gence of the solution be ensured. This will depend on the form of n(z) and the behavior of z8. Similarly, for the x direction, t~x ! 5 1 c E x L~x, u !dx, (27) with L~x, u ! 5 n~x ! @1 1 u2~x !#1/2 , (28) where u(x) 5 1/z8. u(x) has to be a continuous function with continuous first derivative u8(x). We notice that L(u, x) is a nonlinear differential operator. It has the physical meaning of a differential pathway. Also, the condition u(xtp) 5 0 holds (it is not necessarily a bound- ary condition). The point xtp is similar to the turning point in the x 5 x(z) trajectory. Explicit integration of Eq. (27) depends on a law-clad profile. For example, for p . 2, Eq. (27) could diverge, and one would need to ad- just the physical parameters of the waveguide. Equation (27) is one of the forms of the so-called prob- lem of Mayer35 commonly found in terminal control prob- lems. In the present study, t(x) is, moreover, a solution with fixed end points (defining the dimensions of the waveguide); it has approximate analytical or numerical solutions. 6. PARABOLIC PROFILE AND NUMERICAL RESULTS Assume that p 5 2 in Eq. (19). Then Eq. (24) has an im- mediate analytical solution: z~x ! 5 S b̄r n0A2D D arcsinS x xtp D , (29) with slope z8~x ! 5 A2b̄r 2n0xtpAD@1 2 ~x/xtp!2# . (30) Similarly, Eq. (29) determines x 5 x(z) and u 5 u(x). Figure 3 shows the refractive-index-profile distribution in the X–Z plane representing n@x, x(z)# 5 n(x)n@x(z)# for the following physical parameters: r 5 3 mm, D 5 0.00783, n0 5 1.527. To assign a numerical value to b̄, one should take into account Eq. (25). Since b̄ repre- sents the ray-invariant parameter, as defined in the geo- metrical optics approximation, it has to obey 1.50 < b̄ < 1.527, where 1.50 is the refractive index of the clad- ding and 1.527 the refractive index of the core. We have considered b̄ 5 n(xtp) 5 1.515 in the present study.36 For these values, xtp 5 3 mm and zp 5 74.73 mm. Fig- ure 4 displays a three-dimensional plot of the trajectories z 5 z(x) and x 5 x(z) as given in Eqs. (24) and (25). By representing both functions simultaneously, one appreci- ates a complete description of the spatial ray geometry. Similarly, Fig. 5 gives the shape of the slope of the trajec- tory, z8 5 z8(x, z) [see Eq. (30)], showing the constant be- havior in the z direction as expected. Considering z 5 z(x) as the initial policy and z8 as the global constraint, one can determine numerically the op- timized time t(z) and t(x). The results for the time t(z) with the above physical parameters are shown in Fig. 6(a). This function fits a Fig. 3. Refractive-index distribution n@x,x(z)# 5 n(x)n@x(z)# in the X–Z plane for a planar optical waveguide with p 5 2 [see Eq. (19)]. Physical parameters are n0 5 1.527, D 5 0.0078, r 5 3 mm. The refractive-index distribution is shown along both the transverse and the longitudinal trajectories according to Eq. (29). The sinusoidal dependence in the z direction can be seen. Fig. 4. A three-dimensional plot of the two trajectories z(x) and x(z) as given from Eqs. (24) and (25), giving a complete descrip- tion of the spatial ray geometry. The vertical axis is in arbitrary units. z ranges from 2 to 300 mm, and 22 < x < 12 (mm). b̄ 5 1.151 is the ray-invariant parameter. See text for details. Fig. 5. Slope of the trajectory displayed in Fig. 4, z8 5 dz/dx showing, a gradient in the x direction and constant values in the z direction. 878 J. Opt. Soc. Am. A/Vol. 14, No. 4 /April 1997 M. L. Calvo and V. Lakshminarayanan linear variation for points of the sinusoidal trajectory x 5 x(z). This result is comparable to the ray transit time as predicted by the geometrical optics approximation,32 whose distribution is shown in Fig. 6(b). The correspond- ing slope @]t(x, z)#/]z associated with t(z) is given in Fig. 7(a) for the actual trajectory x 5 x(z). The figure shows the gradual change of the slope for different points of the trajectory, exhibiting a skewlike variation. For the sake of comparison, Fig. 7(b) displays the corresponding slope associated with the ray transit time. The shape of these two functions depends dramatically on the physical pa- rameters of the waveguide, with special influence on zp (the value for which the trajectory has zero slope). The linear time distribution t(z) obtained through dynamic programming appears to have a smaller slope than the one corresponding to the transit time. Fig. 6. (a) Linear dependence of the time of the trajectory t 5 t(z) as given by the application of dynamic programming [see Eq. (26)] with x 5 x(z). For simplicity we consider c 5 1. (b) Ray transit time as predicted by the geometrical optics approxi- mation. The inverse of the local speed n(x)/c (c, free-space speed of light) is integrated along the curved ray path x 5 x(z) (c 5 1). Similarly, Fig. 8 gives the variation of the t(x) function. This exhibits linear dependence, as it does for t 5 t(z). Here the slope @]t(x, z)#/]x coincides with the operator L(x, u) as defined in Eq. (28). These results are pre- sented in Fig. 9. The figure shows how the time compen- sates across the transverse section of the waveguide. It is equivalent to the optimum profile given by Snyder and Love37 for the transit time or equalization. We notice that the minimum value Lmin 5 L(0, u 5 const.) equals the ray-invariant parameter b̄, as expected for an opti- mized profile. It is remarkable that the interval Lmax–Lmin 5 0.01181 (for b̄ 5 1.515) is small but suffi- cient for estimating the pulse width: tmax 2 tmin 5 Dt. The variation of the time compensation concerning b̄ for the boundary points of the waveguide is negligible but be- comes much more appreciable for points on the axis, where the highest value is obtained for b̄ 5 1.50 (coincid- ing with the refractive-index distribution of the cladding). Fig. 7. (a) Slope associated with time t 5 t(z) as given by dy- namic programming [see Fig. 6(a)]. (b) Slope of the ray transit time as displayed in Fig. 6(b). M. L. Calvo and V. Lakshminarayanan Vol. 14, No. 4 /April 1997 /J. Opt. Soc. Am. A 879 Thus the maximum spreading of the pulse corresponds to b̄ 5 ncl and the minimum to b̄ 5 n0. Nonetheless, the differences are not dramatic, because we are considering an optimum profile (p 5 2). As the treatment holds for arbitrary n(x) distribution, it could be generalized to any optical media, provided that the convergence of the solu- Fig. 8. Time distribution across the transverse trajectory (cross section of the waveguide) according to dynamic programming [see Eqs. (27) and (28)]. Fig. 9. Slope of the time variation represented in Fig. 8, coin- ciding with Eq. (28), the differential operator L(x, u). This fig- ure indicates how the time compensates across the transverse section of the waveguide. It is equivalent to the optimum profile given by Snyder and Love37 for the transit time or equalization. Minimum equalization corresponds to the waveguide axis (mini- mum local speed). Maximum equalization corresponds to the boundary of the guide (maximum local speed). We represent L(x, u) for four values of b̄: (1) solid curve, b̄ 5 1.500; (2) dot- ted curve, b̄ 5 1.510; (3) dashed curve, b̄ 5 1.515; (4) dashed– dotted curve, b̄ 5 1.527. See text for details. tion is ensured. The extension to arbitrary profile n 5 n(x, y) (two-dimensional problem) is also possible with the restrictions mentioned above. As a proof of con- sistency, one may check that with use of the numerical re- sult for ]t/]x and ]t/]z together with the n(x) profile, Eq. (20) holds. Evaluating thus one finds values of the order of 2.29–2.33 for both sides of Eq. (20). The discrepancies are due to computational round-off errors. 7. DISCUSSION AND CONCLUSIONS Dynamic programming is a powerful technique for solving optimization problems in many fields of science. The op- timal procedure consists of introducing an initial condi- tion that determines the final state of the system. In the optical problem of defining minimum or optimized path- ways (time minimization), dynamic programming pro- vides immediate solutions in the form of equal-time sur- faces, t 5 t(x, z) > const. Its physical meaning is that of wave fronts generated at the initial point and zero time (zero instant). These basic ideas were first noted by Kalaba,10 who demonstrated that the mathematical pro- cedure optimizing pathways under an initial slope of the trajectory is equivalent to solving an eikonal-type equa- tion. The main result of this paper is the development of the transit-time equation based on a dynamic-programming approach in the problem of light propagation in arbitrary optical media. We have proven the simplicity of this technique by studying various cases of waveguides that have gradient- index profiles. The partial differential equation in terms of ]t/]x and ]t/]z for arbitrary refractive-index distribu- tion is a Hamilton–Jacobi type and has Euler–Lagrange implicit solutions. It automatically gives the initial policy—z 5 z(x), trajectory, and z8 5 dz/dx, slope—as global constraints that determine the final state of the system. The integral form of z(x) equals the one from the eikonal equation, with the additional condition that the constant of integration be the ray-invariant b̄. The solutions for t 5 t(z) have behavior equivalent to the ray transit time predicted by the geometrical optics approxi- mation. The numerical estimates carried out demon- strate this statement clearly. Moreover, the slope of the trajectory t 5 t(x) has meaning similar to the optimum profile distribution under the geometrical optics approxi- mation. The slope represents the transit-time variation with the ray invariant b̄. The results indicate that the value of the differential operator L(x, u) at the origin (axial direction) coincides with b̄ 5 1.151 (as in the present study). This curve displays the time equaliza- tion for specific trajectories and may be useful for design- ing specific optical media. It also provides numerical val- ues for the pulse spread. The inverse problem might be analyzed—for example, searching for optimized optical media (fastest propagation of the luminous signal with the lowest distortion). Possible restrictions come from the analytical behavior of the integrals, and one needs to check the convergence of the solution. For example, we analyzed additional results for the law-clad profile p 5 4 (not discussed here for brevity), obtaining divergence for some values of the physical parameters. The method 880 J. Opt. Soc. Am. A/Vol. 14, No. 4 /April 1997 M. L. Calvo and V. Lakshminarayanan is applicable only in the context in which Fermat’s prin- ciple holds. It appears to have an equivalence in quan- tum mechanics. In this context, solutions of the Hamilton–Jacobi-type equations are assured as per Liou- ville’s theorem. Exact solutions are valuable only for specific potentials (having approximately an equivalence with the squared refractive-index distribution). These aspects are now under investigation and will be presented in forthcoming papers. APPENDIX A We look for a compact form of Eq. (11): F S ]e ]x D 2 1 S ]e ]z D 2G 1 F e~x, z ! v~x, z !G 2F S ]v ]x D 2 1 S ]v ]z D 2G 2 2F e~x, z ! v~x, z !GF S ]e ]x D S ]v ]x D 1 S ]e ]z D S ]v ]z D G 5 1. (A1) Dividing both sides of Eq. (A1) by 1/e2, we get F S ] ln e ]x D 2 1 S ] ln e ]z D 2G 1 F S ] ln v ]x D 2 1 S ] ln v ]z D 2G 2 2F] ln e ]x ] ln v ]x 1 ] ln e ]z ] ln v ]z G 5 1 e2 . (A2) Applying relations of vectorial analysis Eq. (A2) produces Ue ¹̄S ln e v D U 2 5 1. (A3) One can object that in Eq. (A3) e/v is not a dimensionless quantity (since it is not mathematically correct to repre- sent a logarithm with dimensions). To avoid this incon- sistency we introduce constant values of e and v—ei and vi—for certain given initial coordinates x 5 xi , z 5 zi , e0 5 e(xi , zi), v0 5 vi(xi , zi). The introduction of these two constants does not alter Eq. (A3). Then Ue ¹̄ lnS ee0 v0v D U2 5 1, (A4) or F ] ]x lnS ee0 v0v D G2 1 F ] ]z lnS ee0 v0v D G2 5 1 e2 , (A5) which is the same as Eq. (12). ACKNOWLEDGMENTS This work was supported in part by the School of Optom- etry of the University of Missouri–St. Louis (UMSL). A UMSL research grant is also gratefully acknowledged. We also thank Jerry Christensen, Dean, School of Optom- etry, USML, for support during the initial stages of this research program. Additional financial assistance from the Spanish Ministry of Education and Science (DGICYT) under project PR94-362 to M. L. Calvo is also acknowl- edged. Partial results were presented at the OSA An- nual Meeting (Portland, Oregon, October 1995) and at the ICO-XVII Congress (Taejon, Korea, August 1996). The authors also thank an anonymous referee for a careful re- view and R. F. Alvarez-Estrada and G. F. Calvo for help- ful suggestions and discussions. M. L. Calvo can be contacted as follows: tel, 341-394- 4684; fax, 341-394-4683; e-mail, mlcalvo@eucmax.sim. ucm.es. REFERENCES AND NOTES 1. R. K. Lagu and R. V. Ramaswamy, ‘‘A variational finite dif- ference method for analyzing channel waveguides with ar- bitrary index profiles,’’ IEEE J. Quantum Electron. QE-22, 968–978 (1986). 2. A. Sharma and P. Bindal, ‘‘Variational analysis of diffused planar and channel waveguides and directional couplers,’’ J. Opt. Soc. Am. A 11, 2244–2248 (1994). 3. A. Sharma, D. V. Kumar, and A. K. Ghatak, ‘‘Tracing rays through graded-index media: a new method,’’ Appl. Opt. 21, 947–987 (1982). 4. B. D. Stone and G. W. Forbes, ‘‘Optimal interpolants for Runge–Kutta ray tracing in inhomogeneous media,’’ J. Opt. Soc. Am. A 7, 248–254 (1990). 5. J. Puchalsky, ‘‘Numerical determination of continuous ray tracing: the four-component method,’’ Appl. Opt. 33, 1900–1906 (1994). 6. K. B. Wolf and G. Krötzen, ‘‘Geometry and dynamics in re- fracting systems,’’ Eur. J. Phys. 16, 14–20 (1995). 7. A. Ghatak, E. Sharma, and J. Kompella, ‘‘Exact ray paths in bent waveguides,’’ Appl. Opt. 27, 3180–3184 (1988). 8. A. Ghatak and E. G. Sauter, ‘‘The harmonic oscillator prob- lem and the parabolic index optical waveguide I. Classical and ray optics analysis,’’ Eur. J. Phys. 10, 136–143 (1989). 9. E. G. Sauter and A. K. Ghatak, ‘‘The harmonic oscillator problem and the parabolic index optical waveguide II. Quantum mechanical and wave optical analysis,’’ Eur. J. Phys. 10, 144–150 (1989). 10. R. Kalaba, ‘‘Dynamic programming, Fermat’s principle, and the eikonal equation,’’ J. Opt. Soc. Am. 51, 1150–1151 (1961). See also V. Lakshminarayanan and S. Varadhara- jan, ‘‘Dynamic programming, Fermat’s principle, and the ei- konal equation—revisited,’’ J. Optimization Theor. Appl. (to be published). 11. R. E. Bellman, Dynamic Programming (Princeton U. Press, Princeton, N.J., 1957). 12. V. Lakshminarayanan, S. Varadharajan, and M. L. Calvo, ‘‘A note on the applicability of dynamic programming to waveguide problems,’’ in Photonics ’96: Proceedings of the International Conference on Fiber Optics and Photonics, J. P. Raina and P. R. Vaya, eds. (Tata McGraw-Hill, New Delhi, 1977), Vol. 1, pp. 209–214. 13. R. E. Bellman and R. Kalaba, ‘‘Dynamic programming, in- variant imbedding and quasi-linearization: comparison and interconnections,’’ in Computing Methods in Optimiza- tion Problems, A. V. Balakrishnan and L. W. Neustadt, eds. (Academic, New York, 1964), pp. 135–145. 14. S. E. Dreyfus, Dynamic Programming and the Calculus of Variation (Academic, New York, 1965). 15. R. E. Bellman and R. Vasudevan, Wave Propagation. An Invariant Imbedding Approach (Reidel, Dordrecht, The Netherlands, 1986). 16. S. E. Dreyfus and A. M. Law, The Art and Theory of Dy- namic Programming (Academic, New York, 1977). 17. R. E. Bellman, R. Kalaba, and G. M. Wing, ‘‘Invariant im- bedding and mathematical physics I. Particle process,’’ J. Math. Phys. 1, 280–308 (1960). 18. R. E. Bellman, Introduction to the Mathematical Theory of Control Processes. Linear Equations and Quadratic Crite- ria (Academic, New York, 1967), Chaps. 3 and 4. 19. R. E. Bellman and S. E. Dreyfus, Applied Dynamic Pro- gramming (Princeton U. Press, Princeton, N.J. 1962), Chap. VIII. 20. M. Born and E. Wolf, Principles of Optics, 6th ed. (Perga- mon, Oxford, 1984), Sec. 4.2. 21. Caratheodory’s theorem establishes the connection between Fermat’s principle and the calculus of variation de- M. L. Calvo and V. Lakshminarayanan Vol. 14, No. 4 /April 1997 /J. Opt. Soc. Am. A 881 fining precise extremals. See C. Caratheodory, Geome- trische Optik (Springer, Berlin, 1937). See also Ref. 20, App. 1, pp. 719–737. 22. M. L. Calvo and V. Lakshminarayanan, ‘‘Dynamic program- ming: an alternative approach to light propagation in ar- bitrary optical media,’’ in 17th Congress of the International Commission for Optics: Optics for Science and New Tech- nology, J. S. Chang, J. H. Lee, S. Y. Lee, and C. H. Nam, eds., Proc. SPIE, 2778, 294–295 (1996). 23. See, for example, Ref. 20, p. 111. 24. R. F. Alvarez-Estrada, Department of Theoretical Physics, Faculty of Physical Sciences, Universidad Complutense, 28040 Madrid, Spain, June 1996 (personal communication). 25. Toda’s potential, also called the Toda lattice, is a well- known solution of Eq. (12) in the field of quantum mechan- ics and nonlinear waves. See M. Toda, ‘‘Vibration of a chain with nonlinear interaction,’’ J. Phys. Soc. Jpn. 22, 431–436 (1967); ‘‘Wave propagation in anharmonic lat- tices,’’ J. Phys. Soc. Jpn. 23, 501–506 (1967); A. Khare and R. K. Bhaduri, ‘‘Exactly solvable noncentral potentials in two and three dimensions,’’ Am. J. Phys. 62, 1008–1014 (1994). 26. R. K. Luneburg, The Mathematical Theory of Optics (U. California Press, Berkeley, Calif., 1964), p. 180. 27. S. Cornbleet, Microwave Optics, The Optics of Microwave Antenna Design (Academic, New York, 1976), pp. 127–132. See also Sec. 2.6. 28. E. T. Whittaker, A Treatise on the Analytical Dynamics of Particles and Rigid Bodies, 4th ed. (Cambridge U. Press, London, 1964), Chap. X, p. 281. 29. L. B. Felsen and M. Marcuwitz, Radiation and Scattering of Waves (Prentice-Hall, Englewood Cliffs, N.J., 1973). 30. Some interesting formalism related to the so-called homogeneous-wave solution can be found in S. Solimeno, B. Crosignani, and P. Di Porto, Guiding, Diffraction and Con- finement of Optical Radiation (Academic, New York, 1986), pp. 63–66. 31. J. Mathews and R. L. Walker, Mathematical Methods of Physics (W. A. Benjamin, New York, 1970), Sec. 8–2. 32. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman & Hall, London, 1983), Sec. 1–9. 33. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman & Hall, London, 1983), Sec. 1–8. 34. Although the Rayleigh–Ritz method is discussed in the lit- erature in relation to variational principles, here we intro- duce applications in the context of dynamic programming. See, for example, Ref. 15, Chap. X, Sec. 3. 35. G. A. Bliss, Lectures on the Calculus of Variation (U. of Chi- cago Press, Chicago, Ill., 1959). 36. The constant of integration b̄ plays an interesting role in the dynamic-programming formulation. It ensures the similarity between the trajectory z(x) given in Eq. (24) and the ray propagation in graded-index-profile waveguides un- der the geometrical optics approximation. Even if this similarity were to be ignored, there should appear values of b̄ > n0, n0 being the refractive index of the core, for which the integral in Eq. (24) diverges. This indicates that only bound trajectories (bound rays) are convergent solutions of Eq. (24). 37. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman & Hall, London, 1983), Sec. 3–2.