Giant piezoelectricity driven by Thouless pump in conjugated polymers

Piezoelectricity of organic polymers has attracted increasing interest because of several advantages they exhibit over traditional inorganic ceramics. While most organic piezoelectrics rely on the presence of intrinsic local dipoles, a highly nonlocal electronic polarization can be foreseen in conjugated polymers, characterised by delocalized and highly responsive ${\pi}$-electrons. These 1D systems represent a physical realization of a Thouless pump, a mechanism of adiabatic charge transport of topological nature which results, as shown in this work, in anomalously large dynamical effective charges, inversely proportional to the band gap energy. A structural (ferroelectric) phase transition further contributes to an enhancement of the piezoelectric response reminiscent of that observed in piezoelectric perovskites close to morphotropic phase boundaries. First-principles Density Functional Theory (DFT) calculations performed in two representative conjugated polymers using hybrid functionals, show that state-of-the-art organic piezoelectric are outperformed by piezoelectric conjugated polymers, mostly thanks to strongly anomalous effective charges of carbon, larger than 5e - ordinary values being of the order of 1e - and reaching the giant value of 30e for band gaps of the order of 1 eV.


INTRODUCTION
Piezoelectricity is a well known phenomenon characterising those materials which have the property to generate a surface charge, and hence an electric tension, when subject to a stress, or conversely to deform elastically in response to an external electric field.Thanks to the possibility they offer to convert mechanical energy into electrical energy and vice-versa, piezoelectric materials are of great interest in various fields and for many applications, from macro-to microscopic electromechanical devices, to energy harvesting and much more [1,2].The most widely used piezoelectric materials are inorganic perovskites, such as lead zirconate titanate (PZT), for their high electromechanical response [3][4][5].However, inorganic materials have very low mechanical flexibility, high fabrication costs and often are toxic because of the lead they contain.These facts motivated an intense research activity aimed at developing and identifying leadfree piezoelectric ceramics [5,6].A promising alternative is to exploit piezoelectric properties of organic materials, a route that has been trodden with some success since early 1990s, mostly focusing on the wide class of organic polymers displaying high flexibility, low fabrication costs and bio-compatibility [7][8][9][10].In this context, the most studied organic piezoelectric is polyvinylidene fluoride (PVDF) [8,11], a saturated polymer derived from polyethylene and comprising molecular units with net (electrical) dipole moments, thus giving rise to ferroelectricity of conformational origin due to the rotation of chains' segments from non-polar to polar isomers.
Despite the intense research efforts and the progress made in the field of organic piezoelectrics, to date inorganic ceramics still display much better piezoelectric * francesco.mauri@uniroma1.itperformance than organic counterparts.The piezoelectric response in PZT and related inorganic materials is strongly enhanced at morphotropic phase boundaries (MPB), marking a composition-driven structural transition between two competing, nearly energetically degenerate phases with distinct symmetries [12][13][14].On general grounds, the properties enhancement close to such phase transition may be traced back to the flattening of freeenergy surfaces, easing polarization extension and/or rotation, as extensively discussed and observed mostly in perovskite oxides [15].Recently, the concept of MPB has been loosely extended to the family of P(VDF-TrFE) copolymers [16,17].Here, the introduction of different TrFE monomers in the semicrystalline PVDF structure has been proposed to lead to an enhanced conformational competition reminiscent of the structural competition realized at MPB, further suggesting an optimal chemical composition for maximizing the piezoelectric coefficient.Even though at the optimal "morphotropic" composition the piezoelectric coefficient roughly doubles the typical values of PVDF, it is still smaller by one order of magnitude compared to characteristic piezoelectric coefficients of inorganic oxide ceramics such as PZT.
Beside saturated polymers as PVDF and P(VDF-TrFE), whose ferroelectric and piezoelectric properties rely on the ordering of built-in molecular dipoles and as such require appropriate poling treatments, a natural alternative choice is represented by conjugated polymers, such as archetypical polyacetylene (PA).Characterised by a delocalised π−orbital along their backbone, these polymers are widely studied for their peculiar electronic properties [18,19].If specific symmetry-lowering effects allowing for piezoelectricity and ferroelectricity are met in a conjugated polymer, one may expect a large polar response of electronic origin due to the redistribution of the responsive π−electronic density along the polymer's backbone.Because of the delocalized nature of conjugated π−state (of Bloch-type in an infinite periodic chain), small changes of atomic positions may lead to a global shift of electrons, revealing the strong nonlocal and ultimately topological character of electronic polarization in these quasi-1D systems [20].Indeed, conjugated polymers have been theoretically put forward as a potential new class of "electronic ferroelectrics" [21].
In order to explore the potential piezoelectricity of conjugated polymers, we adopt the simple model originally introduced by Rice and Mele to study soliton excitations of linearly conjugated diatomic polymers [22], later become a prototypical model for 1D ferroelectricity [20,23] and Thouless adiabatic pumping [24,25].Building on these well-established results, we study how the interplay between these mechanisms affect the system response to a mechanical strain.We find that polar responses, such as Born effective-charge tensors -also known as atomic polar tensors -and piezoelectric coefficients, can be strongly enhanced close to the dimerization phase transition characterized by a bond-length alternation (BLA) between neighboring atoms along the chain.This transition can be controlled by changes in the chemical composition of the polymers or by tuning the electron-phonon (e-ph) interaction, that is expected to be sensitive to the dielectric environment of individual chains.The strength of the piezoelectric effect is found to be determined partly by the MPB-like enhancement -relying on the second-order nature of the transition -, but mostly by a giant effective charge, ultimately due to the topological nature of Thouless' adiabatic charge transport in this class of materials.We verified the model-based predictions and theoretical picture by performing Density Functional Theory (DFT) calculations using a range-separated hybrid functional on two prototypical conjugated polymers derived from polyacetylene, confirming strong enhancement of effective charges and consequently of piezoelectric response, that is found to be larger than the one observed in the celebrated organic piezoelectric PVDF for a wide range of parameters.

A. A model for piezoelectric conjugated polymers
To study the effects of strain on conjugated polymers, we generalise the Rice-Mele model [22] of a 1D linear chain with a unit cell of length a 0 containing two atoms.The only possible strains for a 1D infinite periodic chain are contraction/dilatation of the unit cell, defined by a(ϵ) = a 0 (1 + ϵ), where the adimensional parameter ϵ quantifies the strain along the direction of the chain.In a nearest neighbours tight-binding approximation, the modulation of the hopping energy between atoms at site i and i + 1 can be modelled at linear order in atoms' displacement and strain as t i,i+1 = t(ϵ) + (−1) i+1 δt(ϵ), with Equation ( 1) accounts for the effect of strain on the hopping energy t 0 between equidistant atoms, while Equation (2) describes the variation with respect to t(ϵ) caused by atoms' displacement; at the lowest order, we can assume the same adimensional e-ph coupling parameter β > 0 in describing both effects.The term u(ϵ) is an adimensional fractional coordinate measuring the relative displacement of the two atoms in the unit cell, as defined in the Methods.We indicate with the term ∆ ≥ 0 the on-site energy difference between neighbours.If ∆ = 0, the atoms in the unit cell are equivalent and we recover the well-known SSH model [26], used to describe polyacetylene.However, in order to have a non-trivial polar response, it is necessary to break atoms' equivalence (∆ ̸ = 0) as, e.g., in substituted polyacetylenes (SPA) [27], a class of conjugated polymers formed by inequivalent monomers which can be obtained substituting (one of) the atoms in the C 2 H 2 unit of PA with some element(s) or compound(s).A representation of the model is in the insets of Figure 1a) and 1b).For simplicity, we consider only longitudinal displacements, parallel to the linearchain direction.As we are interested in the equilibrium structures at T = 0 K, we define the optimal displacement u as the one which minimises the total energy, given a set of material-dependent parameters and a strain.For more information on the model and the derivation of the above quantities, see the Methods and the Supplementary Information.
In the absence of strain (ϵ = 0), the properties of the model are well known.In the SSH model (∆ = 0) a finite u ̸ = 0 is allowed by an infinitely small e-ph interaction, conjuring with a Fermi-surface nesting and Peierls electronic instability to produce a dimerized phase with bond-length alternation and the opening of a gap in the energy spectrum E gap = 4|δt| = 4βt 0 u.Breaking the equivalence between atoms (∆ ̸ = 0) also opens a gap which is expected to counteract the Peierls instability and the related bond dimerization.The lowersymmetry structure can be stabilised only by a finite and strong enough e-ph interaction, resulting in a gap Indeed, as shown in Figure 1a), increasing ∆ at fixed e-ph coupling induces a second order structural phase transition at a critical value ∆ c with order parameter u between a distorted (u ̸ = 0) and a higher-symmetric undimerized (u = 0) phase, both displaying an insulating character with a finite band-gap (see Supplementary Information).Such structural transition from the dimerized to undimerized phase can also be induced by decreasing the e-ph coupling at fixed non-zero onsite energy difference, as discussed in the next section and displayed in Figure 1b).In a 2D parametric (∆, δt)-space, where δt depends -through optimal u -on both ∆ and β, the origin of the axes corresponds to a metallic system while the rest of the space to an insulator.Since the milestone paper by Vanderbilt and King-Smith [24], it is known that if this system undergoes an adiabatic evolution along a loop enclosing the origin of this space, a quantised charge |e| is pumped out.This is a prototypical example of the adiabatic charge transport known as Thouless' pump [25,28] realised e.g. in ultracold fermions [29].The crucial ingredient is the presence of the metallic point in the domain enclosed by the loop and, in this sense, it is an example of topological phenomenon (see Supplementary Information).Previous studies on low-dimensional systems showed how topology may reflect onto polar responses, resulting, e.g., in electronic polarization in quasi-1D systems [20] or in large effective charges and piezoelectric coefficients directly related to valley Chern numbers in 2D systems [30,31].In the next two Sections, we use the generalized Rice-Mele model to show how the Thouless-pump topological mechanism conjures with the bond dimerization to give raise to an enhanced piezoelectric effect, and we discuss how the latter can be tuned by acting on the two independent parameters ∆ and β.

B. Morphotropic-like enhancement of the piezoelectric response
In order to have a non-trivial piezoelectric response, the chain must not have points of inversion symmetry, a requirement that is met when the equivalence between atoms is broken in a distorted chain, i.e. both ∆ ̸ = 0 and u ̸ = 0.In this case, the chain becomes also ferroelectric with a net dipole moment per unit cell P [23].The electromechanical response is quantified by the piezoelectric coefficient c piezo , defined as the variation of P due to an applied homogeneous strain ϵ, namely where the derivative of P is decomposed in two contributions.The first one is the so-called clamped ions term which is obtained keeping fixed the relative position of the ions in the unit cell, i.e., for fixed internal fractional coordinate u 0 = u(0).It can be shown (see Methods) that |c c.i. piezo | ≤ |e|β/2π.The second term of Equation (3) takes into account the effect of strain on the internal coordinate u(ϵ) and defines the internal relaxation contribution where we defined the effective charge Z * , namely a measure of how rigidly the electronic charge distribution follows the displacement of the nuclei, as The inclusion of strain in the Rice-Mele model affects explicitly the critical value ∆ c (ϵ) of the phase transition.The behaviour of u(ϵ) as ∆ approaches ∆ c (ϵ) is shown in Figure 1a), following the expected behaviour u(ϵ) ∝ |∆ − ∆ c (ϵ)| 1/2 of the order parameter of second-order phase transitions (see Supplementary Information).It thus follows that: Equation (7) implies that the internal relaxation term diverges as we approach the critical point ∆ c (ϵ) from the distorted phase, in analogy with the MPB mechanism at play in some ferroelectric oxides.Indeed, as the ∆ parameter of the Rice-Mele model accounts for the composition of the system, it allows to continuously tune a morphotropic-like phase transition from the distorted phase (lower symmetry, ferroelectric) to the undistorted one (higher symmetry, paraelectric).On the other hand, at a fixed ∆ ̸ = 0 suppressing the Peierls electronic instability, the second-order phase transition can be driven by the e-ph coupling, as shown in Figure 1b), allowing one to define a critical value β c separating the undimerized and dimerized phases (see Supplementary Information).
Clearly, the order parameter as a function of the variable driving the transition would still follow the expected be- i.e., a diverging internal strain when approaching the critical point from the dimerized phase.

C. Topological contribution to the enhancement
In principle, the diverging behaviour of the internal strain, Equation (7) or (8), guarantees the existence of piezoelectric polymers with arbitrarily high response when close to a morphotropic-like phase boundary, irrespective of the prefactor, namely the effective charge Z * .However, this specific enhancement is a consequence of the second order transition.Numerical evidence in linear acetylenic carbon chains [32] show that quantum anharmonic effects (QAE) may change the order of the structural phase transition, therefore damping the diverging behaviour of c i.r.piezo .A robust enhancement of the piezoelectric coefficient against QAE would depend, therefore, on the strength of the polar response embodied by Z * .Measuring how the electronic charge distribution follows the displacement of the nuclei, effective charge's behaviour in this system is strictly related to the topological charge transport of the Thouless pump [28].A simple geometric argument shows that indeed the effective charge Z * diverges as 1/E gap in the undimerized phase (see Methods and Supplementary Information).As first  ).Near the critical point a huge polar response is also present, quantified by the effective charge Z * ∝ β∆/E 2 gap (see Methods).Approaching the critical point from the undimerized phase, Z * is inversely proportional to the gap Egap(u = 0) ≡ ∆ (c) and linear in β (d).Since the gap is constant as a function of ∆ and linear in β when approaching the critical point from the dimerized phase, Z * displays a linear behaviour in ∆ (c) and it is inversely proportional to β (d).(e) and (f ) The piezoelectric coefficient diverges when approaching the critical point from the dimerized phase.As shown in the inset, the major contribution is due to the internal relaxation term of Equation ( 5).The topological nature of the enhancement guarantees its stability.
noted in Ref. [20], such a remarkable behaviour is a consequence of the topology of the domain of the dipole moment P .As E gap goes to 0, we get closer to the metallic point, where P is not well defined, and even an infinitesimal atomic displacement causes a huge redistribution of the charge density.We contrast this result with the predicted behaviour of polar responses in 2D gapped graphene, where both piezoelectric coefficient and effective charges were found to be independent on the bandgap amplitude [31].Indeed, electron-strain/lattice couplings in 2D hexagonal crystals can be described as gauge fields [30,31] whose effect is to shift the Dirac cone of an amount proportional to the coupling constants, causing the latter to be the only relevant quantities determining the strength of polar responses.In the 1D chain, instead, the e-ph interaction contributes, through dimerization, to the gap opening, thus directly affecting the Thoulesspump topological enhancement of effective charge.We further remark that the absence of a structural transition in gapped graphene causes the piezoelectric response to be mostly due to the clamped-ion contribution, the internal-relaxation one contributing roughly 25% to the total response [31].
The evolution of Z * as a function of ∆/∆ c and β/β c is shown in Figure 1c) and 1d).Even though the system always displays a finite gap preventing the metallic divergence of the effective charge, Z * reaches the giant value of ∼ 30 |e| at the critical points.Such anomalously large effective charges can not be ascribed to the mixed covalent-ionic character of the system only: it can be shown (see Supplementary Information for the simple case of a heteropolar biatomic molecule) that tuning the bond character can lead to a finite enhancement of the effective charges, which are typically only few times the value of the nominal charge [33].Unlike the MPB-related enhancement of the internal strain, also shown in Figures 1c) and 1d), the topological behaviour is expected to be much more stable with respect to QAE, guaranteeing the enhancement of the electromechanical response.The total piezoelectric coefficient, comprising both the clamped ion and internal relaxation contributions, is shown in Fig- ures 1e) and f) as a function of parameters ∆/∆ c and β/β c .Insets highlight how the piezoelectric coefficient is mostly contributed by the internal-relaxation contribution, that is strongly enhanced by the combined effect of diverging internal strain ∂ ū/∂ϵ and anomalously large effective charges.We remark the importance of both mechanisms, since anomalous effective charges alone in general do not guarantee piezoelectric effects if inversion symmetry is kept, as in centrosymmetric CaTiO 3 and SrTiO 3 [33], or if the internal strains are small, as in 2D hexagonal systems and gapped graphene [31].

D. Numerical calculations
We performed ab initio calculations in the framework of Density Functional Theory (DFT) to validate our model predictions, choosing two conjugate polymers representative of the broad class of SPA.One is monofluorinated polyacetylene (MFPA) made by the repetition of the unit CH-CF, obtained by substituting one hydrogen atom of the C 2 H 2 unit of PA with fluorine.The other is polymethineimine (PMI), obtained substituting a CH pair with a nitrogen atom to obtain the unit N-CH.For simplicity, we considered the all-trans structures shown in Figure 2a) and 2b), whose fundamental physical properties are captured by the Rice-Mele model.Even though controlling the fraction of substituted atoms may in principle induce a morphotropic-like transition, this approach poses many challenges both from the computational and experimental side: to our specific purposes, it wouldn't allow to study the phase transition and the associated predicted enhancement of piezoelectric effect by varying with continuity an external parameter (as ∆ in the model).As discussed in sections B and C, the internalstrain and the effective-charge enhancements may be also induced by tuning the parameter β.To achieve this computational task, we take advantage of the effect of screened Coulomb vertex corrections to the dressing of the e-ph coupling [34], leading to an enhancement of e-ph itself especially strong in low-dimensional materials and when phonons at zone boundary are involved [35][36][37][38].Such screened-Coulomb-mediated e-ph enhancement can be captured by hybrid functionals incorporating a fraction of the exact exchange [36,39].Its inclusion has been proven essential for describing the bond-length alternation of trans-polyacetylene [40,41] and related 1D polymers [32,42], whose BLA is typically underestimated by standard local-density or generalized-gradient approximations, indirectly pointing to an enhancement of the eph coupling due to electron-electron interaction.Rangeseparated hybrid (RSH) functionals represent an ideal choice for our purpose, as they have been designed to better account for the screened Coulomb vertex corrections.The latter can be effectively tuned by acting on the longrange (LR) mixing parameter c LR that accounts for the fraction of LR exact exchange in RSH functionals, thus providing a computational knob to continuously vary β.We further remark that the strength of e-ph enhancement due to screened Coulomb effects can be ideally controlled by modifying the screening itself, as proposed for doped graphene [43].Since the optimal mixing parameter c LR is inversely proportional to the scalar dielectric constant of the environment in order to enforce the correct asymptotic potential [44][45][46][47][48], we speculate that controlling the dielectric environment may represent a viable strategy, alternative and complementary to controlling the fraction of substituted atoms, for tuning and optimising the piezoelectric response of conjugated polymers.
Motivated by these reasons, we performed structural optimization of both MFPA and PMI for different values of the LR mixing parameter c LR .For consistency with the model, we considered the coordinates of the C and N atoms along the principal axis of the chain, which we take as the x-axis, to compute the internal coordinate.More details on the effect of c LR on polymers' structures are provided in the Supplementary Information.The evolution of u 0 displayed in Figure 2c) and 2d), clearly hints to the presence of a second order phase transition triggered by c LR for both polymers: the dimerized phase is suppressed by lowering the fraction of mixing, the order parameter showing the expected behaviour as it approaches the second-order phase-transition critical point, in excellent qualitative agreement with model results shown in Figure 1b), d) and f) and confirming a posteriori the direct proportionality between c LR and β.The second-order character of the phase transition is further confirmed by the softening of the corresponding optical phonon found in the higher-symmetry phase when increasing c LR (Figure S6), signalling the onset of a dynamical instability of the undimerized structure.
In Figure 2e) and 2f) the behaviour of Z * and of ∂u/∂ϵ calculated from first principles is shown.For MFPA we took Z * = Z * C F ,xx where C F is the carbon atom bound to the fluorine, while for PMI Z * = Z * C,xx .The full tensors of the effective charges of all the atoms are reported in the Supplementary Information.We highlight the qualitative agreement with the prediction of the model, in particular the huge enhancement of the effective charges around (for MFPA the C bound to the F).In agreement with the model, on the one hand we observe a further hint of the morphotropiclike nature of the transition while on the other hand the huge values of Z * stands out, in particular in the region near the critical points.In (g) and (f ) is shown how the piezoelectric coefficient is greatly enhanced when reaching the critical points from the less symmetric phase.The comparison between c c.i. piezo and c i.r.piezo , in the insets, highlights the internal-relaxation origin of the enhancement.Furthermore, the comparison with the results obtained putting the values of (e) and (f ) in Equation (5) shows that the model very well describes the nature of the enhancement.
the critical point c LR , reaching the strongly anomalous values of ∼ 30 |e| and ∼ 15 |e| in correspondence of c LR for MFPA and PMI, respectively.The covalent character of bonds along the chain prevents a precise definition of the nominal reference value for C, that can be however assumed to be of the order of 1|e|, as the nominal ionic charges for H and F are respectively +1|e| and −1|e|.Effective charges of carbon in both considered chains are strongly anomalous for all considered long-range mix-ing parameters, displaying values between 5|e| and 30|e| even for band gaps exceeding 6 eV.These anomalous values exceed even those reported in oxide ferroelectrics, where effective charges are typically two or three times larger than nominal reference values [33].Figure 2g) and 2h) display the behaviour of the piezoelectric coefficients computed ab initio taking into account also the effects of transverse displacements.In the insets, the different contributions c i.r.piezo and c c.i. piezo are compared, highlighting the internal-relaxation origin of the enhancement.Using the the ab initio values of Z * and ∂u/∂ϵ of Figure 2e) and 2f), we compare the values of c i.r.piezo computed without approximations with those obtained using Equation (5) of the model.The agreement between the two approaches is both qualitatively and quantitatively excellent, notwithstanding the simplifying description provided by the Rice-Mele model, that neglects structural details specific of the two considered polymers as well as transverse displacements.We highlight that despite the behaviour ∂u/∂ϵ ∝ |c LR − c LR | −1/2 , the main contribution to the piezoelectric coefficient is given by the effective charges.The large values attained in a finite range around the second-order critical point and their ultimately topological origin suggest that the piezoelectric effect is robust against quantum and anharmonic effects that may change the order of the phase transition, as predicted in carbyne [32].

DISCUSSION
Finally, we compare the results for the piezoelectric coefficients of MFPA and PMI with those of the best and most widely used piezoelectric PVDF polymer, made by the repetition of the unit CF 2 -CH 2 .The Rice-Mele model fails to capture its main properties, this polymer being not conjugated.Its piezoelectricity indeed derives from the presence of a net dipole moment transverse to the chain, whereas the electromechanical response predicted in conjugated polymers is longitudinal to the chain and ultimately due to the topological-morphotropic enhancement.To have a consistent comparison with available experimental data for PVDF, we computed ab initio the converse piezoelectric coefficient d piezo , which measures the response with respect to an external stress, rather than a strain (see Methods).The results for the converse piezoelectric coefficients of PVDF are compared with those of MFPA and PMI in Figure 3 and are consistent with the values computed in Ref. [49].Even though the calculated d piezo is smaller than reported experimental values, a direct comparison to experiments is hardly drawn because, e.g., of the polymorphic character or low crystallinity of experimental samples, as noticed also in Ref. [49].We remark that piezoelectric response in PVDF is found to be independent on the fraction of exact exchange, confirming the utterly different nature of the electromechanical response in such non-conjugated polymer.We finally mention that mildly anomalous effective charges have been also reported for PVDF [50], the carbon effective charge however not exceeding 1.5|e|, consistently with our results provided in Supplementary Information.On the other hand, both MFPA and PMI display a rather large range of values that are larger than calculated d piezo of PVDF, with up to a six-fold enhancement for PMI close to the dimerization point.The robustness of the enhancement mechanism is confirmed also by the comparison with the piezoelectric coefficients calculated in packed MFPA chains, as discussed in Supplementary Information.Even though the selected prototypical SPAs may not be the most efficient ones for practical realization and engineering of their functional properties, the comparison with first-principles estimate of one of the best available piezoelectric polymer alongside the general validity of the proposed model and consequent robustness of its electromechanical response put forward the broad class of π-conjugated polymers as a promising field for organic piezoelectrics with enhanced functionalities.Additionally, its inverse proportionality to the band gap provides a possible material-design principle for driving the quest of organic polymers with enhanced piezoelectric response.

Details on the model
The Hamiltonian of the chain is the sum of two terms H tot = H L + H e .The lattice contribution H L accounts for the displacement of the atoms with respect to their position in a uniformly spaced chain: where m i and δr i are mass and displacement of atom i, p i its momentum, and K is a spring constant term.In a nearest neighbours tight-binding approximation, the electronic contribution H e reads where the factor 2 accounts for the spin degeneracy, c † i /c i are creation/annihilation operators for electrons, t i,i+1 is a hopping energy while ∆ > 0 accounts for the onsite energy difference between neighbours.The adimensional order parameter u(ϵ) measures the relative displacement between neighbours and is defined as The total energy of the system E tot (u) as a function of the order parameter u has two contributions, namely The first term E L (u) accounts for lattice distortion and is quadratic in u, whereas the electronic contribution E e (u) is linear in u.For a given set of material-dependent parameter t 0 , K, ∆, β and the strain ϵ, we find the optimal displacement u -defined as the one which minimises the total energy -as shown in detail in the Supplementary Information.
We computed the dipole moment per unit cell P using the Berry phase approach [51,52].In the limit E gap ≪ t(ϵ) it holds (see Supplementary Information) with tan θ(ϵ, u(ϵ)) = ∆/2δt(ϵ, u(ϵ)), where we emphasize that the dependence on the strain enters both directly and through u(ϵ).Using Equation ( 2) and ( 13) we have where we notice that |c c.i. piezo | ≤ |e|β/2π, the maximum achievable value being directly proportional to the e-ph coupling constant β.
Plugging Equation (2) in Equation ( 13), we obtain an explicit expression for the effective charge in both the dimerized and undimerized phases: where sin θ = 2∆/E gap .In the undimerized phase u = 0 implies sin θ = 1, hence Z * ∼ β/E gap ∼ β/∆.The diverging behaviour can be further understood using a simple geometric argument that highlights its topological origin.From Equation ( 13), polarization is proportional to an angle θ spanning the parametric (∆, δt)-space along a circumference with radius E gap .As δt ∝ u, it follows that in the undimerized phase E gap dθ ∝ du, hence Z * = ∂P/∂u ∝ ∂θ/∂u ∝ 1/E gap (more details in the Supplementary Information).
The parameters used to produce the results in Figure 1 were obtained fitting the model with the PBE0 [53] total energy profile of carbyne as a function of the relative displacement of the two carbon atoms of its unit cell.In particular, with the value a 0 = 2.534 Å obtained through a cell-relaxation procedure, the fit yields t 0 = 2.239 eV, β = 3.906 and K = 127.9eV/ Å2 , whereas for carbyne it holds ∆ = 0.The behaviour with respect to β/β c was obtained fixing a finite ∆ = 1.05∆ c .We highlight that each carbon atom of carbyne contributes with two electrons to the π-orbital, so it is necessary to put an additional factor of 2 in front of Equation (10) to account for this degeneracy.

Computational details
All DFT calculations were performed using CRYSTAL code [54,55], which employs a basis of local Gaussiantype functions.This approach allows for the simulation of truly isolated systems, as the 1D polymers addressed in our work, using the hybrid functionals, proven to be essential for accurately reproducing the physics of 1D chains, [32,41].The Gaussian-type basis set significantly reduces the computational cost of evaluating real-space integrals and, hence, it allows to drastically reduce the computational cost when compared, e.g., with state-of-the-art plane-waves based DFT codes.We used a triple-ζ-polarised Gaussian-type basis [56] with real space integration tolerances of 10-10-10-15-30 and an energy tolerance of 10 −10 Ha for the total energy convergence.We customised a range-separated LC-ωPBE hybrid exchange-correlation functional [57] varying the value of the long-range (LR) mixing parameter c LR which enters in the definition of the LR part of the functional, namely When c LR = 0 the PBE functional is recovered while if c LR = 1 we have pure Hartree-Fock (HF) exchange.The long-range terms in round brackets depend on the rangeseparation parameter ω that enters in the decomposition of the Coulomb operator 1/r as where erf(•) is the error function; the first and second term in the right-hand side of Eq. 17 account for the short-and long-range part of the Coulomb operator, respectively.All the presented values were obtained with ω = 0.4 a −1 0 .For each c LR , a geometric optimisation was performed and all quantities were computed on the equilibrium configurations (see Supplementary Information).The derivatives ∂u/∂ϵ of the order parameter with respect to the strain were computed with finite differences, performing a fixed-cell optimisation for each strained configuration with cell length a(ϵ) = a 0 (1 ± ϵ) and ϵ = 0.01, while the effective charges were computed as finite differences of polarization.Values of polarization along the chain (parallel to the x-axis) were computed using the Berry phase approach, whereas components transverse to the isolated chain were computed in real space.[58,59].The piezoelectric coefficients c piezo and d piezo of Figure 2g), 2h) and 3, as well as the values of c c.i. piezo and of c i.r.piezo , were computed using the Berryphase approach [60] as implemented in the code [61,62], which accounts also for transverse displacements.The converse piezoelectric coefficient d piezo is linearly related to c piezo through the elastic constants tensor C, namely c piezo = d piezo C. As far as 1D systems are concerned, only a single scalar elastic constant is required, and it can be evaluated as the second derivative of the energy with respect to the strain, i.e.C = ∂ 2 E/∂ϵ 2 .orbital of atom α = A, B is located at r α + R, R = na (n ∈ Z) being the position of the atom's cell along the chain.Without loss of generality, we take r A = −a/4 + δr A and r B = +a/4 + δr B , δr α being the displacement of atom α.For simplicity, we consider only longitudinal displacements, parallel to the linear-chain direction.The basis set of the orbitals {|α, R⟩} is orthonormal and it holds: We define the onsite energy terms of the electronic Hamiltonian H e as We take into account also the energetic contribution due to the overlap between an atom's orbital and the orbitals of its left and right nearest neighbours.In general, the overlap energy term between two electronic orbitals is a function of the distance r between the atoms and is referred to as the hopping energy −t(r), t > 0. For convenience, we indicate with r 1 the distance between atoms in the same cell and with r 2 the distance between neighbouring atoms in adjacent cells, namely: and it holds r 1 +r 2 = a.We can now define the hopping energy between orbitals of atoms in the same cell −t 1 ≡ −t(r 1 ) and the hopping energy between atoms in adjacent cells −t 2 ≡ −t(r 2 ).With the same notation of Equation (S3) we write: As we are interested in the effects of strain ϵ and of atoms' displacement, we define, using Equation (S4) and (S5), an adimensional fractional coordinate u(ϵ): This term quantifies deviations from the equally spaced chain, allowing us to express bond lengths with the following compact expression: With the above definitions, at linear order in atoms' displacement we have which allow us to define the two terms The former term t(ϵ) quantifies the effect of strain on the hopping energy between equidistant atoms while the latter term δt(ϵ) describes the variation with respect to t(ϵ) caused by atoms' relative displacement.In absence of strain (ϵ = 0) we recover the quantities defined in the Rice-Mele model, whereas at linear order in ϵ it holds: where t 0 ≡ t(a 0 /2) and we defined the adimensional parameter β > 0 as This term quantifies the variation of the hopping energy due to a variation of the distance between the atoms and as such it acts as an electron-phonon coupling term.Analogously, we obtain where we define another adimensional e-ph parameter β ′ > 0 as Even though the two e-ph parameters β ′ , β can differ at finite values of the strain, at the lowest order one can safely assume that they coincide and consider β ′ = β.A schematic representation of the model is shown as an inset in Figure S2.

STRUCTURAL PHASE TRANSITION
As we are interested in the structural properties at T = 0 K, we study the total energy per unit cell , where E L (u) and E e (u) are the lattice and electronic contribution, respectively.In particular, we aim to characterise the behaviour of the optimal displacement u, defined as the one which minimises E tot (u) given a set of material-dependent parameters t 0 , K, ∆, β and a strain ϵ.Lattice dynamics being neglected, we write the lattice contribution, which accounts for the displacement of the atoms with respect to their position in a uniformly spaced chain: where we used Equation (S7) and K is an elastic constant term.
To compute the electronic energy per unit cell E e (u), we imagine the linear chain as made of N copies of the unit cell and we adopt periodic boundary conditions.This allows us to define a basis {|α, k⟩} in the reciprocal k-space: where k is defined over the first Brillouin zone, namely and it holds From Equation (S3), (S6) and (S18) we obtain the matrix elements of the electronic Hamiltonian in the reciprocal space: The above results allow to write in the reciprocal space basis a 2 × 2 electronic Hamiltonian matrix H e,k for each k-point: Diagonalising the matrix of Equation (S25) we find the two eigenvalues with the respective eigenvectors, namely the Bloch wave-functions and using Equation (S10), (S11) and S23 we have that For each k-point, Equation (S26) allows to distinguish between a lower and a higher energy level.As shown in Figure S1, the chain present two energy bands.To describe the π−orbitals of conjugated orbitals, we assume that only the lower energy band is filled with electrons and will hereafter refer to it as the occupied band, in contrast with the unoccupied higher energy band.From Equation (S26) we obtain the value of the energy gap E gap between the occupied and unoccupied band: We consider the contribution of all the occupied states to the electronic energy per unit cell and in particular it holds where the factor of 2 accounts for the spin degeneracy.In the limit of an infinite linear chain, namely for N → +∞, the eigenvalues and the eigenvectors' coefficients become continuous functions of k and the sum becomes an integral over the first Brillouin zone: Finally, using Equation (S15), (S17) and (S31), we write the total energy per unit cell as a function of the fractional coordinate u: We can now study the structural properties of the system, encompassed in the optimal displacement u, defined as the one which minimises E tot (u) = E L (u) + E e (u) given a set of material-dependent parameters t 0 , K, ∆, β and a strain

ELECTRONIC POLARIZATION AND THOULESS PUMP
In this section, we derive the expression of the dipole moment per unit cell P as obtained within the modern theory of polarization [51].In this framework, the wave function of the occupied states is required to be continuous at the edges of the Brillouin zone.To satisfy the requirement, we multiply the wave function for the occupied states of Equation (S27) by a phase factor and define with r A = −a/4 + δr A .The dipole moment per unit cell P is defined in term of the Berry phase [63] φ as where e is the electron's charge, the factor 2 at the numerator accounts for the spin degeneracy and the Berry phase φ is defined as We indicated with | ψ occ k ⟩ the periodic part of the wave function |ψ occ k ⟩ and from Equation (S18), (S27) and (S40) it holds where As the scalar product in Equation (S42) is taken over a single unit cell, without loss of generality we consider only the contribution for R = 0 in Equation (S43) and obtain From the gap Equation (S29), we notice that the insulating/metallic character of the system can be visualised in a 2D parametric (∆, δt)−space, where the origin of the axes correspond to a metallic system with E gap = 0, while every other point correspond to an insulating system with E gap ̸ = 0.In this space, we define a parameter θ that allows to identify each point with the polar coordinates (E gap , θ), with the change of coordinates defined by 2∆ = E gap sin θ (S47) 4δt = E gap cos θ. (S48) Applying this change of coordinates in Equation (S46), in the limit E gap ≪ t 0 it can be demonstrated that it holds This property can be appreciated in Figure S3, where we compare, for different values of E gap /t 0 , the behaviour with respect to θ of the dipole moment P , obtained computing Equation (S41) numerically without further approximations.As θ varies, so do the terms ∆ and δt.In particular, a full rotation of 2π implies that the system has returned in its initial state, while the dipole moment P has acquired a quantum of −2|e|.This properties is more general: if the system undergoes an adiabatic evolution along any loop enclosing the origin of the 2D space, a quantised charge is always pumped out.This phenomenon is known as adiabatic charge transport or Thouless' pump [28] and the key ingredient is the presence of the metallic point in the domain enclosed by the loop: in this sense, it is an example of topological phenomenon.This peculiar topology-related property has consequences also in the behaviour of the effective charges of the system.Defined as the variation of polarization due to an atomic displacement, the effective charges are a measure of how rigidly the electronic charge distribution follows the displacement of the nuclei.In the model it holds With the definitions given above one can compute analytically the expression for Z * , however it is possible to appreciate the main feature of its behaviour thanks to a geometric argument in the (∆, δt)-space.The effect of a small displacement on a system with a finite ∆ ̸ = 0 is to span an angle dθ in the space, starting from θ = π/2.Considering that δt ∝ u and that, with a suitable choice of the axis, the points on a circumference in this space correspond to systems with the same E gap , we have that dθ ∝ u/E gap .From Equation (S49), it holds dP ∝ dθ, hence we obtain implying that as we get closer to origin of the axis (E gap goes to 0), even an infinitesimal atomic displacement causes a huge redistribution of the charge density.This geometric argument can be appreciated in Figure S3.

DYNAMICAL CHARGES OF A HETEROPOLAR DIATOMIC MOLECULE
In this Section we compute the dynamical charge of a heteropolar diatomic molecule along the lines discussed in the appendix of Ref. [33], in order to contrast our predicted topological enhancement of effective charge and the contribution of the mixed ionic-covalent character to "anomalous" Born effective charges.Let's consider a diatomic molecule with two monovalent atoms A and B positioned along the x-axis at a distance d AB = R B − R A > 0. In a LCAO tight-binding approach, the electronic Hamiltonian H e reads where t(d AB ) is the hopping energy between the atoms, that will in general depend on some power of the inverse distance.For the sake of clarity and without loss of generality, we assume t(d AB ) ∝ 1/d 2 AB .We can find the occupied electronic orbital |ψ occ e ⟩ diagonalizing H e :  As expected, in the small-gap limit (Egap ≲ t0) we observe the linear behaviour P ∝ −θ predicted in Equation (S49).On the panel, the evolution of P is displayed in the 2D parametric (∆, δt)−space, where each point correspond to a physical realization of the model.An adiabatic evolution along a loop containing the origin of the axes -which correspond to a metallic system -pumps out a quantised charge of −2|e|.In the 2D space on the central panel it is possible to appreciate the geometric argument that yields Equation (S51).Being the polarization proportional to the angle θ spanning the parametric (∆, δt)-space along a circumference with radius Egap, as δt ∝ u, it follows that in the undimerized phase Egapdθ ∝ du, hence Z * = ∂P/∂u ∝ ∂θ/∂u ∝ 1/Egap.Finally, it is interesting to visualize the phase transitions of u0(∆) and u0(β), displayed separately in the two right panels, also as paths in the 2D parametric space.Being ∆ and β the independent parameters and δt ∝ βu, one finds that the transition in ∆ from the dimerized (u0 ̸ = 0, hence δt ̸ = 0) to the undimerized (u0 = 0, hence δt = 0) phase, happens on a circumference with fixed Egap, whereas the evolution with respect to β in the dimerized phase corresponds to a vertical line, the undimerized phase corresponding to a single point at fixed ∆ ̸ = 0 and δt = 0. Arrows on the two paths in the 2D space correspond, respectively, to the arrows in the right panels whereas all the dots of the lower right panel correspond to a single point in the 2D space in (∆ = 1.05∆c, δt = 0).
where Z v A = Z v B = +1|e| are the valence charges of atom A and B, respectively, and r is the position operator, whose matrix elements are well defined in an isolated molecule.The dynamical charge Z * α of atom α = A, B is defined as the derivative of D(d AB ) with respect to the atomic displacement R α , and using all the above definitions we obtain (in units of |e|) where i A = 1 and i B = 2.The above equations tell us that acting on the ionic-covalent character of the bond, accounted for by the term X = ∆/t, one can tune the values of the dynamical charges.Albeit the expression of Eq.S56 apparently reminds the dependence on the gap of the effective charge in the undimerized 1D chain provided in Eq.S51, as it inversely depends on the difference between molecular energy levels, from Eq. S57 it is clear that the effective charge is always limited and shows no diverging behaviour.The dynamical charges of atoms A and B, which obeys the charge-neutrality sum rule Z * A + Z * B = 0, are displayed in Figure S4 as a function of the ioniccovalent character X = ∆/t.The maximum enhancement is indeed found for finite values of X, i.e., arising from the mixed ionic-covalent character of the bond, reaching however a finite value that is roughly ∼ 1.5 the nominal value, consistently with the enhancement reported in [33].We contrast this result with the enhancement of up to 30 times the nominal charge that we find in the 1D chain and highlight the differences with the analytical behaviour found in the model, namely with Equation (S51).Born effective charges of PVDF TABLE S3.Components of the effective charge tensors of the atoms of the unit cell of PVDF.For each component, values reported are obtained as average of those computed with different values of cLR, the highest standard deviation being of the order of 10 −2 .The very weak dependence of the effective charges of PVDF on cLR is a direct consequence of the fact that the structure of PVDF does not depend on cLR, as observed in Figure S8.Our results for the Born effective charges of PVDF are overall consistent with previously reported values calculated within the generalized-gradient approximation [50].The sum rule x y z x -0.603 0.000 0.000 y 0.000 0.012 0.000 z -0.001 0.001 -0.092 x y z x 1.333 0.000 0.000 y 0.000 1.165 0.000 z 0.001 0.000 0.955 x y z x 0.069 0.000 0.000 y 0.000 0.002 -0.048 z 0.000 -0.060 0.034 x y z x -0.434 0.000 0.000 y -0.002 -0.590 -0.292 z 0.000 -0.248 -0.467 x y z x 0.069 0.000 0.000 y 0.000 0.003 0.048 z 0.000 0.060 0.035 S9.On the left panel, piezoelectric coefficients of MFPA in its isolated polymeric chain configuration are compared to those obtained for packed chains.The 3D structure adopted, displayed on the right panel, is similar to the one adopted for PVDF in Ref. [49].The agreement between the results highlight the robustness of the enhancement mechanism.

FIG. 1 .
FIG. 1.(a) and (b)Structural phase transition with order parameter u as a function of the onsite energy difference ∆ and of the e-ph coupling β.When ∆ < ∆c (β > βc) the total energy Etot(u) has a double-well profile with two minima at |u| ̸ = 0, resulting in a distorted chain with bond-length alternation.When ∆ > ∆c (β < βc) the minimum of Etot(u) is at u = 0 and the atoms become equidistant.The behaviour of the order parameter is typical of second order phase transitions, being u ∝ |∆ − ∆c| 1/2 (u ∝ |β − βc| 1/2 ), and the system remains insulating in both phases with gap Egap = 2 ∆ 2 + 4(βt0u) 2 .(c) and (d) Internal strain and effective charge across the phase transition.The inclusion of strain ϵ in the model affects the critical value ∆c(ϵ) (βc(ϵ)) and, as expected in a second order transition, the internal strain displays a diverging behaviour close to the critical point, ∂u/∂ϵ ∝ |∆ − ∆c| −1/2 (∂u/∂ϵ ∝ |β − βc| −1/2 ).Near the critical point a huge polar response is also present, quantified by the effective charge Z * ∝ β∆/E 2 gap (see Methods).Approaching the critical point from the undimerized phase, Z * is inversely proportional to the gap Egap(u = 0) ≡ ∆ (c) and linear in β (d).Since the gap is constant as a function of ∆ and linear in β when approaching the critical point from the dimerized phase, Z * displays a linear behaviour in ∆ (c) and it is inversely proportional to β (d).(e) and (f ) The piezoelectric coefficient diverges when approaching the critical point from the dimerized phase.As shown in the inset, the major contribution is due to the internal relaxation term of Equation (5).The topological nature of the enhancement guarantees its stability.

MFPA∼FIG. 2 .
FIG. 2. (a) Mono-fluorinated polyacetylene (MFPA).(b) Polymethineimine (PMI).In (c) and (d) for MFPA and PMI, respectively, is shown the behaviour of the internal coordinate u0 for different values of the long-range mixing parameter cLR put in the range-separated xc-functional in the DFT calculations.Consistently with the prediction of the model, we observe the behaviour u ≃ |cLR − cLR| 1/2 , with c MFPA LR ≃ 14% and c PMI LR ≃ 26%.In (e) and (f ), for MFPA and PMI respectively, the behaviour of the ∂u0/∂ϵ is reported, along with the values of the effective charge Z * .For each polymer, we chose Z * = Z * C,xx

FIG. 3 .
FIG.3.Comparison of the converse piezoelectric coefficients of MFPA and PMI with respect to PVDF, the current best and most widely used organic piezoelectric.In principle, the mechanism of the morphotropic-topological enhancement allows to outperform the state-of-the-art.

)4
FIG. S1.Electronic bands structure of the model.If Egap = 0, the bands have linear dispersion at the edges of the Brillouin zone, reminding the behaviour of graphene.To describe the properties of conjugated polymers, we assume that only the lower energy bands are occupied with electrons.

− 0 .
FIG. S2.Behaviour of the total energy per unit cellEtot with respect to the fractional coordinate u.If ∆ ≤ ∆c, there are two symmetric minima u ̸ = 0 and the equilibrium structure present a bond length alternation.If ∆ > ∆c, there is only one minimum in u = 0 and the atoms are equidistant.
S53) with x = ∆/ √ ∆ 2 + 4t 2 and ∆ = E B − E A .This allows to define the dipole moment D(d AB ) of the molecule as

1 E gap /t 0 = 1 E
FIG.S3.On the left panel, behaviour of the dipole moment per unit cell P with respect to θ for different values of Egap/t0.As expected, in the small-gap limit (Egap ≲ t0) we observe the linear behaviour P ∝ −θ predicted in Equation (S49).On the panel, the evolution of P is displayed in the 2D parametric (∆, δt)−space, where each point correspond to a physical realization of the model.An adiabatic evolution along a loop containing the origin of the axes -which correspond to a metallic system -pumps out a quantised charge of −2|e|.In the 2D space on the central panel it is possible to appreciate the geometric argument that yields Equation (S51).Being the polarization proportional to the angle θ spanning the parametric (∆, δt)-space along a circumference with radius Egap, as δt ∝ u, it follows that in the undimerized phase Egapdθ ∝ du, hence Z * = ∂P/∂u ∝ ∂θ/∂u ∝ 1/Egap.Finally, it is interesting to visualize the phase transitions of u0(∆) and u0(β), displayed separately in the two right panels, also as paths in the 2D parametric space.Being ∆ and β the independent parameters and δt ∝ βu, one finds that the transition in ∆ from the dimerized (u0 ̸ = 0, hence δt ̸ = 0) to the undimerized (u0 = 0, hence δt = 0) phase, happens on a circumference with fixed Egap, whereas the evolution with respect to β in the dimerized phase corresponds to a vertical line, the undimerized phase corresponding to a single point at fixed ∆ ̸ = 0 and δt = 0. Arrows on the two paths in the 2D space correspond, respectively, to the arrows in the right panels whereas all the dots of the lower right panel correspond to a single point in the 2D space in (∆ = 1.05∆c, δt = 0).
FIG. S4.Dynamical charges of atomsA and B. Tuning the ionic/covalent character of the bond through the ratio ∆/t allows for a maximum enhancement of ∼ 1.5 times the valence charge of 1|e|, in stark contrast with the enhancement of up to 30 times we find in the chain.
Bond angle [deg.] FIG.S7.View of the undimerized structure of (i) mono-fluorinated polyacetylene (MFPA) and (ii) polymethineimine (PMI).In (iii) and (iv) for MFPA and PMI, respectively, we report the values of the bond length between the atoms of the unit cell, calculated with different long-range mixing parameter cLR, whereas in (v) and (vi) we report the values of the bond angles for MFPA and PMI, respectively.For MFPA, we use CH and CF to label carbon atoms bonded to hydrogen and fluorine, respectively.Bond angles are labeled as follows: α and β denote the CFCHCF ( CNC) and CHCFCH ( NCN) bond angles, while γ and δ label the CHCFF ( NCH) and FCFCH ( HCN) bond angles.For MFPA only, we also denote by θ and φ the CFCHH and HCHCF bond angles.As expected, the phase transition is signalled by the bond-length alternation as well as in the bifurcation of γ and δ bond angles.

C H H 1 C H H 2 C F F 1 C F F 2 FIG
FIG. S8.On the left, view of the structure polyvinylidene fluoride (PVDF).On the right, behaviour of the bond lengths for different values of cLR.As expected, the structure of PVDF depends very weakly on the long-range mixing parameter, this polymer being not conjugated.
FIG.S9.On the left panel, piezoelectric coefficients of MFPA in its isolated polymeric chain configuration are compared to those obtained for packed chains.The 3D structure adopted, displayed on the right panel, is similar to the one adopted for PVDF in Ref.[49].The agreement between the results highlight the robustness of the enhancement mechanism.