Orthotropic Piezoelectricity in 2D Nanocellulose

The control of electromechanical responses within bonding regions is essential to face frontier challenges in nanotechnologies, such as molecular electronics and biotechnology. Here, we present Iβ-nanocellulose as a potentially new orthotropic 2D piezoelectric crystal. The predicted in-layer piezoelectricity is originated on a sui-generis hydrogen bonds pattern. Upon this fact and by using a combination of ab-initio and ad-hoc models, we introduce a description of electrical profiles along chemical bonds. Such developments lead to obtain a rationale for modelling the extended piezoelectric effect originated within bond scales. The order of magnitude estimated for the 2D Iβ-nanocellulose piezoelectric response, ~pm V−1, ranks this material at the level of currently used piezoelectric energy generators and new artificial 2D designs. Such finding would be crucial for developing alternative materials to drive emerging nanotechnologies.

the dipole model, or quantum mechanics approaches e.g. to estimate the energy and the structure of atomic systems.
We here define a numerical descriptor 3 (labelled as 't a ') in order to study the effect of external fields on the electromechanical response of bonded atoms. This descriptor accounts for the deviations of the E-field along the bond length relative to the behaviour expected for two positive point charges immersed in an uniform electron density region.
We investigate the dependence between the interatomic position (d) and the angle between E and the interatomic axis (Φ) , i.e. Φ= Φ(d). For ideal behaviours, such dependence is represented by a Heaviside transition between 0º and 180º, in the region between the two extreme atom positions, In the expression above 0 a is the Bohr radius and d 0 represents the point where the 0-180º step transition takes place. In our context, such d 0 -point correspond to the location where E changes the orientation relative to the two-atom bond axis. We define the dimensionless descriptor as 0 a t t a = for finite t -valued functions. This reasoning gives us the following expression for a t , We reported a t values after statistical treatment for a set of N points in each of the studied interatomic interactions. The statistical treatment comprises the fitting of such points to modulated Heaviside functions, eq. S2. The final a t values, were obtained as the average over nearly a thousand points and the standard deviations for each bond was lower than 5 10 − .
The behaviour of the E-vector along the bond axis is indicative of the electrical feasibility of the chemical bonds. In the case of Iβ-nanocellulose (NC) system we restricted the representation of ( ) dependence for three representative hydrogen bonds (HB) pairs and three selected covalent bonds (CB). Deviations from the ideal behaviour are addressed to smooth transitions on the electron density between the pair atoms. The implications of such deviations are highlighted in the main text.

Understanding electrical forces in the ubiquitous water dimer
The water dimer offers a feasible and comprehensive frame to study HB 4-6 . Actually, 2•(H 2 O) is a manageable system to understand the impact of numerical approximations on final numerical outcomes and to rationalize the response of HB relative to coexistent polar CB. Fig. S1 includes the scheme of the geometry we use to perform numerical calculations. The structure is the configuration we obtain after fully relaxing the system with 6-31g (d,p) basis set and using an hybrid functional (HYB-xc) to introduce many-body electron-electron interactions within a reduced mean-field model [7][8][9][10] . The HYB-xc is modelled following the methodology already developed by some of us 11 . Here, HYB-xc is defined as 80% PW91 8 and 20% Fock for exchange contributions. The angle between E vector and HB axis is named Φ. The behaviour of the function is represented in the lower panels in Figure´s S2 and S3. It resembles the profile of two positive poles instead of a dipole model. This argument supports the claim in favour of a modified electrostatic picture for chemical bonds 12 . The use of the index introduced in (Supp. Mat. S1) for the description of ( ) facilitates the understanding of such electrostatic features and establishes a frame to quantify the electrical feasibility of atomic bonds. Table S1 lists the values of the indexes for all the oxygen-hydrogen bonds in the water dimer.
The values are normalized to the index obtained for the archetypal HB of the water dimer. The two-orders of magnitude that differentiate CB indexes from HB index are indicative of the higher electrical inertia of the CB when comparing with HB.

Bond Label
( )  forces and its relation with the selection of the exchange-correlation potential within ab-initio models. In regards with basis sets selection we have applied an accepted reasoning in our field, we increase the basis sets size until we fix the error below a selected limit 11 .

On the Electrical Nature of HB.
The dipole model for HB or polar covalent bonds is the most common approach for treating electrostatic interactions of these bonded systems under a classic (molecular mechanics) framework. However, the consideration of X-H bonds as dipoles is limited to effects occurring at distances longer than the bond distance, otherwise such bonds must be studied under a manybodies quantum mechanics framework.
Then, here we have introduced a novel electrical description of chemical interactions, using an improved DFT approach 11 and ad-hoc developments, Supp. Mat. 1 for further details. This approach permits to derive a mean-field structural (electron density) parameter that characterizes, similar to a classical fashion, the most native electrical features of chemical bonds. This method complements the classical (molecular mechanics) view of hydrogen bonds in terms of a dipole model.

Regarding the limitations of the dipole model for HB within the bond region.
Analysis of E-projections in the main axis of the bonds (Fig. 3 in the main text) show two subtle features to remark (i) the asymptotic behaviors in |E|, coincide with the two atom positions, i.e.  (d OH ). However, as could be seen from Fig. 3c the HB behavior seems closely related to a Heaviside Function, Fig. S5 (a). Deviations encountered in respect with this ideal behavior could be quantified by fitting multiple points (d OH , Φ (d) ) to a limiting Heaviside Function characterized giving rise to a finite t a value (see also Supp. Mat. 1).
Such behaviors numerically represented in Fig. 3 Fig. S6). Under this macroscopic view based for HB dipoles (p 1 ...p 6 in Fig. S6), the existence of a finite PZ effect is straightforward due to existence of finite resultant dipole moment. Considering crystalline symmetries, final contribution will be oriented as: p x = p 5x + p 6x (interchains HB) and p y = p 1y + p 2y (intrachain HB), see Fig. S6. Due to CB constrictions in the x-axis, we explain the atomic origin of the 2D-NC PZ response in terms of contributing p y , equivalently the dipole moments originated in intrachain HB. Figure S6. Representation of far-field HB dipole affecting PZ response in 2D-NC. p 1 , p 2 , p 5 , p 6 represents dipole moments associated to O3-H…O5, O2-H…O6, O6-H…O3' and O6'-H'…O3 respectively.

Interpolation method for the definition of piezoelectric coefficients
Results reported in the main manuscript for the PZ coefficients correspond with the slopes of the linear interpolations done through a set of data (E-field, strain). The selection of a linear interpolation method corresponds with the accepted assumption 13 that we are dealing with a linear response regime (mechanical relaxations in the nearby of the equilibrium position and E values close to electrical threshold).
Tables S2 summarizes the interpolation results for the behaviours represented in Fig. 5a. As could be distinguished, the linear coefficients (estimated PZ coefficient) of both models show the same mean value which is indicative of the robustness of the prediction and the validity of the linear regime. Although the quadratic model evidently shows a better adjustment, the linear approaches are all statistically significant (note that the models' p-values are inferior to 0.05) which also supports the linear framework.