Numerical simulation of light propagation in metal-coated SNOM tips

Presented are the results of numerical simulations accomplished to investigate the propagation of electromagnetic excitations in certain types of metal-coated tapered tips terminating SiO2 multimode optical fibres with a subwavelength output aperture. The numerical simulations were initiated in order to enable better interpretation of previously reported experimental results concerning some features of the mesoscopic effect of spectral modulation observed for a broadband light transmitted by such tips. This effect occurs due to the interference between a small number of waveguide modes exiting a metal-coated tip, and the experimental results indicate a possible mode-selective photon-plasmon coupling in the studied tips. To match the experimental conditions, the tips were modelled for the light wavelength of 800 nm as three-layer systems (with the intermediate adhesion Cr layer and the outer layer of Al or Au). However, due to computational restrictions the end of a tip, only 18 μm long (most significant), was modelled. Numerical simulations yielded the dependences of propagation and attenuation constants on the fibre core radius for the most intensive (both photonic and plasmonic) output modes. The pairs of modes most probably contributing to the observed spectral modulation were identified. Although the simulations did not reveal any explicit mode coupling, the imperfections of real tips can cause mode transformations implying possible involvement of more than two modes. The thin (20 nm) Cr layer plays the main role for plasmonic modes generated on its SiO2 interface, which explains the small outer metal layer influence on the observed modal dispersion.


INTRODUCTION
* Tapered metal-coated tips with subwavelength output aperture terminating an optical fibre are widely used as aperture probes for scanning near-field optical microscopy (SNOM) [1].At certain conditions, the "mesoscopic effect of spectral modulation" (MSM) can be observed with a broadband light transmitted by such a SNOM tip [2].This effect is attributed to the interference of a small number of different photonic modes in the transmitted light.As it has been theoretically shown by Novotny and Hafner [3], the number of different light modes propagating in a metal-coated optical fibre is gradually reduced with decreasing the fibre's subwavelength core diameter.Due to this mode-filtering effect a certain spectral interval can exist for a certain SNOM tip with output aperture diameter d, in which virtually only two modes can be transmitted with significant and comparable amplitudes, resulting in MSM effect.Based on this principle, an experimental technique has been proposed [2] to enable investigation of intermodal dispersion features and plasmonic effects in a metal-coated SNOM tip terminating a multimode optical fibre.
Using this technique, a series of experimental investigations with several different SNOM tips have been performed in Institute of Physics, University of Tartu [4][5][6].The main result is that the intermodal dispersion in a short metal-coated tip region significantly exceeds by absolute value the one in the fibre tail, which can be attributed to a mode-dependent contribution of surface plasmon polaritons: one of the two involved modes couples to plasmons stronger than another one, resulting in a remarkably slower propagation.But some of the results obtained did not have an obvious or trivial interpretation, which motivated us to initiate and carry out the numerical simulations of light propagation in SNOM tips.The main experimental results with lacking interpretation are: (a) on the contrary to previous expectations, replacing Al with Au in the metal coating does not dramatically affect the intermodal dispersion in a SNOM tip, although this value is considered to have a strong plasmonic contribution; (b) experiments with different Cr/Al-coated SNOM tips reveal significantly stronger enhancement of intermodal dispersion in bent-type tips compared to a straight-type tip.
In the current paper, we review our numerical simulations initiated to get an insight on propagation of electromagnetic (EM) excitations in SNOM tips of the types used in the referred experimental works [4][5][6], which could enable better interpretation of the corresponding experimental results.For this reason, the chosen numerical model had to match as much as possible to the experimental conditions and the design of SNOM tips.
In Section 2, the procedure of our numerical simulations is described in mode details.The simulation results are presented and discussed in Section 3, and conclusions are drawn in Section 4.

NUMERICAL SIMULATION PROCEDURE
To match the experimental conditions, the numerical simulations were performed for the light wavelength of 800 nm, for Al-and Au-coated SNOM tips with output aperture diameters d of 200 nm and 150 nm.The SNOM tips were modelled as three-layer systems, as sketched in Fig. 1: first the fused silica core (refractive index n 0 = 1.4533 [7]), then a 20 nm thick chromium adhesion layer (n 1 = 4.0141 + 4.1771i [8]), and finally an infinite layer of aluminium (n 2 = 2.7075 + 8.2713i [8]) or gold (n 2 = 0.1535 + 4.907i [9]).In experiments, the actual outer thickness of the metal layer was 200 nm, which is more than sufficient for it to be virtually completely opaque -thus this layer could safely be modelled as an infinite one.The apex angle of the model SNOM tips was approximately adjusted to a value providing (by an order of magnitude) a match between calculated and experimentally measured transmissions of SNOM tips.The standard approach [3] of dividing a conical tip into many cylindrical disks was used to model a SNOM tip numerically.It is possible to calculate all the EM modes supported for every elementary cylinder by solving the Maxwell equations while assuming no reflections from the ends of the disk.In such a case the task simplifies modelling of an infinite-length round optical fibre.A comprehensive description of the procedure for finding all guided modes of a round optical fiber is presented in [10].We'll give only a short overview in this article.
The electrical field of every mode of a round optical fiber can be expressed as follows: where  and  are the polar coordinates, z is the coordinate along the fibre's optical axis, k z =  + i is the complex propagation constant, and  is the frequency.A similar expression can be written for the magnetic field.The azimuthal dependence of the field E(,) is given by periodic harmonic functions sin(m) and cos(m), where m is the mode's azimuthal order.The radial dependence is given by cylinder functions satisfying the second order Bessel differential equation.The solution of second order differential equation is a linear combination of two independent functions.There are several choices for the two independent functions: Bessel J m (k i ) and Neumann N m functions or Hankel functions of the first (1)   m

H
and the second kind (2)   m H .
The argument of all these functions is k i , where is the radial component of the wave vector of i-th layer, n i is the refractive index of this layer, and k 0 = 2/ is the vacuum wavenumber.The Bessel function J m , which is the only one with a finite value at the center ( = 0), is appropriate for the use in optical fiber core.Similarly, only the (1)   m H function fits for the outer layer (it describes an evanescent wave).
For the intermediate chromium layer, a linear combination of two Hankel functions has to be used.
The continuity requirement for tangential field components (E  , E  , H , H  ) results in a homogeneous linear equation system (described by matrix M) with 4(N -1) unknown field amplitudes, where N is the number of layers (N = 3 in our case).Such a homogeneous system has a solution only if its determinant vanishes.Thus the task is reduced to find the complex propagation constants satisfying the criterion det[M(k z )] = 0.
To solve this task, a C++ program has been developed.For calculation of cylindrical functions of complex argument, the AMOS Fortran library (http://netlib.org/amos/)was ported to C++.Due to the very high density of modes at larger core radiuses, finding the determinant's zeros and matching the modes of different radiuses turned out to be a highly complex task.But finally, alongside with the standard Newton's optimization method, an elegant algorithm developed by Delves and Lyness [11] was applied, which employs the Cauchy's argument principle to calculate number of solutions in complex plane by contour integration.Because of the high mode density and the stability of complex cylindrical functions, the simulation was limited to the maximum core diameter of 18 μm (see Fig. 1), thus only the very end of a SNOM tip could be modelled.Fortunately, this part of a tip also appears to have the most significant effect on the mode structure.
The transmittance T i of a SNOM tip for mode i was calculated as: 0 (0) exp( 2 ( ) ), ( ) where P i (z) is the power of mode i at distance z and L is the total length of the simulated tip [12].To calculate the propagation time T i of mode i in a SNOM tip we integrated the propagation constant as [12]: The optical path difference (OPD) values for a pair of modes ij can be calculated as follows: .

RESULTS AND DISCUSSION
To allow better comparison, the numerical simulation results obtained for SNOM tips with outer metal layer of Al and of Au are presented in separate subsections.

Simulation results obtained for Al-coated SNOM tips
At least 25 different EM modes are still able to propagate in the SNOM tip when the SiO 2 core radius R shrinks to 100 nm, but for most of these modes the transmittance T i is extremely low.As can be concluded from Fig. 3, TE and EH modes propagate mostly inside the silica core, because their attenuation decreases when the core radius R increases.TM and HE modes, on the contrary, evidently propagate on the metal surface, because at higher R values the R-dependence of their attenuation slows down and virtually stops becoming a constant.
It is clear that the modes with highest transmittance are the most relevant for comparison with experimental results.To substantially contribute to the MSM effect a pair of modes should also develop a large OPD while propagating in the SNOM tip.Based on the simulation results we can only roughly estimate the OPD, because our calculations are limited to the very end of the tip.However, as can be seen from Fig. 2, the region R  1.5 μm displays a very high dispersion of β, while already in the region R  5 μm (not shown in figures) β virtually does not depend on R.This indicates that the simulated end of a SNOM tip should play the most important role in developing significant OPD values.
The calculated OPD values between the six strongest modes are shown in Table 1.The observed MSM effect is unlikely to be affected by the interference of two strongest modes (TE01 and EH11) because this pair has a very small OPD (both modes have a similar β dependence on R), and in case of the smaller output aperture EH11 experiences a cutoff.On the other hand, both TE01 and EH11 have a large OPD with HE11, and these are the three most powerful modes.As discussed earlier, TE and EH modes propagate mainly in the core (thus being photonic modes), and HE modes -mainly on the metal surface (plasmonic modes).Thus we can conclude that the strong MSM contribution from the interference between photonic and plasmonic modes is indeed possible.However, taking into account the constantly high attenuation of plasmonic modes along the tip reveals an issue with such a simple explanation, because it makes highly improbable for these modes to propagate distances exceeding 50 μm.An alternative explanation is that the interference occurs between photonic modes with lower attenuation at larger R, e.g., between TE01 (or EH11) and EH12.In this case we have both low (or moderate) attenuation and reasonable OPD.

R (m)
However, we can also imagine a more complex mode excitation scheme involving more than two modes.E.g., the HE11 mode might get excited in the tip region with R  0.4 μm where it has the lowest attenuation of all the modes.Although our simulations do not reveal any interaction between modes, in actual experimental conditions surface roughness or SNOM tip imperfections could probably result in a kind of mode coupling.Bending of the tip could also have an effect on propagation and transformation of modes.Similarly, the TM01 mode might also get excited near the tip end where it has the second lowest attenuation.Calculating OPD theoretically, without having insights on the actual processes of mode transformation, does not seem to be a feasible task.

Simulation results obtained for Au-coated SNOM tips
In this case, for d = 200 nm the ten most powerful modes in the order of transmittance are: TE01, EH11, HE11, EH12, TM01, HE21, HE31, EH22, TE02 and HE41.These are the same ten modes as in the case of Al-coated tip, but in a bit different order.The dependences of propagation and attenuation constants for these modes on the core radius R are shown in   In the case of gold, the attenuation at the end of the tip is higher and modes experience slightly larger cut-off radiuses.The first three strongest modes are the same, but the following modes are mixed up.
The calculated OPD values between the six strongest modes are shown in Table 2.In general, the results are quite similar to the case of Al-coated tip.Again, one possibility to explain the observed MSM effect refers to the interference between plasmonic (HE11) and photonic (TE01 or EH11) modes, though HE11 has constantly high attenuation at larger R. Another possibility is the interference between photonic modes (e.g.TE01 and EH12).As already stated earlier, it is also possible that more than two modes participate, e.g. a photonic mode can excite the HE11 plasmonic mode in the region close to the very end of the tip where this mode has a low attenuation coefficient.

CONCLUSIONS
Although our numerical simulation of SNOM tip endings cannot be directly compared to the experimental OPD values, some conclusions can still be drawn.
First of all we conclude that, in accordance with earlier assumptions, the experimentally demonstrated enhanced levels of intermodal dispersion in metal-coated SNOM tips can be explained by assuming contribution to the observed MSM effect by the interference between photonic and plasmonic modes.However, this also assumes excitation of the plasmonic mode near the end of the tip and thus the involvement of more than two modes.One can speculate that the required mode coupling can be provided by SNOM tip imperfections or surface roughness, but can also be affected by tip bending.The experimental report on stronger intermodal dispersion in bent-type SNOM tips can be interpreted in favour of this mechanism.
On the other hand, an alternative explanation is possible, according to which the two modes contributing to the observed MSM effect are both of photonic origin.In this case the metal coating should not have any measurable effect on the intermodal dispersion in SNOM tips.Thus the reported experimental observation of negligible effect of replacing aluminium with gold in the SNOM tip metal coating seems to confirm this mechanism.
However, according to another conclusion of ours, for different plasmonic modes propagating in our considered SNOM tips the main role is played by a 20 nm thick chromium adhesion layer.Additional simulations of EM field structure (to be in more details presented elsewhere) indicate that typically only a subtle part of the EM energy reaches the outer metal layer, which also can explain its small influence on the intermodal dispersion.
Thus we can finally conclude that although our numerical simulations indicate two possible mechanisms explaining enhanced levels of intermodal dispersion in metal-coated SNOM tips, it is not possible to discern decisively between these two alternatives based on the existing set of experimental data.To overcome this problem, new experiments can be recommended with SNOM tips equipped with much thinner adhesion layer, so the plasmonic modes could fully interact with the outer metal coating.

Fig. 2 .
Fig. 2. Dependences of propagation constants on the SiO 2 core radius for ten strongest modes at the exit of an Al-coated SNOM tip.

Fig. 3 .
Fig. 3. Dependences of attenuation constants on the SiO 2 core radius for ten strongest modes at the exit of an Al-coated SNOM tip.
Figs 4 and 5, respectively.Like in Figs 2 and 3, dashed vertical lines represent output aperture positions for SNOM tips with d = 150 nm and d = 200 nm.At larger R, the dependences are similar to the case of Al-coated tip, though for R  0.2 μm some more pronounced differences appear.

Fig. 4 .
Fig. 4. Dependences of the propagation constants on the SiO 2 core radius for ten strongest modes at the exit of an Au-coated SNOM tip.

Fig. 5 .
Fig. 5. Dependences of the attenuation constants on the SiO 2 core radius for ten strongest modes at the exit of an Au-coated SNOM tip.

Table 1 .
OPD values (given in fs) between the 6 strongest modes in the simulated Al-coated SNOM tip (d = 200 nm) and the transmittance T i of the modes (last column).Highlighted mode combinations can most likely cause the observed MSM effect.

Table 2 .
OPD values (given in fs) between the 6 strongest modes in the simulated Au-coated SNOM tip (d = 200 nm) and the transmittance T i of the modes (last column).Highlighted mode combinations can most likely cause the observed MSM effect