††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# High-harmonic generation enhancement with graphene heterostructures

I. Alonso Calafell  University of Vienna, Faculty of Physics, Vienna Center
for Quantum Science and Technology (VCQ), 1090 Vienna, Austria  L. A. Rozema
University of Vienna, Faculty of Physics, Vienna Center for Quantum Science
and Technology (VCQ), 1090 Vienna, Austria  A. Trenti  University of Vienna,
Faculty of Physics, Vienna Center for Quantum Science and Technology (VCQ),
1090 Vienna, Austria  AIT Austrian Institute of Technology, Center for Digital
Safety & Security, Vienna, Austria  J. Bohn  Department of Physics and
Astronomy, University of Exeter, Exeter, EX4 4QL, UK  E. J. C. Dias  ICFO-
Institut de Ciencies Fotoniques, The Barcelona Institute of Science and
Technology, 08860 Castelldefels (Barcelona), Spain  P. K. Jenke  University of
Vienna, Faculty of Physics, Vienna Center for Quantum Science and Technology
(VCQ), 1090 Vienna, Austria  K.S. Menghrajani  Department of Physics and
Astronomy, University of Exeter, Exeter, EX4 4QL, UK  D. Alcaraz Iranzo  ICFO-
Institut de Ciencies Fotoniques, The Barcelona Institute of Science and
Technology, 08860 Castelldefels (Barcelona), Spain  F. J. García de Abajo
ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of Science and
Technology, 08860 Castelldefels (Barcelona), Spain  ICREA-Institució Catalana
de Recerca i Estudis Avançats, Passeig Lluís Companys 23, 08010 Barcelona,
Spain  F. H. L. Koppens  ICFO-Institut de Ciencies Fotoniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
ICREA-Institució Catalana de Recerca i Estudis Avançats, Passeig Lluís
Companys 23, 08010 Barcelona, Spain  E. Hendry  Department of Physics and
Astronomy, University of Exeter, Exeter, EX4 4QL, UK  P. Walther  University
of Vienna, Faculty of Physics, Vienna Center for Quantum Science and
Technology (VCQ), 1090 Vienna, Austria

(August 27, 2024)

###### Abstract

We investigate high-harmonic generation in graphene heterostructures
consisting of metallic nanoribbons separated from a graphene sheet by either a
few-nanometer layer of aluminum oxide or an atomic monolayer of hexagonal
boron nitride. The nanoribbons amplify the near-field at the graphene layer
relative to the externally applied pumping, thus allowing us to observe third-
and fifth-harmonic generation in the carbon monolayer at modest pump powers in
the mid-infrared. We study the dependence of the nonlinear signals on the
ribbon width and spacer thickness, as well as pump power and polarization, and
demonstrate enhancement factors relative to bare graphene reaching 1600 and
4100 for third- and fifth-harmonic generation, respectively. Our work supports
the use of graphene heterostructures to selectively enhance specific nonlinear
processes of interest, an essential capability for the design of nanoscale
nonlinear devices.

High-harmonic generation (HHG) has been intensely investigated as a route
towards the generation of coherent attosecond radiation in the extreme
ultraviolet and x-ray spectral regions. Atomic gases have thus far been the
most successful among demonstrated systems for HHG McPherson  _et al._ (1987);
Ferray  _et al._ (1988), although they require high vacuum, making them
impractical for the design of integrated devices. As a result, the development
of solid-state HHG systems has become an important challenge, with promising
results demonstrated for various crystalline materials Ghimire  _et al._
(2011); Schubert  _et al._ (2014); Luu  _et al._ (2015). However, the
mechanisms of HHG in solid state appear to be fundamentally different from
that in atomic gases, and are highly sensitive to crystal and polarisation
orientations Langer  _et al._ (2016); You  _et al._ (2017), as well as to the
fundamental optical properties of the materials. In particular, theoretical
models proposed for solid-state HHG have identified both interband transitions
Langer  _et al._ (2016) and intraband electron dynamics You  _et al._ (2017)
as radically different sources for the required anharmonicity.

Graphene constitutes an appealing choice of material in this context because
it features sizeable intraband and interband contributions to the linear Jiang
_et al._ (2018); Soavi  _et al._ (2018); Mikhailov (2016); Rostami  _et al._
(2017) and nonlinear Mikhailov and Ziegler (2008); Marini  _et al._ (2017)
conductivity in the infrared (IR). Moreover, as a two-dimensional material, it
simplifies the conditions needed to achieve phase matching of the signal
generated at different spatial locations (e.g., this is automatically
guaranteed for pumping at normally incidence), which commonly limits the
strength of nonlinear processes in thick crystals. Third- Jiang  _et al._
(2018); Soavi  _et al._ (2018); Kumar  _et al._ (2013); Hong  _et al._ (2013);
Hendry  _et al._ (2010); Baudisch  _et al._ (2018); Beckerleg  _et al._ (2018)
and second-order Constant  _et al._ (2016) nonlinearities have already been
demonstrated in graphene, while higher-order harmonic generation has also been
observed Bowlan  _et al._ (2014); Hafez  _et al._ (2018); Kovalev  _et al._
(2021); Baudisch  _et al._ (2018); Yoshikawa  _et al._ (2017). These results
have stimulated proposals for the exploitation of the nonlinear response of
graphene to implement devices for quantum technology Koppens  _et al._ (2011);
Gullans  _et al._ (2013); Jablan and Chang (2015); Alonso Calafell  _et al._
(2019).

The conical electronic band structure of graphene has been argued to boost the
efficiency of intraband-mediated HHG Mikhailov and Ziegler (2008); Ishikawa
(2010); Al-Naib  _et al._ (2014). However, most measurements of HHG in
graphene have identified a thermal origin and focused on THz frequencies
Bowlan  _et al._ (2014); Hafez  _et al._ (2018); Kovalev  _et al._ (2021).
Recent work has been reported on the observation of mid-IR HHG up to the fifth
harmonic in multilayer graphene Baudisch  _et al._ (2018), as well as mid-IR
HHG from monolayer graphene, though this has required extremely high peak
intensities of $\sim 1\,$TW/cm2 Yoshikawa  _et al._ (2017), at which ablation
of the material is expected to take place Roberts  _et al._ (2011). An
empirical demonstration of the sought-after efficient mid-IR HHG from
monolayer graphene is still awaiting, to the best of our knowledge.

Figure 1: Graphene heterostructure sample and setup. a. Graphene is separated
by a dielectric material (h-BN or a $5\,$nm layer of $\mathrm{Al_{2}O_{3}}$)
of thickness s from an array of gold nanoribbons (width W, gap g) grown on
top. Pump light is normally incident on the sample and the generated harmonics
are collected in transmission. b. The sample shows no third-harmonic signal
when there is only substrate or gold nanoribbons. c. A significant third-
harmonic signal is observed with bare graphene and a strongly enhanced signal
when gold nanoribbons are added on top of the graphene. d. Scanning electron
microscopy (SEM) image of an array. e. Experimental setup of the modified
z-scan measurement. An optical parametric oscillator provides 3900 nm
linearly-polarized femptosecond laser pulses at a 76 MHz repetition rate. The
half-wave plate (HWP), together with the polarizer (pol), allows rotating the
polarization of the light that is then focused down with a short-focal lens
(L1). A second lens (L2) collimates both the initial light and the THG signal
excited in the sample that is moved in the z direction, along the focal point
of the beam. A spectral filter (SF) isolates the nonlinear signal that is then
sent to either a superconducting nanowire single-photon detector (SNSPD) or a
spectrometer to measure THG or FHG, respectively. Figure 2: Characterization
of the nonlinear harmonic signals. a. Spectrum of the THG signal measured with
a Fourier transform spectrometer ($15\,\mathrm{nm}$ resolution), along with
the spectrum of the FHG signal measured with a single-photon spectrometer
($1\,\mathrm{nm}$ resolution). b. Typical THG and FHG signals as a function of
input polarization relative to the ribbon array. As expected, the gold
nanoribbons impose a $\cos^{6}{\theta}$ and $\cos^{10}{\theta}$ dependence for
THG and FHG signals, respectively, maximizing the output when the incident
light is perpendicular to the direction of the nanoribbons. All of the data in
this figure are taken with HBN_70, although the trends are the same with other
parameters. Figure 3: Spacer dependence of the nonlinear enhancement. a. and
b. show the THG and FHG measured power as a function of the input power for
graphene with gold nanoribbons and bare graphene (red squares), respectively,
for a $5\,\mathrm{nm}$ $\mathrm{Al_{2}O_{3}}$ spacer (dark green/blue squares)
and a monolayer h-BN spacer (bright green/blue triangles). For the THG
measurements, we present the results from ALO_500 and HBN_200, while for the
FHG measurements, we show results from ALO_500 and HBN_70. The pump wavelength
is $3900\,$nm for all measurements. In the device with $\mathrm{Al_{2}O_{3}}$
spacer, the THG and FHG signals are enhanced by a factor $26\pm 1$ and $260\pm
30$, respectively. For the h-BN spacer, the enhancement factors are $1600\pm
500$ and $4100\pm 600$, respectively. Error estimates are all $68\,$%
confidence intervals from the fit. Figure 4: Nanoribbon width dependence of
nonlinear enhancement with h-BN spacer. a. and b. are the power dependence of
the THG and FHG signals for bare graphene (red) and graphene with gold
nanoribbons of different widths (yellow through blue). For all the data in
this figure, the pump wavelength is $3900\,$nm, resulting in third- and fifth-
harmonic signals at $1300\,$nm and $780\,$nm, respectively. The interception
points found with the fits, whose slopes are set to 3 and 5, respectively,
enable calculating the enhancement of the signals due to the field
confinement. In THG, the maximum enhancement occurs at a nanoribbon width
$W=200\,\mathrm{nm}$ and for FHG it occurs at $70\,\mathrm{nm}$. c. and d.
summarize $\chi^{(3)}$ and $\chi^{(5)}$ for different nanoribbon widths and
given powers $P=0.4\,\mathrm{mW}$ and $P=10.6\,\mathrm{mW}$, respectively.
Inset: Visible reflection spectra of the $100\,$nm, $200\,$nm, and $500\,$nm
structures, colored light to dark blue as in the color bar.

In this paper, we report on the observation of fifth-harmonic generation (FHG)
from monolayer graphene in the mid-IR spectral range with modest pumping peak
intensities of $<1\,$GW/cm2. This result is made possible by the fabrication
of heterostructures in which graphene is accompanied by metallic components to
boost the external near field actually acting on the carbon sheet. More
precisely, we fabricate heterostructures made of graphene and gold nanoribbons
that are separated by an insulator layer of either a $5\,$nm-thick layer of
aluminium oxide ($\mathrm{Al_{2}O_{3}}$) or a monolayer of hexagonal boron
nitride (h-BN). Similar heterostructures have also been used to launch
acoustic graphene plasmons Iranzo  _et al._ (2018); Lee  _et al._ (2019);
Epstein  _et al._ (2020), observe plasmon-mediated third-harmonic generation
(THG) Alonso Calafell  _et al._ (2021), and enhance HHG at THz frequencies
Deinert  _et al._ (2020). Here, the nanoribbons act as mid-IR antennas,
serving to enhance the incident pump field in the graphene layer and,
consequently, boost the resulting harmonic generation. Importantly, we
corroborate that ribbons exhibit no observable THG or FHG signals without the
graphene layer. Upon investigation of the role of the insulator spacer
thickness $s$ and the ribbon width $W$, we find that different ribbon widths
are obtained to optimize the THG and FHG efficiencies, a surprising result
that we attribute to an increased absorption of the fifth harmonic field by
plasmonic modes in the metal ribbons at the generated optical wavelengths. Our
findings are pivotal for the design of high-efficiency nanoscale nonlinear
frequency-conversion devices.

A schematic of our heterostructure design is shown in Fig. 1a. The samples
consist of a $\mathrm{SiO_{2}}$-$\mathrm{Si}$-$\mathrm{SiO_{2}}$ substrate, on
which a monolayer of graphene is grown by chemical vapor deposition. This is
then covered by a dielectric: either a monolayer of h-BN or a $5\,$nm film of
$\mathrm{Al_{2}O_{3}}$. On top of that, arrays of gold nanoribbons with
different widths $W$ and fixed inter-ribbon gaps of $g=50\,$nm are etched
using electron beam lithography and lift off (see Refs. Iranzo  _et al._
(2018); Alonso Calafell  _et al._ (2021) for more details). Fig. 1b shows a
typical scanning electron microscopy (SEM) image of the resulting high-quality
metal edges of a representative sample. We will refer to the different
structures as follows: the samples with the $5\,$nm $\mathrm{Al_{2}O_{3}}$
spacer and ribbon width $W$ are called ALO_W, while we refer to the samples
with the monolayer h-BN spacer and ribbon width $W$ as HBN_W. For example,
ALO_200 is the heterostructure with a $5\,$nm $\mathrm{Al_{2}O_{3}}$ spacer,
$g=50\,$nm and $W=200\,$nm.

To measure the nonlinear emission produced by our samples we use a modified
$z$-scan setup, in which the nonlinear signal (either FHG or THG) is measured
while the sample is moved along the z axis through the focus of the laser beam
(see Fig. 1e). Our pump beam is a linearly-polarized pulsed laser with a $\sim
260\,$fs pulse duration, a central wavelength of $3900\,\mathrm{nm}$
($0.225\,$eV), and a $76\,$MHz repetition rate. This beam is generated by an
optical parametric oscillator (OPO), fed by a mode-locked femtosecond
Ti:saphhire laser. We use a half-wave plate (HWP) to rotate the polarization
of the incoming beam to that set by the polarizer (pol). A lens with a
$5.95\,\mathrm{mm}$ focal length focuses the pump beam down to a waist of
$\sim 13.5\,\mathrm{\mu m}$. When the sample is moved parallel to the pump
beam (along the $z$ axis), the nonlinear emission occurs most efficiently
where the fluence is maximum (i.e. at the focal point). Afterwards, a lens
with an $11\,\mathrm{mm}$ focal length collimates the beam, which is then sent
through a set of spectral filters that separate the nonlinear emission from
the pump beam. The resulting nonlinear signal is coupled into a multimode
fiber, which can be sent to one of two detectors. The THG signal (at $\sim
1300\,\mathrm{nm}$) is measured using a large-area SNSP detector from
PhotonSpot with a $\sim 60\,\%$ detection efficiency at this wavelength. We
combine this detector with a Gemini interometer from NIREOS to measure the
spectrum of the THG. To acquire the FHG signal, which is centered around
$780\,\mathrm{nm}$, we connect the fiber to a single-photon sensitive silicon
Andor spectrometer with a resolution of $1\,\mathrm{nm}$ and a detection
efficiency of $\sim 10\,\%$.

To verify the origin of our nonlinear signals, we perform $z$-scans on various
regions of our samples. As shown in Fig. 1b, neither the substrate (green
diamonds) nor the gold nanoribbons without the graphene layer (green circles)
display any measureable nonlinear response. However, both graphene (red
triangles) and graphene with nanoribbons (red squares) show a significant
nonlinear signal at the focal point of the pump beam (see Fig. 1c). These
control measurements demonstrate that only the graphene layer contributes to
any measurable nonlinearity in our setup. Moreover, they already show the very
large nonlinear enhancement provided by the gold nanoribbons.

To further verify that our observed nonlinear signals are indeed associated
with HHG, we measure their spectra. As displayed in Fig. 2a, we find a THG
signal at $\sim 1300\,\mathrm{nm}$ and a FHG signal at $\sim
780\,\mathrm{nm}$. Both of them are at the expected wavelengths for our
$3900\,\mathrm{nm}$ pump beam. In addition to these two signals, we observe
broadband white light generation due to thermal photoluminescence from the
graphene layer. This has been previously reported in Ref. Lui  _et al._
(2010), and, as we discuss in the Appendix, the model presented therein fits
our observed spectrum.

Owing to the geometry of the gold nanoribbons, the pump field enhancement is
strongly dependent on the polarization of the input light. As shown in Fig.
2b, the enhancement and the resulting nonlinear signal are maximized when the
light polarization is perpendicular to the direction of the nanoribbons.
However, when the polarization is rotated away from such perpendicular
direction, the nonlinear signal decreases as $(\cos{\theta})^{6}$ and
$(\cos{\theta})^{10}$, for the THG and FHG signals, respectively, reaching a
minimum when the polarization is oriented along the nanoribbons. At this
point, the signal strength is even lower than that of planar graphene without
gold nanoribbons because of the screening produced by currents induced in the
metal. In the Appendix, Fig. 5 shows similar measurements performed on
graphene without the nanoribbons, which demonstrate only a small (few percent)
polarisation dependence of the nonlinear signals, attributed to birefringence
of the silicon substrate.

To better understand the nature of the enhancement, we first study its
dependence on the material and thickness $s$ of the dielectric spacer between
the graphene and gold nanoribbons, and then proceed to characterize the
dependence of the third- and fifth-harmonic signals on the nanoribbon width.
To quantify the effect of the spacer material and thickness, in Fig. 3a we
show the enhancement of the THG signals for both ALO_500 (dark green squares)
and HBN_200 (bright green triangles). These are compared to the signal from
bare graphene (red squares). We fit the nonlinear signals as a function of
input power (straight lines) assuming a third order dependence and with the
interception point in the logarithmic plot (which determines the magnitude of
the third-order signal) as the only fitting parameter. The enhancement factor
(EF) is then given by the difference between these interception points. We
find an EF of $26\pm 1$ with the $5\,\mathrm{nm}$ $\mathrm{Al_{2}O_{3}}$
spacer and an EF of $1600\pm 500$ with the monlayer h-BN spacer. In Fig. 3b,
we show a similar analysis for the FHG signals, assigning a fifth-order
dependence of the signal. Here, we present measurements from ALO_500 (dark
blue squares) and HBN_70 (bright blue squares). We find an EF of $260\pm 30$
for the $5\,\mathrm{nm}$ $\mathrm{Al_{2}O_{3}}$ spacer and an EF of $4100\pm
600$ with the h-BN spacer. We discuss the width dependence in detail below.

Our power-scaling measurements indicate that we remain in the perturbative
regime for both the THG and FHG measurements, in which the power of the
$n^{\mathrm{th}}$-order harmonic scales with $p_{\mathrm{pump}}^{n}$, where
$p_{\mathrm{pump}}$ is the average pump power Wildenauer (1987); Li  _et al._
(1989). Deviation from the expected power dependencies is only observed for
the highest incident powers used here. As pointed out in Ref. Alonso Calafell
_et al._ (2021), this is caused by an increased electron temperature Yu  _et
al._ (2017), which reduces the nonlinearity. Efficient HHG for higher-
harmonics in gas media typically takes place in the non-perturbative regime,
wherein the power of all harmonics scales as $p_{\mathrm{pump}}^{2}$ Higuchi
_et al._ (2014). This limit has also been reached in solid-state HHG Schubert
_et al._ (2014); Liu  _et al._ (2017), including in graphene Yoshikawa  _et
al._ (2017). We were unable to reach this limit either in extended graphene or
in our heterostructures without damaging the samples. In particular, we found
that for, average pump powers above $\sim 30\,$mW ($\sim 1\,$GW/cm2), the
power of the THG and FHG signals began to permanently decrease after several
seconds, presumably due to laser-induced structural modifications, as found in
Ref. Roberts  _et al._ (2011).

To quantify the effect of the nanoribbon width on the THG and FHG signal
enhancement, in Fig. 4 we compile a series of measurements performed for
devices with the monolayer h-BN spacer. For the $\mathrm{Al_{2}O_{3}}$
samples, we were only able to measure a FHG signal using ALO_500. In our
attempts to measure FHG in the $\mathrm{Al_{2}O_{3}}$ samples with other
widths, we observed laser-induced damage of the nanoribbons before detecting
FHG. As such, we will not discuss the width dependence in that sample. In
panels a and b, we show the third- and fifth-harmonic signals as a function of
input power for nanoribbon widths of $20\,$nm to $500\,$nm with the monolayer
h-BN spacer (yellow to blue) and bare graphene (red squares). Again, the
points represent experimental data and the straight lines are fits, where only
the interception points are free parameters (one per line) and the slopes are
set to 3 and 5, respectively. Due to the long acquisition times and slow laser
power fluctuations, we were unable to acquire reliable FHG power-scaling plots
for HBN_20, HBN_50 and HBN_500. Note also that it was necessary to carry out
these measurements in planar graphene using higher powers, which resulted in
visible damage to these samples for the highest powers used.

In panels c and d, we calculate the effective third-order and fifth-order
susceptibilities of the graphene heterostructures Jiang  _et al._ (2018); Boyd
(2003), and find the maxima to be $\chi^{(3)}=5.6\cdot 10^{-6}$ and
$\chi^{(5)}=1.8\cdot 10^{-14}\,\mathrm{esu}$. These values are several orders
of magnitude larger than for bare graphene, highlighting the benefit of such
nanostructures in enhancing the nonlinear properties of 2D materials. Looking
directly at the conversion efficiencies of the different processes (defined as
the ratio of the generated TH or FH power over the input power), we find a
maximum THG conversion efficiency of $\sim 2\times 10^{-8}$ with HBN_200 at a
pump power of $2.5$ mW, and an FHG conversion efficiency of $\sim 4\times
10^{-12}$ in HBN_70 with a pump power of $10$ mW. While these conversion
efficiencies are lower than those that have been achieved in the THz regime
Bowlan  _et al._ (2014); Hafez  _et al._ (2018); Kovalev  _et al._ (2021),
they are significantly higher than in previous work within the mid-IR due to
the enhancement produced by the heterostructures in our work. For example,
Ref. Baudisch  _et al._ (2018) found THG and FHG conversion efficiencies of
$3\times 10^{-10}$ and $9\times 10^{-14}$, respectively, in five-layer
graphene.

Interestingly, we observe that the optimal widths for harmonic generation is
different for THG and FHG: the THG signal is strongest for ribbons with a
width of $200\,\mathrm{nm}$, an effect which arises due to a maximum in field
enhancement for larger ribbons Alonso Calafell  _et al._ (2021). Somewhat
unexpectedly, the FHG signal is strongest for the $70\,\mathrm{nm}$ ribbons.
As the properties of pump beam are nominally identical (other than the power),
this discrepancy most likely originates in the width-dependent additional
absorption of the generated beams, as the metal nanoribbons in our
heterostructures display plasmonic modes in the visible range. The width
dependence of this plasmonic resonance is illustrated in the inset of Fig. 4d,
which shows the measured reflection spectra for samples with ribbons of widths
$100\,$nm, $200\,$nm, and $500\,$nm. We observe that the peak of this
plasmonic mode red shifts from about $675\,$nm for the $100\,$nm structure to
$750\,$nm for the $500\,$nm structure. The encroaching plasmon resonance
attenuates the FHG signal (generated at $780\,$nm) more for the wider ribbons,
while the THG signal does not undergo this additional absorption. This
plasmon-quenching effect accounts for the weaker-than-expected FHG for wider
ribbons.

In conclusion, we have demonstrated that graphene-based heterostructures can
be used to enhance nonlinear conversion efficiencies in the mid-IR spectral
range at relatively low pump powers. This has enabled us to observe FHG in the
perturbative regime, wherein the generated nonlinear signal scales with the
input pump power to the power of five. Such enhancement also allows us to
observe third- and fifth-harmonic generation in graphene at modest pump powers
in the mid-IR. We study the dependence of the nonlinear signals on the ribbon
width and spacer thickness, as well as pump power and polarization, and
demonstrate enhancement factors for third- and fifth-harmonic generation of up
to 1600 and 4100, respectively, compared to bare graphene.

We have shown that the ribbon width is one of the key parameters to adjust in
order to optimise nonlinear conversion. For THG, we observe the largest
nonlinear signals for ribbon widths of $\sim 200$ to $500\,$nm, an effect
stemming from the optimal enhancement of the pumping field acting on the
graphene. For FHG, we observe a maximum nonlinear conversion for smaller
ribbon widths of $70\,$nm. For larger ribbon widths, the FHG is likely
quenched due to spectral overlap with the plasmon resonance of the ribbons,
which shifts towards the FHG emission wavelength of $780\,$nm for larger
ribbons. This demonstrates the fine tuning required in order to maximise a
structure for a specific generation wavelength. This information is crucial in
the design of future plasmonic, nanoscale HHG devices, particularly those
intended to generate visible and near-IR wavelengths.

## acknowledgments

We acknowledge support from the Austrian Science Fund (FWF) through BeyondC
(F7113) and Research Group 5 (FG5), the Air Force Office of Scientific
Research (AFOSR) via PhoQuGraph (FA8655-20-1-7030), the Engineering and
Physical Sciences Research Council (EPSRC) of the United Kingdom via the EPSRC
Centre for Doctoral Training in Metamaterials (Grant No. EP/L015331/1), the
Spanish MICINN (PID2020-112625GB-I00 and SEV2015-0522), the Catalan CERCA, and
Fundació Privada Cellex, the European Union’s Horizon 2020 research and
innovation program under the Marie Sklodowska‐Curie grant agreement No 801110
and the Austrian Federal Ministry of Education, Science and Research (BMBWF).

## Appendix

###  .1 Extraction of the third-order and fifth-order susceptibilities—

Experimentally, $\chi^{(3)}$ and $\chi^{(5)}$ have been estimated following
the procedure reported in Ref Jiang  _et al._ (2018). The measured average
power is proportional to the squared modulus of the electric field through the
following relation:

$P(\omega)=\frac{1}{8}\left(\frac{\pi}{\ln{2}}\right)^{3/2}f\tau
W^{2}n\epsilon_{0}c\frac{|E(\omega)|^{2}}{2},$ (1)

where $f$ is the laser repetition rate, $\tau$ is the temporal pulse width,
$W$ is the Gaussian beam diameter, $n$ is the refractive index, and
$\epsilon_{0}$ and $c$ are the permittivity and speed of light in vacuum. In
the experiment here reported, $f=76\,\mathrm{MHz}$, $\tau=260\,\mathrm{fs}$
(FWHM), $W=27\,\mathrm{\mu m}$, and the refractive index $n=2.4$ is considered
to be constant.

The generated THG and FHG electric fields in a graphene monolayer of thickness
$d_{gr}=0.33\,\mathrm{nm}$ are related to the electric field at the
fundamental pump frequency $E(\omega_{0})$ through the following relations
Boyd (2003); Shen (1984):

$\displaystyle
E(3\omega_{0})=i\frac{\omega_{0}}{8nc}\chi^{(3)}d_{gr}E(\omega_{0})^{3},$ (2a)
$\displaystyle
E(5\omega_{0})=i\frac{\omega_{0}}{8nc}\chi^{(5)}d_{gr}E(\omega_{0})^{5}.$ (2b)

From Eqs. (2), it is thus possible to calculate the effective $\chi^{(3)}$ and
$\chi^{(5)}$ susceptibilities. Note that all the quantities in Eqs. (1) and
(2) are given in SI units, while the $\chi^{(3)}$ and $\chi^{(5)}$ values in
Fig. 3(c-d) are given in the electrostatic system of units (esu). The
expressions

$\displaystyle\chi^{(3)}(\mathrm{SI})=\frac{4\pi}{(3*10^{4})^{2}}\,\,\chi^{(3)}(\mathrm{esu}),$
(3a)
$\displaystyle\chi^{(5)}(\mathrm{SI})=\frac{4\pi}{(3*10^{4})^{4}}\,\,\chi^{(5)}(\mathrm{esu}),$
(3b)

provide the relative conversion factors.

Figure 5: Polarization dependence of the transmission through the substrate.
Red points show the normalized polarization-dependent transmission through the
substrate. Green points are the normalized THG power as a function of input
polarization. Blue points are the cube of the transmitted pump power in the
red points. All measurements are performed with a carrier wavelength of
$3900\,\mathrm{nm}$. Figure 6: Ultrafast thermal emission from planar
graphene. Spectral emission of monolayer graphene on a Si wafer, pumped by
$260\,$fs pulses with a peak intensity of $100\,$MW/cm2. The experimental data
(orange) shows not only the FHG signal around $780\,$nm, but also a broadband
background. This latter is modelled as ultrafast thermal emission by hot
carriers in graphene (blue). The inset shows the corresponding time-dependent
electron temperature.

###  .2 Polarization dependence

In this section, we discuss the polarization dependence of THG from bare
graphene. The upper panel of Fig. 5 shows the transmission of the pump beam
through the sample as a function of the orientation of its linear
polarization. In this experiment, far-field light with a wavelength of
$3900\,\mathrm{nm}$ is normally incident on the sample from the side of the
substrate, composed of a stack of two layers of silica and one layer of
silicon (see Fig. 1). The silica layers have a thickness of
$285\,\mathrm{nm}$, while the silicon layer has a thickness of
$500\,\mathrm{\mu m}$ (see Ref. Iranzo  _et al._ (2018) for more detail of the
sample fabrication). The measurements reported in Fig. 5 are carried out on an
area of the sample containing only graphene (i.e., there are no gold
nanoribbons), just like in the violet triangles in Fig. 1a of the main text.
The red points in Fig.5 show a periodic pump transmission through the sample
as a function of the incoming linear polarization angle. The measured relative
difference in transmission is about 5%. In the lower panel, the green points
show the polarization-dependent normalized THG power together with the
normalized cubed transmitted power (blue points). The quantitative agreement
between the two explains the polarization dependence of the THG signal.

The physical origin of the polarization effect could be explained by a slight
birefringence in the silicon Ref. Krüger  _et al._ (2015). Another possible
explanation is the dependence of the Fresnel reflection coefficient on
polarization for off-normal incidence Saleh and Teich (2019). Indeed, for an
interface between air and silica, an incidence angle of $\sim
7\,\mathrm{\degree}$ relative to the interface normal is enough to induce a
relative difference in reflectance by $5\,\mathrm{\%}$ between $s$ and $p$
polarizations.

###  .3 Thermal radiation background

In this section, the physical origin of the background signal observed in the
FHG spectrum of Fig. 2 is identified as ultrafast photoluminescence from
graphene. The significant broadband light emission from graphene under
excitation by femtosecond laser pulses has been studied in Ref. Lui  _et al._
(2010) in great detail and is employed here to model our measurements. The
ultrafast heating of the electron gas is modelled with a two-temperature
model, comprising the electronic temperature $T_{el}$ and the temperature of
strongly coupled optical phonons $T_{op}$:

$\displaystyle\frac{dT_{el}(t)}{dt}$ $\displaystyle=$
$\displaystyle\frac{I(t)-\Gamma\left(T_{el},T_{op}\right)}{c_{e}\left(T_{el}\right)},$
(4) $\displaystyle\frac{dT_{op}(t)}{dt}$ $\displaystyle=$
$\displaystyle\frac{\Gamma\left(T_{el},T_{op}\right)}{c_{op}\left(T_{op}\right)}-\frac{T_{op}(t)-T_{0}}{\tau_{op}},$
(5)

with a phonon decay time $\tau_{op}=1\,$ps and a specific heat capacity of the
phonons $c_{op}\left(T_{op}\right)$ and electrons $c_{el}\left(T_{el}\right)$.

We define the electron-phonon energy exchange rate $\Gamma$ as

$\begin{split}\Gamma\left(T_{el},T_{op}\right)&=\alpha\cdot(n_{op}(T_{op})+1)\int
E\cdot(E-E_{op})\cdot
f_{el}(E,T_{el})\cdot[1-f_{el}((E-E_{op}),T_{el})]\text{d}E\\\
&-\alpha\cdot(n_{op}(T_{op}))\int E\cdot(E-E_{op})\cdot
f_{el}(E-E_{op},T_{el})\cdot[1-f_{el}((E),T_{el})]dE,\end{split}$ (6)

where $E_{op}=0.18$ eV is the phonon energy,
$n\left(T_{op}\right)=\left[\exp\left(E_{op}/kT_{op}\right)-1\right]^{-1}$
represents the phonon population at temperature $T_{op}$, and
$f\left(E,T_{el}\right)=\left[\exp\left(E/kT_{el}\right)+1\right]^{-1}$ is the
Fermi-Dirac distribution for electrons at temperature $T_{el}$. The phonon
coupling constant $\alpha$ corresponds to

$\alpha=3\dfrac{9\left(\frac{\partial t}{\partial b}\right)^{2}}{\pi\rho
E_{op}\hbar^{3}v_{F}^{4}},$ (7)

with $\frac{\partial t}{\partial b}=45\,$eV/nm and $\rho=7.6\times
10^{-7}\leavevmode\nobreak\ \text{kg}/\text{m}^{2}$ as the density of graphene
[Ref2].

The calculated electron temperature in graphene for $260\,$fs pulses and an
incident peak intensity of $100\,$MW/cm2 is plotted in the inset of Fig. S6.
To get the corresponding spectral radiant fluency, we use Plancks’s law:

$\mathcal{F}\left(\lambda,T_{el}\right)=\varepsilon(\lambda)\frac{2hc^{2}}{\lambda^{5}}\int\left[\exp\left(\frac{hc/\lambda}{kT_{el}(t)}\right)-1\right]^{-1}dt.$
(8)

The emissivity $\varepsilon(\lambda)$ of graphene is calculated via the
absorption of light coming from air, taking into account the wavelength-
dependent modulation introduced by the cavity formed by the $285\,$nm layer of
glass between the graphene and Si. Assuming that the illuminated graphene area
acts like a thermal Lambert emitter, we estimate our lens to collect up to an
angle $\theta=\arctan(\frac{D_{eff}/2}{f})$, where $f=11\,$mm is the focus
length and the effective diameter is scaled to match the experimental signal
via $D_{eff}=30\%\times 5.5\,$mm. This is the only scaling parameter, which
has been tweaked to match the thermal estimate with the background that we
observe in the experiment, as shown in Fig. 6. It compensates alignment
issues, as the apparatus is designed to measure HHG, but not white light
generation. It also compensates for filter-reflection losses, coupling into
the spectrometer, and similar effects. It does not compensate for wavelength-
dependent performance, which we assume to be most influenced by the rough
estimate of our absorbed mid-IR fluence ($\sim 1.5\,$%), but also by the
collection lens having a chromatic alignment and wavelength-dependent
efficiency. Taking all of these aspects into account, our estimate matches
sufficiently well to identify ultrafast thermal emission of hot electrons as
the dominant source of the background observed in the FHG measurements.

## References

  * McPherson  _et al._ (1987) A. McPherson, G. Gibson, H. Jara, U. Johann, T. S. Luk, I. A. McIntyre, K. Boyer, and C. K. Rhodes, J. Opt. Soc. Am. B 4, 595 (1987).
  * Ferray  _et al._ (1988) M. Ferray, A. l’Huillier, X. F. Li, L. A. Lompre, G. Mainfray, and C. Manus, Journal of Physics B: Atomic, Molecular and Optical Physics 21, L31 (1988).
  * Ghimire  _et al._ (2011) S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro, and D. A. Reis, Nature physics 7, 138 (2011).
  * Schubert  _et al._ (2014) O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, _et al._ , Nature photonics 8, 119 (2014).
  * Luu  _et al._ (2015) T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T. Hassan, and E. Goulielmakis, Nature 521, 498 (2015).
  * Langer  _et al._ (2016) F. Langer, M. Hohenleutner, C. P. Schmid, C. Pöllmann, P. Nagler, T. Korn, C. Schüller, M. Sherwin, U. Huttner, J. Steiner, _et al._ , Nature 533, 225 (2016).
  * You  _et al._ (2017) Y. S. You, D. A. Reis, and S. Ghimire, Nature physics 13, 345 (2017).
  * Jiang  _et al._ (2018) T. Jiang, D. Huang, J. Cheng, X. Fan, Z. Zhang, Y. Shan, Y. Yi, Y. Dai, L. Shi, K. Liu, _et al._ , Nature Photonics 112, 430 (2018).
  * Soavi  _et al._ (2018) G. Soavi, G. Wang, H. Rostami, D. G. Purdie, D. De Fazio, T. Ma, B. Luo, J. Wang, A. K. Ott, D. Yoon, S. A. Bourelle, J. E. Muench, I. Goykhman, S. Dal Conte, M. Celebrano, A. Tomadin, M. Polini, G. Cerullo, and A. C. Ferrari, Nature nanotechnology 13, 583 (2018).
  * Mikhailov (2016) S. A. Mikhailov, Phys. Rev. B 93, 085403 (2016).
  * Rostami  _et al._ (2017) H. Rostami, M. I. Katsnelson, and M. Polini, Phys. Rev. B 95, 035416 (2017).
  * Mikhailov and Ziegler (2008) S. A. Mikhailov and K. Ziegler, Journal of Physics: Condensed Matter 20, 384204 (2008).
  * Marini  _et al._ (2017) A. Marini, J. D. Cox, and F. J. García de Abajo, Phys. Rev. B 95, 125408 (2017).
  * Kumar  _et al._ (2013) N. Kumar, J. Kumar, C. Gerstenkorn, R. Wang, H.-Y. Chiu, A. L. Smirl, and H. Zhao, Physical Review B 87, 121406 (2013).
  * Hong  _et al._ (2013) S.-Y. Hong, J. I. Dadap, N. Petrone, P.-C. Yeh, J. Hone, and R. M. Osgood Jr, Physical Review X 3, 021014 (2013).
  * Hendry  _et al._ (2010) E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, and S. A. Mikhailov, Phys. Rev. Lett. 105, 097401 (2010).
  * Baudisch  _et al._ (2018) M. Baudisch, A. Marini, J. D. Cox, T. Zhu, F. Silva, S. Teichmann, M. Massicotte, F. H. L. Koppens, L. S. Levitov, F. J. García de Abajo, and J. Biegert, Nat. Commun. 9, 1018 (2018).
  * Beckerleg  _et al._ (2018) C. Beckerleg, T. J. Constant, I. Zeimpekis, S. M. Hornett, C. Craig, D. W. Hewak, and E. Hendry, Applied Physics Letters 112, 011102 (2018).
  * Constant  _et al._ (2016) T. J. Constant, S. M. Hornett, D. E. Chang, and E. Hendry, Nat. Phys. 12, 124 (2016).
  * Bowlan  _et al._ (2014) P. Bowlan, E. Martinez-Moreno, K. Reimann, T. Elsaesser, and M. Woerner, Physical Review B 89, 041408 (2014).
  * Hafez  _et al._ (2018) H. A. Hafez, S. Kovalev, J.-C. Deinert, Z. Mics, B. Green, N. Awari, M. Chen, S. Germanskiy, U. Lehnert, J. Teichert, _et al._ , Nature 561, 507 (2018).
  * Kovalev  _et al._ (2021) S. Kovalev, H. A. Hafez, K.-J. Tielrooij, J.-C. Deinert, I. Ilyakov, N. Awari, D. Alcaraz, K. Soundarapandian, D. Saleta, S. Germanskiy, _et al._ , Science advances 7, eabf9809 (2021).
  * Yoshikawa  _et al._ (2017) N. Yoshikawa, T. Tamaya, and K. Tanaka, Science 356, 736 (2017).
  * Koppens  _et al._ (2011) F. H. Koppens, D. E. Chang, and F. J. García de Abajo, Nano Lett. 11, 3370 (2011).
  * Gullans  _et al._ (2013) M. Gullans, D. E. Chang, F. H. L. Koppens, F. J. García de Abajo, and M. D. Lukin, Phys. Rev. Lett. 111, 247401 (2013).
  * Jablan and Chang (2015) M. Jablan and D. E. Chang, Phys. Rev. Lett. 114, 236801 (2015).
  * Alonso Calafell  _et al._ (2019) I. Alonso Calafell, J. Cox, M. Radonjić, J. Saavedra, F. García de Abajo, L. Rozema, and P. Walther, npj Quantum Information 5, 1 (2019).
  * Ishikawa (2010) K. L. Ishikawa, Physical Review B 82, 201402 (2010).
  * Al-Naib  _et al._ (2014) I. Al-Naib, J. E. Sipe, and M. M. Dignam, Phys. Rev. B 90, 245423 (2014).
  * Roberts  _et al._ (2011) A. Roberts, D. Cormode, C. Reynolds, T. Newhouse-Illige, B. J. LeRoy, and A. S. Sandhu, Applied Physics Letters 99, 051912 (2011).
  * Iranzo  _et al._ (2018) D. A. Iranzo, S. Nanot, E. J. C. Dias, I. Epstein, C. Peng, D. K. Efetov, M. B. Lundeberg, R. Parret, J. Osmond, J.-Y. Hong, J. Kong, D. R. Englund, N. M. R. Peres, and F. H. L. Koppens, Science 360, 291 (2018).
  * Lee  _et al._ (2019) I.-H. Lee, D. Yoo, P. Avouris, T. Low, and S.-H. Oh, Nature nanotechnology 14, 313 (2019).
  * Epstein  _et al._ (2020) I. Epstein, D. Alcaraz, Z. Huang, V.-V. Pusapati, J.-P. Hugonin, A. Kumar, X. M. Deputy, T. Khodkov, T. G. Rappoport, J.-Y. Hong, _et al._ , Science 368, 1219 (2020).
  * Alonso Calafell  _et al._ (2021) I. Alonso Calafell, L. A. Rozema, D. Alcaraz Iranzo, A. Trenti, P. K. Jenke, J. D. Cox, A. Kumar, H. Bieliaiev, S. Nanot, C. Peng, _et al._ , Nature Nanotechnology 16, 318 (2021).
  * Deinert  _et al._ (2020) J.-C. Deinert, D. Alcaraz Iranzo, R. Perez, X. Jia, H. A. Hafez, I. Ilyakov, N. Awari, M. Chen, M. Bawatna, A. N. Ponomaryov, _et al._ , ACS nano 15, 1145 (2020).
  * Lui  _et al._ (2010) C. H. Lui, K. F. Mak, J. Shan, T. F. Heinz, _et al._ , Physical review letters 105, 127404 (2010).
  * Wildenauer (1987) J. Wildenauer, Journal of applied physics 62, 41 (1987).
  * Li  _et al._ (1989) X. Li, A. l’Huillier, M. Ferray, L. Lompré, and G. Mainfray, Physical Review A 39, 5751 (1989).
  * Yu  _et al._ (2017) R. Yu, A. Manjavacas, and F. J. Garcia de Abajo, Nature communications 8, 1 (2017).
  * Higuchi  _et al._ (2014) T. Higuchi, M. I. Stockman, and P. Hommelhoff, Physical review letters 113, 213901 (2014).
  * Liu  _et al._ (2017) H. Liu, Y. Li, Y. S. You, S. Ghimire, T. F. Heinz, and D. A. Reis, Nature Physics 13, 262 (2017).
  * Boyd (2003) R. W. Boyd,  _Nonlinear optics_ (Elsevier, 2003).
  * Shen (1984) Y.-R. Shen, New York, Wiley-Interscience, 1984, 575 p. (1984).
  * Krüger  _et al._ (2015) C. Krüger, D. Heinert, A. Khalaidovski, J. Steinlechner, R. Nawrodt, R. Schnabel, and H. Lück, Classical and Quantum Gravity 33, 015012 (2015).
  * Saleh and Teich (2019) B. E. Saleh and M. C. Teich,  _Fundamentals of photonics_ (John Wiley & Sons, 2019).


# Neural information coding for efficient spike-based image denoising

###### Abstract

In recent years, Deep Convolutional Neural Networks (DCNNs) have outreached
the performance of classical algorithms for image restoration tasks. However
most of these methods are not suited for computational efficiency and are
therefore too expensive to be executed on embedded and mobile devices. In this
work we investigate Spiking Neural Networks (SNNs) for Gaussian denoising,
with the goal of approaching the performance of conventional DCNN while
reducing the computational load. We propose a formal analysis of the
information conversion processing carried out by the Leaky Integrate and Fire
(LIF) neurons and we compare its performance with the classical rate-coding
mechanism. The neural coding schemes are then evaluated through experiments in
terms of denoising performance and computation efficiency for a state-of-the-
art deep convolutional neural network. Our results show that SNNs with LIF
neurons can provide competitive denoising performance but at a reduced
computational cost.

Index Terms—  Image denoising, Spiking Neural Networks, Neural information
coding, Neuromorphic computing.

##  1 Introduction and Related Work

Smartphone cameras, because of their reduced size and high pixel count, are
intrinsically more susceptible to noise than conventional digital cameras.
Image denoising algorithms are therefore intensively used in smartphones to
recover image quality by reducing the amount of noise of the raw image. Image
denoising performance have increased during the last few years and recent
methods based on Deep Convolutional Neural Networks (DCNNs) have provided very
high scores [1] to the point of outreaching classical spatial and patch-based
algorithms [2]. However, deploying AI-based algorithms on embedded devices
poses many problems. The limited amount of memory available, power consumption
and thermal dissipation are indeed critical for embedded battery powered
platforms. The field of neuromorphic engineering, especially SNNs, is emerging
as a new paradigm for the design of low-power and real-time information
processing hardware [3]. The spike information coding used by SNNs enables
sparse and event-based computation through the network. The combination of
these properties may lead to more energy efficient hardware implementations of
neural networks, allowing state-of-the-art AI algorithms to be executed on
mobile platforms with a reduced power budget [4]. However, to achieve these
energy gains while simultaneously reaching the level of performance of DCNNs,
SNNs must be able to encode analog data with high precision using very compact
codes, i.e. spike trains. In recent years, several training and conversion
methods have been proposed to improve the accuracy of SNNs on large-scale
machine learning tasks. To take advantage of better performance provided by
supervised learning, several methods have been developed to convert ANNs,
trained using standard schemes like backpropagation, into SNNs for event-
driven inference [5] [6]. The ANN-SNN conversion is based on the idea that
firing rates of spiking neurons should match the activations of analog
neurons. Rate-based conversion methods have achieved significant results over
the last few years, thus reducing the accuracy gap with ANNs [7]. However,
rate-based conversion methods have a major drawback since they require a large
amount of timesteps to precisely match the activations of analog neurons. A
new approach, called surrogate gradient learning, has been proposed to train
SNNs directly in the spike domain using standard supervised learning
algorithms [8]. Recent studies reported competitive performance on a series of
static and dynamic datasets using surrogate gradient training [9]. In this
paper, we will extend these previous works and study the trade off between
accuracy and efficiency of SNNs for the specific and uncovered case of image
denoising. This task is challenging for two reasons. First, as denoising is a
regression task, the network has to predict a continuous value (i.e. the noise
amplitude) for each pixel of the image. Moreover, state of the art results
have been obtained with very deep networks (17 layers or more). In Section 2
we study two spike coding approaches and we formalize the trade-off between
performance and activation sparsity in SNN. In Section 3 we propose, for the
first time, an image denoising solution based on SNN. The network trained
directly in the spike domain provides a level of performance close to the
state of the art CNN based solution. In the last section, we conclude the
paper and we discuss future work.

##  2 Neural coding

Neural coding schemes convert input pixels into spikes that are transmitted to
spiking neurons for information processing. Specifically, we are interested in
the two complementary neural coding procedures called _encoding_ and
_decoding_ that map analog values into train of spikes and its inverse. Two
types of encoding and decoding schemes are studied: the coding with
information conversion and the rate-coding. The following section presents an
analysis of the neural conversion with LIF neurons. The comparison between the
two coding scheme will be discussed in Sec. 2.2

###  2.1 Neural coding with LIF neurons

In this first coding scheme, that we call _Neural information conversion_ , a
LIF neuron, located in the first layer (encoding) of the network is fed with a
constant input $x$ through $T$ timesteps. We are interested in finding the
encoding function that defines the relation between $x$ and $z(t)$, the
spiking output. Let us first recall the equations that govern the dynamic of a
LIF neuron in the discrete case [10]:

$V[n]=V[n-1]+\frac{1}{\tau}(-(V[n-1]-V_{reset})+x[n])$ (1)

Whenever $V[n]\geq V_{th}$, where $V_{th}$ is the threshold voltage, the
neuron emits a binary spike. The spike train representation of the analog
input $x$, is thus encoded in the following function:

$z(t)=\sum_{j=1}^{T}\delta(t-t_{j})$ (2)

Where $\delta(t)$ is the Dirac delta function and $t_{j}$ are the spike-times
indexed by $j$.

In the following we also consider that the membrane potential is completely
discharged after a spike emission, $V_{reset}=0$. Let us find the value of $x$
that makes the neuron fires at each timestep, thus producing a constant firing
rate of 1. The conditions that produce this firing pattern are shown below:

$\begin{cases}V[n-1]=0\\\ V[n]\geq V_{th}\\\
x[n]=x,\text{$\forall\,n$}\end{cases}$ (3)

Since a spike has to be generated at the current timestep $n$, the membrane
potential is greater or equal to $V_{th}$. Moreover, with a firing rate of 1,
the neuron resets its membrane potential at each timestep, after the spike
emission. Therefore the membrane potential at the previous timestep, $V[n-1]$,
equals 0. Replacing conditions 3 into Eq. 1 we obtain:

$x\geq V_{th}.\tau$ (4)

When the input $x$ is greater than or equal to $V_{th}.\tau$, the LIF neuron
fires a spike at each timestep. Following the same reasoning, let us find the
values of $x$ that produce a firing rate of 0.5. In such a case, the neuron
periodically alternates between two states: [charge, fire&reset, charge,
fire&reset, …]. The conditions that lead to this behaviour are shown below:

$\begin{cases}V[n-1]=x/\tau\\\ V[n]\geq V_{th}\\\
x[n]=x,\text{$\forall\,n$}\end{cases}$ (5)

Replacing conditions 5 into Eq. 1 we obtain:

$V[n]=V_{th}=x/\tau+\frac{1}{\tau}(x-x/\tau)$ (6)

The values of $x$ that produces a firing rate of 0.5 are defined by:

$\frac{V_{th}}{2/\tau-1/\tau^{2}}\leq x\leq V_{th}\,.\tau$ (7)

The same approach can be used to determine the production conditions of the
other firing rates of a LIF model as depicted in figure 1.
Before proceeding in the analysis, it is now interesting to focus on the spike
patterns generated by a LIF neuron in more details. We may wonder, for
example, if all spiking patterns of a LIF neuron are allowed, as well as the
number of different firing rates. Let us start with a simple example where the
neuron codes information over $T=8$ timesteps. Note that a sequence that leads
to a generation of a spike (in the case of a constant input) must have the
following format: [charge] during $k$ timesteps, followed by [fire&reset].
Table 1, shows the spiking pattern generated by a LIF neuron, when simulated
for $T=8$ timesteps.

Timesteps | $f_{r}$ | $k$
---|---|---
0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |  |
1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1.0 | 0
0 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 0.5 | 1
0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0.25 | 2
0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0.25 | 3
0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0.125 | 4
0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0.125 | 5
0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0.125 | 6
0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0.125 | 7
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 8
Table 1:  Output codes of a LIF neuron stimulated with a constant value. $T$
represents the number of timesteps (here $T=8$). $f_{r}$ is the firing rate
and $k$ is the number of timesteps used for the charge phase before the first
spike. The value 1 means that a spike has been generated on the corresponding
timestep.

We can observe that the number of output codes does not depend on the input
$x$ but only on the value of $T$ (simulation timesteps). In fact, only $T+1$
output codes can be generated by a LIF neuron, when stimulated with a constant
input. We have so far characterized the _encoding_ function, the relation
between the analog input $x$ and the spike pattern, $z(t)$, generated at the
output of the neuron. The inverse process called _decoding_ aims at mapping
the inverse function, that is the reconstruction of the analog value
($\hat{x}$) from a spiking input. To do so, we use the firing-rate of the
neuron as the information carrying variable and express the decoded analog
output as follows, where $f_{r}$ denotes the firing rate:

$\hat{x}=f_{r}=\frac{1}{T}\sum^{T}z(t)$ (8)

The output codes obtained by simulating a LIF neuron are plotted in Fig. 1 as
a function of the input value $x$.

Fig. 1: Firing-rate as a function of the input $x$ for $T=8$
($V_{th}=1.0,\tau=2.0$). The thresholds for $f_{r}=1.0$ and $f_{r}=0.5$ and
$f_{r}=0.25$ are also shown.

As it can be observed, the conversion operated by a LIF neuron is highly non-
uniform as it provides more quantization steps at amplitudes near $V_{th}$
than at higher amplitudes. However, the quantization step sizes decrease while
approaching $V_{th}$. At the opposite, codes that carry a high $f_{r}$, have
large quantization step sizes. As an example and as shown in Fig. 1, a $f_{r}$
of 0.5 will be generated by the neuron when $x\in[1.33,2.0]$ .

###  2.2 Comparison between neural conversion and rate coding

To assess the performance of the neural conversion scheme, we use a set of 12
natural images [11] (Set12), to measure the accuracy of the quantizer. The
Peak-Signal-to-Noise ratio (PSNR) is used as quality criterion. The pixel
intensities of the test images are normalized in the interval $[1.0,2.0]$ to
match the conversion range of the LIF neuron shown in Fig. 1. The normalized
pixel intensities are fed (without noise) into a LIF neuron membrane for $T$
timesteps. The spikes generated by each neuron are then collected and an
estimate of each pixel value is computed using equation 8. We compare the
neural coding of LIF neurons with the rate coding, a well known scheme that
has been extensively used in the SNN community for coding dense information in
the spike domain [12] [13] [14]. In the rate coding scheme we assume that
spike trains are independent realizations of Poisson processes with rates
$r_{i}$, where the pixel intensity $x_{i}$ is the firing probability,
normalized between [0,1], at each timestep. Eq. 8 is also used to decode the
spiking output. The average PSNR on Set12 is shown on the left side of Fig. 2.

Fig. 2: PSNR (left) and $\theta$ (right) as a function of the number of
timesteps for the LIF neural conversion and rate coding schemes on the Set12
dataset. Each point of the curves represents an average over the number of
pixels of all the dataset images. Here $V_{th}=1$ and $\tau=2$ for the LIF
neurons.

As we can observe, the PSNR increases quickly for the LIF coding scheme and
saturates at $T\sim 10$. Adding more timesteps to the conversion does not
improve the image reconstruction. From the previous analysis, presented in
Sec. 2.1, we know that increasing $T$ adds more quantization intervals and
therefore also increases the number of output codes that can be generated by a
LIF neuron. However, as it can be observed from Fig. 1, the sizes of newly
added quantized intervals decrease fast and vanish when $x$ approaches
$V_{th}$. As a result, the non-uniform quantization scheme that emerges from
the LIF neuron does not allow decreasing the distortion between the original
and the quantized signals by adding more quantization bits, i.e. increasing
$T$. On the other hand, the rate conversion scheme does not set any limits on
the accuracy since the PSNR increases with a logarithmic shape as a function
of $T$. Let us now study a property of great interest for a neural information
coding scheme: the activity of the spiking neurons. This property is key for
reducing computation costs, thus the energy consumption in networks of spiking
neurons. The activity of a neural network is defined as the average number of
spikes generated by each neuron during $T$ timesteps. Referring to the
previous image quantization example, the activity is defined as follows:

$\theta=\frac{\displaystyle\sum_{i=0}^{n}\sum_{j=0}^{m}z_{i,j}(t)}{n\times m}$
(9)

Where $(n,m)$ is the size of the input image. Summing over $z$, which is a
$T\times n\times m$ binary matrix, results in the total number of spikes
generated by all the LIF neurons. The activity, $\theta$, of the rate coding
and LIF conversion scheme can be seen on the right side of Fig. 2. This figure
shows that, like in each rate conversion scheme, the number of spikes
increases almost linearly with $T$. However, as we can observe from the PSNR
curve shown in Fig. 2, the amount of information carried by each new spike in
the neural conversion scheme saturates above $T\sim 10$. The process of coding
dense information into stream of spikes is key for SNN and has been pointed
out as one of the main reasons for the current performance gap between SNN and
ANN. In the next section we investigate how these rules and properties show up
in larger and complex networks of spiking neurons.

##  3 Image denoising with spiking neurons

Our study of image denoising with spiking neurons is based on the DnCNN
network proposed in [1]. We focus the study on Gaussian denoising with a
certain noise level. The considered network is composed of 17 convolutional
layers. Activation functions are replaced with LIF neurons ($V_{th}=1$,
$\tau=2$). The input of the network is a noisy observation $y=x+v$ where $v$
is additive white Gaussian noise with zero mean and standard deviation
$\sigma$. As proposed in [1] we follow a residual learning formulation, to
train a mapping $R(y)\sim v$ and recover the clean image from $x=y-R(y)$. The
averaged mean square error between the desired and the estimated residuals is
used as the loss function. Training data for gaussian gray image denoising are
generated using the same method proposed in [1]. A dataset composed of 12
images (not included in the train set) is used for testing. The surrogate
gradient [8] learning is used to train the SNN networks. Denoising results and
neuron activity are shown in Fig. 3. As we can observe the network performance
(PSNR) follows the same trend observed for information coding shown in Fig. 2.
The LIF conversion scheme can provide competitive denoising performance with
few timesteps, but PSNR saturates when $T>10$. Rate coding could theoretically
achieve the same PSNR as DnCNN but at the cost of hundreds of timesteps.

Fig. 3: Denoising PSNR (left) and $\theta$ (right) as a function of the number
of conversion timesteps for the LIF neural conversion scheme and the rate
coding. The dotted line represent the performance of DnCNN for a noise level
$\sigma=25$.

Fig. 4 illustrates the visual results of the coding methods on the C.man
image. As it can be seen with only 7 timesteps, the LIF conversion method
provides better images compared with rate coding. The latter scheme would
require a large amount of timesteps to encode analog values with the precision
needed for the denoising task.

(a) Clean image

(b) Noisy image (20.18 dB)

(c) LIF (28.35 dB)

(d) Rate coding (21.59 dB)

Fig. 4: Denoising results of the image “C.man” with noise level 25. Denoised
images are shown in Fig. c and d for LIF conversion and rate coding, both with
$T=7$.

##  4 Conclusion

In this paper we have presented an analysis, based on information coding, for
SNNs and its application for image denoising. Our analysis showed that
information coding at the neuron level can explain the performance at the
network level. As future work, we aim at using our approach to guide the
design of spiking neural models. Our objective is to encode information with
both low latency and high precision for further hardware neuromorphic
implementation.

## References

  * [1] Kai Zhang, Wangmeng Zuo, Yunjin Chen, Deyu Meng, and Lei Zhang,  “Beyond a Gaussian Denoiser: Residual Learning of Deep CNN for Image Denoising,”  IEEE Transactions on Image Processing, vol. 26, no. 7, pp. 3142–3155, July 2017,  arXiv: 1608.03981.
  * [2] Shuhang Gu, Lei Zhang, Wangmeng Zuo, and Xiangchu Feng,  “Weighted Nuclear Norm Minimization with Application to Image Denoising,”  in 2014 IEEE Conference on Computer Vision and Pattern Recognition, Columbus, OH, USA, June 2014, pp. 2862–2869, IEEE.
  * [3] Nassim Abderrahmane, Edgar Lemaire, and Benoît Miramond,  “Design space exploration of hardware spiking neurons for embedded artificial intelligence,”  Neural Networks, vol. 121, pp. 366–386, 2020.
  * [4] Javier Mendez, Kay Bierzynski, MP Cuéllar, and Diego P Morales,  “Edge intelligence: Concepts, architectures, applications and future directions,”  ACM Transactions on Embedded Computing Systems (TECS), 2022.
  * [5] Peter U. Diehl, Daniel Neil, Jonathan Binas, Matthew Cook, Shih-Chii Liu, and Michael Pfeiffer,  “Fast-classifying, high-accuracy spiking deep networks through weight and threshold balancing,”  in 2015 International Joint Conference on Neural Networks (IJCNN), July 2015, pp. 1–8,  ISSN: 2161-4407.
  * [6] Bodo Rueckauer, Iulia-Alexandra Lungu, Yuhuang Hu, Michael Pfeiffer, and Shih-Chii Liu,  “Conversion of Continuous-Valued Deep Networks to Efficient Event-Driven Networks for Image Classification,”  Frontiers in Neuroscience, vol. 11, pp. 682, 2017.
  * [7] Abhronil Sengupta, Yuting Ye, Robert Wang, Chiao Liu, and Kaushik Roy,  “Going Deeper in Spiking Neural Networks: VGG and Residual Architectures,”  Frontiers in Neuroscience, vol. 13, pp. 95, 2019.
  * [8] Emre O. Neftci, Hesham Mostafa, and Friedemann Zenke,  “Surrogate Gradient Learning in Spiking Neural Networks: Bringing the Power of Gradient-Based Optimization to Spiking Neural Networks,”  IEEE Signal Processing Magazine, vol. 36, no. 6, pp. 51–63, Nov. 2019.
  * [9] Stanisław Woźniak, Angeliki Pantazi, Thomas Bohnstingl, and Evangelos Eleftheriou,  “Deep learning incorporating biologically inspired neural dynamics and in-memory computing,”  Nature Machine Intelligence, vol. 2, no. 6, pp. 325–336, June 2020\.
  * [10] Effrosyni Doutsi, Lionel Fillatre, Marc Antonini, and Panagiotis Tsakalides,  “Dynamic Image Quantization Using Leaky Integrate-and-Fire Neurons,”  IEEE Transactions on Image Processing, vol. 30, pp. 4305–4315, 2021\.
  * [11] S. Roth and M.J. Black,  “Fields of Experts: a framework for learning image priors,”  in 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), June 2005, vol. 2, pp. 860–867 vol. 2,  ISSN: 1063-6919.
  * [12] Wenzhe Guo, Mohammed E. Fouda, Ahmed M. Eltawil, and Khaled Nabil Salama,  “Neural Coding in Spiking Neural Networks: A Comparative Study for Robust Neuromorphic Systems,”  Frontiers in Neuroscience, vol. 15, pp. 638474, Mar. 2021.
  * [13] Nassim Abderrahmane and Benoit Miramond,  “Neural coding: adapting spike generation for embedded hardware classification,”  in 2020 International Joint Conference on Neural Networks (IJCNN), Glasgow, United Kingdom, July 2020, pp. 1–8, IEEE.
  * [14] Romain Brette,  “Philosophy of the Spike: Rate-Based vs. Spike-Based Theories of the Brain,”  Frontiers in Systems Neuroscience, vol. 9, pp. 151, 2015.