Frequency response in step index plastic optical fibers obtained from the generalized power flow equation

We present a method to obtain the frequency response of step index (SI) plastic optical fibers (POFs) based on the power flow equation generalized to incorporate the temporal dimension where the fibre diffusion and attenuation are functions of the propagation angle. To solve this equation we propose a fast implementation of the finite-difference method in matrix form. Our method is validated by comparing model predictions to experimental data. In addition, the model provides the space-time evolution of the angular power distribution when it is transmitted throughout the fibre which gives a detailed picture of the POFs capabilities for information transmission. Model predictions show that angular diffusion has a strong impact on temporal pulse widening with propagation. © 2009 Optical Society of America OCIS codes: (060.0060) Fiber optics and optical communications; (060.2270) Fiber characterization; (060.2300) Fiber measurements; (060.2310) Fiber optics.


Introduction
Transmission properties in multimode fibers, such as attenuation and bandwidth, show a nonlinear dependence with fibre length whose origin is generally attributed to mode coupling.For example, the concatenation factor, that describes the dependence of fibre bandwidth with length, is an empirical parameter widely used when estimating the fiber maximum span for a given data rate.Plastic optical fibers, with highly multimode transmission and strong mode coupling, show a complex relationship between bandwidth and length different from ray-theory predictions.There are several proposals in the literature to obtain bandwidth from other parameters more easily measured, such as the numerical aperture (NA) [1], but we have shown that these simple models cannot give a complete description of the changes in bandwidth with length [2].
In another work, we devised a method based on Gloge's power flow equation and on experimental far field patterns (FFPs) to obtain the angular diffusion and attenuation functions characteristic of a given fiber [3].These functions provide an account of the fiber spatial behavior and, along with the power flow equation, can be used to predict output power angular distributions at any fiber length and for any launching condition [4].This description was incomplete because the temporal dependence was not explicit in the equation and frequency response and bandwidth could not be calculated.
Here, we propose the use of Gloge's generalization of the differential power flow equation to obtain the pulse temporal spread with propagation [5].The temporal dimension is introduced into the equation which can then be solved in the frequency domain.This approach was used in [6] to obtain the impulse response solving the differential equation by the Crank-Nicholson scheme.In this paper we present a fast matrix approach of the finite-difference method to solve the power flow equation for the given angular diffusion and attenuation functions.Thus, the frequency responses for a given fibre can be obtained at a range of lengths to derive the bandwidth dependence with distance.The method is verified by comparing its estimates to experimental measurements of frequency responses for the same fibers and conditions as the FFPs used to determine the angular diffusion and attenuation functions.
In this paper, we describe our method, based on Gloge's differential equation, and the fast procedure devised to solve it.Then, we briefly describe the experimental set-up and methods to obtain the frequency responses and the FFPs versus length for the tested POFs.Third, we present the results, comparing the experimental and predicted frequency responses and the bandwidth dependence on fibre length.Afterwards, we discuss the model predictions that can be applied to understand POF behavior in real links, and finally, we summarize the conclusions.

Matricial approach proposed to solve the space-time power flow equation
We use Gloge's power flow equation to describe the evolution of the modal power distribution as it is transmitted throughout a POF where different modes are characterized by their inner propagation angle with respect to fiber axis (θ ), which can be taken as a continuous variable [3].We make no assumptions about the angular diffusion and attenuation which are described as functions of θ , d(θ ) and α(θ ) respectively.Following the procedure described in [5] to introduce the temporal dimension, the partial derivative of the optical power, P(θ , z,t), with respect to z produces the two terms on the left hand of the equation, giving: and, given that ∂t ∂ z = n c cos θ , we get the following equation: Then, we take the Fourier transform at both sides of Eq. ( 2) and use the Fourier derivation property to obtain the following simplified equation: where p(θ , z, ω) is the Fourier transform of P(θ , z,t).
To solve this differential equation we implement a finite-difference method where we use a forward difference for the first z derivative, and a first and second-order central differences for the first and second angular derivatives respectively.Thus, the power at angle θ and distance z + Δz is obtained as the linear combination of the power at the same angle and the two adjacent angles (θ + Δθ , θ − Δθ ) for a distance z as shows the following equation: Equation ( 4) can be expressed in a more compact representation in matrix form.In fact, the differential changes in the angular power distribution at each Δz step are given by a simple matrix product.Thus, given the angular power at an initial length z 1 , the power distribution at a longer length z 2 can be calculated with the following matrix equation: where p is a vector where each component k is the power at the discretized propagation angle θ = k • Δθ , and m = z 2 −z 1 Δz is an integer that can be found for any pair of lengths, z 2 > z 1 providing we choose a small Δz.
A is a diagonal matrix that describes power propagation without diffusion.Its elements are obtained from Eq. (4) as which is the first order approximation of that it is the exact solution of the Eq.(3) in the absence of diffusion.We have used this later expression to ensure that the matrix elements are always positive which solves stability problems for large α values and permits the use of greater Δz steps.Notice that A is the only frequency dependent term in Eq. (5).Iteration over the values of ω gives the complete spatial and temporal evolution of the optical power in the fiber.The matrix D is a tri-diagonal matrix which accounts for diffusion along the fiber.Its elements for k > 0 are: These matrix elements describe power diffusion through a differential length of the fiber indicating the fraction of power that flows out from a given angle and the fraction that drifts to this angle from the adjacent ones.
To derive the undetermined value at k = 0, corresponding to θ = 0, we need to consider only the last term of Eq. ( 1) which is the following: and expressing the derivatives by the central differences, the following expression is obtained, As P(0 − Δθ , z) = P(0 + Δθ , z) by the even nature of P(θ ), Eq. ( 11) is reduced to and the boundary condition at θ = 0 for a Δz step is given by The other boundary condition at the maximum k = N is given by Eq. ( 8), but the value of D N,N−1 must compensate for the absence of the term D N,N+1 such that the sum of terms is zero, resulting in The physical meaning of this condition is that there are no losses due to diffusion and it holds for all rows of matrix D.
(C) 2009 OSA 16 February 2009 / Vol.17, No. 4 / OPTICS EXPRESS 2853 Matrix (A(ω)+D) carries all space-time information concerning power propagation through the fiber and thus, gives a complete description of the fiber as a transmission system.In fact, for ω = 0, the solution of Eq. ( 1) is the radial profile of the FFP at a given length L, P(θ , z = L).
The key of the method we propose to solve Eq. ( 3) is to take advantage of the sparse nature of this matrix.Therefore, calculating multiple matrix powers is more efficient than performing the same number of iterations, particularly when using MatLab ® .Even more, it is not necessary to re-calculate these matrices when changing the initial condition to obtain the space-time output power distributions, as they only depend on the fiber diffusion and attenuation.The values of Δz and Δθ that are critical for convergence have been determined according to the required precision.In the calculations presented here we have used Δz = 0.001 m and Δθ = 0.005 rad obtaining accurate results.The execution time for these values and a typical simulation of a 150 m fibre length to obtain the power distribution at 5 m steps is 0.03 s, more than 70 times faster than the method used in [3] , based on the MatLab ® partial derivative equation solver.
Once the initial condition in vector form is multiplied by the system matrix, it is possible to obtain the frequency response at a given length z = L for each output angle, p(θ , z, ω), as well as the global frequency response by integrating the total power for all angles at a given frequency: Its inverse Fourier transform gives the pulse temporal spreading at this length, h(L,t).An important aspect revealed by the triple dependence of power with propagation angle, length and time is that to obtain the frequency response at one given length L, it is not enough to know the frequency response at any shorter length, H(L 0 < L, ω).To compute the total acquired delay at a given angle and fibre length it is necessary to know the previous path followed by the power reaching that angle, which implies to know either p(θ , L 0 < L, ω) or p(θ , L 0 < L,t).Thus, to be able to calculate the frequency response at any given length, it is necessary to know the angular power distribution right at the fibre input: P(θ , z = 0,t = 0) or p(θ , z = 0, ω), where there is no propagation acquired temporal delay.In fact, previous experimental results suggest that the input distribution has a strong impact on bandwidth changing the balance of diffusion and differential attenuation [2].

Experimental methods to obtain POF frequency response and FFP
We measured the frequency responses and FFPs versus fiber length for two samples of each of three PMMA fibers of 1mm diameter from different manufacturers: ESKA-PREMIER GH4001 (GH) from Mitsubishi, HFBR-RUS100 (HFB) from Agilent, and PGU-FB1000 (PGU) from Toray.The GH and PGU fibers have a NA of 0.5, corresponding to a 19.5º inner critical angle.The HFB fiber has a NA of 0.47 which implies an 18.5º inner critical angle.For each sample both the FFPs and frequency responses were measured under the same launching conditions, taken sequentially starting from a long fiber (175m-100m) down to 10m.We measured the frequency response at each fibre length by feeding pure sinusoidal waveforms of different frequencies to an AlGaInP laser diode (LD Sanyo DL-3147-021).The laser source emits a maximum power of 5 mW at 645 nm and has a typical divergence of 30º in the perpendicular plane, and of 7.5º in the parallel plane.Power was launched directly into the fibre using a connector, which limits the launching NA to near 0.19 (under-filled launch) but is close to the conditions achieved in real links.The receptor is based on a 1mm diameter photodiode (FDS010) with a 50Ω load resistance and whose bandwidth is only 200 MHz.The fibre output is free-space coupled to the receptor using another connector.Most of the power is captured due to the large area of the detector, avoiding spatial filtering that could modify the measured frequency responses.The receptor output is amplified using a 40 dB amplifier (Mini-Circuits ZKL-1R5) with a band-pass from 10 MHz to 1.5 GHz.A wideband Infinium DCA 86100A oscilloscope from Agilent is connected to the output of the amplifier and captures the received signal whose amplitude is directly related with the frequency response of the system at that frequency.Device control and data acquisition are performed by the computer through the GPIB.
Data is processed to make our method more robust and to extend the bandwidth measurements far above the system bandwidth (up to 1 GHz).The frequency response of a short segment of fiber (75 cm) is measured to be used as a reference to characterize the effect of the electrical components.The detailed procedure to obtain the frequency responses is described in [8].Simultaneously, the FFP images were captured and their radial profiles were extracted using the set-up and procedures described in [3,9].The angular attenuation and diffusion functions used in this paper to characterize each fibre type were obtained from these experimental FFPs in [3].

Results
In this section, we demonstrate how our method provides a suitable model for POFs that can be applied to reproduce experimental data and to predict fibre behaviour for previous findings that up to now had not been thoroughly explained.

On the shape of the frequency response of POFs
The impulse and frequency responses of plastic optical fibres have been usually modelled in the literature as Gaussian functions [10], but our previous measurements of the frequency response never exhibited a Gaussian shape [2,8].Thus, we first compare the shape of both our experimental and predicted frequency responses to Gaussian functions whose bandwidths are the same as those obtained from our measured frequency responses.Fig. 1 shows the data for one of the samples of the HFB fibre, as an example, at three different lengths: 15 meters on the left graph, 50 meters in the middle and 100 meters on the right.The red lines are Gaussian functions with the same bandwidth as the corresponding experimental frequency responses, shown as black dots.The green lines show the model predicted frequency responses obtained introducing directly the angular distribution estimated from the transmitter characteristics as the input power distribution.The diffusion and attenuation functions, d(θ ) and α(θ ) used in the model calculations were those previously estimated for the three different POFs from experimental FFPs [3].All three graphs in Fig. 1 reveal that the shape of both the experimental and predicted frequency responses is far from Gaussian.Particularly, the high frequency fall-off is much steeper for the Gaussian functions than for both the experimental and the model predicted frequency responses.This discrepancy occurs even for the frequency responses at 100 m where the measured FFPs were near the steady state distribution.On the other hand, the green lines, showing the model predictions for 50 and 100 meters, closely tailor our experimental curves both at high and low frequencies.
For the 15 m fibre shown in the leftmost graph of Fig. 1, however, model prediction and experimental data, although similar in shape, are shifted by an offset, with the simulated frequency response being higher than the measured one.These model discrepancies from our measurements were found in all fibre samples at the shorter measured lengths, from 10 meters and up to 25-30 meters.This result is not unexpected and, in fact, it agrees with our previous conclusions derived from the analysis of our measurements of FFPs for short fibres [11] .Those measurements suggested that optical power suffers a strong initial diffusion at the fibre input, which is higher than the subsequent propagation diffusion, and whose pattern strongly depends on the fibre type.Thus, we propose to introduce this strong initial diffusion into our present framework as an independent effect by means of another matrix, which we will call injection matrix J from now on, and will be different for each fibre type.Thus, to obtain an estimate of the J matrix, we propose an arrangement of our experimental radial profiles for 1.25 meter fi- bres obtained launching a He-Ne laser beam at different angles.Thus, the profile for each input angle is introduced as the corresponding matrix column to describe how the power injected at one given input angle is spread over other angles at the fibre input.The J matrices for the three fibre types are shown in Fig. 2 as colour images.Dark red shows the highest power level and dark blue, the lowest one.Horizontal and vertical axes show the input and output inner propagation angles in degrees, which are the rows and columns of J respectively.Therefore, the matrix product of the injection matrix and the transmitter angular power distribution p(z = 0 − ) gives the vector describing the angular power distribution just after entering the fibre, p(z = 0 + ): With these input power distributions, we obtain new model predictions shown as blue lines that present a better agreement to the experimental data at short lengths than either with the Gaussian or with our model without the injection matrix.At the longer presented lengths the predictions using directly the transmitted input power and those including the injection matrix practically coincide and therefore, both reproduce the experimental measurements.This finding is consistent with the fact that the shape of the initial distribution is not critical after the power has been propagating throughout the fibre for several tens of meters bearing the effects of diffusion and differential attenuation.The modularity and flexibility of our method allowed a straightforward inclusion of the injection matrix.

Bandwidth dependence on fibre length
POF bandwidth as a function of fibre length in logarithmic coordinates has been usually fitted with straight lines.The slope of these lines is known as the concatenation factor and is related to mode coupling.In the absence of mode coupling and differential attenuation, the concatenation factor is one and therefore, a slope higher or lower than one is evidence of diffusive non-linear effects [2]. Figure 3 shows bandwidth versus length for the three tested fibres in a two-axis logarithmic scale.Bandwidths obtained from experimental frequency responses are shown as symbols of different colours for each of the two samples of the same fibre type to illustrate the slight inter-sample variability.The blue solid line represents the model prediction extended to fibre lengths longer than our measurements.The three graphs show that the model captures the general tendency of the experimental data, showing an initial decrease with a slope slightly higher than one, and a subsequent shallower dependence at longer lengths.We found that the experimental data for the GH and HFB fibres can be fitted by single straight lines of slopes 1.2 Fig. 2. Image representation of the injection matrix for the GH, HFB and PGU fibers (left, center and right graphs) obtained from our previous measurements in [11].Horizontal and vertical axis are input and output angles respectively.and 1.1 respectively.The PGU bandwidths have been measured beyond the transition length and thus, its data for both samples can be better fitted with two lines intercepting at 60 meters, with a slope of 1.3 for the first line and 0.8 for the second.These two slopes are consistent to those found in our previous measurements of different samples of the PGU fibre under similar launching conditions [2].The model predictions show that the other two fibres exhibit a similar behaviour which is not evident from the experimental measurements as they do not reach the transition length.We argue that the slope higher than one is caused by our underfilled launch that, combined with diffusion, degrades performance and thus, provokes a faster decrease in bandwidth.This effect is even greater for the PGU fibre, with a higher slope for the first segment, suggesting a very intense diffusion.This finding is consistent with the values of its diffusion function d(θ ), which are higher than for the other fibre types.The shallower segment, with a slope of 0.8, shows that the strong diffusion of the PGU fibre makes it reach the equilibrium at shorter lengths than for the other fibre types.Using the model, we found the same behaviour for the other fibre types, but with the change in slope appearing at longer lengths as Fig. 3 shows.None of our measurements, however, reveal slopes below 0.8 even though we measured up to fibre lengths where the FFP was near the steady state [3].However, previous bandwidth measurements using a scrambler at the emitter end produced slopes down to 0.6 [2].We can reproduce these lower slopes using our method to obtain the bandwidths predicted introducing in the model an overfilled launch.

Insight on the space-time power distribution
To have further knowledge of fibre behaviour, the model can be used to obtain the output power distribution as a joint function of output angle and time at any fiber length by taking the inverse Fourier transform of the calculated frequency responses for each output angle.In Fig. 4, the upper leftmost graph shows the power at the output of a 150 m PGU fiber as an example to explain the image representation that we have adopted to help the visualization of the fiber behavior.Time is shown on the horizontal axis in nanoseconds and output angle in degrees on the vertical axis.In the image, each row represents the temporal pulse arriving at a given output angle.The integrated power over the output angle results in the temporal pulse spread which is shown normalized in the graph below the image (lower leftmost graph).The pulse spread is very asymmetric with power rising first very steeply to its maximum and then, slowly decreasing with a long tail extending up to 30 ns.The image columns are the radial profiles of the spatial  power distribution at fixed times.The integrated power over time gives the radial profile of the FFP represented on the upper rightmost graph as normalized power on the horizontal axis versus angle in degrees on the vertical axis.On the image, the superimposed solid blue line joins the angular positions at which the maximum power reaches the fiber end at each temporal delay.The dashed magenta line shows the delay, δ (θ ), obtained without diffusion which is given by the ray-theory inverse cosine law: where L is the fiber length and n the refraction index of PMMA.This precise knowledge of how the angular power distribution evolves with time and fibre length helps to understand the improvement of the frequency response through spatial filtering that we previously found in [12].

Discussion
We can state that, on the basis of our experimental results, Gaussian functions can be discarded to fit POF frequency responses.However, the graphs in Fig. 1 show that the shape of the measured frequency responses is well reproduced by our model predictions.When inspecting the whole set of data, there are some discontinuities in the experimental data which are not followed by the model predictions.These discontinuities arise from localized defects or strains in the fiber which, when removed by the cut-off procedure, let the remaining fiber recover its normal predictable behavior.Although these effects are usually unknown and practically impossible to detect prior the measurements, they can be modelled using the present matrix framework and included in the method to assess their influence on fibre properties.
In addition, we have shown how the initial angular power distribution is critical to predict the temporal behaviour for short and middle length fibres.Based on previous and present results, we argue that optical power suffers strong diffusion when injected into a POF that conditions its subsequent propagation.Moreover, we have estimated these initial power distributions from our measurements for very short fibres and obtained a better agreement than when using the transmitter distribution directly.Although data for short fibres are very sensitive to defects, curvatures and strains which explain their higher variability and wider deviations from the model, our method permits to reproduce the variability in the measurements by introducing these localized effects and evaluating their impact on the results.
Figure 3 shows bandwidth versus length relationship for both our experimental data and model predictions demonstrating a good model performance in predicting the bandwidth de-  pendence on length.Bandwidth-length dependence is slightly higher than linear for both the GH and the HFB fibres up to 100 meters.We argue that the slope higher than one is caused by underfilled launch combined with diffusion which makes bandwidth decrease steeper with length.The shallower segment at longer lengths indicates that the fibre has reached its transition length.None of our measurements, however, reveal slopes below 0.8 and, only in extreme non-practical conditions, the model simulations present slopes near the theoretical limit of 0.5.The image representation in Fig. 4 helps to visualize the power distribution over space and time at the fiber output and to understand the relationship between the angular power distribution and the pulse spreading through modal diffusion.The image shows that optical power that exits the fiber over a cone from 0º to 8º is concentrated over a relatively narrow time slot.Above this angular range, pulses have a wider time spread and their peaks increase with the output angle as shows the blue line in the image.At these angles, the power peaks are reached at lower times than for the cosine prediction indicating noticeable shorter delays than those that would be obtained in the absence of diffusion.In other words, diffusion improves fiber transmission capability.The horizontal time shift at the lowest angles does not affect the fibre behavior as it is an overall delay.Therefore, the power exiting the fiber at the highest angles is also that with the longest delays which suggests an easy way to improve fiber capacity by spatial filtering out of the tail at the higher angles as was shown before [12,13].As most power is confined in a range of lower angles, filtering out the power at the highest angles will imply small power loss while producing a narrower overall impulse response.

Conclusions
We propose a fast and robust method that reproduces experimental measurements of different POF parameters and that allows to extend model predictions where it is difficult or unpractical to measure some fibre propagation properties.The method provides the space-time optical power distribution with length from which angular power distribution, attenuation, bandwidth and pulse spreading can be derived.In addition, it offers a flexible tool to study the effects of using different devices, such as scramblers, tappers, etc, or the impairments occasioned by defects and imperfections as they can be modelled as matrices to be introduced in our framework.Using the information provided by the space-time power distribution fiber transmission characteristics can be enhanced using an appropriate spatial filter, as we had shown experimentally before.Therefore, a good fiber characterization can be applied to optimize fiber performance in POF links.

Fig. 1 .
Fig.1.Frequency responses for one of the samples of the HFB fiber at three lengths (15m on the left, 50m in the center, and 100m on the right).Experimental measurements are shown in black dots.Gaussian functions with the same bandwidths as the corresponding experimental curve are shown as red lines.Model predictions without the injection matrix are shown as green lines, while the blue lines show the predictions with injection matrix.

Fig. 3 .
Fig. 3. Bandwidth versus fibre length for the three fibre types: GH on the left, HFB in the centre and PGU on the right.Symbols of different colours represent data for different samples of the same fibre.The blue lines show the model predictions extended up to 750 m.

Fig. 4 .
Fig.4.The graph on the upper left is the image representation of the space-time power distribution at the output of 150 m of the PGU fiber.Below, the lower leftmost graph shows the overall pulse spread obtained as the integral of power over output angle.The power integral over time renders the radial profile of the FFP shown on the upper rightmost graph.
(ω) are obtained by sampling the angular frequency ω as required for a precise calculation of the inverse discrete Fourier transform of p(θ , z, ω) to obtain P(θ , z,t).