[1]Hanying Deng [3]Nicolae C. Panoiu [4] Fangwei Ye

[1]School of Optoelectronic Engineering, Guangdong Polytechnic Normal
University, Guangzhou, 510665, China

2]Department of Physics, Changzhi University, Shanxi, 046011, China

[3]Department of Electronic and Electrical Engineering, University College
London, Torrington Place, London, WC1E 7JE, United Kingdom

[4]School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai,
200240, China

# Topologically enhanced nonlinear optical response of graphene nanoribbon
heterojunctions

<EMAIL_ADDRESS>Zhihao Qu  Yingji He  Changming Huang
<EMAIL_ADDRESS><EMAIL_ADDRESS>*  [  *  *

###### Abstract

We study the nonlinear optical properties of heterojunctions made of graphene
nanoribbons (GNRs) consisting of two segments with either the same or
different topological properties. By utilizing a quantum mechanical approach
that incorporates distant-neighbor interactions, we demonstrate that the
presence of topological interface states significantly enhances the second-
and third-order nonlinear optical response of GNR heterojunctions that are
created by merging two topologically inequivalent GNRs. Specifically, GNR
heterojunctions with topological interface states display third-order harmonic
hyperpolarizabilities that are more than two orders of magnitude larger than
those of their similarly sized counterparts without topological interface
states, whereas the second-order harmonic hyperpolarizabilities exhibit a more
than ten-fold contrast between heterojunctions with and without topological
interface states. Additionally, we find that the topological state at the
interface between two topologically distinct GNRs can induce a noticeable red-
shift of the quantum plasmon frequency of the heterojunctions. Our results
reveal a general and profound connection between the existence of topological
states and an enhanced nonlinear optical response of graphene nanostructures
and possible other photonic systems.

##  1 Introduction

The topology of matter has become a subject of extensive research interest due
to its ability to explain many intriguing physical phenomena, such as
topological insulators [1, 2], the quantum spin Hall effect [3, 4], and
Majorana fermions [5] and their analogs. In particular, interfaces between two
materials of different topological classes give rise to highly robust
topological states [6, 7, 8]. Since the first discovery of topological states
in condensed matter physics, these states have been observed in a variety of
systems spanning different areas of science, including mechanics [9, 10, 11],
acoustics [12, 13, 14], and optics [15, 16, 17, 18, 19]. Recently, localized
topological states have been achieved at the interface between two segments of
graphene nanoribbons (GNRs) with different topological phases [20, 21, 22].
The topological phases of GNRs are determined by their width, edge, and end
termination [20, 21]. With recent advances in bottom-up synthesis from
precursor molecules, atomically precise GNRs with various widths, well-defined
edges, and end terminations can be readily produced [23, 24, 25, 26],
providing a rich variety of graphene systems suitable for further development
of possible applications.

Despite the extensive study and rapid progress in topological physics, most of
existing studies have focused on the linear regime. However, recent research
suggests that the interplay between topology and nonlinearity can give rise to
important novel phenomena [27, 28, 29, 30]. For instance, a topological source
of quantum light has been realized in a two-dimensional array of ring
resonators [27], whereas four-wave mixing of topological edge plasmons has
been observed in graphene metasurfaces [28]. In addition, strong third-
harmonic generation has been experimentally demonstrated in a SSH-like left-
handed transmission line [29]. Moreover, GNRs inherit the intrinsically large
nonlinearity from extended graphene sheets and produce large field enhancement
via plasmon excitation [31, 32, 33]. Consequently, GNR heterojunctions
composed of two topologically inequivalent GNRs provide an ideal and versatile
platform to investigate the interplay between topology and optical
nonlinearity at the nanoscale.

When the size of graphene nanostructures is less than $\sim$ 10 nm, the widely
used classical nonlinear conductivity is inadequate to describe their
nonlinear optical properties. This is so because nonlocal and finite-size
effects play a crucial role in determining their optical response,
necessitating full-quantum calculations [34, 35, 36, 37, 38]. One widely used
quantum mechanical method to describe the electronic states of graphene
nanostructures is the tight-binding model [39, 40]. However, the standard
tight-binding model considers only the coupling of nearest-neighbor atoms, so
that the accuracy of its predictions can be improved. Recently, we proposed a
distant-neighbor quantum mechanical (DNQM) approach to compute the linear and
nonlinear optical properties of graphene nanostructures [37, 38]. Compared to
the tight-binding model, the DNQM calculations are more accurate, as it
includes the coupling of the $\pi$-orbital electron of each atom with the core
potential of all atoms in the graphene nanostructure. In particular, our DNQM
method has been applied to calculate the second- and third-order nonlinear
optical response of certain graphene nanostructures [38].

In this work, we use DNQM calculations to investigate the nonlinear optical
properties of GNR heterojunctions composed of two segments of GNRs with either
the same or different topological phases. We show that due to the emergence of
topological interface states, the nonlinear optical response of GNR
heterojunctions consisting of two topologically distinct GNRs is significantly
enhanced. More specifically, compared to GNR heterojunctions without
topological interface states, the third-order nonlinear hyperpolarizability of
GNR heterojunctions with topological interface states can increase by more
than two orders of magnitude. Additionally, second-harmonic generation in GNR
heterojunctions with topological interface states is shown to be more than 10
times larger than that in topologically trivial heterojunctions with similar
size. We also find that the resonance frequency of quantum plasmons in GNR
heterojunctions is red-shifted when the topological interface states are
generated.

##  2 System configuration

The GNR heterojunctions we considered are shown in Fig. 1. They consist of two
armchair-edged GNRs ($W$-AGNRs) that may or may not be topologically
equivalent. Here, $W$ represents the number of rows of carbon atoms that form
the width of GNRs. The topology of GNRs is determined by their width, edge,
and end termination, and can be characterized by a $Z_{2}$ invariant [20, 21].
This invariant is an integer and takes a value of $Z_{2}=1$ ($Z_{2}=0$) for
topologically nontrivial (trivial) ribbons with a unit cell that possesses
spatial inversion and/or mirror symmetry. We first focus on the AGNRs with
zigzag, zigzag′, and bearded terminations, as they are the most commonly
synthesized forms of GNRs. For odd $W$, the values of $Z_{2}$ of $W$-AGNRs are
given by [20]:

$\displaystyle
Z_{2}=\frac{1+(-1)^{[\frac{W}{3}]+[\frac{W+1}{2}]}}{2}\quad\mathrm{zigzag-
terminated},$ (1a) $\displaystyle
Z_{2}=\frac{1-(-1)^{[\frac{W}{3}]+[\frac{W+1}{2}]}}{2}\quad\mathrm{zigzag^{\prime}-terminated}.$
(1b)

Here, $[x]$ is the floor function that takes the maximum integer smaller than
or equal to the real number $x$. For even $W$, the values of $Z_{2}$ of AGNRs
are [20]:

$\displaystyle
Z_{2}=\frac{1-(-1)^{[\frac{W}{3}]+[\frac{W+1}{2}]}}{2}\quad\mathrm{zigzag-
terminated},$ (2a) $\displaystyle
Z_{2}=\frac{1-(-1)^{[\frac{W}{3}]}}{2}\quad\mathrm{bearded-terminated}.$ (2b)

Knowing the value of the $Z_{2}$ invariant of $W$-AGNRs with various
terminations and widths, we construct four heterojunctions lying on the
$x$-$y$ plane, at which two AGNRs of same or distinct topology are joined
together, as shown in Fig. 1.

Figure 1: Armchair GNR heterojunctions formed with two GNRs of (a) W=9 zigzag
and $W=7$ zigzag′ terminations, (b) $W=9$ zigzag and $W=7$ zigzag
terminations, (c) $W=6$ and $W=4$ bearded terminations, and (d) $W=6$ zigzag
and $W=4$ zigzag terminations. $W$ is the number of rows of carbon atoms along
the lateral direction. The two GNR segments in (a) and (c) are topologically
inequivalent, whereas those in (b) and (d) are topologically equivalent. The
dashed red rectangles indicates a unit cell.

The two heterojunctions presented in Figs. 1(a) and 1(b) consist of a zigzag-
terminated $W=9$ AGNR with $Z_{2}=1$ and a zigzag′-terminationed $W=7$ AGNR
with $Z_{2}=0$ [Fig. 1(a)], and a zigzag-terminated $W=7$ AGNR with $Z_{2}=1$
[Fig. 1(b)]; we labeled them as 9/7-Zigzag-8D and 9/7-Zigzag-8S, respectively.
Note that the two heterojunctions contain the same number of carbon atoms,
$N_{c}=256$, and have the same length of $L=6.7$ nm. We also consider two
heterojunctions containing the same number of carbon atoms, $N_{c}=240$, and
having the same side length of $L=10.1$ nm by connecting two bearded- or
zigzag-terminated AGNRs with even rows of carbon atoms, $W=4$ and 6, and
denote them as 6/4-Bearded-12D and 6/4-Zigzag-12S, as shown in Figs. 1(c) and
1(d), respectively. The two AGNR segments in Fig. 1(c) are topologically
inequivalent, while those in Fig. 1(d) are topologically equivalent.

##  3 Energy spectra and charge density distributions

Figure 2: Energy spectra of (a) 9/7-Zigzag-8D and (c) 6/4-Bearded-12D
heterojunctions. A topological interface state exists at the interface of (b)
9/7-Zigzag-8D and (d) 6/4-Bearded-12D heterojunctions. The topological
interface states in 9/7-Zigzag-8D and 6/4-Bearded-12D heterojunctions are
marked by the red dots in (a) and (c), respectively.

We first investigate energy spectra and charge density distributions of the
AGNR heterojunctions using the DNQM method (details about the DNQM
calculations can be found in Appendix A). Figures 2(a) and 2(c) depict the
energy spectra of the 9/7-Zigzag-8D and 6/4-Bearded-12D heterojunctions,
respectively, with modes arranged in increasing order of their enenergy. We
also calculate their corresponding charge density distributions. In accordance
with the bulk-edge correspondence principle, our computations demonstrate that
a topologically induced interface state arises at the 9/7-Zigzag-8D
heterojunction, localized at the interface between the two topologically
distinct AGNR segments, as shown in Fig. 2(b). Remarkably, the energy of the
topological interface state in the 9/7-Zigzag-8D heterojunction corresponds to
the HOMO-1 (with HOMO meaning the highest-occupied molecular orbital), having
a value of $E=-154$ eV. The red dot in the energy spectrum in Fig. 2(a)
indicates this energy level.

Similarly, Fig. 2(d) illustrates a localized interface state of the
6/4-Bearded-12D heterojunction, due to the fact that the two GNR segments
possess different topological phases. We also note that the energy of the
topological interface state of the 6/4-Bearded-12D heterojunction corresponds
to the HOMO, with a value of -118.3 eV, as indicated by the red dot in the
energy spectrum in Fig. 2(c). By contrast, no localized interface states exist
at GNR heterojunctions formed with two topologically equivalent GNR segments
(9/7-Zigzag-8S and 6/4-Zigzag-12S).

##  4 Linear and nonlinear optical response

We next investigate the linear and nonlinear optical response of the AGNR
heterojunctions composed of two GNRs that are topologically equivalent and
inequivalent. We first consider the 9/7-Zigzag-8D and 9/7-Zigzag-8S
heterojunctions. As discussed above, the 9/7-Zigzag-8D heterojunction has a
topological state located at the interface separating its two topologically
inequivalent AGNR segments, whereas no localized interface states exist in the
9/7-Zigzag-8S heterojunction. Additionally, the two AGNR heterojunctions have
the same length and contain the same number of carbon atoms. Notice that,
because its symmetry to reflection w.r.t. the $x-z$ plane, the second-order
nonlinear optical response of the 9/7-Zigzag-8D heterojunction is not allowed
for $y$-polarized incident fields. Thus, we assume that the incident electric
field is $x$-polarized. The results of the calculations for the linear
polarizability and second- and third-order hyperpolarizabilities of
9/7-Zigzag-8D and 9/7-Zigzag-8S heterojunctions are summarized in Fig. 3 (see
calculation details in Appendix B).

Figure 3: Real and imaginary parts of [(a), (b)] linear polarizability, [(c),
(d)] second-order, and [(e), (f)] third-order hyperpolarizabilities of
9/7-Zigzag-8D and 9/7-Zigzag-8S heterojunctions, respectively. The red (green)
curves indicate the imaginary (real) part of the polarizabilities. The purple
dots denote the resonance peaks. The polarizabilities are given in atomic
units (a.u.) with $e=\hbar=m_{e}=a_{0}=1$.

In Figs. 3(a) and 3(b), the real and imaginary parts of the linear
polarizability spectra, $\alpha_{xx}(\omega)$, of 9/7-Zigzag-8D and
9/7-Zigzag-8S heterojunctions are shown, respectively. Quantum plasmons exist
when the peak of the imaginary parts of the polarizability occurs at a zero of
the real part[36, 37, 41]. From Figs. 3(a) and 3(b), we can see that the first
resonance frequency of quantum plasmon in the topologically nontrivial
heterojunction is substantially smaller than that in the topologically trivial
heterojunction. Specifically, the 9/7-Zigzag-8D and 9/7-Zigzag-8S
heterojunctions have the first resonance peak located at 0.73 eV and 1.44 eV,
respectively. This feature can be explained by the molecular orbital
associated with its topological interface state participating in the
transition. Thus, the first resonance peak in the linear polarizabilities of
the 9/7-Zigzag-8D heterojunction is located at ${\hbar}\omega=0.73$ eV, which
corresponds to the transition from the HOMO-1 state (which as discussed above
is the topological state) to the LUMO+2.

Figures 3(c) and 3(d) show the real and imaginary parts of the second-order
hyperpolarizability, $\beta_{xxx}(2\omega)$, corresponding to the second-
harmonic generation (SHG), of 9/7-Zigzag-8D and 9/7-Zigzag-8S heterojunctions,
respectively. Remarkably, the SHG hyperpolarizability spectra of these two
heterojunctions are significantly different. Firstly, it can be seen that
$\beta_{xxx}(2\omega)$ of topologically nontrivial 9/7-Zigzag-8D
heterojunction is more than one order of magnitude larger than that of
topologically trivial 9/7-Zigzag-8S heterojunction. Additionally, similar to
the linear case, the resonance frequency of the second-order nonlinear plasmon
of 9/7-Zigzag-8D heterojunction is red-shifted in comparison to that of
9/7-Zigzag-8S heterojunction. These two interesting characteristics of the
topologically nontrivial 9/7-Zigzag-8D heterojunction can be ascribed to its
topological interface state participating in the transitions. In particular,
as shown in Fig. 3(c), the most pronounced peak in the spectrum of SHG
hyperpolarizabilities of 9/7-Zigzag-8D heterojunction is located at
${\hbar}\omega=0.027$ eV, which corresponds to the transition from the HOMO-1,
corresponding to the topological interface state, to the LUMO, indicating that
the presence of the topological interface state can dramatically enhance the
second-order nonlinear response of the 9/7-Zigzag-8D heterojunction and induce
a red-shift of the resonance frequency of its quantum plasmons.

In Figs. 3(e) and 3(f), we show the spectra of the third-harmonic generation
(THG) hyperpolarizability, $\gamma$${}_{xxxx}(3\omega)$, of 9/7-Zigzag-8D and
9/7-Zigzag-8S heterojunctions, respectively. Similar to the case of the
second-order nonlinear optical response, due to the presence of the
topological interface state, the THG of the topologically nontrivial
9/7-Zigzag-8D heterojunction is dramatically enhanced and its plasmon
frequency is red-shifted. Interestingly, we find that the hyperpolarizability
$\gamma$${}_{xxxx}(3\omega)$ of the topologically nontrivial 9/7-Zigzag-8D
heterojunction is enhanced by more than two orders of magnitude in comparison
to that of the topologically trivial 9/7-Zigzag-8S heterojunction.

Figure 4: Real and imaginary parts of [(a), (b)] linear, [(c), (d)] second-
order, and [(e),(f)] third-order hyperpolarizability of 6/4-Bearded-12D and
6/4-Zigzag-12S heterojunctions. [(g), (h)] Evolution of the peak value of (g)
SHG and (h) THG hyperpolarizabilities with the side length of 6/4-Bearded-12D
(blue curves) and 6/4-Zigzag-12S (purple curves) heterojunctions,
respectively. The red (green) curves indicate the imaginary (real) part of the
hyperpolarizabilities. The purple dots denote the resonance peaks.

We now consider the 6/4-Bearded-12D and 6/4-Zigzag-12S heterojunctions and use
the DNQM approach to investigate their linear and nonlinear optical response,
the corresponding results being summarized in Fig. 4. As explained above, the
6/4-Bearded-12D heterojunction possesses a topological interface state, while
no topological interface states exist in the 6/4-Zigzag-12S heterojunction.
The real and imaginary parts of linear polarizability, $\alpha_{xx}(\omega)$,
of 6/4-Bearded-12D and 6/4-Zigzag-12S heterojunctions are presented in Figs.
4(a) and 4(b), respectively. Again, due to the emergence of the topological
interface state, the most pronounced peak in the spectrum of
$\alpha_{xx}(\omega)$ of 6/4-Bearded-12D heterojunction is significantly red-
shifted in comparison to that of 6/4-Zigzag-12S heterojunction. The nonlinear
hyperpolarizabilities associated with SHG ($\beta_{xxx}(2\omega)$) and THG
($\gamma$${}_{xxxx}(3\omega)$) are presented in Figs. 4(c,d) and 4(e,f),
respectively, for the 6/4-Bearded-12D and 6/4-Zigzag-12S heterojunctions.
Comparing these spectra, we can draw similar conclusions to those revealed by
Fig. 3. More specifically, due to the existence of the topological interface
state, both the SHG and THG hyperpolarizabilities of the 6/4-Bearded-12D
heterojunction exceed by more than one order of magnitude those of the
topologically trivial 6/4-Zigzag-12S heterojunction. Furthermore, the presence
of topological interface states also leads to an obvious red-shift of the
nonlinear plasmon frequency of 6/4-Bearded-12D heterojunction.

In Figs. 4(g) and 4(h), we present an overview of the dependence of the
maximum SHG and THG hyperpolarizabilities on the side length of the
6/4-Bearded-12D and 6/4-Zigzag-12S heterojunctions, respectively. We can see
that the SHG and THG hyperpolarizabilities of 6/4-Bearded-12D heterojunction
that possesses a topological interface state are more than one order of
magnitude larger than those of the topologically trivial 6/4-Zigzag-12S
heterojunction of similar side length. As expected, the magnitudes of the SHG
and THG hyperpolarizabilities of both 6/4-Bearded-12D and 6/4-Zigzag-12S
heterojunctions increase with the increase of their side length. Interestingly
enough, we note that the rate of the increase of the nonlinear
hyperpolarizabilities with side length for the topologically nontrivial
6/4-Bearded-12D heterojunction is larger than that for the topologically
trivial 6/4-Zigzag-12S heterojunction.

Figure 5: GNR heterojunctions formed with (a) a zigzag′-terminationed $W=7$
AGNR and a $P=5$ cove-edged GNR and (b) a zigzag-terminated $W=7$ AGNR and a
$P=5$ cove-edged GNR. The two GNR segments in (a) and (b) are topologically
inequivalent and equivalent, respectively. Charge density distributions of (c)
7AGNR/5cove-8D and (d) 7AGNR/5cove-8S. Real and imaginary parts of [(e), (f)]
linear polarizability and [(g), (h)] second-order, and [(i), (j)] third-order
hyperpolarizability of 7AGNR/5cove-8D and 7AGNR/5cove-8S heterojunctions. The
red (green) curves indicate the imaginary (real) part of the polarizabilities.
The purple dots denote the resonance peaks.

Finally, we focus on a more complex GNR heterojunction composed of an AGNR and
a cove-edged GNR, as schematically shown in Figs. 5(a) and 5(b). Cove-edged
GNRs have been recently synthesized by bottom-up process from precursor
molecules [42]. Their topological properties are determined by the ribbon
width, edge shape, and end termination [21]. Figures 5(a) and 5(b) present two
heterojunctions consisting of a $P=5$ cove-edged GNR with tilting angle of
$30^{\text{o}}$ ($Z_{2}=1$) and either a zigzag′-terminationed $W=7$ AGNR
($Z_{2}=0$) or a zigzag-terminated $W=7$ AGNR ($Z_{2}=1$), which are denoted
by 7AGNR/5cove-8D and 7AGNR/5cove-8S, respectively. Here, $P$ denotes the
number of zigzag chains forming the width of the cove-edged GNR. The two
heterojunctions contain the same number of carbon atoms, $N_{c}=256$. Again,
by bulk-boundary correspondence principle, a topological interface state
occurs for the 7AGNR/5cove-8D heterojunction and is localized at the interface
between the two GNRs of distinct topology, as shown in Fig. 5(c). By contrast,
for the 7AGNR/5cove-8S heterojunction comprised of two topologically
equivalent GNRs, no localized interface states exist and the corresponding
charge profile is extended throughout the heterojunction, as shown in Fig.
5(d).

We show in Figs. 5(e-j) the frequency dependence of linear polarizability, and
SHG and THG hyperpolarizabilities of the 7AGNR/5cove-8D and 7AGNR/5cove-8S
heterojunctions, respectively. As expected, due to the occurrence of the
topological interface state, the resonance frequency of both the linear and
nonlinear plasmons of 7AGNR/5cove-8D heterojunction is significantly red-
shifted, as compared with that of 7AGNR/5cove-8S heterojunction. Furthermore,
the magnitudes of the second- and third-order nonlinear hyperpolarizabilities
of the topologically nontrivial 7AGNR/5cove-8D heterojunction are
significantly larger than those of the topologically trivial 7AGNR/5cove-8D
heterojunction, as indicted by the red dots in Figs. 5(g-j), respectively.

##  5 Conclusion

In conclusion, we have utilized the DNQM approach to investigate the optical
response of GNR heterojunctions, comprising of two topologically equivalent or
inequivalent ribbons. Our findings reveal a significant enhancement in the
nonlinear optical response of GNR heterojunctions that possess topological
interface states as compared to the trivial ones. Specifically, we have
observed an increase in the third-order hyperpolarizabilities of GNR
heterojunctions with topological interface states by over two orders of
magnitude in comparison to those without such states, despite having the same
number of carbon atoms. Moreover, due to the presence of topological interface
states, the second-order hyperpolarizabilities of GNR heterojunctions with
topologically distinct GNRs are over one order of magnitude larger than those
of topologically trivial heterojunctions. Furthermore, we have noticed a large
red-shift in the resonance frequency of quantum plasmons in GNR
heterojunctions that contain topological states. Remarkably, these results
were established for several heterojunctions with different configurations,
which suggests that there is a general and deep connection between the
existence of topological states and enhanced nonlinear optical response in
GNRs and other photonic systems.

## Appendix A: Calculation of energy spectra and charge density distributions
of GNR heterojunctions

We use a DNQM approach to calculate the energy spectra and charge density
distributions of GNR heterojunctions shown in Fig.1 and Fig.5. Our
calculations include the interactions between $\pi$ electrons and all cores of
carbon atoms in the GNR heterojunctions. It goes beyond the tight-binding
approximation used successfully to determine the linear and nonlinear optical
response of graphene nanoflakes [31, 32, 34], as in our case the coupling
constants are calculated from first principles. Specifically, we assume that
the $\pi$ electrons reside in the $2p_{z}$ carbon orbitals oriented
perpendicular to the GNR heterojunction plane and populate on average with one
such electron per carbon site. Thus, we employ the Clementi orbital [43] as a
basis function on each atom, which takes the form of

$\displaystyle\psi_{2p_{z}}=R_{2p}(r)Y_{2p_{z}}(\theta,\varphi),$ (A1)

with

$\displaystyle
R_{2p}(r)=\frac{1}{\sqrt{6}}\frac{Zr}{na_{0}}Z^{\frac{3}{2}}e^{-\frac{Zr}{na_{0}}},$
(A2)

and

$\displaystyle
Y_{2p_{z}}(\theta,\varphi)=\sqrt{\frac{3}{4\pi}}\cos\theta=\sqrt{\frac{3}{4\pi}}\frac{z}{r}.$
(A3)

Here, $a_{0}$ is the Bohr radius, $n=2$ is the orbital number, and $Z=3.136$
is the effective nuclear charge for the $2p_{z}$ orbital of a carbon atom[43].

The Hamiltonian operator for a single electron in a GNR heterojunction is

$\displaystyle\widehat{H}=-\frac{\hbar}{2m}\nabla_{\vec{r}}^{2}-\sum_{\alpha=1}^{N}\frac{Z_{\text{eff}}e^{2}}{\vec{r}-\vec{r}_{0\alpha}},$
(A4)

where the value of effective core charge, $Z_{\text{eff}}=0.637$, has been
adjusted so that the computed HOMO-LUMO gap of the simplest GNR structure
(i.e. benzene) is 6 eV [40].

The energy spectra and charge density distributions are calculated using the
the Schrödinger equation

$\displaystyle\widehat{H}\psi(\vec{r})=E\psi(\vec{r}),$ (A5)

where $E$ and $\psi(\vec{r})$ are the eigenenergy and eigenfunction,
respectively. This eigenfunction is expanded in the atom-centered basis
functions as

$\displaystyle\psi(\vec{r})=\sum_{q=1}^{N}c_{q}\psi_{q}(\vec{r})=\sum_{q=1}^{N}c_{q}\psi_{2p_{z}}(\vec{r}-\vec{r}_{oq}).$
(A6)

From Eq. (A5) and Eq. (A6), we can have

$\displaystyle\sum_{q=1}^{N}\widehat{H}c_{q}\psi_{q}=E\sum_{q=1}^{N}c_{q}\psi_{q}.$
(A7)

We then use $\psi_{s}^{*}(\vec{r})$, $s=1,...,N$, to multiply Eq. (A7), and
integrate over the entire space:

$\displaystyle\sum_{q=1}^{N}c_{q}\int\psi_{s}^{*}\widehat{H}\psi_{q}d\vec{r}=E\sum_{q=1}^{N}c_{q}\int\psi_{s}^{*}\psi_{q}d\vec{r}.$
(A8)

Now, we define the matrix elements

$\displaystyle\widehat{H}_{sq}$
$\displaystyle=\int\psi_{s}^{*}\hat{H}\psi_{q}d\vec{r},$ (A9a)
$\displaystyle\widehat{S}_{sq}$
$\displaystyle=\int\psi_{s}^{*}\psi_{q}d\vec{r},$ (A9b)

so that the Schrödinger equation is reduced to a finite generalized eigenvalue
problem with eigenvectors, $\widehat{c}$, and eigenvalues, $E$

$\displaystyle\widehat{H}\hat{c}=E\widehat{S}\hat{c},$ (A10)

where $\hat{c}=(c_{1},c_{2},\ldots,c_{N})^{T}$ is an $N$-dimensional column
vector. We can obtain the energy spectra and charge density distributions of a
GNR heterojunction containing $N$ carbon atoms by numerically solving Eq.
(A10).

## Appendix B: Quantum perturbative approach to the linear and nonlinear
optical polarizabilities of GNR heterojunctions

We calculate linear and nonlinear optical polarizabilities of GNR
heterojunctions using a well-known quantum perturbative approach [44]. The
linear polarizability, $\alpha_{ij}(\omega)$, and the nonlinear
hyperpolarizabilities corresponding to second-harmonic generation,
$\beta_{ijk}(2\omega)$, and third-harmonic generation,
$\gamma_{hijk}(3\omega$) are given by [44]

$\displaystyle\alpha_{ij}(\omega)=\sum_{g,m}^{transition}\Big{(}\frac{\mu_{gm}^{i}\mu_{mg}^{j}}{\omega_{mg}-\omega}+\frac{\mu_{gm}^{j}\mu_{mg}^{i}}{\omega_{mg}^{*}+\omega}\Big{)},$
(B1)

$\displaystyle\beta_{ijk}(2\omega)=P_{l}\sum_{g,m,n}^{transition}\Big{[}\frac{\mu_{gn}^{i}\mu_{nm}^{j}\mu_{mg}^{k}}{(\omega_{ng}-2\omega)(\omega_{mg}-\omega)}$
$\displaystyle+\frac{\mu_{gn}^{j}\mu_{nm}^{i}\mu_{mg}^{k}}{(\omega_{ng}^{*}+\omega)(\omega_{mg}-\omega)}$
$\displaystyle+\frac{\mu_{gn}^{j}\mu_{nm}^{k}\mu_{mg}^{i}}{(\omega_{ng}^{*}+\omega)(\omega_{mg}^{*}+2\omega)}\Big{]},$
(B2)

$\displaystyle\gamma_{hijk}((3\omega)=$ $\displaystyle
P_{l}\sum_{g,m,n,v}^{transition}\Big{[}\frac{\mu_{gv}^{h}\mu_{vn}^{i}\mu_{nm}^{j}\mu_{mg}^{k}}{(\omega_{vg}-3\omega)(\omega_{ng}-2\omega)(\omega_{mg}-\omega)}$
$\displaystyle+\frac{\mu_{gv}^{i}\mu_{vn}^{h}\mu_{nm}^{j}\mu_{mg}^{k}}{(\omega_{vg}^{*}+\omega)(\omega_{ng}-2\omega)(\omega_{mg}-\omega)}+\frac{\mu_{gv}^{i}\mu_{vn}^{j}\mu_{nm}^{h}\mu_{mg}^{k}}{(\omega_{vg}^{*}+\omega)(\omega_{ng}^{*}+2\omega)(\omega_{mg}-\omega)}$
$\displaystyle+\frac{\mu_{gv}^{i}\mu_{vn}^{j}\mu_{nm}^{k}\mu_{mg}^{h}}{(\omega_{vg}^{*}+\omega)(\omega_{ng}^{*}+2\omega)(\omega_{mg}^{*}+3\omega)}\Big{]},$
(B3)

where $\vec{\mu}=e\vec{r}$ is the dipole moment operator,
$\mu_{gm}=\int\psi_{g}^{*}\vec{\mu}\psi_{m}d\vec{r}$ is the transition dipole
moment, $\omega_{mg}=\frac{E_{m}-E_{g}}{\hbar}-i\eta$,
${\psi_{m}}=\sum_{l}c_{ml}{\psi_{l}}$, and $E_{m}$ and $c_{m}=\sum_{l}c_{ml}$
are the eigenenergy and engenvector of the eigenstate $m$, respectively. In
addition, $g$, $m$, $n$ and $v$ are labels used to distinguish between levels
of transitions, and $\eta=0.1$ eV is related to the lifetime of excited
states. $P_{l}$ is a intrinsic permutation operator that defines the
expression to be summed over all permutations of the Cartesian indices: $h$,
$i$, $j$, and $k$.

## References

  *   * Qi and Zhang [2011] Qi, X.-L., Zhang, S.-C.: Topological insulators and superconductors. Rev. Mod. Phys. 83, 1057–1110 (2011) https://doi.org/10.1103/RevModPhys.83.1057
  * Hasan and Kane [2010] Hasan, M.Z., Kane, C.L.: Colloquium: Topological insulators. Rev. Mod. Phys. 82, 3045–3067 (2010) https://doi.org/10.1103/RevModPhys.82.3045
  * Kane and Mele [2005a] Kane, C.L., Mele, E.J.: ${Z}_{2}$ topological order and the quantum spin hall effect. Phys. Rev. Lett. 95, 146802 (2005) https://doi.org/10.1103/PhysRevLett.95.146802
  * Kane and Mele [2005b] Kane, C.L., Mele, E.J.: Quantum spin hall effect in graphene. Phys. Rev. Lett. 95, 226801 (2005) https://doi.org/10.1103/PhysRevLett.95.226801
  * Kitaev [2001] Kitaev, A.Y.: Unpaired majorana fermions in quantum wires. Physics-Uspekhi 44, 131 (2001) https://doi.org/10.1070/1063-7869/44/10S/S29
  * Konig et al. [2007] Konig, M., Wiedmann, S., Brune, C., Roth, A., Buhmann, H., Molenkamp, L.W., Qi, X.-L., Zhang, S.-C.: Quantum spin hall insulator state in hgte quantum wells. Science 318(5851), 766–770 (2007)
  * Grusdt et al. [2013] Grusdt, F., Höning, M., Fleischhauer, M.: Topological edge states in the one-dimensional superlattice bose-hubbard model. Phys. Rev. Lett. 110, 260405 (2013) https://doi.org/10.1103/PhysRevLett.110.260405
  * Wang and Zhang [2017] Wang, J., Zhang, S.-C.: Topological states of condensed matter. Nature materials 16(11), 1062–1067 (2017)
  * Jawed et al. [2015] Jawed, M.K., Dieleman, P., Audoly, B., Reis, P.M.: Untangling the mechanics and topology in the frictional response of long overhand elastic knots. Phys. Rev. Lett. 115, 118302 (2015) https://doi.org/10.1103/PhysRevLett.115.118302
  * Huber [2016] Huber, S.D.: Topological mechanics. Nature Physics 12(7), 621–623 (2016)
  * Lee et al. [2018] Lee, C.H., Li, G., Jin, G., Liu, Y., Zhang, X.: Topological dynamics of gyroscopic and floquet lattices from newton’s laws. Phys. Rev. B 97, 085110 (2018) https://doi.org/10.1103/PhysRevB.97.085110
  * He et al. [2016] He, C., Ni, X., Ge, H., Sun, X.-C., Chen, Y.-B., Lu, M.-H., Liu, X.-P., Chen, Y.-F.: Acoustic topological insulator and robust one-way sound transport. Nature physics 12(12), 1124–1129 (2016)
  * Ye et al. [2022] Ye, L., Qiu, C., Xiao, M., Li, T., Du, J., Ke, M., Liu, Z.: Topological dislocation modes in three-dimensional acoustic topological insulators. Nature Communications 13(1), 508 (2022)
  * Yang et al. [2015] Yang, Z., Gao, F., Shi, X., Lin, X., Gao, Z., Chong, Y., Zhang, B.: Topological acoustics. Phys. Rev. Lett. 114, 114301 (2015) https://doi.org/10.1103/PhysRevLett.114.114301
  * Lu et al. [2014] Lu, L., Joannopoulos, J.D., S, M.: Topological photonics. Nature photonics 8(11), 821–829 (2014)
  * Shaolin et al. [2017] Shaolin, K., Bing, W., Hua, L., Kai, W., Peixiang, L.: Topological edge modes in non-hermitian plasmonic waveguide arrays. Optics Express 25(10), 11132–11143 (2017)
  * Ye et al. [2016] Ye, F., Deng, H., Panoiu, N.C., Chen, X.: Topological surface plasmons in superlattices with changing sign of the average permittivity. Optics Letters 41(18), 4281 (2016)
  * Ozawa et al. [2019] Ozawa, T., Price, H.M., Amo, A., Goldman, N., Hafezi, M., Lu, L., Rechtsman, M.C., Schuster, D., Simon, J., Zilberberg, O., Carusotto, I.: Topological photonics. Rev. Mod. Phys. 91, 015006 (2019) https://doi.org/10.1103/RevModPhys.91.015006
  * Deng et al. [2019] Deng, H., Chen, Y., Huang, C., Ye, F.: Topological interface modes in photonic superlattices containing negative-index materials. Europhysics Letters 124(6), 64001 (2019)
  * Cao et al. [2017] Cao, T., Zhao, F., Louie, S.G.: Topological phases in graphene nanoribbons: Junction states, spin centers, and quantum spin chains. Phys. Rev. Lett. 119, 076401 (2017) https://doi.org/10.1103/PhysRevLett.119.076401
  * Lee et al. [2018] Lee, Y.L., Zhao, F., Cao, T., Ihm, J., Louie, S.G.: Topological phases in cove-edged and chevron graphene nanoribbons: Geometric structures, $Z_{2}$ invariants, and junction states. Nano Lett. 18, 7247–7253 (2018)
  * Zhao et al. [2021] Zhao, F., Cao, T., Louie, S.G.: Topological phases in graphene nanoribbons tuned by electric fields. Phys. Rev. Lett. 127, 166401 (2021) https://doi.org/10.1103/PhysRevLett.127.166401
  * Cai et al. [2014] Cai, J., Pignedoli, C.A., Talirz, L., Ruffieux, P., Söde, H., Liang, L., Meunier, V., Berger, R., Li, R., Feng, X., et al.: Graphene nanoribbon heterojunctions. Nat. Nanotechnol. 9(11), 896–900 (2014)
  * Ruffieux et al. [2012] Ruffieux, P., Cai, J., Plumb, N.C., Patthey, L., Prezzi, D., Ferretti, A., Molinari, E., Feng, X., Müllen, K., Pignedoli, C.A., et al.: Electronic structure of atomically precise graphene nanoribbons. Acs Nano 6(8), 6930–6935 (2012)
  * Liu et al. [2020] Liu, Z., Chen, Z., Wang, C., Wang, H.I., Wuttke, M., Wang, X.-Y., Bonn, M., Chi, L., Narita, A., Müllen, K.: Bottom-up, on-surface-synthesized armchair graphene nanoribbons for ultra-high-power micro-supercapacitors. J. Am. Chem. Soc. 142(42), 17881–17886 (2020)
  * Narita et al. [2015] Narita, A., Feng, X., Müllen, K.: Bottom-up synthesis of chemically precise graphene nanoribbons. The Chemical Record 15(1), 295–309 (2015)
  * Mittal et al. [2018] Mittal, S., Goldschmidt, E.A., Hafezi, M.: A topological source of quantum light. Nature 561(7724), 502–506 (2018)
  * You et al. [2020] You, J.W., Lan, Z., Panoiu, N.C.: Four-wave mixing of topological edge plasmons in graphene metasurfaces. Sci. Adv. 6(13), 3910 (2020)
  * Wang et al. [2019] Wang, Y., Lang, L.-J., Lee, C.H., Zhang, B., Chong, Y.: Topologically enhanced harmonic generation in a nonlinear transmission line metamaterial. Nat. Commun. 10(1), 1102 (2019)
  * Smirnova et al. [2019] Smirnova, D., Kruk, S., Leykam, D., Melik-Gaykazyan, E., Choi, D.-Y., Kivshar, Y.: Third-harmonic generation in photonic topological metasurfaces. Phys. Rev. Lett. 123(10), 103901 (2019)
  * Cox et al. [2016] Cox, J.D., Silveiro, I., Abajo, F.J.: Quantum effects in the nonlinear response of graphene plasmons. ACS nano 10(2), 1995–2003 (2016)
  * Cox and Garcia de Abajo [2019] Cox, J.D., Abajo, F.J.: Nonlinear graphene nanoplasmonics. Acc. Chem. Res. 52(9), 2536–2547 (2019)
  * Hendry et al. [2010] Hendry, E., Hale, P.J., Moger, J., Savchenko, A., Mikhailov, S.A.: Coherent nonlinear optical response of graphene. Phys. Rev. Lett. 105(9), 097401 (2010)
  * Cox and Javier García de Abajo [2014] Cox, J.D., Abajo, F.: Electrically tunable nonlinear plasmonics in graphene nanoislands. Nat. Commun. 5(1), 5725 (2014)
  * Manjavacas et al. [2013] Manjavacas, A., Thongrattanasiri, S., Abajo, F.J.G.: Plasmons driven by single electrons in graphene nanoislands. Nanophotonics 2(2), 139–151 (2013)
  * Manrique et al. [2017] Manrique, D.Z., You, J.W., Deng, H., Ye, F., Panoiu, N.C.: Quantum plasmon engineering with interacting graphene nanoflakes. J. Phys. Chem. C 121(49), 27597–27602 (2017)
  * Deng et al. [2018] Deng, H., Manrique, D.Z., Chen, X., Panoiu, N.C., Ye, F.: Quantum mechanical analysis of nonlinear optical response of interacting graphene nanoflakes. APL Photon. 3(1), 016102 (2018)
  * Deng et al. [2021] Deng, H., Huang, C., He, Y., Ye, F.: Quantum plasmon enhanced nonlinear wave mixing in graphene nanoflakes. Chin. Phys. B 30(4), 044213 (2021)
  * Cox and Garcia de Abajo [2015] Cox, J.D., Abajo, F.J.: Plasmon-enhanced nonlinear wave mixing in nanostructured graphene. Acs Photon. 2(2), 306–312 (2015)
  * Ezawa [2007] Ezawa, M.: Metallic graphene nanodisks: Electronic and magnetic properties. Phys. Rev. B 76(24), 245415 (2007)
  * Tame et al. [2013] Tame, M.S., Mcenery, K.R., Özdemir, Ş., Lee, J., Maier, S.A., Kim, M.S.: Quantum plasmonics. Nat. Phys. 9(6), 329–340 (2013)
  * Liu et al. [2015] Liu, J., Li, B.-W., Tan, Y.-Z., Giannakopoulos, A., Sanchez-Sanchez, C., Beljonne, D., Ruffieux, P., Fasel, R., Feng, X., Müllen, K.: Toward cove-edged low band gap graphene nanoribbons. J. Am. Chem. Soc. 137(18), 6097–6103 (2015)
  * Clementi and Raimondi [1963] Clementi, E., Raimondi, D.L.: Atomic screening constants from SCF functions. J. Chem. Phys. 38, 2686–2689 (1963)
  * Boyd [2003] Boyd, R.W. (ed.): Nonlinear Optics, 3rd Ed. Taylor $\&$ Francis, (2003)


§.§ Detection and Equalization

Our metric of optimality is the  maximum  a posteriori probability.  Given the multipath induced by the DiMC, to achieve good performance, sequence detection will be necessary in many cases.  Under the assumption that each possible symbol is equally likely, this strategy reduces to maximum likelihood sequence detection (MLSD).  As famously shown in the seminal work,  complexity reduction is achievable if the memory/multipath is of finite duration and the Viterbi detector can be employed <cit.>. In a DiMC system, the LTI-Poisson channel model clearly reveals the inherent multipath and channel memory. Recalling that the transmission vector $\boldsymbol{x}$ is related to $\boldsymbol{s}$ through the employed modulation scheme, and given data is transmitted in blocks with block length $S$, the MLSD operates by performing

\begin{equation}
\label{eq:MLSD}
\begin{split}
\hat{\boldsymbol{s}} &= \argmax_{\boldsymbol{s}} P(\boldsymbol{y} | \boldsymbol{x}) \\
&= \argmax_{\boldsymbol{s}} \prod_{n=1}^{S N} \frac{ \Tilde{\lambda}[n]^{y[n]}     e^{-\Tilde{\lambda}[n] } } {(y[n])!} \\
&= \argmax_{\boldsymbol{s}} \sum_{n=1}^{S N} y[n] \ln (\Tilde{\lambda}[n]) - \Tilde{\lambda}[n].
\end{split}
\end{equation}

Here, according to the candidate symbol sequence $\boldsymbol{s}$ and the channel model presented in Section <ref>, the $SN$-vector $\boldsymbol{\Tilde{\lambda}}$ is expressed as

\begin{equation}
\Tilde{\lambda}[n] = \sum_{n=1}^{S N} \bar{x}[n-k+1] h[n].
\end{equation}

We note that the MLSD presented in (<ref>) is a generalized expression for an arbitrary (single-molecule) modulation scheme. That said, as mentioned in Section <ref>, we focus on OOK-based BCSK as the default modulation scheme. For the special case where $L=1$ (no-ISI) and $N=1$ (single sample per bit), the maximum likelihood (ML) detector for BCSK is of the threshold form and can be written as

\begin{equation}
\label{eq:threshold_FTD}
\hat{s}[k] = y[k] \mathop{\gtrless}_{0}^{1} \gamma.
\end{equation}

We refer to this specific scenario as the  fixed threshold detector (FTD) throughout the paper. For the no-ISI scenario, the optimal $\gamma$ value can be found by

\begin{equation}
\gamma = \frac{2M h[1]}{\ln \big( \frac{2M h[1] + \lambda_s}{\lambda_s} \big)}.
\end{equation}

The FTD is a very simple detector with very low computational complexity, which is a desirable trait for a nano- or micro-scale transceiver. However, it is also fragile against ISI (even when $\gamma$ is optimized considering ISI <cit.>) and is sub-optimal for the case where $L>1$. Motivated by this, several studies consider adaptively changing the threshold and perform

\begin{equation}
\label{eq:threshold_general}
\hat{s}[k] = y[k] \mathop{\gtrless}_{0}^{1} \gamma[k]
\end{equation}

to decode the symbol. We call this family of strategies adaptive threshold detectors (ATD) in this paper. The simplest form of ATDs selects $\gamma[k] = y[k-1]$ for $N=1$. For the general case where $N \geq 1$, it can be expressed as

\begin{equation}
\label{eq:ATD}
\hat{s}[k] = \sum_{n = (k-1)N +1}^{kN} y[k] \mathop{\gtrless}_{0}^{1} \gamma[k]
\end{equation}

where $\gamma[k] = \sum_{n = (k-2)N +1}^{(k-1)N} y[k]$, which is equivalent to the energy collected within the previous symbol duration <cit.>. This strategy is motivated by the fact that consequent symbols see comparable ISI energy, hence comparing the energy of consequent symbols provide information about the increase/decrease of concentration, hence the transmitted symbol. It is shown in <cit.> that this detector outperforms FTD in higher data rates whilst falling short when the communication rate is slower.

Note that the strategy of <cit.> does not take channel properties into account. Using the channel coefficient vector $\boldsymbol{h}$, <cit.> devises a more sophisticated ATD for BCSK,  denoted the memory-limited decision aided decoder (MLDA). Given the strong ISI, the MLDA also employs decision feedback, that is we seek the following detector structure:

\begin{equation}
\label{eq:MLDFE}
\hat{s}[k] = \argmax_{j \in \{0, 1\}} P(\boldsymbol{y}_k | s[k] = j, \hat{s}[k-L],\cdots, \hat{s}[k-1]),
\end{equation}

where $\boldsymbol{y}_k = \begin{bmatrix} y[(k-1)N +1] & \dots & y[kN] \end{bmatrix}^T$. In MLDA, the current symbol is detected under the assumption that the previous $L$ symbols were correctly detected. We observe that a maximum a posteriori sequence detector was also considered in <cit.>. For $N=1$, <cit.> shows that given perfect information about ISI and noise, the symbol-by-symbol maximum likelihood detector results in a threshold rule, where the threshold is given as:

\begin{equation}
\label{eq:MLDA_1}
\gamma[k] = \frac{2M h[1]}{\ln \big( \frac{2Mh[1] + \mathcal{I}[k] + N \lambda_s}{\mathcal{I}[k] + N \lambda_s} \big)}.
\end{equation}

Here, $\mathcal{I}[k]$ is the number of received molecules due to ISI. In practice, $\mathcal{I}[k]$ cannot be perfectly known, but can be estimated using past decisions and $\boldsymbol{h}$ by computing,

\begin{equation}
\hat{\mathcal{I}}[k] = \sum_{i=2}^{L-1} h[i] \hat{x}[k-i+1].
\end{equation}

We note that due to complexity reasons, a nano-transceiver might not hold all $L$ previously detected symbols in memory. In that case, <cit.> considers storing only $L' < L$ previously detected symbols. For the general case where $N \geq 1$ (i.e. multiple samples per symbol), MLDA's symbol-by-symbol maximum likelihood detection becomes

\begin{equation}
\label{MLDA_N}
\sum_{q=1}^{N} y[(k-1)N+q] \ln \left(\frac{2M h[q] + \hat{\mathcal{I}}_q[k] + \lambda_s}{\hat{\mathcal{I}}_q[k]  + \lambda_s}\right) \mathop{\gtrless}_{0}^{1}   2M \sum_{q=1}^{N} h[q],
\end{equation}

where $\hat{\mathcal{I}}_q[k]$ is the estimated mean ISI contribution on the $q^{th}$ sample of the intended symbol.

MLDA provides strong ISI mitigation, hence an improved error performance. However, as it uses previous decisions for ISI estimation, it might be prone to error propagation when channel conditions are not favorable. Furthermore, it requires a more complex receiver than the conventional FTD, which might be undesirable for some applications. Motivated by this, rather than adaptively changing $\gamma[k]$ and keeping the transmission power constant for a bit-$1$ (i.e. $2M$), <cit.> consider the FTD at the receiving end, but vary the number of emitted molecules depending on earlier transmissions (i.e. pre-equalization). By <cit.>, this strategy is called adaptive transmission rate and constant thresholding (ATRaCT). It should be noted that ATRaCT requires the availability of CSI (i.e., the vector $\boldsymbol{h}$) at the transmitter. That said, given CSI at the transmitter, ATRaCT is able to provide a strong performance improvement over FTD while still conserving the simplicity of the receiver nano-machine. Furthermore, it is shown in <cit.> that the scheme is more robust to distance variations than FTD, suggesting applicability when $d$ is not known exactly or is estimated erroneously.

Using the Viterbi decoder <cit.>, the MLSD is of complexity $\mathcal{O}(2^L S)$. This exponential complexity can become prohibitive when the data rate is increased, as one needs to consider a larger $L$ to capture the same portion of $f_{hit}(t)$. Motivated by this concern, decision-feedback (DFE) and linear minimum mean squared error (linear MMSE, LMSSE) equalizers are derived in <cit.>. We note that although they are classified under detectors, methods like MLDA also implicitly incorporate a DFE-type equalization within their set of operations. Using a similar approach, DFE is combined with Wald's sequential probability ratio test (SPRT, <cit.>) in <cit.> to devise a symbol-by-symbol detector. We will extend our discussion on SPRT approaches when we consider the impact of sychronization errors.
Overall, it can be inferred that despite its potential for error propagation, the DFE is a good and relatively low-complexity alternative to be paired with any type of DiMC detector.

Another new class of receiver side equalizers employ discrete time derivative pre-processing. First proposed by <cit.>; therein it is observed that
applying a time derivative to the received sample vector $\boldsymbol{y}$ is shown to decrease the time dispersion of the signal. These initial endeavors are generalized for an arbitrary derivative order $m$ in <cit.>, and it has been shown that higher order derivatives are able to outperform $m=1$ or no derivatives ($m=0$). However, it is also shown in <cit.> that increasing $m$ arbitrarily does not lead to a monotonically improving error performance, as $m$ governs a trade-off between ISI mitigation and noise amplification. Herein, we will characterize said noise amplification phenomenon.

$\frac{\partial^m f_{hit}(t)}{\partial t^m}$ for $m = 0$, $1$, and $2$ for $r_0 = \SI{15}{\micro\meter}$, $r_r = \SI{5}{\micro\meter}$, $D = 80 \frac{\SI{}{\micro\meter\squared}}{\SI{}{\second}}$. $m=0$ corresponds to the $f_{hit}(t)$ function itself in (<ref>).

Since the derivative operation is linear, it can be represented by a matrix $\boldsymbol{D}$, where

\begin{equation} \label{eq:Dmatrix}
\boldsymbol{D} = \begin{bmatrix}
-1 & 1 & 0  & \cdots & 0 & 0\\
0 & -1 & 1  & \cdots & 0 & 0\\
\vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\
\vdots & \vdots &  & -1 & 1 & 0 \\
0 & 0 & \cdots  & 0 & -1 & 1 \\
0 & 0 & \cdots  & 0 & 0 & -1 \\
\end{bmatrix}_{SN \times SN}.
\end{equation}

In addition, we can re-write the Gaussian approximation in (<ref>) in vector form as

\begin{equation} \label{eq:matrix_form}
\begin{split}
\boldsymbol{y} = (\boldsymbol{H}\boldsymbol{x} + \lambda_s \boldsymbol{j}) + \boldsymbol{\eta},
\end{split}
\end{equation}

where $\boldsymbol{j}$ is an vector of ones, $\boldsymbol{H}$ is a Toeplitz matrix of the coefficient vector $\boldsymbol{h}$, and $\boldsymbol{\eta} \sim \mathcal{N} (0,\textrm{diag}\{\boldsymbol{H}\boldsymbol{x}\} + \lambda_s \boldsymbol{I})$, following the definition in (<ref>). Denoting $\textrm{diag}\{\boldsymbol{H}\boldsymbol{x}\} + \lambda_s \boldsymbol{I}$ as $\boldsymbol{\Sigma}$, applying the $m^{th}$ order derivative can be represented by

\begin{equation} \label{eq:y_derivatived}
\begin{split}
\boldsymbol{y}_{(m)} &= \boldsymbol{D}^m \boldsymbol{y} \\
&= \boldsymbol{D}^m \boldsymbol{H}\boldsymbol{x} + \boldsymbol{D}^m \lambda_s \boldsymbol{j} + \boldsymbol{D}^m \boldsymbol{\eta}. \\
&= \boldsymbol{D}^m \boldsymbol{H}\boldsymbol{x} + \boldsymbol{D}^m \boldsymbol{\eta}.
\end{split}
\end{equation}

Note that since $\eta \sim \boldsymbol{\Sigma}$, we have $\boldsymbol{D}^m \boldsymbol{\eta} \sim \mathcal{N} (0,\boldsymbol{D}^m \boldsymbol{\Sigma} (\boldsymbol{D}^T)^m ) $, which results in noise amplification.

Overall, <cit.> show that the $\boldsymbol{D}^m$ operator provides powerful ISI mitigation due to compression of the right-tail of the $f_{hit}(t)$ function. As a consequence of this, if the transceiver is able to transmit using large transmission powers, it can circumvent the noise amplification issue and achieve an order of magnitude higher data rates than conventional schemes while preserving reliable communication. In addition, the $\boldsymbol{D}^m$ operator provides:

* Flexibility: Note that the $m^{th}$ order derivative is an equalizer-like operator, and similar to the DFE, can be used as a separate block before the detection step. This makes operators like $\boldsymbol{D}^m$ and DFE compatible with different kinds of detectors, and provide flexibility in design.

* A more feasible banded MLSD: Recall that the MLSD is of complexity $\mathcal{O}(2^L S)$ using the Viterbi decoder. One might consider a sub-optimal banded MLSD, by considering a shorter memory $L' < L$, <cit.>. Note that this version of the Viterbi decoder is of complexity $\mathcal{O}(2^{L'} S)$. However, in high data rate DiMC systems, this approach heavily underestimates ISI, resulting in poor error performance. That said, due to the aggressive ISI mitigation provided by the $\boldsymbol{D}^m$ operator, it is shown in <cit.> that such a banded MLSD yields operable error performances. Therefore, the $\boldsymbol{D}^m$ operator makes MLSD-like detectors considerably more feasible for high data rate DiMC systems.

* Very low computational complexity: The $m^{th}$ order derivative operator is computationally very cheap. Even though its description is made through matrix multiplications for clarity of argument, essentially each derivative operation requires a single discrete time shift and a vector subtraction operation. Hence, recalling $S$ as the block length and $N$ as the samples per bit parameter, the $m^{th}$ order derivative operator is of complexity $\mathcal{O}(mSN)$ per block, which is linear in block length. Note that unlike DFE where the complexity is tied to the considered memory window, this complexity is regardless of the channel memory $L$.

§.§ On CSI Estimation and Transceiver Design

Many detectors and equalizers mentioned in Section <ref> rely on channel state information (CSI) at the receiver end, which is generally considered to be either the $f_{hit}(t)$ function itself, or the channel coefficient vector $\boldsymbol{h}$. In many studies throughout the literature, receiver CSI is assumed to be readily available, whereas in reality, the $\boldsymbol{h}$ vector needs to be estimated using practical algorithms <cit.>. Herein, we discuss recent endeavors in channel estimation for DiMC systems and non-coherent transceiver design strategies.

DiMC channel estimation is studied in <cit.> in the context of joint parameter estimation, where the authors provide the Cramer-Rao Bound (CRB), ML estimation, and more practical estimators for jointly estimating $r_0$, $D$, molecule degradation rate, and number of emitted molecules for a single-shot emission. Note that estimating these parameters indirectly yields the $f_{hit}(t)$ function, if the form of $f_{hit}(t)$ is analytically known. However, the form is currently known for a small subset of topologies, under idealized conditions. Fortunately, for digital DiMC, the only information the transceivers need is the $\boldsymbol{h}$ vector to characterize the channel (assuming $\lambda_s = 0$ for the sake of discussion). Motivated by this, <cit.> proposes DiMC channel estimation strategies using pilot symbols. The CRB, ML estimation, and the linear least squares estimator are provided in the study, as well as strategies to select good pilot symbol sequences. In case the transceiver has access to prior information about $\boldsymbol{h}$, MAP estimation and corresponding Bayesian CRB are also derived. The endeavors of this study are generalized to molecular MIMO channels by <cit.>, where the estimation is done for channel coefficient matrix $\boldsymbol{H}_{\text{MIMO}}$. Note that the estimators of $\boldsymbol{H}_{\text{MIMO}}$ need to also take ILI.

As they rely on introducing pilot symbols in front of a data block, the strategies presented in <cit.> inherently reduce the effective data rate. Motivated by decreasing the number of pilot bits, the study is extended by <cit.> where semi-blind channel estimation strategies are considered. In this context, expectation-maximization and decision-directed methods are proposed.

Instead of estimating $\boldsymbol{h}$ and employing the aforementioned detectors/equalizers in Subsection <ref>, one may also consider non-coherent detectors that circumvent the need of channel estimation. For instance, some modulation schemes are compatible with simple $\argmax$ demodulators, i.e. maximum count detectors (MCD). Examples include the basic type modulation (MoSK), pulse-position modulation (PPM), and the basic spatial modulation (MSSK <cit.>). We note that since threshold optimization strategies require the CSI at the receiver end, the conventional demodulator of BCSK (FTD presented in (<ref>)) cannot be realized non-coherently. However, the simple ATD proposed in <cit.> is capable of non-coherent operation, and might be compatible for simple nano- or micro-scale use due to its low complexity. In addition, neural network-based non-coherent detection is proposed in <cit.> for DiMC systems. It is shown in <cit.>, that the proposed detectors outperform Viterbi-based MLSD detection under imperfect CSI. That said, the method of<cit.> involves training a neural network, which can be a substantial effort, we note that their implementation is more suitable for macro-scale testbeds such as <cit.>.

§.§ Synchronization

Like the availability of CSI, perfect synchronization between the transmitter and the receiver is also generally assumed for DiMC system design. However, achieving synchronization is a particularly challenging task in a DiMC system, especially considering the low molecule budget and limited computation capability of molecular transceivers. To this end, we discuss DiMC transceiver clock synchronization strategies and asynchronous receiver design herein.

In <cit.>, clock synchronization for a one-way DiMC with positive flow is considered. In the presented scheme, the messenger molecule is assumed to hold information about the transmission instant, which the receiver exploits to estimate the delay. In addition, a two-way message exchange is considered in <cit.>. In <cit.>, delay and clock frequency imperfections are jointly considered. When the information molecule holds the transmission instant, it is shown by <cit.> that clock recovery can be performed using a single symbol transmission in a single-input multiple-output (SIMO) system, using the arrival time differences at different antennas.

The above methods are inspired by conventional radio frequency-based (RF) communications, where the timing information is typically embedded within a packet. That said, encoding the transmission instant information within the molecule's chemical structure is a complex task for a nano-machine to handle. Therefore, it is of particular interest to discuss methods that do not assume such an encoding. In <cit.>, the DiMC synchronization problem is tackled using two types of molecules: one for synchronization and the other for information transmission. For $K$-ary MoSK modulation, a blind synchronization approach (where the data is treated as a random variable) is presented in <cit.>. A dual-molecule clock synchronization scheme for mobile molecular communication is also presented in <cit.>.

Even without considering timing information within messenger molecule structure, the aforementioned synchronization methods rely on the use of multiple molecule types and or complex biochemical processes, which might be undesirable in terms of nano- to micro-scale transceiver complexity. Motivated by this, <cit.> considers synchronization for single-molecule type intensity modulations.

In addition to delay estimation approaches, the synchronization problem can also be tackled through developing asynchronous detection schemes. The simple asynchronous detector (ADS) and the asynchronous detector with decision feedback (ADDF) are presented in <cit.>. For each BCSK symbol, ADS picks the sample with largest molecule count and compares the count with a fixed threshold, which can be written as

\begin{equation}
\label{eq:ADS}
\hat{s}[k] = \mathop{\max}_{q \in \{1, \dots, N \}} y[(k-1)N+q] \mathop{\gtrless}_{0}^{1} \gamma.
\end{equation}

Note that ADS resembles the conventional FTD presented in (<ref>), with the main difference being in the considered received signal (i.e. ADS considers the maximum sample, whereas (<ref>) considers the sum of samples). Using the same principle with an MLDA-like DFE mechanism, ADDF considers to mitigate ISI before employing ADS, which can be described by

\begin{equation}
\label{eq:ADDF}
\hat{s}[k] = \mathop{\max}_{q \in \{1, \dots, N \}} \big\{y[(k-1)N+q] - \hat{\mathcal{I}}_q[k] \big\} \mathop{\gtrless}_{0}^{1} \gamma,
\end{equation}

where the estimated ISI mean $\hat{\mathcal{I}}_q[k]$ can be found in a similar manner to MLDA.

We note that the maximum operations in (<ref>) and (<ref>) do not predict which sample holds the maximum according to the available CSI: they just take the maximum sample assuming the delay is $\tau = 0$. This way, even though the fixed-maximum sample method misses the peak for $\tau \neq 0$, ADS and ADDF are able to recover the maximum concentration as long as the peak stays within the symbol duration.  In <cit.>, both ADS and ADDF are shown to outperform conventional schemes in terms of robustness against synchronization errors.

Another approach by which to increase robustness to synchronization errors is to adopt the principles of a sequential probability ratio test <cit.>. In this scenario, the decision rule $\delta(\cdot)$ can be expressed as

\begin{align}
\label{eq:DFE_SPRT}
\delta(L_q(\mathbf{y}_q^k| \mathbf{x}))=\begin{cases}
\hat{s}[k] = 0, &L_q(\mathbf{y}_q^k | \mathbf{x}) \leq A_{\text{SPRT}}\\
\hat{s}[k] = 1, &L_q(\mathbf{y}_q^k | \mathbf{x}) \geq B_{\text{SPRT}}\\
\mbox{sample}, &\mbox{else.}
\end{cases}
\end{align}

where $q < N$, $L_q(\cdot)$ is the computed likelihood from the collected $\mathbf{y}_q^k$ vector, and $\mathbf{y}_q^k = \begin{bmatrix} y[(k-1)N+1] & \dots & y[(k-1)N+q] \end{bmatrix}^T$.
The constants $A_{\text{SPRT}}$ and $B_{\text{SPRT}}$ are given by (as in <cit.>)

\begin{equation}
A_{\text{SPRT}}=\frac{1-P_D}{1-P_{FA}}, \; B_{\text{SPRT}} = \frac{P_D}{P_{FA}}
\label{eq:sprt_thresholds}
\end{equation}

where $P_{FA}$ is the false alarm rate and $P_D$ is the detection rate. The values of $P_{FA}$ and $P_{D}$ are set manually as targets. If a decision has not been made by $m=N$, a truncation rule based on the minimum distance rule is applied to prevent sampling into the subsequent channel symbol. In <cit.>, the DFE-SPRT is proposed wherein the likelihood function is generalized to allow for decision feedback as in (<ref>).

In <cit.>, the DFE-SPRT is shown to provide strong robustness against synchronization errors, at the cost of a moderate increase in computational complexity over ADDF and MLDA. Furthermore, the synchronization error is estimated using a single molecule type by the maximum log-likelihood estimator. The CRB of the estimation is derived and the estimator is shown to yield close-to-CRB performance when the data rate is low. As described in (<ref>), the DFE-SPRT operates in iterations, where at each iteration, the receiver takes an additional sample and computes a likelihood ratio test (LRT) with the available vector of samples. Depending on the confidence, the symbol might be detected, or an additional sample is collected. We note that since the DiMC system is subject to ISI and the synchronization error estimation is non-perfect, the LRT is mismatched. It is shown in <cit.> that under these imperfections, there exists an optimal number of samples and this number decreases with increasing mismatch.

An interesting phenomenon is that for a DiMC system using BCSK and FTD, a small clock mismatch between the transmitter and the receiver actually improves error performance <cit.>. It is noted in <cit.> that the molecules that arrive within the first $\tau$ duration of the symbol interval are mostly due to ISI. Hence, assuming the receiver clock lags the transmitter clock by $\tau$, such a delay acts as an ISI mitigator and improves error performance. Of course, a large mis-synchronization results in bad performance, since the intended received signal power also decreases. It is shown in <cit.> that this effect is applicable to PPM and MCPM as well. Furthermore, it is demonstrated in <cit.> that higher order modulation schemes yield increased robustness against more severe synchronization errors (i.e. large $\tau$). We note that while satisfying the same bit duration ($t_b$) constraint, a higher order modulation scheme is able to transmit at a larger symbol duration (see Section <ref>). This makes $\tau$ comparatively smaller with respect to $t_{sym}$, alleviating the debilitating effect of the mismatch.

Overall, the synchronization problem is a particularly challenging one for nano- to micro-scale DiMC systems, due to their low-complexity requirements and the DiMC channels' properties. Although the sub-field is still relatively less explored, we believe that detector and equalizer designs that seek to provide robustness against synchronization errors may have greater practical utility than those that do not.  This remains an open avenue for research.

In addition, single-molecule type clock recovery strategies can be devised, as they are in line with the low-complexity requirements of DiMC transceivers.