Translator Disclaimer
19 December 2017 Electro-optical study of nanoscale Al-Si-truncated conical photodetector with subwavelength aperture
Author Affiliations +
A type of silicon photodiode has been designed and simulated to probe the optical near field and detect evanescent waves. These waves convey subwavelength resolution. This photodiode consists of a truncated conical shaped, silicon Schottky diode having a subwavelength aperture of 150 nm. Electrical and electro-optical simulations have been conducted. These results are promising toward the fabrication of a new generation of photodetector devices.



Near-field scanning optical microscopy (NSOM) has recently made significant progress in the manufacture of tips using processed methods,1,2 and the simulation of tip designs for future manufacture.3 Started three decades ago, by Lewis et al.,4 the technique enables us to overcome one of the fundamental limits of classical optics, the Abbe’s diffraction limit. The basic concept relies on detecting the evanescent field of light penetrating through a subwavelength aperture. Over time, a variety of useful applications have been developed, making near-field optics an exciting field of applied investigations: single molecule emission,5 spectroscopy,6 nanorecording, nanolithography, integrated circuits failure analysis, and biological inspection. To be sure, the key component of the NSOM remains the probes. Because of its high conversion efficiency relative to fiber probes, a semiconductor photodiode appears to be an attractive solution to detect the optical near field.7 In the past, pioneering works8 were performed on a micromachined pyramidal-shaped silicon Schottky diode with submicrometer aperture. The small lateral dimensions of this device required special care to prevent spurious photo-induced signals due to light absorption at device edges. Subsequently, improved tips were produced using the same kind of process; however, they were used for a scanning tunneling microscope purpose, rather than for an NSOM.9 To optimize the electro-optic performance of a silicon Schottky diode with submicrometer aperture, we developed a finite-element method (FEM) simulation using Comsol multiphysics. In this article, we present an analysis of the forecasted behaviors when running electro-optical simulations, as presented in the schematic flowchart of Fig. 1. Several Comsol modules, such as optics and semiconductor, are both coupled through the semiconductor and the optoelectronics interface. To enable a complete analysis of the behaviors, several optical input parameters (wavelength, aperture, and contact height) were varied; electrical input parameters (voltage, power, and concentration of dopants ND and NA) were varied as well to simulate the I–V curves.

Fig. 1

Schematic flow of the Comsol multiphysics used for the study of the photodetector device.


The simulation was performed with the Comsol multiphysics, which uses the FEM. The simulation couples the semiconductor and the wave optics modules. The device was modeled as a truncated cone. The half angle of the cone used was 33 deg. The radius of the upper active surface is 75 nm giving an active detector area of 17.7×103  μm2. The total height is 1.54  μm, and the lower radius is <1  μm.

The material chosen for the interior of the device was silicon, n-doped to a level Nd=1016  cm3 using analytic doping. Trap-assisted recombination was activated, with default values. The sides of the device serve as the Schottky junction. This is realized using the metal contact boundary condition with the rectifying junction option selected, together with “thermionic emission.” The material chosen for the model was aluminum, which sets the work function, and the external bias was set to the appropriate voltage. The Schottky junction extends only along the upper 1  μm of the device. The rest of the exterior is given an insulating boundary condition. This defines the boundary conditions relative to the semiconductor module.

To simulate photodetection, “optical transitions” was activated in the interior. This adds a photogeneration source term to the convection–diffusion equations, which is proportional to the intensity of the electromagnetic field. The precise form of this term is determined by Fermi’s golden rule, together with an empirical model for indirect absorption in silicon, based on an empirical model of Green and Keevers.10

The optical excitation is effected at the upper surface using a port boundary condition with the relevant value of the electromagnetic intensity. The lower surface is given an open boundary condition using the scattering boundary condition selection. Finally, with regard to the electric potential the upper surface is given an insulating boundary condition, while the lower surface is defined as a metal contact boundary condition.


Theory and Simulation Model


Evanescent Waves

In the work of Bethe,11 it was demonstrated that when a plane wave impinges on a flat surface with a round aperture of dimension significantly smaller than the wavelength, λ, the transmitted radiation field has the same form as if the aperture itself were a radiating dipole. In particular, for large distances from the aperture, the intensity falls off with inverse square of the distance to the center

Eq. (1)


Although Bethe’s solution addressed the inconsistencies of the scalar field diffraction by solving the full vector formalism of Maxwell’s equations, we will use the angular spectrum of the scalar theory.12 One assumes a propagating wave function ψ of wavelength λ satisfying Helmholtz’s (stationary) equation

Eq. (2)

(2+k2)ψ=0,where  k=2π/λ.

Time-independent plane wave solutions have the form

Eq. (3)


A general solution is a superposition of plane wave solutions

Eq. (4)


The constraint on k2=k2 allows us to solve one component of the wave number in terms of the other two, and write the whole expression as an integral over the remaining two components

Eq. (5)


Eq. (6)

defines the wave vector in the transverse direction associated with the transverse position vector r. The expression for the field Eq. (5) becomes

Eq. (7)

Note that the choice of sign + reflects forward propagation of the wave in the z-direction. Typically in the case of a propagating field, it is assumed that kx and ky are bounded by the condition k<k such that the root is real and the exponent term is complex. However if kk, the root is imaginary and the exponent term describes now a decay, i.e., the evanescent field. According to Eq. (7), for z=0, the propagating wave is

Eq. (8)

where ψ^(kx,ky) is the Fourier transform of the field over the plane z=0, defined by

Eq. (9)

If a flat image, a screen with an aperture is located at z=0, ψ(x,y,0) will describe the magnitude and phase of the field in space, and ψ^(kx,ky) its spatial frequency spectrum (k   space).

In light of the comments above, we arrive at the following connection between superresolution and evanescent waves: spatial frequencies greater than the wave number describe fine details with characteristic size smaller than the wavelength. These high-frequency elements of the spectrum will not excite propagating modes but rather will lead to evanescent, decaying modes. Subwavelength apertures, protrusions, etc. belong to the category of fine detail, which have high-frequency content. To see this explicitly, consider the canonical example of a circular aperture of radius a. For an incoming plane wave at normal incidence, the value of the field at z=0—just beyond the screen—is taken to be a constant value within the aperture, and zero outside

Eq. (10)


In each case the spatial frequency distribution is given by the Fourier transform, in accordance with Eq. (9). For the circular aperture, the Fourier transform of the circ-function is the well-known Airy function

Eq. (11)

where J1(x) is the Bessel function of first order.

These profiles are shown in Fig. 2 for two radius values: 0.6 and 0.2  μm.

Fig. 2

Spatial frequency profiles for circular aperture wave functions. Apertures radius: a1=0.6  μm, a2=0.2  μm. The corresponding first minima of the profiles are located at k1=z1/a1=6.4  μm1 and k2=z1/a1=19  μm1. Illumination wavelength is λ=0.5  μm (k=2π/λ=12  μm1). z1=1.22π is the first zero of the first Bessel function.


Both profiles are characterized by a central lobe or disc bounded by the first zero, and smaller side lobes of decaying amplitude. The location of the lobes is determined by the zeros of ψ^ function, i.e., for k1=z1/a1=6.4  μm1 and k2=z1/a2=19  μm1, where z1=1.22π is the first zero of the first Bessel function J1(z).

Considering now an incident wavelength λ=0.5  μm (k=2π/λ=12  μm1) such k1<k<k2. According to Fig. 2, the incident wave will pass through the 0.6-μm radius aperture as a propagating wave while it will pass through the 0.2-μm radius aperture as a decaying wave. Conversely, apertures of radius below z1/k=0.32  μm will lead to a spatial decay of the incident wave.

The intensity is proportional to the square of the field in Eq. (8); the total intensity is obtained by integrating over the plane z=0. By Parseval’s theorem, this is equal to the square of the integral of the Fourier transform distribution Eq. (11) over all spatial frequencies. For the circular aperture, the central peak contains 84% of the total intensity. In both cases, the overwhelming majority of spatial frequencies are concentrated within the central peak, and thus it is justified to focus on its location as characterizing the distribution of the frequencies. For the circular aperture, the central peak corner is defined by the first zero (z1) of the Bessel function

Eq. (12)

so most of the spatial frequencies k   lie within the cutoff limit k*

Eq. (13)


So, the smaller the aperture the larger is the cutoff frequency.

To understand the connection with evanescent waves, we consider the extent of the Fourier Transform relative to the cutoff frequency signaling the transition between propagating modes and decaying modes. When the width of the main peak increases beyond this cutoff, the frequency distribution will include decaying modes. Reciprocally, the more the aperture size decreases, the more the decaying modes are excited.

For our simulated photodetector, the radius a is 75 nm so the cutoff spatial frequency for the decaying mode behavior will be k*=z1/a=51.1  μm1, i.e., a cutoff wavelength λ* of 123 nm.


Subwavelength Imaging

We next present a brief discussion of the impact of evanescent waves on image resolution. Returning to Eq. (7), which relates the field at r=(r,z) to its spatial frequency components at z=0, we can define the factor

Eq. (14)

as the coherent optical transfer function. The coherent point-spread function is obtained by substituting Eq. (9) back (or simply by Fourier transforming with respect to rr)

Eq. (15)


Eq. (16)


The OTF Eq. (14) changes from a propagating phase to an exponential decay factor, which suppresses the frequencies for above k

Eq. (17)


The OTF thus acts as a low-pass filter directly in frequency space. The larger the distance z the more drastically it cuts off the high frequencies and the more the resulting image will be blurred.

For a subwavelength aperture of given radius a, it is possible to estimate the maximum distance at which detection is practical. Equivalently, for a given detection distance, it is possible to estimate the minimum-sized aperture detectable.

The starting point is Eq. (17) and the characterization of the frequency distribution by the central peak contained within the cutoff limit k*   is defined in Eq. (12). Suppression of frequency k is significant when

Eq. (18)

for some constant ϵ. Choosing suppression to a level ϵ=1 as indicative of signal loss is somewhat restrictive; instead we take as our criterion

Eq. (19)


Eq. (20)

For a given distance z, apertures with

Eq. (21)

will be detectable. Substituting for k* using Eq. (12), one obtains the following estimation for minimum detectable aperture size:

Eq. (22)


At large distances this reduces to the usual diffraction criterion

Eq. (23)


At small distances, however this becomes

Eq. (24)


So that resolution is bounded only by the minimum separation of the detector from the sample, and is independent of the wavelength.

For completeness, we include an expression for maximum estimated distance to detect an aperture of radius a

Eq. (25)

which reduces to Eq. (29) in the limit of small aperture size. The independence of the resolution from the wavelength in the limit of close proximity, which is predicted by Eq. (28) may seem surprising; it should be less so when one considers that when z is close to 0, the expression for the PSF in Eq. (16) approaches a delta function.

Stated otherwise, the notion of wavelength is fundamentally a property of propagating, oscillating fields; evanescent modes are closer in nature to electrostatic fields. By measuring the sample at increasingly smaller distances, the exponential suppression can be mitigated, and the detector approaches sampling the precise field profile of the source.


Electrical Device Modeling

The device couples a variety of physical phenomena to achieve detection of incoming radiation: the incoming electromagnetic wave is described by Maxwell’s equations, the electronic transport and charge distribution within the silicon is determined by the drift-diffusion equations for electrons and holes in a semiconductor coupled to Poisson’s equation, and these are augmented by statistical mechanics to describe the effect of doping on the carrier charge population. The interaction between the radiation and the silicon—photo-detection—involves indirect and direct interband transitions. This interaction occurs primarily within the Schottky junction that develops between the silicon bulk and the aluminum surface; the electronic properties of this junction reflect primarily the phenomenon of thermionic emission between the two materials. We will present here the differential equations in the bulk, and the boundary conditions, which are applied to the numerical algorithm to simulate these phenomena.

When a metal is brought in contact with an extrinsic semiconductor, a similar situation results, provided that the Fermi level of the metal is lower than that of the semiconductor, for n-type doping (in our case study). Electrons then diffuse from the semiconductor to the metal creating a depletion layer. In the simulation, the Schottky barrier junction between the aluminum and the silicon is realized as a boundary condition on the potential. This is possible because for a junction between a metal and a nondegenerately doped semiconductor the built-in potential falls, to good approximation, entirely on the side of the semiconductor. The metal, with its high density of states, can provide a space charge of almost arbitrarily high density. A very thin layer of charge is thus sufficient to neutralize an impinging field. The potential drop, which is proportional to the width, is thus much smaller on the metal side. At a junction between a metal and nondegenerately doped semiconductor, the metal acts essentially as an extremely highly doped semiconductor. For n-type-doped silicon, for instance, the metal acts as if it were p-type, accepting electrons and accumulating them in a thin surface layer as described above. We can thus make the approximation that the potential ϕ takes the value it has in the bulk of the metal, on the boundary of the semiconductor, and falls to the bulk semiconductor value entirely within the semiconductor. As a quantitative analysis, one starts with the decomposition of the chemical potential into an internal component and an external component

Eq. (26)


The “external” component is an electrostatic contribution due to space charge accumulation, in this case the depletion layer

Eq. (27)

where ϕ   is the semiconductor potential at the surface contact and measured as the band bending.

The “internal” component is determined by the band structure—in particular the conduction band edge Ec and the valence band edge Ev—together with the Fermi level. For a nondegenerated semiconductor, for instance, the electrons are well described as an ideal gas obeying Boltzmann statistics; the latter contribution is the diffusion potential. This gives

Eq. (28)

where μint is in fact the position of semiconductor Fermi level relative to the vacuum level, i.e., the semiconductor work function ϕs. n is the electron density and Nc is the density of states of the conduction band. This is measured relative to the ground state of a free electron in the vacuum. The sign signifies the fact that the value of the Fermi level lies beneath the zero energy-vacuum level. For a metal, the chemical potential is reduced to its internal component, i.e., the work function ϕm. Assuming the metal-semiconductor contact at equilibrium, we have

Eq. (29)


Eq. (30)


Identifying the difference of the total chemical potential with the external bias, μtotal,mμtotal,s=qV gives the following boundary condition for the semiconductor potential ϕ:

Eq. (31)


One should note that in Comsol, the electrical potential is defined relatively to the metal vacuum level. In the simulation, we thus used the relation

Eq. (32)


The Schottky-junction is a heterojunction; there is a discontinuity in the band structure. As described above, the formation of a depletion layer leads to band bending; a potential barrier between the semiconductor and the metal is the result. Ideally, for an n-type semiconductor, the energy barrier height is defined13

Eq. (33)

where χ is the electron affinity of the semiconductor and qϕm is the metal work function. Assuming an n-type silicon and χ=4.05  eV in contact with aluminum, and by taking qϕm=4.72  eV, we obtain an effective barrier height13

Eq. (34)

qϕBn=0.67  eV.

The electronic behavior can be analyzed using the model of thermionic emission from the semiconductor to the metal. In this model, the current is found by using Boltzmann statistics to compute the average velocity as a thermal sum over all energies higher than the barrier. The barrier changes with the external bias V and thus the velocity and hence the current I depends on the external bias. The result is given by

Eq. (35)

where A is the cross-section area of the junction. The saturation current density is given by

Eq. (36)

where ϕBn is the height of the zero-bias potential barrier. The constant C* is

Eq. (37)

with m* the effective mass of the charge carriers in the semiconductor. For free electrons, C* value is found to be the Richardson’s constant

Eq. (38)

For n-type silicon of crystallographic type 111, the value of the ratio is C*/C=2.2.13


Electro-Optical Results and Interpretation


Photodetector Structure Simulation

The NSOM photodetector device is simulated as a truncated conical structure of 1.54-μm height, 1-μm bottom radius, and a 75-nm top radius (Fig. 3). The material is defined as n-type silicon crystal. Since Comsol software is based on finite elements calculation, a high density mesh of vertices is necessary to define the shape with enough resolution. To form a Schottky diode, a 1-μm-high surface of the structure from the top is defined as metal (Al) contact, where the bias voltage is applied respective to the grounded bottom.

Fig. 3

NSOM 3-D truncated conical-shaped photodetector device structure. The main parameters are: height is 1.54  μm, top radius is 75 nm, and bottom radius is <1  μm.


In the simulated device, the contact area defined as the truncated cone surface is

Eq. (39)

A=2.741  μm2.


Preliminary Electrical Characterization

As a first attempt to check the physical consistency of the device simulation as a Schottky diode, the I–V curve is simulated without illumination at 300 K (Fig. 4). The current range is limited to 500 nA to keep a reasonable current density (lower than 10  A/cm2). An exponential fit perfectly matches the I–V curve with an exponential factor of 27  V1. This value is different from the calculated value from Eq. (35), i.e., 38.5  V1 and may be related to the nonplanar geometry of the device.

Fig. 4

Simulated I–V curve without illumination at 300 K and exponential fit.


Then, we simulated the distribution of the electron current density in Fig. 5 at forward bias of 0.25 V. The scale of the current density is coherent with current values as shown in Fig. 4.

Fig. 5

Electron current density (A/m2) without illumination for a forward bias of 0.25 V.


Since a diode photodetector is generally biased at reverse voltage, we investigated the distribution of the minority carriers (holes) concentration at 0.5  V relative to the grounded bottom electrode, without illumination, as shown in Fig. 6. As observed there, the holes are accumulated at the top and at the 1-μm-height surface contact of the device. At the surface, max of the log concentration is 11.

Fig. 6

Simulated distribution of the holes concentration without illumination for a reverse bias of 0.5  V.


Further analysis was conducted for a reverse bias of 0.5  V using several electrical simulations: distribution of the electron concentration showing the expected depletion region below the surface (Fig. 7) and the silicon electric potential ϕ (Fig. 8) is defined by Eq. (32).

Fig. 7

Simulated distribution of the electron concentration showing the extension of the depletion area for a reverse bias of 0.5  V without illumination.


Fig. 8

Simulated distribution of the electric potential for a reverse bias of 0.5  V. Metal work function for the Schottky contact is 4.72 eV.



Influence of the Illumination

Simulation of light illumination applied at the top aperture of the truncated conical shape is shown in Fig. 9. A wavelength λ of 550 nm was used to study the simulated distribution of and the minority carrier’s (holes) concentration under illumination power of 0.1 W as shown in Fig. 10. Also, in Fig. 10 we can see the small increases in the holes (minority carriers) concentration due to illumination. The max log of concentration is 11.5 at the top, when compared with 11 without illumination (Fig. 6).

Fig. 9

Illumination coupled to the top of the truncated conical shape. Wavelength is 550 nm.


Fig. 10

Distribution of holes concentration showing changes in depletion size area. Bias voltage 0.5  V under illumination power 0.1 W and λ=550  nm.



Evanescent Decay of the Radiation Field

As discussed above, radiation transmitted through a subwavelength aperture is expected to decay with distance depending on the distance z relative to the aperture radius a. Indeed, at large distance (za), the field behaves like a classic propagating wave and decays as a square law expressed in Eq. (1). But, at short distance za the radiation will decay as an exponential law as mentioned by the OTF expression [Eq. (17)]. In addition, optical absorption by the silicon will also lead to an exponential decay of the intensity within the device, according to the Beer–Lambert rule. Figure 11 shows the behavior of the electric field as a function of the distance z along the central axis of the device for several wavelengths (50 and 300 nm) in air and silicon media. It is observed that the 50 nm does not present decay in air medium, since it is not influenced from the diffraction phenomena from the 150-nm aperture. However, it presents a clear decay in the silicon medium due to Beer–Lambert law (for such wavelengths the absorption is very small, about 10 nm). For 300 nm, there is a similar decay in air and silicon, mainly due to the diffraction from the aperture of 150 nm (evanescent wave), so Beer–Lambert law is less influencing. The difference in the electrical field amplitudes between the wavelengths seems to be caused by the difference in the energy, since the energy is a function of the wavelength (1/λ), and of the squared electrical field. Under such conditions, the 50 nm has bigger amplitude than the 300 nm in the air.

Fig. 11

Decay of the incident electric field for two different wavelengths: 50 and 300 nm, respectively, below and above the diffraction limit (for a circular aperture of 75 nm radius).




We have presented a type of silicon photodiode, which has been designed and simulated to probe the optical near field. This photodiode consists of a submicron aluminum/silicon Schottky diode (150-nm diameter aperture) on top of a truncated conical shape. Electrical and electro-optical simulations have been conducted. These results present a promising device as a function of several parameters such as wavelengths. Since the 75-nm radius of the circular aperture is still far from the well-known 10-nm quantum range, quantum effects have been ignored in this analysis, and the focus was on NSOM. Still further investigations, such as smaller varying sizes of the aperture, are needed to optimize the performance of such a challenging nanoscale photodetector.



F. Lu, M. Jin and M. A. Belkin, “Tip-enhanced infrared nanospectroscopy via molecular expansion force detection,” Nat. Photonics, 8 (4), 307 –312 (2014). NPAHBY 1749-4885 Google Scholar


L. Stern et al., “Near field phase mapping exploiting intrinsic oscillations of aperture NSOM probe,” Opt. Express, 19 (13), 12014 –12020 (2011). OPEXFF 1094-4087 Google Scholar


S. Arora et al., “Design of MEMS based microcantilever using Comsol multiphysics,” Int. J. Appl. Eng. Res., 7 (11), 2012 (2012). Google Scholar


A. Lewis et al., “Development of a 500 Å spatial resolution light microscope. Part I. Light is efficiently transmitted through lambda/16 diameter apertures,” Ultramicroscopy, 13 (3), 227 –231 (1984). ULTRD6 0304-3991 Google Scholar


R. X. Bian et al., “Single molecule emission characteristics in near-field microscopy,” Phys. Rev. Lett., 75 (26), 4772 –4775 (1995). PRLTAO 0031-9007 Google Scholar


J. K. Trautman et al., “Near-field spectroscopy of single molecules at room temperature,” Nature, 369 (6475), 40 –42 (1994). Google Scholar


G. Kolb, K. Karrai and G. Abstreiter, “Optical photodetector for nearfield optics,” Appl. Phys. Lett., 65 (24), 3090 –3092 (1994). APPLAB 0003-6951 Google Scholar


R. A. Chelly et al., “Pyramid-shaped silicon photodetector with subwavelength aperture,” IEEE Trans. Electron Devices, 49 (6), 986 –990 (2002). IETDAI 0018-9383 Google Scholar


I. J. Vera Marún and R. Jansen, “Multiterminal semiconductor/ferromagnet probes for spin-filter scanning tunneling microscopy,” J. Appl. Phys., 105 (7), 07D520 (2009). JAPIAU 0021-8979 Google Scholar


M. A. Green and M. J. Keevers, “Optical properties of intrinsic silicon at 300 K,” Prog. Photovoltaics Res. Appl., 3 189 –192 (1995). PPHOED 1062-7995 Google Scholar


H. A. Bethe, “Theory of diffraction by small holes,” Phys. Rev., 66 (7–8), 163 –182 (1944). Google Scholar


J. W. Goodman, Introduction to Fourier Optics, 2nd Ed.McGraw-Hill, New York (1996). Google Scholar


S. M. Sze and K. K. Ng, Physics of Semiconductor Devices, 3rd ed.John Wiley and Sons, Inc., Hoboken, New Jersey (2007). Google Scholar


Matityahu Karelits is pursuing his BSc degree in applied physics /electro-optics engineering at the Lev Academic Center (formerly JCT—Jerusalem College of Technology).

Yaakov Mandelbaum studied mathematics and physics at the undergraduate and graduate levels at the University of Pennsylvania, MIT, and the Hebrew University of Jerusalem. After working in industry as a physicist and electro-optics team leader, he joined the Lev Academic Center in 2013 as a lecturer and researcher in the Department of Applied Physics and Electro-Optics. In parallel, he is currently pursuing his PhD at Bar-Ilan University.

Avraham Chelly received his PhD in solid state physics from the Université de Haute Alsace, Mulhouse, France, in 1997. During his postdoc at the Hebrew University of Jerusalem, Israel, he was involved in the development and fabrication of a silicon nano-detector for near-field optics. In 2004, he moved to Bar-Ilan University where he established the Advanced Semiconductor Devices Laboratory. He is giving lectures on semiconductor processing and devices.

Avi Karsenty received his PhD in applied physics/material science (microelectronics/electro-optics) from the Hebrew University of Jerusalem in 2003. His research focuses on nanoscale electro-optics coupled-devices. After 22 years in high-tech industries, part of which as an engineer and a manager for 16 years with Intel, he is today the physics/electro-optics engineering department’s head. He is IEEE senior member and OSA senior member, and has received 35 awards in engineering/physics.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Matityahu Karelits, Ya'akov M. Mandelbaum, Avraham R. Chelly, and Avi Karsenty "Electro-optical study of nanoscale Al-Si-truncated conical photodetector with subwavelength aperture," Journal of Nanophotonics 11(4), 046021 (19 December 2017).
Received: 4 October 2017; Accepted: 28 November 2017; Published: 19 December 2017


Back to Top