Novel Model for the Angle and Skewness Dependent Transmission Behavior of Step-Index Polymer Optical Fiber

Abstract

Step-index polymer optical fibers (SI-POFs) are deployed in both sensing and data transmission systems. The optical transmission behavior of these fibers is complex and affected by intrinsic influences like modal dispersion, scattering and attenuation as well as extrinsic influences like the launching condition and the angular sensitivity of the receiver. Since a proper modeling of the transmission behavior is important in order to evaluate the suitability of the fiber for a specific application, we present a novel model for step-index multi-mode fibers (SI-MMFs) which considers all the previously mentioned impacts. Furthermore, the model differentiates scattering and attenuation for propagating rays not only by their propagating angle ${\theta }_z$ but also by the skewness ${\theta }_{\varphi }$. It is therefore possible to distinguish between guided, tunneling and refracted modes. The model uses scatter and attenuation data from previously published measurements of an SI-POF and computes the impulse response of the transmission system which is transferred to the frequency domain to derive the amplitude and phase response. A possible application for SI-POF is the length or strain measurement of the fiber by measuring the phase of a harmonically modulated signal. These sensors rely on a linear relation between the length of the fiber and the phase of the modulated signal. We demonstrate the application of the model by simulating the length measurement error that occurs for these sensors by obtaining the phase response for the corresponding optical transmission system. Furthermore, we will demonstrate the flexibility of the model by varying several influences including the excitation of different mode categories and evaluate the impact on the measurement error. Finally, we compare the simulated length error derived from the model to real data obtained from a cutback measurement. An implementation of the model, which was used for all simulations in this paper, is publicly available.

1. Introduction

Multi-mode fibers (MMFs) like polymer optical fibers (POFs) differ from single-mode fibers (SMFs) and few-mode fibers by their potential modeling approaches. Due to the limited number of modes, the latter ones can be modeled by solving the wave equation for each possible mode individually. However, this approach is unfeasible for MMFs because of their large number of modes. For example, the typical step-index polymer optical fiber (SI-POF) can guide several million modes. Common modeling approaches for MMFs include ray tracing simulations and the power flow equation (PFE) [1,2,3]. Both approaches are based on optical powers, which makes it impossible to consider interference between different modes, which occurs due to the periodicity of the electro-magnetic field. While there are approaches based on the propagation of the elctro-magnetic field in MMFs [4,5], working with optical powers is sufficient for POF since it is known not to maintain the polarization for fibers longer than ≈50 cm [6]. The PFE has been used to derive the transfer function of MMFs depending on several parameters [7]. Furthermore, the necessary coefficients describing the mode conversion and attenuation have been measured for SI-POF as well [8]. However, the separate measurement of attenuation and mode conversion caused by scattering is not without problems. The measurement of one effect is always influenced by the other since none of them can be avoided completely. This is especially true for POF since it faces significantly more scattering than a glass optical fiber (GOF).

An additional drawback of both ray tracing and the PFE is the neglect of a ray’s skewness, which makes it impossible to distinguish between tunneling and refracted modes as proposed by Snyder and Love [9]. Figure 1a shows two rays with the same propagation angle ${\theta }_z$. The red ray is a meridional ray since it crosses the optical axis. It can be seen that the trajectory of the ray changes significantly when its origin is shifted towards the edge of the fiber (black ray), even though ${\theta }_z$ remains constant. The resulting ray is called a skew ray. The difference between both rays can be described by the skewness ${\theta }_{\varphi }$, which is the angle between the projection of a ray’s trajectory on the fiber’s front surface and the tangent at the point of reflection and is depicted in Figure 1b. Closely related to ${\theta }_z$ and ${\theta }_{\varphi }$ is a third angle $\alpha$ (Figure 1a), which is the angle between the ray and the normal to the surface at the point of reflection. $\alpha$ can be expressed depending on ${\theta }_z$ and ${\theta }_{\varphi }$:

$$cos\left(\alpha \right)=sin\left({\theta }_z\right)·sin\left({\theta }_{\varphi }\right).$$

(1)

Due to the relation between the three angles, it is sufficient to consider only two of them to distinguish the different mode categories. For the further discussion, we choose ${\theta }_z$ and $\alpha$. Figure 2 shows the possible mode categories with respect to the chosen angles. It is well known that a mode can only be guided if ${\theta }_z$ is smaller than the maximum angle ${\theta }_c$, for which total internal reflection occurs:

$${\theta }_c=arccos\left(\frac{n_{clad}}{n_{core}}\right).$$

(2)

$n_{core}$ is the refractive index of the core of the fiber and $n_{clad}$ is the refractive index of the cladding. If ${\theta }_z$ exceeds ${\theta }_c$, the mode category depends on $\alpha$. Similar to ${\theta }_c$, a minimum angle ${\alpha }_{min}$ is defined:

$${\alpha }_{min}=arcsin\left(\frac{n_{clad}}{n_{core}}\right)=90{}^{\circ }-{\theta }_c.$$

(3)

If $\alpha$ exceeds ${\alpha }_{min}$, the mode is a tunneling mode and suffers therefore from a higher attenuation than a guided mode. If $\alpha$ is smaller than ${\alpha }_{min}$, the mode is a refracted mode. The attenuation of refracted modes is even higher than for tunneling modes, which is why they only contribute to the transmission behavior of very short fibers.

The number of reflections at the core-cladding interface, which a ray faces while propagating along the fiber, depends on the skewness as can be seen from Figure 1a. As a consequence, the scattering experienced by a mode depends on the skewness as well. The transmission behavior of a step-index multi-mode fiber (SI-MMF) therefore does not only depend on the angular power distribution of the light source over ${\theta }_z$, but also on the power distribution over ${\theta }_{\varphi }$. We propose a new fiber model for SI-MMF that is capable of considering the dependency of both attenuation and scattering on both angles. This allows us to investigate the influence of the angular power distribution and the different mode categories on specific applications by evaluating the phase and amplitude response obtained from the impulse response of the simulated optical transmission system.

We have previously proposed an optical strain sensor based on an SI-POF [10,11]. An intensity-modulated LED is used to excite the fiber with a fixed modulation frequency $f_{mod}$ and the strain of the fiber can be monitored by measuring the modulation phase at the end of the fiber as shown in Figure 3.

When the length of the fiber is changed by $\Delta L$, the modulation phase experiences a phase shift $\Delta \varphi$. In an ideal fiber, which is free from dispersion, scattering and other parasitic influences, the modulation phase at a specific position $\varphi \left(z\right)$ has a linear dependency on the position z. As a result, the shift of the modulation phase $\Delta \varphi$ caused by the elongation of an ideal fiber is proportional to the length change:

$$\Delta {\varphi }_{ideal}\propto \Delta L.$$

(4)

However, in a real fiber, the development of the modulation phase along the fiber is not linear to the position z and therefore not linear to the length change $\Delta L$ as well. A measured phase change $\Delta \varphi$ can correspond to a different length change $\Delta L$ at a specific position z than it does at another. The modulation phase at the end of the fiber is exposed to additional influences when the fiber is actually strained as we have recently investigated [12]. However, in order to avoid a measurement error, it is crucial to know how the real development of the modulation phase along the fiber differs from the ideal modulation phase.

We will demonstrate the application of the model by simulating the nonlinear development of the modulation phase along the fiber. For this purpose, we simulate the impulse responses at discrete positions along the sensing fiber. The impulse responses are then transferred to the frequency domain via a Fourier transform. The argument of the obtained frequency responses are the phase responses, which allows us to predict the development of the modulation phase along the sensing fiber at arbitrary modulation frequencies. Since the sensor is based on an SI-POF, the model employs scatter data we have obtained from a previously performed scatter and attenuation measurement of an SI-POF described in [13]. The model itself is independent from the actual fiber type and can model any kind of SI-MMF as long as the required scatter and attenuation data for the specific fiber type are available.

In order to study the the impact on the development of the modulation phase, we will vary the following simulation parameters:

• Launching condition,
• Mode categories,
• Angular sensitiviy of the receiver,
• Modulation frequency.

Finally, we will compare the simulated development of the modulation phase along the fiber to real values obtained from a cutback measurement.

2. Materials and Methods

2.1. Fiber Model

The model computes the impulse response of an optical transmission system including the light source, the fiber and the optical receiver. For the evaluation of the suitability of the transmission system for a specific application, the impulse response has to be transferred to the frequency domain by a Fourier transform. A brief introduction to the required system theoretical principles can be found in Appendix A, which also includes the definition of an optical impulse response in order to consider the dimensions of the involved physical quantities properly.

2.1.1. General Description of the Model

The proposed fiber model is based on a description of the power distribution at discrete, equidistant positions z along the fiber depending on ${\theta }_z$, ${\theta }_{\varphi }$ and the time t. Each position is represented by a matrix M$\left[n\right]$, with n being the position index. A matrix consists of a two-dimensional array with $M_z\times M_p$ matrix cells MC$[m_z,m_p]$ with each cell being responsible for a specific combination of a ${\theta }_z$ range (index $m_z$) and a ${\theta }_{\varphi }$ range (index $m_p$). As indicated in Figure 4, each matrix cell contains two elements:

• Cell power (CP),
• Cell scatter matrix (CSM).

The cell power (CP) describes the power over time with which a cell is being excited and can occur in two different shapes. The first one is an Dirac-shaped power pulse as described by Equation (A3) and is used to excite the cells of the first matrix ($n=0$). The second one is called cell impulse response (CIR) and describes a discrete power over time with the time index k. It is used for the excitation of all subsequent matrices ($n\ne 0$). The cell scatter matrix describes how much of the power of the current cell is transferred into each cell of the next matrix M$[n+1]$. The model does not distinguish between attenuation and scattering, since both are included jointly in the scatter matrices. While the cell scatter matrices are unique to each cell of a matrix, the total set of CSMs are identical for each matrix and therefore independent from n. As a result, there is only one set of cell scatter matrices that is referenced by each matrix.

The model requires two different types of data as input:

• Distribution of the initial excitation,
• Cell scatter matrices.

The initial power distribution describes the power distribution over ${\theta }_z$ and ${\theta }_{\varphi }$ at $n=0$ and $t=0$ and depends on the angular and positional power distribution of the light source. We have derived the mathematical equations to describe the power distribution of a light source over ${\theta }_z$ and ${\theta }_{\varphi }$ previously [14]. Furthermore, the angular conversion and Fresnel losses, which occur when the light enters the fiber, have to be considered. Both were covered in [14] as well. The discrete power $P_0$ $[m_z,m_p]$, exciting each cell, can be derived from the continuous power distribution of the total excitation $\frac{{\partial }^{2}P}{\partial {\theta }_z\partial {\theta }_{\varphi }}$:

$$P_0[m_z,m_p]={\int }_{{\theta }_{z,min}\left[m_z\right]}^{{\theta }_{z,max}\left[m_z\right]}{\int }_{{\theta }_{\varphi ,min}\left[m_p\right]}^{{\theta }_{\varphi ,max}\left[m_p\right]}\frac{{\partial }^{2}P}{\partial {\theta }_z\partial {\theta }_{\varphi }}\mathrm{d}{\theta }_z\mathrm{d}{\theta }_{\varphi }$$

(5)

with ${\theta }_{z,min}\left[m_z\right]$, ${\theta }_{z,max}\left[m_z\right]$, ${\theta }_{\varphi ,min}\left[m_p\right]$ and ${\theta }_{\varphi ,max}\left[m_p\right]$ being the angular boundaries of the excited matrix cell. As discussed in Appendix A, the transmission system has to be excited with a Dirac-shaped light pulse of a certain energy $E_0$. As a consequence, each matrix cell has to be excited with a Dirac-impulse containing an energy $E_0[m_z,m_p]$ which is proportional to $P_0[m_z,m_p]$:

$$E_0[m_z,m_p]=E_0·\frac{P_0[m_z,m_p]}{P_0}$$

(6)

$P_0$ is the total power exciting the system:

$$P_0={\int }_{{\theta }_{z,min}}^{{\theta }_{z,max}}{\int }_{{\theta }_{\varphi ,min}}^{{\theta }_{\varphi ,max}}\frac{{\partial }^{2}P}{\partial {\theta }_z\partial {\theta }_{\varphi }}\mathrm{d}{\theta }_z\mathrm{d}{\theta }_{\varphi }$$

(7)

with ${\theta }_{z,min}$, ${\theta }_{z,max}$, ${\theta }_{\varphi ,min}$ and ${\theta }_{\varphi ,max}$ being the boundaries of the total angular range considered by the model. As indicated by Equation (5), the model works with powers and not power densities. For reasons discussed in Section 2.1.3, the angular step widths of the cells are not equidistant. Working with absolute powers makes the model therefore more transparent and the power exchange between the cells during each matrix transition easier to model. The dimensions of the matrix $M_z\times M_p$ represent the number of discrete angular steps over ${\theta }_z$ and ${\theta }_{\varphi }$ and can freely be chosen. The same holds for the angular boundaries of each cell and the total angular ranges that are considered by the model. The only requirement is a valid set of scatter matrices for the chosen parameters. The cell scatter matrices are obtained from measurements of a fiber of a certain length L as described in [13]. Hence, the length $\Delta L$ between two matrices M$\left[n\right]$ and M$[n+1]$ is also defined by the cell scatter matrices.

2.1.2. Propagation Algorithm

The propagation algorithm that describes the evolution of the power distribution along the fiber is depicted in Algorithm 1. The first step is the initial excitation of the cell powers of the first matrix $\mathrm{M}\left[0\right]$. Each cell is excited with a Dirac-impulse (Equation A3) containing the energy $E_0[m_z,m_p]$:

$$P_{\delta }[m_z,m_p]\left(t\right)=E_0[m_z,m_p]·\delta \left(t\right).$$

(8)

Once this is accomplished, the main loop, which cycles through the matrices until the desired fiber length has been reached, starts. Whenever a subsequent matrix $\mathrm{M}[n+1]$ is initialized, the cell impulse response $\mathrm{CIR}\left[k\right]$ of each cell is initialized with a given number of steps K and a minimum and maximum time that can occur at the position in the fiber that is related to the current matrix. The power of each interval of $\mathrm{CIR}\left[k\right]$ is set to 0. The task of the next loop is to iterate all matrix cells of the current matrix M$\left[n\right]$. In case of the first matrix ($n=0$), the excitation power of the considered cell $P_0[m_z,m_p]$ is transferred to the next matrix according to the scatter matrix of the cell. In all other cases ($n\ne 0$), the power of each of the K time intervals of the cell impulse response has to be transferred to the next matrix $\mathrm{M}[n+1]$ according to the scatter matrix of the current cell. A detailed description of the scattering process is given in Section 2.1.3. When all required matrix transitions are completed, Fresnel losses occurring at the end of the fiber are applied to the cell impulse responses. By applying angular filters, specific receiver characteristics can be considered. The next step is to sum up the cell impulse responses of all matrix cells of the final matrix to obtain the total impulse response $P_g\left[k\right]$ that describes the whole transmission system. Finally, the obtained impulse response is normalized with respect to the exciting energy $E_0$ to obtain the normalized impulse response $g\left[k\right]$ as described in Appendix A.

Algorithm 1 Simulation of the power propagation in the fiber.

1: Load the initial distribution of the excitation. 2: for each Matrix $\mathrm{M}\left[n\right]$ do 3:     Prepare the impulse responses. 4:     for each Matrix cell $\mathrm{MC}[m_z,m_p]$ do 5:         if $0=n$ then 6:            Spread $P_0[m_z,m_p]$ to the next matrix. 7:         else 8:            for each Time interval of the cell impulse response $\mathrm{CIR}\left[k\right]$ do 9:                Spread the power to the next matrix. 10:            end for 11:         end if 12:     end for 13: end for 14: Apply the Fresnel losses. 15: Apply the receiver characteristics. 16: Sum up all CIRs. 17: Normalize $P_g\left[k\right]$.

2.1.3. Details of the Available Implementation

An implementation of the proposed model is available on GitHub [15] along with the required scatter matrices of a standard SI-POF with a diameter of 1 $\mathrm{m}$$\mathrm{m}$. Instructions on how to build and use the model are provided in the file README.md. The model applies scattering and attenuation to a propagating light wave depending on the used cell scatter matrices. As a result, some properties of the model like the angular boundaries of the matrix cells (${\theta }_{z,min}\left[m_z\right]$, ${\theta }_{z,max}\left[m_z\right]$, ${\theta }_{\varphi ,min}\left[m_p\right]$ and ${\theta }_{\varphi ,max}\left[m_p\right]$) and the distance between two matrices $\Delta L$ are defined by the CSMs. The properties of the CSMs in turn depend on the measurement setup, with which the scatter data were obtained. In this section, we will give a brief overview of the measurement setup, which was used to obtain the CSMs used in this paper and the subsequent adjustments of the model.

The scatter data we use for the model were obtained from a standard SI-POF with a diameter of 1 $\mathrm{m}$$\mathrm{m}$ and a length of $L_{Fiber}=16$ cm. As a result, the distance $\Delta L$ between two discrete propagation lengths represented by $\mathrm{M}\left[n\right]$ and $\mathrm{M}[n+1]$ is $\Delta L=16$ cm as well. A detailed description of the measurement setup can be found in [13].

Angular Bounds of the Cells

During the measurements, the fiber was excited with a collimated laser beam under different angles ${\theta }_{z0}$ between the optical axis and the collimated beam outside the fiber. Since the measurements were carried out for equidistant steps for the angle ${\theta }_{z0}$, the steps for ${\theta }_z$ inside the fiber were not equidistant. The central angle ${\theta }_{z,c}$ of a matrix cell with the index $m_z$ is therefore

$$\begin{array}{cc}\hfill {\theta }_{z,c}\left[m_z\right]& =arcsin\left(\frac{sin\left({\theta }_{z0}\left[m_z\right]\right)}{n_{core}}\right)\hfill \\ & =arcsin\left(\frac{sin\left(\frac{m_z·{\theta }_{z0,max}}{M_z-1}\right)}{n_{core}}\right)\hfill \end{array}$$

(9)

with $0\le m_z<M_z$. The measurement was performed with $M_z=86$ steps and a maximum angle ${\theta }_{z0,max}$ of 85°, resulting in a step size of 1°. However, in order to increase the resolution of the scatter matrices, the data were interpolated for a step size of 0.1°, resulting in a scatter matrix with $M_z=851$ rows. The cell bounds for ${\theta }_z$ can be expressed as

$${\theta }_{z,min}\left[m_z\right]=\frac{{\theta }_{z,c}[m_z-1]+{\theta }_{z,c}\left[m_z\right]}{2}$$

(10)

and

$${\theta }_{z,max}\left[m_z\right]=\frac{{\theta }_{z,c}\left[m_z\right]+{\theta }_{z,c}[m_z+1]}{2}$$

(11)

There are two exceptions at the boundaries with

$${\theta }_{z,min}\left[0\right]=0{}^{\circ }$$

(12)

and

$${\theta }_{z,max}[M_z-1]=arcsin\left(\frac{sin\left({\theta }_{z0,max}\right)}{n_{core}}\right)$$

(13)

During the measurements, the skewness of the excited modes was varied by partly covering the fiber’s front surface with a movable blade as shown in Figure 5a. The skewness of a ray can be influenced by a lateral shift x, which is the distance between the incident ray and center of the fiber as indicated in Figure 5b. The blade covered at least the upper half of the front surface ($x=0$) and was moved downwards in 5 $\mathsf{\mu }$ steps until the full front surface was covered ($x=r_{fiber}$), with $r_{fiber}$ being the radius of the core of the fiber. For $r_{fiber}=0.49\mathrm{m}\mathrm{m}$ follows $\frac{r_{fiber}}{\Delta x}=\frac{0.49\mathrm{m}\mathrm{m}}{5\mathsf{\mu }\mathrm{m}}=98$ measurements, since the measurement of the completely covered front surface was unnecessary for obvious reasons. As a result, the obtained cell scatter matrices have $M_p=98$ columns with $0\le m_p<M_p$. As mentioned before, $M_p$ is the number of discrete angular steps of the model for ${\theta }_{\varphi }$ and is given by the provided cell scatter matrices. During the measurements, the lower limit of the excited ${\theta }_{\varphi }$ range was independent from the position of the blade and was always $0{}^{\circ }$. The upper limit can be expressed depending on x or $m_p$ respectively:

$$\begin{array}{c}{\theta }_{\varphi ,max}\left(x\right)=\frac{\pi }{2}-arcsin\left(\frac{x}{r_{fiber}}\right)\\ {\theta }_{\varphi ,max}\left[m_p\right]=\frac{\pi }{2}-arcsin\left(\frac{m_p}{M_p}\right)\end{array}$$

(14)

The scatter data for the ${\theta }_{\varphi }$ range of a specific cell was derived from the difference of two measurements with neighboring blade positions. The maximum value of the ${\theta }_{\varphi }$ range of a cell can be obtained from Equation (14) as well. The minimum value can be expressed as:

$$\begin{array}{cc}\hfill {\theta }_{\varphi },min\left[m_p\right]& =\frac{\pi }{2}-arcsin\left(\frac{m_p+1}{M_p}\right)\hfill \end{array}$$

(15)

Cell Impulse Response

The model uses cell impulse responses ($\mathrm{CIR}\left[k\right]$), representing a power over K time intervals, to excite all matrices but the first. The minimum and maximum time covered by a particular $\mathrm{CIR}\left[k\right]$ depends on the position index n of the matrix but is independent of the particular cell of the matrix or $m_z$ and $m_p$, respectively:

$$t_{min}\left[n\right]=\frac{n·\Delta L·n_{g,core}}{c_0},$$

(16)

$$t_{max}\left[n\right]=\frac{n·\Delta L·n_{g,core}}{c_0·cos\left({\theta }_{z,max}\right)}.$$

(17)

$c_0$ is the speed of light in vacuum. The influence of chromatic dispersion on the pulse broadening in SI-POF is negligible in comparison to the influence of modal dispersion. However, chromatic dispersion does affect the group velocity of modulated signals. Therefore, we consider the refractive group index $n_{g,core}$. The temporal boundaries of ($\mathrm{CIR}\left[k\right]$) are set to the extremal values that can occur in the current matrix. As a result, a power of any ray which is scattered into a specific cell of a matrix can be recognized by ($\mathrm{CIR}\left[k\right]$), independent of its previous path.

A simplified example with three matrices is shown in Figure 6 to illustrate the application of the CIRs in the model. Since this example involves CIRs of different matrix cells and of different matrices, we need to expand the corresponding notation to ($\mathrm{CIR}[n,m_z,m_p,k]$. The example uses only three rays with the angles ${\theta }_{z0}=0^{\circ }$ (green) , 15${}^{\circ }$ (red) and 30${}^{\circ }$ (yellow). Therefore, the number of discrete steps over ${\theta }_z$ is $M_z=3$ and ${\theta }_{z0,max}=30{}^{\circ }$. We choose the number of time intervals per CIR to be $K=3$ so every angular step has its own corresponding time interval. Since the skewness of a ray has no influence on its transit time, we neglect the skewness for this example and choose $M_p=1$ and $m_p=0$, respectively.

The cells of the first matrix $\mathrm{M}\left[0\right]$ are excited with Dirac-shaped power pulses $P_{\delta }[m_z,0]\left(t\right)$. As indicated in Figure 6, the angular step width between the three rays is equidistant outside the fiber (${\theta }_{z0}$), but not inside the fiber due to the angular conversion at the entry of the fiber. Furthermore, the transit time of an ideal ray propagating under a constant ${\theta }_z$ is not proportional to ${\theta }_z$ but to $\frac{1}{cos\left({\theta }_z\right)}$. As a result, the times $t_c[n,m_z]$, at which an ideal ray arrives at a matrix $\mathrm{M}\left[n\right]$ has no linear dependency on $m_z$:

$$t_c[n,m_z]=\frac{n·\Delta L·n_{g,core}}{c_0·cos(arcsin\left(\frac{sin\left(\frac{m_z·{\theta }_{z0,max}}{M_z-1}\right)}{n_{core}}\right))}.$$

(18)

Hence, the distances between the arrival times of the different ideal rays are not equidistant at each matrix. This can be seen from the colored dashed lines, which mark the arrival times for each ray. Consequently, our implementation of a cell impulse response supports non equidistant intervals. The ideal ray propagating under the angle of ${\theta }_z=0{}^{\circ }$ reaches $\mathrm{CIR}[1,0,0,k]$ at the minimum possible time $t_c[1,0]=t_{min}\left[1\right]$ and its energy contributes to the first interval $k=0$. The thick green bar indicates the arrival of the ray in the CIR and the light green area marks the resulting power in the corresponding interval. The second ideal ray (${\theta }_z=15{}^{\circ }$) arrives $\mathrm{CIR}[1,1,0,k]$ at $t_c[1,1]$ and the third ideal ray reaches its corresponding CIR at $t_c[1,2]=t_{max}\left[1\right]$. The boundaries of the temporal steps of the CIRs are centered between the arrival times $t_c[n,m_z]$, resulting in the temporal step width $\Delta t_k[n,k]$. As long as a ray is not scattered during its transition from one matrix to another, the ideal transit time $\Delta t_c\left[m_z\right]$ can be calculated from the difference of the arrival times at two neighboring matrices under the same $m_z$:

$$\Delta t_c\left[m_z\right]=t_c[n+1,m_z]-t_c[n,m_z].$$

(19)

An ideal, unscattered ray will always arrive at a particular matrix $\mathrm{M}\left[n\right]$ at the time $t_c[n,m_z]$. During the transition from $\mathrm{M}\left[1\right]$ to $\mathrm{M}\left[2\right]$, this is the case for the green and the yellow ray. The red ray is however partly scattered during the transition. While half of its energy is maintained in the original angle and therefore arrives $\mathrm{M}\left[2\right]$ at $t_c[2,1]$, the other half is scattered to $m_z=2$. In this case, we assume that the scattering occurs in the middle of the transition and the arrival time at the next matrix is calculated from the mean transit times of the original angle and the angle to which the energy is scattered. The scattered energy contributes to the CIR of its destination angle $\mathrm{CIR}[2,2,0,k]$ and to the interval in which boundaries the time of arrival lies.

The temporal step widths of the impulse response intervals were derived in the absence of scattering. However, even when scattering is considered, the largest part of the scattered power will maintain its original angle ${\theta }_z$. As a consequence, the non equidistant spacing of the $\mathrm{CIR}\left[k\right]$ is still justified. A minor drawback of the non equidistant impulse responses is the lack of a standard procedure for the transfer to the frequency domain. However, as we show in Appendix B, the Discrete Time Fourier Transformation (DTFT) can be adjusted to work with unevenly sampled data.

Spreading the Power of a Cell Impulse Response to the Next Matrix

During the transition from a matrix $\mathrm{M}\left[n\right]$ to the next matrix $\mathrm{M}[n+1]$, the power contained by the matrix cells of $\mathrm{M}\left[n\right]$ is scattered to the matrix cells of $\mathrm{M}[n+1]$ according to the scatter matrices of the originating matrix cells. Since the power of each matrix cell depends on the time and is therefore stored in a cell impulse response, the scattering process is carried out for the power of each time interval of a cell impulse response. We have discussed the determination of the arrival time of a scattered power at the destination matrix cell and the subsequent assignment of the corresponding time interval of the destination cell impulse response in the previous section. In the following, we use a simplified example to focus on the scattering process and the involved power conversions. So far, we have denoted matrix cells as $\mathrm{MC}[m_z,m_p]$. Since we discuss matrix cells of different matrices in this section, we expand the notation to $\mathrm{MC}[n,m_z,m_p]$ to consider the position index n as well. We consider a matrix with the angular dimensions $M_z=4$ and $M_p=4$ and cell impulse responses with $K=4$ time intervals. As shown in Figure 7, we examine the scattering of the power contained by the time intervals of a single cell impulse response $\mathrm{CIR}[n,1,2,k]$. If the example was about an ideal fiber free from scattering, $\mathrm{CIR}[n,1,2,k]$ would only have a power greater than 0 for $k=1$, since it belongs to a matrix cell for $m_z=1$ and $M_z=K$. Since the CIR has some power for $k=2$ as well, the power of this time interval must have been scattered into this CIR from a CIR with a different index $m_z$ during a previous matrix transition.

The first step is the derivation of the time interval energy $E[n,1,2,k]$ for each time interval of $\mathrm{CIR}[n,1,2,k]$ which can be achieved by multiplying the power of each time interval by the corresponding temporal step width $\Delta t_k[n,k]$. The cell scatter matrix $\mathrm{CSM}[1,2,m_{z,d}]$ defines how each computed time interval energy is distributed over the destination matrix cells of the next matrix $\mathrm{M}[n+1]$ using the angular index $m_{z,d}$. As of now, the model only supports scattering over ${\theta }_z$ which is explained in Section 2.1.4. For an ideal fiber, $\mathrm{CSM}[1,2,m_{z,d}]$ would have the value 1 for $m_{z,d}=1$ and 0 for all other indices. In the example, 80% of the power maintains the original angle ($\mathrm{CSM}[1,2,1]=0.8$) with only a little power being scattered to different angles. The sum of all values of $\mathrm{CSM}[1,2,m_{z,d}]$ is <1, which means that the total power decreases during the matrix transition due to attenuation. The next step is the calculation of the time interval energies $E[n+1,m_{z,d},2,k],$ which are scattered into particular matrix cells of the next matrix from the original time interval energy $E[n,1,2,k]$. The value of the time interval energy for each destination matrix cell $\mathrm{MC}[n+1,m_{z,d},2]$ can be obtained by multiplying the current time interval energy $E[n,1,2,k]$ with the corresponding factor of the cell scatter matrix $\mathrm{CSM}[1,2,m_{z,d}]$.

The time of arrival and subsequently the index k of the destination cell impulse response is obtained as depicted in Figure 6. The power added to the destination time interval of $\mathrm{CIR}[n+1,m_z,m_p,k]$ is computed by dividing the arriving energy by the temporal step width of the destination time interval. Due to scattering, ${\theta }_z$ can change during the transition, which is indicated for $\mathrm{CIR}[n+1,3,2,k]$. In this case, power was scattered to a matrix cell with a higher ${\theta }_z$. As a consequence, the index k at which the time interval energies arrive can be higher than the original index, from which the time interval energies came from. Making it easier to keep track of the index changes, the areas corresponding to the time interval energies are colored yellow (original index $k=1$) and green (original index $k=2$). The index k of the destination impulse response can also be smaller than the original index k when a time interval energy is scattered to a lower angle ${\theta }_z$ as it is depicted for $\mathrm{CIR}[n+1,0,2,k]$.

Multithreading

Depending on the size of the matrices, the length of the fiber and the number of intervals of the impulse responses, the computation of the model can be demanding on both time and memory. Due to the parallel nature of the matrix cells, the model is ideally suited for multithreading. The current implementation supports up to $M_z$ threads. Provided a CPU with enough cores, the time necessary to compute the model can be significantly reduced. While the complete model is available as a Java implementation, the time-consuming part of the model was rewritten using the compute unified device architecture (CUDA) and integrated into the regular model. CUDA is a parallel computing platform that grants access to the graphics processing unit (GPU) of Nvidia graphics cards (Nvidia Corp., Santa Clara, CA, USA), which can offer up to a few thousand cores. While not as powerful as the cores of a regular CPU, the sheer amount of available cores makes CUDA attractive for algorithms that can be split into a large number of threads. Given that the graphics card has enough memory and depending on the regular CPU, the computation time can be further reduced.

2.1.4. Advantages and Drawbacks of the Model

The main advantage of the proposed model over existing approaches is the proper consideration of tunneling and refractive modes and their impact on the amplitude and phase response of an SI-MMF. Tunneling and refractive modes can be disabled separately in the model to investigate the impact of each mode category individually. The model can deal with light sources with arbitrary angular distributions and geometries provided their power distribution can be expressed depending on ${\theta }_z$ and ${\theta }_{\varphi }$. Equations for point sources and a two-dimensional source have already been derived [14]. Additional influences like Fresnel losses and the angular sensitivity of a given receiver can easily be considered.

A Java implementation, which was built with a focus on the modularity of the different components, exists for the model. A variety of different scenarios can be modeled and adjusted without the requirement of extensive programming skills. The implementation supports multithreading and benefits from the large number of cores of modern multicore CPUs. Working with the CUDA implementation offers the potential of an even larger performance gain, depending on the available memory and the number of cores of the GPU.

The source code of the model is open source, released under the GNU General Public License Version 3 (GPLv3) and under active development. While the features discussed in this paper are already implemented, the model is regularly expanded. Due to the licensing, everyone interested is allowed to study the code and expand it by one’s needs.

One of the drawbacks of the model is that the exact modeling of a fiber is only possible if the length is an integer multiple of $\Delta L$. However, both amplitude and phase response change rather slowly over the length compared to $\Delta L$, which allows good results for intermediary values by interpolation of the exact values. Another disadvantage is a consequence of the currently used scatter matrices, which were obtained from far field measurements. Even though the scatter matrices are available for all combinations of ${\theta }_z$ and ${\theta }_{\varphi }$, they only describe the power distribution over ${\theta }_z$ and neglect the impact on ${\theta }_{\varphi }$. Due to the lack of more detailed data, the model currently keeps the power distribution over ${\theta }_{\varphi }$ constant.

2.2. Simulation of the Phase Response

We demonstrate the application of the proposed model by evaluating the impact of all influences considered by the model on an optical strain or length sensor, which is based on the phase measurement of an optical signal whose power is harmonically modulated (see Figure 3). The sensor relies on a linear relationship between the length L of a fiber and the modulation phase at the end of the fiber [10,11,16,17,18]. The requirement for this linear relation applies to every measurement technique based on phase measurements including optical frequency domain reflectometry (OFDR). The ideal phase for which this condition holds corresponds to a fiber that allows the light to propagate under the angle of ${\theta }_z=0{}^{\circ }$ only:

$${\varphi }_{ideal}(L,{\omega }_{mod})=\frac{L·n_{core}·{\omega }_{mod}}{c_0}.$$

(20)

It can be seen that ${\varphi }_{ideal}(L,{\omega }_{mod})$ is linear to both the length of the fiber and the angular modulation frequency ${\omega }_{mod}$. However, influences like modal and chromatic dispersion, scattering, the angular power distribution of the light source and the angular receiver sensitivity influence the actual phase of the modulation at the end of the fiber. The proposed model considers all these influences when computing the impulse response $g\left[k\right]$ of a fiber of a particular length. As described in Appendix A, the impulse response can be transferred to the frequency domain leading to the transfer function $H\left({\omega }_{mod}\right)$ of the transmission system. The complex argument of the transfer function is the phase response $\varphi \left({\omega }_{mod}\right)$ which describes the modulation phase at the end of the fiber of a particular length depending on the modulation frequency. In order to evaluate the development of the modulation phase along the fiber, the model has to compute the impulse response of the fiber for each discrete position up to the maximum fiber length.

For the evaluation of the measurement error of the length sensor, which is caused by the deviation of the phase response of a particular length L from the ideal modulation phase for the same length, we use the relative error which expresses the relative deviation of a measured or simulated length $L_{meas}$ from the real length L:

$$E_r=\frac{L_{meas}}{L}-1=\frac{\varphi (L,{\omega }_{mod})}{{\varphi }_{ideal}(L,{\omega }_{mod})}-1.$$

(21)

It is important to mention that the real phase response $\varphi (L,{\omega }_{mod})$ is neither proportional to the length L nor to angular modulation frequency ${\omega }_{mod}$. The resulting relative error $E_r$ is therefore not constant over ${\omega }_{mod}$ or L. For demonstrative purposes, we compute the relative error of the phase based length measurement of an SI-POF with a length up to $L_{Fiber}=12\mathrm{m}$. We predict the impact of the launching condition, receiver characteristics and the modulation frequency $f_{mod}$ on the relative error. Furthermore, we investigate the influence of each mode category. All simulations will be performed with the scatter and attenuation matrices obtained for an SI-POF Asahi TC-1000 (Asahi Kasei Corp., Tokyo, Japan) as described in [13]. For the refractive index and the refractive group index of the core, we use the values obtained by Beadie et al. [19] for a wavelength of $650\mathrm{n}\mathrm{m}$: $n_{core}=1.49$ and $n_{g,core}=1.51$.

2.2.1. Light Source

We use two different light sources for the simulations. The first one ($L_U$) has a uniform power distribution over a one-dimensional cross section of ${\theta }_{z0}$. Since we consider the light source to be axially symmetric around the optical axis, it has a uniform power distribution over the solid angle $\Omega$, which is spanned by ${\theta }_{z0}$. If $P_T$ is the power that would be emitted in the half space $2\pi$, the power P, which is emitted in the solid angle $\Omega$, can be expressed as

$$P=\frac{\Omega }{2\pi }·P_T.$$

(22)

The power distribution over $\Omega$ is

$$\frac{\mathrm{d}P}{\mathrm{d}\Omega }=\frac{P_T}{2\pi }.$$

(23)

The solid angle can be expressed with respect to ${\theta }_{z0}$:

$$\Omega =2\pi ·(1-cos\left({\theta }_{z0}\right)).$$

(24)

This leads to

$$\frac{\mathrm{d}\Omega }{\mathrm{d}{\theta }_{z0}}=2\pi ·sin\left({\theta }_{z0}\right)$$

(25)

for the derivative. The power distribution of $L_U$ can therefore be expressed with respect to ${\theta }_{z0}$ as well:

$$\frac{\mathrm{d}P}{\mathrm{d}{\theta }_{z0}}=\frac{\mathrm{d}P}{\mathrm{d}\Omega }\frac{\mathrm{d}\Omega }{\mathrm{d}{\theta }_{z0}}=sin\left({\theta }_{z0}\right)·P_T.$$

(26)

Since we want $L_U$ to excite the full half space, we define a maximum angle of ${\theta }_{z0,max}=90{}^{\circ }$. Figure 8 shows the power distribution of $L_U$ normalized to its own maximum. While the definition of $L_U$ is not close to a realistic light source, it allows a uniform excitation of all possible modes, which is useful for evaluating their impact on the simulations. The source has a circular area with the same diameter as the core of the fiber and is placed right in front of the fiber. We assume that there is no gap between the light source and the front surface of the fiber. As a consequence, all rays enter the fiber at the same time. Nevertheless, we do consider Fresnel losses at the entry and the angular conversion according to Snell’s law. The power is equally distributed over the area.

The second light source $L_G$ has the same geometry and near field, but a non-uniform distribution over the angle and matches the angular distribution at the outputs of a Y-coupler that we used in the following experiment. The cross section of the distribution over ${\theta }_{z0}$ has a Gaussian shape and can be expressed as

$$\frac{\mathrm{d}{P}_{2D}}{\mathrm{d}{\theta }_{z0}}=exp\left(-\frac{{\theta }_{z0}^2}{2{\sigma }^2}\right)·P_T$$

(27)

with $\sigma =10.41{}^{\circ }$. The rotational symmetric power distribution around the optical axis can be obtained by integrating over the angle $\beta$ around the optical axis:

$$\frac{\mathrm{d}P}{\mathrm{d}{\theta }_{z0}}={\int }_0^{2\pi }exp\left(-\frac{{\theta }_{z0}^2}{2{\sigma }^2}\right)·P_T·\frac{1}{2\pi }sin\left({\theta }_{z0}\right)\mathrm{d}\beta =exp\left(-\frac{{\theta }_{z0}^2}{2{\sigma }^2}\right)·sin\left({\theta }_{z0}\right)·P_T.$$

(28)

The rotational symmetric power distribution of $L_G$ is shown in Figure 8 as well. Just like $L_U$, $L_G$ is normalized to its own maximum in Figure 8. If we consider the maximum angle of acceptance to be

$${\theta }_{z0,ac}=arcsin\left(\sqrt{n_{core}^2-n_{clad}^2}\right)=arcsin\left(\sqrt{{1.49}^2-{1.42}^2}\right)=26.83{}^{\circ },$$

(29)

it can clearly be seen that most of the power of $L_U$ will excite unguided modes, whereas the majority of the power of $L_G$ will end up in guided modes.

2.2.2. Receiver Characteristics

We consider two different receivers, both silicon PIN photodiodes, that will be used in the simulations and in the measurement. The first one is an S5052 [20] from Hamamatsu (Hamamatsu Photonics, Hamamatsu City, Japan) with a relatively narrow angular sensitivity. The second one is a BPW34 [21] from Vishay (Vishay Intertechnology, Inc., Malvern, PA, USA) that accepts larger angles. The angular sensitivities of both photodiodes taken from the data sheets were fitted as functions of the third order as shown in Equation (30) and are depicted in Figure 9. The corresponding coefficients are shown in

Table 1. Parameters of the angular receiver sensitivity.

Receiver	a (${\left({}^{\circ }\right)}^{-3}$)	b (${\left({}^{\circ }\right)}^{-2}$)	c (${\left({}^{\circ }\right)}^{-1}$)	d	${\theta }_{z0,rmax}$ $\left({}^{\circ }\right)$
S5052	$-1.67\times {10}^{-5}$	$-6\times {10}^{-4}$	$-3.33\times {10}^{-4}$	1	30
BPW34	0	$-1.52\times {10}^{-4}$	$2.44\times {10}^{-3}$	1	90

:

$$S_r\left({\theta }_{z0}\right)=\left\{\begin{array}{cc}a·{\theta }_{z0}^3+b·{\theta }_{z0}^2+c·{\theta }_{z0}+d,\hfill & \mathrm{for}0\le {\theta }_{z0}\le {\theta }_{z0,rmax},\hfill \\ 0,\hfill & \mathrm{for}{\theta }_{z0}>{\theta }_{z0,rmax}.\hfill \end{array}\right.$$

(30)

2.2.3. Modulation Frequency of the Transmitter

Since the phase response is a function of ${\omega }_{mod}$, we will vary the frequency with which the power of the transmitter is being modulated and evaluate the impact.

2.2.4. Mode Categories

One of the features of the proposed fiber model is the inclusion of leaky modes. For the investigation of the influence of each mode category, we will perform the simulations for three different combinations:

• Guided modes,
• Guided and tunneling modes,
• Guided, tunneling and refracted modes.

2.2.5. Influence of Reflections

Another influence on the development of the modulation phase along the fiber is optical reflection, as shown in Figure 10. When the incident light with the power $P_0$ enters the fiber, it is affected by Fresnel losses. The power which enters the fiber ($P_i$) can be expressed by the transmittance T or the reflectance R of the end surfaces of the fiber respectively:

$$P_i=P_o·T=P_0·\left(1-R\right).$$

(31)

If we assume a propagation angle of ${\theta }_z=0{}^{\circ }$, the reflectance of the end surfaces of the fiber can be expressed as

$$R={\left(\frac{n_{core}-n_0}{n_{core}+n_0}\right)}^2={\left(\frac{1.49-1.0}{1.49+1.0}\right)}^2\approx 3.9\%$$

(32)

and the transmittance as

$$T=1-R\approx 96.1\%.$$

(33)

Therefore, most of the exciting power enters the fiber and only a small part is reflected. For the same reason, the largest part of $P_i$ is transmitted at the end of the fiber:

$$P_{trans}=P_i·T.$$

(34)

However, a small part is reflected back at the end surface of the fiber. The reflected light travels back to the beginning of the fiber, is reflected again and propagates along the fiber again. The largest part of this twofold reflected light will be transmitted through the end surface, which we will call first rebound:

$$P_{reb1}=P_i·R^2·T.$$

(35)

The first rebound superimposes the transmitted light $P_{trans}$. When the fiber is excited with a light source whose power is modulated with a fixed frequency, the first rebound can interfere with the transmitted light. It is important to mention that the interference occurs between the modulated optical powers and not between the electromagnetic fields of the light. The reflected light will travel the fiber again and lead to the second rebound

$$P_{reb1}=P_i·R^4·T$$

(36)

and so on. The relation between the powers of the first rebound and the transmitted light can be expressed as:

$$\frac{P_{reb1}}{P_{trans}}=R^2\approx 1.5\times {10}^{-3}.$$

(37)

The second rebound faces four reflections:

$$\frac{P_{reb2}}{P_{trans}}=R^4\approx 2.25\times {10}^{-6}.$$

(38)

From this consideration, one can see that the power of the second rebound is so small that its impact on the transmission behavior of a fiber is negligible. The impact of the first rebound can easily be integrated in our model by calculating a second impulse response with three times the fiber length that is attenuated by the discussed reflectivities. For the estimation of the power of the impulse rebounds, we have neglected the attenuation of the fiber. The fiber model, however, does consider absorption and scattering for the impulse rebound as well as the angular dependency of the Fresnel losses.

2.3. Experimental Determination of the Phase Response

To verify the results of the model, we created a cutback measurement setup which allowed us to investigate the development of the phase of a harmonically modulated signal along a fiber. This allowed us also to evaluate the error that affects a strain or length sensor that uses an SI-POF as medium. The measurement setup is depicted in Figure 11. We used the output of a vector network analyzer (VNA) to modulate the power of an LED having a wavelength of 650 $\mathrm{n}$$\mathrm{m}$ with a fixed frequency. A Y-coupler acted as a mode scrambler and homogenized the slightly irregular far field of the LED, resulting in the angular power distribution of $L_G$ as described before at both outputs. One output was connected to the fiber under test, the other output of the Y-coupler remained unconnected. The fiber had a total length of $L_{Fiber}=12\mathrm{m}$ and was coiled with a diameter of $d_{Fiber}=20\mathrm{c}\mathrm{m}$. In previous tests, the diameter had been varied between $10\mathrm{c}\mathrm{m}$ and $40\mathrm{c}\mathrm{m}$ and had been found to have a negligible impact on the development of the modulation phase. Only the last meter of the fiber was straight with a receiver attached to it. Two receivers with the described PIN photodiodes (S5052/BPW34) were used to convert the optical power into an electrical voltage which was returned to the VNA. The fiber was being cut back in 16 $\mathrm{c}$$\mathrm{m}$ steps to match the step width of the model and the phase was measured with both receivers for each position. The fiber was only being cut along the straight section. When this was not possible anymore, one loop of the coiled section was unwinded. The fiber was being cut with a guided high quality blade but not polished. Polishing was avoided since it is a time-consuming process, which would have allowed the phase measurement to be affected by a possible temperature drift of the receivers. As a trade-off, the far field at the end of the fiber was affected by the relatively poor quality of the end surface which caused fluctuations in the measurement of the phase along the fiber.

The absolute modulation phase, which was measured by the receiver for each cutback step, did not only depend on the light propagation in the fiber. The electrical wiring, the transmitter and the receiver introduced a phase delay as well which influenced the measurement. Therefore, a reference measurement of the modulation phase was performed directly at the end of the mode scrambler. It was performed for each modulation frequency and for both receivers. Furthermore, the VNA could only determine the modulation phase within one modulation wavelength ${\lambda }_{mod}$ along the fiber. The development of the absolute phase, which was required to obtain the relative error Equation (21), was therefore reconstructed by adding up the measured phase differences between two cuts and subtracting the phase from the reference measurement.

3. Results and Discussion

3.1. Basic Validation of the Model

Before the model was used to simulate the behavior of an optical transmission system, several basic checks had been performed on the model to ensure its validity. The model offers the possibility to bypass the cell scatter matrices. In this case, the impulse response is only affected by chromatic and modal dispersion but not by attenuation and scattering. It was verified that the total energy carried by the impulse responses of the transmission system is maintained during the propagation. Furthermore, the power of each time interval k of the cell impulse response of each matrix cell is constant over the fiber length. Additionally, ray tracing simulations were set up to verify these ideal impulse responses using LightTools (version 8.4, Synopsys, Mountain View, CA, USA). The results obtained from LightTools matched the computed impulse responses of our model within the limits of accuracy given by the discrete nature of both approaches.

In order to classify the accuracy of the used scatter data, the attenuation of the Asahi TC-1000 was extrapolated from simulations performed with the model. The attenuation was found to be ≈250 dB/km, which is higher than the $150\mathrm{d}\mathrm{B}/\mathrm{k}\mathrm{m}$ stated by the manufacturer [22]. However, one has to keep in mind that the used scatter data were obtained from a $16\mathrm{c}\mathrm{m}$ piece of fiber which can introduce a relatively large measurement error. Under these conditions, we consider the obtained attenuation to be in an adequate range.

3.2. Example Simulation

For the verification of the model, we simulated the relative error of a fiber length measurement based on the phase measurement of the modulated optical power according to Equation (21). However, this required further processing of the impulse responses computed by the model. To get an impression of a simulated impulse response and the corresponding amplitude and phase response, we present the results of a example simulation in this section. Figure 12a shows an impulse response computed by the model for an SI-POF with a length of $L_{Fiber}=12\mathrm{m}$. The fiber was excited with the Gaussian power distribution $L_G$ and the receiver characteristics of a BPW34 photodiode were considered. Furthermore, only guided modes were excited. Figure 12b shows the corresponding amplitude response up to a maximum modulation frequency of 240 $\mathrm{M}$$\mathrm{Hz}$ obtained by the DTFT as described in Appendix B. The computed impulse response $P_g\left[k\right]$ was normalized by the energy with which the guided modes were excited. As it can be seen, the attenuation for the modulation frequency of 240 $\mathrm{M}$$\mathrm{Hz}$ is ≈1 dB higher than for a continuous wave (CW) signal. The computed phase response $\varphi \left({\omega }_{mod}\right)$ is not depicted since the deviation from the ideal phase response is so small that it is not directly observable. Therefore, Figure 12c shows the group delay, which is the negative derivation of the phase response with respect to the angular modulation frequency:

$${\tau }_g\left(\omega \right)=-\frac{\mathrm{d}\varphi \left({\omega }_{mod}\right)}{\mathrm{d}{\omega }_{mod}}.$$

(39)

As it can be seen, the group delay changes only a little within the depicted frequency range.

3.3. Results of the Simulations

We used the relative error as given by Equation (21) to determine the accuracy of the fiber length obtained from a phase measurement, when influences like modal dispersion are not considered. Therefore, the influence of impacts like modal and chromatic dispersion, the launching condition and the angular receiver sensitivity on the accuracy of a phase based length or strain sensors can be evaluated. Unless stated otherwise, all simulations were evaluated at a modulation frequency of $f_{mod}=240\mathrm{M}\mathrm{Hz}$. Since we were working with an LED, $240\mathrm{M}\mathrm{Hz}$ was close to the maximum achievable modulation frequency.

3.3.1. Influence of Light Source, Receiver and Mode Categories

Figure 13 shows the relative error up to a length of $L=12\mathrm{m}$ simulated with both launching conditions $L_U$ (uniform distribution) and $L_G$ (Gaussian distribution) for the photodiode BPW34, which has a large acceptance angle. The simulation was performed with guided modes only (G), guided and tunneling modes (G + T) and finally with all possible modes including refracted modes (G + T + R). Inside the observed range, the relative error is positive in all scenarios but decreasing with the length. A positive relative error means that the measured length is longer than the real length. The light source $L_U$ couples more power into higher order modes than $L_G$ and the BPW34 detects modes up to relatively high angles (see Figure 9). Since the transit time of each ray depends on ${\theta }_z$, the relative errors obtained for $L_U$ are larger than the ones for $L_G$. If we compare the results of the three different mode categories, we can see that it does make a difference if we consider only guided modes (G) or tunneling modes additionally (G + T). The difference is more pronounced for the excitation with $L_U$ as for $L_G$ since the latter one is hardly exciting leaky rays. The further consideration of refracted modes (G + T + R) shows almost no influence on the relative error since they carry only little power, no matter which light source we choose. In fact, the influence of the refracted modes is not visible at all for the excitation with $L_G$. Overall, we can see that the model predicts only small errors in the phase measurements that vary slowly with the fiber length. The optical strain or length sensor would therefore measure a strain or length that is larger than the actual strain or length. Depending on the light source, the error ranges between 2.4% and 2.8% for $L_G$ and 2.7% and 5.7% for $L_U$.

The simulations for the same conditions but for the photodiode S5052 are depicted in Figure 14. In comparison to the previous scenario, all results show a smaller relative error, which is caused by the limited angular sensitivity of this photodiode (Figure 9). The relation between the different mode categories is similar. Refracted modes play a negligible role, but this time the influence of the tunneling modes is very limited as well, since they are hardly detected by the receiver. This effect is especially pronounced for $L_G$ where all three mode categories can barely be distinguished, but also for the excitation with $L_U$, the difference between guided and unguided modes became smaller. The error ranges between 2.3% and 2.55% for $L_G$ and 2.65% and 3.55% for $L_U$.

It might be confusing to see that, for the excitation with $L_U$, the relative error simulated for guided modes can actually get slightly larger than for the simulations that consider leaky modes as well. In the given scenario, this happens at $L_{Fiber}⪆10\mathrm{m}$. Since the simulation with all possible modes carries more power in higher order modes, it leads to a broader impulse response and is therefore often expected to result in a larger measurement error. However, we have to take into account that the signal is modulated with a frequency of $f_{mod}=240\mathrm{M}\mathrm{Hz}$. The phase cannot be obtained by simply considering the mean transit time of all excited modes, which would lead to a smaller relative error for the simulation of only guided modes. Instead, it is important to consider the modulation phase of each mode whose periodicity depends on its transit time and on the modulation frequency. In other words, the phase delay, which an input signal faces, does not only depend on the mean transit time of the impulse response but also on the shape of the impulse response and the modulation frequency of the input signal.

3.3.2. Influence of the Modulation Frequency

Furthermore, we evaluated the influence of the modulation frequency on the simulations. Figure 15 shows the relative error for the excitation with $L_G$ and both photodiodes BPW34 and S5052 for the modulation frequencies $80\mathrm{M}\mathrm{Hz}$, $160\mathrm{M}\mathrm{Hz}$ and $240\mathrm{M}\mathrm{Hz}$. The modulation frequencies were chosen from a range that could be experimentally evaluated. The simulations considered all mode categories. The relative error for $f_{mod}=240\mathrm{M}\mathrm{Hz}$ was already shown in the previous diagrams. As we can see, for both photodiodes, the decrease of the relative error with the fiber length is less pronounced for lower modulation frequencies. For the BPW34, the relative error at the maximum fiber length of 12 $\mathrm{m}$ increases from $2.4\%$ at 240 $\mathrm{M}$$\mathrm{Hz}$ to $2.5\%$ at 80 $\mathrm{M}$$\mathrm{Hz}$. The effect is a little less pronounced for the S5052, where the relative error changes from $2.3\%$ (240 $\mathrm{M}$$\mathrm{Hz}$) to $2.35\%$ at 80 $\mathrm{M}$$\mathrm{Hz}$.

3.3.3. Influence of Reflections

Figure 16 shows the same scenario but considers the first rebound as described in Section 2.2.5. For the modulation frequency of 80 $\mathrm{M}$$\mathrm{Hz}$, it can clearly be seen that the reflections cause an oscillating behavior of the relative error. The absolute error oscillates with a constant amplitude over the length of the fiber. Since the relative error is related to the fiber length, its amplitude has a maximum at the beginning of the fiber and decreases with the length of the fiber. The rebound has to travel two times the fiber length before it interferes with the transmitted light; therefore, the spatial oscillation period is

$${\lambda }_{osc}=\frac{{\lambda }_{mod}}{2}=\frac{c_0}{n_{core}·f_{mod}·2}$$

(40)

with ${\lambda }_{mod}$ being the modulation wavelength. For $f_{mod}=80\mathrm{M}\mathrm{Hz},$ this leads to ${\lambda }_{osc}\approx 1.26\mathrm{m},$ which can be seen in the depiction. The oscillation period for $f_{mod}=160\mathrm{M}\mathrm{Hz}$ is $\approx 0.63\mathrm{m}$ and is harder to observe in the depiction due to the relatively large simulation step width of 16 $\mathrm{c}$$\mathrm{m}$. The oscillation at the modulation frequency of $240\mathrm{M}\mathrm{Hz}$ can hardly be recognized at all since the oscillation period has decreased to ≈0.42 m. For the same reason, the maximum relative error for $f_{mod}=160\mathrm{M}\mathrm{Hz}$ and $f_{mod}=240\mathrm{M}\mathrm{Hz}$ is less pronounced than for $f_{mod}=80\mathrm{M}\mathrm{Hz}$. Apart from the oscillations, we can see that the reflections do not influence the overall development of the relative error. For a better comparability of the different curves, we define three different errors:

• Peak error is the maximum error that occurs at the first simulated position ($0.16$ $\mathrm{m}$),
• Max. error is the maximum occuring error at the beginning of the fiber neglecting the peak,
• Min. error is the error at the length of 12 $\mathrm{m}$.

The peak error has its maximum for $f_{mod}=80\mathrm{M}\mathrm{Hz}$. The real peak value at $L=0\mathrm{m}$ is independent from the modulation frequency. However, since the first simulated position is $L=0.16\mathrm{m}$ and the oscillation period is shorter for higher frequencies, the peak error decreases for higher modulation frequencies.

3.4. Results of the Experiment

The experimental results were obtained with the measurement setup described in Section 2.3 and are shown in Figure 17. For a better comparability the depiction shows the same error range as Figure 16. Furthermore, the simulated and measured peak, maximum and minimum errors are compared in

Table 2. Measured and simulated relative error.

Receiver	Modulation Frequency	Peak Error Measured/Simulated	Max. Error Measured/Simulated	Min. Error Measured/Simulated
BPW34	80 $\mathrm{M}$$\mathrm{Hz}$	4.3%/3.45%	2.95%/2.85%	2.85%/2.5%
BPW34	160 $\mathrm{M}$$\mathrm{Hz}$	3.4%/3.25%	2.95%/2.85%	2.8%/2.45%
BPW34	240 $\mathrm{M}$$\mathrm{Hz}$	3.4%/3.0%	2.95%/2.85%	2.65%/2.4%
S5052	80 $\mathrm{M}$$\mathrm{Hz}$	2.9%/3.15%	2.4%/2.55%	2.45%/2.35%
S5052	160 $\mathrm{M}$$\mathrm{Hz}$	2.35%/2.95%	2.4%/2.55%	2.4%/2.35%
S5052	240 $\mathrm{M}$$\mathrm{Hz}$	2.35%/2.7%	2.4%/2.55%	2.35%/2.3%

. As it can be seen, the overall development of the relative error is predicted relatively well by the model. Since the data were slightly smoothed, the oscillations for 160 $\mathrm{M}$$\mathrm{Hz}$ and 240 $\mathrm{M}$$\mathrm{Hz}$ are less pronounced than for the simulations. The simulated maximum error for the BPW34 is only ≈0.1% smaller than the measured maximum error. For the S5052, the predicted maximum error is ≈0.15% larger.

The reflection induced peak errors differ a little more. For the BPW34, the simulated peak error at 80 $\mathrm{M}$$\mathrm{Hz}$ is 0.85% too small. The simulated peak errors for the other modulation frequencies are too small as well (0.15% at 160 $\mathrm{M}$$\mathrm{Hz}$ and 0.4% at 240 $\mathrm{M}$$\mathrm{Hz}$). A possible explanation is that the simulation only considers reflections at the end surfaces of the fiber. It is likely that an additional reflection occurs at the surface of the receiver, which is even stronger due to the high refractive index of silicon. The simulated peak errors for the S5052 are larger than the measured errors (0.25% at 80 $\mathrm{M}$$\mathrm{Hz}$, 0.6% at 160 $\mathrm{M}$$\mathrm{Hz}$ and 0.35% at 240 $\mathrm{M}$$\mathrm{Hz}$). This supports the just stated theory if we consider the fact that this photodiode has a convex lens attached to it which can reduce the amount of power that is reflected back into the fiber. Since the simulated peak errors are larger than the measured peak errors, one could conclude that the reflection at the unpolished fiber’s end surface might not be relevant at all.

The development of the relative error over the length of the fiber shows the expected dependency on the modulation frequency. The higher the frequency, the steeper the decrease of the error. In general, the decrease predicted by the model is more pronounced than for the measurements. Compared to the measurements, the minimum error predicted by the model is ≈0.35% smaller for the BPW34 and ≈0.05% smaller for the S5052. This behavior could be caused by the used scatter data, which were obtained from a 16 $\mathrm{c}$$\mathrm{m}$ piece of the Asahi TC-1000. Obtaining the scatter data of such a short fiber is a complicated task since the accuracy of the results depends on a variety of different influences as described in [13]. As the decrease of the relative error happens faster in the simulation than it does in the measurement, it could be concluded that the scattering of the used data is slightly too strong.

Like in the simulation, the oscillations can best be observed at the lowest modulation frequency of 80 $\mathrm{M}$$\mathrm{Hz}$ and their periods as well as the positions of the minima and maxima being in good agreement with the theoretical results. The curves for the S5052 are less smooth than for the BPW34, which can be explained by the smaller angular detection range. Due to the missing polishing of the end surface after each cut, the angular power distribution is partially randomized before hitting the receiver. Due to the different optical path lengths of the modes, the resulting phase depends on the sum of all detected modes. Since the receiver detects a random set of modes at each cut, the detected phase dithers within a certain range. The larger the angular detection range of the receiver, the smaller the variance of the set of detected modes. Another effect is also related to the small angular range of the S5052. Between the fiber length of 8 $\mathrm{m}$ and 12 $\mathrm{m}$, the relative error varies slowly. In fact, the relative error for the S5052 increases for all frequencies between $11.5$ $\mathrm{m}$ and 12 $\mathrm{m}$, which cannot be explained by the mode mixing at the end surface of the fiber. However, during the cutback measurement, the geometry of the fiber is altered whenever another loop is unwinded. This could have a small impact on the angular power distribution inside the fiber and therefore cause a slowly developing drift of the phase. Again, the BPW34 is not affected due to the large angular detection range.

4. Conclusions

We have presented a new propagation model for SI-MMF that enables the prediction of the transmission behavior and therefore allows the evaluation of the fiber’s behavior for both analog sensing and digital data transmission applications via the phase and amplitude response obtained from the computed impulse response. It considers the launching condition and the angular sensitivity of the receiver as well as the angle- and skewness dependent attenuation and scattering inside the fiber. The model allows the selective consideration of the different mode categories (guided, tunneling and refracted modes) to study their impact. An implementation of the model that takes full advantage of modern computer hardware is publicly available.

Simulations performed with the model lead to the conclusion that, while the influence of refracted modes on the phase response is negligible, the impact of tunneling modes depends strongly on the launching condition and on the angular sensitivity of the applied receiver. For realistic light sources, however, their impact is very limited as well.

Furthermore, we have demonstrated the application of the model by investigating the relative error that occurs in optical length or strain sensors that are based on the phase measurement of a harmonically modulated optical signal. For a fiber up to a length of 12 $\mathrm{m}$, the resulting length measurement error lies in the range of 2% to 3%, if the discussed influences like modal dispersion or the refractive group index are not properly considered, and changes only slowly with the length of the fiber as long as reflections are neglected. They can, however, be covered by the model as well by adding a second impulse response that is caused by the reflections. The determination of the proper reflectivities is, however, a demanding task since the reflectivity of the photodiode has to be taken into account as well.

A cutback measurement has been performed with different receivers and modulation frequencies to verify the model. The results are generally in good agreement with the simulations and the deviations can be explained. If we neglect reflections, the maximum deviation of the relative error of the measurement from the model occurs for the minimum error of the BPW34 at $80\mathrm{M}\mathrm{Hz}$, where the measured error is ≈0.35% larger than predicted. While receivers with a broader angular acceptance range show a larger overall error, they are less prone to fluctuations of the modal power distribution caused by unprepared end surfaces or the geometry of the fiber.