Resonant Transmission Line Method for Unconventional Fibers

Abstract

We provide a very general review of the resonant transmission line method for optical fiber problems. The method has been found to work seamlessly for a variety of difficult problems including elliptical and eccentric core fibers as well as “holey” photonic crystal fibers. This new version has been shown to offer great versatility with respect to cases of unconventional, inhomogeneous index profiles.

1. Introduction

The method of resonant transmission lines (RTL) has evolved from some original observations around the generic theme of transverse resonance in direct association with the standard theory of the solutions of the general telegrapher’s equation in transmission lines. The particular advantage of the method is that it offers the possibility of reducing even full three-dimensional problems into one its dimensional equivalent, at least wherever a full variable separation is possible due to the particular symmetries of a given coordinate system. A full historical review of the use of transmission lines for modeling electromagnetic problems would require too much space but we can provide some important milestones in chronological order in brief.

The transmission line theme can be traced as far back as in the old treatise by Schelkunoff [1]. Further physical justification is offered by a seminal paper by Marcuvitz and Schwinger [2] where we find a first version of the transverse resonance condition as a means to satisfy certain boundary conditions via appropriate impedance matching. Extensive use of the telegrapher’s equation is also found in Ramachandran [3], Gallawa [4], Clarricoats and Oliner [5,6], Yoneyama [7], Borneman and Arndt [8], Tao [9], Shigesawa and Tsuji [10] as well as Dahl et al. [11]. More recently, Moshonas et al. [12] also used transmission line methods for the study of photovoltaic cells. An important comparison between transverse stationary modes and their similarity with quantum wells is mentioned by Bialynicki-Birula [13] and has led to the application of the RTL method by some of the authors in Schrödinger type and general Sturm–Liouville problems [14,15,16]. An early application of similar ideas in optical fibers can be traced in a paper by Yeh and Lindgren [17]. Tamir also discussed guided waves based on impedance matching in Reference [18], while later Carlin and Zmuda [19] applied transfer matrix methods for inhomogeneous fibers. More recently, Mencarelli and Rossi [20] used transverse resonance for studying multilayered photonic crystals.

The modern RTL method was originally introduced by Papageorgiou and Boucouvalas [21,22,23,24,25,26] in simple cases of standard step index and radial index profiles. Lately, it was also applied successfully in the case of superlattices with periodic potentials [27] which are similar to a certain class of photonic crystal fibers (PCFs) or Bragg fibers in general. Several other cases of unusual refractive index profiles have been successfully treated with the same method [28,29,30,31]. In the present work, we expose more recent developments in an effort to expand the validity of the RTL method in all unconventional fiber models.

We shall heretofore refer to all the cylindrical optical fibers as the class of conventional optical fibers (COFs) for as long as they can be separated in a set of n very thin successive cylindrical layers of average radius r with uniform refractive index values given as $n_i$. Indices of successive layers will be allowed in general to be different for each step varying as $n_i=n(r_i)$. The resulting, total index profile n(r) variation from the center of the fiber up to the limit outer air medium, completely defines the propagation properties inside the fiber.

We shall also separate another class of unconventional optical fibers (UOFs) as the ones in which at least in some of its successive thin cylindrical layers present an additional variation of the total index profile along its radial coordinate φ in the form η(r,φ). Such cases include elliptic core, non-symmetric or eccentric core fibers, and in general, all cases of not strictly circular core fibers. In these cases, any discretization into thin cylindrical layers that “cuts” through the core and cladding, results in a variation of the refractive index along φ. In Figure 1, we show such a case of a thin, circular cylindrical layer cutting an elliptic core fiber. As is evident in the schematic, the variation is analogous to the arc length of the circular sector cut by the ellipse for each discretization step.

Additionally, photonic fibers made of silica with a set of small air holes around their centers will have many cylindrical layers of varying refractive indexes along θ. In Figure 2, we show an example of a PCF with a hexagonal lattice of 6 rows of air holes.

Again, radial discretization with thin, circular cylindrical layers results in alternating cuts through air holes and silica of a PCF as depicted in the schematic. The main aim in the present work is to develop the original RTL method so as to be able to include all such UOF cases via an appropriate transformation to mathematically equivalent COF cases.

In the following sections, we first introduce the exact mathematical formulation of the discretization process for Maxwell’s equations in the context of the RTL method in the following section. We also show how to perform an algebraic decoupling of the resulting transmission line equations and the transfer function method for computing the actual fields from the eigenvalues obtained. In Section 3, the same tools are applied in four more difficult cases of UOF with either, asymmetry or eccentricity of the central core. We also perform the same analysis for elliptic and rectangular cores as well as for the PCF case.

2. Mathematical Equivalence of Homogeneous Circular Cylindrical Layers to Electric Transmission Lines

The basis for the application of the previously introduced RTL method is the radial discretization of all cylindrical fibres via a separation into a succession of thin cylindrical layers, each one with its own constant refractive index n. These layers can be made to extend outside of the cladding in order to take into consideration the effect of the surrounding air (n = 1). Each thin cylindrical layer could have thickness δr proportional to each average radius r which means that given discrete steps as $\delta r=r_2-r_1$ with $r=\frac{r_1+r_2}{2}$ one has

$$\frac{r_2-r_1}{r_2+r_1}=\frac{c}{2}=>\{\begin{array}{c}\frac{1+\raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{1-\raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{r}_1=r_2 \left(out\right)\\ \frac{1-\raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{1+\raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{r}_2=r_1\left(in\right)\end{array}$$

(1)

For any such circular cylindrical layer Maxwell equations (for a constant wavelength, i.e., constant frequency “ω”) can be written in their standard form as

$$\{\begin{array}{c}\nabla X\overrightarrow{E}=-j\omega {\mu }_0\overrightarrow{H}\\ \nabla X\overrightarrow{H}=j\omega {\epsilon }_0{n}^2\left(r\right)\overrightarrow{E}\end{array}$$

(2)

Taking into consideration the relations $\omega {\mu }_0=k_0{z}_0$ and $\omega {\epsilon }_0=\frac{k_0}{z_0}$ where $k_0=\frac{\omega }{c}, z_0=120\mathsf{\pi }$ and replacing $z_0\overrightarrow{H} \mathrm{with} \overrightarrow{H}$ in order $\overrightarrow{E} \mathrm{and} \overrightarrow{H}$ to have the same units $\left(\mathrm{V}/\mathrm{m}\right)$, Maxwell equations become then

$$\{\begin{array}{c}\nabla X\overrightarrow{E}=-j{k}_0\overrightarrow{H}\\ \nabla X\overrightarrow{H}=j{k}_0{n}^2\left(r\right)\overrightarrow{E}\end{array}$$

(3)

In circular cylindrical geometry of coordinates (r, φ, z) the following set of three partial differential equations can be derived by the first vector Maxwell equation as

$$\{\begin{array}{c}\frac{1}{r}\frac{\partial E_z}{\partial \phi }-\frac{\partial E_{\phi }}{\partial z}=-j{k}_0{H}_r\\ \frac{\partial E_r}{\partial z}-\frac{\partial E_z}{\partial r}=-j{k}_0{H}_{\phi }\\ \frac{1}{r}\frac{\partial \left(r{E}_{\phi }\right)}{\partial r}-\frac{1}{r}\frac{\partial E_r}{\partial \phi }=-j{k}_0{H}_z\end{array}$$

(4)

Applying a Fourier Transform along “z” and “φ” with wave numbers “β” and “l”, where $l$ is integer (because along “φ” we have Fourier series of period 2π), the set (4) becomes:

$$\{\begin{array}{c}\frac{jl}{r}\overline{E_z}-j\beta \overline{E_{\phi }}=-j{k}_0\overline{H_r}\\ j\beta \overline{E_r}-\frac{\partial \overline{E_z}}{\partial r}=-j{k}_0\overline{H_{\phi }}\\ \frac{1}{r}\frac{\partial \left(r\overline{E_{\phi }}\right)}{\partial r}-\frac{jl}{r}\overline{E_r}=-j{k}_0\overline{H_z}\end{array}$$

(5)

In (5) we use new variables $\overline{E_r}, \overline{E_{\phi }}, \overline{E_z}, \overline{H_r}, \overline{H_{\phi }}, \overline{H_z}$ to denote the Fourier Transforms of the respective electromagnetic field components. Furthermore, replacing β and r by their reduced variables according to the following relations:

$$\{\begin{array}{l}\frac{\beta }{k_0}=>\beta \\ r{k}_0=>r\end{array}$$

Then (5) takes the form

$$\{\begin{array}{c}\frac{jl}{r}\overline{E_z}-j\beta \overline{E_{\phi }}=-j\overline{H_r}\\ j\beta \overline{E_r}-\frac{\partial \overline{E_z}}{\partial r}=-j\overline{H_{\phi }}\\ \frac{1}{r}\frac{\partial \left(r\overline{E_{\phi }}\right)}{\partial r}-\frac{jl}{r}\overline{E_r}=-j\overline{H_z}\end{array}$$

(6)

Following a similar approach, the second Maxwell vector Equation (3) can be written in the form

$$\{\begin{array}{c}\frac{jl}{r}\overline{H_z}-j\beta \overline{H_{\phi }}=j{n}^2\left(r\right)\overline{E_r}\\ j\beta \overline{H_r}-\frac{\partial \overline{H_z}}{\partial r}=j{n}^2\left(r\right)\overline{E_{\phi }}\\ \frac{1}{r}\frac{\partial \left(r\overline{H_{\phi }}\right)}{\partial r}-\frac{jl}{r}\overline{H_r}=j{n}^2\left(r\right)\overline{ E_z}\end{array}$$

(7)

Furthermore, following a cumbersome analysis, it is possible to prove that the system of Equations (6) and (7) can be transformed in a set of four differential Equation (8), relating the equivalent “voltage’ and “current” functions $V_M, I_M, V_E, I_E$ defined as follows:

$$\begin{gathered} V_M=\frac{l\overline{H_{\phi }}+\beta r\overline{H_z}}{jF} \\ {I}_M=\frac{r\overline{H_r}}{j}=\frac{\beta r\overline{{\overline{{{\rm E}}_{\phi }}}_{}}-l\overline{{{\rm E}}_z}}{j} \\ {V}_E=\frac{l\overline{{{\rm E}}_{\phi }}+\beta r\overline{E_z}}{F} \\ {I}_E=n^{2}r\overline{E_r}=l\overline{H_z}-\beta r\overline{H_{\phi }} \end{gathered}$$

where we use the notation $F=\frac{{\left(\beta r\right)}^2+l^2}{r}$

$$\{\begin{array}{c}\frac{\partial V_M}{\partial r}=-\frac{{\gamma }^2}{jF}{I}_M-jM{I}_E\\ \frac{\partial I_M}{\partial r}=-jF{V}_M\\ \frac{\partial V_E}{\partial r}=-\frac{{\gamma }^2}{j{n}^{2}F}{I}_E-jM{I}_M\\ \frac{\partial I_E}{\partial r}=-j{n}^{2}F{V}_E\end{array}$$

(8)

In (8) we introduced the total propagation factor ${\gamma }^2=\frac{l^2}{r^2}+{\beta }^2-n^2$ and the auxiliary function $M=\frac{2l\beta }{\left[{\left(\beta r\right)}^2+l^2\right]F}$.

At this point it is noticed that $V_M,I_M,V_E,I_E$ are continuous functions at the boundaries because the tangential components of electric and magnetic fields $\overrightarrow{H_{\phi }} \overrightarrow{ H_z}$ and $\overrightarrow{E_{\phi } } \overrightarrow{E_z}$ on the cylindrical surface are continuous functions passing the boundaries’ of the cylindrical layer. Using the previous relations, the Fourier Transforms of the electromagnetic field components along (r, $l$, β) can be expressed as functions of their equivalent “voltages” and “currents” functions with the auxiliary relations

$$\begin{gathered} \overline{H_r}=\raisebox{1ex}{$j{I}_M$}\!\left/ \!\raisebox{-1ex}{$r$}\right., \overline{E_r}=\frac{I_E}{n^{2}r} \\ \overline{H_{\phi }}=jl{V}_{M }/r-\frac{\beta }{F}{I}_E \\ \overline{E_{\phi }}=l{V}_E/r+j\frac{\beta }{F}{I}_{{\rm M}} \\ \overline{H_z}=\frac{l}{Fr}{I}_E+j\beta V_M \\ \overline{{{\rm E}}_z}=-j\frac{l}{Fr}{I}_M+\beta V_E \end{gathered}$$

It becomes evident by inspection that the final Equation (8) represents two coupled electric transmission lines.

2.1. Decoupling the Transmission Line Equations

The prescribed set of Equation (8) constitutes a homogeneous set of ordinary differential equations of r and considering that the all vectors [$V_M,I_M,V_E,I_E$] can be turned into exponential functions of r given by $V_M=V_M{e}^{\xi r}, I_M=I_M{e}^{\xi r}, V_E=V_E{e}^{\xi r}, I_E=I_E{e}^{\xi r}$, where $V_M,I_M,V_E,I_E$ are constants, i.e., not functions of r. Thus, the system (8) can be transformed in an algebraic set of the following four equations

$$\{\begin{array}{c}\xi V_M=-\frac{{\gamma }^2}{jF}{I}_M-jM{I}_E\\ \xi I_M=-jF{V}_M\\ \xi V_E=-\frac{{\gamma }^2}{j{n}^{2}F}{I}_E-jM{I}_M\\ \xi I_E=-j{n}^{2}F{V}_E\end{array}$$

(9)

Replacing $I_M=-\frac{jF}{\xi }{V}_M$, $I_E=-\frac{j{n}^{2}F}{\xi }{V}_E$, we obtain a set of two homogeneous equations

$$\{\begin{array}{c}\xi V_M=\frac{{\gamma }^2}{\xi }{V}_M-M\frac{n^{2}F}{\xi }{V}_E\\ \begin{array}{c}\xi V_E=\frac{{\gamma }^2}{n^2\xi }{V}_E-M\frac{F}{\xi }{V}_M\end{array}\end{array}\mathrm{or}\{\begin{array}{c}{\xi }^2{V}_M={\gamma }^2{V}_M-n^{2}MF{V}_E\\ \begin{array}{c}{\xi }^2{V}_E={\gamma }^2{V}_E-MF\end{array}\end{array}$$

This then leads to the eigenvalue equations

$$\{\begin{array}{c}\left({\xi }^2-{\gamma }^2\right)V_M+n^{2}MF{V}_E=0\\ MF{V}_{{\rm M}}+\left({\xi }^2-{\gamma }^2\right)V_{{\rm E}}=0\end{array}$$

From the standard form of the eigenvalue problem we obtain through the determinant differential equations as follows

$${\left({\xi }^2-{\gamma }^2\right)}^2-n^2{M}^2{F}^2=0,\mathrm{or}{\xi }^2={\gamma }^2\pm nMF$$

(10)

Hence the system has two eigenvalues and two mutually excluded or “normal” eigenvectors. The eigenvectors will be found by replacing ${\xi }^2$ by its value. Thus, for ${\xi }^2={\gamma }^2-nMF$, $n^{2}MF{V}_E-nMF{V}_M=0=>V_M=n{V}_E$ and the eigenvector is $V_S=V_M+n{V}_E$ While for ${\xi }^2={\gamma }^2+nMF$, $V_M=-n{V}_E$ and the eigenvector becomes $V_d=V_M-n{V}_E$. Their respective “current” eigenvectors are related as $\frac{I_M}{I_E}=\frac{V_M}{n^2{V}_E}=\frac{1}{n}$, thus $I_M=\frac{I_E}{n}$ and $I_s=I_M+\frac{I_E}{n}$, $I_d=I_M -\frac{I_E}{n}$. Since the auxiliary M function has the sign of l, the set $\left(V_s, I_s\right)$, for $l=-l$ becomes equal to the set $\left(V_d, I_d\right)$. Thus, we can consider as a unique solution for the set $\left(V_s, I_s\right)$ and the integer ‘$l$ ’varies $-\infty to+\infty$, and of course ${\xi }^2={\gamma }^2-nMF$:

$$\{\begin{array}{l}\frac{\partial V_s}{\partial r}=-\frac{{\xi }^2}{jF}{I}_s\\ \frac{\partial I_s}{\partial r}=-jF{I}_s\end{array}$$

(11)

Furthermore, $V_s, I_s$ should be continuous functions at their boundaries although n(r) varies from layer to layer. This is achieved via the adjustment $V_s=V_M+n{V}_E=2V_M$ and $I_s=I_M+\frac{I_E}{n}=2I_M$, which are continuous functions of r by definition.

$$\{\begin{array}{l}\frac{\partial V_M}{\partial r}=-\frac{{\xi }^2}{jF}{I}_M\\ \frac{\partial I_M}{\partial r}=-jF{I}_M\end{array}$$

(12)

Another option for achieving continuity is to consider the functions $V_{ss}=\frac{V_M}{n}+V_E$ and $I_{ss}=n{I}_M+I_E$. In this case, $V_{ss}=2V_E$ and $I_{ss}=2I_E$ are also continuous. Thus

$$\{\begin{array}{l}\frac{\partial V_E}{\partial r}=-\frac{{\xi }^2}{j{n}^{2}F}{I}_E\\ \frac{\partial I_E}{\partial r}=-jF{n}^2{I}_E\end{array}$$

(13)

Thus, the set of two coupled transmission lines (9) is equivalent to two independent transmission lines (12) and (13).

The two waves represented by the equations of transmission lines (12) and (13), are geometrically normal because the first is related to the magnetic field and the second to the electric field that are geometrically normal for transmitted EM waves. This property is an inherent property of EM modes in optical fibers related to birefringence phenomena. However, the β respective values, for any mode, are always found to be very close and can be considered as practically equal.

2.2. Equivalent Circuits for Cylindrical Layers, Boundary Conditions, and Birefringence

Taking into consideration the transmission line theory, it can be proved that each layer of infinitesimal thickness δr is equivalent to a T-circuit as the one shown in Figure 3

$$\{\begin{array}{c}{Z}_B=\frac{\xi }{jF}\tanh \left[\raisebox{1ex}{$\left(\xi \delta r\right)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right]\\ Z_p=\frac{\xi }{jF\sinh \left(\xi \delta r\right)}\end{array}$$

For $\xi \delta r\ll 1$ the impedances can be approximated by the equivalent relations

$$\{\begin{array}{c}{Z}_B=\frac{{\xi }^2\left(\raisebox{1ex}{$\delta r$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right) }{jF}\\ Z_p=\frac{1}{jF\delta r}\end{array}$$

(14)

If ${\xi }^2>0$, both $Z_B, Z_p$ are “capacitive” reactances, for ${\xi }^2<0$ however $Z_B$ becomes “inductive” reactance. For $\left(V_E, I_E\right)$ the approximate respective impedances of the T-circuit are given as

$$\{\begin{array}{c}{Z^{\prime }}_B=\frac{{\xi }^2\left(\raisebox{1ex}{$\delta r$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right) }{j{n}^{2}F}\\ {Z^{\prime }}_p=\frac{1}{j{n}^{2}F\delta r}\end{array}$$

(15)

As previously stated, the functions $\left(V_M, I_M\right)$ of each layer are continuous at the cylindrical boundaries of the layer, thus if we divide the fiber (including a sufficient number of air layers) in successive thin layers and replace them by their equivalent T-circuits, an overall lossless transmission line is formed with only reactive elements. For given “l”, the “β” values that lead to the resonance of the overall transmission line are the eigenvalues of the whole optical fiber.

When a transmission line is in resonance, at any arbitrary point $r_0$ of the line, the sum of reactive impedances arising from the successive T-circuits on the left and right sides of $r_0$ should be equal to zero, thus the equation giving the eigenvalues of the transmission line is the following:

$$\{{\dot{Z}}_{L.r_0}+{\dot{Z}}_{R.r_0}=0$$

(16)

Equation (16) provides the eigenvalues “β” for a given “l”, where ${\dot{Z}}_{L.r_0}, {\dot{Z}}_{R.r_0}$ are the overall reactive impedances of successive T-circuits on the left and right of $r_0$, using Equations (14) or (15). The value of $r_0$ is usually given by the core radius. For the same “$l$” the Equations (14) and (15) give usually slightly different values of ‘β’. This phenomenon is called “Birefringence”. For circular step index fibers, the birefringence is negligible; however, for elliptic fibers and fibers of any other non-circular cores, the birefringence phenomenon could be not negligible.

In order to calculate the overall reactive impedances on the left and right of $r_0$ we should find the impedances for $r\to 0$ and for $r\to \infty$. As we proceed to 0 or to $\infty$, the remaining piece of transmission line becomes “homogeneous”, i.e., its overall reactive impedance is equal to its characteristic impedance given by $Z=\frac{\xi }{jF} \left(or \frac{\xi }{j{n}^{2}F}\right).$ Then we must have

$$\begin{gathered} r\to \infty : F\to {\beta }^{2}r, MF\to 0, \xi \to \sqrt{{\beta }^2-n^2}\to Z_{r\to \infty }=0 \\ r\to 0: F\to \frac{l^2}{r}, \xi \to \frac{l}{r}\to Z_{r\to 0}=\frac{1}{j|l|} \left(or \frac{1}{j{n}^2|l|}\right) \end{gathered}$$

For l = 0 $Z_{r\to 0}=\infty$(open circuit at the center of the equivalent transmission line) It is useful to notice that there is an equivalence between our formulation and the classic formulation modes of optical fibers. In particular, for $l$ = 0, the modes (V_M,I_M) are the TM modes, while the modes (V_E,I_E) are the TE modes. For $l$ > 0, the modes (V_M,I_M) are the HE modes, while the modes (V_E,I_E) are their HE birefringence modes. For $l$ <0 the modes (V_M,I_M) are the EH modes, while the modes (V_E,I_E) are their EH birefringence modes. For any given $l$, using the resonance technique the β values of the two birefringence modes can be calculated. The Equation (16) is given as a MATLAB function in Appendix A.

Let us consider for example a step-index fiber of n₁ = 1.54, n₂ = 1.47 the V_M, V_E, fundamental modes for V = 3.3, can be calculated and their β/k₀ values are respectively 1.518934962534846 and 1.518340184686295, hence their birefringence is equal to 0.0004947 or 0.0391%. The β/k₀ value for the equivalent mode V_eq was also calculated and was equal to 1.518638548412019 (that is approximately equal to the mean value of the previous β/k₀ values), while the β/k₀ value calculated classically by Bessel functions is equal to 1.518642063686336. These β values are very close differing only 0.0002315%.

In the following Figure 4, the normalized birefringence of the step-index fibers for n₁ = 1.54, n₂ = 1.47, and of n₁ = 1.475, n₂ = 1.47 as functions of V are shown.

We notice that for any V, the normalized birefringence is almost proportional to Δn² = (n₁ − n₂)², thus the birefringence of step-index fibers of very small Δn is negligible. For instance, for a value of V = 2.4, and Δn = 1.54 − 1.47 = 0.07, the birefringence is found to be 0.168 × 0.0049 = 0.0008232 or ~0.055% on the average β, while for Δn = 1.475 − 1.47 = 0.005, the birefringence becomes 0.168 × 0.000025 = 0.000042 or ~0.0028% on the average β. What is remarkable is that our method is sensitive and calculates it.

2.3. Calculating “Voltages” V_M, V_E and “Currents” I_M, I_E and Resulting Fields

For any given $l$, using the resonance technique the β values of the two birefringence modes can be calculated. These β values are practically the same, thus we can consider them as equal or we can consider as the proper value of β the mean value of the two modes. Taking V_M = 1 at the center point of the fiber (r = 0), the respective value of I_M at the same point can be calculated by the respective terminal impedance. Using the matrix relations between input–output for the equivalent successive T-circuits, the values of $V_M$ and $I_M$ at the rest thin cylindrical layers can be calculated. In fact, from the general theory of the telegrapher’s equation we know that the inputs and outputs are associated via a transfer matrix as follows

$$\begin{array}{l}[V_{out}\underset{}{}{I}_{out}]=\left(\begin{array}{cc}\cosh (\xi (r)\cdot \delta r_{})& Z(r)\cdot \sinh (\xi (r)\cdot \delta r)\\ \sinh (\xi (r)\cdot \delta r)/Z(r)& \cosh (\xi (r)\cdot \delta r)\end{array}\right)\left[\begin{array}{l}{V}_{in}\\ I_{in}\end{array}\right]\\ \approx \left(\begin{array}{cc}1& Z(r)\cdot (\xi (r)\cdot \delta r)\\ (\xi (r)\cdot \delta r)/Z(r)& 1\end{array}\right)\left[\begin{array}{l}{V}_{in}\\ I_{in}\end{array}\right]\\ =\left(\begin{array}{cc}1& {\xi }^2(r)\cdot \delta r/jF(r)\\ jF(r)\cdot \delta r& 1\end{array}\right)\left[\begin{array}{l}{V}_{in}\\ I_{in}\end{array}\right]\end{array}$$

(17)

In Equation (17), the characteristic impedance should be taken as $Z\left(r\right)=\xi \left(r\right)/jF\left(r\right)$ to fit with the previous analysis. Using the relations $n{V}_E=V_M$ and $n{I}_M=I_E$ the respective values of their birefringence partners can also be calculated for every thin cylindrical layer r_i. Finally, we obtain the actual fields via the relations

$$\begin{gathered} \overline{H_r}=\raisebox{1ex}{$j{I}_M$}\!\left/ \!\raisebox{-1ex}{$r$}\right., \overline{E_r}=\frac{I_E}{n^{2}r} \\ \overline{H_{\phi }}=jl{V}_{M }/r-\frac{\beta }{F}{I}_E \\ \overline{E_{\phi }}=l{V}_E/r+j\frac{\beta }{F}{I}_{{\rm M}} \\ \overline{H_z}=\frac{l}{Fr}{I}_E+j\beta V_M \\ \overline{{{\rm E}}_z}=-j\frac{l}{Fr}{I}_M+\beta V_E \end{gathered}$$

A very useful field component for optical fibers is the value of the overall electric field at any thin cylindrical layer of average radius r that can be calculated by the formula:

$${|{\overrightarrow{E}}^{}(r)|}^2=|{\overline{E}}_r^{}{|}^2+|{\overline{E}}_{\varphi }^{}{|}^2+|{\overline{E}}_z^{}{|}^2$$

After some algebra, this leads to the formula

$$|\overrightarrow{E}(r){|}^2=\left[\frac{{(\beta r)}^2+l^2}{{(nr)}^2}\right]|V_M{|}^2+\left[\frac{1}{{(\beta r)}^2+l^2}+\frac{1}{{(nr)}^2}\right]|I_M{|}^2$$

In the next section, we extend our analysis in certain UOF cases.

3. Unconventional Fibers

The refractive index n(r,φ) of the fiber with a UOF profile in general can be described as a function of both r and φ. Each cylindrical layer of an average radius r is considered to have a local value η (φ) for $r=\frac{r_1+r_2}{2}$. Again, we make use of the generic form of Maxwell equations as in Equation (2) of the previous section. Fourier Transforming the first vector equation gives the same set of equations as in (6), since there is no difference with the UOF case, however the second Maxwell vector equation due to the presence of the general n function should now be written as

$$\{\begin{array}{c}\frac{jl}{r}\overline{H_z}-j\beta \overline{H_{\phi }}=jn{\left(l\right)}^2\otimes \overline{E_r}\\ j\beta \overline{H_r}-\frac{\partial \overline{H_z}}{\partial r}=jn{\left(l\right)}^2\otimes \overline{E_{\phi }}\\ \frac{1}{r}\frac{\partial \left(r\overline{H_{\phi }}\right)}{\partial r}-\frac{jl}{r}\overline{H_r}=jn{\left(l\right)}^2\otimes \overline{ E_z}\end{array}$$

(18)

The symbol ⊗ means convolution arising by the product of two functions of the variable φ. In the following paragraphs it will be shown how to escape this mathematical difficulty for the usual unconventional optical fibers.

3.1. UOF with Non-Circular, Non-Symmetric, or Eccentric Cores

For unconventional fibers of non-circular cores there are a set of circular layers where the refractive index varies between the inner and outer core and cladding values, respectively. In any such case, the function $n{\left(\phi \right)}^2$ is a sum of a steady component $n^2$ and a periodic function of φ of period 2π thus can be written as a Fourier series $n{\left(\phi \right)}^2=n^2+{{\displaystyle \sum }}_{-\infty }^{+\infty }{N}_k\exp \left(jk\phi \right)$. Taking into consideration that the convolution of the product of an exponential function $\exp \left(jk\phi \right)$ with any function $\mathrm{A}\left(\phi \right)$ of a Fourier Transform $A\left(l\right)$ is equal to $A\left(l+k\right)$, i.e., the convolution generates “harmonics”. The function $n{\left(\phi \right)}^2$ is in a set of cylindrical thin layers, a sum of step functions alternating between the values $n_1^2 \mathrm{and}{n}_2^2$, where $n_{1 } \mathrm{and} n_2$ are refractive indexes of core and cladding. Considering that in optical fibers the refractive indices of core and cladding are very close, one effectively has that $(n_{1 }-n_2)/n_{1 }\ll 1$. As a result, any harmonic factors N_k of the function $n{\left(\phi \right)}^2$ are negligible in comparison to its steady component $n^2$ and can be omitted.

As an example, the harmonics become maximal for equal alternation steps. In this case, the first harmonic, that has the maximum value of all harmonics, is equal to $A_1=\frac{2\left(n_1^2- n_2^2\right)}{\pi },$ while the steady component $n^2$ equals $(n_1^2+n_2^2)/2$ and we may make an approximation as $A_1/n^2$ $\simeq \left(4/\pi \right)(n_{1 }-n_2)/n_{1 }\ll 1$.

Thus, for optical fibers we can always assume that $n{\left(l\right)}^2$≃$n^2.$ Then the system (18) will become equivalent to the following

$$\{\begin{array}{c}\frac{jl}{r}\overline{H_z}-j\beta \overline{H_{\phi }}=j{n}^2\overline{E_r}\\ j\beta \overline{H_r}-\frac{\partial \overline{H_z}}{\partial r}=j{n}^2\overline{E_{\phi }}\\ \frac{1}{r}\frac{\partial \left(r\overline{H_{\phi }}\right)}{\partial r}-\frac{jl}{r}\overline{H_r}=j{n}^2\overline{ E_z}\end{array}$$

(19)

We can then follow the analysis that we did with the conventional fibers, where n² is the average value of the η²(φ) of each layer along φ in the [0, 2π] interval.

3.2. Application to Elliptic Core Fibers

The method was applied in the calculation of fundamental modes of a fiber of elliptic core of a and b major and minor semi-axis, respectively, with refractive index n₁ = 1.54, and a cladding value of n₂ = 1.47 (Figure 5) for various wavelengths (defined by various V factor values $V=2\pi b\sqrt{n_1^2-n_2^2}$) and four ellipticity ratios a/b = 1.1, 1.3, 1.5, and 2.0. Results are presented in tabulated format (

Table 1. Comparison of ellipse for the fundamental modes.

a/b = 1.1	Mathieu	RTL	Differences	Relative Diff. (0/00)
V	b11, Nο	b11, Νο
1.5	1.487454917000000	1.48753837672558	−0.00008345972558	0.056109079
1.7	1.493245000000000	1.49345098615188	−0.00020598615188	0.137945315
1.9	1.498457700000000	1.49872279558264	−0.00026509558264	0.17691229
2.1	1.503029750000000	1.50331246585078	−0.00028271585078	0.188097309
2.3	1.506994250000000	1.50727053366314	−0.00027628366314	0.183334252
2.5	1.510418500000000	1.51067588097959	−0.00025738097959	0.170403752
2.7	1.513376830000000	1.51360949624267	−0.00023266624267	0.153739794
2.9	1.515936970000000	1.51614477492828	−0.00020780492828	0.13708019
3.1	1.518160870000000	1.51834484072766	−0.00018397072766	0.121179996
3.3	1.520101500000000	1.52026268667918	−0.000161186679180	0.106036787

,

Table 2. Comparison of ellipse for the fundamental modes.

a/b = 1.3	Mathieu	RTL
V	b11, Nο	b11, Νο	Differences	Relative Diff. (0/00)
1.5	1.491188512000000	1.491027765079550	0.000160746920450	0.107797853
1.7	1.497028990000000	1.496897637511800	0.000131352488200	0.087742114
1.9	1.502119714000000	1.501986592174880	0.000133121825120	0.088622647
2.1	1.506471523000000	1.506335918376860	0.000135604623140	0.090014727
2.3	1.510205927000000	1.510039049600890	0.000166877399110	0.110499764
2.5	1.513423500000000	1.513196002523000	0.000227497477000	0.150319773
2.7	1.516170122000000	1.515897403413190	0.000272718586810	0.179873342
2.9	1.518513480000000	1.518220317018870	0.000293162981130	0.193059189
3.1	1.520539300000000	1.520228501721160	0.000310798278840	0.204400030
3.3	1.522298190000000	1.521974104375010	0.000324085624990	0.212892341

,

Table 3. Comparison of ellipse for the fundamental modes.

a/b = 1.5	Mathieu	RTL	Differences	Relative Diff. (0/00)
V	b11, Nο	b11, Νο
1.5	1.494250610000000	1.493636836360350	0.000613773639650	0.410756827
1.7	1.499922346000000	1.499340905733040	0.000581440266960	0.387646913
1.9	1.504818670000000	1.504203108321470	0.000615561678530	0.409060368
2.1	1.509039170000000	1.508315227121150	0.000723942878850	0.479737633
2.3	1.512533150000000	1.511793195339390	0.000739954660610	0.489215500
2.5	1.515493100000000	1.514745797603970	0.000747302396030	0.493108412
2.7	1.518011570000000	1.517265908951310	0.000745661048690	0.491209068
2.9	1.520166130000000	1.519429880525980	0.000736249474020	0.48432172
3.1	1.522019120000000	1.521299532774000	0.000719587226000	0.472784617
3.3	1.523620800000000	1.522924688766490	0.000696111233510	0.456879582

and

Table 4. Comparison of ellipse for the fundamental modes.

a/b = 2	Mathieu	RTL	Differences	Relative Diff. (0/00)
V	b11, Nο	b11, Νο
1.5	1.499390500000000	1.497675992184210	0.001714507815790	1.143469840
1.7	1.504727250000000	1.502885066432660	0.001842183567340	1.224264110
1.9	1.509110880000000	1.507251064452310	0.001859815547690	1.232391584
2.1	1.512712190000000	1.510915134944720	0.001797055055280	1.187968913
2.3	1.515667864000000	1.514006539095880	0.001661324904120	1.096100896
2.5	1.518085750000000	1.516632675098700	0.001453074901300	0.957175773
2.7	1.520047200000000	1.518879733445780	0.001167466554220	0.768046252
2.9	1.521631209000000	1.520816115931370	0.000815093068630	0.535670578
3.1	1.522869900000000	1.522496060692140	0.000373839307860	0.245483418
3.3	1.523799100000000	1.523962746538420	−0.000163646538420	0.107393775

) compared with previous results calculated with Mathieu functions together with differences and relative differences showing a deviation which goes as only 0.01/0.123% on the average. Note the results in Table 1, Table 2, Table 3 and Table 4 are not normalised the same way as in Figure 4, they are simply the odd modes, b₁₁ = β_o/k_o. We can see that for small ellipticities especially, the results compare very well for all V values quite well with the Mathieu Functions results. As the ellipticity becomes very large as in Table 4, the results begin to differ. However, the accuracy of the Mathieu Functions used in this case is not known precisely, so the trend is correct and the actual difference could be debated.

The steady component of the refractive index for the calculations for each radius r is defined as n₁ for r < b, n₂ for r > a, and as $(n_1{\phi }_1+n_2{\phi }_{2 })/\pi$ when b < r < b, where φ₁, φ₂ are the arcs of the circle of radius r, inside and outside the ellipse in the upper semi ellipse.

In the following Figure 6, the β diagram of the fundamental even mode of an elliptic fiber with semi axis ratio aa/bb = 2, n₁ = 1.54 and n₂ = 1.47 and variable factor defined by: $V=bb·k_0·\sqrt{n_1^2-n_2^2}$) is shown for completeness in agreement with [32].

3.3. Application to a Rectangular Core Fiber

The method was applied also in the calculation of fundamental modes of a fiber with an rectangular core, Figure 7, with aa and bb semi-sides, with refractive index n₁ = 1.54, and a cladding of refractive index n₂ = 1.47, for various wavelengths, defined by various V factor values $V=bb·k_0·\sqrt{n_1^2-n_2^2}$ and four ratios a/b = 1.1, 1.3, 1.5, and 2.0.

Birefringence results between the elliptical and rectangular waveguides are compared in tabular form in

Table 5. Comparison of birefringence for elliptical and orthogonal core fibers.

	Rectangular Core		Elliptic Core
a/b = 1.1	Fundamental Mode Values	Birefringence (TR)	Birefringence (TR)
V
1.5	1.492539945916010	0.000346704531850	0.000303585935370
1.7	1.498357447022790	0.000513384576600	0.000575417316090
1.9	1.503350570368440	0.000576764589930	0.000712035787010
2.1	1.507590076049970	0.000580793732230	0.000757892695190
2.3	1.511183551577450	0.000554099561380	0.000750121835720
2.5	1.514237341469060	0.000513428535340	0.000713597265680
2.7	1.516844499544290	0.000468033918970	0.000663528284040
2.9	1.519082616576080	0.000422786387590	0.000608733773240
3.1	1.521015087494710	0.000380083447450	0.000554096605350
3.3	1.522693317372840	0.000340959781570	0.000502161305990

,

Table 6. Comparison of birefringence for elliptical and orthogonal core fibers.

	Rectangular Core		Elliptic Core
a/b = 1.3	Fundamental Mode Values	Birefringence (TR)	Birefringence (TR)
V
1.5	1.495847758922450	0.000344065006270	0.000325722889410
1.7	1.501459058153440	0.000451639480470	0.000520210503450
1.9	1.506179273220210	0.000480411129310	0.000601842058430
2.1	1.510134432440230	0.000468166987430	0.000615857816190
2.3	1.513456882268610	0.000436878544030	0.000593498171210
2.5	1.516262818410440	0.000398422704380	0.000553682765940
2.7	1.518648025192190	0.000358904227450	0.000507155757960
2.9	1.520689537973830	0.000321272665480	0.000459732120240
3.1	1.522448794667530	0.000286796641510	0.000414381305640
3.3	1.523974761485380	0.000255871394240	0.000372472456950

,

Table 7. Comparison of birefringence for elliptical and orthogonal core fibers.

	Rectangular Core		Elliptic Core
a/b = 1.5	Fundamental Mode Values	Birefringence	Birefringence (TR)
V
1.5	1.498120450723510	0.000278361237480	0.000295449461120
1.7	1.503463400911670	0.000357573590440	0.000442306154260
1.9	1.507908503417930	0.000376605061620	0.000497613358440
2.1	1.511610102671980	0.000365648210230	0.000501375052490
2.3	1.514709614765410	0.000341283334480	0.000478771792790
2.5	1.517324041962060	0.000312143747180	0.000444266137270
2.7	1.519546684126380	0.000282536889750	0.000405759261530
2.9	1.521451022056960	0.000254483058070	0.000367378187930
3.1	1.523094834995030	0.000228815613970	0.000331142666700
3.3	1.524523718254910	0.000205765806060	0.000297919827130

and

Table 8. Comparison of birefringence for elliptical and orthogonal core fibers.

	Rectangular Core		Elliptic Core
a/b = 2	Fundamental Mode Values	Birefringence (TR)	Birefringence (TR)
V
1.5	1.501154370588750	0.000115100455570	0.00017851280204
1.7	1.505905455821760	0.000185956233950	0.00028097261437
1.9	1.509839664490200	0.000217496051150	0.00032331073887
2.1	1.513123501262850	0.000227441962700	0.00033221109868
2.3	1.515890438928040	0.000225620967080	0.00033221109868
2.5	1.518243818452060	0.000217498908240	0.00030589412246
2.7	1.520263082984900	0.000206163759050	0.00028459934024
2.9	1.522009559228890	0.000193389032200	0.00026227629912
3.1	1.523531014266760	0.000180203858540	0.00024038749997
3.3	1.524865056067340	0.000167204354840	0.00021967555642

with equal semi axis. The constant component of the refractive index for calculating with each radius r is defined as

$$n=\{\begin{array}{c}{n}_1,r<b\\ \frac{1}{\pi }\left(n_1{\phi }_1+n_2{\phi }_2\right)\\ n_2,r>b\end{array}$$

where φ₁, φ₂ being the arcs of the circle of radius r, inside and outside the rectangle in the upper half plane.

For the results in the Table 5, Table 6, Table 7 and Table 8 we have defined birefringence as simply the difference in the orthogonal fundamental modes normalized to k₀, for the rectangular and elliptical waveguides. We can see the elliptical waveguide has greater birefringence in all cases of ellipticity and for all wavelengths (represented by the V value).

To make a direct comparison, the average refractive indexes as functions of r, for an elliptic core fiber and for a rectangular core fiber of the same aa and bb and aa/bb = 2, are shown in Figure 8. In Appendix B, the MATLAB codes for the average refractive indexes of an elliptic and a rectangular core fiber are given. In Figure 8 we observe that the average refractive index of the elliptical waveguide is dropping off faster than the rectangular and this facilitates the explanation of the stronger birefringence for this case.

3.4. The PCF Case

Ιn the case of a holey fiber, we separate the whole fiber circular cross-section into a set of thin cylindrical layers, variable η along φ extending beyond the cladding to take into account the surrounding air with η = 1. Each layer’s thickness is δr = r₁ − r₂. We can then approximate n(r,φ) as n(φ) for the average r <r> = r + δr/2. The refractive index can be written as a Fourier series, i.e., as $n{\left(\phi \right)}^2={\langle n\rangle }^2+{{\displaystyle \sum }}_{-\infty }^{+\infty }{N}_k\exp \left(jl\phi \right)$. Taking into account the properties of the Fourier Transform we see that $FT(\exp (jl\varphi )\cdot f(\phi ))=f(l+l^{\prime })$ so that the expressions in the second terms of Equation (18) spread around a spectrum of harmonics. This is also to be understood as a result of successive scatterings from the bored air holes. We can now use the natural geometry of the usual hexagonal lattice to see that for each set of holes we can have either 6k harmonics. For the fundamental harmonic of $l^{\prime }$ = 1, the derived harmonics passing through a layer of 6k holes should be 6k + 1. Thus, for the fundamental wave crossing the successive layers it “sees” a different set of periodic rectangle functions that will be shown rigorously to contribute a different number of harmonics (7, 13, 19, …).

For a common harmonic to pass through, one must then take an integer product which leads to higher and higher harmonics, thus cutting out the entire spectrum apart from the last highest frequency. We conclude that for holey optical fibers, the approximation for any of its cylindrical thin layers, $\eta {(r,l)}^2\approx {\eta }^2(r)={\eta }^2$ suffices for further analysis of the resulting equations. Thus, the original system (18) becomes

$$\{\begin{array}{c}\frac{jl}{r}\overline{H_z}-j\beta \overline{H_{\phi }}=j{n}^2\overline{E_r}\\ j\beta \overline{H_r}-\frac{\partial \overline{H_z}}{\partial r}=j{n}^2\overline{E_{\phi }}\\ \frac{1}{r}\frac{\partial \left(r\overline{H_{\phi }}\right)}{\partial r}-\frac{jl}{r}\overline{H_r}=j{n}^2\overline{ E_z}\end{array}$$

(20)

For the usual hexagonal pattern of holes, we may utilize elementary analytical geometry to derive the two separate regions where the refractive index alternates between the air refractive index n = 1 value and the higher value of the crystal material. We assume that along each separate layer a large circle corresponding to each cylindrical shell of radius r from the center of the fiber to the center of a smaller hole of radius r << r₀ is cut while moving clockwise along the large circle.

Prescribing a set of circles of successive radii r for each of which we can find the air holes (in 1/6 angle of the PCF) which are cut by the particular radius each time. Each arc is computed inside its respective air hole and the total sum of them divided by π/3 expresses the average squared refractive index. As a matter of fact, the square of the refractive index in this sum is equal to one, while the refractive index in the rest arc is the square of the silica refractive index. Hence the average refractive index can be easily calculated along r. In Figure 9 and Figure 10, we show the average refractive index and the electric field of a hexagonal PCF, of n = 1.46 with a twelve layers lattice, as functions of the reduced distance from the center of the fiber for the fundamental mode. Figure 8 was generated by a MATLAB code for air-hole diameter equal to 0.8 of the air-hole distance and the air hole diameter was 3.14 times the transmitted wave length (V = 1.607847). We also notice the parametrization used as $=\left(\Lambda -d/2\right)\times 2\times \pi \sqrt{n_1^2-n_2^2}$), n₁ for the silica refractive index, n₂ ≃ minimum refractive index = 1.123, Λ for the reduced air hole distance, d for the reduced air hole diameter, and $\Lambda -d/2$ = reduced inner core of PCF.

For completeness we also show Figure 11, which is the β/V diagram of the even fundamental for a PCF. This shows how easily we can produce some very useful results for a series of unconventional fibers using the transmission line method we have developed.

In Appendix B, the MATLAB code for the average refractive index of a hexagonal photonic crystal fiber is also given.

In Appendix B, the MATLAB code for the average refractive index of a hexagonal photonic crystal fiber is also given.

4. Conclusions

The presented resonance technique can be used for the study of unconventional fibers, i.e., fibers with cores of any shape, as long as the difference between core and cladding refractive indices is sufficiently small which holds true for almost all the monomode and holey fibers. The unconventional case is proven reducible to the same technique of conventional fibers, where for each mode order $l$ we can approximate by a set of two, independent and non-homogeneous, resonant transmission lines (RTLs), each one representing one mode of the birefringence.

The simulation of unconventional fibers with RTLs gives a new, simple, and effective method for computing the eigenvalues of the RTLs representing the various modes of the holey fibers. Furthermore, for each eigenvalue, the average values of E.M. fields for every thin cylindrical layer of radius r of the unconventional fiber is directly computable from the relevant eigenfunctions of the RTLs.