Broadband Tunable Infrared Light Emission from Metal-Oxide-Semiconductor Tunnel Junctions in Silicon Photonics

Broadband near-infrared light emitting tunnel junctions are demonstrated with efficient coupling to a silicon photonic waveguide. The metal oxide semiconductor devices show long hybrid photonic–plasmonic mode propagation lengths of approximately 10 μm and thus can be integrated into an overcoupled resonant cavity with quality factor Q ≈ 49, allowing for tens of picowatt near-infrared light emission coupled directly into a waveguide. The electron inelastic tunneling transition rate and the cavity mode density are modeled, and the transverse magnetic (TM) hybrid mode excitation rate is derived. The results coincide well with polarization resolved experiments. Additionally, current-stressed devices are shown to emit unpolarized light due to radiative recombination inside the silicon electrode.


Supplementary information Fabrication
We fabricated the tunnel devices from a 340 nm thick silicon-on-insulator (SOI), which in a ≈ first step were globally doped to approximately via a two-step ion  phosphorus ≈ (5 ⋅ 10 18 ) cm -3 implantation and activation and diffusion annealing.Then, using e-beam lithography, the silicon was fully etched with an inductively coupled plasma etch (ICP) to electrically separate the different devices and to form the photonic waveguides.A following partial etch formed the thin ridges for the silicon contact.Next, using optical lithography resist and a lift-off, nickel was deposited via electron beam evaporation onto the silicon and subsequently annealed in forming gas to form a silicide contact.After that, a conformal SiO 2 insulating layer was deposited using plasma-enhanced atomic layer deposition (PE-ALD), which was then locally opened down to the silicon using electron-beam (e-beam) lithography and reactive ion etching (RIE), in order to define the tunnel junction areas.A very short (2 s) dip in buffered hydrofluoric acid (1:7 solution) was then used to remove the residual oxide over the tunnel area (ALD and native oxide) before the samples were directly transferred to the PE-ALD to deposit the tunnel barriers.Since an oxygen plasma was used in the PE-ALD which oxidizes the silicon surface, there is always an offset of tunnel oxide thickness, as illustrated in Supplement Figure 1.The deposited thicknesses were measured via ellipsometry on a pure silicon reference sample, and the final thickness used was .The last step was  SiO2 = 2 nm then the e-beam lithography and e-beam evaporation of 1 nm Titanium and 100 nm ≈ ≈ Gold to define the second electrode.

Cutback Measurement
To characterize the propagation losses of the hybrid plasmonic mode, we performed a cutback measurement, where several devices with different hybrid propagation lengths were measured.After setup and grating coupler loss subtraction, the resulting transmission losses are averaged from 1480 nm to 1560 nm to minimize resonant effects.In addition to the propagation length, the coupling efficiency between the photonic and hybrid WG can also be extracted as %, as the intersection of the fitted line at zero hybrid length   = 74 +6 -8 corresponds to twice the coupling loss.

Emission Measurement Setup and Calibration
The emission is coupled to a silicon photonic waveguide and guided to the edge of the diced and mechanically polished chip, where it is collected by a high numerical aperture (NA) NIR microscope objective (50x, NA=0.65).A c-coated linear polarizer filters the emission and an anti-reflective coated objective lens then focuses the light to a Czerny-Turner spectrometer with a 45 lines per millimeter, 1.75 μm blazed grating.There, a thermoelectrically cooled InGaAs focal plane array is used for alignment and power measurements, and a liquid nitrogen cooled linear InGaAs detector is used to record the spectra.
Supplement Figure 3: Sketch of the used emission measurement setup, where either the Nirvana or Pylon detector was used for power or spectral measurements respectively.
To measure the emitted power via the InGaAs focal plane array, the polished edge facet is focused on to the detector with a mirror installed in the spectrometer's grating turret.Each measurement takes one emission image and one background image for dark subtraction.The recorded pixel-intensities are summed up over a limited area where the emission is recorded.This summed intensity corresponds to the total number of photons being detected due to the detector's unity gain.The same is done on a second area on the signal image, right below the emitting area.From this second dark area, the standard deviation of the pixel values is calculated to estimate the dark current noise in the measurement, which is shown in the error bars in the main text.
To estimate the actual emitted power into the waveguide as well as the source efficiency, the following values for the losses are considered, where the dipole-photon coupling factor is only used in the efficiency estimation.

Efficiency factor Value Comments
Dipole excitation to photon The efficiency values in the main text thus consist of the vacuum source efficiency ( ) ≈ 10 -6 times the transversal, i.e. heterostructure, LDOS enhancement factor ( ) ≈ 300 For the spectral measurements, the spectrometer's wavelength was calibrated using a Hg/Ne calibration source.The relative system response calibration was performed with a known Halogen lamp, which was beforehand measured on a NIST traceable NIR spectrometer.

1D mode enhancement
The 1-dimensional mode density enhancement is calculated by the angular spectrum representation of the dipole's emitted field.We calculate the total field at the dipoles position by applying the appropriate boundary conditions from the layer stacks permittivities 2 .The method is based upon and explained in detail in work from Parzefall et.al. 3 We calculate the optical mode density enhancement of the heterostructure via: )

Influence of Silicon Thickness
While the tunnel oxide thickness is very limited in variability to maintain direct tunnelling, we could potentially vary the thickness of the silicon device layer.Supplement Figure 5 shows the calculations of the emission efficiency into the fundamental TM0 mode as well as the total emission into the TM0 for a heterostructure with varying silicon thickness.
Supplement Figure 5: Left: For the heterostructure Au-SiO2-Si the emission efficiency into the TM0 mode is plotted for varying silicon thickness and wavelengths.Right: For the same structure and variations, the total emission into the TM0 is plotted, instead of the efficiency as in the left plot.

Transverse multimodes
To better understand the influence of the width of the tunnel section, we perform 3D FEM simulations with a single dipole emitter in the center of the thickness of the oxide of a hybrid waveguide at 1400 nm.The resulting optical field in the hybrid mode WG is overlapped with the guided eigenmodes, showcasing the excitation of the different modes for varying waveguide widths.The overall emitted power for a single dipole emitter does not change, however in the transverse multi-mode case, the emitted power is distributed relatively evenly between the different modes.The coupling into the respective mode is calculated as Where is the field overlap and the Poynting vector component along x, y, or z.

𝜒 𝑃 𝑥,𝑦,𝑧
As is evident in Supplement Figure 6b, with increasing width the number of transverse modes increases very quickly.As a result, a wide (25 μm), highly transversally multi-moded LETJ can support around 120 different modes, all with effective refractive indices between 1.5 and 3.5.In addition, the lateral dipole position influences which modes can be excited, as illustrated in Supplement Figure 6c, where only the even transverse modes are excited due to the central position of the simulatied excitation-dipole.The effect of non-homogenous tunneling path distribution should be more pronounced the smaller the cavity area becomes, as fewer and fewer tunnel paths will contribute to the emission.Reinforcing this hypothesis is the recorded spectrum of a transversally very large tunnel device, illustrated in Supplement Figure 7. There, we assume that in total more tunnel paths are actively emitting and thus effectively more excitation dipoles are distributed inside the cavity.This leads to a better agreement in the cavity mode intensities with the simulated model where evenly distributed dipoles are assumed.However, the wide waveguide structure now supports more than 120 transversally higher order modes, all of which will be excited.As we only collect a significant amount of light from the first few modes, all having similar effective refractive indices -and thus similar longitudinal resonance frequencies -we still record the cavity resonances, albeit with reduced quality factors due to higher order mode averaging.

Cavity Model
To accurately model the resonances of the measured emission spectrum, we consider the emitter inside a Fabry-Pérot cavity with frequency-dependent mode index, loss, and reflectivities.
We assume that the inelastic tunneling processes occur incoherently and can be regarded as spontaneously emitting dipole sources distributed along the cavity's length .The  spontaneous emission enhancement due to the increased density of optical states can be calculated with the field enhancement at the location of the dipole. 4As we have many dipolelike emitters throughout the cavity, we need to average the field enhancement inside the cavity over the entire length to model the total emission. 5pplement Figure 8: Sketch of the Fabry-Pérot cavity model used.
Inside the cavity, we have a hybrid plasmonic mode with an effective refractive index and  eff resulting k-vector , where is the radial frequency and the vacuum speed of light.The emitted field leaving the cavity can be calculated as  0  t the sum of all the fields inside, originating from a dipole at position To get the total  0 .emission, we integrate the field enhancement over all different dipole positions 0 <   < .
We define a propagator and the mirror reflection and transmission () =  - ′   - ′′  coefficients and respectively.Forming a geometric series of the reflected fields inside  1,2  1,2 the cavity, we can define the coefficient And considering the propagation from , we can write out the emitted field as The total field enhancement than becomes: Since the reflectors are strongly frequency dependent, they are simulated with 3D FDTD, where the geometry is modelled as close to the fabricated device as possible.

Tunneling calculations
Both elastic and inelastic tunneling are treated as a first order perturbation on the Hamiltonian of the electrode/barrier system 6 .This assumption is valid as long as the electrons in both electrodes are only weakly coupled, which is the case for the barrier thicknesses considered in this work 7,8 .The total system Hamiltonian can thus be written as 2,9  =  0 +  el +  inel Where we have and as the elastic and inelastic perturbative terms.    Supplement Figure 9: left: the case of the elastic tunneling rate calculation is sketched.Right: the inelastic tunneling rate calculation considers all transition where the electrons in the left and right electrode have energies differing by .ℎ We now want to determine the wavefunctions separately for both electrodes. L,R

Elastic Tunneling Rate
Due to the weak coupling, we model each electrode's potential separately: , Where is the barrier thickness.As a simplification, we assume the barrier to be rectangular. A more detailed calculation could be made with a discretization of the barrier potential 10 , which would however require a fully numerical calculation of the transmission rates.
With rectangular potentials, the electrons' wavenumbers become ℏ Where and are the electron's wavenumber in the left electrode and the left barrier is the effective electron mass inside the barrier, which for SiO 2 is assumed  b as 11 .is the left electrode's barrier height.We apply the bias voltage  b,SiO2,parabolic = 0.3  0  L on the right electrode and get for its wavenumber with the assumed rectangular barrier ℏ Following Bardeen's tunneling theory 6,7 , we get the elastic transition matrix element () Where are the unperturbed wavefunctions on the left and right side ( . L,R →∞) For our rectangular barrier system, we can solve the time-independent Schrödinger equation to get the following wavefunctions Inserting (Error!Reference source not found.)into (Error!Reference source not found.)leaves us with the elastic tunnelling transition matrix element for a rectangular barrier By integrating over all energies and including the electrodes, respective density of states  L/R , we end up with the total elastic electron transition rate .

Inelastic tunnelling rate
For our considered case of weak coupling, we can treat the inelastic tunnelling transition in the same way as the elastic transition as a perturbation to the system Hamiltonian.Here we use the Light-Matter interaction Hamiltonian With as the vector potential operator and as the momentum operator.The   = -ℏ / inelastic transition matrix element is the expectation value when applying the (,ℎ) momentum operator: We assume that the vector potential inside the rectangular tunnel barrier is varying much slower than the electron wavefunction, which allows us to evaluate this transition matrix element at any point inside the barrier.This gives us the evaluated element ′′  - ′′  Following Fermi's golden rule we can then calculate the total spectral inelastic tunneling rate by integrating over all allowed energies for the photon energy is the vacuum permittivity, the electron charge  0 We assume the gold to have a constant electronic density of states (eDOS) over the considered electron energies 3 , while we take a square-root dependence of the eDOS at the band-edges of silicon due to the parabolic dispersion of electrons.
Parameters for the tunnel rate calculations:

Emission at different bias voltages
When measuring the TM tunnel emission under different voltages, shown in Supplement Figure 12, no discernible shift in the spectra is visible, apart from an overall stronger emission.When looking at the ratios of the cavity peak intensities no clear trend is visible, as is expected from the flat vacuum source spectrum shown in Figure 2(a) in the main text, where for biases above 1.25 V no significant change in the vacuum source spectrum is visible in our system.

Electrical characteristics of the tunnel junctions
In Supplement Figure 13 an exemplary IV curve of a MOS tunnel junction is shown.We can also calculate the theoretical direct tunneling current following a Wentzel-Kramers-Brillouin approximation 12 :

Supplement Figure 1 :
left: The fabrication steps are illustrated.right: Ellipsometrically measured tunnel oxide thickness on silicon, dependent on cycle numbers in the PE-ALD at 200 °C.With the use of an oxygen plasma in the growth, there is always some offset of tunnel oxide thickness, i.e. the tunnel oxide consists of both plasma-grown and ALD-grown oxide.
prop = 9.6 +1.9 -1.4 μm Supplement Figure2: left: an SEM image of the cutback measurement structure is shown.right: the extracted device losses for various hybrid mode section lengths are shown, including the linear fit in the semi logarithmic plot to extract the losses.The loss uncertainties are taken from the uncertainty of the fit and signify the 95% confidence bounds.
where is the vacuum wavevector, is the in plane wavevector component, and is the  0  ||   out of plane wavevector component in material (i.e.SiO 2 ).and are the upwards and   ↑   ↓  downwards traveling complex e-field amplitudes, and is the permittivity of material . Applied to our MOS structure, we can calculate the dipole emission enhancement as shown in Supplement Figure4, where the TM0 and TM1 modes are visible as resonances.The integrated area under the TMx mode then corresponds to the value of in the main text,    0 which is mentioned to be around of the total emission.25 % Supplement Figure4: Calculated dipole radiation enhancement purely due to the MOS heterostructure.The normalized wavenumber is used, and replicates the of the supported modes.The  =  || / 0   strong plasmonic resonance for frequencies above 550 THz is clearly visible as a broadening of the TM0 mode's resonance and its increase in wavenumber.

Supplement Figure 6
: a: The transverse fundamental mode TM00 and the first order mode TM01 are illustrated.b: Resulting from 2D eigenmode analysis, the real part of the effective mode indices at vacuum wavelength of 1400 nm are plotted versus an increasing waveguide width.c: Originating from a central dipole excitation, coupling into the various modes is shown for varying WG width.

Supplement Figure 7 :
(a): Identical sub-figure from the main text, showcasing the different contributions to the simulated emission spectrum on the right.(b): colorized SEM image of wide, transversally multimoded device, whose TM spectrum is shown in (c).The same longitudinal cavity modes as for the single-moded device can be observed from the emission under a bias voltage .  = 2.5  (  ≈ (7.04 ± 0.08)   2 )

Supplement Figure 12 :
Left: The TM emission of one device under different bias voltages is shown.1-5 label the peak numbers used in the right plot, and the gray areas show the averaging windows used to calculate the peak intensity.Right: The peak intensity ratios between the peaks in the left plot are shown, with the error bars signify one standard deviation of the noise over the averaging window.

⋅ 1 ℎ)
(  B - ox 2 ) ⋅ exp ( -4 2 * e  SiO2 ⋅ (  B - ox 2 ) ⋅ 1 ℎ ) ) -(  2 2ℎ  SiO2 ⋅  B )  B exp(-4 2 * e  SiO2  B ⋅ With being the electron charge, the Planck constant, the tunnel barrier thickness,  ℎ  SiO2  B the tunnel barrier height, the voltage over the oxide, and the reduced electron mass  ox  * e inside the tunnel barrier.The oxide voltages for the tunneling model are calculated by solving the Poisson equation.A reasonable fit of the tunnel currents through the thin oxide can be made for barrier 2 nm potentials between electrode and tunnel barrier of for positive biases and   = 1.45 eV  b = 2 for negative biases.The difference in barrier potentials can be explained due to the eV tunneling occurring either from the silicon conduction band into the gold (positive biases) or from the gold with a large work function into the silicon (negative biases).Supplement Figure 13: Measured and calculated direct tunnel currents are shown.