# Large shift current via in-gap and charge-neutral exciton excitations in BN
nanotubes and single BN layer

Yi-Shiuan Huang1 Yang-Hao Chan2,3<EMAIL_ADDRESS>Guang-Yu Guo1,2
<EMAIL_ADDRESS>1Department of Physics and Center for Theoretical
Physics, National Taiwan University, Taipei 10617, Taiwan  2Physics Division,
National Center for Theoretical Sciences, Taipei 10617, Taiwan  3Institute of
Atomic and Molecular Sciences, Academia Sinica, Taipei 10617, Taiwan

###### Abstract

We perform ab initio many-body calculations to investigate the exciton shift
current in small diameter zigzag BN nanotubes and also single BN sheet, using
the GW plus Bethe-Salpeter equation (GW-BSE) method with the newly developed
efficient algorithms. Our GW-BSE calculations reveal a giant in-gap peak in
the shift current spectrum in all the studied BN systems due to the excitation
of the A exciton. The peak value of the excitonic shift current is more than
three times larger than that of the quasiparticle shift current, and is
attributed to the gigantic enhancement of the optical dipole matrix element by
the A exciton resonance. The effective exciton shift current conductivity is
nearly ten times larger than the largest shift conductivity observed in
ferroelectric semiconductors. Importantly, the direction of the shift current
in the BN nanotubes is found to be independent of the tube chirality ($n,0$)
(or diameter), contrary to the simple rule of
$sgn(J_{\text{shift}})=\text{mod}(n,3)$ predicted by previous model
Hamiltonian studies. Finally, our ab initio calculations also show that the
exciton excitation energies decrease significantly with the decreasing
diameter due to the curvature-induced orbital rehybridization in small
diameter zigzag BN nanotubes.

##  I INTRODUCTION

Shift current is one of the primary mechanisms for the bulk photovoltaic
effect (BPVE) (also known as the photogalvanic effect), which can generate DC
photocurrent in noncentrosymmetric crystals due to their second order optical
response to light irradiation. The shift of the real space Wannier charge
center of the excited electron is responsible for the shift current [1, 2, 3,
4]. In contrast to the conventional photovoltaic effect, shift current is a
bulk phenomenon that does not require a p-n junction to separate the optically
generated electron-hole pairs for a DC photocurrent. Consequently, the BPVE
can be exploited to generate the above band gap photovoltage [5] and thus to
fabricate high power-conversion efficiency solar cells [6, 7]. Therefore,
there has been a resurge of interest in the BPVE in recent years.

An exciton is a type of collective excitation formed by a bound electron-hole
pair interacting via Coulomb interaction, and can be optically generated.
Although an exciton itself is charge-neutral, the process of exciton
excitation can generate a shift current called exciton shift current due to
Wannier charge center shift of the electron and hole in noncentrosymmetric
crystals [8]. Furthermore, the exciton shift current can have exotic sub-
bandgap peaks owing to the exciton binding energy [8, 9]. Recently, the sub-
bandgap exciton shift current was observed in the noncentrosymmetric
semiconductor CdS by using THz emission spectroscopy at the low temperature of
2 K[10]. The amplitude of the exciton shift current is comparable to the shift
current driven by free electron-hole pair excitation.

BN nanotubes (BN-NTs) are formed by rolling up a single hexagonal BN sheet
(see Fig. 1) along a specific chiral vector $(n,m)$. [11] There are three
types of BN-NTs, namely armchair $(n,n)$ nanotubes, zigzag $(n,0)$ nanotubes
(Fig. 1), and chiral $(n,m)$ nanotubes where $n\neq m$ [11]. However, there is
no second-order nonlinear optical (NLO) response (e.g., second-harmonic
generation and BPVE) in armchair BN-NTs [12]. Moreover, previous experiments
[13] indicate that among the grown BN-NTs, the zigzag structure is usually
favored. Therefore, we focus on zigzag BN-NTs in this paper.

BN-NTs exhibit large many-body interaction effects owing to their one-
dimension (1D) and wide bandgap nature, as have been demonstrated by previous
ab initio many-body theory studies at the level of GW plus Bethe-Salpeter
equation (GW-BSE) [14, 15]. Furthermore, huge diameter-dependent A exciton
peaks have been observed in BN-NTs [17, 16], being consistent with the
theoretical predictions [14, 15]. Therefore, we expect zigzag BN-NTs to be
ideal candidates for observing the large subbandgap exciton shift current due
to their following two properties [10]. First, exciton absorption spectra in
BN-NTs are well separated from the continuum of free electron-hole excitation
owing to their significantly renormalized optical spectra and large exciton
binding energy [14, 15]. Second, strongly bounded excitons with a large
binding energy impede thermal dissociation into free electron-hole pairs [14,
15].

Figure 1: Crystal structures of the zigzag BN-NT (6,0), and the hexagonal BN
sheet. (a) Side and (b) top views of the (6,0) BN-NT. (c) Top view of the
single BN sheet. The solid box in (a) indicates the unit cell. Dashed lines in
(a) and (c) denote the mirror planes.

However, ab initio studies of the BPVE in BN-NTs at GW-BSE level have not been
reported, mainly because such ab initio calculations are technically and
computationally challenging, as will be explained below. Nevertheless, simple
model Hamiltonian calculations without [18] and with [19] the inclusion of the
the excitonic effect have been performed, predicting that zigzag BN-NTs have
shift current along the tube axis with the direction determined by the chiral
index ($n$,0). In particular, the shift current direction would follow the
simple rule of sgn($J_{\text{sh}})=\text{mod}$($n$,3) [18, 19]. According to
this rule, in a bundle of zigzag BN-NTs, 1/3 of them would have the shift
current along the positive tube axis, another 1/3 would possess the current
flowing in the opposite direction and the other 1/3 would have the zero
current. That is to say, a bundle of zigzag BN-NTs would have either zero or
very small net shift current, i.e., zigzag BN-NTs would not be suitable for
applications in photovoltaic solar cells and nonlinear optoelectronic devices.
On the other hand, the prediction of this simple rule is theoretically
surprising. A symmetry analysis (see Sec. II. A below) indicates that the
shift current would flow along the $y$-axis in a single BN sheet [see Fig.
1(c)], which is also the tube axis when the BN sheet is rolled up to form a
zigzag BN-NT. Consequently, for zigzag BN-NTs with a large diameter, the shift
current would always be along the positive tube axis since the curvature
effect would be very small [20, 12], i.e., the above mentioned rule would not
occur. In this context, it is important to perform the state-of-the-art ab
initio calculations to investigate this detrimental prediction.

In this work, therefore, we perform the state-of-the-art ab initio GW-BSE
calculations of the shift current in small diameter zigzag BN-NTs
[$(5,0),(6,0),(7,0),(8,0)$] as well as single BN layer. Because of possible
large curvature effect on the optical properties [20, 12], the shift current
in small BN-NTs could depend significantly on the tube diameter. On the other
hand, the optical properties of large diameter BN-NTs would be rather similar
to that of the single BN sheet [20, 12], and thus are not considered here.
Among other things, our ab initio calculations reveal that the shift current
due to the exciton excitations will be dramatically enhanced. Furthermore, the
direction of the shift current calculated without and with the excitonic
effect included, is always along the $c$-axis, i.e., being independent on the
tube index ($n,0$) (or tube diameter). Therefore, our work demonstrates that
the zigzag BN-NT bundles are promising materials for high power conversion
efficiency solar cells as well as high sensitivity photodectors.

We would like to comment that extending the existing ab initio GW-BSE approach
to calculate the exciton shift current is not straightforward [9]. In fact,
there has been only one reported fully ab initio study on the exciton shift
current [9]. In Ref. [21], to account for the excitonic effect, Fei et al.
used the linear optical coefficients derived from their ab initio GW-BSE
calculations. However, since the Coulomb interaction between the electron and
hole in an exciton is not explicitly taken into account, such calculations are
still based on an independent-particle approximation [21] and hence are not a
fully ab initio GW-BSE approach. In Ref. [9], to take the strong excitonic
effect in two-dimensional (2D) materials into account, a time-dependent
adiabatic GW (TD-aGW) approach was developed and used to calculate the exciton
shift current. Although the ab initio TD-aGW method can properly include the
excitonic effect on the optical responses [9], the extremely high
computational cost of conducting real-time propagation prevents it from
studying complex structures such as BN-NTs. In this work, we thus develop a
computationally efficient approach that combines the GW-BSE and sum-over-state
formalism derived from the perturbative density-matrix approach within the
mean-field approximation [22] to calculate the exciton shift current.

The rest of this paper is organized as follows. In Sec. II, we introduce the
crystal structure of the BN-NTs and the BN sheet, followed by a brief
description of the theories and computational details used in this work. In
particular, the computationally efficient approach mentioed above for the
exciton shift current calculations will be outlined. The main results are
presented in Sec. III. In Sec. III A, we present the calculated electronic and
optical properties of the BN sheet. The distinguished features in the optical
spectra are analysed in terms of the electronic band structure and interband
optical transition matrix elements. In Sec. III B, the electronic properties
of the BN-NTs are reported, which will be used to understand the calculated
optical absorption and shift current spectra in subsections that follows. In
Secs. III C and III D, the calculated optical absorption and shift current
spectra of the BN-NTs are presented, respectively. Finally, the conclusions
drawn from this work are summarized in Sec. IV.

##  II CRYSTAL STRUCTURES AND COMPUTATIONAL METHODS

###  II.1 Symmetry and shift conductivity tensor

Single hexagonal BN sheet is a noncentrosymmetric crystal with space group
$P\overline{6}m2$ and point group $D_{3h}$. A symmetry analysis would show
that the single BN sheet has three nonzero shift conductivity tensor elements
$\sigma^{xxy}$, $\sigma^{yxx}$, and $\sigma^{yyy}$, where the first Cartesian
index is the direction of the current and the second and third indice are the
polarization directions of the external field. Furthermore,
$\sigma^{xxy}=\sigma^{yxx}=-\sigma^{yyy}$ [12, 23], i.e., there is only one
inequivalent nonzero element. As an example, let us consider the difference
between in-plane shift conductivity tensor elements $\sigma^{xxx}$ and
$\sigma^{yyy}$ in the single BN sheet [see Fig. 1(c)]. Figure 1(c) indicates
that the atomic structures on the left and right regions to the mirror plane
(denoted by the dashed line) (the $y$-axis) are symmetric. Therefore, the
shift current along the direction normal to the mirror plane (the $x$-axis)
would be zero if the second order combination of external fields is even under
the mirror symmetry, i.e., $\sigma^{xxx}=0$. On the other hand, such mirror
symmetry does not exist in the direction along the mirror plane (the
$y$-axis), and thus the conductivity element $\sigma^{yyy}$ would not be zero.

When a hexagonal BN sheet is wrapped up to form a zigzag BN-NT with chiral
index ($n,0$), the $y$-axis would become the tube axis [the $c$-axis, see Fig.
1(a)] and the point group of the resultant BN-NT ($n$,0) is $C_{2nc}$. Thus,
the conductivity tensor element $\sigma^{ccc}$ is nonzero. On the other hand,
there is no azimuthal shift current in the zigzag BN-NTs because it
corresponds to the direction of the $x$-axis in the single BN sheet (see Fig.
1).

###  II.2 Density functional theory calculations

Ab initio calculations, based the density functional theory (DFT) with the
local density approximation (LDA), are performed to determine the ground state
properties of the considered zigzag BN-NTs and also the single BN sheet. A
supercell geometry is adopted to simulate a BN-NT in which the nanotubes are
arranged in a square array with a minimum distance of 12 Å between the
neighboring nanotubes. A slab supercell method is used to model the single BN
sheet, and the inter-sheet distance used is over 16 Å. In the structural
optimization calculations, the accurate projector-augmented wave (PAW) method
plus the conjugate gradient approach, as implemented in the VASP package [24,
25], is used to determine the atomic positions and lattice constants of BN-
NTs. A large plane-wave cutoff energy of 450 eV is adopted. Theoretical
equilibrium structures are obtained when the forces acting on all the atoms
and the uniaxial stress were less than 0.005 eV/Å and 2.0 kBar, respectively.

The ground electronic structure calculations are performed by using the plane-
wave pseudopotential method as implemented in the Quantum Espresso package
[26]. The optimized norm-conserving Vanderbilt pseudopotential [27] is
exploited here. The $1\times 1\times 32$ and $18\times 18\times 1$ Monkhorst-
Pack $k$-grids [28] are used to evaluate the Brillouin zone (BZ) integrals for
the BN-NTs and the single BN sheet, respectively. The energy cutoff for the
plane-wave basis set is 50 Ry. The resultant electronic structures are used in
the subsequent GW-BSE calculations, as described below.

###  II.3 Quasiparticle band structure calculations

The present GW-BSE calculations are performed via the BerkeleyGW package [29,
31, 30]. The quasiparticle energy bands are calculated by solving the Dyson
equation,

$[-\frac{1}{2}\nabla^{2}+V_{\mathrm{ion}}+V_{\mathrm{H}}+\Sigma(E^{\mathrm{QP}}_{n\bf{k}})]\psi^{\mathrm{QP}}_{n\bf{k}}=E^{\mathrm{QP}}_{n\bf{k}}\psi^{\mathrm{QP}}_{n\bf{k}},$
(1)

where $\Sigma$, $E^{\mathrm{QP}}_{n\bf{k}}$, and
$\psi^{\mathrm{QP}}_{n\bf{k}}$ are the self-energy operator, the energy, and
the wave function of the quasiparticles within the $G_{0}W_{0}$ approximation,
respectively [31].

In the present one shot $G_{0}W_{0}$ calculations [32], a nonuniform neck
subsampling (NNS) $k$-grid of $1\times 1\times 8$ ($18\times 18\times 1$) with
a subsampling of 10 points in the mini-Brillouin zone, 600 (1600) bands, and a
dielectric cutoff energy of 50 (50) Ry for the BN-NTs (the single BN sheet)
are used. Truncation of the Coulomb interactions between the BN-NT (the BN
sheet) and its periodic images is implemented [33]. The dynamic dielectric
matrix is computed within the independent particle approximation (IPA) and
Hybertsen-Louie generalized plasmon pole model [29].

###  II.4 Exciton excitation calculations

The exciton wavefunction can be expressed as the linear combination of the
quasiparticle electron-hole pairs,

$\Psi_{s}(\textbf{r}_{e},\textbf{r}_{h})=\sum\limits_{\textbf{k},m,n}A^{s}_{mn\textbf{k}}\psi_{\textbf{k},n}(\textbf{r}_{e})\psi^{*}_{\textbf{k},m}(\textbf{r}_{h}),$
(2)

where band index $m$ ($n$) sums over the valence (conduction) bands only.
Exciton envelope function of the $s$-th exciton state $A^{s}_{mn\textbf{k}}$
can be obtained by solving the BSE,

$(E^{\mathrm{QP}}_{n\textbf{k}}-E^{\mathrm{QP}}_{m\textbf{k}})A^{s}_{mn\textbf{k}}+\sum\limits_{m^{\prime}n^{\prime}\textbf{k'}}\braket{mn\textbf{k}}{K^{eh}}{m^{\prime}n^{\prime}\textbf{k'}}=\Omega^{s}A^{s}_{mn\textbf{k}}$
(3)

where $\Omega^{s}$ and $E^{\mathrm{QP}}_{m\textbf{k}}$
($E^{\mathrm{QP}}_{n\textbf{k}}$) are the $s$-th exciton excitation energy and
the valence (conduction) band quasiparticle excitation energies, respectively.
$K^{eh}$ is the electron-hole interaction kernel, which includes an exchange
repulsive bare Coulomb term and a direct electron-hole attractive screened
Coulomb term [31, 30]. The imaginary part of the dielectric function including
the excitonic effect can be expressed as

$\displaystyle\varepsilon^{"}_{\mathrm{BSE}}$ $\displaystyle=\frac{g_{s}\pi
e^{2}}{\epsilon_{0}N_{k}V_{c}}\sum\limits_{s}|\sum\limits_{mn\textbf{k}}A^{s}_{mn\textbf{k}}\textbf{x}_{mn\textbf{k}}\cdot\textbf{e}|^{2}\delta(\omega-E^{s})$
$\displaystyle=\frac{g_{s}\pi
e^{2}}{\epsilon_{0}N_{k}V_{c}}\sum\limits_{s}|\textbf{X}_{s}\cdot\textbf{e}|^{2}\delta(\omega-E^{s})$
(4)

where $\textbf{x}_{mn\textbf{k}}\cdot\textbf{e}$ is the dipole matrix element
along the polarization direction e and $g_{s}$ accounts for the spin
degeneracy. $N_{k}$ is the number of $k$ points and $V_{c}$ is the unit cell
volume. The exciton dipole matrix element $\textbf{X}_{s}$ is defined as

$\textbf{X}_{s}\equiv\sum\limits_{mn\textbf{k}}A^{s}_{nm\textbf{k}}\textbf{x}_{mn\textbf{k}}.$
(5)

When neglecting the excitonic effect, optical excitations are given by direct
transitions between quasiparticle electron-hole pairs. Within the IPA
approximation, the imaginary part of the dielectric function reduces to [30],

$\varepsilon^{"}_{\mathrm{IPA}}=\frac{g_{s}\pi
e^{2}}{\epsilon_{0}N_{k}V_{c}}\sum\limits_{mn\textbf{k}}|\textbf{x}_{mn\textbf{k}}\cdot\textbf{e}|^{2}\delta(\omega-(E^{\mathrm{QP}}_{n\textbf{k}}-E^{\mathrm{QP}}_{m\textbf{k}})).$
(6)

Equation (4) has the same structure with Eq. (6) except that the dipole matrix
element is replaced by the exciton dipole matrix element. Our BSE calculations
for the BN-NTs (the single BN sheet) are computed on a dense $k$-grid of
$1\times 1\times 64$ ($72\times 72\times 1$) with the dielectric cutoff of 10
(10) Ry.

###  II.5 Exciton shift current calculations

The shift current density along the $a$-axis is given by [4]

$J^{a}_{sh}(\omega)=2\sum\limits_{bc}\sigma^{abc}(0;\omega,-\omega)E^{b}(\omega)E^{c}(-\omega),$
(7)

where $\sigma^{abc}$ is the third-rank conductivity tensor, $b$ and $c$ denote
the polarization directions of the electric fields of the incident light. In
the IPA, $\sigma^{abc}$ can be expressed as [4, 34]

$\sigma^{abc}(0;\omega,-\omega)=-\frac{i\pi
e^{3}}{4N_{k}V_{c}}\sum\limits_{n,m,\textbf{k}}f_{nm}\\\
(\textbf{x}^{b}_{mn}[\textbf{x}^{c}_{nm}]_{;k_{a}}+\textbf{x}^{c}_{mn}[\textbf{x}^{b}_{nm}]_{;k_{a}})[\delta(\omega_{mn}-\omega)+\delta(\omega_{nm}-\omega)],$
(8)

where $\textbf{x}^{a}_{nm}$ and $[\textbf{x}^{a}_{nm}]_{;\textbf{k}}$ are the
dipole matrix element and its generalized momentum derivative, respectively.
The generalized momentum derivative of quantity $O_{mn}(\textbf{k})$ is given
by

$[O^{b}_{nm}(\textbf{k})]_{;k_{a}}=\partial_{ka}O^{b}_{nm}(\textbf{k})-i[\xi^{a}_{nn}(\textbf{k})-\xi^{a}_{mm}(\textbf{k})]O^{b}_{nm}(\textbf{k})$
(9)

where $\xi^{a}_{nm}(\textbf{k})$ is the $a$ component of the Berry connection.
[4] $f_{nm}=f_{n}-f_{m}$ and $\hbar\omega_{nm}=E_{n}-E_{m}$ where $f_{n}$
denotes the occupation factor and $E_{n}$ is the energy of the $n$th band at
the k point. Note that here band indices $m$ and $n$ should sum over all the
states.

As mentioned before, to include the excitonic effect, we use the efficient
perturbative density-matrix approach to calculate $\sigma^{abc}$ [see Eq.
(B1a) in Ref. [22]]. In this case, $\sigma^{abc}$ can be written as

$\displaystyle\sigma^{abc}(0,\omega,-\omega)$
$\displaystyle=-\frac{g_{s}e^{3}}{2V}\sum_{ss^{\prime}}\left[\frac{X_{s}^{b}\Pi_{ss^{\prime}}^{a}X_{s^{\prime}}^{c*}}{\left(\hbar\omega-
E_{s}+i\eta\right)\left(\hbar\omega-
E_{s^{\prime}}-i\eta\right)}+\frac{V_{s}^{a*}X_{ss^{\prime}}^{b}X_{s^{\prime}}^{c}}{-E_{s}\left(\hbar\omega-
E_{s^{\prime}}+i\eta\right)}+\frac{V_{s}^{a}X_{ss^{\prime}}^{b*}X_{s^{\prime}}^{c*}}{-E_{s}\left(-\hbar\omega-
E_{s^{\prime}}-i\eta\right)}\right]$ $\displaystyle+\left(b\leftrightarrow
c,\omega\leftrightarrow-\omega\right),$ (10)

where we define
$V^{a}_{s}\equiv\sum_{cv\mathbf{k}}A^{s*}_{cv\mathbf{k}}v^{a}_{cv\mathbf{k}}$
with velocity matrix elements $v^{a}_{cv\mathbf{k}}$, and inter-exciton
coupling matrix elements, $X_{ss^{\prime}}$ and $\Pi_{ss^{\prime}}$ are
defined as,

$\Pi_{ss^{\prime}}^{a}\equiv\sum_{cvc^{\prime}\mathbf{k}}A_{cv\mathbf{k}}^{s}v_{c^{\prime}c\mathbf{k}}^{a}A_{c^{\prime}v\mathbf{k}}^{s^{\prime}*}-\sum_{cvv^{\prime}\mathbf{k}}A_{cv^{\prime}\mathbf{k}}^{s}v_{v^{\prime}v\mathbf{k}}^{a}A_{cv\mathbf{k}}^{s^{\prime}*},$

and $X_{ss^{\prime}}^{b}=Y_{ss^{\prime}}^{b}+Q_{ss^{\prime}}^{b}$. We note
that $c$ and $v$ denote that the summation runs over the conduction and
valence band, respectively. Here, $\textbf{X}_{ss^{\prime}}$ is an ill-defined
operator. However, it can be separated into the well-defined
$\textbf{Y}_{ss^{\prime}}$ (interband part) and $\textbf{Q}_{ss^{\prime}}$
(intraband part) operators, where

$Q_{ss^{\prime}}^{b}=i\sum_{c^{\prime}v^{\prime}\mathbf{k^{\prime}}}A_{c^{\prime}v^{\prime}\mathbf{k}^{\prime}}^{s*}\left(A_{c^{\prime}v^{\prime}\mathbf{k}^{\prime}}^{s^{\prime}}\right)_{;k_{b}^{\prime}},$

and

$Y_{ss^{\prime}}^{b}=\sum_{\begin{subarray}{c}c^{\prime}\neq c_{1}\\\
v^{\prime}\neq
v_{1}\end{subarray}}\sum_{c_{1}v_{1}\mathbf{k}}A_{c^{\prime}v^{\prime}\mathbf{k}^{\prime}}^{s*}\left[A_{c_{1}v^{\prime}\mathbf{k}^{\prime}}^{s^{\prime}}x_{c^{\prime}c_{1}}^{b}-A_{c^{\prime}v_{1}\mathbf{k}^{\prime}}^{s^{\prime}}x_{v_{1}v^{\prime}}^{b}\right].$

We develop a post-processing program to implement this formalism, which is
then used to calculate the excitonic shift current conductivity using the
outputs from the BerkeleyGW package.

It is not straightforward to evaluate the intraband portion of the position
operator $\textbf{Q}_{ss^{\prime}}$. $\textbf{Q}_{ss^{\prime}}$ involves the
numerical derivative of the exciton envelope function $A^{s}_{nm\bf{k}}$ with
respect to k. Furthermore, $A^{s}_{mn\bf{k}}$ has an arbitrary $k$-dependent
gauge. To solve the problem, we use the locally smooth gauge adopted in Ref.
[9] to calculate the intraband position operator. The idea of the locally
smooth gauge is to rotate the wave functions at neighboring $k$ points in such
a way that the overlap of connected wave functions is Hermitian [9, 35, 36].

In this work, the $\delta$-function is approximated by a Gaussian function
with a 0.1 eV broadening. For the BN NTs, we use the effective unit cell
volume rather than the supercell volume. The effective unit cell volume is
given by $V_{c}=\pi[(D/2+d/2)^{2}-(D/2-d/2)^{2}]T=\pi DdT$, where $D$ is the
tube diameter, $T$ is the length of translational vector and $d$ is the
effective thickness of the nanotube walls, which is set to the interlayer
distance (3.28 Å) of $h$-BN [20]. For the single BN sheet, the effective unit
cell volume is $V_{c}=A_{c}d$ where $A_{c}$ and $d$ are the area of the unit
cell and effective thickness of the sheet, respectively.

Figure 2: Quasiparticle band structure of the single BN sheet from our GW
(solid red lines) and LDA (blue dashed lines) calculations. Valence band
maximum is set at 0 eV, which is at the $K$ point indicated by the pink
diamond. The conduction band minima of the GW and LDA band structures are
denoted by the green star and orange square, respectively.

##  III RESULTS AND DISCUSSION

###  III.1 Electronic and optical properties of the single BN sheet

Both the GW and LDA quasiparticle band structures of the single BN sheet are
displayed in Fig. 2. Our GW quasiparticle band structure agrees well with the
previous calculations [37, 38, 39, 43, 40, 41, 42]. The single BN sheet shares
the similar structure with graphene, the honeycomb lattice. However, unlike
graphene which is a semimetal with $\pi$ and $\pi^{*}$ bands degenerate at the
$K$ point, here the $\pi$ and $\pi^{*}$ bands are well separated due to the
ionicity difference between B and N [33]. Hence, the single BN sheet has a
large band gap (Fig. 2). Table I shows that the single BN sheet is an
insulator with a $K^{v}\rightarrow K^{c}$ direct band gap of 4.65 eV from the
LDA calculation. However, the single BN sheet becomes an insulator with a
larger $K^{v}\rightarrow\Gamma^{c}$ indirect band gap of 7.66 eV from the GW
calculation. The large GW quasiparticle correction to the LDA band gap can be
attributed to the weak dielectric screening, a consequence of 2D nature and
the wide band gap of the BN sheet. The GW quasiparticle correction of the
$K^{v}\rightarrow K^{c}$ gap (3.26 eV) is slightly (0.37 eV) larger than that
of $K^{v}\rightarrow\Gamma^{c}$ (2.89 eV). This makes the single BN sheet
transform from the direct band gap ($K^{v}\rightarrow K^{c}$) in the LDA to
the indirect bandgap ($K^{v}\rightarrow\Gamma^{c}$) at the GW level. The
$\pi^{*}$ band, arising from $2p_{z}$ orbitals, has the out-of-plane charge
density and hence the weaker dielectric screening. The weaker dielectric
screening of the $\pi^{*}$ band explains why the quasiparticle correction of
$K^{v}\rightarrow K^{c}$ gap is smaller than that of
$K^{v}\rightarrow\Gamma^{c}$ gap [44]. Consequently, the GW quasiparticle
correction not only modifies the band gap, but also changes the dispersion of
the bands. This indicates the significance of GW quasiparticle calculations,
and thus the complex screening effect cannot be fully taken into account
through a simple scissor correction. We notice that our GW band widths are in
nearly perfect agreement with the experimental values (Table I).

Table 1: Transition energies from the valence band maximum ($v$) to conduction
band minimum ($c$) ($\Delta E$), band gaps ($E_{g}$), band widths, the onset
of the continuum optical absorption ($E_{abs}$), exciton excitation energy
(i.e., optical band gap) ($\Omega$) and exciton binding energy ($E_{b}$) of
single BN sheet. The available experimental data (Exp.) are also listed for
comparison.

|  $\Delta E$ ($E_{g}$) (eV) |  $E_{abs}$ (eV) | Band width (eV) |  $\Omega$ (eV) |  $E_{b}$ (eV)
---|---|---|---|---|---
| $K^{v}\rightarrow K^{c}$ | $K^{v}\rightarrow\Gamma^{c}$ |  |  $\pi$ Band |  $\sigma_{1}$ Band |  $\sigma_{2}$ Band | Exciton A (B)
LDA | 4.65 (4.65) | 4.77 | 4.52 | 5.45 | 5.80 | 7.65 | … | …
GW | 7.91 | 7.66 (7.66) | 7.72 | 5.93 | 6.41 | 8.20 | … | …
GW+BSE | … | … | … | … | … | … | 5.91 (6.76) | 1.81 (0.96)
Exp. | … | … | … | 5.8a | 6.5a | 8.2a | 6.05b, 6.03c, 6.1d | …

aReference [45]; bReference [46]; cReference [47]; dReference [48]

Figure 3(a) shows the calculated imaginary (absorptive) part ($\varepsilon"$)
of the dielectric function of the single BN sheet. We notice that the
$\varepsilon"$ spectrum looks very similar to that of the previous LDA
calculation [20], albeit with the onset of the absorption and peak positions
at much higher photon energies due to the GW quasiparticle corrections. For
optical excitations below 10 eV with the in-plane polarization of the electric
field, only transitions between the $\pi$ and $\pi^{*}$ bands are dipole
allowed. The oscillator strength of the dielectric function is mainly
contributed by the $\pi$ and $\pi^{*}$ bands. The onset energy of the
absorption spectrum at the GW-IPA level corresponds to the direct transition
at the $K$ point. On the other hand, the optical transition near the $M$ point
contributes to the largest peak at $\sim$9.0 eV. In other words, the
oscillator strength of the dielectric function stems primarily from the
optical transitions between two flat $\pi$ and $\pi^{*}$ bands along the high
symmetry $K-M$ line.

Remarkably, when the excitonic effect is included at the GW+BSE level, two
prominent absorption peaks [labelled A and B in Fig. 3(a)] occur in the
quasiparticle band gap and they are due to the excitation of the first two
bright excitions. Unlike bulk semiconductors such as CdS in which the exciton
peak usually appears a shoulder on the absorption edge [10], the A and B
exciton peaks in the single BN sheet are completely detached from the
quasiparticle absorption edge due to their gigantic binding energies of 1.81
and 0.96 eV (Table I), respectively. Note that the calculated large excitation
energy of the A exciton agrees well with the measured optical band gap (see
Table I).

The large exciton binding energies are the consequence of the weak screening
in 2D wide bandgap materials such as the single BN sheet, as mentioned in the
previous paragraph. Figure 4 shows the exciton A wavefunction
($|\psi(\textbf{r}_{e},\textbf{r}_{h})|^{2}$). Clearly, the exciton is
strongly localized, and the electron probability distribution beyond the
nearest neighbor is lower than $30\%$. Consequently, the exciton envelope
function of the A exciton is widespread in the BZ [37], resulting in the
mixing between the interband electron-hole excitations along the high-symmetry
$K-M$ line. The mixing strongly redistributes the oscillator strength of the
spectrum [see Eq. (4)]. The in-gap exciton states A and B carry almost all the
oscillator strength, and the absorption spectrum almost diminishes above the
quasiparticle band gap. The magnitude of the A exciton peak is more than three
times larger than the largest quasiparticle absorption peak at $\sim$9.0 eV at
the GW-IPA level.

Figure 3: (a) Imaginary part of the dielectric function and (b) shift current
conductivity of the single BN sheet from both GW+BSE (solid red line) and
GW+IPA (dashed blue line) calculations. The first two bright exciton peaks are
labelled by A and B. The green vertical line indicates the onset of the
continuum optical absorption (see Table I). Figure 4: The electron probability
distribution $|\psi(\textbf{r}_{e},\textbf{r}_{h})|^{2}$ [see Eq. (2)] of the
doubly degenerate A exciton in the single BN sheet with the hole position
($\textbf{r}_{h}$) fixed at the N atom (the pink sphere).

Figure 3(b) displays the calculated shift current conductivity of the BN
sheet. Clearly, the shift current conductivity spectra are rather similar to
the corresponding absorption spectra. In particular, when the excitonic effect
is included, two huge peaks occur within the quasiparticle band gap (red
line). Comparing the shift current conductivity to the imaginary part of the
dielectric function indicates that the two peaks are due to the A and B
exciton excitations. The height of the maximal peak at 5.91 eV is nearly eight
times larger than that of the shoulder at 8.15 eV from the GW-IPA calculation,
i.e., the effect of the exciton excitation not only create two in-gap peaks
but also greatly enhance the shift current. The significant increase in the
shift current is primarily due to the in-gap A exciton excitation, and is
mainly caused by the significant enhancement of the optical dipole matrix
element [see Eq. (4) and (10)] due to strong overlap of wavefunctions of the
electron and hole within the electron-hole pair.

###  III.2 Electronic properties of BN nanotubes

Figure 5: Quasiparticle bandgaps of the zigzag BN-NTs vs the tube diameter
from both LDA and GW calculations. Bandgaps from the previous LDA calculation
[20] are denoted by red open squares. The horizontal red and blue lines
indicate the quasiparticle band gaps of the single BN sheet from the LDA and
GW calculations, respectively.

Our LDA and GW calculations indicate that the zigzag BN-NTs considered here
are direct bandgap insulators with $\Gamma^{v}\rightarrow\Gamma^{c}$
transition, being consistent with previous LDA calculations (see, e.g., Ref.
[20]). The calculated band gaps of the zigzag BN-NTs are listed in Table 2 and
also displayed in Fig. 5. Figure 5 shows that for small diameter BN-NTs, both
GW and LDA band gaps get reduced dramatically as the diameter decreases.
Interestingly, the GW correction to the LDA band gap ($\Delta E_{g}$) is
almost independent of the diameter, being about 3.0 eV (Table 2 and Fig. 5).
The lowering of the band gap relative to the BN sheet already occurs at the
LDA level [20], and was attributed to the curvature-induced orbital
rehybridization. The $\pi^{*}$ and $\sigma^{*}$ orbitals are orthogonal to
each other in the flat BN sheet. However, such orthogonality does not hold in
the BN-NTs where the local curvature exists. Consequently, the $\pi^{*}$ and
$\sigma^{*}$ orbitals can hybridize and form ring-like charge distribution, as
shown in Fig. 6. The ring-like charge distribution can effectively reduce the
ionicity and hence lowers the energies of $\pi^{*}$ bands [50, 49, 51],
thereby reducing the quasiparticle bandgaps of the BN-NTs. As the diameter
increases, the effect of curvature-induced orbital rehybridization decreases
and the band gap converges to that of single BN sheet.

Figure 6: Charge density distribution (in units of $10^{-3}$ Å-3) of the
states in the vicinity of the conduction band minimum of the zigzag
[($5,0),(6,0),(7,0),(8,0)$] BN-NTs. B and N atoms are indicated by green and
gray spheres, respectively.

###  III.3 Absorption spectra of BN nanotubes

Figure 7: (a) The imaginary (absorptive) part of the dielectric function and
(b) shift current conductivity of the single-walled zigzag
[$(5,0),(6,0),(7,0),(8,0)$] BN-NTs and the single BN sheet from both the
GW+BSE (solid red lines) and GW+IPA (dashed blue lines) calculations. The
prominent peak labelled A, is due to the first bright exciton (exciton A).
Both the polarization direction of light electric field [(a) and (b)] and
shift current (b) are parallel to the tube axis ($y$-axis) in the BN-NTs (BN
sheet). For comparison, the experimental dielectric function of the multi-
walled BN-NT [53] is also displayed in (a). Table 2: Diameter ($D$), band gap
($E_{g}$), the onset of the continuum optical absorption ($E_{abs}$), exciton
excitation energy ($\Omega$), exciton binding energy ($E_{b}$) and the
difference between the LDA and GW band gaps ($\Delta E^{GW-LDA}_{g}$) of the
zigzag BN-NTs and also single BN sheet (Table I). For comparison, the measured
optical band gaps ($\Omega$) of large diameter BN-NTs are also listed.

|  |  |  |  | Exciton A
---|---|---|---|---|---
|  $D$ (Å) |  $E_{g}^{LDA}$ (eV) |  $E_{g}^{GW}$ (eV) |  $\Delta E^{GW-LDA}_{g}$ (eV) |  $E_{abs}$ (eV) |  $E_{b}$ (eV) |  $\Omega$ (eV)
BN-NT |  |  |  |  |  |  |
(5,0) | 4.11 | 2.13 | 5.19 | 3.06 | 7.06 | 2.23 | 4.83
(6,0) | 4.87 | 2.70 | 5.70 | 3.00 | 7.37 | 2.32 | 5.05
(7,0) | 5.65 | 3.31 | 6.31 | 3.00 | 7.61 | 2.27 | 5.34
(8,0) | 6.43 | 3.51 | 6.51 | 3.00 | 7.71 | 2.28 | 5.43
Exp.a | $<10$ | … | … | … | … | … | 5.51
Exp.b | 15-30 | … | … | … | … | … | 5.8 $\pm$ 0.2
Exp.c | 50 $\pm$ 10 | … | … | … | … | … | 5.82 $\pm$ 0.01
Exp.c | 600 $\pm$ 100 | … | … | … | … | … | 5.90 $\pm$ 0.01
sheet | $\infty$ | 4.65 | 7.66 | 3.01 | 7.72 | 1.81 | 5.91

aReference [16]; bReference [52]; cReference [17]

We show the calculated absorptive part ($\varepsilon"$) of the dielectric
function of the considered BN-NTs and also single BN sheet in Fig. 7(a). In
the BN-NTs, optical transitions between $\pi$ and $\pi^{*}$ bands dominate the
absorption spectra. First, compared with the single BN sheet, the energies of
the $\pi^{*}$ bands of the BN-NTs are lowered owing to the curvature-induced
orbital rehybridization, as discussed above. Consequently, the energies of the
$\varepsilon^{"}$ continuum onset at the GW-IPA level decrease with decreasing
diameter, as shown in Fig. 7 (a) and Table 2. Second, contrary to the moderate
diameter ($D$ >10 Å) BN-NTs where the $\varepsilon"$ spectrum consists of a
single distinct peak [20], the absorption spectra of the small diameter BN-NTs
considered here feature multiple peaks at the GW-IPA level. For example, the
$\varepsilon^{"}$ spectrum of BN-NT (5,0) possesses three distinct peaks at
7.4, 7.9, and 8.7 eV, respectively.

BN-NTs are 1D wide bandgap insulators. Consequently, due to much reduced
dielectric screening in these 1D wide bandgap materials, large exciton binding
energies arise. For small diameter BN-NT (5,0), in particular, the exciton
binding energy $E_{b}$ of peak A is 2.23 eV, which is larger than that (1.81
eV) of the 2D single BN sheet (Table I) and also that (0.63 eV) of the typical
2D transition metal dichalcogenide semiconductor $\text{MoS}_{2}$ [54]. It is
also much larger than that of 3D bulk $h$-BN (0.72 eV) [55] and 2H-MoS2 (0.05)
[56]. Note that the exciton binding energy of typical bulk semiconductor CdS
is as small as 0.03 eV [10]. The large $E_{b}$ results in the emergence of the
huge in-gap absorption peak which is completely detached from the onset of the
continuum absorption $\varepsilon^{"}$ spectrum of the BN-NT. Most of the
oscillator strength of $\varepsilon^{"}$ spectra is carried by the A exciton
state of the BN-NTs.

Table 2 shows the excitation energies ($\Omega$) of exciton A in BN-NTs. Our
calculated $\Omega$ converges to that of the single BN sheet and agrees well
with the experimental values [52, 17]. We also compare our calculated
$\varepsilon^{"}$ with the experimental one derived from electron energy-loss
spectra. Experimental $\varepsilon^{"}$ of medium diameter multi-walled BN-NTs
with small momentum transfer $q=0.1$ Å-1 is plotted in Fig. 7 (a) (orange
dotted line) [53]. The complex main peak at 5.5 eV in the experimental
$\varepsilon^{"}$ spectrum may consist of exciton A peaks from different walls
of the multi-walled BN-NT. The substantial broadening of the complex primary
peak could be due to the consequence of the finite exciton lifetime. The
shoulder from around 6.0 to 8.0 eV may contribute by minor peaks other than
exciton A peak.

In the low-dimensional BN systems, exciton A is strongly localized, as shown
in Fig. 4. The localized nature of exciton A suppresses the curvature effect
on the exciton binding energy as the tube diameter increases, and leads to the
fast convergence of the exciton binding energy [see Table 2] [14]. The red
shift of the A peak can be attributed to the curvature-induced orbital
rehybridization that lowers the energies of $\pi^{*}$ bands and reduces the
onset of the continuum $\varepsilon^{"}$ at the GW-IPA level.

###  III.4 Exciton shift current of BN nanotubes

In Fig. 7 (b), calculated shift current conductivity spectra of the BN-NTs are
plotted. Figure 7 shows that the shift current spectra are very similar to the
corresponding $\varepsilon^{"}$ spectra. This correlation is facilitated by
the dipole matrix elements present in both the expression for
$\varepsilon^{"}$ [Eq. (4)] and the expression for the shift current [Eq.
(9)]. In particular, as for the $\varepsilon^{"}$ spectrum, the gigantic
excitonic effect results in the occurrence of the huge in-gap shift current
peak from the GW-BSE calculation [Fig. 7 (b)]. The prominent shift current
peak in the BN-NTs and single BN sheet is due to the excitation of the A
exciton. Furthermore, the excitonic effect substantially increases the largest
shift current peak of all the BN-NTs by nearly a factor of three. This large
enhancement of the shift current by the A exciton excitation is mainly due to
the large enhancement of the dipole matrix element caused by the strong
overlap of the wavefunctions of the electron and hole in the A exciton state.

The optical responses of BN-NTs in the low-energy region is dominated by
$\pi$-$\pi^{*}$ transitions. As mentioned before, the curvature-induced
orbital rehybridization lowers the energy of the $\pi^{*}$ bands of BN-NTs. As
a result, the onset of the shift current spectra at the GW-IPA level decreases
with decreasing diameter, as shown in Fig. 7 (b). Similar to $\varepsilon^{"}$
spectra, shift current spectra of small diameter BN-NTs feature multiple peaks
at the GW-IPA level.

Figure 8:  The k-resolved shift vector ($R$) (a, b) amd shift current
conductivity ($\sigma$) (c, d) of the single BN sheet from the GW-IPA
calculation. The green lines in (a) and (b) indicate the 1D BZ of the ($6,0$)
and ($8,0$) BN-NTs, respectively.

Figure 7 (b) shows that the shift current spectra of the ($n,0$) BN-NTs from
both ab initio GW-IPA and GW-BSE calculations do not change sign as the tube
diameter (or $n$) increases, i.e., the direction of the calculated shift
current does not follow the simple rule of
$sgn(J_{\text{shift}})=\text{mod}(n,3)$ predicted by the previous tight-
binding model calculations [19, 18]. To better understand this important
finding, we compute the k-resolved shift conductivity and shift vector, which
is given by $R^{yyy}_{\text{GW-
IPA}}(\textbf{k})={\text{Im}}[\textbf{x}^{y}_{mn}[\textbf{x}^{y}_{nm}]_{;k_{y}}]/\left|\textbf{x}^{y}_{mn}\right|^{2}$
[57], of the single BN sheet, as displayed in Fig. 8 for two incident photon
energies of 8.0 eV (optical absorption edge) and 9.0 eV (the peak position) at
the GW-IPA level. Figure 8 indicates that the k-resolved shift current and
shift vector of the single BN sheet are positive in the whole BZ for both
selected photon energies. At the absorption edge (8.0 eV), the shift current
and vector are localized on the $K$ and $K^{\prime}$ points. As the energy
increases to the peak position (9.0 eV), the weights of the shift current and
shift vector extend toward the $M$ point and eventually toward the
$K$-$\Gamma$ symmetry line. The 1D BZ of one BN-NT consists of a few discrete
lines on the 2D BZ of the BN sheet [11], as displayed for the ($6,0$) and
($8,0$) BN-NTs in Fig. 8. Consequently, the direction of the shift current of
the zigzag ($n,0$) BN-NTs would remain unchanged when the tube diameter [or
the chirality ($n,0$)] varies, contrary to the prediction of the tight-binding
model calculations [19, 18].

We notice that previous tight-binding calculations also predicted that the
direction of the electric polarization of the zigzag BN-NTs would follow the
same rule of $sgn(J_{\text{shift}})=\text{mod}(n,3)$ as the shift current [18,
58]. This is not surprising since the electric polarization and shift current
are closely connected [59]. However, a subsequent $ab$ $initio$ calculation
showed that the electric polarization of the zigzag BN-NTs grows monotonically
as the tube diameter increases [60], i.e., the direction of the electric
polarization does not depend on the chiral index ($n$). Therefore, both the
previous [60] and present ab initio calculations demonstrate the importance of
the full ab initio calculations for the electric polarization and shift
current in the BN-NTs.

Furthermore, a close examination of the simple tight-binding model Hamiltonian
suggests that the erroneous prediction of the tight-binding model calculations
[18, 19] would stem from the fact that only $\pi$-bands in the single BN sheet
and hence the transitions between the lowest two azimuthal subbands in zigzag
BN-NTs were included, whereas the predominant exciton peak in BN-NTs is
composed of a coherent supposition of transitions from several different
subband pairs [15]. Consequently, the rehybridization effect, induced by the
curvature of small diameter BN-NTs, is missing in the simple tight binding
method [18, 19]. Indeed, previous DFT calculations showed that the
rehybridization effect plays a crucial role in the strong bandgap
renormalization for small diameter zigzag BN-NTs [20].

Let us now compare the magnitude of the excitonic photocurrents in the
considered BN systems with the observed shift currents in well-known
materials. First, the observed photocurrent in ferroelectric BaTiO3 above the
absorption edge under light intensity $I=0.5$ mW/cm2 and sample width $w=0.15$
cm is around 5$\times 10^{-13}$ A [61]. We find that under the same
conditions, the A exciton shift current of a single BN sheet would reach
$5.0\times 10^{-12}$ A, and that from a thin film consisting of, e.g., the
($5,0$) BN-NT array with an intertubular distance ($d$) of 3.28 Å would be as
large as $\sim$10 mA. Second, the effective conductivity of the A exciton
shift current in the considered BN systems (Fig. 7) is about two-orders of
magnitude larger than the excitonic shift conductivity ($\sim 0.2$ $\mu$A/V2)
observed in semiconductor CdS [10], and also more than five times larger than
the largest observed shift conductivity ($\sim 10.0$ $\mu$A/V2) in
ferroelectric SbSI [62].

##  IV CONCLUSIONS

Using the newly developed computationally efficient algorithms, we have
performed the state-of-the-art ab initio GW-BSE culations to investigate the
exciton shift current as well as the electronic and optical properties of the
zigzag [$(5,0),(6,0),(7,0),(8,0)$] BN-NTs and also the single BN sheet. First
of all, we find a giant in-gap peak in both the shift current and optical
absorption spectra in all the studied BN systems due to the excitation of the
A exciton. This excitonic peak is nearly three times higher than that in the
continuum due to the excitation of the free electron-hole pairs (Fig. 7), and
may be attributed to the gigantic enhancement of the optical dipole matrix
element by the A exciton resonance. Second, our ab initio calculations show
that the exciton excitation energies and also the onset of the continuum
spectra decrease significantly with the decreasing diameter due to the
curvature-induced orbital rehybridization in small diameter zigzag BN
nanotubes. This orbital rehybridization lowers the ionicity of the $\pi^{*}$
bands by the formation of ring-like charge distribution inside the BN-NTs
(Fig. 6), thus reducing the energy of the $\pi^{*}$ bands and hence the
bandgap. The quasiparticle GW correction to the bandgap is almost independent
of the tube diameter.

Third, we find that the direction of the shift current in the BN-NTs is
independent of the tube chirality ($n,0$) (or diameter), contrary to the
simple detrimental rule of $sgn(J_{\text{shift}})=\text{mod}(n,3)$ reported by
previous model Hamiltonian studies [18, 19]. Importantly, this implies that in
a bundle of aligned zigzag BN-NTs, the contributions of the BN-NTs to the
shift current would be additive rather than cancelling each other as the
simple rule suggests [18, 19]. Finally, the effective exciton shift current
conductivity is nearly ten times larger than the largest shift conductivity
observed in ferroelectric semiconductors [62]. Furthermore, strongly bounded
excitons with a large binding energy in the considered BN systems would impede
thermal dissociation into free electron-hole pairs. All these properties would
make the zigzag BN-NTs excellent candidates for experimentally studying the
exciton shift current and also for promising applications in nano-scale
optoelectronic devices.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the support from the National Science and
Technology Council and the National Center for Theoretical Sciences of The
R.O.C. The authors also thank the National Center for High-performance
Computing (NCHC) in Taiwan for providing computational and storage resources.

## References

  * [1] W. Kraut and R. v. Baltz, Anomalous bulk photovoltaic effect in ferroelectrics: A quadratic response theory: A quadratic response theory, Phys. Rev. B 19, 1548 (1979).
  * [2] B. I. Sturman and V. M. Fridkin, The Photovoltaic and Photorefractive Effects in Noncentrosymmetric Materials, (Gordon and Breach Science Publishers, 1992).
  * [3] C. Aversa and J. E. Sipe, Nonlinear optical susceptibilities of semiconductors: Results with a length-gauge analysis, Phys. Rev. B 52, 14636 (1995).
  * [4] J. E. Sipe and A. I. Shkrebtii, Second-order optical response in semiconductors, Phys. Rev. B 61, 5337 (2000).
  * [5] A. Bhatnagar, A. R. Chaudhuri, Y. H. Kim, D. Hesse and M. Alexe, Role of domain walls in the abnormal photovoltaic effect in BiFeO3, Nature Commun. 61, 5337 (2000).
  * [6] M. A. Green and S. P. Bremner, Energy conversion approaches and materials for high-efficiency photovoltaics, Nat. Mater. 16, 23 (2016).
  * [7] T.-Y. Cai, S.-C. Liu, S. Ju, C.-Y. Liu and G.-Y. Guo, Multiferroic double perovskites ScFe1-xCrxO3 (1/6 $\leq x\leq$ 5/6) for highly efficient photovoltaics and spintronics, Phys. Rev. Appl. 8, 034034 (2017).
  * [8] T. Morimoto and N. Nagaosa, Topological aspects of nonlinear excitonic processes in noncentrosymmetric crystals, Phys. Rev. B 94, 035117 (2016).
  * [9] Y. H. Chan, D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Giant exciton-enhanced shift currents and direct current conduction with subbandgap photo excitations produced by many-electron interactions, Proc. Natl. Acad. Sci. USA 118, e1906938118 (2021).
  * [10] M. Sotome, M. Nakamura, T. Morimoto, Y. Zhang, G.-Y. Guo, M. Kawasaki, N. Nagaosa, Y. Tokura, and N. Ogawa, Terahertz emission spectroscopy of ultrafast exciton shift current in the noncentrosymmetric semiconductor CdS, Phys. Rev. B 103, L241111 (2021).
  * [11] R. Saito, G. Dresselhaus, and M. S. Dresselhaus, Physical Properties of Carbon Nanotubes, (Imperial College, London, 1998).
  * [12] G. Y. Guo,and J. C. Lin, Second-harmonic generation and linear electro-optical coefficients of BN nanotubes, Phys. Rev. B 72, 075416 (2005).
  * [13] R. S. Lee, J. Gavillet, M. Lamy de la Chapelle, A. Loiseau, J. L. Cochon, D. Pigache, J. Thibault, and F. Willaime, Catalyst-free synthesis of boron nitride single-wall nanotubes with a preferred zig-zag configuration, Phys. Rev. B 64, 121405(R) (2001).
  * [14] L. Wirtz, A. Marini, and A. Rubio, Excitons in Boron Nitride Nanotubes: Dimensionality Effects, Phys. Rev. Lett. 96, 126104 (2006).
  * [15] C. H. Park, C. D. Spataru, and S. G. Louie, Excitons and Many-Electron Effects in the Optical Response of Single-Walled Boron Nitride Nanotubes, Phys. Rev. Lett. 96, 126105 (2006).
  * [16] L. H. Li, Y. Chen, M. Y. Lin, A. M. Glushenkov, B. M. Cheng, and J. Yu, Single deep ultraviolet light emission from boron nitride nanotube film, Appl. Phys. Lett. 97, 141104 (2010).
  * [17] J. Yu, D. Yu, Y. Chen, H. Chen, M. Y. Lin, B. M. Cheng, J. Li, and W. H. Duan, Narrowed bandgaps and stronger excitonic effects from small boron nitride nanotubes, Chem. Phys. Lett. 476, 240-243 (2009).
  * [18] P. Král, E. J. Mele, and D. Tománek, Photogalvanic Effects in Heteropolar Nanotubes, Phys. Rev. Lett. 85, 1512 (2000).
  * [19] S. Konabe, Exciton effect on shift current in single-walled boron-nitride nanotubes, Phys. Rev. B 103, 075402 (2021).
  * [20] G. Y. Guo and J. C. Lin, Systematic $ab$ $initio$ study of the optical properties of BN nanotubes, Phys. Rev. B 71, 165402 (2005).
  * [21] R. Fei, L. Z. Tan, and A. M. Rappe, Shift-current bulk photovoltaic effect influenced by quasiparticle and exciton, Phys. Rev. B 101, 045104 (2020).
  * [22] A. Taghizadeh, and T. G. Pedersen, Gauge invariance of excitonic linear and nonlinear optical response, Phys. Rev. B 97, 205432 (2018).
  * [23] R. W. Boyd, Nonlinear Optics, (Elsevier, Amsterdam, 2003).
  * [24] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47, 558 (1993).
  * [25] G. Kresse and J. Furthmuller, Efficient interactive schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * [26] P. Giannozzi, et al., QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter 21, 395502 (2009).
  * [27] D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013).
  * [28] H. J. Monkhorst and J. D. Pack, Special points for Brillonin-zone integrations, Phys. Rev. B 13, 5188 (1976).
  * [29] M. S. Hybertsen and S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies, Phys. Rev. B 34, 5390 (1986).
  * [30] M. Rohlfing and S. G. Louie, Electron-hole excitations and optical spectra from first principles, Phys. Rev. B 62, 5188 (2000).
  * [31] J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, S. G. Louie, BerkeleyGW: A massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures, Comp. Phys. Commun. 183, 1269-1289 (2012).
  * [32] F. H. d. Jornada, D. Y. Qiu, and S. G. Louie, Nonuniform sampling schemes of the Brillouin zone for many-electron perturbation-theory calculations in reduced dimensionality, Phys. Rev. B 95, 035109 (2017).
  * [33] S. Ismail-Beigi, Truncation of periodic image interactions for confined systems, Phys. Rev. B 73, 233103 (2006).
  * [34] J. Ibanez-Azpiroz, S. S. Tsirkin, and I. Souza, Ab initio calculation of the shift photocurrent by Wannier interpolation, Phys. Rev. B 97, 245143 (2018).
  * [35] I. Souza, J. Íñiguez, and D. Vanderbilt, Dynamics of Berry-phase polarization in time-dependent electric fields, Phys. Rev. B 69, 085106 (2004).
  * [36] K. S. Virk and J. E. Sipe, Semiconductor optics in length gauge: A general numerical approach, Phys. Rev. B 76, 035213 (2007).
  * [37] F. Zhang, C. S. Ong, J. W. Ruan, M. Wu, X. Q. Shi, Z. K. Tang, and S. G. Louie, Intervalley Excitonic Hybridization, Optical Selection Rules, and Imperfect Circular Dichroism in Monolayer h-BN, Phys. Rev. Lett. 128, 047402 (2022).
  * [38] H. Ahin, S. Cahangirov, M. Topsakal, E. Bekaroglu, E. Akturk, R. T. Senger, and S. Ciraci, Monolayer honeycomb structures of group-IV elements and III-V binary compounds: First-principles calculations, Phys. Rev. B 80, 155453 (2009).
  * [39] F. A. Rasmussen, P. S. Schmidt, K. T. Winther, and K. S. Thygesen, Efficient many-body calculations for two-dimensional materials using exact limits for the screened potential: Band gaps of MoS2, h-BN, and phosphorene, Phys. Rev. B 94, 155406 (2016).
  * [40] I. Guilhon, M. Marques, L. K. Teles, M. Palummo, O. Pulci, Silvana Botti, and F. Bechstedt, Out-of-plane excitons in two-dimensional crystals, Phys. Rev. B 99, 161201 (2019).
  * [41] H. Mishra and S. Bhattacharya, Giant exciton-phonon coupling and zero-point renormalization in hexagonal monolayer boron nitride, Phys. Rev. B 99, 165201 (2019).
  * [42] F. Ferreira, A. J. Chaves, N. M. R. Peres, and R. M. Ribeiro, Excitons in hexagonal boron nitride single-layer: a new platform for polaritonics in the ultraviolet, J. Opt. Soc. Am. B 36, 674-683 (2019).
  * [43] T. Galvani, F. Palean, H. P. C. Miranda, A. Molina-Sanchez, L. Wirtz, S. Latil, H. Amara and F. Ducastelle, Excitons in boron nitride single layer, Phys. Rev. B 94, 125303 (2016).
  * [44] H. C. Hsueh, G. Y. Guo, and S. G. Louie, Excitonic effects in the optical properties of a SiC sheet and nanotubes, Phys. Rev. B 84, 085404 (2011).
  * [45] A. Nagashima, N. Tejima, Y. Gamou, T. Kawai, and C. Oshima, Electronic dispersion relations of monolayer hexagonal boron nitride formed on the Ni(111) surface, Phys. Rev. B 51, 4606 (1995).
  * [46] C. Elias, P. Valvin, T. Pelini, A. Summerfield, C. J. Mellor, T. S. Cheng, L. Eaves, C. T. Foxon, P. H. Beton, S. V. Novikov, B. Gil and G. Cassabois, Direct band-gap crossover in epitaxial monolayer boron nitride, Nat. Commun. 10, 2639 (2019).
  * [47] X. Li, H. Qiu, X. F. Liu, J. Yin, and W. Guo, Wettability of Supported Monolayer Hexagonal Boron Nitride in Air, Adv. Funct. Mater. 27, 1603181 (2017).
  * [48] R. J. P. Roman, F. J. R. Costa Costa, A. Zobelli, C. Elias, P. Valvin, G. Cassabois, B. Gil, A. Summerfield, T. S. Cheng, C. J. Mellor, P. H. Beton, S. V. Novikov and L. F. Zagonel, Band gap measurements of monolayer h-BN and insights into carbon-related point defects, 2D Mater. 8, 044001 (2021).
  * [49] C. Thomsen, S. Reich, J. Maultzsch, Carbon Nanotubes: Basic Concepts and Physical Properties, (WILEY‐VCH, 2004).
  * [50] X. Blase, A. Rubio, S. G. Louie, and M. L. Cohen, Stability and Band Gap Constancy of Boron Nitride Nanotubes, Europhys. Lett. 28, 335 (1994).
  * [51] A. Rubio, J. L. Corkill, and M. L. Cohen, Theory of graphitic boron nitride nanotubes, Phys. Rev. B 49, 5081 (1994).
  * [52] R. Arenal, O. Stéphan, M. Kociak, D. Taverna, A. Loiseau, and C. Colliex, Electron Energy Loss Spectroscopy Measurement of the Optical Gaps on Individual Boron Nitride Single-Walled and Multiwalled Nanotubes, Phys. Rev. Lett. 95, 127601 (2005).
  * [53] G. G. Fuentes, E. Borowiak-Palen, T. Pichler, X. Liu, A. Graff, G. Behr, R. J. Kalenczuk, M. Knupfer, and J. Fink, Electronic structure of multiwall boron nitride nanotubes, Phys. Rev. B 67, 035429 (2003).
  * [54] D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Screening and many-body effects in two-dimensional crystals: Monolayer MoS2, Phys. Rev. B 93, 235435 (2016).
  * [55] B. Arnaud, S. Lebegue, P. Rabiller and M. Alouani, Huge excitonic effects in layered hexagonal boron nitride, Phys. Rev. Lett. 96, 026402 (2006)
  * [56] A. R. Beal, J. C. Knights, and W. Y. Liang, Transmission spectra of some transition metal dichalcogenides. II. Group VIA: Trigonal prismatic coordination, J. Phys. C: Solid State Phys. 5, 3540 (2006)
  * [57] A. Strasser, H. Wang, and X. Qian, Nonlinear Optical and Photocurrent Responses in Janus MoSSe Monolayer and MoS$2$-MoSSe van der Waals Heterostructure, Nano Lett. 22, 4145–4152 (2022).
  * [58] E. J. Mele and P. Král, Electric Polarization of Heteropolar Nanotubes as a Geometric Phase, Phys. Rev. Lett. 88, 056803 (2002).
  * [59] B. M. Fregoso, T. Morimoto and J. E. Moore, Quantitative relationship between polarization differences and the zero-averaged shift photocurrent, Phys. Rev. B 96, 075421 (2017)
  * [60] S. M. Nakhmanson, A. Calzolari, V. Meunier, J. Bernholc, and M. B. Nardelli, Spontaneous polarization and piezoelectricity in boron nitride nanotubes, Phys. Rev. B 67, 235406 (2003).
  * [61] W. T. H. Koch, R. Munser, W. Ruppel and P. Wurfel, Anomalous photovoltage in BaTiO3, Ferroelectrics 13, 305 (1976)
  * [62] M. Sotome, M. Nakamura, J. Fujioka, M. Ogino, Y. Kaneko, T. Morimoto, Y. Zhang, M. Kawasaki, N. Nagaosa, Y. Tokura, and N. Ogawa, Spectral dynamics of shift current in ferroelectric semiconductor SbSI, Proc. Natl. Acad. Sci. USA 116, 1929 (2019).


11institutetext: Universidade Federal de Ouro Preto, Ouro Preto 35402-163,
Brasil
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>
22institutetext: Cloudwalk, Inc, São Paulo, 05425-070 São Paulo, Brasil
22email<EMAIL_ADDRESS>
# FeatGeNN: Improving Model Performance for Tabular Data with Correlation-
based Feature Extraction††thanks: The author would like to acknowledge FAPEMIG
(Fundação de Amparo à Pesquisa do Estado de Minas Gerais), CAPES (Coordenação
de Aperfeiçoamento de Pessoal de Nível Superior), UFOP (Universidade Federal
de Ouro Preto) and Cloudwalk, Inc, for the financial support which has been
instrumental in the successful execution of our research endeavors.

Sammuel Ramos Silva 1122 0000-0002-2274-1390 Rodrigo Silva 11
0000-0003-2547-3835

###### Abstract

Automated Feature Engineering (AutoFE) has become an important task for any
machine learning project, as it can help improve model performance and gain
more information for statistical analysis. However, most current approaches
for AutoFE rely on manual feature creation or use methods that can generate a
large number of features, which can be computationally intensive and lead to
overfitting. To address these challenges, we propose a novel convolutional
method called FeatGeNN that extracts and creates new features using
correlation as a pooling function. Unlike traditional pooling functions like
max-pooling, correlation-based pooling considers the linear relationship
between the features in the data matrix, making it more suitable for tabular
data. We evaluate our method on various benchmark datasets and demonstrate
that FeatGeNN outperforms existing AutoFE approaches regarding model
performance. Our results suggest that correlation-based pooling can be a
promising alternative to max-pooling for AutoFE in tabular data applications.

###### Keywords:

Automated Feature Engineering feature creation correlation-based pooling
tabular data machine learning.

##  1 Introduction

Creating effective features is a crucial aspect of machine-learning projects.
Essentially, it involves deriving new features from existing data to train a
model or extract more information for statistical analysis. Discovering novel
features from raw datasets is often the key to improving model performance
[1].

Traditionally, feature creation is a manual process that heavily relies on an
analyst’s domain knowledge and programming skills. However, this approach can
be limiting, as an analyst’s intuition and expertise often influence the
features created. To overcome these limitations, researchers have been
exploring the field of Automated Feature Engineering (AutoFE). AutoFE aims to
automate the feature creation process, enabling the discovery of more complex
and effective features without relying solely on human input.

Automated feature engineering methods involve applying transformations to raw
data to create new features. One commonly used technique is the expansion-
reduction method [4], which generates a large number of features and then
applies a feature selection algorithm to reduce their dimensionality. During
the expansion phase, various transformations, such as logarithmic, max/min, or
sum, are applied to the raw data. In the reduction phase, a feature selection
method is utilized to identify the most effective set of features, which can
significantly enhance a model’s performance.

The possible number of transformation operations that can be performed on
already-transformed features is practically infinite, which leads to an
exponential increase in the feature space. This issue can cause a problem in
reducing the number of feature evaluations required. To address this issue,
researchers have proposed adaptive methods for AutoFE. For instance, Khurana
et al. [5] introduced a Q-learning agent capable of performing feature
transformation search, achieving higher performance but still generating a
large number of features. In another study [6], a Multi-Layer Perceptron (MLP)
was trained to suggest the best transformations for each raw feature,
resolving the problem of excessive feature generation. More recently, DIFER
[7], a gradient-based method for differentiable AutoFE, has demonstrated
superior performance and computational efficiency compared to other
approaches, although it still requires significant computation

In recent years, the use of deep neural networks (DNNs) has become
increasingly widespread across a range of fields, such as computer vision and
natural language processing [19, 18]. Typically, these models extract new
features by feeding input features into the hidden layers of a DNN. While this
approach is effective in capturing complex interactions between implicit and
explicit features, it may not always generate useful new features due to a
lack of relevant interactions in the dataset [9]. Moreover, most existing
works use max-pooling in the pooling layer, which may not be optimal for
tabular data because it does not preserve the order and context of features in
the data matrix. Additionally, max-pooling is intended to identify the most
significant features within an image, which may not always be relevant or
effective for tabular data.

To address the limitations of existing AutoFE methods, we propose FeatGeNN, a
convolutional approach that leverages correlation as a pooling function to
extract and generate new features. FeatGeNN first applies convolutional
filters to the raw data to extract high-level representations. Then, instead
of using traditional pooling functions like max or average pooling, it
computes the correlation between the extracted features, which helps to
identify the most informative features. The selected features are then passed
through a multi-layer perceptron (MLP) to create the final set of new
features. Preliminary results indicate that FeatGeNN outperforms existing
AutoFE methods in both the number of generated features and model performance,
demonstrating its potential as a potent tool for creating features in machine
learning.

##  2 Related work

The main goal of feature engineering is to transform raw data into new
features that can better express the problem to be solved. Training a model
with the generated features can increase the performance of the model.
However, the process of feature engineering can be limited by the expertise,
programming skills, and intuition of the person working with the data. For
this reason, AutoFE approaches have recently gained attention.

The authors of [7] propose a differentiable AutoML model that efficiently
extracts low and high-order features. The model includes three steps:
Initialization, Optimizer Training, and Feature Evaluation. In initialization,
features are constructed randomly and evaluated using a machine-learning model
in the validation set. In training the optimizer, a tree-like structure is
created with an encoder, a predictor, and a decoder, called a parse tree. The
encoder maps the post-order traversal string to a continuous space, the
predictor is a 5-layer MLP that maps the representation to the score computed
by a machine learning model, and the decoder maps the embedding to the
discrete feature space. In the final step of feature evolution, the best n
features are selected and optimized using a gradient-based approach.

In [4], the authors present an algorithm that uses mathematical functions to
generate new features for relational databases. The algorithm begins by
identifying the entities that make up the database and defines a set of
mathematical functions that are applied at both the entity level and the
relational level. The proposed approach first enumerates all possible
transformations on all features and then directly selects features based on
their impact on model performance. However, due to the potentially large
number of features generated, it is necessary to perform feature selection and
dimensionality reduction to avoid overfitting and improve the interpretability
of the model.

In [6], the authors propose a novel model for feature engineering in
classification tasks that can generalize the effects of different feature
transformations across multiple datasets. The model uses an MLP for each
transformation to predict whether it can produce more useful features than the
original set. The Quantile Sketch Array (QSA) achieves a fixed-size
representation of feature values to handle features and data of different
lengths. The QSA uses Quantile Data Sketch to represent feature values
associated with a class label.

The authors of [12] have proposed an RNN-based approach to address the feature
explosion problem in feature engineering and support higher-order
transformations. Their architecture uses an RNN to generate transformation
rules with a maximum order for each raw feature within a fixed time limit. For
datasets with multiple raw features, the authors use multiple RNNs as
controllers to generate transformation rules for each feature. The transformed
features are evaluated using a machine learning algorithm and the controller
is trained using policy gradients. The model includes two special unary
transformations: ”delete” and ”terminate”, which remove a feature and
terminate the current transformation, respectively, to determine the most
appropriate transformation order.

In [5] they propose a heuristic model for automating feature engineering in
supervised learning problems. Their model is based on a tree structure, where
the raw dataset is the root, each node is a transformed dataset, and the edges
represent the transformation functions. The goal is to find the node with the
highest score, reducing the feature construction problem to a search problem.

The authors present three exploration strategies to traverse the tree. The
first is depth-first traversal,” in which a random transformation is applied
to the root and the algorithm then explores a branch until there is no further
improvement. Then it chooses another node with the highest score and starts
the process again. The second is the ”Global Traversal”, where a global search
is performed to find the most promising node out of all the nodes explored so
far. The third is ”Balanced Traversal”, in which the algorithm chooses either
an exploration or exploitation strategy at each step based on a time or node
budget. To handle the explosive growth of columns, feature selection is
required as they grow. Cognito allows the selection of features after each
transformation to clean up the dataset and ensure a manageable size. In
addition, at the end of the model execution, the algorithm performs another
feature selection for all columns in the dataset, including the newly created
columns.

AutoFeat is a method presented in [11] that generates and selects non-linear
input features from raw inputs. The method applies a series of transformations
to the raw input and combines pairs of features in an alternating multi-step
process to generate new features. However, this leads to an exponential
increase in the size of the feature space, so a subsampling procedure is
performed before computing new features. The authors have shown that two or
three steps of the feature technique are usually sufficient to generate new
features.

After feature engineering, the new dataset has a higher number of features
than the original dataset. To reduce the dimensionality, the authors developed
a feature selection procedure. First, they remove new features that are highly
correlated with the original or simpler features. Then they apply a wrapper
method with L1-regular linear models to select the most informative and non-
redundant features from the dataset. In the end, only a few dozen features are
retained and used after the feature creation and selection process.

In [27], autolearn is proposed, a learning model based on regression between
pairs of features and aimed at discovering patterns and their variations in
the data. The method selects a small number of new features to achieve the
desired performance. The proposed method consists of four phases: Pre-
processing to reduce dimensionality, where the authors perform feature
selection based on information gain (IG); Mining of correlated features to
define and search for pairwise correlated features, where the distance
correlation [8] is calculated to determine if there is an interesting
predictive relationship between a pair of features; Feature generation, where
regularized regression algorithms are used to search for associations between
features and generate new features; and Feature selection, where features that
do not add new information to the dataset are discarded.

The authors of [3] have proposed a novel model that achieves both memorization
and generalization by simultaneously training a linear model component and a
neural network component. The model consists of two components: The Wide
component, which is a generalized linear model of the form $yW^{t}xb$, where
$y$ denotes prediction, $x$ denotes features, $w$ denotes model parameters and
$b$ denotes bias. The input features can be either raw or transformed, the
most important transformation being the cross-product transformation; and the
Deep component, which is a feed-forward neural network. For categorical
features, an embedding is created, which is then added to the dataset and fed
into the network.

The authors of [2] have proposed a model for predicting CTR that can handle
interactions between low and high-order features by introducing a
factorization-machine (FM) based neural network. The model consists of two
parts: the FM component, which generates low-order features and can generate
interactions between 1st and 2nd-order features with low computational cost,
and the deep component, a feed-forward neural network that learns interactions
between higher-order features. The input to the network is a high-dimensional
vector of sparse data containing categorical and continuous variables as well
as grouped fields.

The FGCNN is another approach proposed in CTR for prediction[9]. This model
consists of two components, namely the Feature Generation and the Deep
Classifier. The Feature Generation component uses the mechanisms inherent in
the Convolutional Neural Network (CNN) and the Multilayer Perceptron (MLP) to
identify relevant local and global patterns in the data and generate new
features. The Deep Classifier component then uses the extended feature space
to learn and make predictions.

Our work introduces a CNN-based model with correlation-pooling for extracting
high-order features and improving model performance. Unlike traditional
pooling functions such as max-pooling, which focus on selecting the maximum
value within a pooling region, correlation-pooling considers the linear
relationships between features in the data matrix. It measures the correlation
coefficient between the features and aggregates them based on their
correlation values to capture the interdependencies and patterns in the data.
By incorporating correlation-based pooling into the feature extraction
process, FeatGeNN can effectively extract high-order features that reflect the
underlying relationships among input variables. Our proposed method achieves
competitive results on a range of problems, suggesting that correlation-based
pooling is a promising technique for working with tabular data in neural
networks.

##  3 Proposed Approach

In this section, we describe the proposed Feature Generation with Evolutionary
Convolutional Neural Networks (FeatGeNN) model in detail.

###  3.1 Problem Formulation

Given an dataset $D=\left\langle F,tg\right\rangle$, where
$F=\left\\{f_{1},f_{2},...,f_{n}\right\\}$ are the raw features and $tg$ the
target vector. We denote as $L_{E}^{M}(D,tg)$ the performance of the machine
learning model $M$ that is learned from $D$ and measured by the evaluation
metric $E$ (e.g. accuracy). In addition, we transform a raw set of features
$D$ into $D_{new}$ by applying a set of transformation functions
$T=\left\\{t_{1},t_{2},...,t_{n}\right\\}$.

Formally, the goal of the AutoFE is to search the optimal transformed feature
set $D^{*}$ where $L_{E}^{M}(D^{*},tg)$ is maximized.

###  3.2 FeatGeNN Model

Figure 1: The FeatGeNN process.

In this study, we use a convolutional neural network to extract features that
can improve the performance of a machine learning model (i.e., Random Forest).
As explained earlier, using an MLP alone to generate new features would not
result in a good set of new features. The reason for this is the relationship
between the number of informative interactions between features and the total
number of features in the feature space. Also, using a CNN alone might not
lead to good performance because a CNN only considers local interactions and
does not consider many important global interactions between features [9].

To overcome this problem, we use an architecture that combines the MLP with
the CNN. The FeatGeNN model includes two main blocks, namely local feature
extraction and global feature generation (Figure 1). The first block attempts
to identify the most informative interactions between local features, while
the second block generates new features from the features extracted by the
local feature extraction block and combines them globally.

The Local Feature Extraction block includes two main operations, namely
Pooling and Convolution. Among these operations, the pooling operation plays a
crucial role in reducing dimensionality and preserving the most informative
features for subsequent layers. In previous work on feature generation for
tabular data with CNN, max-pooling was mainly used. However, we found that
using max-pooling for tabular data may not give the desired result because the
model may not compare closely related features, thus affecting the features
generated by the model. Therefore, we propose the use of correlation-pooling
to address this issue.

Figure 2: Correlation-Pooling process.

In correlation pooling, the variant of pooling used in our Local Feature
Extraction block, uses Pearson correlation[29] to group features that are
highly correlated. By grouping these features, correlation-pooling can
preserve the relationship between closely related features and thus improve
the quality of the features extracted by the CNN model. This is in contrast to
max-pooling, which preserves only the most dominant feature in a group and may
ignore other relevant features that are closely related. Therefore, by
incorporating Pearson correlation in the pooling operation, correlation-
pooling can effectively circumvent the limitation of max-pooling and help
generate more informative features for subsequent layers in the CNN model. The
Pearson correlation coefficient can be formulated as follows for our problem:

$r=\frac{n\sum xy-(\sum x)(\sum y)}{\sqrt{[n\sum x^{2}-(\sum x)^{2}][n\sum
y^{2}-(\sum y)^{2}]}}$ (1)

where x and y represent the values of the features X and Y respectively, and
$X,Y\in F$, where F is the set of all features. The variable n denotes the
number of samples in the dataset D.

To avoid having to run the Pearson algorithm twice, we have introduced an
iterative calculation of the Pearson coefficient. This means that at the
current stage of model development, we compute the Pearson coefficient r to
perform the pooling operation for the subsequent evolutionary generation of
the model. To reduce the computations required, we also added a threshold to
limit the number of data sent to the correlation calculation, i.e., a model
can only use 70% of the data to calculate the correlation value for the
features.

While the Pearson correlation is a statistical measure that describes the
linear relationship between two variables, it is not suitable for analyzing
relationships between more than two characteristics. To overcome this
limitation, we use the multivariate correlation matrix, which consists of
pairwise Pearson correlation coefficients between all pairs of variables. This
matrix allows us to analyze relationships between multiple variables and
identify the most highly correlated variables. The overall correlation value
for the feature f can be formulated as follows:

$CS_{f}=\frac{\sum_{k}^{N}r_{fk}}{N}$ (2)

where $CS_{f}$ is the correlation score for the feature f, $r_{fk}$ represent
the person correlation score for the feature tuple (f,k) and N the total
number of feature in the dataset.

In the Global Feature Generation block, an MLP is utilized to merge the
features extracted from the Local Feature Extraction block and generate novel
features. These novel features are then appended to the original dataset and
used in the machine-learning model.

###  3.3 Evolution Process

In this work, we adopt an evolution process for conducting AutoFE, as depicted
in Figure 3. This process involves three distinct steps: (1) Feature
Selection, (2) Population Initialization, and (3) Feature Evolution.

Figure 3: Feature Evolution process.

The first step of our proposed approach is to reduce the combination of
uncorrelated and redundant features using a feature selection method. We used
the Maximum Relevance-Minimum Redundancy (MRMR) [16] method for this purpose.
By minimizing the combination of such features, we aim to reduce the
introduction of noise into the model and improve the quality of the features
generated by the CNN model.

During the population initialization step, we generate a population of the CNN
model $POP$ that is evolved in the Feature Evolution step. To evaluate this
initial population, we use a machine learning model on the dataset resulting
from step (1). Specifically, we take the set of features F* from the Feature
Selection step and input them into the CNN model p (where $p\in POP$) to
generate n new features f. These newly created features are concatenated with
the original dataset D to create a new dataset $D*\left\\{F\cup f\right\\}$,
which is then evaluated by the machine learning model $L^{m}$ to obtain a
score $S_{p}$.

In the trait evolution step, a genetic algorithm [17] is used to evolve the
population and identify the most effective traits to improve the performance
score obtained by $L^{m}$. During each epoch of the genetic algorithm, for
each model $p$ that is not part of the elite group E (where $Eis\in POP$), a
crossover is performed between its weights and those of a second model
$p^{\prime}$, which is selected using a round-robin tournament [28]. Following
the crossover process, the offspring generated by this operation can be
subjected to mutation. The features produced by the offspring are then
evaluated, as described in the initialization of the population initialization
step. If the score obtained by $L^{m}$ is better than the current score for p
or if depreciation is allowed, the offspring replaces the current model p and
the score is updated.

##  4 Results

In this section, we aim to answer the following research questions:

  * •

RQ1: How effective is correlation-pooling compared to Max-Pooling?

  * •

RQ2: Study of the impact of the number of data on the correlation pooling
computation?

  * •

RQ3: How effective is the proposed FeatGeNN approach? (Comparison with
literature)

###  4.1 Experimental Setup

To evaluate the performance of the FeatGeNN model, on classification problems,
6 classification datasets from the UCI repository, which were used in the
state-of-the-art methods [7] [12], were selected. The description of each
dataset in terms of the Number of Features and Number of Samples is presented
in Table 1.

Table 1: Statistics of the benchmarks used to perform the evaluation of the FeatGeNN features. Datasets | Samples | Features
---|---|---
SpamBase | 4601 | 57
Megawatt1 | 253 | 37
Ionosphere | 351 | 34
SpectF | 267 | 44
Credit_Default | 30000 | 25
German Credit | 1001 | 24

In our experiments, we use the f1-score as the evaluation measure, which is
also commonly used in the related works [12] and [7]. The threshold for
questions RQ1 and RQ2 was set at 80% of the available data in the dataset. To
ensure robustness and reliability, we use 5-fold cross-validation, in which
the dataset is divided into five subsets or folds and the evaluation is
performed five times, with each fold serving once as a test set. This approach
helps mitigate the effects of data variability and provides a more
comprehensive assessment of the model’s performance. As for the chosen
algorithm, we use Random Forest as the base method in all our experiments.
Random Forest is a popular and widely used ensemble learning method known for
its robustness and ability to handle different types of data.

###  4.2 Effectiveness of Correlation-Pooling vs. Max-Pooling (RQ1)

In this subsection, this experiment aims to answer: Can our FeatGeNN with
Correlation-Pooling achieve competitive results compared to the version with
Max-Pooling? Table 2 shows the comparison results in terms of F1 score. The
results show that the FeatGeNN with correlation-pooling outperforms the
version with max-pooling in most datasets. The only exceptions are the
Megawatt1 and Credit_Default datasets, where the results are very similar.
This result can be attributed to the fact that correlation-pooling takes into
account the relationships between features when generating new features, which
contributes to its relatively better performance.

Table 2: Comparing FeatGeNN performance with Correlation-Pooling and Max-Pooling. The * denotes the version of the FeatGeNN that was executed with Correlation-Pooling. The results are the average score, and the standard deviation, after 30 runs Dataset | Base | FeatGeNN | FeatGeNN*
---|---|---|---
SpamBase | 0.9102 | 0.9422 (0.011) |  0.9530 (0.016)
Megawatt1 | 0.8890 | 0.9148 (0.002) | 0.9151 (0.002)
Ionosphere | 0.9233 | 0.9587 (0.012) |  0.9667 (0.004)
SpectF | 0.7750 | 0.8682 (0.018) |  0.8776 (0.013)
Credit_Default | 0.8037 | 0.8092 (0.003) | 0.8095 (0.003)
German Credit | 0.7401 | 0.7775 (0.006) |  0.7814 (0.002)

###  4.3 Impact of the Number of Data on the Correlation Pooling Computation
(RQ2)

In this subsection, our experiment aims to answer the question: What is the
influence of the number of available data on the Correlation-Pooling
computation?. Figure 4 shows the performance of three versions of FeatGeNN:
FeatGeNN (using all available data), FeatGeNN (using 60% of the data), and
FeatGeNN* (using 30% of the data).

The results show that, as expected, the performance of the model varies with
the amount of data used to compute the correlation-pooling. On average, the
version with access to the entire dataset achieves a performance improvement
of 0.76% and 1.38% compared to the FeatGeNN and FeatGeNN* versions,
respectively. Compared to the version that used 80% of the available data, the
result after 30 epochs is very similar, although the version with more data
performs better in fewer epochs. These results indicate that the performance
of FeatGeNN is still competitive with the original version, even though the
performance decreases slightly with less available data.

Figure 4: The performance of the different versions of FeatGeNN is compared in
terms of the amount of data used for computation. In the image, the symbol
represents the version that used 60% of the available data, the * symbol
represents the version that used 30% of the data, and the symbol ° represents
the version that used 100% of the data. The FeatGeNN without symbol stands for
the version that used 80% of the available data.

###  4.4 Effectiveness of FeatGeNN (RQ3)

In this subsection, this experiment aims to answer: Can our FeatGeNN with
Correlation-Pooling achieve competitive results when compared to the state-of-
the-art models?. We compare FeatGeNN on 6 datasets with state-of-the-art
methods, including (a) Base: Raw dataset without any transformation; (b)
Random: randomly apply a transformation to each raw feature; (c) DFS [4]; (d)
AutoFeat [11]; (e) LFE [6]; (f) NFS [12]; and (g) DIFER [7].

Table 3: Comparison between FeatGeNN with other methods from the literature, reported on [7]. * reports the average and standard deviation across 30 runs, while the FeatGeNN column reports the maximum value across the same runs. Dataset | Base | Random | DFS | AutoFeat | NFS | DIFER | FeatGeNN* | FeatGeNN
---|---|---|---|---|---|---|---|---
SpamBase | 0.9102 | 0.9237 | 0.9102 | 0.9237 | 0.9296 | 0.9339 | 0.9530 (0.016) | 0.9644
Megawatt1 | 0.8890 | 0.8973 | 0.8773 | 0.8893 | 0.9130 | 0.9171 | 0.9151 (0.002) | 0.9171
Ionosphere | 0.9233 | 0.9344 | 0.9175 | 0.9117 | 0.9516 | 0.9770 | 0.9644 (0.012) | 0.9713
SpectF | 0.7750 | 0.8277 | 0.7906 | 0.8161 | 0.8501 | 0.8612 | 0.8776 (0.013) | 0.8802
Credit_Default | 0.8037 | 0.8060 | 0.8059 | 0.8060 | 0.8049 | 0.8096 | 0.8095 (0.003) | 0.8102
German Credit | 0.7410 | 0.7550 | 0.7490 | 0.7600 | 0.7818 | 0.7770 | 0.7814 (0.002) | 0.7827

Table 3 shows the comparative results of FeatGeNN relative to existing methods
(results reported in [7]). From Table 3 we can observe that in the
classification tasks, the comparison shows that FeatGeNN, performs the best
for the SpamBase, Credit_Default, German Credit, and SpectF benchmarks, the
second best for the Ionosphere benchmark and achieves the same result as the
DIFER method for the Megawatt1 benchmark. Although DIFER achieves the best
performance in the Ionosphere benchmark, they only achieve 0.58% more than the
best result obtained by our proposed method.

Regarding the number of features, Table 4 shows that FeatGeNN excels in
producing fewer features for the Megawatt1, SpectF, and Credit_Default
datasets compared to other methods. For the remaining datasets, FeatGeNN
achieves comparable results with the same number of features.

Compared to the performances of Base and Random, FeatGeNN achieved an average
improvement of 5.89% considering all datasets, which demonstrates the
potential of the features generated by our proposed model.

Table 4: Comparison between FeatGeNN, DIFER, AutoFeat, and Random ($\ast$ the results reported on [7]). Dataset | Random | AutoFeat* | NFS* | DIFER* | FeatGeNN
---|---|---|---|---|---
SpamBase | 1 | 46 | 57 | 1 | 1
Megawatt1 | 8 | 48 | 37 | 29 | 8
Ionosphere | 1 | 52 | 34 | 1 | 1
SpectF | 8 | 37 | 44 | 9 | 8
Credit_Default | 4 | 30 | 25 | 5 | 4
German Credit | 1 | 22 | 24 | 1 | 1

##  5 Conclusion

In this study, we presented a novel approach for generating new features in
tabular data that combines feature selection and feature generation to improve
the performance of predictive models. Our proposed method uses a CNN
architecture to effectively capture local features during convolution
operations (Local Feature Extraction), thereby reducing the number of
combinations required in the MLP phase (Global Feature Generation). In
addition, we integrated a correlation-pooling operation as a dimensionality
reduction step. Our approach demonstrates efficient feature learning and
achieves competitive results compared to the architecture used by Max-Pooling
and state-of-the-art methods.

As a direction for future research, we intend to explore information theory
methods as possible alternatives for pooling operations. This could further
increase the effectiveness of our approach to learning new features.

## References

  * [1] Domingos, Pedro M.. “A few useful things to know about machine learning.” Communications of the ACM 55 (2012): 78 - 87.
  * [2] Guo, H., Tang, R., Ye, Y., Li, Z., and He, X.. ”DeepFM: A Factorization-Machine based Neural Network for CTR Prediction.” 2017 arXiv. https://doi.org/10.48550/arXiv.1703.04247
  * [3] Cheng, H., Koc, L., Harmsen, J., Shaked, T., Chandra, T., Aradhye, H., Anderson, G., Corrado, G., Chai, W., Ispir, M., Anil, R., Haque, Z., Hong, L., Jain, V., Liu, X., and Shah, H.. ”Wide & Deep Learning for Recommender Systems.”. 2016 arXiv. https://doi.org/10.48550/arXiv.1606.07792
  * [4] Kanter, James Max and Kalyan Veeramachaneni. “Deep feature synthesis: Towards automating data science endeavors.” 2015 IEEE International Conference on Data Science and Advanced Analytics (DSAA) (2015): 1-10.
  * [5] Udayan, Khurana, Turaga Deepak, Samulowitz Horst and Parthasrathy Srinivasan. “Cognito: Automated Feature Engineering for Supervised Learning.” (2016).
  * [6] Nargesian, Fatemeh, Horst Samulowitz, Udayan Khurana, Elias Boutros Khalil and Deepak S. Turaga. “Learning Feature Engineering for Classification.” IJCAI (2017).
  * [7] Zhu, Guanghui, Zhuoer Xu, Xu Guo, C. Yuan and Yihua Huang. “DIFER: Differentiable Automated Feature Engineering.” ArXiv abs/2010.08784 (2020): n. pag.
  * [8] Gábor J. Székely. Maria L. Rizzo. Nail K. Bakirov. ”Measuring and testing dependence by correlation of distances.” Ann. Statist. 35 (6) 2769 - 2794, December 2007. https://doi.org/10.1214/009053607000000505
  * [9] Liu, Bin, Ruiming Tang, Yingzhi Chen, Jinkai Yu, Huifeng Guo and Yuzhou Zhang. “Feature Generation by Convolutional Neural Network for Click-Through Rate Prediction.” The World Wide Web Conference (2019): n. pag.
  * [10] Cheng, Weiwei, Gjergji Kasneci, Thore Graepel, David H. Stern and Ralf Herbrich. “Automated feature generation from structured knowledge.” CIKM ’11 (2011).
  * [11] Horn, F., Robert T. Pack and Michael Rieger. “The autofeat Python Library for Automatic Feature Engineering and Selection.” ArXiv abs/1901.07329 (2019): n. pag.
  * [12] Chen, Xiangning, Bo Qiao, Weiyi Zhang, Wei Wu, Murali Chintalapati, Dongmei Zhang, Qingwei Lin, Chuan Luo, Xudong Li, Hongyu Zhang, Yong Xu, Yingnong Dang, Kaixin Sui and Xu Zhang. “Neural Feature Search: A Neural Architecture for Automated Feature Engineering.” 2019 IEEE International Conference on Data Mining (ICDM) (2019): 71-80.
  * [13] Zoph, Barret and Quoc V. Le. “Neural Architecture Search with Reinforcement Learning.” ArXiv abs/1611.01578 (2017): n. pag.
  * [14] Atzmüller, Martin and Eric Sternberg. “Mixed-Initiative Feature Engineering Using Knowledge Graphs.” Proceedings of the Knowledge Capture Conference (2017): n. pag.
  * [15] Lam, Hoang Thanh, Johann-Michael Thiebaut, Mathieu Sinn, Bei Chen, Tiep Mai and Oznur Kirmemis Alkan. “One button machine for automating feature engineering in relational databases.” ArXiv abs/1706.00327 (2017): n. pag.
  * [16] Peng, Hanchuan, Fuhui Long and C. Ding. “Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy.” IEEE Transactions on Pattern Analysis and Machine Intelligence 27 (2005): 1226-1238.
  * [17] Katoch, Sourabh, Sumit Singh Chauhan and Vijay Kumar. “A review on genetic algorithm: past, present, and future.” Multimedia Tools and Applications 80 (2021): 8091 - 8126.
  * [18] Bahdanau, Dzmitry, Kyunghyun Cho and Yoshua Bengio. “Neural Machine Translation by Jointly Learning to Align and Translate.” CoRR abs/1409.0473 (2015): n. pag.
  * [19] Iandola, Forrest N., Matthew W. Moskewicz, Khalid Ashraf, Song Han, William J. Dally and Kurt Keutzer. “SqueezeNet: AlexNet-level accuracy with 50x fewer parameters and ¡1MB model size.” ArXiv abs/1602.07360 (2016): n. pag.
  * [20] Cheng, Heng-Tze, Levent Koc, Jeremiah Harmsen, Tal Shaked, Tushar Chandra, Hrishikesh B. Aradhye, Glen Anderson, Gregory S. Corrado, Wei Chai, Mustafa Ispir, Rohan Anil, Zakaria Haque, Lichan Hong, Vihan Jain, Xiaobing Liu and Hemal Shah. “Wide & Deep Learning for Recommender Systems.” Proceedings of the 1st Workshop on Deep Learning for Recommender Systems (2016): n. pag.
  * [21] Shcherbakov, Maxim V., Adriaan Brebels, Anton Tyukov, Timur Janovsky and Valeriy Anatol. “A Survey of Forecast Error Measures.” (2013).
  * [22] Hosmer, David W. and Stanley Lemeshow. “Applied Logistic Regression.” (1989).
  * [23] Cortes, Corinna and Vladimir Naumovich Vapnik. “Support-Vector Networks.” Machine Learning 20 (2004): 273-297.
  * [24] Chen, Tianqi and Carlos Guestrin. “XGBoost: A Scalable Tree Boosting System.” Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (2016): n. pag.
  * [25] Tibshirani, Robert. “Regression Shrinkage and Selection via the Lasso.” Journal of the royal statistical society series b-methodological 58 (1996): 267-288.
  * [26] Smola, Alex and Bernhard Schölkopf. “A tutorial on support vector regression.” Statistics and Computing 14 (2004): 199-222.
  * [27] A. Kaul, S. Maheshwary and V. Pudi, ”AutoLearn — Automated Feature Generation and Selection,” 2017 IEEE International Conference on Data Mining (ICDM), 2017, pp. 217-226, doi: 10.1109/ICDM.2017.31.
  * [28] D. E. Goldberg and K. Deb, ”A comparative analysis of selection schemes used in genetic algorithms,” in Foundations of Genetic Algorithms, G. J. E. Rawlins, Ed. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 1991, pp. 69-93.
  * [29] Pearson, K. (1895). ”Note on regression and inheritance in the case of two parents.” Proceedings of the Royal Society of London, 58: 240-242.