# Time-Resolved Open-Circuit Conductive Atomic Force Microscopy for
Quantitative Analysis of Nanowire Piezoelectricity and Triboelectricity

Yonatan Calahorra,1∗ Wonjong Kim,2 Jelena Vukajlovic-Plestina,2
Anna Fontcuberta i Morral,23 Sohini Kar-Narayan1

1Department of Materials Science and Metallurgy
University of Cambridge, CB3 0FS, Cambrdige, UK
2Laboratory of Semiconductor Materials, Institute of Materials, School of
Engineering
Ecole polytechnique fédérale de Lausanne (EPFL), 1015 Lausanne, Switzerland
3Institute of Physics, School of Basic Sciences
Ecole polytechnique fédérale de Lausanne (EPFL), 1015 Lausanne, Switzerland

∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
> Piezoelectric nanowires are promising materials for sensing, actuation and
> energy harvesting, due to their enhanced properties at the nanoscale.
> However, quantitative characterization of piezoelectricity in nanomaterials
> is challenging due to practical limitations and the onset of additional
> electromechanical phenomena, such as the triboelectric and piezotronic
> effects. Here, we present an open-circuit conductive atomic force microscopy
> (cAFM) methodology for quantitative extraction of the direct axial
> piezoelectric coefficients of nanowires. We show, both theoretically and
> experimentally, that the standard short-circuit cAFM mode is inadequate for
> piezoelectric characterization of nanowires, and that such measurements are
> governed by competing mechanisms. We introduce an alternative open-circuit
> configuration, and employ time-resolved electromechanical measurements, to
> distinguish between electrical generation mechanisms and extract the
> piezoelectric coefficients. This method was applied to nanowires of GaAs, an
> important semiconductor, with relatively low piezoelectric coefficients. The
> results obtained for GaAs piezoelectric coefficient, $\sim$0.4-1 pm/V, are
> in good agreement with existing knowledge and theory. Our method represents
> a significant advance in understanding the coexistence of different
> electromechanical effects, and in quantitative piezoelectric nanoscale
> characterization. The easy implementation will enable better understanding
> of electromechanics at the nanoscale.

## Introduction

The several decades long interest in semiconductor nanowires (?) (NWs) has
brought focus to a topical niche - that of piezoelectric semiconductor NWs,
pioneered by Wang et al. (?, ?, ?), with potential applications in sensing,
energy harvesting and logic (?, ?, ?). Three distinct electromechanical
effects are manifested strongly in semiconductor nanowires: the high aspect
ratio allows large elastic deformations, therefore enhancing the piezoelectric
effect (?), describing changes in surface polarisation due to applied strain;
while the increased surface-to-volume ratio enhances interfacial effects such
as triboelectricity (?, ?), relating to surface charge transfer upon contact
with a dissimilar material; in addition, the combination of semiconducting and
piezoelectric properties results in a unique electromechanical phenomenon
known as the piezotronic effect (?, ?), whereby the height of a semiconductor
energy barrier for charge carrier transport is changed due to mechanical
pressure.
When considering electromechanical current/voltage generation from
semiconductor NWs, both single NWs and NW ensembles or arrays have been
considered (?, ?, ?, ?). In particular, an abundance of conductive atomic
force microscopy (cAFM) measurements have been reported on piezoelectric
semiconductor NWs (mostly ZnO and III-N). In cAFM, the tip is scanned along
the surface, while the current is recorded simultaneously, usually under
applied bias. This mode of operation is also useful for other physical
mechanism characterization, such as photovoltaics (?). However,
straightforward application of cAFM for electromechanical characterization of
piezoelectric NWs is challenging. The two common variants for piezoelectric NW
characterization are the short-circuit cAFM mode: recording the current
generated by deforming a NW with an AFM tip (see Fig. 1) (?, ?, ?, ?); and the
loaded configuration (resiscope mode): recording the voltage developed across
a resistor in parallel to the deformed NW (?, ?, ?, ?, ?).
Although this is a widely studied topic and a commonly conducted experiment,
very little attention has been given to the complete set of mechanically
induced current flow mechanisms, namely piezoelectric, triboelectric and
piezotronic. The combination of these effects is to be expected due to the
nature of measurement involving dynamic forces and contact characteristics
experienced by the NW-tip system, as the tip is scanned in and out of contact
with NWs. Indeed, we are aware of only two reports attempting a quantitative
analysis of the measured voltage or current (?, ?), to extract or compare the
piezoelectric coefficients with theory - albeit with limited success.
The contrast between abundant experimental demonstrations of electrical
generation from piezoelectric NWs and lack of quantitative discussion, is
closely related to an uncertainty in the physical origins of the measured
signals. In particular, the initial report on cAFM based generation from ZnO
NW (?), was followed by alternative analysis of the results (?, ?), and a
rebuttal from Wang (?). Furthermore, flat ferroelectric samples have been
studied by similar approaches (?, ?, ?). In particular, a direct piezoresponse
force microscopy (PFM) method, introduced by Gomez et al. (?), was used to
quantitatively extract direct piezoelectric coefficients; however, the nature
of that method precludes its application to nanomaterials, or materials of
lower piezoelectric coefficients, as discussed below.
Herein, we report a new methodology to perform this experiment, enabling
reliable extraction of the piezoelectric coefficients using cAFM. The
prevailing method for quantitative analysis of piezoelectricity is PFM, which
is generally not suitable for non-planar samples (?), and characterises
converse piezoelectricity (electrical to mechanical). We performed a detailed
set of cAFM experiments at different conditions to isolate and distinguish the
various current mechanisms in vertical pressure based cAFM measurements. Using
scanning and ramp-mode cAFM in conjunction, we show that all three current
mechanisms highlighted above are present in the case of GaAs NWs. We also show
that the piezotronic and triboelectric currents could dominate the
measurements to the point of obscuring piezoelectric data. We distinguish
between short- and open-circuit configurations, and their practical
applications in relation to cAFM based measurements. We suggest an open-
circuit configuration and a time-resolved measurement methodology to extract
the piezoelectric coefficients demonstrating good agreement to known values of
GaAs. Obtaining these numbers for GaAs, a weak piezoelectric material,
indicates promise to use with other piezoelectric materials. Our method offers
a route for quantitative analysis of the direct piezoelectric effect. It is
easily implemented, and can be extended to other piezoelectric nanostructures,
as well as to characterize triboelectricity.

### Piezoelectricity, triboelectricity and piezotronics

Research into nanoscale electromechanical phenomena has evolved from mere
piezoelectricity to include a wide body of work related to the unique
combination of piezoelectricity and semiconducting properties, coined as the
piezotronic effect (?, ?), as well as to triboelectricity or contact
electricfication (?). The work presented here aims to distinguish these
effects in cAFM.

### Piezoelectricity

Piezoelectricity, first described by Jacques and Pierre Curie, is the linear
interrelation between electrical polarization and applied stress in non-
centrosymmetric materials (?). There are four types of piezoelectric
coefficients, correlating the electric displacement D, electric field E,
stress field T and strain field S (?)

$\displaystyle d_{ij}$ $\displaystyle=\left(\frac{\partial D_{i}}{\partial
T_{j}}\right)^{E}=\left(\frac{\partial S_{j}}{\partial E_{i}}\right)^{T}$ (1)
$\displaystyle e_{ij}$ $\displaystyle=\left(\frac{\partial D_{i}}{\partial
S_{j}}\right)^{E}=-\left(\frac{\partial T_{j}}{\partial E_{i}}\right)^{S}$ (2)
$\displaystyle g_{ij}$ $\displaystyle=-\left(\frac{\partial E_{i}}{\partial
T_{j}}\right)^{D}=\left(\frac{\partial S_{j}}{\partial D_{i}}\right)^{T}$ (3)
$\displaystyle h_{ij}$ $\displaystyle=-\left(\frac{\partial E_{i}}{\partial
S_{j}}\right)^{D}=-\left(\frac{\partial T_{j}}{\partial D_{i}}\right)^{S}$ (4)

indices i,j follow Voigt notation, e.g. $d_{ij}$ describes the $i$ component
of polarization in response to stress applied in the $j$ direction, or
alternatively strain in the $j$ direction in response to an electric field
along the $i$ direction. The two alternative descriptions describe the direct
(left-hand side equality) and converse (right-hand side) piezoelectric effect.
The superscripts indicate differentiation under constant or zero fields. This
has a bearing on the application of the relation to experimental conditions,
depending on the boundary conditions. For example, a clamped piezoelectric
structure operates under zero strain condition and will develop stress
following the application of an electric signal. Alternatively, a structure
clamped on one side operates under zero or constant stress. This is the
operational condition of PFM, where an electric field is applied to the
sample, and the strain is monitored - a direct measurement of the out of plane
converse piezoelectric strain coefficient $d_{jj}$ (alternatively named the
piezoelectric charge coefficient when considering the direct effect). The
application of these equations to the direct piezoelectric effect will be
discussed below in detail.
The coefficients are superpositioned with the linear mechanical or electrical
equations to yield the complete linear electromechanical relations, for
example

$\displaystyle D$ $\displaystyle=\epsilon_{T}E+dT$ (5) $\displaystyle S$
$\displaystyle=d^{t}E+s_{E}T$ (6)

where $\epsilon$ is the dielectric permittivity of the material, subscript $T$
represents constant stress, $d$ is the piezoelectric charge tensor, and
$d^{t}$ is the transpose of $d$ (the interchangeably used piezoelectric strain
tensor). $s_{E}$ is the elastic compliance at constant electric field.
Notably, the piezoelectric charge coefficient, $d$, is used in short-circuit
configuration (constant or zero electric field), where the piezoelectric
element is treated as a current source. In an open-circuit configuration the
piezoelectric voltage coefficient, $g$, is used - where the piezoelectric
element is treated as a voltage source.
Piezoelectricity has many practical applications, including in sensing and
actuation, however, there is a common thread to them: piezoelectricity is a
decaying phenomenon. As such, only changes in stress yield useful electrical
signals, and the main avenue of applications is in the AC regime.

### Triboelectricity

Triboelectricity or contact electrification is a result of charge transfer
between materials with different work-functions, where charge is transferred
to equate the electrochemical potential upon contact (?). Triboelectricity is
not limited to semiconductors or metals, or crystalline materials, and is a
property of soft and biomaterials as well (?). The interest in nanoscale
triboelectricity has peaked in the past decade due to the enhanced surface
area of nanomaterials, significantly increasing the efficiency of
triboelectric power generation (?, ?, ?). Similar to piezoelectricity,
triboelectricity is also instantaneous. It occurs upon contact between the two
materials, and decays with the interfacial electrochemical potential reaching
equilibrium.

### Piezotronics

The research into semiconductor nanowire piezoelectricity has brought into
focus the combination of semiconducting and piezoelectric properties in a
single material - the piezotronic effect. (?, ?, ?) Briefly, the main
characteristic of a semiconductor is tunable conductivity through application
of an electric field, as well as through physical contact with a different
material. For example, in a pn junction or a metal-semiconductor contact, a
depletion region in the semiconductor is formed in proximity to the interface.
When it comes to piezoelectricity, depletion, or lower conductivity, is
associated with higher piezoelectric responses. The piezotronic effect
therefore describes the effect strain has on a semiconductor device, in
particular the change of an energy barrier height due to changes in the
interface charge, induced by the application of stress. (?) The adequate
physical treatment of this phenomena is not significantly different than that
of non-ideal metal-semiconductor contacts having interface states within the
band-gap, as suggested by Bardeen (?, ?). The piezotronic effect takes the
form

$\Delta\phi_{B}=\phi_{B,PT}-\phi_{B,0}\sim\ P_{piezo}$ (7)

where the difference between the unperturbed barrier height ($\phi_{B,0}$) and
the strained barrier height ($\phi_{B,PT}$) is proportional to the
polarization ($P_{piezo}$).
There are two important distinctions between the piezotronic effect and the
piezoelectric effect, from which it originates: i) the piezotronic effect does
not decay, and the barrier height change (Eq. 7) holds as long as the strain
is held; ii) the piezotronic effect results in an exponential relation between
the device current and the applied strain/stress, unlike piezoelectricity
which is a linear effect by nature. This is due to the exponential relation
between the current and the barrier height. This fact presents tremendous
promise for piezotronic pressure and strain sensors, as we have recently
demonstrated using GaAs NWs similar to the ones examined in this work (?).

### Electromechanical cAFM measurements of semiconductor NWs

The prevalent cAFM measurement of vertically aligned (as-grown) NWs is
bending-induced current generation. This is due to two reasons, related to the
geometry of the AFM tip/NW system: i) ideally, contact mode piezoresponse
force microscopy (PFM (?, ?)) would have been desirable to measure the
piezoelectric properties of the NW, and specifically the axial piezoelectric
coefficient. However, NWs are brittle when considering contact mode operation,
and therefore such reports are scarce, and require mechanical constraints (?).
We have previously addressed this issue by developing a non-destructive PFM
operation mode (?) and applying it to horizontal and vertical III-V NWs (?).
This limitation, together with the need to characterise the direct
piezoelectric effect (rather than the converse effect in the case of PFM), has
directed work towards non-contact mode electromechanical cAFM; ii) while
slender NWs are easily bent, their axial stiffness is considerably larger. The
bending (cantilever) and compression (axial) force constants of a typical NW
are given by

$\displaystyle
k_{bending}=\frac{3YI}{L^{3}}\longrightarrow\frac{Ya^{4}}{4L^{3}}$ (8)
$\displaystyle k_{axial}=\frac{YA}{L}\longrightarrow\frac{Ya^{2}}{L}$ (9)

where $Y$ is Young’s modulus, $I$ is the second moment of inertia, $L$ the NW
length, $A$ is NW cross section area, and $a$ is the edge of a hypothetical
square cross-sectioned NW. For a moderate aspect ratio square cross-section
NW, 1 $\mu$m long and 100 nm wide, with a Young’s modulus of 100 GPa, we get
$k_{axial}=1000$ N/m, while $k_{bending}=2.5$ N/m. It is therefore clear that
bending NWs is generally much easier than compressing them, and hence the
common cAFM NW characterization methodology measures current during bending
while the tip is scanned across the NWs.
The common description for carrying out NW deformation-induced current
generation measurements is as follows (?, ?): as the tip is scanned across the
sample and a NW is bent, the compressed and stretched sides of the NW develop
a potential difference. A current flows when the tip is in contact with one of
these sides. Interestingly, the common observation is that only one of the
voltage polarities is associated with measurable current flow. This is
attributed to the rectifying electrical properties of the tip-NW contact (?).
Notably, this is in contrast to the common observation in NW ensemble device
measurements, where two current/voltage peaks of opposite sign are associated
with straining and relaxing the device (?, ?). Moreover, the role of the
rectifying contact is deemed even more fundamental. It was found that the
existence of a rectifying contact is critical for electrical generation by the
strain/compression of the piezoelectric semiconductor nanowire, to prevent a
non-efficient current route (i.e. involving power dissipation within the NW)
developing between the polarised regions of the NW (?, ?, ?, ?, ?). The
underlying reasoning could be that NWs with good electrical contacts are in
fact highly doped, and their piezoelectric properties deteriorate. These
points have been a part of the discussion regarding the origins of the
measured signal (?, ?).

## Results and discussion

### Pressure-dependent contact current

Figure 1: (a) Measurement schematic, showing the AFM cantilever descending
atop the NW, and the current measured; (b) The sample used in these
measurements viewed by SEM in an angle of 20∘. Scale bar is 1 $\mu$m.

Figure 1a shows a schematic of the experimental procedure, where the AFM tip
is scanned in PeakForce mode (?) along the surface, while recording the
resulting current, referred to as the PF-TUNA mode. Figure 1b shows a scanning
electron microscope (SEM) image of the examined NW sample. Nanowire growth
time was controlled to maintain a relatively short length, rendering the NWs
more stable under AFM scanning. Noticeably, alongside the NWs, parasitic
growth dominates the silicon surface, on which the NWs were grown by molecular
beam epitaxy (MBE).

Figure 2: Scanning peak-force influence on measured short-circuit current map
(a) Single NW height channel; (b,c,d) contact current maps obtained with 22.5,
45, 62.5 nN peak-force. The active area increased with applied force.

Figure 2 shows a set of PF-TUNA current maps obtained from a NW (topography in
black and white), with increasing peak-force (the maximum mechanical load, see
Methods Section) values of 20-60 nN. The electrical response area increased
with the mechanical stimulus, however the highest current was measured at 45
nN peak force. This indicates that the relation between the measured current
and applied force is not straightforward. Furthermore, the current was not
correlated to the top of the NW, but to a region along its side, found to be
consistent throughout the measurements. This result could be explained with an
increased deformation of this side, however the measured current was found to
be uncorrelated with the deformation maps (Supporting Information Figure S1),
and this result recurred throughout this work. Moreover, this observation was
reinforced through measurements obtained from the sample in general: zero-bias
current is observed throughout the sample, upon contact with the surface, and
the parasitic GaAs growth (see Figure S2). Bearing in mind that the parasitic
growth is not expected to produce a significant piezoelectric signal, and is
irregularly deformed, we conclude that the current observed is not
piezoelectric per se. The following set of experiments explore the two
additional electromechanical current mechanisms.

### Bias-dependent contact current

To further explore the nature of the measured current, a small bias was
applied to the sample during scanning. Figure 3 shows current maps obtained
with lower/higher peak force setting (top and bottom) and three biasing
conditions: negative, zero and positive. The bias used was $\pm 10$ mV. The
middle column of Fig. 3, is similar to Fig. 2, where increased peak force
results in broader current “hot-spots”. However, the bias sign was found to
dictate the measured current sign. This is another indication that the
measured hot-spot current is not piezoelectric. The increase in current hot-
spot area with pressure could indicate that the dominant mechanism is
improvement of the contact, and not piezoelectricity. We note that, following
this experiment, we have found that there is an electrical ground imbalance in
the instrument of about 3.5 mV, which probably contributes to the current
measured in zero bias. I-V measurements taken from the current hot-spots and
the low-current areas on top of the NW show a complementary picture where the
low-current areas show a rectified I-V curve, while the hot-spots show better
conductivity. This is most likely related to the geometry of the NW and tip,
and a local barrier reduction due to Schottky effect (Supporting Information
Figure S3). In these experiments, there where no conclusive findings relating
the current to applied pressure, i.e. the piezotronic effect.
Furthermore, since GaAs is considered here, there might be an optoelectronic
contribution to the current through carrier generation by the AFM laser, as
demonstrated recently by Alekseev et al. (?). It is therefore of interest to
eliminate steady-state (ohmic/piezotronic) contributions to the current.

Figure 3: Bias influence on measured current map. Contact current maps
obtained with 42 or (d-f) 62 nN peak-force, and sample-tip bias of -10/0/10 mV
(a,d/b,e/c,f). The bias had a significant effect on measured current.

### Open-circuit current measurements

In order to reduce unwanted steady-state current contributions, we introduced
a standard microscope glass slide between the sample and the AFM stage; this
is schematically shown in Fig. 4a. Notably, the existence of a “steady-state”
electrical contact (i.e. ohmic or rectifying contact) is not necessary for the
observation of generated current/voltage in piezoelectric generators (?, ?).
This is obvious when considering the issue outside the context of
piezoelectric semiconductors: the best piezoelectric materials are insulating,
and therefore do not have good electrical contacts. If so, the central role of
the Schottky contacts in cAFM measurements of electromechanically induced
current (as mentioned above) implies that this current is indeed not purely
due to piezoelectricity. Therefore, we set out to examine cAFM-induced current
in an open-circuit configuration.
Figure 4b,c provides complementary results to Figure 3, obtained under open-
circuit conditions. The current was found to be independent of applied bias,
in a striking contrast to the standard configuration. This implies that steady
state currents were dominating the measurements before. Nonetheless, the
measured current was still not negligible, and maintained the hot-spot
characteristics. This result is intriguing: i) if ohmic/piezotronic
contributions are eliminated, then triboelectric and piezoelectric currents
are possible; ii) considering the uncorrelated nature of the deformation maps
and the current maps, the reasonable conclusion is that triboelectricity is
the dominant effect in the generation of current hot spots in the open-circuit
regime. This result is in good agreement with recent reports of the co-
existence of triboelectricity and photovoltaic current generation in a III-V
material probed by cAFM (?). The question therefore remains: is it possible to
record the piezoelectric generation arising from a single NW? (?, ?, ?)

Figure 4: a) Open-circuit cAFM measurement schematic; b,c) consecutive current
maps obtained from the sample while applying -50 [b], and 50 mV [c]. These
images encompass both parasitic growth and nanowires which contribute to the
signal. The isolation rendered the measurement not sensitive to the bias.

### Time-resolved current measurements

Figure 5: Mechanical ramping induced open-circuit measured current: time-
resolved. (a) Topography (b) current map and (c) deformation map of a NW later
used to obtain feedback-controlled time-resolved current from the high (d) and
low (e) current locations; (d,e) a series of time-resolved current
measurements obtained by mechanical ramping on-top the same position, to the
same maximal force, with varying ramping rate. Note the current axis for (e)
and (d) are not identical. low- and high-current positions yield distinct
time-resolved current characteristics, interpreted as piezoelectric and
triboelectric contributions.

In order to further investigate electromechanical current generation, we moved
on to explore the time-resolved characteristics obtained through mechanical
ramping in open-circuit configuration. This has several advantages compared to
the scanning mode: i) the currents in scanning mode are averaged over
significant periods of time (e.g. the entire contact time or the entire peak-
force period). Therefore time-resolved measurements provide additional
information; ii) ramping is less destructive than scanning and larger forces
can be used, allowing stronger indentations. Here we drove the cantilever into
the sample up to a force of ∼ 1300 nN, which is roughly 10 times larger than
the peak-force used during scanning.
Figure 5 shows time-resolved current measurements obtained from a current hot-
spot and from the top of the NW, where generally low currents were found. The
top part of the image shows the scanning results, demonstrating once more the
lack of correlation between deformation and current. Low- and high-current
spots were chosen, and a series of mechanical ramps with a preset maximal
force were executed, while changing the ramping rate from 20 nm/sec (left-hand
side) to 505 nm/sec (right-hand side). The tip was completely out of touch
with the sample after every measurement, as evident from the force curve
plateau.
The two ramping points gave rise to distinct traits: the current hot-spots
show a current peak measured directly upon tip-sample contact, and a
subsequent decay of the current. Conversely, the low-current NW top, gave rise
to a current peak which is correlated with the trend of the applied force -
most evident in faster (shorter) ramps. These findings support the assumption
that triboelectric current is the dominant mechanism explaining the current
hot-spots, considering the instantaneous onset of current and the lack of
correlation with deformation (green arrows in Fig. 5). Possibly, the geometry
of the tip and sample results in spatial locations which favour the onset of
triboelectric current, for example where sharp features result in a locally
enhanced electric field supporting charge transfer.
The current measured on top of the NW, with good correlation to the applied
force (Fig. 5e), is an indication of piezoelectric generation. This follows
the elimination of DC currents and triboelectric currents, as the dominant
mechanisms in this position. However, there are still ambiguities in
understanding this current: firstly, measurement of GaAs under the red AFM
laser induces an inherent contribution from carrier generation, (?) and it is
possible that the contact current here is affected by that mechanism. Indeed,
overlaying the currents in Fig. 5d,e indicates that a possible background
current is present during the measurement, apart from the peaks (see
supporting Figure S4). Secondly, the lack of an opposite current peak is
intriguing: piezoelectricity is linear, and therefore it might be expected
that both positive and negative currents will appear upon compressing and
releasing the NW.

### Extraction of piezoelectric coefficients

Figure 6: Electronic circuits depicting the different configurations for
piezoelectric generation measurement a) the piezoelectric element as a current
or voltage source, with the NW capacitance and resistance connected in the
appropriate manner; b) the loaded resiscope mode, where neither open-circuit
voltage nor short-circuit current are measured, with $R_{L}$ as the load and
the op-amp as a voltage follower; c) the open-circuit configuration, realizing
a voltage differentiator, with the glass slide capacitance $c_{glass}$, and
the feedback resistor and capacitor, $R_{A},\ c_{p}$; d) the short-circuit
configuration.

For extracting the direct piezoelectric coefficient, mirroring (converse) PFM,
it is necessary to quantify this measurement. We examine the different
configurations in their electronic context in Figure 6. The operation of a
generator is influenced by its load, hence the piezoelectric element acts as a
voltage source under open-circuit conditions, and as a current source under
short-circuit conditions (Fig. 6a). Analysis of the loaded configuration is
not straightforward, and indeed, the operation efficiency of piezoelectric
energy harvesters (as any electrical generator) is affected by the load (?, ?)
(Fig. 6b). Inversely, the open-circuit and short-circuit configurations are
straightforward for (theoretical) analysis, and correspond to the
piezoelectric element acting as a voltage and current source correspondingly
(Fig. 6c,d). In particular, our implementation of the open-circuit
configuration, with a capacitor in series to the device, is a voltage
differentiator, where the measured current maintains

$i_{measured}=c_{glass}\frac{dv}{dt}$ (10)

The other components shown in Fig. 6 are the NW internal resistance
($r_{NW}$), which is expected to be high in our case, and capacitance
($c_{NW}$), which is determined by the surface area of a single NW and
therefore very low ($\sim{}10^{-18}$ F). $R_{A}$ is the current to voltage
converter feedback resistor, and $c_{p}$ is the feedback/stray capacitance.
According to AFM manufacturer $R_{A}\simeq 5$ G$\Omega$ and $c_{p}\sim{}$ fF -
accounting for about 10 kHz operation bandwidth 111private communication with
Bruker. The glass capacitor used here is comprised of an AFM sample holder
plate (1.25 cm in diameter), and the AFM stage, with a 1 mm glass slide in
between. This result is $c_{glass}=5.3$ pF, assuming relative permittivity of
5 (?). Therefore $c_{glass}$ is significantly larger than the other capacitors
in the circuit, and comes into play. Furthermore, the time-resolved experiment
does not involve rapid changes to the equivalent circuits, since the tip is in
consistent contact with the NW, therefore there are no significant changes
expected in the contribution of the op-amp current bias to the signal (?).

### Short-circuit configuration

Short-circuit across the NW (3-axis) indicates that (for simplicity,
considering axial fields and displacements alone)

$\displaystyle E_{3}$ $\displaystyle=0$ (11) $\displaystyle D_{3}$
$\displaystyle=P_{3}=d_{33}T_{3}$ (12)

since there are no fields outside the nanowire, the boundary conditions for
the displacement field $D$ hold $D=\pm Q/A$, where $A$ is the effective
electrode area, and $Q$ the free interfacial charge, neutralising the
polarisation charge, brought by the deformation. Ideally, the current measured
in the experiment is the movement of charges to neutralise the polarisation -
a generated current (Fig. 6d).
The dimensions of the NW measured in Fig. 5 are $L_{NW}=210$ nm and
$r_{NW}=40$ nm. We have tried to take the significant uncertainty in the
diameter of the NW (lateral AFM measurement), of about $\pm 20$ nm, into
account. The piezoelectric charge coefficient of GaAs is $d_{33}\simeq
1.5-2.5$ pm/V (pC/N) (?). If so the application of 1300 nN should result in

$Q=d_{33}\cdot F_{app}=2.5\cdot 10^{-12}\cdot 1300\cdot 10^{-9}=3.25\cdot
10^{-18}\ \mathrm{C}$ (13)

Let us assume this force is reached within 0.1 sec, corresponding to the
fastest ramps in Fig. 5. In that case, the current measured will be $Q/\Delta
t=3.25\cdot 10^{-17}$ A. This current is in fact much lower than the noise
level in our system ($\sim$100 fA). Moreover, even if we consider a material
with $d_{33}=100$ pC/N, and a force 100 times stronger, the current will still
be roughly within the noise limit. Indeed, Gomez et al. circumvented this
limitation by using an external current amplifier in their direct-PFM
application (?). They measured relatively strong piezoelectric materials,
having a flat topography.
We therefore conclude that the short circuit AFM configuration, where the
piezoelectric charge coefficient dominates (?), is mostly not suitable for the
measurement of piezoelectric current generation. When considering that it is
in fact a single NW measurement, this result is not surprising, and
corresponds with previous criticism of single NW current generation (?, ?).
This also corroborates to our observation and conclusion above that other
current mechanisms prevail in this configuration, and are probably what has
been measured in previous work utilizing short-circuit cAFM. The piezoelectric
charge (strain) coefficient, $d$, is directly measurable under the application
of constant electric field and monitoring of strain, i.e. PFM mode.

### Open-circuit configuration

We move forward to the open-circuit configuration, where the boundary
conditions hold $D_{boundary}=0$ (?, ?). The relevant coefficient under zero
displacement field is $g_{33}$, which is the ratio of electric field to
applied stress - the piezoelectric voltage constant. As discussed above, in
our configuration the voltage is not measured directly but the current in the
circuit is recorded. This current is linearly related to the voltage time
derivative (Eq. 10).

Figure 7: Integrated current during force application as a function of
cantilever ramping rate (marks), and power law fitting of the data to extract
the rate independent part. The results were obtained from several NWs, as well
as a control experiment measured on the oxide covered surface. Doted lines are
fittings according to Eq. 14. Tip 1 and Tip 2 are both Adama AD-40-AS, with
force constants of 38 and 27 N/m, correspondingly.

To extract the voltage related to the application of the force we integrate
the current. We focus on the charging part of the time-resolved force-current
curves (ramp approach). Figure 7 shows that charge for several NWs, as a
function of cantilever ramping rate, obtained by numerically integrating the
current up to the point where peak-force is reached (Fig. 5). As discussed
earlier, there still might be parasitic mechanisms contributing to the
measurement. Assuming these are linear with the contact time, the
“piezoelectric” component should be mostly dependent upon the force, rather
than the time. The integrated current is therefore fitted to

$Q=a/rate+q_{PE}$ (14)

where the $1/rate$ term is proportional to contact duration, and therefore
$q_{PE}$ is the time independent charge. The results for NW #1 indicate that
for the 1300 nN force, the charge accumulated on the capacitor was 2.68 pC. As
mentioned earlier, $c_{glass}=5.3$ pf. If so, the voltage drop on the
capacitor is ${q_{PE}/c_{glass}}$ which yields 0.505 V. Using the NW geometry
mentioned above, we find the associated electric field is

$E_{piezo}=V_{piezo}/L_{NW}$ (15)

which yields $E_{1300nN}=2.4$ V/m. The associated stress is given by $T=F/A$,
yielding $T_{1300nN}=258$ N/m2. Calculating the piezoelectric voltage
coefficient, $g_{33}$ ($E/T$), yields $g_{33,GaAs}=0.0091$ Vm/N. The charge
and voltage coefficients hold

$d=g\cdot\epsilon_{0}\epsilon_{r}$ (16)

using $\epsilon_{r,GaAs}=13$, (?) we get a value for the charge coefficient
$d_{33,GaAs}=1.05$ pC/N. A similar analysis for two other NWs (NW 3 and the
higher value for NW 2, see Fig. 7) yields $d_{33,GaAs}=0.4-0.55$ pC/N. These
values are in good agreement with known values and theoretical calculations
for GaAs of 1.5-2.5 pC/N. We note that the experimental results are somewhat
lower than theory. A possible reason is that the force could contribute to
bending the NW and not only compressing it, effectively reducing the apparent
mechanical to electrical coupling, resulting on a lower coefficient
calculated. Furthermore, we see the uncertainty in diameter as the main
contributor to the experimental error. The coefficients are proportional to
the diameter square, and hence the associated error in the coefficient is
double that of the diameter, which we asses as 20% - yielding $\pm\ 40\%$
uncertainty in the coefficients.
The time-resolved nature of this analysis determines the limits of the method,
to be maintained for a valid analysis: 1) the ramping rate needs to be fast
enough for resolving the contribution of the piezoelectric charging, and
diminishing the contributions of other mechanisms (Fig. 7); 2) the ramping
rate needs to be slow enough for the system to respond to the changes. In our
case with $C\simeq 5$ pF, and the NW resistance assessed in
M${\Omega}$-G${\Omega}$, we have a time constant of $\tau=RC\simeq 10\ \mu$S -
$10$ mS. This indicates that even our fastest ramps (100-200 mS long) are slow
enough for the piezoelectric NW to charge the capacitor. This also indicates
that the rapid peak-force operation ($\sim$1-2 kHz), might be too fast to
measure the developed voltage, explaining the lack of measured signal in
scanning mode.
We performed several control experiments to validate the measurements
(Supporting Information Figure S5): i) the open-circuit time-resolved
procedure on the native oxide layer atop the Si growth substrate away from the
NW growth area (Fig. 7), as well as on an unrelated sample with a top ITO
electrode; ii) short-circuit time-resolved measurements of the ITO electrode;
iii) short-circuit time-resolved measurements of a different GaAs NW sample,
where the Si substrate is undoped. It was found that the silicon oxide layer
gave rise to an electromechanical signal, which was lower, though comparable,
to the integrated current during NW measurements. Although native silicon
oxide is not piezoelectric, it was recently demonstrated that the AFM
electromechanical apparatus is sensitive to flexoelectricity - an electrical
response brought by internal strain gradient and vice-versa. Abdollahi et al.
have shown that AFM tips induce non-uniform electrical fields to layers during
PFM operation, which results in the onset of the converse flexoelectric effect
- manifested as an apparent piezoelectric signal (?). Furthermore,
flexoelectricity can contribute to tip-enhanced photovoltage generation (?).
It is very likely that in our case, where the direct piezoelectric effect is
examined, flexoelectricity is responsible for the signals obtained from the
oxide dielectric. The sharp tips used in this work (10 nm tip radius)
corroborate with the quantitative analysis presented by Abdollahi and co-
workers. Notably, piezoelectricity and flexoelectricity can co-exist in
piezoelectric materials. In that sense, vertical NWs are advantageous compared
to films, as the deformation of a rod would be more uniform than that of a
layer, when subjected to a vertical load by an indenter of similar dimensions.
Measurements on top of an ITO electrode in the open-circuit configuration
(Supporting Information Figure S5) gave rise to the decaying curves, similar
to those measured around the current hot-spots, attributed to
triboelectricity. In the short-circuit configuration, the results were
entirely different: the current was 10-100 times larger, and in the opposite
polarity. Finally, measurements of a similar GaAs NW grown on a non-conductive
Si substrate gave a flat reading. These controls show that the current
measured in the differentiator configuration is not an artefact of the system.
This method bears great implications for piezoelectric analysis at the
nanoscale. As mentioned above, only a handful of reports attempt the
extraction of piezoelectric coefficients from the measured data. On the one
hand, using the short-circuit configuration, the generated currents are too
low, and using the measurements generated by other mechanisms yield faulty
results. On the other hand, the loaded configuration (resiscope mode), is
neither open- nor short-circuit, making analysis more complex (?). This leaves
our open-circuit configuration as a promising alternative for the analysis of
piezoelectricity in nanoscale structures, providing a applicable route to
decouple various electromechanical effects.

## Summary

To conclude, we provide a comprehensive theoretical and experimental analysis
for AFM-based piezoelectric generation. We have demonstrated for the first
time that three electromechanical current mechanisms, piezotronic,
piezoelectric and triboelectric, are present in the common short-circuit
configuration, and that triboelectricity and piezotronic effects dominate the
measured current. These findings help settle controversy related to this type
of measurement (?, ?, ?), by showing that indeed piezoelectric current
generation is unlikely the effect measured, while other electromechanical
effects are present. We introduced an open-circuit configuration,
circumventing fundamental limitations we found in the commonly used short-
circuit configuration. Through time-resolved measurements of the generated
current (effectively voltage differentiation), we extracted the piezoelectric
voltage coefficient, which is related to the piezoelectric charge coefficient.
This methodology and analysis procedure will allow characterization of
nanoscale piezoelectric materials and advance the efforts of improving
piezoelectric actuators, sensors and energy harvesters. It can be generalized
for examining triboelectric, flexoelectric and piezotronic devices, depending
on the dominant mechanism.

## Materials and Methods

### Nanowire Growth and Morphology

Ga-catalyzed GaAs NWs were grown on doped silicon by molecular beam epitaxy
(MBE), in a process similar to previous reports (?, ?, ?). See Supporting
Information Section S1 for full details.

### AFM characterization

For AFM characterization, the sample back side was sputtered with gold, and
mounted on a circular metallic AFM holder using silver paint. Atomic force
microscopy was carried out using a Bruker Dimension Icon microscope, using the
PF-TUNA mode, through scanning, I-V and I-Z (ramping) operation. In PF-TUNA
scanning mode, the tip is periodically oscillated above the surface such that
it form intermittent contact with the sample, up to a specified ford, i.e. the
peak-force. The tips used were MESP-RC-V2 (coated silicon) by Bruker, and
Adama AD-40-AS (diamond tips, to prevent tip damage; all time-resolved
measurements reported were done using these tips). To realise an open-circuit
configuration an optical microscope glass slide was introduced between the
conductive AFM stage and the sample holder. The general experimental procedure
was: i) locating a NW through a $\sim$10-15 $\mu$m scan, zooming-in, switching
to ramp mode in the desired location, where the current was recorded during
ramping with controlled rate and force.

### Analysis

Growth analysis was done using SEM (Zeiss MERLIN). AFM results were analysed
using NanoScope Analysis software, and MATLAB was used for calculations and
presentation of results.

## References

  * 1. L. Güniat, P. Caroff, A. Fontcuberta i Morral, Vapor phase growth of semiconductor nanowires: Key developments and open questions, Chemical reviews 119, 8958 (2019).
  * 2. M.-H. Zhao, Z.-L. Wang, S. X. Mao, Piezoelectric characterization of individual zinc oxide nanobelt probed by piezoresponse force microscope, Nano Letters 4, 587 (2004).
  * 3. Z. L. Wang, J. Song, Piezoelectric nanogenerators based on zinc oxide nanowire arrays, Science 312, 242 (2006).
  * 4. X. Wang, et al., Piezoelectric field effect transistor and nanoforce sensor based on a single ZnO nanowire, Nano letters 6, 2768 (2006).
  * 5. Z. L. Wang, Piezopotential gated nanowire devices: Piezotronics and piezo-phototronics, Nano Today 5, 540 (2010).
  * 6. Z. L. Wang, Triboelectric nanogenerators as new energy technology for self-powered systems and as active mechanical and chemical sensors, ACS nano 7, 9533 (2013).
  * 7. J. Briscoe, S. Dunn, Piezoelectric nanogenerators–a review of nanostructured piezoelectric energy harvesters, Nano Energy 14, 15 (2015).
  * 8. G. Zhu, et al., Toward large-scale energy harvesting by a nanoparticle-enhanced triboelectric nanogenerator, Nano letters 13, 847 (2013).
  * 9. Y. S. Choi, Q. Jing, A. Datta, C. Boughey, S. Kar-Narayan, A triboelectric generator based on self-poled Nylon-11 nanowires fabricated by gas-flow assisted template wetting, Energy & Environmental Science 10, 2180 (2017).
  * 10. T. Frömling, R. Yu, M. Mintken, R. Adelung, J. Rödel, Piezotronic sensors, MRS Bulletin 43, 941 (2018).
  * 11. C.-H. Wang, et al., Effects of free carriers on piezoelectric nanogenerators and piezotronic devices made of GaN nanowire arrays, Small 10, 4718 (2014).
  * 12. N. Gogneau, et al., Gan nanowires for piezoelectric generators, physica status solidi (RRL)-Rapid Research Letters 8, 414 (2014).
  * 13. N. Gogneau, et al., From single iii-nitride nanowires to piezoelectric generators: New route for powering nomad electronics, Semiconductor Science and Technology 31, 103002 (2016).
  * 14. G. Liu, et al., Nanogenerators based on vertically aligned InN nanowires, Nanoscale 8, 2097 (2016).
  * 15. D. Mikulik, et al., Conductive-probe atomic force microscopy as a characterization tool for nanowire-based solar cells, Nano Energy 41, 566 (2017).
  * 16. W. Su, Y. Chen, C. Hsiao, L. Tu, Generation of electricity in GaN nanorods induced by piezoelectric effect, Applied Physics Letters 90, 063110 (2007).
  * 17. P. A. Alekseev, et al., Piezoelectric current generation in wurtzite GaAs nanowires, physica status solidi (RRL)–Rapid Research Letters 12, 1700358 (2018).
  * 18. X. Wang, et al., Electricity Generation based on One-Dimensional Group-III Nitride Nanomaterials, Advanced materials 22, 2155 (2010).
  * 19. N. Jegenyes, et al., High Piezoelectric Conversion Properties of Axial InGaN/GaN Nanowires, Nanomaterials 8, 367 (2018).
  * 20. M. Alexe, S. Senz, M. A. Schubert, D. Hesse, U. Gösele, Energy harvesting using nanowires?, Advanced Materials 20, 4021 (2008).
  * 21. M. A. Schubert, S. Senz, M. Alexe, D. Hesse, U. Gösele, Finite element method calculations of ZnO nanowires for nanogenerators, Applied physics letters 92, 122904 (2008).
  * 22. Z. L. Wang, Energy Harvesting Using Piezoelectric Nanowires–A Correspondence on “Energy Harvesting Using Nanowires?” by Alexe et al., Advanced Materials 21, 1311 (2009).
  * 23. S. Hong, et al., Charge gradient microscopy, Proceedings of the National Academy of Sciences 111, 6566 (2014).
  * 24. A. Gomez, M. Gich, A. Carretero-Genevrier, T. Puig, X. Obradors, Piezo-generated charge mapping revealed through direct piezoelectric force microscopy, Nature communications 8, 1113 (2017).
  * 25. O. Kwon, et al., Direct Probing of Polarization Charge at Nanoscale Level, Advanced Materials 30, 1703675 (2018).
  * 26. Y. Calahorra, M. Smith, A. Datta, H. Benisty, S. Kar-Narayan, Mapping piezoelectric response in nanomaterials using a dedicated non-destructive scanning probe technique, Nanoscale 9, 19290 (2017).
  * 27. J. Curie, P. Curie, Phénomènes électriques des cristaux hémièdres à faces inclinées, Journal de Physique theorique et appliquee 1, 245 (1882).
  * 28. D. Damjanovic, Ferroelectric, dielectric and piezoelectric properties of ferroelectric thin films and ceramics, Reports on Progress in Physics 61, 1267 (1998).
  * 29. Q. Jing, S. Kar-Narayan, Nanostructured polymer-based piezoelectric and triboelectric materials and devices for energy harvesting applications, Journal of Physics D: Applied Physics 51, 303001 (2018).
  * 30. J. Zhou, et al., Flexible piezotronic strain sensor, Nano letters 8, 3035 (2008).
  * 31. P. Keil, T. Frömling, A. Klein, J. Rödel, N. Novak, Piezotronic effect at schottky barrier of a metal-zno single crystal interface, Journal of Applied Physics 121, 155701 (2017).
  * 32. J. Bardeen, Surface states and rectification at a metal semi-conductor contact, Physical Review 71, 717 (1947).
  * 33. S. M. Sze, K. K. Ng, Physics of semiconductor devices (John wiley & sons, 2006).
  * 34. Y. Calahorra, et al., Highly sensitive piezotronic pressure sensors based on undoped GaAs nanowire ensembles, Journal of Physics D: Applied Physics 52, 294002 (2019).
  * 35. A. Gruverman, S. V. Kalinin, Piezoresponse force microscopy and recent advances in nanoscale studies of ferroelectrics, Journal of materials science 41, 107 (2006).
  * 36. L. Collins, Y. Liu, O. Ovchinnikova, R. Proksch, Quantitative Electromechanical Atomic Force Microscopy, arXiv preprint arXiv:1904.06776 (2019).
  * 37. M. Minary-Jolandan, R. A. Bernal, I. Kuljanishvili, V. Parpoil, H. D. Espinosa, Individual gan nanowires exhibit strong piezoelectricity in 3d, Nano letters 12, 970 (2012).
  * 38. Y. Calahorra, et al., Exploring piezoelectric properties of III–V nanowires using piezo-response force microscopy, Semiconductor Science and Technology 32, 074006 (2017).
  * 39. M. R. Hasan, S.-H. Baek, K. S. Seong, J. H. Kim, I.-K. Park, Hierarchical ZnO nanorods on Si micropillar arrays for performance enhancement of piezoelectric nanogenerators, ACS applied materials & interfaces 7, 5768 (2015).
  * 40. M. Johar, et al., Stable and high piezoelectric output of GaN nanowire-based lead-free piezoelectric nanogenerator by suppression of internal screening, Nanomaterials 8, 437 (2018).
  * 41. R. Yang, Y. Qin, L. Dai, Z. L. Wang, Power generation with laterally packaged piezoelectric fine wires, Nature nanotechnology 4, 34 (2009).
  * 42. B. Pittenger, N. Erina, C. Su, Quantitative Mechanical Property Mapping at the Nanoscale with PeakForce QNM, Bruker.
  * 43. V. A. Sharov, et al., InP/Si Heterostructure for High-Current Hybrid Triboelectric/Photovoltaic Generation, ACS Applied Energy Materials (2019).
  * 44. S. Crossley, S. Kar-Narayan, Energy harvesting performance of piezoelectric ceramic and polymer nanowires, Nanotechnology 26, 344001 (2015).
  * 45. ECCOSTOCK dielectric materials chart.
  * 46. Y. Calahorra, C. Ou, C. Boughey, S. Kar-Narayan, Semiconductors and Semimetals (Elsevier, 2018), vol. 98, pp. 445–478.
  * 47. Y. Li, S. Hu, Z. Liu, L. Chen, Effect of electrical boundary conditions on ferroelectric domain structures in thin films, Applied physics letters 81, 427 (2002).
  * 48. E. A. Eliseev, S. V. Kalinin, A. N. Morozovska, Finite size effects in ferroelectric-semiconductor thin films under open-circuit electric boundary conditions, Journal of Applied Physics 117, 034102 (2015).
  * 49. GaAs basic properties, Ioffe Institute.
  * 50. A. Abdollahi, N. Domingo, I. Arias, G. Catalan, Converse flexoelectricity yields large piezoresponse force microscopy signals in non-piezoelectric materials, Nature communications 10, 1266 (2019).
  * 51. M.-M. Yang, D. J. Kim, M. Alexe, Flexo-photovoltaic effect, Science 360, 904 (2018).
  * 52. C. Colombo, D. Spirkoska, M. Frimmer, G. Abstreiter, A. F. i Morral, Ga-assisted catalyst-free growth mechanism of GaAs nanowires by molecular beam epitaxy, Physical Review B 77, 155326 (2008).
  * 53. F. Matteini, G. Tütüncüoglu, H. Potts, F. Jabeen, A. Fontcuberta i Morral, Wetting of Ga on SiOx and Its Impact on GaAs nanowire growth, Crystal Growth & Design 15, 3105 (2015).
  * 54. F. Matteini, et al., Impact of the Ga droplet wetting, morphology, and pinholes on the orientation of GaAs nanowires, Crystal Growth & Design 16, 5781 (2016).

## Acknowledgments

Y.C. and S.K-N. are grateful for support from ERC Starting Grant (Grant No.
ERC–2014–STG–639526, NANOGEN), as well as Henry Royce Institute - Cambridge
Equipment grant EP/P024947/1 and the Centre of Advanced Materials for
Integrated Energy Systems ”CAM-IES” grant EP/P007767/1. W.K., J.V.P. and
A.F.i.M. thank SNF for funding through the NCCR QSIT and H2020 via project
Indeed and Nanoembrace.

## Supplementary materials

Supplementary Text
Figs. S1 to S5

### Section S1. Experimental: NW growth

GaAs nanowires were obtained by the Ga-assisted growth method on silicon (?,
?, ?). For this, four inch $<$111$>$ p-doped Si wafers with a resistivity of
$<$ 0.03 $\Omega\cdot$cm (doping $10^{18}$ cm${}^{-3})$ were diced into 35
$\times$ 35 mm2 square chips to fit in the MBE sample holder. The native oxide
of the silicon was first removed by dipping the sample in a buffered HF
solution (7:1) for 5 min. Subsequently, a 1.1 nm oxide was grown using 200 W,
200 sccm O2 plasma power in TEPLATM GigaBatch for 30 seconds. The final oxide
thickness was measured by spectroscopic ellipsometry. The prepared Si chips
were introduced into the ultra-high vacuum (UHV) environment MBE machine (DCA
P600) and subsequently annealed at 500 ∘ C for 2 hours in UHV to ensure a
pristine surface free of water and organic molecules. They were then
transferred to the growth chamber where they were de-gassed at 770 ∘ C for 30
min to further remove any possible surface contaminants. The sample was then
brought to the growth temperatue, 634∘ C. Then, Ga was pre-deposited for 10
min at beam equivalent pressure (BEP) of 1.4$\times 10^{-7}$ Torr
(corresponding to a nominal deposition rate of 1.1 $\mathring{A}$/s). After
than, the As4 source was opened for 20 min at BEP of 2$\times 10^{-6}$ Torr.
No dopants were introduced and hence the NWs were nominally undoped. Identical
nanowires were used in the publication by Calahorra et al. (?).

Fig. S1. (a) Deformation map and (b) Contact current map of the NW in Figure 2
in the main text. Demonstrating the lack of correlation between the high
current areas and deformation. The deformation and current spatial
characteristics remained similar throughout the measurements, even following
sample rotation. They seem to be mostly linked to tip geometry.

Figure S2: (a) Topography and (b) Contact current map of typical large-area
view of the sample. Demonstrating that current is not exclusive to NWs.

Fig. S3. I-V spectroscopy obtained from a NW, with initially low-current
characteristics. a) NW current map; b-d) I-V at increasing contact deflection
setpoints 0.2,0.3,0.4 V; e) back to 0.25 V setpoint after the (b-d)
measurements, recreating the higher resistance. These results show an
increasing voltage with pressure, however it seems like contact is improving
rather than a piezotronic effect. Blue and red curves are trace and retrace,
correspondingly.

Fig. S4. Overlaid time-resolved current measurements from main text,
indicating the possible existence of a parasitic current during the
measurement.

Fig. S5. Control time-resolved current recordings: open-circuit of an ITO
electrode (top five panels), with different ramp rates and max. force; short-
circuit of an ITO electrode (bottom left), and a GaAs NW grown on undoped Si
(bottom right). The controls show the richness of the measured current
characteristics on different samples, and that there is measurement artefact –
in particular in the measurement of the silicon oxide.


# Accurate Spectral Collocation Solutions to some Bratu’s Type Boundary Value
Problems.

Călin-Ioan Gheorghiu Tiberiu Popoviciu Institute of Numerical Analysis, Str.
Fantanele 57, Cluj-Napoca, Romania, email<EMAIL_ADDRESS>
###### Abstract

We solve by Chebyshev spectral collocation some genuinely nonlinear Liouville-
Bratu-Gelfand type, 1D and a 2D boundary value problems. The problems are
formulated on the square domain $[-1,1]\times[-1,1]$ and the boundary
condition attached is a homogeneous Dirichlet one. We pay a particular
attention to the bifurcation branch on which a solution is searched and try to
estimate empirically the attraction basin for each bifurcation variety.

The first eigenvector approximating the corresponding the first eigenfunction
of the linear problem is used as an initial guess in solving the nonlinear
algebraic system of Chebyshev collocation to find the “small”solution. For the
same value of the bifurcation parameter we use another initial guess, namely
lowest basis function (1 point approximation), to find the “big”solution. The
Newton-Kantorovich method solves very fast the nonlinear algebraic system in
no more than eight iterations. Beyond being exact, the method is numerically
stable, robust and easy to implement. Actually, the MATLAB code essentially
contains three programming lines. It by far surpasses in simplicity and
accuracy various methods used to solve some well-known problems.

We end up by providing some numerical and graphical outcomes in order to
underline the validity and the effectiveness of our method, i.e., norms of
Newton updates in solving the algebraic systems and the decreasing rate of
Chebyshev coefficients of solution.

Keywords: Bratu’s problem, nonlinear eigenvalue problem, spectral collocation,
accuracy, Chebfun.

MSC2010: 34C23, 35B32, 65N35, 65N25, 65N12, 65H10

##  1 Introduction

The boundary value problem attached to the Liouville-Bratu-Gelfand operators
in two dimensions is classical in mathematics and physics. Nevertheless,
computational methods for highly accurate solutions are still of much current
interest, particularly in applications. Thus for technical point of view, the
exothermic chemical reaction of bulk materials held in a storage is a common
hazard.

Consider the distribution of temperature, $u,$ within such a material, which
is increased by the release of energy from reaction. We have the model (see
for instance [9])

$u_{t}=D\Delta u+\alpha e^{u},\ \ x\in\Omega,\ \ t>0,$ (1)

where $D$ is the heat conductivity and$\alpha$ is a positive constant. The
rate $\alpha e^{u}$ is actually a simplification of the more usual Arrhenius
law. Let us suppose that $u$ is held fixed, say $=u_{b},$ at the boundary of
the domain $\Omega.$ This is the most representative example in the class of
problems we are going to consider.

Various numerical methods, finite elements, finite difference and finite
volumes for space discretization along with various finite difference schemes
to march in time, have been designed in order to solve initial-boundary value
problems attached to the equation (1). We will mainly confine ourselves to one
and two dimensional steady equations coming from (1). We will also solve some
nonlinear problems where the non linearity in (1) is of the form
$f(\varepsilon,u):=f(\varepsilon,e^{u})$ with
$f:{\mathbb{R}}_{+}\times{\mathbb{R}}\rightarrow{\mathbb{R}}$ continuous and
differentiable with respect to both arguments.

Our main aim is to show that the Chebyshev spectral collocation (ChC) along
with Newton-Kantorovich produces highly accurate solutions to these problems
irrespective to the non linearity and corner singularity of the domain. The
method is robust, very efficient and easy implementable. For a given value of
bifurcation parameter we compute both solutions, the “small”and the “big”one.
For 2D case, up to our knowledge this a novelty.

The strategy in question emerges from the work of John Boyd [1]. However,
recently using ChC along with Newton-Kantorovich we have successfully solved a
third-order singular and nonlinear boundary value problem on the half-line
(see [8]). Actually, we have used o Chebfun code in order to implement both
algorithms. This problem possesses an integral invariant and thus we could
guarantee an accuracy tending to the machine precision. For some similar
nonlinear boundary value problems we have used, beyond collocation method,
Galerkin type methods in order to solve them. Thus we quote some of them [3],
[6] and [7].

In a recent paper [12] the authors propose a variational method in order to
solve Bratu’s problem for two dimensions in an adequate approximate space of
bi quadratic spline functions. Their numerical results look correct but the
authors do not provide very convincing arguments in support of their accuracy.

This paper divides naturally into two parts. The first Section 2 is devoted to
the 1D Bratu’s problem. We find the “small”and the “big”solution and study
their stability. At least graphically we show that the Newton-Kantorovich
process is of order larger than two. We also observe that for both solutions
their Chebyshev coefficients decrease exponentially. A simple Chebfun code has
been very efficient. In the second Section 3, in principle, we have done
almost the same for the 2D Bratu’s problem. However, in order to discretize
the 2D Bratu’s equation we have used the classical Chebyshev differentiation
matrices. Moreover, we believe that we are actually showing for the first time
how the Chebyshev coefficients of solutions to boundary or eigenvalues
problems decrease.

Some concluding remarks are also supplied.

##  2 1D Bratu’s problem

Suppose now $\Omega:=[-L,L]$ in one dimension and consider the steady-state
problem

$u_{xx}+\lambda e^{u}=0,\ u\left(-L\right)=u\left(L\right)=0,\ L>0,$ (2)

where $\lambda:=\alpha/D$. The solution of this problem is even, so normally
one solves the half-problem on $[0,1]$ under mixed boundary conditions
$u^{\prime}\left(0\right)=u\left(L\right)=0.$ Gelfand in [5] solved Bratu’s
problem as an initial value problem with conditions
$u^{\prime}\left(0\right)=0,$ $u\left(0\right)=A,$ $A>$ by shooting. He found
a solution of the form

$w\left(x\right)=2\ln\left[\frac{e^{A/2}}{\cosh\left(Bx\right)}\right],\
B:=\left(\lambda e^{A}/2\right).$ (3)

For a fixed $L,$ by enforcing the boundary condition at $x:=L$ in solution (3)
one can find the following transcendental relation between $A$ and $\lambda$

$\frac{e^{A/2}}{cosh\left(\sqrt{\lambda e^{A}/2}L\right)}=1.$ (4)

It provides the bifurcation pattern for problem (2) (see again Gelfand [5] for
other comments on global solutions to Bratu’s problem (2) and their
invariant).

###  2.1 Bifurcation of positive solutions

Figure 1: Bifurcation diagram for 1D Bratu’s problem.

In order to visualize the bifurcation pattern (4) we use the roots method for
class Chebfun2 (see [4] Part II Functions of two variables (Chebfun2)). Thus
we have got Fig. 1 with the critical point (which is defined to be where
$\lambda(A)$ is a maximum) of coordinates

$\left(A_{*},\lambda_{*}\right)=\left(1.187331536443172,0.8786312538512331\right).$

Specifically, solutions are unique for $\lambda\leq 0$ and for a single
positive value $\lambda_{*}$ called the Frank-Kamenetskii critical value; do
not exist for $\lambda>\lambda_{*}$; and two solutions exist for
$0<\lambda<\lambda_{*}.$ For $\lambda_{*}$ Boyd reported in [2] the value
$0.8784576797812903015$ which is fairly closed to that reported above. We have
to mention that for $L:=1/2$ the critical value $\lambda_{*}=3.51360308$ and
for $L:=2$ we have got $\lambda_{*}=0.2196644$ which means this value
decreases as the length of the integration interval increases.

###  2.2 Chebfun solutions to 1D Bratu’s problem and their linear stability.

Figure 2: Solutions $u_{1}$ and $u_{2}$ (left panel) and their representation
in the phase plane (right panels).

For the same value of $\lambda$, i.e., $\lambda:=0.25$ we compute, using a
simple Chebfun code, the solutions corresponding to the two branches of
bifurcation diagram. The first one, denoted by $u_{1},$ results from Chebfun’s
default initial guess (the zero function). This is the “small”solution. The
other, the “big”one denoted by $u_{2},$ results from the alternative initial
guess $6(1-x^{2})$, i.e., the simplest basis function for Galerkin spectral
method. They both are represented in the left panel of Fig. 2 along with their
trajectories in the phase plane (right panels).

Figure 3: Norms of Newton updates in solving for $u_{1}$ and $u_{2}$ (left
panel) and the Chebyshev coefficients of these solutions (right panel).

In order to evaluate the speed of the convergence process we display in the
left panels of Fig. 3 the evolution of Newton’s updates in computing both
solutions. In the right panel, in a log-linear plot, we display the decreasing
rate of the even Chebyshev coefficients of solutions. The odd coefficients are
lower than $O(10^{-15}).$ A brief look to Fig. 3 leads us to the conclusion
that the computation process performed by Chebfun code is much more expensive
for the second solution than for the first one (the number of Newton’s
iterations and the number of Chebyshev coefficients of solutions) is
considerably higher for $u_{2}$ than for $u_{1}$.

In addition, we will show in what follows that the first “small”solution is
stable but the second one “big”is not. Let’s consider now the unsteady problem
corresponding to (2). It reads

$u_{t}=u_{xx}+\lambda e^{u},\ x\in\left(-L,\ L\right),\ u\left(\pm
L\right)=0,$ (5)

where $\lambda$ is now fixed in the interval $(0,\ \lambda_{\ast}).$

Linearizing about $u=u_{i},\ i=1,2$ we have

$v_{t}=v_{xx}+\left(\lambda e^{u_{i}}\right)v,\ x\in\left(-L,\ L\right),\
u\left(\pm L\right)=0,\ i=1,2.$

Now the eigenvalue problems

$-\phi_{xx}-\left(\lambda e^{u_{i}}\right)\phi=\mu\phi,\ x\in\left(-L,\
L\right),\ \phi\left(\pm L\right)=0,\ i=1,2,$

have the eigenvalues

$\mu_{i,n}=\left(\frac{n\pi}{2L}\right)^{2}-\lambda e^{u_{i}},\ i=1,2.$

It easy to verify that for $L:=1,\ \mu_{1,n}>0$ and thus solution $u_{1}$is
stable, but $\mu_{2,n}<0$ which means that solution $u_{2}$ is unstable.

It is important to observe that in the text [14] the authors obtain the
bifurcation diagram (4) by integrating successively the problem (2) using a
Chebfun code.

##  3 2D Bratu’s problem

The two-dimensional Bratu’s problem is an elliptic partial differential
equation with homogeneous Dirichlet boundary conditions. This problem is given
by

$\left\\{\begin{array}[]{c}\Delta u+\lambda e^{u}=0,\ \textrm{in}\
\Omega:=\left(-L,\ L\right)\times\left(-L,\ L\right),\ L>0,\\\ u=0,\
\textrm{on}\ \partial\Omega.\end{array}\right.$ (6)

The problem is a nonlinear eigenvalue one that is commonly used as a test
problem for many numerical methods.

###  3.1 Bifurcation of positive solutions

In his paper [1] John Boyd finds in a very ingenious and at the same time a
simple way an approximate relationship between the eigenparameter $\lambda$
and the maximum norm of solution denoted by $A$. He uses the lowest basis
function that is a polynomial and vanishes at the four sides of the square,
namely

$u\left(x,y\right):=A\left(1-x^{2}\right)\left(1-y^{2}\right),\
A\in\mathbb{R}.$ (7)

Then he uses a type of the weighted residual argument and finds the one-point
analytic approximation

$\lambda\approx 3.2Ae^{-0.64A}.$ (8)

of the bifurcation diagram. He improves this relationship working with a
three-point collocation approximation. In any case, the bifurcation diagram
remains perfectly analogous to that in Fig. 2.1.

#### Corner singularities and symmetries.

In the same paper J. Boyd also observe that the solution to problem (6) has
singularities at the four corners of the domain proportional to

$r^{2}\log\left(r\right),$

where $r$ is the radial coordinate of a local polar coordinate system centered
on a given corner. This is a very weak singularity in the sense that second-
order finite differences still give an error that is asymptotically of order
$h^{2}$, where $h$ is the grid spacing, as noted in Haidvogel and Zang [10].

It is clear that the solution to Bratu’s equation is symmetric about the
origin in both $x$ and $y$. Consequently, only the even-degree Chebyshev
polynomials are needed to construct the appropriate two-dimensional basis
functions. For a given resolution, this double parity reduces the size of the
basis set by a factor of four. Anyway, we will not exploit this symmetry.
However, if we rotate the square through an angle of $90^{o}$ neither the
equation nor the boundary conditions are altered. Thus, in the language of
group theory, the solution is also “invariant under the rotation group
$C_{4}$”which is the group of rotations through any integer multiple of
$90^{o}$.

###  3.2 ChC solutions to 2D Bratu’s problem.

In order to find this solution we have used the ChC method based on the
collocation differentiation matrices from [15].

This procedure starts with the second order Chebyshev differentiation matrix
$\mathbf{D}^{\left(2\right)}$ of order $N.$ In order to implement the
homogeneous Dirichlet boundary conditions we discard the first and the last
rows and columns of this matrix and obtain the square matrix of order $N-2$
denoted by $\widetilde{\mathbf{D}}^{\left(2\right).}$ Now in order to obtain
the differentiation matrices corresponding to the second order partial
derivatives $\partial^{2}u/\partial x^{2}$ and $\partial^{2}u/\partial y^{2}$
we take the following Kronecker products

$\mathbf{D20}=kron(\mathbf{eye}(N-2,N-2),\widetilde{\mathbf{D}}^{\left(2\right)});\mathbf{D02}=kron(\widetilde{\mathbf{D}}^{\left(2\right)},\mathbf{eye}(N-2,N-2)).$

The discrete Laplacian is now simply $\mathbf{\Delta=D20+D02}$, and the MATLAB
code eigs is now used to find the set of eigenvalues and eigenfunctions.
Actually using the code reshape we transform eigenvector $U$ of length $N^{2}$
into the eigensolution matrix $\mathbf{UU}$ of dimension $N\times N.$

When we solve the 2D Bratu’s problem we add to the above Laplacian the matrix
$\lambda exp(\mathbf{U})$ and call the routine MATLAB fsolve.

#### Fist eigenfunction as an initial guess for Newton-Kantorovich algorithm.

The first eigenfunction of the eigenproblem

$\left\\{\begin{array}[]{c}\Delta u+\lambda u=0,\ \textrm{in}\
\Omega:=\left(-L,\ L\right)\times\left(-L,\ L\right),\\\ u=0,\ \textrm{on}\
\partial\Omega.\end{array}\right.$ (9)

will provide the initial guess for the following iterative process designed to
solve the nonlinear algebraic system.

In order to solve the eigenproblem (9), i.e., to discretize the Laplacian
operator, we use ChC. Subsequently we make use of MATLAB code eig. We find at
least the first ten eigenvalues with a precision of order
$O\left(10^{-10}\right).$ It is known that domains with symmetries will often
have eigenvalues with multiplicities greater than one. The authors of [13]
have observed that, for instance for the square of side $\pi,$ the
eigenproblem (9) has eigenvalues

$\lambda_{m,n}=m^{2}+n^{2},\ m,n=1,2,\ldots.$

Whenever $m\neq n$, $\lambda_{m,n}$ will have multiplicity at least two. A
simple shift of the domain $[-1,1]\times[-1,1]$ into $[0,\pi]\times[0,\pi]$
has enabled us to compare our computed eigenvalues with the above
$\lambda_{m,n}$. Thus we prove the accuracy of our computation.

A contour plot which displays the isolines of the matrix solution to problem
(6) is reported in Fig. 4. It confirm the symmetries discussed above.

Figure 4: Izolines of the “small”solution to problem (6) produced by ChC with
$N:=16$ when $\lambda:=0.5.$

In Fig. 5 we display a sixteen square-point approximation to problem (6).
Actually, in order to have a more suggestive image we interpolate the solution
to a finer grid and then plot. We use the MATLAB routine interp2 with the
option ’spline’.

Figure 5: ChC “small”solution to 2D Bratu’s problem for $\lambda=0.5$ when the
order of approximation of ChC is $N:=20.$

The Chebyshev coefficients of the first eigenfunction and of the solution to
problem (2) are displayed in the left respectively right panels of Fig. 6 in a
log-linear plot. They both confirm the assertion from Paragraph 3.1, i.e., the
coefficients of the odd-degree Chebyshev polynomials must vanish. Actually,
they are found in strips of points of order $10^{-15}$ or less. Comparing the
two panels in this figure we also can conclude that the coefficients of the
first eigenfunction fall off exponentially rather than algebraically with $N$
but the non linearity in Bratu’s equation reduces the rate of this fall. These
coefficients are obtained using fast Chebyshev transform (FCT).

Figure 6: Absolute values of the Chebyshev coefficients for the first
eigenfunction (left panel) and the Chebyshev coefficients of solution to
problem (6) (right panel). In both cases the order of approximation of ChC is
$N:=16.$

#### The lowest basis function as an initial guess for Newton-Kantorovich
algorithm.

Let’s consider now the simplest basis function (7) as an initial guess for
Newton-Kantorovich algorithm. We want to approximate a solution situated on
the right branch of the bifurcation diagram.

Figure 7: ChC “big”solution to 2D Bratu’s problem for $\lambda=0.5$ when the
order of approximation of ChC is $N:=16.$

Now the highest point of solution attains

$u_{max}=4.677164395529806e+00$

which is much higher than in the previous case when we have obtained only

$u_{max}=1.865174060688610e-01.$

It is also easy to see that the allure of the “big”solution is much sharper
than in the previous case.

A contour plot which displays the isolines of the matrix solution in this case
is reported in the left panel of Fig. 8. It confirms the symmetries discussed
above.

Figure 8: Izolines of the “big”solution to problem (6) (left panel) and in a
log-linear plot the Chebyshev coefficients of solution (right panel). All
these results are produced by ChC when the order of approximation of ChC is
$N:=16.$

In Table 1 we report some output parameters that fsolve provides when they
solve the nonlinear algebraic system generated by ChC algorithm. Comparing its
two lines we notice that the computation effort in the second case is double
compared to those for the first case. In terms of elapsed time the difference
is not so large.

Table 1: fsolve outcomes for both initial guesses. Initial guess | Nr. iter. | Func-count. | First-order optimality
---|---|---|---
First eigenfunction of (9) | 4 | 985 | 2.73e-09
Lowest basis function (7) | 8 | 1773 | 2.05e-10

###  3.3 ChC solutions to some L-B-G problems.

Let’s consider in this section first the so called Gelfand’s perturbed problem

$\left\\{\begin{array}[]{c}\Delta u+\lambda e^{\left(\frac{u}{1+\varepsilon
u}\right)}=0,\ \textrm{in}\ \Omega:=\left(-L,\ L\right)\times\left(-L,\
L\right),\ 0<\varepsilon<<1,\\\ u=0,\ \textrm{on}\
\partial\Omega.\end{array}\right.$ (10)

We apply to this problem the same ChC technique along with Newton-Kantorovich
and we do not notice major differences from the Bratu’s problem regardless of
the values of parameter $\varepsilon$. The same conclusion holds when we
consider in the equation (1) the non linearity $f(u):=\cosh(u)$ or
$f(u):=\sinh(u).$ For the sake of brevity, we will not present results for
these cases.

##  4 Concluding remarks

All in all, we have solved some 2D nonlinear Liouville-Bratu-Gelfand problems
with exponential accuracy in spite of corner singularities of the domain. The
exponential accuracy is proved by the rate the Chebyshev coefficients of
solutions decay. We have found out both, the “small”solution as well as the
“big”solution. This latter result is rarely reported.

The Chebyshev pseudospectral method along with the Newton-Kantorovich
algorithm is so efficient that it produce remarkably accurate results with a
very modest approximation orders. Actually we have worked with $10<N<20$ and
the iterative algorithm converged in just a few iterations. Also we have to
remark that the Chebfun, used here for 1D problems, as well as MATLAB codes
based on classical ChC are extremely simple, robust and efficient.

## References

  * [1] Boyd, J.P.: An Analytical and Numerical Study of the Two-Dimensional Bratu Equation. J. Sci. Comput. 1, 183–201 (1986)
  * [2] Boyd, J.P.: Chebyshev polynomial expansions for simultaneous approximation of two branches of a function with application to the one-dimensional Bratu equation. Appl. Math. Comput. 143, 189–-200 (2003)
  * [3] Boyd, J.P., Gheorghiu, C.I.: All roots spectral methods: Constraints, floating point arithmetic and root exclusion. Applied Mathematics Letters 67, 28-–32 (2017)
  * [4] Driscoll, T.A., Hale, N., Trefethen, L.N.: Chebfun Guide. 1st Ed. For Chebfun version 5. (2014) https://www.chebfun.org.
  * [5] Gelfand, I.M.: Some problems in the theory of quasi-linear equations. AMS Trans. Ser. 29, 295–381 (1963)
  * [6] Gheorghiu, C.I., Trif, D.: The numerical approximation to positive solution for some reaction diffusion problems. PU. M. A. 11, 243–253 (2000)
  * [7] Gheorghiu, C.I., Trif, D.: Direct and Indirect Approximations to Positive Solution for a Nonlinear Reaction-Diffusion Problem I. Direct (Variational) Rev. Anal. Numér. Théor. Approx.31 61–-69 (2002)
  * [8] Gheorghiu, C.I., A third-order nonlinear BVP on the half-line. https://www.chebfun.org/examples/ode-nonlin/GulfStream.html
  * [9] Grindrod, P.: The Theory and Applications of Reaction-Diffusion Equation, 2nd edn. Patterns and Waves. Clarendon Press, Oxford (1996)
  * [10] Haidvogel, D. B., Zang, T.: The efficient solution of Poisson’s equation in two dimensions via Chebyshev approximation, J. Comp. Phys. 30, 167–180 (1979).
  * [11] Huang, S-Y., Wang, S-H.: Proof of a Conjecture for the One-Dimensional Perturbed Gelfand Problem from Combustion Theory. Arch. Rational Mech. Anal. 222, 769–-825 (2016). https://doi.org/: 10.1007/s00205-016-1011-1
  * [12] Kouibia, A., Pasadas, M., Akhrif, R.: A variational method for solving two-dimensional Bratu’s problem. Numer. Algor. 84, 1589–-1599 (2020). https://doi.org/10.1007/s11075-020-00957-y
  * [13] Kuttler, J. R., Sigillito, V. G.: Eigenvalues of the Laplacian in Two Dimensions. SIAM Rev. 26, 163–193 (1984)
  * [14] Trefethen, L.N., Birkisson, Á, Driscoll, T. A.: Exploring ODEs. SIAM, Philadelphia (2018)
  * [15] Weideman, J.A.C., Reddy, S. C.: A MATLAB Differentiation Matrix Suite. ACM T. Math. Software. 26, 465–519 (2000)