Comprehensive Modeling of Multimode Fiber Sensors for Refractive Index Measurement and Experimental Validation

We propose and develop a comprehensive model for estimating the refractive index (RI) response over three potential sensing zones in a multimode fiber. The model has been developed based on a combined ray optics, Gaussian beam, and wave optics analysis coupled to the consideration of the injected interrogating lightwave characteristics and validated experimentally through the realization of three sensors with different lengths of stripped cladding sections as the sensing region. The experimental results highly corroborate and validate the simulation output from the model for the three RI sensing zones. The sensors can be employed over a very wide dynamic RI range from 1.316 to over 1.608 at a wavelength of 1550 nm, with the best resolution of 2.2406 × 10−5 RI unit (RIU) obtained in Zone II for a 1-cm sensor length.

region 14 , the refractometer can be used to detect methane (CH 4 ) concentration. The cryptophane-based molecular traps will absorb or entrap the CH 4 molecules, and reversibly produce a bond in the bulk polymeric material that will induce variation in the RI of the sensitized region as a function of CH 4 concentration.
Various operating principles have been exploited for fiber RI sensing such as employing Fresnel reflection at the end face of a single-mode fiber (SMF) 8 , surface plasmon resonance (SPR) [15][16][17][18][19] , tapered optical fibers [20][21][22][23] , multimode interference (MMI) 20,[24][25][26] , fiber Brag gratings (FBGs) 27,28 , long-period gratings (LPGs) 29 , etc. A further technique that is simple, fast to realize, and cost effective is by employing a totally or partially stripped-cladding core as the sensing region in a multimode fiber (MMF) [30][31][32][33][34] . Previous work has reported typical dynamic ranges from 1.33-1.55 RIU with resolutions between 6.48 × 10 −3 RIU and 9.68 × 10 −4 RIU for different fiber types 30,31 . The de-cladded section is generally substituted by a sensing medium which interacts with lightwaves propagated in the fiber. Through preliminary experimental results 34 , we have demonstrated and reported the classification of three operating zones in the RI response in a stripped-cladding step index MMF subject to variation in the external medium's RI in contact with the sensing element, according to three different sensing mechanisms. This three-zone phenomenon is theoretically reviewed, explained and experimentally validated in this paper, the outcome of which can potentially be employed for detecting RI variation of any external medium subject to or undergoing perturbation. More advantageously, it can enable RI measurement over a very wide dynamic range, virtually unlimited by the RI of the MMF core since detection is always or still possible even when the external medium's RI is greater than that of the core.
In previous work 34 , Zone I has been defined for a range of variation of the external medium's RI up to the cladding index of the MMF. The variation of the power guided in the fiber, induced by the change in RI is due to evanescent wave absorption (EWA). In Zone II, where the induced RI variation is situated between the cladding RI (n cl ) and the core RI (n co ), i.e. greater than n cl but less than n co , we have demonstrated that there are two concurrent optical (or in this case sensing) phenomena occurring caused by EWA and the reduction of the number of propagation modes due to the modification of the critical angle as a consequence of the RI change in the sensing or external medium. The final RI response is Zone III where the RI of the external medium is higher than n co . The power guided for this condition is due to the reflection of the "lost" lightwaves, which have been transmitted to the external medium interface, back into the fiber core. This is the most logical hypothesis to explain the existence of the RI response for indices higher than the core index (equally demonstrated experimentally by Mukherjee et al. 35 ), since under this operating condition, there is no (further) EWA as no further lightwave propagation exists by total internal reflection (TIR). The explanation of the RI sensing behavior is backed up by experimental evidence, via the transmission characteristics of the propagating lightwaves through the fiber.
In this paper, a comprehensive study is carried out on step-index MMFs for RI sensing where the light source is a single-mode DFB laser diode pigtailed to a single-mode fiber (SMF). A microscope objective (MO) is employed to inject the required lightwave modes into the MMF under investigation. From this launch condition, we propose an effective model for this system and define the three sensing zones by stripping the MMF cladding. The model is adapted from the combination of the analytical equations for wave optics, beam optics (Gaussian beam), and ray optics. The motivation in proposing this model is to overcome the complexity in the simulation of propagating lightwaves in MMFs by wave optics, especially in the case of simulation by the finite element method (FEM) in 3D since the mesh size must be smaller than the wavelength. Furthermore, the dimension of MMFs will contribute to a very high number of unknowns, implying a very significant computational load which is difficult to resolve by generic computing. Hence, the simplest way to explain the propagation principle and the power loss phenomenon in MMFs as RI sensors is by ray optics. Nevertheless, EWA and the beam distribution injected into MMFs, which cannot be determined using ray optic analysis, can be compensated by concurrently applying the analytical equations of the evanescent wave, and employing the Gaussian beam equations to define the injected beam power distribution in an MMF. Using the model developed, we demonstrate accurate simulation results which are corroborated by experimental data.

Results
Propagation in multimode optical fibers. The laser beam is propagated in an optical fiber by TIR which occurs as the beam injected into the fiber core is incident at the core-cladding interface at an angle higher than the critical angle. This critical angle, θ c , can be explained by Snell's law as a function of the RI contrast (n co and n cl , respectively) as follows: For an SMF, only one mode propagates in the fiber. However, in an MMF numerous modes can propagate. The propagating modes are modes that are incident at the core-cladding interface with an angle between θ c and 90°. Theoretically, in wave optics the number of modes (M) for a very large number of modes in a step-index MMF can be estimated by 36 where V is an all-important parameter for determining whether a fiber is single or multi mode. For V < 2.405, the fiber will be single-mode while an MMF has a V-parameter greater than 2.405. The value of V relates to n co , n cl , the fiber core radius (a), and the wavelength (λ) of the injected beam as where − n n co cl 2 2 is also known as the numerical aperture (NA), which defines the acceptance angle within which the injected beam via the MO can be propagated by the fiber or radiated by the fiber. Hence, to obtain all the modes over the entire possible acceptance angle, the injected beam must be adapted to the NA of the step-index MMF.
Experimental set-up. During the measurements, a differential probe configuration is employed, as shown in Fig. 1, with one MMF serving as the reference fiber and the other MMF the sensing fiber to compensate common-mode noise produced by both MMFs. The sensor, based on the ratiometric intensities from the two MMFs, thus measures a transmitted power for the reference MMF and the power in the sensing MMF arm using two identical Ge-type Thorlabs PDA50B photodetectors (PDs) with a peak spectral response within 800-1800 nm. A single-mode fiber-pigtailed DFB laser diode from Modulight, Inc. emitting at 1550 nm and driven by precision current and temperature controllers is employed to interrogate the fibers. The laser output beam is divided by a single-mode fiber coupler to obtain two equal or symmetrical beams which are first collimated and then injected into the reference and sensing MMF arms of the sensor via two identical MOs with an NA of 0.65 in order to transmit the launched beam over all acceptance angles in the fibers. The MMFs are plastic-clad silica (PCS) fibers exhibiting an NA of 0.48 ± 0.02, with 200 µm core diameter and 230 µm cladding diameter. The output power detected from both MMF arms are then transmitted to a dedicated computer via a 14-bit 2-MHz Agilent data acquisition system (DAQ) for further processing.
The MMF sensing region is realized by thermally removing a certain length of the fiber cladding (as shown in the inset of Fig. 1), followed by a simple procedure of cleaning the fiber core using acetone and isopropanol solutions. Although the buffer still contains some filaments of the coating material, both the core and cladding sections are smooth (i.e. continous) and undamaged, and hence will not influence or affect the propagation of the launched lightwaves within the sensing region. Three different sensing lengths of 1 cm, 2.5 cm, and 4 cm of stripped fiber cladding have been prepared. This sensing area is thus sensitive to variations in the RI value through variations of the power transmitted as a function of the medium's RI. As previously mentioned, 3 different conditions or sensing phenomena can intervene in the MMF RI sensor. These detection mechanisms are further elaborated below.
Sensing mechanism in Zone I. Consider the case of the sensing medium having an RI value (n sm ) less than the cladding RI. Under this condition, all the modes propagate in the fiber core up to the sensing zone by the phenomenon of TIR, at which point evanescent waves will be generated along the core-cladding interface, as shown in Fig. 2. Figure 1. Experimental schematic set-up of MMF-based sensor for measuring RI under differential configuration to compensate common-mode noise. The inset shows an enlarged image of the stripped-cladding core of the MMF used in this work, with the junctions between the core and cladding, and between the cladding and buffer/coating clearly illustrated. There is, in particular, no spurious transitions between the core and cladding (i.e. no obvious filaments left over from the cladding removal process). www.nature.com/scientificreports www.nature.com/scientificreports/ When the incident angle (θ) is larger than θ c , almost all power is confined within the core while the previously transmitted part of the beam propagates as evanescent waves penetrating into the cladding. Hence, the resulting power transmitted to the fiber end while propagating over a stripped cladding of length, L, can be predicted by 37 where γ is the evanescent wave's attenuation coefficient and P 0 the initial power propagated by TIR in the MMF. For a weakly guiding fiber, i.e. for a fiber with only a small RI difference between the cladding and core, γ is a function of the bulk absorption coefficient (α) and the fractional power in the cladding and core (r), i.e. the ratio of the power in the cladding to the total power, given by γ = r.α. Now, r can be estimated by a simple equation related to the number of modes (M) using r = 4/(3 M ) 37 . Nevertheless, since the MMF is generally a non-weakly guiding fiber, γ thus has to be re-defined to account for the number of modes with a transmission coefficient (T) penetrating into the cladding over a number of reflections per unit length (N) as a function of the critical angle of the sensing medium (θ csm = sin −1 (n sm /n co )) in replacement of the cladding and the incident angle (θ), for θ > θ c 37 such that Further, by adapting Equations (5), (6) and (7) to solve Equation (4) for non-weakly guiding fibers, the power over a sensing length, L, can be calculated as a function of the angle over all the acceptance angles (θ, π/2) as In this case, due to the comparatively low power transmission loss in the MMF over a length of ~80 cm, the absorption effect by the fiber cladding is assumed to be negligible (attenuation ~0.028 dB/m from manufacturer's datasheet). Hence, the initial power can be taken as the combined power of all the rays launched into the fiber under TIR as shown in Fig. 3.
The single-mode Gaussian beam from the SMF is collimated by an F220FC-1550 collimator to a beam diameter (φ = 2w′) 38 of 2.1 mm (red beam in Fig. 3(b) to the left of the lens) while the diameter of the total beam power can be calculated by D = πw′. This collimated beam (consisting of both the red and yellow beams) is next focused by an MO with an NA of 0.65 to match all possible acceptance angles for propagation in the MMF. The beam waist (w 0 ) can be estimated by Equation (9) which relates λ, the focal length of the MO (f), and the spot radius of the collimated beam (w′) in the form 36 As illustrated in Fig. 3(b), when the total lens area of the MO is completely illuminated by the input collimated beam, the injection angle β′ = β 0 , can be obtained from the relationship, NA = n air sin β 0 . However, since the input beam does not illuminate the entire lens area due to the lens diameter (d lens ≈ 5.5 mm) being larger than the total beam diameter (D = πw′ = 3.2987 mm), β′ cannot, consequently, be calculated directly using the NA equation which gives a β 0 value of ~40.54°. By simplifying the MO to a one-lens system, f can first be obtained using Equation (10) relating β 0 to the lens radius (r lens = 0.5d lens ), after which β′ can then be obtained for a given D value via Equation (11) In ray optics, mode propagation in the MMF is represented by individual rays at specific angles injected into and, subsequently, propagating in the fiber, and can thus be simply analyzed by considering the beam with an incident θ at the core-cladding interface. Rays which have incident angles less than θ c of the fiber will not be propagated in the core, when compared to rays with angles greater than θ c . These latter rays will be guided in the core by TIR. To simplify the propagation model, the incident beam characteristics illustrated in Fig. 3(b) and (c) are analyzed. The beam with a propagating angle β in the MMF will arrive at the core-cladding interface at a distance z i , calculated using where β (=90° − θ) is related to β′ by Snell's law via n air sin β′ = n co sin β and i = 1 or 2. The incident Gaussian beam can be represented as a set of rays with different individual optical power values exhibiting a Gaussian distribution, where the outer rays have a lower power than the rays towards the center fiber axis. However, since the onset of TIR starts at z i , the diffraction effects of the Gaussian beam are annulled at z i 39 , thus the power distribution for each angle of the individual ray propagated along the MMF follows on from the ray power distribution at point z i , as will be explained below. Consequently, the sensing region must be placed after the point z i . Then, by using Equation (12) and considering that the beam is injected into the center of the fiber core with a = 100 µm, for a 1 cm sensing length, the beam will be incident at the interface of the sensing area at least once up to θ = 89.427°. Hence, almost all of the rays injected will reach the sensing region. Further, for a Gaussian beam injected from air and expanded in the fiber core (silica with n co = 1.444), the beam spot radius at z i (w zi ) can be calculated by 40 The optical intensity at z i (I zi ) and the optical power injected into the fiber (P inj ) can respectively then be given by 36,38 Subsequently, as illustrated in Fig. 3(b), when θ is lower than θ c of the fiber, the total power arrives at the core-cladding interface at the point z 1 . Only the optical power between θ and θ c will exit the core while the remaining power will still be guided in the fiber core by TIR represented as the beam with a radius a′ which can be calculated by simple trigonometry as i c Substituting the new radius for a′ in Equation (15) then allows P 0 , the power guided by TIR for θ < θ c of the fiber, to be obtained as Figure 3. Schematic of beam launching conditions into MMF where the red beam represents the portion of the beam which transmits ~86% of the total power (or 1/e 2 beam intensity) and which is used to determine the spot size radius (w 0 ), while the yellow beam represents over 99% of the total beam power. (a) Collimated beam from lead-in SMF with Gaussian power distribution is focused by MO into the MMF for more specific launch parameters when (b) θ < θ c of the fiber and when (c) θ > θ c of the fiber.
However, if θ > θ c , all the incident beams arriving in the sensing region are reflected by TIR. Thus, a′ = a and P inj = P 0 . Equation (17) must then employ z i = z 2 (see Fig. 3(c)), and w zi takes on w z2 . Since the surface of the MO lens is not entirely illuminated by the collimated beam (d lens ≈ 5.5 mm, and D ≈ 3.2987 mm), θ can be obtained by applying Equations (10, 11) and Snell's law. Here, θ = 71.57° is higher than θ c (70.60°), thus all injected beams are guided by TIR. Consequently, the previous definition of Zone I 34 has to be re-evaluated, i.e. the limit of Zone I is modified to the RI value (n eq ) which is equivalent to the angle of the incident beam (via θ = sin −1 (n eq /n co )). Further, w z2 = 67.898 µm is obtained via Equation (13) for the point z 2 = 300.14 µm and w 0 = 1.51 µm obtained by Equation (9). Then, to solve Equation (8), P 0 in Equation (17) can be discretized into individual ray powers between θ and 90° (i.e. P 0 (θ, π/2)) at the point z 2 by varying a' in discrete steps (from 100 µm to 0) corresponding to the θ value for each ray. Substituting these step variations of a' into Equation (17) then allows the decremental power difference, P 0 (θ, π/2), to be calculated for each individual ray with a corresponding set of T(θ, π/2) and N(θ, π/2) parameters. Finally, these powers are summed over θ − 90° in Equation (8), to obtain the total transmitted power (P L ) subject to EWA for a given sensing length L.
Sensing mechanisms in Zone II. For the second operating condition or Zone II response, where the RI of the sensing medium (n sm ) falls between the cladding RI and the core RI (i.e., n cl < n sm < n co ), two optical phenomena concurrently influence the optical power loss 34 (see Fig. 4): (1) reduction of the number of propagation modes due to the modification of the critical angle (from θ c to θ csm ) upon variation of n sm and (2) EWA since TIR is always present in this operating zone, identical to the explanation of losses described for Zone I.
The reduction in the number of propagation modes occurs for increasing n sm since the critical angle also increases as θ csm = sin −1 (n sm /n co ). This leads to a decrease in the guided power by TIR (P 0 ) as the beam radius a' at point z 2 decreases (see Equation (17)), subsequently reducing the sectional area at this point as illustrated by Fig. 5. The beam radius, a′, is calculated by the following formula: The model for Zone II is subsequently obtained for different values of P 0 (θ csm , π/2) as a result of the reduction of the number of modes by modification of the critical angle via Equation (17) while the remaining power for rays from θ csm to 90° which are subject to EWA can be estimated using Equation (8) for a given α and sensing length L as in Zone I.
Sensing mechanism in Zone III. The last condition is Zone III where n sm > n co . Here, propagation by TIR is no longer supported although a very small portion of optical power can still be guided due to the phenomenon of external reflection. This can be explained by Fresnel equations for a beam incident at an interface between two media of different RI values as illustrated in Fig. 6.
Since n sm > n co , most of the optical power of the ray will be transmitted to the exterior. Furthermore, since the sensing medium is much thicker than the core diameter, the transmitted ray will not be back-reflected into the core. However, a small proportion of power is still reflected back into the core by the external reflection mechanism at the core-sensing medium interface. This proportion of back-reflected power can be obtained by calculating the reflectivity (R) for the light containing both P-polarization (r p ) and S-polarization (r s ) components using the Fresnel equations below 36  Figure 4. Illustration of sensing mechanisms in Zone II due to the combination of EWA and mode power loss by the modification of the critical angle (θ csm ). Under initial conditions, the blue line represents the mode propagated in the MMF by TIR at critical angle; however, when n sm increases, θ c is modified to θ csm and this mode will be transmitted through the medium and lost to the exterior. The remaining rays incident at an angle greater than the new θ csm will be guided in the MMF by TIR and concurrently subject to EWA. The reflection at the sensing medium-air interface is negligible since the medium thickness is more than 20 times the core diameter in addition to its strong absorption at 1550 nm. The dotted blue and red lines represent very weak reflections obtained by Fresnel equations at the core-medium boundary which are neglected for simplifying the model for Zone II. with θ i , θ t and θ r the incident, transmitted and reflected angles, respectively. Thus, for a given R, the power of the individual ray guided in the core for Zone III as a function of the total number of reflections (NL) at the core-sensing medium interface can be estimated by Figure 5. Evolution of optical power and intensity in the MMF RI sensor for various values of n sm by modification of the critical angle, θ csm , for n cl < n sm < n co . (a) increasing n sm will increase θ csm which reduces the acceptance angle of the propagating beam, (b) power in transversal plane decreases for increasing n sm , and (c) illustrates decreasing optical intensity over different a' obtained by Equation (14).

Figure 6.
Optical power guided by external reflection in Zone III when n sm > n co . The majority of the power is transmitted to the exterior while a small portion is reflected back into the core. This mechanism again arises when the ray intersects the interface between the sensing region and the external medium. This phenomenon is more prominent for rays which have incident angles close to 90°. Simulation and experimental results. A suite of simulations is then performed based on the above models and validated against experimental data for the three different sensing Zones. Figure 7 illustrates the high level of agreement between simulations and the experimentally-measured data. The sensing medium is composed of a combination of glycerol-water mixture which is adjusted to obtain different values of RI from 1.3164 to 1.4571 at a wavelength of 1550 nm 5 . For RI values beyond 1.4571, calibrated RI liquids were employed. Furthermore, different glycerol-water concentration levels exhibit varying absorption coefficients (α) with that of water being experimentally determined in a controlled environment to be 11.49 cm −1 , and that of glycerol being 11.10 cm −1 . Assuming a linear relationship between the imaginary part or extinction coefficient of the mixture's RI (k tot ) and its volume fraction, which is a function of the glycerol (f g ) and water (f w ) fractions (with f g + f w = 1), the extinction coefficients of glycerol (k g ) and water (k w ), we can estimate 41 k tot = f g k g + f w k w . The composite absorption coefficient of the solution can then be obtained via 41 The sensitivity curves of the RI sensors are next obtained for each length of sensing region and plotted in Fig. 8. These curves are obtained by differentiating the fitted curves through the experimental data in Fig. 7 based on our models for each zone.
Here, the sensitivity for each sensor has been obtained by first fitting specific response curves to the 3 respective sensing Zones in Fig. 8 followed by their differentiation with respect to RI.

Discussion
Three complete models for the three different sensing zones of an MMF refractometer for RI measurement have, for the first time, been developed as a function of the injected lightwave characteristics into the fiber. The mathematical models for Zone I and Zone II are similar as a consequence of direct influences from EWA. Although the principal model for EWA in Equation (8) has been employed for both Zone I and Zone II, they are subject to different initial powers P 0 (in Zone I, P 0 = P 0 (θ, π/2) for θ > θ c or P 0 = P 0 (θ c , π/2) if θ < θ c , while in Zone II, P 0 = P 0 (θ csm , π/2)).
For Zone I response, P 0 is conserved for each of the rays over the entire RI response of the sensing medium as all the rays within the sensing region are guided by TIR. Hence, if θ is less than θ c of the MMF, the injected rays will propagate in the sensing region, with a critical angle of ~70.60° (equivalent to a cladding RI of 1.362), as predicted by Equation (1). However, in this work, since θ = 71.57°, corresponding to an equivalent RI value of ~1.370, the Zone I response is extended from our preliminary results 34 to this new RI value. Subsequently, all incident rays from MO will be guided in the MMF core by TIR. The simulated results plotted in Fig. 7 demonstrate good agreement with the experimental data, with the blue, red, and black continuous lines corresponding to simulation for a 1 cm, 2.5 cm, and 4 cm stripped fiber cladding, respectively, while the symbols ' ' , ' ' , and 'o' represent the respective experimental measurements. Longer sensing lengths will incur higher losses, hence the power measured decreases when the RI increases due to the increasing transmission (T) penetrating into the cladding, i.e. more power is absorbed at higher RI, as described by Equations (7) and (8). For Zone II, where n cl < n sm < n co (n co = 1.444), the variation of the optical power guided for different values of n sm is due to two optical phenomena as illustrated in Fig. 4. Here, the rays arrive in the sensing region guided by TIR since the incident beam angle is 71.54° and thus, the beginning of Zone II is equivalent to an RI value of ~1.370. The green line in Fig. 7 represents the simulated initial power variation P 0 (θ csm , π/2) as a function of RI due uniquely to the modification of the critical angle and does not reflect the influence of EWA. However, the contribution from EWA in the sensing region will lead to a smaller additional decrease in the total guided power along the fiber. This is clearly demonstrated through simulations in the form of continuous blue, red, and black lines for L = 1 cm, 2.5 cm, and 4 cm, respectively, as predicted by Equation (8). Now, identical to the EWA phenomenon in Zone I, the longer sensing region in Zone II is subject to higher absorption (Fig. 7). The high agreement between the simulated results and experimental measurements clearly demonstrates that the higher power losses occurring in Zone II are due to the contribution of both modification of the critical angle and EWA. The contribution of losses due to modification of the critical angle is also found to be more dominant than EWA in Zone II. Further, the power variation is observed to increase slowly at the beginning of Zone II toward the end of the Zone where it then decreases sharply. This sharp decrease is due to the power distribution of the Gaussian beam which increases sharply from the sides or wings of the distribution (i.e. top hat radius) towards the peak (i.e. circular aperture) of the Gaussian curve, but which decreases at the peak area (i.e. circular aperture area) 38 . Hence, the optical power response in Fig. 7 will exhibit a flat response or inflexion point at exactly the boundary between Zone II and Zone III. The last operating regime of the MMF refractometer is Zone III which can be employed for probing a medium with n sm > n co . Under this condition, there is no propagation by TIR. Nevertheless, the incorporation of the Fresnel relations in Equations (19) and (20) into Equation (21) postulates the existence of guided power in the core by the phenomenon of external reflection, in particular, from the rays which have incident angles θ in the sensing medium close to 90°. Although the guided power is small, this will increase for further increases in the value of RI beyond that of n co . This is validated experimentally in Fig. 7 through the measurement of increasing power at the fiber output end as n sm increases beyond that of the core. Complementary to this, Equation (21) further correctly predicts the higher power guided over shorter lengths of the sensing region in the MMF since there are fewer reflections (NL) which, in turn, reduce the transmitted power or rays to the exterior through the sensing medium. Nonetheless, the respective discrepancies between the simulations and the experimental results for the different sensing lengths in Zone III could be due to non-ideal conditions during the experimental study, such as the existence of very small bends in the MMF which can alter the optical power distribution and/or modify the beam quality factor (M 2 ) by inducing changes to the MMF index profile in the bending area 42 .
The sensitivity curves plotted in Fig. 8 illustrate the best sensitivity being achieved in Zone II for the shortest sensing length (i.e. 1-cm stripped fiber cladding) as the shortest sensor is subject to the least EWA. Consequently it suffers higher losses through modification of the critical angle (mode loss mechanism) as predicted by the continuous green line in Fig. 7. Toward the end of the Zone II response, there is more power variation for a small RI variation. The three sensitivity curves first increase sharply from the middle of Zone II, and then decline less sharply toward almost the end of this Zone before decreasing back toward zero at the core index (1.4444) which is the minimum point of the optical power response (see Fig. 7). The respective inflexion points of the sensitivity curves in Fig. 8 occur before the end of Zone II and correspond to the beginning of the decreasing gradient of the optical power response as described above with respect to the circular aperture area of the Gaussian beam. A www.nature.com/scientificreports www.nature.com/scientificreports/ zero sensitivity value could potentially be obtained when the power response in Fig. 7 occurred over very small RI variations (i.e. tending toward 0 or ΔRI → 0) of the sensing medium. However, since the practical RI variations induced in this work cannot be infinitesimal, the sensitivity of the three sensing lengths obtained at the end of Zone II cannot reach zero value. Nevertheless, the sensitivity curves as plotted in Fig. 8 decrease toward zero when the measured RIs approach the end of Zone II (i.e. close to the core index).
Conversely, for Zone I, the longest sensor has better sensitivity since the principally EWA contributions to the sensing mechanism are cumulative over the entire sensing length. There is thus more absorption by the longest sensor resulting in the largest power variation as a function of RI. In Zone III, on the contrary, the shortest sensor is again more sensitive since less power is lost to the exterior, and consequently more power is guided in the fiber core. Hence, according to Equations (19), (20) and (21), the guided power increases with increasing RI, with the increasing power being sharper at the beginning of Zone III, and subsequently declines less sharply with increasing RI. Thus, the sensitivity decreases with increasing RI in the sensing medium since the rate of power variation with RI decreases. However, although this sensor theoretically has virtually unlimited dynamic range for operation over Zone III, its performance could be limited to only a certain RI range when the sensitivity approaches the noise level.
Based on the sensitivity curves in Fig. 8, the sensor resolution has been determined with respect to the measurement noise level using the 6σ-definition (6 times RMS noise corresponding to ~99.7% confidence level) 43 with only ~0.3% of the samples lying outside of this distribution. The best resolution achieved is 2.2447 × 10 −5 RIU by the 1-cm sensor in Zone II. It is also in this Zone that the 2.5-cm and 4-cm sensors have the best relative resolutions of 2.9847 × 10 −5 RIU and 3.2517 × 10 −5 RIU, respectively, compared to the other two Zones. For Zone I, the best resolution is achieved by the 4-cm long sensor with a minimum detection level of 1.6116 × 10 −3 RIU while the 1-cm and 2.5-cm sensors are capable of resolutions of 5.5905 × 10 −3 RIU and 1.7528 × 10 −3 RIU, respectively. The achievable sensor resolution in Zone I is not very high due to the induced multiplicative noise from multiple reflections in the sensing region as well as the relatively low sensitivity in this zone. For the Zone III response, the noise level is relatively low since most of the injected power, including the noise from multiple reflections, are transmitted to the exterior. In this Zone, the normalized noise level ranges from approximately 1.2 × 10 −4 (a.u.) − 1.6 × 10 −4 (a.u.) and is typically dominated by the measured photodetector noise normalized to 8.66 × 10 −5 (a.u.). The minimum RI resolution that can be detected in Zone III are 1.0031 × 10 −3 RIU, 1.8070 × 10 −3 RIU, and 3.1920 × 10 −3 RIU for the 1-cm, 2.5-cm, and 4-cm sensing lengths, respectively. The low resolution obtained in Zone III can simply be understood by the low sensitivity in this zone as a consequence of higher losses arising from external reflection, as explained previously.

conclusions
We have proposed three models based on a combination of analytical wave optics to obtain the EWA equation, Gaussian beam optics to describe the injected power distribution, and ray optics to explain the principle of optical mode losses in an MMF configured for refractive index measurements. These models have been adapted to consider the three different sensing mechanisms as a function of the relative cladding and core RIs. Nonetheless, the models for Zone I and Zone II are fundamentally similar, whereby both Zones are subject to EWA as the fundamental loss mechanism. However, Zone II involves the additional phenomenon of critical angle modification, which modifies the model employed through the use of different values of P 0 (θ csm , π/2) as a function of RI variation. Further, since the incident beam angle in the fiber is higher than θ c , the boundary between Zone I and Zone II is no longer the cladding RI value (1.362), but the RI which corresponds to the incident beam angle (RI~1.370). Finally, the model for Zone III exploits Fresnel relations, where the rays propagating in the sensing region exhibit different power variations as a function of the ray angle with respect to their initial P 0 (θ, π/2) for an acceptance angle carried over from Zone II.
The experimental measurements performed are found to validate the simulation results derived from our models to describe the three different optical sensing mechanisms in the MMF refractometer. The results confirm that in Zone I, the sensing mechanism is uniquely via EWA which induces the largest losses in the sensor with the longest sensing region. For Zone II, the best sensor resolution of 2.2406 × 10 −5 RIU is achieved for the 1-cm sensor. The sharp power decrease occurring in Zone II is a consequence of the losses induced by modification of the critical angle for a Gaussian beam, where most power is concentrated at the center axis (top hat area), corresponding to an incident angle close to 90° (i.e. close to the core RI). However, at the beginning of Zone II, the losses are relatively small due to the weaker power distribution at the edges of the Gaussian beam which corresponds to an RI approaching that of the cladding, such that the losses are dominated by EWA. Last, but not least, in Zone III, when the external reflection mechanism intervenes, only a relatively small initial guided power exists, which subsequently increases as the external medium's RI increases due to the increasing reflectivity of the rays back into the fiber core. As predicted by our model, the longest sensor will guide less optical power since more reflections are induced by a longer sensing region, resulting in more rays being transmitted toward the exterior.
Future work will undertake modeling and analysis of the MMF refractometer described above for practical in-situ applications, for example, for detecting dissolved methane in aqueous environments. This would involve employing a thin PDMS film incorporating Cryptophane-A molecular traps as the sensing region whose bulk RI (~1.40 at 1550 nm) varies with varying methane concentration. We anticipate a potentially achievable measurement resolution of 1.5592 × 10 −4 RIU at 1550 nm by the 1-cm MMF sensor. This translates to ~28.35 nM of dissolved methane for a specified sensitivity of 5.5 × 10 −6 RIU/nM 44 . This performance would also be improved several-fold by exploiting lock-in detection techniques.

Methods
Simulation method. Simulations are realized for 3 sensing lengths of the stripped-cladding MMF through the use of the equations in the Results Section which are adapted to efficiently and accurately describe the sensing SCIENTIFIC REpoRtS | (2018) 8:5912 | DOI:10.1038/s41598-018-24153-0 mechanisms induced by the three optical phenomena in the fiber. The simulation results illustrate the sensor response in terms of the normalized optical power as a function of RI. The simulation procedure has been carried out according to the following flow chart shown in Fig. 9.
// Start program // Enter input parameters such as a, n air , n co , n cl , λ, NA, φ, D, r lens , n w , n gly , and β 0 // Incident angle (θ = 90° − β) obtained by determining f in Equation (10), β′ in Equation (11) and from Snell's law // If θ < θ c , Zone I limit is cladding RI value, but if θ > θ c , Zone I limit is RI equivalent of θ // To obtain power distribution, calculate spot radius with Equation (16) then determine power with Equation (17), and discretize this power for individual rays.
Experimental system. The interrogating laser beam from a Modulight, Inc. laser diode is injected into the MMFs as shown in Fig. 3. Three-axis MDE122 Martock translation stages (50 nm resolution) from Elliot Scientific are used for the launching of the laser beam into the MMF. The beam injection condition is optimized using adjustable-gain transimpedance photodetectors from Thorlabs Inc., when the highest transmitted power is detected at the fiber output end. The experimental set-up was mounted and stripped over at least 5 times to demonstrate the consistency and repeatability of the results. Measurement procedure. The stripped-cladding area is used to measure the variation of RI in liquids.
A combination of pure distilled water and pure glycerol from Sigma-Aldrich is used to obtain a range of RIs  5 which can be accessed from http://faculty.chem.queensu.ca/ people/faculty/loock/publications.htm. According to the data, a plot of RI as a function of the concentration of glycerol in water at 25 °C is traced, from which we then obtained the relationship between RI and the glycerol concentration levels. A homogeneous sensing medium (water + glycerol) is then obtained with the aid of a magnetic stirrer operating for ten minutes for each mixing process. This sensing medium is next calibrated against a reference Hanna Instruments optical refractometer incorporating automatic temperature compensation from 10 °C to 40 °C at an operating wavelength of 589 nm. The data from the reference refractometer are subsequently adjusted through the relationship obtained in Hoyt 45 to determine the mixture's RI at 1550 nm wavelength. For RI values beyond 1.4571, calibrated oils from Cargille Laboratories have been employed in our experimental measurements.
Data availability. Both simulation and experimental data are publicly available in the "supplementary information" file.