1. Introduction

A paraxial optical system (POS) is a system in which rays propagate forward with a small inclination angle and a small transverse distance about the optic axis. Generally, ray propagation obeys the Fermat principle, and solving this equation is a complicated process in most cases. As the paraxial approximation can greatly simplify the theoretical analysis of ray dynamics problems when the light source and ray trajectory is close to the optic axis, the POS theory is widely applied for optical design, such as image formation [1,2], optical transformation [3,4], paraxial beam analysis [5–7] and aberration analysis [8].

The mathematical tool to analysis the state and transformation of a beam in POS is the ray transfer matrix, also called ABCD matrix, which was first introduced by Kogelnik and Li to study the incident and outgoing laser beams and optical resonators [9]. It is a very important technique in treating problems of propagation and diffraction [6,7,10–13], transformation [5,14–16] and self-consistent solution [17] for the input and output rays of the optical elements. Locating a light spot on the ray trajectory and mark its propagation direction is to apply a two-dimensional state vector ${\textbf u} = \left( {\begin{array}{{c}} {y(x)}\\ {\theta (x)} \end{array}} \right)$, where y represents its transverse coordinate and θ represents the inclination angle of the ray about optic axis (x axis). For example, the state vectors for input and output are shown in Fig. 1. When light propagation through a linear optical element, the relationship between input ${{\textbf u}_0}\textrm{ = }\left( {\begin{array}{{c}} {{y_0}}\\ {{\theta_0}} \end{array}} \right)$ and output ${{\textbf u}_1}\textrm{ = }\left( {\begin{array}{{c}} {{y_1}}\\ {{\theta_1}} \end{array}} \right)$ is described by the following linear algebra equation

 

Fig. 1. Ray propagation through a linear optical element. The state vector of a light spot is given by Eq. (1). The subscripts 0 and 1 respectively represent the input and output rays. The states between the input and output rays are linked through the ABCD matrix of the optical element.

Download Full Size | PDF

As seen in Eq. (1), the matrix elements are determined by

From Fig. 1, we see the conditions ${\theta _0} = 0$ and ${y_0} = 0$ represent two special rays. The former is parallel to the optic axis, while the incident point of the latter is located on the optic axis.

If the light propagates through a family of optical elements, the total transfer matrix is expressed as a multiplication of all ABCD matrices

As a convenient tool for the calculation of light propagation in POS, the restriction of ABCD matrix also exists. First, as shown in Eq. (3), it is usually inadvisable for an optical medium with continuous refractive index function, especially the input and output present a more complicated (probably nonlinear) relationship. Although quite a large portion of nonlinear systems can be linearized by keeping the first-order term and drop the high-order ones [18], the linearization process does not always go swimmingly, because the Jacobi matrix $\frac{{\partial {{\textbf u}_1}}}{{\partial {{\textbf u}_0}}}|{_0} $ of input-output mapping may be singular. As the optical path is reversible, the linearization should ensure it is an invertible linear transformation. This is a ubiquitous problem for the gradient-refractive-index (GRIN) medium [19–23]. If the Jacobi matrix $\frac{{\partial {{\textbf u}_1}}}{{\partial {{\textbf u}_0}}}|{_0} $ is not invertible, the POS can no longer be treated as a linear system and the ABCD matrix does not exist. One has to consider the high-order terms. Since analyzing the high-order terms in series expansion is laborious, it is necessary to seek for a uniform equation including the refractive index instead of expanding it as a power series [22,23]. Second, calculating the elements of a single matrix, according to Eq. (2), should generally consider two special rays. If an optical system consists of several elements, the calculation of the whole ABCD matrix will be more costly. Moreover, this procedure is only effective for linear system. Third, the ABCD matrix regards an optical element as a black box, from which we only know the relationship between the input and output, but it does not show the details that how the ray evolves inside. Due to the complexity and sensitivity of optical systems, how to control the fluctuation of light path is a potential challenge for optic design. This is not only shown from the input and output, but also related to the inner structure. Thus tracking the ray trajectory and analyzing the stability of a POS is a necessary work [24–26]. However, to the best of our knowledge, this has not been paid much attention for a continuous optical system. To solve these problems, the crucial point is to extend the evolution equation of a state vector from a family of discrete optical elements to a continuous form. It is reasonable to predict that Eq. (1), as an algebra equation, is no longer a proper description in a continuous system. In order to describe the light propagation in a general POS, one can predict that an ordinary differential equation (ODE) model, as a replacement of Eqs. (1) and (3), may be more suitable for continuous cases.

From the perspective of Lagrangian optics [27], if x axis is analogous to a generalized time axis, and ray trajectory to a generalized time evolution pattern, one can regard y as a generalized coordinate and θ ≈ tan θ = y’ as a generalized velocity. Therefore, the state vector u(x) represents a point in phase space of a POS, while the ray propagation problem is ascribed to how a state vector evolves in the phase space of a geometric optical system [28–30]. Therefore, the problem here is to develop an equation which can describe phase evolution under the paraxial approximation, and as an extension of ABCD matrix, it should be expressed as a consistent form for both continuous and layered medium. The first component of phase vector y(x) describes the ray trajectory, and the whole vector u(x), as a parameter equation of y(x) and θ(x), describes the phase trajectory. All the information about ray propagation is included in the phase vector, which is the goal we need to achieve.

This paper is organized as follows: In Sec. 2, we derive the state evolution formula (SEF) of a POS based on the optical Lagrange equation. In Sec. 3, we put forward some application examples to solve the phase vector by using SEF. In Sec. 4, we briefly introduce the concept of stability of a POS based on SEF, and show the process to check whether an optical system is stable. Sec. 5 summarizes this paper.

2. State evolution formula

For convenience of the following mathematical treatment, there should be a statement before the presentation of our theory. For all functions appearing in this paper, every jump point is considered as continuous, and its value takes the average of left and right limits. The algorithm of these functions conforms to that of continuous ones. This regulation is in accordance with the Dirichlet convergence theorem, and these functions can be regarded as a limitation of continuous function sequences. For example, the Heaviside step function is defined as

Although H(x) is discontinuous at x = 0, it can be viewed as a limitation of continuous function sequence ${H_n}(x) = \frac{1}{2} + \frac{1}{\pi }\arctan nx,\;n \to \infty $. The variation tendency of H_n (x) is shown in Fig. 2. Since H_n (x) is continuous, the order of operation about limitation and integral can be exchanged. Therefore we have $\int_{ - \infty }^{ + \infty } {H(x)\delta (x)dx} = \mathop {\lim }\limits_{n \to \infty } \int_{ - \infty }^{ + \infty } {{H_n}(x)\delta (x)dx} = \frac{1}{2}$, where $\delta (x)$ is the Dirac function. Similar regulation holds for other discontinuous functions. This procedure will greatly facilitate the analysis of the layered medium with discontinuous refractive index functions.

 

Fig. 2. Variation tendency of continuous function sequence H_n (x), n = 1, 5, 10 and ∞.

Download Full Size | PDF

We begin our discussion from the optical Lagrange equation [27], which describes the geometric trajectory of a ray propagating in a medium with refractive index function $n = n(x,y)$

Undoubtedly, it can be also derived from the optical Hamilton equations [31]. Equation (6) can be further simplified through performing paraxial approximation ${y^{{\prime}^2}} \ll 1$

The error of this approximation can be estimated. If we assume that |θ|≤10° (≈ 0.174) in the paraxial region, and then ${y^{{\prime}^2}} \le {\tan ^2}{10^ \circ } \approx 0.03 < < 1$, the error is limited within 3%. Therefore it can get a relatively exact result.

Since n may be relevant to y, Eq. (7) is generally a nonlinear equation. In the previous work which assumes the medium is axial-symmetric (or concentric), the linearization procedure is to expand $n = n(x,y)$ as a power series, $n = {n_0}(x) + {n_1}(x){y^2} + {n_2}(x){y^4} + O(6)$ [22]. However, the first-order approximation may lose the GRIN properties when ${n_0}(x) \equiv C$ and ${n_1}(x) \equiv 0$ because of the singularity of linearization. To avoid this problem and seek for a general equation for a distributed refractive index function, here we define two parameters ${N_1}(x,y) = \frac{{\partial \ln n}}{{\partial x}}$ and ${N_2}(x,y) = \frac{{\partial \ln n}}{{\partial y}}$, called axial and transverse logarithmic refractive index coefficients (axial LRIC and transverse LRIC) respectively, expressing the variation of logarithmic refractive index in unit length. Then we replace each total derivative about x with an over-dotted sign. Note that $\dot{y} = y^{\prime} = \theta $ under the paraxial approximation, one can cast Eq. (7) into a phase-space form

Equation (8) is the main result of this paper, giving the state evolution formula (SEF) of a paraxial geometric optical system. Compared with Eqs. (1) and (3), it tells how a phase vector on the ray trajectory evolves continuously in a medium with a distributed refractive index. Therefore, it extends the discrete ABCD matrix method to a continuous case. Through this formula, we are capable to track an arbitrary light spot on the ray trajectory, and do not need to consider the special rays to solve the ABCD matrix of a black box. Moreover, one can predict that even for a POS that cannot be linearized, the relation between input and output still keeps the form of Eq. (8). This overcomes the difficulty of expanding the refractive index as a power series, and also avoids the singularity problem during the linearization process [22,23]. It is worth to mention that Eq. (8) is applicable for an arbitrary refractive index function, and this eliminates the restriction that the medium must be axial-symmetric [22]. By transforming the jump function into a limitation of a continuous function sequence, it takes a consistent form for continuous and discontinuous cases, which can be treated uniformly rather than separately. The initial condition of Eq. (8), namely input state, which gives the position and direction of the light source (or equivalent light source), is

Equations (8) and (9) construct a Cauchy problem of ODE, and its solution ${\textbf u} = {\textbf u}(x,{x_0},{{\textbf u}_0})$ can uniquely determine the evolution of phase vector and ray trajectory.

In the following sections, we will introduce some application examples of this formula, including solving the phase vector and analysis the stability of a POS.

3. Application: solution of state evolution formula

3.1 Homogeneous medium

For a homogenous medium, $n(x,y) = C$, the LRICs ${N_1} = {N_2} = 0$, from which Eq. (8) takes the form

The solution of Eq. (10) is derived through a direct integral under the initial condition Eq. (9)

It can be seen that M₁ is in accordance with the ABCD matrix of a homogenous medium [5,9]. This verifies the validity of the SEF.

3.2 Axial inhomogeneous medium

We consider a medium with axial inhomogeneous and smooth refractive index function, $n(x,y) = n(x)$, which leads to ${N_1} = \frac{{\dot{n}(x)}}{{n(x)}},{N_2} = 0$. Then Eq. (7) becomes

The second component of Eq. (13) can be integrated under initial condition Eq. (9), $\theta = \frac{{n({x_0})}}{{n(x)}}{\theta _0}$, which gives the solution

This is the evolution law of a ray propagating in an axial inhomogeneous medium. Moreover, we can also rewrite Eq. (14) as a matrix form

Therefore M₂ is the ABCD matrix of an axial inhomogeneous medium. In particular, if $n(x) \equiv n({x_0})$, Eq. (14) will be transformed into Eq. (11). This is in accordance with the homogeneous case.

3.3 Spherical interface

Despite the state evolution formula is proposed from continuous medium, it can be shown that Eq. (8) also goes for discontinuous cases. As stated above, this only needs to suitably rewrite the refractive index function by using the Heaviside function in Eq. (4). For example, the ray propagation through a spherical interface is shown in Fig. 3. The interface is between two homogenous medium with refractive indices n₁ and n₂ respectively.

 

Fig. 3. Schematic diagram of a ray passing through a spherical interface with radius R. The two homogenous medium with refractive indices n₁ and n₂ are astride the interface.

Download Full Size | PDF

The refractive index function of a spherical interface with radius R, in the paraxial region, takes the form

Since the media before and after the interface are both homogenous, we are interested in the transition of the state vector astride the interface. Here the superscripts – and + are used to denote the states of the ray trajectory behind and after the interface respectively. On the left side of the interface, one can rewrite the initial condition as $\left( {\begin{array}{{c}} {y(x_0^ - )}\\ {\theta (x_0^ - )} \end{array}} \right) = \left( {\begin{array}{{c}} {{y_0}}\\ {{\theta_0}} \end{array}} \right)$. Obviously, the ray trajectory satisfies the continuous condition $y(x_0^ + ) = y(x_0^ - ) = {y_0}$. In order to obtain the relationship between $\theta (x_0^ + )$ and $\theta (x_0^ - )$, one can perform integral $\int_{{x_0} - \varepsilon }^{{x_0} + \varepsilon } {( \cdots )\textrm{d}x} $ in the second component of Eq. (17), and let $\varepsilon \to {0^\textrm{ + }}$. We get $\theta (x_0^ + ) - \theta (x_0^ - ) ={-} 2\frac{{{n_2} - {n_1}}}{{{n_2} + {n_1}}}\left( {\frac{{{y_0}}}{R}\textrm{ + }\frac{{\theta (x_0^ + ) + \theta (x_0^ - )}}{2}} \right)$, which can be further simplified as

One can also prove that Eq. (18) is valid. The transition of the phase vector on the interface can be presented as a matrix form

This result conforms to the ABCD matrix of a spherical interface. If $R \to \infty $, it becomes a planar interface, and Eq. (18) gives ${n_2}\theta (x_0^ + ) = {n_1}\theta (x_0^ - )$, which is exactly the Snell’s law.

3.4 Transverse inhomogeneous medium

Now we consider a nonlinear element whose refractive index function takes a transverse inhomogeneous form, $n(x,y) = n(y)$. Then ${N_1} = 0,{N_2} = \frac{{\textrm{d}\ln n}}{{\textrm{d}y}}$, and Eq. (8) becomes

The second equation of Eq. (20) can be transformed as $\frac{{\textrm{d}(\ln y)}}{{\textrm{d}y}} = \dot{\theta } = \frac{{\textrm{d}\theta }}{{\textrm{d}y}}\dot{y} = \theta \frac{{\textrm{d}\theta }}{{\textrm{d}y}} = \frac{1}{2}\frac{{\textrm{d}({\theta ^2})}}{{\textrm{d}y}}$. From the initial condition Eq. (9) one can get the first integral of Eq. (20)

Then y can be solved from the first component of Eq. (20)

The integral can be calculated once n(y) is given. The sign is determined by the direction of light source. Since the phase vectors of input and output present a nonlinear relation, it is impossible to apply the ABCD matrix in such a case.

From the above examples, it can be seen that the phase vector of a POS can be directly solved from the SEF without considering the special rays. This is not only applicable for linear and discrete ABCD systems [5,7,12], but also for nonlinear and continuous systems. Instead of calculating four matrix elements respectively, solving SEF is relatively simple as it is a one-step process to get the evolution of phase vector. In addition, the relationship between input and output in nonlinear systems cannot be described by the ABCD matrix, while the SEF can cope with such situations.

4. Application: stability of a paraxial optical system

A POS with an arbitrary refractive index function, in general, is a nonlinear dynamic system, as seen in Eq. (8). A perturbation imposed upon the light source of a geometric optical device may cause a deviation on the ray trajectory by changing its position ${y_0}$ and direction ${\theta _0}$. If the perturbation on the state of light source is small, the deflection of ray trajectory is also small, we say that the system is stable, and otherwise it is unstable, as shown in the schematic Fig. 4. An unstable system is usually low-performance because it may lead to the ray trajectory drastically deflects from the principle optical axis. This may break the paraxial properties and make the system be out of control, even under a small perturbation. Therefore, the stability should be fully considered in the optic design.

 

Fig. 4. Schematic diagram for the stability of a POS. If a small perturbation is imposed upon the light source, making its position and direction be changed, the ray trajectory does not deflect drastically, the system is stable, and otherwise it is unstable.

Download Full Size | PDF

Previous researches have discussed the stability of a discrete and periodic optical system [24,25]. Periodic POS with purely real ray matrices are classified as either stable or unstable from the semi-trace criterion, depending on the absolute value of the semi-trace of its ABCD matrix is greater or less than 1 [9,24]. However, a general case should consider a continuous and a non-periodic system, which has not been previously investigated. Here we define the stability of a general POS from the SEF. Assume the perturbation is fixed in the reference plane of light source $x = {x_0}$, and if not, one can use an equivalent light source for substitution. The unperturbed state of light source u₀ is described by Eq. (9), under which the state vector of light spot u is the solution of Eq. (8). If a perturbation occurs upon the light source, shifting the initial state ${{\textbf u}_0}$ to ${\tilde{{\textbf u}}_0}$, the state vector of ray will correspondingly deflect from ${\textbf u}$ to $\tilde{{\textbf u}}$, as shown in Fig. 4. One can define the deviation vector $\delta {\textbf u} = \tilde{{\textbf u}} - {\textbf u}$ to describe the deflection after the perturbation. Note that both ${\textbf u}$ and $\tilde{{\textbf u}}$ are the solutions of Eq. (8), only with different initial conditions. Substituting them into Eq. (8) and subtracting, we get the following ODE

The initial condition of Eq. (23) is the deviation vector of light source

It can be seen that $\delta {\textbf u} \equiv {\textbf 0}$ is always the solution of Eq. (23). The physical significance of this trivial solution is that the light trace has no perturbation, and thus no deviation. Therefore, if light source is small perturbed, the deviation vector $\delta {\textbf u}$ just fluctuates around zero instead of drastic deflection, the system is stable. This can be concluded to the following mathematical statement

In other words, if the deviation of source $||{\delta {{\textbf u}_0}} ||$ is restricted into a small range ($< M$), the deviation of whole trajectory $||{\delta {\textbf u}} ||$ can also reach a sufficiently small extent ($< \varepsilon $), we can say the system is stable. Equation (25) is the strict definition of a stable POS. Furthermore, if the deviation vector satisfies

Since the deviation vector is the variation of phase vector, performing variation can quickly judge the stability if the solution is known. For example, the deviation vector in a homogenous medium can be obtained from Eq. (11)

As seen, even a small perturbation $\delta {\theta _0}$ upon the direction of light source will lead to a boundless deviation on the ray trajectory when x − x₀ is large enough. So the homogenous medium with an infinite length, theoretically, is an unstable optical system. Fortunately, the homogeneous medium always has a finite length in practice, and x − x₀ is bounded. If propagation length is sufficiently small, the small deviation will not be enlarged to a great extent, and the ray will not stray from the paraxial region. In such a case the perturbation can be ignored.

Another example is the axial inhomogeneous medium, of which the deviation vector can be obtained from Eq. (14)

As seen, if axial LRIC ${N_1} < 0$, the medium becomes rarer along x axis (an extreme case is $\frac{{n({x_0})}}{{n(x)}} \to \infty$), the system is unstable, since the factor $\frac{{n({x_0})}}{{n(x)}}$ will enlarge the direction deviation of light source $\delta {\theta _0}$. Therefore, to obtain a stable system, the medium should become denser along the propagation direction. According to Eq. (15), this conclusion is in accordance with semi-trace criterion of periodic systems [9,24]. However, for a continuous and non-periodic system, it is only a necessary condition and cannot ensure its stability.

Explicitly, we set $n(x)\sim {x^p},p > 0$ and $n({x_0}) = 1$. If p > 1, the system is stable as the integral is bounded. However, it is not asymptotically stable as Eq. (26) is not satisfied. If p ≤ 1, the integral is boundless when $x \to \infty$, yielding an unstable system. The ray trajectories for p = 2, p = 1 and p = 0.8 are depicted in Fig. 5. Here, we set the initial condition as ${{\textbf u}_0} = \left( {\begin{array}{{c}} {0.1}\\ {0.1} \end{array}} \right)$ and ${\tilde{{\textbf u}}_0} = \left( {\begin{array}{{c}} {0.11}\\ {0.11} \end{array}} \right)$, which ensures that the transverse distance y and the inclination angle θ are in the paraxial region. As seen, the deviation between two trajectories will tend to a constant value in a stable system (Fig. 5(a)), but be dramatically enlarged with the ray propagation in an unstable system (Fig. 5(b) and (c)). Through slightly changing the position and direction of the light source, we also numerically simulated the ray trajectories in the cases of Fig. 5(a) - (c), and the related results are correspondingly depicted in Fig. 5(d)–5(f). The analytical and numerical results are mostly matched.

 

Fig. 5. Ray trajectories before and after perturbation in an axial inhomogeneous medium. (a) ∼ (c) are analytical results, and (d) ∼ (f) are light path simulation results with unperturbed and perturbed sources. ${x_0} = 1$, ${{\textbf u}_0} = \left( {\begin{array}{{c}} {0.1}\\ {0.1} \end{array}} \right)$, and ${\tilde{{\textbf u}}_0} = \left( {\begin{array}{{c}} {0.11}\\ {0.11} \end{array}} \right)$. (a) (d) p = 2, stable; (b) (e) p = 1, unstable; (c) (f) p = 0.8, unstable. The analytical and numerical results are well matched.

Download Full Size | PDF

5. Conclusion

In conclusion, we put forward a systematic theory of the phase vector evolution of a ray propagating in POS. By performing the paraxial approximation based on optical Lagrange equation, we formulated an ODE model, called SEF, to describe the evolution of phase point and ray trajectory. Compared with traditional ABCD matrix method, SEF is more effective and versatile since it is suitable for optical medium with both continuous and discontinuous refractive index functions, even the input and output present a nonlinear relationship. Compared with the ray dynamic equation proposed previously, it eliminates the restriction that the medium is axial symmetric, and avoids the singularity problem of linearization. It also proves correct from the existing results of discrete geometric optical elements. Based on SEF, we propose a strict definition for stability of a POS and obtain the stability criterion, through which one can judge whether a POS can mostly hold its light path when the source is under small perturbation. Our study may give a reference model for the theoretical analysis of ray dynamics in optical devices and systems.