Relationship between local coordinates and thermal conductivity in amorphous carbon

Amorphous carbon is a carbon allotrope, characterized by a complex network structure.^1–3 Depending on the fabrication conditions, such as temperature and deposition rate, various microstructures can be obtained from high-density and sp³ bond-rich structures to low-density and sp² bond-rich structures.^4–9 This complexity translates into characteristic mechanical and chemical properties that are advantageous for numerous applications, ranging from protective coating of magnetic storage devices to anodes of sodium-ion batteries.^1,10,11 One of the most important physical properties is thermal conductivity. It has been reported that the thermal conductivity of amorphous carbon spans a wide range of values, from 0.2 W/m K (which is close to that of polymers) to 3.5 W/m K.^12–16 Although it is believed that structural variations drive the thermal conductivity variation, the impact of the sp²/sp³ ratio and relevant atomic configurations on thermal conductivity remains elusive.

Theoretical analysis is useful for addressing the above issues and for exploring the relationship between carbon nanostructures and its thermal conductivity. The thermal conductivity of amorphous carbon has been theoretically simulated using classical molecular dynamics (MD).^17–19 It has been reported that classical MD tends to overestimate the thermal conductivity, compared with experimentally obtained values. One reason for this is the lack of the quantum effect in classical MD.¹⁹ Another problem is the accuracy of the empirical potentials used in the classical MD. Typical interatomic potentials for carbon (Tersoff,²⁰ REBO,²¹ and AIREBO²²) yield amorphous structures that differ in the proportion of sp²/sp³ bonds and mechanical properties.²³ 

By using the ab initio calculations, we can avoid the problem of empirical potentials. Previous studies have analyzed the sp²/sp³ bond ratio and its effect on the electronic and vibrational states of amorphous carbon for various densities by melt-quench simulation using ab initio molecular dynamics (AIMD).^24–28 An evaluation of the interatomic force by ab initio calculations is expected to improve the accuracy of thermal conductivity calculations; however, ab initio simulations of the thermal conductivity of amorphous carbon have not yet been performed.

One reason for the lack of these results is the high computational cost of ab initio calculations. Thermal conductivity calculations using the Green–Kubo method require MD calculations longer than several hundreds of picoseconds for obtaining a sufficient number of ensembles for statistical averaging. Heat flux calculations by non-equilibrium MD require a large system for introducing a reasonable temperature gradient. It is challenging to complete these simulations in a reasonable time using ab initio calculations, which are much more expensive than classical MD calculations.

In this study, we performed ab initio thermal conductivity analysis of amorphous carbon using another Kubo formula-based method, the Allen–Feldman theory.^29–31 We report a proportional relationship between the density of an amorphous structure and its thermal conductivity. We also report that the density strongly correlates with the sp/sp²/sp³ bond ratio. Moreover, topological structural properties, such as the ring structure of carbon atoms, also correlate with the density. Using persistent homology analysis, the correlation between these structural features and density can be quantified, providing a predictive model of thermal conductivity.

In this study, we generated the amorphous carbon structures and determined the force constants by ab initio calculations based on the density functional theory using VASP.^32–35 The PAW^36,37 pseudopotential was used for C atoms. The exchange-correlation functional was described using the local density approximation (LDA). The k-space was sampled only at the Г point. We sampled the density between 2.0 and 3.4 g/cm³ in steps of 0.2 by changing the size of the cubic cell containing 216 atoms. We prepared two data sets, each of which included eight samples with different densities. We set the initial structure to the diamond lattice for data set 1 and to the cubic lattice for data set 2. The melt-quench procedure was the same for the two data sets. All the time steps were 1 fs, and the ensemble was set to NVT with a Nosé–Hoover thermostat.^38,39 First, liquid-state carbon was simulated by increasing the temperature from 300 to 9000 K in 2 ps and by equilibrating at 9000 K for 3 ps. Subsequently, the resultant liquid was quenched to 5000 K in 2 ps and was equilibrated at 5000 K for 2 ps. The amorphous model was obtained by cooling the sample from 5000 K to 300 K in 2 ps. We performed another 2-ps-long equilibration after this melt-quench procedure and sampled the last snapshot as the structural model of amorphous carbon. The structure was relaxed before the thermal conductivity calculation until the force on the atoms dropped below 0.01 eV/Å.

Heat conduction in crystals is described by the phonon gas model, in which quanta of collective vibrations (phonons) carry heat. By combining the phonon gas model and the Boltzmann transport theory, we can calculate the thermal conductivity, including the effects of phonon scattering by other phonons, impurities, and boundaries. However, for highly disordered amorphous materials, the assumption of the heat transport by phonons with well-defined group velocities is no longer valid. The Allen–Feldman theory enables the thermal conductivity computation in such cases. This theory is briefly explained below. For a detailed description, see Refs. 29 and 30.

Based on the Kubo formula, the diagonal component of the thermal conductivity tensor is defined as

In the above, $μ=x,y,z$⁠; V is the system’s volume; T is the temperature; $β=1/kBT$⁠; and $Sμ$ is the $μ$ component of the heat flux operator $S$⁠. In the harmonic approximation, the heat flux operator is decomposed into the contributions of respective modes. The decomposed heat flux operator $Smn$ is defined by the eigenvectors of vibrational modes and the dynamical matrix $Dαβij$ as follows:

In the above, $m,n$ are the indices of the vibrational modes; $am†$ $(an)$ corresponds to the generation (annihilation) operator of a vibrational mode; $α,β=x,y,z$⁠; $ωm$ is the vibrational frequency of the mode m; $γ,γ′$ are indices that distinguish the atoms in a supercell; $Rγ$ is the position vector of atom $γ$⁠; and $eα(γ,m)$ is the polarization of the eigenvector of the vibrational mode m at atom $γ$ along the $α$ direction. Note that the definition of the dynamical matrix $Dαβγγ′$ is

where E is the total energy, $mγ$ is the mass of the atom, and $uγα$ is the displacement of the atom.

The thermal conductivity is expressed as

In the above, $Cm(T)=ℏωm2T[−∂nm∂ωm]$ is the specific heat of mode m. To make the correspondence between this equation and the Boltzmann transport theory clear, mode diffusivity $Dm$ is defined as

Using this definition, Eq. (6) can be expressed as

The mode diffusivity measures the contribution of each mode to thermal conductivity. In this study, the dynamical matrix required for calculating the mode diffusivity was obtained from the combination of phonopy⁴⁰ and VASP package. The mode diffusivity and thermal conductivity from the dynamical matrix were calculated using a custom Python package. The code is available on GitHub.⁴¹ 

Persistent homology^42–45 is an emerging technique for topological data analysis.^46,47 It has been used for analyzing the complex atomic structures of SiO₂ glass,^48–51 metallic glass,^48,52,53 and amorphous ice.⁵⁴ The topological characteristics of atomic configurations, such as rings and voids, were extracted by filtration, as schematically shown in Fig. 1.

In the filtration process, spheres were placed at the respective atoms. Subsequently, the radius of each sphere was gradually increased. In the context of persistent homology, the process of the radius increase is often referred to as “time.” With increasing radius, the spheres start to intersect with one another. When the edges form a closed ring, a “cycle” is obtained. The ring becomes fully covered by spheres as the radius further increases. This is interpreted as the cycle converting into another class of topological feature called a “boundary.” The above-mentioned definition of the first persistent homology can be extended to the higher-order persistent homologies. In the second persistent homology, the formation and coverage of polyhedral void during the filtration procedure corresponds to the cycle-to-boundary transition. The topology of data is represented by the pairs of birth and death times at which a cycle appears and is converted into a boundary. The two-dimensional visualization of birth and death time pairs is called a persistence diagram (PD). The inverse analysis technique⁵⁵ provides a characteristic structure that generates a specific birth–death pair in the PD.

We used the HomCloud code⁵⁶ for analyzing the persistent homologies of the amorphous models. PDs were evaluated for the optimized structures used in this study. We focused on the first and second homologies, which represent the characteristics of the rings and voids in the amorphous phase. For principal component analysis (PCA) and ridge regression, the diagrams were converted into vectors using the persistent image method.^57,58 When constructing persistent image vectors, we choose linear weighting for the distance from the diagonal line of the PD. The standard deviation of the Gaussian distribution used in the persistent image method was 0.06. We focused on the square region [0.0, 6.0] _birth× [0.0, 6.0]_death, which was divided into a 300 × 300 mesh.

The structural models with the densities of 2.0, 2.4, 2.8, and 3.2 g/cm³ obtained by the melt-quench simulations are shown in Fig. 2(a). The color of the atoms corresponds to the coordination number, which indicates that the ratio of the coordination number changes with density. The thermal conductivities of the 16 structures with various densities are shown in Fig. 2(b). A proportional relationship between thermal conductivity and density is evident.

In the Allen–Feldman theory, the mode diffusivity in Eq. (7) represents the contribution of each vibrational mode to thermal conductivity. Figure 3 shows that the diffusivity of the vibrational mode in the 200–500 cm⁻¹ range increases as the density increases.

As was pointed out by Allen and Feldman, the vibrational modes of amorphous materials can be categorized into three groups: (1) low-frequency and long-wavelength modes that retain elastic wave characteristics (propagons), (2) modes that are delocalized but have no clear periodicity in their eigenvectors (diffusons), and (3) localized modes (locons). The results of the inverse participation ratio in Fig. 4 show that modes up to 500 cm⁻¹ are not localized. As shown in Fig. 3, the mode diffusivity decays rapidly to zero below 200 cm⁻¹, especially for high-density systems. This occurs owing to the small size of these systems and indicates an underestimation of the contribution of propagons. In general, the density of states of the propagon mode follows the Debye approximation that scales like ω⁻², whereby the mode diffusivity follows a continuous distribution.⁵⁹ Small systems cannot reproduce such densities of states and, thus, underestimate the contribution of propagons. In fact, the thermal conductivity of a-C with a density of 3.0 g/cm³ was experimentally reported to be 2.2–2.7 W/m.K,¹³ but in our calculations, it was only approximately 0.9 W/m K. This was attributed to the underestimation of the contribution of propagons. Because of this underestimation, it was difficult to fully reproduce the experimental values for the present system size, but we believe that the density dependence of diffuson-derived thermal conductivity in the 200–500 cm⁻¹ region was accurately evaluated.

To investigate the origin of this relationship between thermal conductivity and density, we analyzed the characteristics of the structures. We first calculated the distribution of the nearest-neighbor distance and bond angle. The results are shown in Fig. 5. Similar to the previously reported results,^24–28 the mean nearest-neighbor distance increased in proportion to the density, while that of the bond angle decreased.

The percentage of atoms with sp, sp², and sp³ bonds determined from the coordination number of each C atom are summarized in Fig. 6. In Fig. 6, the threshold value for determining the presence or absence of bonds was set to 1.8 Å. We assigned the C atoms with coordination numbers of two, three, and four to sp, sp², and sp³ bonded atoms, respectively. In Fig. 6, the results for all 16 samples from data sets 1 and 2 are plotted. At densities above 2.4 g/cm³, there are almost no sp-bonded atoms, and sp² and sp³ bonds are dominant. The ratio of sp³-bonded atoms is linearly proportional to the density, whereas sp² is inversely proportional. Below 2.4 g/cm³, the change in the ratio of sp² and sp³ bonds becomes less visible and that in the ratio of the sp bonds becomes more abundant with decreasing density.

The relationship between these changes in connectivity and local structures was explored using ring statistics (Fig. 7). We used the definition by Matsumoto et al.⁶⁰ for identifying rings. For both data sets 1 and 2, the number of rings with more than eight vertices increases with decreasing density. These results can be explained by the following mechanism: the rings with five, six, and seven vertices are formed by a mixture of sp² and sp³ bonds, and the increase in the number of rings with more vertices is attributed to the addition of sp bonds. This implies that the change in the topological features induced by the change in the coordination number is a key factor determining the density dependence of thermal conductivity.

The above information about nearest-neighbor distances, bond angles, coordination numbers, and ring statistics are encoded in the persistent homology. The first and second PDs for the lowest- and highest-density samples are shown in Fig. 8. In the first persistence diagram, the region around the birth time ∼0.6 Ų shows a strong density dependence. The value of 0.6 Ų corresponds to the square of the half of the nearest-neighbor distance in amorphous carbon (∼0.77² Ų). In the second persistence diagram, although the distribution is not as clear as that for the first persistence diagrams, the overall distribution tends to be concentrated in the region with smaller birth-death times as the density increases.

In our previous work, we showed that for amorphous Si, the local structure affecting thermal conductivity can be identified from the PD.⁶¹ Therefore, we performed a similar analysis for amorphous carbon. Persistent image vectors were generated from the PDs. We determined the regions in the diagrams that govern the density dependence, using the PCA of the persistent image vectors. The results are shown in Fig. 9. The different-density structures are well separated by the first and second principal components in both the first and second PDs [Fig. 9(a)]. In the first-order (second-order) persistence diagrams, the first principal component explains 68.3% (27.7%) of the data variance, while the second principal component explains 13.5% (12.1%) of the data variance. Figure 9(b) shows the projection of the PCA coefficients onto the PD, revealing regions that characterize different-density structures. The PCA coefficients show that the first homology is characterized by the death time in the region with the birth time of ∼0.6 Ų, which is consistent with the results in Fig. 8. In the second homology, a smaller birth-death time correlates with a higher density structure, as in the PD. This can be intuitively understood by the smaller size of the polyhedron in the amorphous structure at higher densities.

To understand the characteristic death time dependence for the first homology, we determined the local structure corresponding to birth–death pairs in regions with large absolute values of PCA coefficients by a stable volume.⁶² Figure 10 shows a histogram of the number of vertices in stable volume cycles. The number of vertices for the birth–death pair in the red region in Fig. 9(c) was small, while that in the blue region in Fig. 9(c) was large. In addition, the number of birth–death pairs in the blue region itself was very small in the high-density structure. Although the definition of the cycle determined by persistent homology is different from that of the ring structure identified in ring statistics, this result is consistent with the trend observed in ring statistics: large rings increase with decreasing density.

Based on these results, persistent homology is expected to be a useful descriptor for capturing both geometrical and topological characteristics and can be used for predicting physical properties. Indeed, we succeeded in constructing a linear regression model for thermal conductivity, using the strong correlation between the persistent image vector and the amorphous structure. Figure 11 shows the results of the fivefold cross validation of the ridge regression model using the persistent image vector. The model’s predictions for the training and validation data from the five runs are plotted on the same axis. Although the accuracy of the predictions was limited by the amount of data in this study, the average root mean square error (RMSE) of thermal conductivity was approximately 0.05 W/m K. These results indicate that persistent homology is an excellent descriptor that integrates information about the nearest-neighbor distance, bond angle, coordination number, and ring statistic, linking between complex structures and physical properties.

We performed ab initio thermal conductivity analysis of amorphous carbon using the Allen–Feldman theory. The obtained thermal conductivity values showed a proportional relationship with the density of carbon. We demonstrated that the density dependence was owing to the contribution from diffusons in the 200–500 cm⁻¹ range. Note that the contribution of propagons was underestimated in this work, owing to the small size effect, and an accurate evaluation of the density dependence of the propagons’ contribution remains to be addressed.

To understand the density dependence of the diffuson contributions, we analyzed the structures at various densities. Similar to the results reported in previous studies, the mean nearest-neighbor distance increased in proportion to the density, while the mean bond angle decreased. The density also strongly correlated with the fraction of sp, sp², and sp³ bonds, which was reflected in the topological features extracted from ring statistics. This implies that the change in the topological features induced by the change in the coordination number is a key determinant of the density dependence of thermal conductivity.

We demonstrated that the above geometrical and topological information are integrated into persistent homology. The PCA of persistent image vectors generated from PDs identified the characteristic local structures at different densities. The number of vertices was small for the cycles in high-density structures, whereas it was large for low-density structures. This result was consistent with the trend observed for ring statistics. Furthermore, we succeeded in constructing a linear regression model for thermal conductivity using the above strong correlation between the persistent image vector and the amorphous structure. This indicates that persistent homology is a promising descriptor for capturing the correlation between structural and physical properties of complex systems.

This study was supported by JST, PRESTO Grant Nos. JPMJPR17I7, JPMJPR17I5, JPMJPR19I4, JPMJPR1923, JPMJPR2198 and MEXT KAKENHI 21H01816, 19H02544 19H00834, 20H05884, Japan. The calculations were performed using a computer facility at the Research Center for Computational Science (Okazaki, Japan).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

 Emi Minamitani received her bachelor’s, master’s, and Ph.D. degrees from the Department of Applied Physics at Osaka University, Japan in 2005, 2007, and 2010, respectively. She was a postdoctoral researcher at RIKEN from 2011 to 2013 and an assistant professor and lecturer at the University of Tokyo from 2013 to 2015 and 2015 to 2019. She is currently an associate professor at the Department of Theoretical and Computational Molecular Science at the Institute for Molecular Science. Her research focuses on theoretical studies of nanoscale physics, including surface magnetism, electron–phonon interactions, and thermal properties. She received the L’Oréal–UNESCO Japan National Fellowship for Women in Science in 2008 and the Young Scientists’ Prize by the Minister of Education, Culture, Sports, Science and Technology in 2019.

 Takuma Shiga is a lecturer at the Department of Mechanical Engineering, School of Engineering, University of Tokyo. He received his bachelor’s and master’s degrees in Physics in 2008 and 2010 from the Tokyo University of Science and received his Ph.D. degree in Engineering in 2013 from the University of Tokyo. His research focuses on heat conduction in materials from first-principles and thermal control utilizing phononic effects. He is a recipient of the Young Scientist’s Prize by the Minister of Education, Culture, Sports, Science, and Technology and the Academic award of the Heat Transfer Society of Japan.

 Makoto Kashiwagi is an assistant professor at the Department of Chemistry and Biological Science, Aoyama Gakuin University. He received a B. Eng in 2008 from the National Institution for Academic Degrees and Quality Enhancement of Higher Education and received M. Eng and Ph.D. in Engineering in 2010 and 2015 from the Kyushu Institute of Technology. He was a postdoctoral researcher at the University of Tokyo from 2015 to 2018 and a junior researcher at Waseda University from 2018 to 2019. He moved to the Aoyama Gakuin University in 2019. He focuses on the thermal transport phenomena in disordered materials and thermophysical properties of nanostructured materials.

 Ippei Obayashi is a professor at the Cyber-Physical Engineering Research Core (Cypher), Okayama University. He received a B. Sc in 2004, M. Sc in 2006, and Ph.D. in mathematics in 2010 from Kyoto University. He was a postdoctoral researcher in applied mathematics at Kyoto University from 2010 to 2015, an assistant professor at the Advanced Institute for Materials Research (AIMR), Tohoku University from 2015 to 2018, and a research scientist at the Center for Advanced Intelligence Project (AIP), RIKEN, from 2018 to 2021. He moved to Okayama University in 2021. He is an applied mathematician interested in topological data analysis and dynamical systems. Currently, he mainly focuses on the mathematical theory, software, and applications of persistent homology. He is the main developer of HomCloud, a persistent homology-based data analysis software. He received the Young Researchers Prize of the 6th Hiroshi Fujiwara Prize on Mathematical Science for scientific research on persistent homology and development of HomCloud.