Statistical Modelling of Photonic Crystal Fibre Based Surface Plasmon Resonance Sensors Resonant Peak Wavelength for Tolerance Studies

We report a statistical approach to model the resonant peak wavelength (RPW) equation(s) of a photonic crystal fibre (PCF)-based surface plasmon resonance (SPR) sensors in terms of the PCF structural parameters (air-hole diameter, pitch, core diameter and gold layer thickness) at various tolerance levels. Design of experiments (statistical tool) is used to investigate the role played by the PCF structural parameters for sensing performance evaluation—RPW, across three tolerance levels (±2%, ±5% and ±10%). Pitch of the hollow-core PCF was discovered to be the major influencing parameter for the sensing performance (RPW) of the PCF-based SPR sensor while the inner metal (gold) layer thickness and core diameter are the least contributing parameters. This novel statistical method to derive the sensing performance parameter(s) of the PCF-based SPR sensors can be applied effectively and efficiently in the designing, characterisation, tolerance analysis not only at the research level, but also in optical fibre sensor fabrication industry to improve efficiency and lower cost.


Introduction
Hollow-core photonic crystal fibres (HC-PCFs) are types of photonic crystal fibres (PCFs) that offer enormous promise for sensing applications because of the direct interaction of light and the gaseous or liquid analyte within the hollow fibre core. The periodic cladding structure of these PCFs gives them a photonic band-gap that restricts light inside the hollow-core or core filled with analytes having lower refractive index (RI). A key advantage of PCFs is that optical characteristics, such as confinement loss, dispersion and non-linearity can be influenced by varying their structural parameters [1][2][3][4]. The air-holes in the cladding can also be filled with biological and chemical analytes in gaseous or liquid forms. These are the sensing samples through which several chemical, physical and biological entities, such as temperature, pressure, liquid levels and composition (compounds and micro-organisms) are being sensed, quantified, calibrated for measurements and monitoring [5]. The sensing performance of HC-PCF-based surface plasmon resonance (SPR) sensors across these samples can be described in terms of the resonant peak wavelength (RPW) shifts, confinement loss from each sensing analyte, linearity and resolution of the sensors.
Recent studies have been reported on filling of air-holes for different SPR sensing applications. A HC-PCF-based SPR sensor was proposed by Luan and Yao [6] for SPR sensing of analytes within an RI range of 1.33 to 1.5 with sensitivity of 1800 nm/RIU. Liu et al. [7] presented a new kind of SPR sensor based on silver-coated hollow fibre structure for the detection of liquids with high RI. The highest sensitivity achieved from diameter d = 1.4 µm, pitch Λ = 2 µm, gold layer thickness t g = 50 nm, selected air-hole diameter d s = d − t g = 1.35 µm, hollow-core diameter d c = 1.8 µm.
The schematic cross-section of the proposed HC-PCF-based SPR sensor is illustrated in Figure 1. The sensing performance of the proposed sensor was analysed by FEM through COMSOL Multiphysics software and the obtained results are analysed for a RI range of 1.43 to 1.49.  Figure 1 shows the air-holes are arranged in a hexagonal lattice. The entire fibre cross-section radius was considered to be 9 µm. A perfectly matched layer (PML) of thickness 2 µm was added to the sensor structure design to absorb incident radiations without producing any back reflections [4]. The fused silica was used as the background material in this HC-PCF-based SPR sensor and its wavelength-dependent RI equation was calculated by the Sellmeier equation [27] n SiO 2 (λ) = 1 + a 1 λ 2 where n SiO 2 is the RI of the fused silica, λ is the incident light wavelength and it is measured in µm. The Sellmeier coefficients are a 1 , a 2 , a 3 , b 1 , b 2 , b 3 and their values are 0.691663, 0.4079426, 0.8974794, 0.004 679 148 µm 2 , 0.013 512 06 µm 2 , and 97.934 000 25 µm 2 , respectively. Gold is used as the plasmonic metal layer and its dielectric function is derived from the following Drude model: where λ p is the plasma wavelength of 0.168 26 µm and λ c is the collision wavelength of 8.9342 µm for gold (Au) [28]. The confinement loss obtained at the RPW is derived as, where α loss is the confinement loss while (n eff ) is the imaginary part of the effective refractive index of the fundamental core mode. The analyte RI range investigated is from 1.43 to 1.49. The sensor could be applied for detection on lymphocyte (n a = 1.43 ± 0.05) [29], monocytes (n a = 1.43 ± 0.04) [29], aqueous enzymatic extracted algae oil (n a = 1.445) [30], adsorbed single-stranded DNA layers (n a = 1.46) [31], poly-acrylic acid (n a = 1.47) [32] and poly-allylamine hydrochloride (n a = 1.49) [32].

FEM Simulation Results and SPR Discussion
The sensing performance of this HC-PCF-based SPR sensor was evaluated in terms of the sensitivity, which is proportional to the wavelength change of resonance peaks for two different analytes (∆λ RPW ). The sensitivity is calculated with the formula where ∆n a is the minimum RI change of analyte. Figure 2 shows the confinement loss spectra characteristics against wavelength range from 0.95 µm to 1.6 µm for sensing analyte with different RIs varying from 1.43 to 1.49. Results show that an increase in the RI of the analyte causes the RPW to steadily shift to lower wavelengths. As a result of the RPW shifting to shorter wavelengths when the RI of the analyte is changed from 1.43 to 1.49, it can be noticed that the confinement peak loss in the core mode steadily reduces too as shown in Table 1   Results show that for analytes with RI = 1.45, the surface plasmons were excited from the gold metal layers at a wavelength of 1244 nm. This is due to the coupling between the fundamental core mode and plasmonic mode at that wavelength. During this process, maximum incident light energy is transferred from the fundamental core mode to the plasmonic mode. Figure 3 shows the electric field distributions of the fundamental core mode for three wavelengths; 1200 nm, 1244 nm (RPW), 1300 nm and plasmonic mode (or surface plasmon mode) for the phase-matching point or RPW (1244 nm). The characteristics of fundamental core mode and plasmonic mode and confinement loss at phase-matching point or RPW, are illustrated in Figure 4. The dotted line in Figure 4 indicates the point where the fundamental core mode strongly couples with the plasmonic mode. This point as shown in the diagram corresponds to the peak of the confinement loss.

Modelling of RPW Equation(s) Using DoE for HC-PCF Sensor Tolerance Studies
Factorial design and Response Surface Methodology are the two types of DoE methods used in this research work for the statistical model derivation of RPW equation(s) with respect to the structural parameters of the HC-PCF-based SPR sensor across different tolerance levels.

DoE Factorial Design
Factorial design is expressed in the form of a first order polynomial (linear) equation [33]: where Y is the response (RPW), k is the number of parameters (k = 4 PCF structural parameters that correspond to air-hole diameter (x 1 ), pitch (x 2 ), gold layer thickness (x 3 ) and core diameter (x 4 )), is the residual error. β 0 , β i are offset constant and linear term coefficient, respectively. The offset constant which is the initial constant of the linear model Equation (5) is the average value of the response (RPW). The method of least squares is used to estimate the β coefficients [33]. Defining Y as a vector of length n = 2 k made of the RPW values calculated from n runs of FEM simulations as response, Equation (5) in matrix form is written as: where D is the matrix of size n × (k + 1) constructed with the offset constant coefficient and different levels of parameters. β is a vector of length k + 1 that consists of the offset constant β 0 and linear term coefficients β i to be calculated. To solve for β, Equation (6) Table 2. The Minitab statistical software is used to perform the statistical and numerical analysis for k = 4 parameters corresponding to n = 2 k = 2 4 = 16 simulations (see Table 3).  Table 3. The high (+2%) and low (−2%) levels for each structural parameter ( Table 2) are used for the FEM simulations run (first column in Table 3) with air-hole diameter varied for every one (2 0 = 1) of them (second column in Table 3), pitch varied for every two (2 1 = 2) of them (third column in Table 3), gold layer thickness varied for every four (2 2 = 4) of them (fourth column in Table 3) and core diameter varied for every eight (2 3 = 8) of them (fifth column in Table 3).
Factorial design is used for the preparation of the matrix D in Equation (6). Using the structure parameters for the D matrix and the values for RPW (response Y vector from last column) in Table 3, Equation (6) can be written as: Using the calculated value of β for the 2-level linear model (5), we derive: 4 (8) where x 1 , x 2 , x 3 and x 4 represent the 4 HC-PCF-based SPR sensor parameters under consideration: air-hole diameter, pitch, gold layer thickness and core diameter, respectively. Residual analysis was carried out to confirm the effectiveness of the derived linear model. In Figure 5a,b, the residual errors (blue coloured circled marks) are plotted against the estimated RPW values and number of FEM simulation runs. The statistical model equation (Equation (8)) derived for the RPW using DoE factorial method provides only the direct linear relationship of the 4 structural parameters influencing the sensing performance. For the derivation of the statistical model for that RPW that includes the nonlinear terms and inter-connected structural parameters influencing terms, we have to utilise the Response Surface Methodology available in DoE.

DoE Response Surface Methodology
The Response Surface Methodology (RSM) was introduced by Box and Wilson [26] to derive the nonlinear statistical model equations to estimate the response of interest that will include the higher-order and interplay effects of the parameters. The RSM statistical method used in our work to perform this DoE procedure is known as central composite design (CCD) [34,35]. This statistical procedure utilises a centre point to design a quadratic model for the RPW response of the HC-PCF-based SPR sensor. The model is expressed in the form of a second-order quadratic equation [34]: where β 0 , β i , β ii and β ij represent constant, linear, intra-quadratic and inter-quadratic terms coefficients, respectively. Thus, for the CCD procedure, the vector β length becomes k+2 C 2 .
In addition, Y becomes a vector of length k+2 C 2 to be constructed with the FEM simulations RPW (response) values. Consequently, D becomes a matrix of size n × k+2 C 2 constructed with the offset constant coefficient and different levels of parameters in sequence of linear, intra-quadratic and inter-quadratic parameters values. Also for CCD, n = 2 k + 2k + 1.
In the CCD procedure, the 2 levels (maximum and minimum) used for factorial analysis are increased to 5 levels. The centre point is introduced as well as the embedded fractional factorial levels.
Using the data in Table 4, statistical analysis and numerical analysis for k = 4 parameters corresponding to n = 2 4 + 2 × 4 + 1 = 25 scenarios are tested using the Minitab statistical software. The variations of +1% and −1% levels for the structural parameters are same as described in Table 3 for factorial design. As required, more FEM simulations runs (rows 17-25) are carried out using the variations of +2%, −2% and 0% levels (Table 4) for the structural parameters as depicted in Table 5.
Using the structure parameters for the D matrix and the values for RPW (response Y vector from last column) in Table 5, Equation (6) can be written as: 1240  1241  1248  1250  1237  1238  1245  1248  1240  1240  1248  1250  1237  1238  1245  1248  1241  1244  1234  1252  1246  1241  1244  1243 where the D matrix values are calculated as 1(offset), d, Λ, t g , d c , d 2 , Λ 2 , t 2 g , d 2 c , dΛ, dt g , dd c , Λt g , Λd c and t g d c . Table 5 shows the RPWs from 25 different FEM simulation runs with the designed structural parameters for an analyte RI of n a = 1.45. Using the calculated value of β for the 5-level quadratic model Equation (9), we derive: The relationship between the DoE RPW model obtained in Equation (11) and the RPW-FEM simulations for the different structural parameters are shown in Figure 6.  Figure 6. Effects of four factors (air-hole diameter, pitch, gold layer thickness and core diameter of analyte) on RPW curves generated from quadratic CCD ±2% model. The dotted, solid and dashed curves show the RPW calculated from the ±2% DoE model (11) with three parameters (expected target parameter on x-axis) fixed at their own tolerance levels, −2%, 0% and +2%, respectively. The RPW marks ( −→ −2%, −→ 0% and −→ +2%) are calculated from corresponding FEM simulations.
The 15 RPW values calculated from FEM simulations are plotted on the corresponding RPW model (11) curves (5 marks for each characteristic curve) in Figure 6. It shows that the prediction from the derived DoE Equation (11) is a very good estimate. The dotted, solid and dashed curves show the RPWs calculated from DoE model (11) with ±2% variations for the PCF parameter represented by the x-axis label and with three parameters fixed at their own tolerance levels, −2%, 0% and +2%, respectively (see Table 4). The , and marks show the corresponding RPW values from the FEM simulation results for DoE model accuracy verification of those dotted, solid and dashed lines.
From the subplots, except the one for the pitch, the curves for the RPW in the air-hole diameter, gold layer thickness and core diameter plots are separated from each other without any crossing. The RPW lines in the air-hole diameter subplot are having more curvature, which shows that air-hole diameter has a slight quadratic effect influence for the RPW estimation for the ±2% DoE model. For the gold layer thickness and core diameter subplots, the equally spaced straight lines show that these two parameters are dominated only by the linear effects on the derived RPW for the ±2% DoE model.
We can conclude that pitch plays the most important role (nonlinear influence) in determining the RPW for the ±2% DoE model. However, there is possibility for the characteristics from the air-hole diameter plot to have crosspoint(s) at higher tolerance levels as weak quadratic nonlinear effect is depicted at ±2% variations in the values for the structural parameters. To investigate further on the tolerance studies, we have performed the DoE statistical modelling for two higher level variations, respectively, for ±5% and ±10% ranges. •

DoE Model with ±5%
The DoE quadratic model of a higher tolerance (±5%) is investigated to further find the possible best fit model for fabrication purposes. The full details for the FEM simulation of the HC-PCF-SPR sensor with (±5%) fabrication tolerance are provided in the following and the DoE equation derived is of same functional form like the ±2% variations, but with a different set of the constant and coefficients.
The DoE modelling calculation of the HC-PCF-SPR sensor with ±5% tolerance is presented in detail. Table 6 shows the 5 levels design of four parameters with ±5% fabrication tolerance. Using the data in Table 6, statistical analysis and numerical analysis for 25 scenarios are tested using the Minitab statistical software.
Using the structure parameters for the D matrix and for RPW (response Y vector from last column) in Table 7, Equation (6) can be written similar to (10). Table 7 shows the RPWs from 25 different FEM simulation runs with the designed structural parameters for an analyte RI of n a = 1.45. For the calculated value of β for the 5-level quadratic model (9), we derive:  Figure 7. Effects of four factors (air-hole diameter, pitch, gold layer thickness and core diameter of analyte) on RPW curves generated from quadratic CCD ±5% model. The dotted, solid and dashed curves show the RPW calculated from the ±5% DoE model (12) with three parameters (expected target parameter on x-axis) fixed at their own tolerance levels, −5%, 0% and +5%, respectively. The RPW marks ( −→ −5%, −→ 0% and −→ +5%) are calculated from corresponding FEM simulations.
The 15 RPW values calculated from FEM simulations are plotted on the corresponding RPW model (12) curves (5 marks for each characteristics curve) in Figure 7. It shows that the prediction from the derived DoE Equation (12) is a very good estimate. The dotted, solid and dashed curves show the RPWs calculated from DoE model (12) with ±5% variations for the PCF parameter represented by the x-axis label and with three parameters fixed at their own tolerance levels, −5%, 0% and +5%, respectively (see Table 6). The , and marks show the corresponding RPW values from the FEM simulation results for DoE model accuracy verification of those dotted, solid and dashed lines.
From the subplots, except the air-hole diameter and pitch in this case, the curves for the RPW in the gold layer thickness and core diameter plots are still separated from each other without any crossing. The RPW lines in the air-hole diameter subplot are having more curvature and dashed line (+5%) and solid line (0%) will have a crossing for lower air-hole diameter which shows that air-hole diameter has stronger quadratic effect influence for the RPW estimation in the ±5% DoE model compared with the ±2% DoE model.
We can summarise that besides pitch, air-hole diameter plays also one of the important roles (nonlinear influence) in determining the RPW for the ±5% DoE model. The possibility for the characteristics from the air-hole diameter plot to have crosspoint(s) at higher tolerance levels as weak quadratic nonlinear effect can be expected at ±5% variations in the values for the structural parameters. However, this needs further confirmation. To investigate further on the tolerance studies, we have performed the DoE statistical modelling for another higher level of variation (±10% range). •

DoE Model with ±10%
The DoE quadratic model of a higher tolerance (±10%) is investigated to further find the possible best fit model for fabrication purposes. Similar to the previous case for the ±5% tolerance study, the full details for the FEM simulation of the HC-PCF-SPR sensor with (±10%) fabrication tolerance are provided in the following and the DoE equation derived is of same functional form like the previous ±2% and ±5% fabrication tolerances, but with a different set of the constant and coefficients. The DoE modelling calculation of the HC-PCF-SPR sensor with ±10% tolerance is presented in detail. Table 8 shows the 5-level design of four parameters with ±10% fabrication tolerance. Using the data in Table 8, statistical analysis and numerical analysis for 25 scenarios are tested using the Minitab statistical software.
Using the structure parameters for the D matrix and the values for RPW (response Y vector from last column) in Table 9, Equation (6) can be written similar to (10). Table 9 shows the RPWs from 25 different FEM simulation runs with the designed structural parameters for an analyte RI of n a = 1.45. For the calculated value of β for the 5-level quadratic model (9), we derive: The relationship between the DoE RPW models obtained in Equation (13) Figure 8. Effects of four factors (air-hole diameter, pitch, gold layer thickness and core diameter of analyte) on RPW curves generated from quadratic CCD ±10% model. The dotted, solid and dashed curves show the RPW calculated from the ±2% DoE model (13) with three parameters (expected target parameter on x-axis) fixed at their own tolerance levels, −10%, 0% and +10%, respectively. The RPW marks ( −→ −10%, −→ 0% and −→ +10%) are calculated from corresponding FEM simulations.
The 15 RPW values calculated from FEM simulations are plotted on the corresponding RPW model (13) curves (5 marks for each characteristics curve) in Figure 8. It shows that the prediction from the derived DoE Equation (13) is a very good estimate. The dotted, solid and dashed curves show the RPWs calculated from DoE model (13) with ±10% variations for the PCF parameter represented by the x-axis label and with three parameters fixed at their own tolerance levels, −10%, 0% and +10%, respectively (see Table 8). The , and marks show the corresponding RPW values from the FEM simulation results for DoE model accuracy verification of those dotted, solid and dashed lines.
From the subplots, except the air-hole diameter and pitch in this case, the curves for the RPW in the gold layer thickness and core diameter plots are separated from each other without any crossing as in the previous case. The RPW lines in the air-hole diameter subplot are having more curvature and dashed line (+5%) and solid line (0%) now have two crossings in lower air-hole diameter region which confirms that the air-hole diameter has strong quadratic effect influence for the RPW estimation in the ±10% DoE model compared with the ±2% and ±5% DoE model.
We can conclude that pitch and air-hole diameter not only play important roles (nonlinear influence) in determining the RPW for the ±5% DoE model, but also they both exhibit the characteristics in the subplots to have crosspoints at higher tolerance levels as quadratic nonlinear effect can be expected at ±10% variations in the values for the structural parameters. The level of quadratic nonlinear effect caused by air-hole diameter increases as the tolerance level increases. Moreover, the crosspoints and the strong quadratic nonlinear effect from air-hole diameter on the RPW cannot be easily seen/inferred for a DoE model with low tolerance.
Comparison between the coefficients of the three DoE models with different fabrication tolerances are listed in Table 10. High absolute values for β 1 , β 11 , β 2 and β 22 (first-and second-order coefficients for air-hole diameter and pitch, respectively) can be found for all three models. This matches with the conclusion that those two structure parameters play the dominating roles for the determination of RPW. It should be noted that for quadratic coefficient of pitch (β 22 ), it kept almost unchanged for all three tolerances with high absolute values, whereas the second-order coefficient of air-hole diameter (β 11 ) increases drastically in the negative direction. These characteristics of the coefficients confirm the conclusions arrived from the behaviour noticed from the graphs for the quadratic effect behaviour of pitch and air-hole diameter on the RPW. With respect to the coefficients β 4 and β 44 corresponding to the first-order and second-order terms for the gold layer thickness, even though, both the absolute values and changes are significant, the role played by this parameter has barely any quadratic effect on the determination of the RPW as the unit for the metal coating (nm) is three orders lower than the other three parameters (µm). Therefore, the quadratic effects for parameters of gold layer thickness and core diameter with small coefficients are not crucial for the determination of the RPW, which again confirms the findings from the graphs.  (11)), ±5% (Equation (12)) and ±10% (Equation (13) Figure 9 illustrates the RPW response surface (3D) and contour (2D) plots for the interaction of any two among the four structural parameters of the proposed HC-PCF SPR sensor, corresponding to air-hole diameter, pitch, gold layer thickness and core diameter tolerance levels of ±2%, ±5% and ±10%. From these plots, it is noticeable that the ±10% DoE model can provide overall characteristics of the RPW variations with any two structural parameters and it can represent most of the parametric ranges that are provided by the ±2% and ±5% DoE models. Thus, Equation (13) corresponding to ±10% DoE model is proved to be the general RPW equation that also includes ±2% and ±5% tolerance level cases. However, it should be noted that finer detail differences can still be found between the surfaces for different tolerance levels due to the model residual error defined in Equation (9). Particularly, this is observed in ±5% and ±10% surfaces in Figure 9c and ±2% and ±5% surfaces in Figure 9f. Hence, the proposed statistical modelling of the PCF-SPR sensor response (RPW in this work) for the structural parametric range to be developed is based on requirements such as the overall general characteristics for larger range variations or finer details for particular small range variations. . Response surface (3D) and contour (2D) plots of RPW with respect to the interaction of any two among the four structural parameters (air-hole diameter, pitch, gold layer thickness and core diameter) for ±2%, ±5% and ±10% tolerance levels. Different levels of transparency for the surface and contour plots are set for identification as 0% (for ±2%), 25% (for ±5%), 50% (for ±10%).
Results from the DoE models for three tolerance levels, ±2% (Equation (11)), ±5% (Equation (12)) and ±10% (Equation (13)) as illustrated in Figures 6-8, show that DoE model can present good predictions and exhibit linear or quadratic influences on the determination of the RPW of HC-PCF-based SPR sensors. For the proposed HC-PCF sensor, pitch has been discovered to have the most quadratic effect in determining the RPW of this sensor among all three tolerance levels. Gold layer thickness and core diameter were discovered to have very steady linear effects as tolerance increases. This means that during the fabrication of HC-PCF-based SPR sensors, structural variations of the gold layer thickness and core diameter from tolerances ±2%, ±5% to ±10% will not affect the stability of RPWs and, hence, the sensing performances of the sensor. On the other hand, the air-hole diameter will cause comparatively significant quadratic effect on the RPW of HC-PCF-based SPR sensors for high fabrication tolerance level (±10%), moderate quadratic effect with typical fabrication tolerance level (±5%) and the quadratic effect can hardly be found for low fabrication tolerance level (±2%).

Conclusions
In this paper, we presented a statistical approach of determining the influence of the optical fibre structural design parameters on the sensing performance of HC-PCF-based SPR sensors at three tolerance levels (±2%, ±5% and ±10%). Two statistical DoE methods, factorial design and RSM design, were presented in detail to derive the model equation(s) for the RPW in terms of the structural parameters of the HC-PCF-based SPR sensor for the tolerance studies. The whole work begun with FEM simulation for the sensor design which could achieve a maximum sensitivity of 5600 nm/RIU for analytes detection with RI range of 1.43 to 1.49. Then, factorial design and CCD were demonstrated for the modelling and calculations procedures. The DoE models based on FEM simulation results for RPW estimation with three tolerance levels were analysed and compared in detail.
Results from the investigation on three models show that structural parameters (airhole diameter, pitch, gold layer thickness and core diameter) have linear or quadratic effects on the RPW determination for the sensing process. A clear understanding of the individual parameters influence in determining the RPW changes was presented with the DoE methodology combined with a few FEM simulation results which would not be feasible with numerical computations only. For the proposed HC-PCF sensor, pitch was the most significant parameter for the quadratic effect influence on RPW and gold layer thickness and core diameter are the least contributing parameters. The significance of quadratic effect influence for air-hole diameter was found to be proportional to the increase of fabrication tolerance level. Due to the amount of time and computational resource saved by using the proposed research methodology on tolerance study, we envisage that this novel statistical model derivation of sensing performance parameters can be applied in fibre sensor fabrication industry to improve efficiency and lower cost.

Conflicts of Interest:
The authors declare no conflict of interest.