This is the accepted version of the article: Cappai A., Antidormi A., Bosin A., Narducci D., Colombo L., Melis C.. Impact of synthetic conditions on the anisotropic thermal conductivity of poly(3,4-ethylenedioxythiophene) (PEDOT): A molecular dynamics investigation. Physical Review Materials, (2020). 4. 035401: - . 10.1103/PhysRevMaterials.4.035401. Available at: https://dx.doi.org/10.1103/PhysRevMaterials.4.035401 1 PHYSICAL REVIEW MATERIALS 00, 005400 (2020) 2 Impact of synthetic conditions on the anisotropic thermal conductivity of 3 poly(3,4-ethylenedioxythiophene) (PEDOT): A molecular dynamics investigation 4 Antonio Cappai ,1 Aleandro Antidormi ,2 Andrea Bosin,1 Dario Narducci ,3 Luciano Colombo ,1 and Claudio Melis 1,* 1 5 Department of Physics, University of Cagliari, Cittadella Universitaria, I-09042 Monserrato (CA), Italy 2 6 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain 3 7 Department of Materials Science, University of Milano-Bicocca, via Cozzi 55, 20125 Milano, Italy 9 (Received 23 December 2019; accepted 10 February 2020; published xxxxxx) 10 In this work we study the effect of different synthetic conditions on thermal transport properties of poly(3,4- 11 ethylenedioxythiophene) (PEDOT) by focusing in particular on the role of proton scavengers. To this aim, 12 different PEDOT samples were generated in silico using a novel computational algorithm based on a combination 13 of first-principles density functional theory and classical molecular dynamics simulations. The corresponding 14 thermal conductivities were then estimated using the approach to equilibrium molecular dynamics methodology. 15 The results show that the initial synthetic conditions strongly affect the corresponding thermal conductivities, 16 which display variations up to a factor of ∼2 depending on the proton scavenger. By decomposing the thermal 17 conductivity tensor along the direction of maximum chain alignment and the corresponding perpendicular 18 directions, we attribute the thermal conductivity differences to the variations in the average polymer chain length 19 λave . A dependence of the thermal conductivity with the polydispersity index was finally observed, suggesting a 20 possible role of intercrystallite chains in enhancing thermal transport properties. By means of the Green-Kubo 21 modal analysis, we eventually characterize the vibrational modes involved in PEDOT thermal transport and 22 investigate how they are related to the thermal conductivity values of the samples. 23 DOI: 10.1103/PhysRevMaterials.00.005400 24 I. INTRODUCTION where T is the average operating temperature. As mentioned, 48 at the state of the art, current values of ZT in organic polymers 49 25 In the ongoing quest for nonfossil energy supplies, ther- are disappointingly low to make possible large-scale energy 50 26 moelectric devices are a promising opportunity. Thermal en- production by organic-based TEGs [3,7,8]. Hence, several 51 27 ergy is in fact produced almost ubiquitously, and in large investigations have been devoted to devise new strategies to 52 28 amounts, in all conventional power supply stations, but a very achieve at least a meaningful increase of ZT [9–12]. 53 29 large fraction (usually referred to as waste heat) is almost It is generally recognized that the thermal and electrical 54 30 completely lost. In this scenario, the discovery of efficient properties of organic polymers are usually strongly dependent 55 31 thermoelectric materials could promote a sizable replacement on their synthesis and on the resulting morphological features 56 32 of conventional devices for energy production with more [13]. This state of affairs justifies the large variability of 57 33 environment-friendly thermoelectric ones [1]. This situation thermal and electrical conductivities reported in the literature. 58 34 has motivated intense experimental and theoretical efforts In particular, the chain length has been identified as one of 59 35 aimed at increasing the efficiency of thermoelectric generators the key factors in determining the overall transport properties 60 36 (TEGs), as well as at searching for new materials with optimal [14]. 61 37 physical features. In this respect, organic conductive polymers Among organic polymers, poly(3,4-ethylenedioxythio- 62 38 are possible candidates as active constituents of micro-TEGs phene) (PEDOT) is especially promising because of its 63 39 [2,3]. In fact, even if in terms of conversion efficiency they electrical properties (having a hole mobility as high as 3 64 40 are not comparable with their inorganic counterparts, organic cm2 V−1 s−1 [15]), air stability [16,17], and low thermal 65 41 polymers are under investigation because of their interesting conductivity. According to several experimental and 66 42 properties such as flexibility and biocompatibility, very useful computational investigations, thermal conductivity values 67 43 in some specific applications where efficiency of conversion have been reported as small as 0.3 W m−1 K−1 for poly(3,4- 68 44 is not the main goal. The efficiency of conversion can be ethylenedioxythiophene):tosylate (PEDOT:tos) samples 69 45 quantified by means of the figure of merit ZT [4], defined by [18–21]. This good feature is paralleled by a simple synthesis 70 46 a combination of the electrical conductivity σ , the Seebeck process which allows large-scale manufacturing [16,22]. 71 47 coefficient S, and the thermal conductivity κ [5,6] as In a previous work [23] we demonstrated how the poly- 72 σ S2 merization conditions, which are dictated by the choice of 73 ZT = T, (1) proton scavengers, greatly affect the micromorphology of the 74 κ resulting PEDOT samples, with efficient proton scavengers 75 leading to an increase of the maximum and average chain 76 lengths. Moreover, we analyzed the impact of the chain- 77 * Corresponding author: cmelis@dsf.unica.it length distribution on the degree of crystallinity. However, 78 2475-9953/2020/00(0)/005400(8) 005400-1 ©2020 American Physical Society ANTONIO CAPPAI et al. PHYSICAL REVIEW MATERIALS 00, 005400 (2020) 79 the actual interplay between the synthetic process and the TABLE I. Polymerization conditions analyzed in this work and 80 resulting thermoelectric properties is still a matter of debate. eigenvalues ε' and ε p1,p2 associated with the three axes of the inertia 81 Motivated by this, in this paper we explore the impact of tensor calculated for each sample. Eigenvalues are expressed in 82 the proton scavengers on thermal transport. By combining nondimensional units. 83 molecular dynamics simulations and the Green-Kubo modal 84 analysis [24] we compute the main component of thermal Sample ID Scavenger A Scavenger B ε' ε p1 ε p2 ↔ 85 conductivity tensor κ . A modal analysis is then performed in pyr Pyridine Pyridine 0.82 0.76 0.42 86 order to identify the impact of morphology on the vibrational pyr+wat Pyridine Water 0.81 0.70 0.49 87 modes involved. tos Tosylate Tosylate 0.83 0.73 0.43 88 We provide evidence that PEDOT thermal conductivity pyr+tos Pyridine Tosylate 0.77 0.73 0.50 89 is strongly affected by the initial synthetic conditions, i.e., tos+wat Tosylate Water 0.78 0.65 0.57 90 the choice of the proton scavengers employed in the poly- wat Water Water 0.70 0.66 0.64 91 merization process. By analyzing the thermal conductivity 92 components along the direction of maximum chain alignment, 93 we attribute such a difference to the variations in the average GAUSSIAN 16 software package [26]. Then, a simulation box 128 94 chain length λave . Furthermore, by estimating the spatial ex- containing N 3,4-ethylenedioxythiophene (EDOT) monomers 129 95 tension of the different vibrational modes involved in PEDOT is generated and aged in an NVT ensemble by means of classi- 130 96 thermal transport, we observe large differences depending on cal molecular dynamics, using the LAMMPS software package 131 97 the proton scavengers used. [27]. Interactions have been described using the AMBER force 132 field [28]. A cutoff radius of 10 Å is used for Lennard-Jones 133 98 II. METHOD interactions, while electrostatic long-range interactions (i.e., 134 exceeding the cutoff radius) were computed using a particle- 135 99 A. Sample preparation particle particle-mesh (pppm) solver with an accuracy of 10−4 136 100 We simulated the entire polymerization process starting kcal mol−1 Å−1 . 137 101 from initial monomers and taking properly into account the in- The formation of a chemical bond between two molecules 138 102 fluence of different combinations of the main reactants used in is described as a Markov-like mechanism: using the $G(n) 139 103 real experimental conditions. The computational protocol em- calculated in the previous step, a bond creation probability 140 104 ployed and the stipulated reaction mechanism are described in p(n) is defined as 141 105 full detail in our previous work [23]. Here we briefly outline p(n) = min(1, exp($G(n)/(kB T ))), (5) 106 the main features of the two-step algorithm. The Gibbs free 107 energy variation $G(n) (where n is the number of monomeric where kB is the Boltzmann constant and T the operating tem- 142 108 units composing the oligomer) is first calculated for each perature, to quantitatively describe how likely is the formation 143 109 polymerization reaction: of a bond between two units closer than a given distance 144 0→∗ threshold (3 Å in our case). If the bond is created, the force 145 $G(n) = $G0 (n) + $G(n) + $G(n)sol , (2) field underlying molecular dynamics (MD) simulations is 146 110 where $G0 (n) is the gas-phase Gibbs free energy of reaction updated by introducing the new bonds, angles, and dihedrals 147 & using an approach similar to the methodology developed by 148 HXn−1 H+ + HXH+ + Am + Bm Gissinger et al. [29]. By developing a specific polymerization 149 !HXn H + AHm+1 + BHm +1 & (3) code applying this protocol, PEDOT samples resulted to be 150 characterized by different chain-length distributions. The av- 151 111 where HXn−1 H and HXn H are polymer chains of length n − 1 erage chain length ranged from 2 up to 19 monomeric units. 152 & 112 and n, respectively, and Am and Bm are the proton scavengers, In order to evaluate the impact of the chain-length distri- 153 113 either neutral (m, m = 0) or anionic (m, m& = −1). $G0→∗ & bution on thermal transport, the protocol described in the next 154 114 is the thermodynamic correction accounting for the transition section has been used. 155 115 from gas-phase state to standard condition [25], which can be 116 calculated as ! " B. Thermal conduction 156 RT $G0→∗ = $nRT ln , (4) The calculation of thermal conductivity in our samples 157 pV ∗ was performed using the approach to equilibrium molecular 158 117 where $n is the difference between the sum of moles of dynamics (AEMD) method [30–32]: a steplike temperature 159 118 products and the sum of moles of reactants, p is the gas-phase profile was initially imposed within the system along a given 160 119 pressure (1 atm in the present case), 1/V ∗ is the solute con- direction, next evolving according to the heat equation 161 120 centration (assumed to be 1 mol L−1 ), T the polymerization ∂T → ↔ → 121 temperature, and R the ideal gas constant. ρc = ∇ ·( κ · ∇ T ), (6) ∂t 122 Finally, $Gsol is the solvation Gibbs free energy account- ↔ 123 ing for solute-solvent interactions for each reactant and prod- where κ is the thermal conductivity tensor, ρ the mass density, 162 124 uct. The scavengers investigated in this work are summarized and c the specific heat. 163 125 in Table I. In order to create the steplike profile, two regions, one 164 126 At this stage a fully ab initio calculation is performed for each half of the sample, were defined and, respectively, 165 127 via density functional theory (DFT), as implemented in the equilibrated at the target temperature T1 = 350 K and T2 = 166 005400-2 IMPACT OF SYNTHETIC CONDITIONS ON THE … PHYSICAL REVIEW MATERIALS 00, 005400 (2020) modal analysis (GKMA) was employed [33]. Relying both on 199 long MD runs and on the diagonalization of the dynamical 200 matrix [34], this methodology represents a valuable tool to 201 investigate thermal properties in systems with broken transla- 202 tion invariance, where the validity of concepts like phononic 203 group velocity and mean free paths is questionable. For this 204 reason, GKMA has already been successfully employed in the 205 study of amorphous and paracrystalline materials [33]. 206 In GKMA, each Cartesian component κµν of the thermal 207 ↔ conductivity tensor κ is calculated from the equilibrium 208 current-current autocorrelation function, according the Green- 209 Kubo relation 210 $ ∞ * κµν = (qµ (t ) · qν (0))dt, (8) kB T 2 0 where * is the system volume, kB is the Boltzmann constant, 211 T is the system temperature, and the angular brackets denote 212 FIG. 1. Simulation box for AEMD simulation. The simulations an ensemble average of the product of the heat fluxes µ, ν 213 cell (L side) is divided in a hot region (left, red), thermostatted at Cartesian components qµ,ν (t ) at a given time t. The compo- 214 T1 = 350 K, and a cold region (right, blue) at T2 = 250 K to establish nents are obtained from the vibrational eigenmodes extracted 215 a steplike temperature profile. from the dynamical matrix, defined as [33] 216 1 ∂Flα 167 250 K using a Nosé-Hoover thermostat (Fig. 1). During Dlα, jβ = − √ , (9) 168 the equilibration of the left (right) region, the atoms of the ml m j ∂r jβ 169 opposite right (left) region were kept fixed in order to avoid where ml is the mass of the lth atom, Flα the α Cartesian 217 170 any heat exchange during the thermostatting period. After component of the force acting on the lth atom, and r jβ the 218 171 the initial equilibration, the following unconstrained dynamics β Cartesian component of the displacement of the jth atom. 219 172 was generated in an NVE ensemble by recording the average The eigenmodes of the system, extracted from the diago- 220 173 temperatures (T1i (t )) and (T2i (t )) of both regions. Since the nalization of the dynamical matrix given in Eq. (9), yield a set 221 174 samples display no structural order, the initial temperature el,s of 3N vectors (with N the number of atoms in the sample). 222 175 profile was separately imposed along the i = x, y, z directions, For each s eigenmode, a heat flux operator Qs (t ) is defined 223 176 corresponding to the orthogonal edges of the periodically according to 224 177 repeated cubic simulation box (L side). The evolution of the average temperature difference N % N & 1 # # 178 179 $T i (t ) = (T2i (t )) − (T1i (t )) between the two regions of the Qs (t ) = El ẋl (s, t ) + (Flk ẋl (s, t )rlk ) , (10) 180 system was fitted onto the solution of Eq. (6), which reads * l=1 k=1 #+∞ ! " where * is the total volume of the sample, El is the total 8$T i (0) κii π 2 2 225 $T i (t ) = exp − (2n + 1) t , (7) energy of the lth atom, Flk is the force acting on the lth atom 226 π 2 (2n + 1)2 ρcL 2 n=0 when the kth atom is displaced, and rlk is the distance vector 227 181 by using the fact that off-diagonal components of the thermal between the lth and kth atoms. Each Cartesian component µ 228 182 ↔ conductivity tensor κ were found to be sufficiently small of the heat flux operator Qµ,s (t ) is related to the corresponding 229 183 to be negligible. Because of this, only diagonal components qµ (t ) heat flux in Eq. (8) by the relation 230 184 κii enter in Eq. (7). We remark that in the fitting procedure 3N # 185 the sum was truncated at n = 20. We checked that further qµ (t ) = Qµ,s (t ). (11) 186 increases of n do not affect κii . s=1 187 Simulations were carried out on cubic samples of 118 Å 188 side in periodic boundary conditions along each direction: The contribution of the sth mode to the velocity of the lth 231 189 each sample contains as many as ∼105 atoms. Equations of atom is represented by ẋl (s, t ) and it is obtained by projecting 232 190 motion have been integrated by the velocity Verlet algorithm the total velocity vl of the lth atom onto the eigenvectors of 233 191 with a time step of 0.5 fs. Since the temperature-difference the matrix given in Eq. (9) according to 234 192 profiles appear to be very noisy when the average temperature xl (s, t ) = (vl · el,s )el,s . (12) 193 difference between the regions is below ∼5 K, an ensemble 194 average (ten trajectories for each curve) was performed for During an NVE run, Qs (t ) is calculated for a simulation time 235 195 each sample by varying the initial equilibration periods. as long as 200 ps; the overall thermal conductivity κ, defined 236 (as is common in the literature [33]) in terms of the trace of 237 196 C. Modal analysis the thermal conductivity tensor, 238 197 In order to unveil the effect of the chain-length distribution 1 ↔ κxx + κyy + κzz on the thermal properties of the samples, the Green-Kubo κ= Tr( κ ) = , (13) 198 3 3 005400-3 ANTONIO CAPPAI et al. PHYSICAL REVIEW MATERIALS 00, 005400 (2020) Pyridine-based samples Tosylate-based samples Water-based samples 120 120 120 <ΔTx(t)> <ΔTx(t)> <ΔTx(t)> <ΔTy(t)> <ΔTy(t)> <ΔTy(t)> <ΔTz(t)> <ΔTz(t)> <ΔTz(t)> 100 100 100 80 80 80 <ΔTi> (K) <ΔTi> (K) <ΔTi> (K) 60 60 60 40 40 40 20 20 20 0 0 0 0 20 40 60 80 100 120 140 0 20 40 60 80 100 120 140 0 20 40 60 80 100 120 140 t (ps) t (ps) t (ps) FIG. 2. Time evolution of the average ($T i (t )) along the three directions i = x, y, z of a pyr (left), tos (center), and wat sample (right): smooth lines represent the result of the fitting on simulated data. 239 can then be computed within the the Green-Kubo formalism, given mode. It is calculated as 266 240 referring to Eqs. (8) and (11), as *N 2 l=1 nC,l el,s $ ∞ '# 3N #3N ( n̄s = *N 2 , (17) l=1 el,s * κ= Qs (t ) · Qs& (0) dt. (14) 3kB T 2 0 s=1 s& =1 where nC,l is the coordination number of the lth atom, and e2l,s 267 241 The contribution of each mode can be then explicitly calcu- is the squared modulus of its displacement vector. 268 242 lated by relying on the total heat flux operator in Eq. (11) 243 evaluated at t = 0 in Eq. (14), which can be recast in the form III. RESULTS AND DISCUSSIONS 269 3N $ 3N A. Dependence of thermal conductivity κ 270 * # ∞ # κ= (Qs (t ) · q(0))dt = κs , (15) from the adopted scavenger 271 3kB T 2 s=1 0 s=1 Figure 2 shows the time evolution of the average tempera- 272 244 where κs is the contribution of the sth mode to the thermal ture difference ($T i (t )) during the simulation in the case of 273 245 conductivity. pyr (left), tos (center), and wat (right) samples. 274 246 In order to investigate the number of monomers involved For pyr and tos samples, a different ($T i (t )) evolution 275 247 in each vibrational sth mode, a participation ratio (PR) was depending on the i direction is observed, while in the case 276 248 defined as of wat samples the evolution is isotropic. We argue that the 277 ) *N 2 +2 actual choice of the proton scavenger for the polymerization 278 1 l=1 el,s process largely affects the heat propagation, thus generating 279 PR = *N 4 , (16) anisotropy which we attributed to the fact that pyr and tos sam- N l=1 el,s 280 ples shows significantly higher crystalline volume fractions 281 249 where el,s are derived from the diagonalization of the dynam- with respect to wat samples [23], as well as higher average 282 250 ical matrix. chain lengths λave . 283 251 The definition of the PR makes it possible to determine the In order to better clarify the observed anisotropies occur- 284 252 spatial extension of each of the s modes involved in thermal ring for relatively high λave values, we estimated the compo- 285 ↔ 253 transport. In fact, if the eigenmode analyzed involves a large nents of the thermal conductivity tensor κ along the direction 286 254 group of atoms, i.e., if the spatial extension of the eigenmode of maximum chain alignment (κ' component) and along the 287 255 is large, PR ∼1, while in the case of localized modes the PR two corresponding orthogonal directions (κ p1,p2 components). 288 256 is very close to zero and decays as a function of the number They were obtained as follows. 289 257 of atoms, N, involved [PR = O(N −1 )] [35]. Starting from this Starting from a polymerized sample, we defined an orien- 290 258 different behavior, Allen et al. [35] proposed, by studying the tation vector for each unit as the direction identified by the 291 259 paradigmatic case of amorphous silicon, a classification of position of two specific carbon atoms of the thiophenic ring 292 260 extended modes in diffusons, for which the PR quickly decays [Fig. 3(a)]. These M normalized vectors were then translated 293 261 to zero, and propagons, characterized by a larger and almost in order to point out from a specific origin [Fig. 3(b)]. This 294 262 N-insensitive value of PR. representation is equivalent to M points within a sphere [yel- 295 263 Finally, in order to better characterize the vibrational low in Fig. 3(b)] of unitary radius. In systems characterized 296 264 modes, we define [33] the average coordination number as the by long-range order, it is expected to observe all such points 297 265 weighted number of bonds with first neighbors involved in a concentrated around two poles of the sphere, corresponding 298 005400-4 IMPACT OF SYNTHETIC CONDITIONS ON THE … PHYSICAL REVIEW MATERIALS 00, 005400 (2020) (a) (b) (c) FIG. 3. Schematic representation of the procedure used to estimate the direction of maximum chain alignment: (a) polymeric configuration in real space: the black arrows represents the specific C-C orientational vectors; (b) polymeric configuration in the orientational vector space where unitary fictitious masses are attributed to each point lying on the unitary sphere surface; and (c) representation of the inertia tensor in the orientational vector space: the blue arrow represents the direction of maximum alignment. 299 to the direction of maximum chain alignment. In amorphous conductivity component along these three main directions. 321 300 structure, instead, the distributions of these points within the The off-diagonal entries resulted in being negligible with 322 301 sphere must be uniform, since no special directions are present respect to the diagonal κ' , κ p1 , κ p2 components, which are 323 302 in the sample. plotted for each sample in Fig. 4. The dashed band represents 324 303 In order to determine the direction along which the largest the κ previously reported for PEDOT:tos samples [18–21]. 325 304 number of chains is aligned, a fictional unitary mass was The results provide robust evidence of the role of chain 326 305 attributed to each of the M points and all the components of length in determining the magnitude of the resulting thermal 327 306 the corresponding inertia tensor were computed. The mean conductivity: it is in fact shown that samples characterized 328 307 direction of maximum alignment can then be found by esti- by longer chains (pyr and pyr+tos samples) exhibit higher 329 308 mating the eigenvalues of the real non-negative corresponding in-chain thermal conductivities compared to wat samples, in 330 309 density matrix and looking for its column vector associated which κ' reaches its minimum, resulting in being indistin- 331 310 with the lowest eigenvalue [Fig. 3(c)]. guishable from the two κ p1 , κ p2 components. 332 311 The analysis of the eigenvalues (Table I) associated with Referring to Fig. 4, we observed a good agreement between 333 312 each density map clearly suggests the existence of a long- the results for κ' obtained from our calculations, κ' = (0.29 ± 334 313 range order in pyr samples. At variance, wat samples appear 0.01) Wm−1 K−1 , and the values reported in the literature 335 314 to be mainly amorphous. This observation is consistent with [21], κ' = (0.286 ± 0.001) Wm−1 K−1 , for PEDOT:tos sam- 336 315 the results reported in Fig. 2 and previously [23], where we ples. 337 316 provided evidence about the increase in crystallinity of pyr It is interesting to observe that even if out-of-chain com- 338 317 samples compared to wat samples. ponents are negligible with respect to the κ' component 339 318 After the calculations of the three axes of the inertia in high anisotropic samples, a slight reduction is observed 340 319 ellipsoid, the inertia and thermal conductivity tensors were moving from pyr samples to the wat ones. We argue that in 341 320 simultaneously diagonalized in order to obtain the thermal amorphouslike samples the average distance between chains 342 is greater than the corresponding distance in more crystalline 343 samples, resulting in a decrease of the interaction which, 344 0.4 Thermal conductivity (W K m ) -1 -1 0.35 0.4 κ|| 0.3 0.35 κp 0.25 κ|| 0.3 κ(WK-1m-1) κp1 0.2 κp2 0.25 0.15 0.2 0.1 0.15 0.05 0.1 py py to py to w 0.05 at s s+ r r+ r+ 0 2 4 6 8 10 12 14 16 18 20 w to w at at s λave(m.u.) FIG. 4. Thermal conductivity components along the direction of maximum alignment, κ' , and the two orthogonal directions, κ p1,p2 . FIG. 5. Average κ' and κ p values of the samples analyzed plotted The dashed band displays the range of experimental κ previously as a function of average chain length λ (in monomeric units): κ p is reported for PEDOT:tos samples [18–21]. defined as the average value of κ p1 and κ p2 . 005400-5 ANTONIO CAPPAI et al. PHYSICAL REVIEW MATERIALS 00, 005400 (2020) 0.4 30 Pyridine Water 0.35 Weighted coordination number 25 0.3 κ|| (WK-1m-1) 20 0.25 15 0.2 10 0.15 0.1 5 1 1.1 1.2 1.3 1.4 1.5 1.6 PDI 0 0 20 40 60 80 100 FIG. 6. Average κ' values of the samples analyzed plotted as a Frequency(THz) function of PDI. FIG. 8. Average monomeric units involved in vibrational modes as a function of the oscillation frequency for wat (purple) and pyr 345 in turn, is detrimental for thermal transport. As observed in (green) samples. 346 Ref. [23] an increase in crystallite fraction from 0.8% to 347 49% is related to a shift of the peaks of the sulfur-sulfur 348 radial distribution function of (0.06 ± 0.01) nm, suggesting a length on the overall thermal conductivity. In particular, along 371 349 reduction of the average interchain distance. the direction of maximum alignment number the increase of 372 350 Figure 5 shows the average κ' and κ p values of the analyzed κ is almost monotonic. 373 351 samples as a function of λave . A sharp increase of κ' as a 352 function of the average chain length can be observed while B. Modal analysis results 374 353 the average κ p shows a much weaker dependence on λave . We 354 further observe that samples characterized by longer chains In order to provide additional characterization of the im- 375 355 (pyr and pyr+tos samples) exhibit higher parallel thermal pact of morphological features on the vibrational modes (i.e., 376 356 conductivities compared to wat samples, in which κ' reaches heat carriers) involved, we performed a Green-Kubo modal 377 357 its minimum, resulting in being indistinguishable from the two analysis. In Fig. 8 the average weighted coordination number 378 358 perpendicular components. ns [see Eq. (17)] for pyr and wat samples is plotted as a 379 359 In addition, κ' depends also on the polydispersity index function of the frequency ν of the vibrational eigenmodes. 380 360 (PDI) (Fig. 6), suggesting a possible role for the intercrys- Simulations indicate that the number of monomers in- 381 361 tallite chain in enhancing the thermal transport properties volved in each vibrational mode in pyr samples is significantly 382 362 of PEDOT, in strict analogy with the effects on electronic higher compared to that of wat samples. This is quite relevant 383 363 transport observed by Noriega et al. [14] and Kang and Snyder in explaining the higher thermal conductivity corresponding 384 364 [36] in their works. to pyr samples, since the number of units involved in each 385 365 The strong dependence of thermal conductivity on the vibrational mode is a direct measure of the spatial extent of the 386 366 chain length is finally confirmed by estimating the thermal mode. We conclude that in pyr samples even low-frequency 387 367 conductivity for different samples corresponding to different vibrational modes are characterized by an overall extension 388 368 stages of the polymerization procedure (in the case of pyr of about 12–15 monomeric units. Therefore, they are more 389 369 samples), in which λave increases with time. The result, shown effective in the thermal energy propagation if compared to the 390 370 in Fig. 7, clearly shows the impact of increasing the chain case of wat samples, in which the extent of each vibrational 391 0.7 0.28 Pyridine κ|| Water 0.26 κp 0.6 0.24 0.5 0.22 κ(WK m ) -1 0.4 0.2 -1 PR 0.18 0.3 0.16 0.2 0.14 0.12 0.1 0.1 1 2 3 4 5 6 7 8 9 0 0 5 10 15 20 5 # 10 MD steps ν (THz) FIG. 7. Evolution of thermal conductivity components of the pyr FIG. 9. Estimated PR as a function of the frequency for wat sample, κ' and κ p , during the polymerization process. (purple) and pyr (green) PEDOT samples. 005400-6 IMPACT OF SYNTHETIC CONDITIONS ON THE … PHYSICAL REVIEW MATERIALS 00, 005400 (2020) 1.3 increase in the number of modes related with a pure diffu- 423 Pyridine 1.2 Water sion of the thermal energy (i.e., diffusons: vibrational modes 424 immediately scattered and lacking a meaningful mean free 425 1.1 path [33,35,37]) for which the participation ratio PR quickly 426 1 decays [see Eq. (16)]. 427 Normalized κ 0.9 0.8 IV. CONCLUSIONS 428 0.7 In this work we have investigated the effect of different 429 0.6 synthetic procedures on PEDOT thermal transport properties. 430 0.5 To this aim we generated realistic PEDOT samples by means 431 of a novel computational algorithm based on a combination of 432 0.4 first-principles density functional theory and classical molecu- 433 0.3 lar dynamics simulations. We then estimated their correspond- 434 0 2 4 6 8 10 ing thermal conductivities by using the AEMD methodology. 435 Frequency(THz) We observe a clean dependence of PEDOT thermal trans- 436 FIG. 10. Normalized cumulative conductivity as function of the port properties on the initial synthetic conditions, i.e., the 437 frequency for wat (purple) and pyr (green) samples. choice of the proton scavengers employed in the polymer- 438 ization process. By decomposing the thermal conductivity 439 parallel and perpendicular components with respect to the 440 392 mode is limited to two to three monomeric units in all the direction of maximum chain alignment, we prove that the 441 393 range of frequencies investigated. main contribution to the overall κ is in the chain direction, 442 394 By further analyzing the corresponding participation ratio while thermal conductivities in the perpendicular directions 443 395 as defined in Eq. (16) and reported in Fig. 9, we observe that are significantly lower and less influenced by the average 444 396 in pyr samples the low-frequency modes are more spatially chain length. 445 397 extended compared to the case of pure wat samples. The We further characterized the vibrational modes involved in 446 398 participation ratio decay is in fact much slower than in the case PEDOT thermal transport in terms of their spatial extension by 447 399 of wat samples, determining a larger efficiency in propagating computing the average number of monomers involved in each 448 400 thermal energy. In the case of wat samples, on the other vibrational eigenmode and the contributions of each mode to 449 401 hand, the faster decay of PR indicates that such vibrational the thermal conductivity. We remark that, while in the case of 450 402 eigenmodes are spatially more confined. This property is pyr samples the spatial extension of the modes is relatively 451 403 detrimental in terms of thermal transport since eigenmodes high even at low frequencies (0–2 THz), wat samples are 452 404 above the frequency threshold of ∼5 THz cannot efficiently characterized by vibrational modes involving few monomeric 453 405 sustain extended heat currents along all the sample (Fig. 10). units (two or three). Such a difference can be directly related 454 406 Our picture is confirmed by Fig. 10 in which the nor- to the different thermal conductivity values: in particular the 455 407 malized cumulative thermal conductivity κs [defined as in higher thermal conductivity value estimated for pyr samples 456 408 Eq. (15)] is shown as a function of the frequency ν of the is directly caused by the higher chain lengths resulting in a 457 409 vibrational modes in the range 0–10 THz. In both cases de- higher delocalization of the corresponding vibrational modes. 458 410 picted in Fig. 10 the cumulative thermal conductivity quickly The present scenario unveils the critical influence of proton 459 411 converges to its final value at around 2 THz; however, the scavengers in determining a high degree of anisotropy in PE- 460 412 approach is different: while in the case of pyr samples the con- DOT thermal transport by affecting the overall average chain 461 413 vergence value is reached almost monotonically, in the case length and leading to a spatial limitation of low-frequency 462 414 of wat samples, after an initial sharp increase, negative con- vibrational modes crucially responsible for thermal conduc- 463 415 tributions are present in the range 1–2 THz, leading to a de- tivity. 464 416 crease of the cumulative thermal conductivity. This behavior 417 is not in contrast with the definition of κs (ν) since the terms 418 involved in Eq. (15) can also be nonpositive, revealing the ACKNOWLEDGMENT 465 419 presence of a sharp decrease in the number of modes involved 420 in heat propagation (the so-called propagons: low-frequency A.C. gratefully acknowledges Regione Autonoma della 466 421 modes extending for several monomeric units before scatter- Sardegna for financial support under project “Programma 467 422 ing [33,35,37]). At the same time, we observe a corresponding Operativo Regionale FSE 2014-2020”. 468 [1] J. D. König, K. Bartholomé, H. Böttner, D. Jänsch, M. Klein [3] D. Beretta, N. Neophytou, J. M. Hodges, M. G. Kanatzidis, Altstedde, M. Köhne, J. Nurnus, A. Roch, K. Tarantik, and G. D. Narducci, M. Martin-Gonzalez, M. Beekman, B. Balke, G. Hirn, BINE-Themeninfo I/2016, 1 (2016). Cerretti, W. Tremel et al., Mater. Sci. Eng. R 138, 210 (2019). [2] J. L. Blackburn, A. J. Ferguson, C. Cho, and J. C. Grunlan, Adv. [4] M. V. Vedernikov and E. K. Iordanishvili, in 17th Mater. 30, 1704386 (2018). International Conference on Thermoelectrics: Proceedings 005400-7 ANTONIO CAPPAI et al. PHYSICAL REVIEW MATERIALS 00, 005400 (2020) ICT98, Nagoya, Japan (IEEE, New York, 1998), [26] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, pp. 37–42. M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. [5] Z. U. Khan, O. Bubnova, M. J. Jafari, R. Brooke, X. Liu, Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, R. Gabrielsson, T. Ederth, D. R. Evans, J. W. Andreasen, M. J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Fahlman et al., J. Mater. Chem. C 3, 10616 (2015). Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. [6] K. Wijeratne, M. Vagin, R. Brooke, and X. Crispin, J. Mater. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, Chem. A 5, 19619 (2017). B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. [7] C. B. Vining, Nat. Mater. 8, 83 (2009). Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, [8] A. Shakouri, Annu. Rev. Mater. Res. 41, 399 (2011). M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, [9] D. A. Mengistie, C.-H. Chen, K. M. Boopathi, F. W. Pranoto, L.- T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. J. Li, and C.-W. Chu, ACS Appl. Mater. Interfaces 7, 94 (2015). Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, [10] F. Kong, C. Liu, H. Song, J. Xu, Y. Huang, H. Zhu, and J. Wang, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, Synth. Met. 185, 31 (2013). V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. [11] Z. Zhu, C. Liu, F. Jiang, J. Xu, and E. Liu, Synth. Met. 225, 31 Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. (2017). Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. [12] M. Culebras, C. Gómez, and A. Cantarero, J. Mater. Chem. A Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, 2, 10109 (2014). J. B. Foresman, and D. J. Fox, GAUSSIAN16 Revision C.01 [13] A. Crnjar, C. Melis, and L. Colombo, Phys. Rev. Mater. 2, (gaussian Inc., Wallingford, CT, 2016). 015603 (2018). [27] S. Plimpton, J. Comput. Phys. 117, 1 (1995). [14] R. Noriega, J. Rivnay, K. Vandewal, F. P. Koch, N. Stingelin, P. [28] D. A. Case, I. Y. Ben-Shalom, S. R. Brozell, D. S. Cerutti, T. E. Smith, M. F. Toney, and A. Salleo, Nat. Mater. 12, 1038 (2013). Cheatham, III, V. W. D. Cruzeiro, T. A. Darden, R. E. Duke, D. [15] M. Neophytou, J. Griffiths, J. Fraser, M. Kirkus, H. Chen, Ghoreishi, M. K. Gilson, H. Gohlke, A. W. Goetz, D. Greene, C. B. Nielsen, and I. McCulloch, J. Mater. Chem. C 5, 4940 R. Harris, N. Homeyer, S. Izadi, A. Kovalenko, T. Kurtzman, (2017). T. S. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. [16] A. Elschner, S. Kirchmeyer, W. Lovenich, U. Merker, and K. Luo, D. J. Mermelstein, K. M. Merz, Y. Miao, G. Monard, C. Reuter, PEDOT: Principles and Applications of an Intrinsically Nguyen, H. Nguyen, I. Omelyan, A. Onufriev, F. Pan, R. Qi, Conductive Polymer (CRC Press, Boca Raton, FL, 2010). D. R. Roe, A. Roitberg, C. Sagui, S. Schott-Verdugo, J. Shen, [17] K. Sun, S. Zhang, P. Li, Y. Xia, X. Zhang, D. Du, F. H. C. L. Simmerling, J. Smith, R. Salomon-Ferrer, J. Swails, R. C. Isikgor, and J. Ouyang, J. Mater. Sci.: Mater. Electron. 26, 4438 Walker, J. Wang, H. Wei, R. M. Wolf, X. Wu, L. Xiao, D. M. (2015). York, and P. A. Kollman, AMBER (University of California, San [18] L. Cigarini, A. Ruini, A. Catellani, and A. Calzolari, Phys. Francisco, 2018). Chem. Chem. Phys. 20, 5021 (2018). [29] J. R. Gissinger, B. D. Jensen, and K. E. Wise, Polymer 128, 211 [19] O. Bubnova, Z. U. Khan, A. Malti, S. Braun, M. Fahlman, M. (2017). Berggren, and X. Crispin, Nat. Mater. 10, 429 (2011). [30] E. Lampin, P. L. Palla, P.-A. Francioso, and F. Cleri, J. Appl. [20] H. Kwok, J. Electron. Mater. 42, 355 (2013). Phys. 114, 033525 (2013). [21] K. Zhang, J. Qiu, and S. Wang, Nanoscale 8, 8033 (2016). [31] C. Melis, R. Dettori, S. Vandermeulen, and L. Colombo, Eur. [22] W. Hong, Y. Xu, G. Lu, C. Li, and G. Shi, Electrochem. Phys. J. B 87, 96 (2014). Commun. 10, 1555 (2008). [32] G. Fugallo and L. Colombo, Phys. Scr. 93, 043002 (2018). [23] A. Cappai, A. Antidormi, A. Bosin, D. Galliani, D. Narducci, [33] A. Antidormi, X. Cartoixà, and L. Colombo, Phys. Rev. Mater. and C. Melis, Phys. Chem. Chem. Phys. 21, 8580 (2019). 2, 056001 (2018). [24] H. R. Seyf, K. Gordiz, F. DeAngelis, and A. Henry, J. Appl. [34] W. Lv and A. Henry, New J. Phys. 18, 013028 (2016). Phys. 125, 081101 (2019). [35] P. B. Allen, J. L. Feldman, J. Fabian, and F. Wooten, Philos. [25] G. C. Shields and P. G. Seybold, Computational Approaches Mag. B 79, 1715 (1999). for the Prediction of pKa Values (CRC Press, Boca Raton, FL, [36] S. Dongmin Kang and G. Snyder, Nat. Mater. 16, 252 (2017). 2013). [37] H. R. Seyf and A. Henry, J. Appl. Phys. 120, 025101 (2016). 005400-8