Carbon nanotubes and graphene are promising materials for thermal management applications due to their high thermal conductivities. However, their thermal properties are anisotropic, and the radial or cross-plane direction thermal conductivity is low. A 3D Carbon nanotube (CNT)-graphene structure has previously been proposed to address this limitation, and direct molecular dynamics simulations have been used to predict the associated thermal conductivity. In this work, by recognizing that thermal resistance primarily comes from CNT-graphene junctions, a simple network model of thermal transport in pillared graphene structure is developed. Using non-equilibrium molecular dynamics, the resistance across an individual CNT-graphene junction with *sp*^{2} covalent bonds is found to be around $6\xd710\u221211$ m^{2}K/W, which is significantly lower than typical values reported for planar interfaces between dissimilar materials. In contrast, the resistance across a van der Waals junction is about $4\xd710\u22128$ m^{2}K/W. Interestingly, when the CNT pillar length is small, the interfacial resistance of the *sp*^{2} covalent junction is found to decrease as the CNT pillar length decreases, suggesting the presence of coherence effects. To explain this intriguing trend, the junction thermal resistance is decomposed into interfacial region and boundary components, and it is found that while the boundary resistance has little dependence on the pillar length, the interfacial region resistance decreases as the pillar length decreases. This is explained by calculating the local phonon density of states (LDOS) of different regions near the boundary. The LDOS overlap between the interfacial region and the center region of CNT increases as the pillar length decreases, leading to the decrease of interfacial region resistance. The junction resistance *R _{j}* is eventually used in the network model to estimate the effective thermal conductivity, and the results agree well with direct MD simulation data, demonstrating the effectiveness of our model.

## I. INTRODUCTION

Over the past few decades, the size of electronic devices has decreased to nanoscale dimensions, and thermal management of these devices has become a significant issue to ensure high performance and reliability. Materials with high thermal conductivity are needed to effectively dissipate thermal energy. Carbon nanotubes (CNTs)^{1–3} and graphene^{4–7} are promising high thermal conductivity materials, and at the same time, they offer low mass density. However, both suffer from anisotropy in thermal transport. The thermal conductivity of graphene stacks or graphite in the cross-plane direction is two or more orders of magnitude lower than that of the in-plane direction.^{8} CNT bundles also show similar behavior in the radial direction.^{9} Hence, the idea of constructing three-dimensional interconnected structures, with CNT, graphene or both as building blocks, such as pillared graphene,^{10} carbon nanotubes with intramolecular junctions,^{11} and CNT networks,^{12} has been proposed. These structures offer potential utility in thermal management of electronic devices.

The thermal transport through pillared graphene with *sp*^{2} covalent bonds at junctions, as shown in Fig. 1(a), has been modeled recently by Varshney *et al.*^{13} using non-equilibrium molecular dynamics (NEMD). They calculated the thermal conductivity of both in-plane and cross-plane directions to fall in the range of about 1–20 W/(mK). Their results suggest that the hierarchical structures can help to reduce thermal anisotropy. Lee *et al.* used a wave packet method to explore phonon energy transmission at the CNT-graphene junction.^{14} Park and Prakash presented reverse NEMD (RNEMD) simulations of thermal transport in two pillared graphene architectures with different symmetries.^{15} These MD simulations were done on the entire 3D network and can only give the thermal conductivity of one particular architecture at a time. On the other hand, the thermal transport across an individual junction is not well understood.

In this paper, we assume that the thermal resistance of the complex 3D architecture is dominated by the CNT-graphene junction resistance, while the resistances of the graphene sheets and CNT pillars are relatively small in comparison. Therefore, a simple thermal network model is proposed, shown in Fig. 1(b) as an effective and efficient thermal transport design tool for such pillared graphene networks. NEMD simulation is first used to predict the thermal resistance across an individual CNT-graphene junction with respect to different structural parameters. Trends between junction resistance *R _{j}* and geometric parameters are found and are then explained using the local density of states (LDOS) concepts. Finally, the junction resistance

*R*is used in the network model to predict the thermal conductivity of the 3D pillared graphene structure, and the results are compared to direct MD simulations. The thermal conductivity of a network model including CNT resistance is also calculated as a comparison.

_{j}## II. MD SIMULATION METHOD

In the previous predictions of thermal transport in a pillared graphene structure,^{13,15} the thermal conductivity in the cross-plane direction was found to depend strongly on the CNT-graphene junction, whereas the in-plane conductivity was less sensitive. Hence, in this work, we focus on the cross-plane thermal properties. Because of the high thermal conductivity of both CNTs and graphene, the assumption that the thermal resistance in the cross-plane direction mainly comes from junction resistance is made. Therefore, a network model is constructed, as shown in Fig. 1(b), composed of junction resistances in parallel and series. Once the junction resistance is calculated, the thermal network can be used to easily find the total resistance of the structure, from which the thermal conductivity can be derived. This approach can save significant computation time and can give insights to thermal transport in such complex structures.

In order to obtain the junction thermal resistance, the NEMD^{16} is performed on just a small cell of the pillared graphene structure. The small cell is composed of three (6, 6) CNT pillars and two graphene sheets, as shown schematically in Fig. 1(c). The CNT is oriented in the *z* direction, while the graphene sheets are in the *x-y* plane. Three important geometric parameters exist in the simulation domain: the length of CNT pillar between two graphene sheets called pillar length (*PL*), the minimum graphene length between two CNT pillars called minimum inter-pillar distance (*MIPD*),^{13} and the distance from the thermal reservoir to the junction (here we call it “reservoir distance,” or *RD*). Periodic boundary conditions are imposed in the *x* and *y* directions, and fixed boundary conditions are used for the cross-plane (*z*) direction. Interactions between carbon atoms are described by the AIREBO potential,^{13,17,18} and the potential parameters are from Ref. 17. The same structure but with van der Waals interaction at the junction is also investigated. The Lennard-Jones 12–6 potential with $\u03f5=4.7483\xd710\u221222$ J and $\sigma =0.3407$ nm (Ref. 19) is used in this case.

In our simulation, the *PL* varies from 0.4 to 100 nm, the *MIPD* varies from 1 to 5 nm, and the *RD* varies from 1 to 100 nm. The thermal reservoirs (heat source and sink regions) have lengths of 5 nm each, as marked in Fig. 1(c). At first, an NPT ensemble is run for 2 ns to relax the structure and allow it to reach thermal equilibrium at 300 K. Then it is switched to the NVE ensemble. The simulation time step is 0.4 fs. A heat flux of 1 eV/ps is then added to the source and extracted from the sink for 20 ns. After the system reaches steady state, the temperature of each small cell is obtained by averaging over 10^{7} time steps over the last 4 ns. Such temperature profiles are used to obtain the thermal transport results.

## III. RESULTS AND DISCUSSION

### A. Temperature distribution and junction resistance

#### 1. *sp*^{2} covalent junctions

Figure 2 shows a typical temperature profile of the middle CNT pillar and two graphene sheets with *sp*^{2} covalent junctions. The temperature distributions in both CNT and graphene are nearly linear, while at the junction an abrupt temperature jump exists. Hence after linear fitting, the temperature jump $\Delta Tj$ at the junction can be calculated. The junction thermal resistance *R _{j}* is given by

Here, *q* is the heat flow rate. $Ac=8.57\xd710\u221219$ m^{2} is the cross-sectional area of a CNT, approximated as a circular ring of thickness 0.335 nm and of radius (average of inner and outer radius) which equals to the radius of CNT.

The junction resistance *R _{j}* is predicted as a function of

*PL*, and the results are shown in Figure 3(a).

*R*ranges from $4.1\xd710\u221211$ to $7.2\xd710\u221211$ m

_{j}^{2}K/W. As a comparison, we also calculate the diffusive limit of $Rj,d$ by using a junction with very long CNT and large graphene, and the value is $7.5\xd710\u221211$ m

^{2}K/W. These values are very low compared to typical thermal boundary resistances across other dissimilar materials, which are often in the range of $10\u221210$–$10\u22128$ m

^{2}K/W.

^{20–23}This fact can probably be attributed to no mismatch in atomic mass, and small mismatch in bonding and phonon dispersion across the CNT-graphene junction. Also,

*R*for finite

_{j}*PL*is smaller than the diffusive limit $Rj,d$, and it decreases as

*PL*decreases, in contrast to the common case where

*R*is a constant. This behavior indicates that the two junctions are coupled but not independent. Similar coherent behaviors have been observed in superlattices with small periods.

_{j}^{24}It will be discussed in more depth in Section III B. Using the temperature drop of CNT bulk region in Eq. (1), we can calculate the CNT resistance in a similar way. It is typically in the range of $2.5\xd710\u221211$ to $3.3\xd710\u221211$ m

^{2}K/W for $1.4<PL<5.3$ nm, and this weak dependence on

*PL*is due to the ballistic transport in CNT. Since each CNT is associated with two junctions, the CNT resistance is indeed smaller than the total resistance of two junctions. Therefore, our hypothesis that neglecting the CNT resistance for estimation of the network thermal conductivity is reasonable, noting that the network thermal conductivity will be overestimated a bit by doing so. Nevertheless, we will compare the final results with and without CNT resistance in Section III D.

Figure 3(b) shows the dependence of *R _{j}* on

*MIPD*and

*RD*, and the dependencies are very weak, indicating that our calculations are converged on these parameters.

#### 2. van der Waals junctions

A typical temperature profile for a van der Waals junction is shown in Fig. 4, and the calculated *R _{j}* is $4.0\xd710\u22128$ m

^{2}K/W, which is much higher than that of the covalent junction and is the result of weak bonding. This is consistent with previous experimental observations of bonding effects on thermal interfacial transport.

^{25}Also, it has little dependence on

*PL*, indicating the absence of coherence effects. Moreover, this resistance becomes orders of magnitude higher than the CNT bulk resistance, hence in this case neglecting CNT and graphene resistance will not affect the accuracy of our network model.

### B. Decomposition of the junction resistance

#### 1. *sp*^{2} covalent junctions

The dependence of *R _{j}* of

*sp*

^{2}covalent junction on

*PL*is intriguing. Under the well-known acoustic and diffuse mismatch theories, interfacial thermal resistance can be attributed to the mismatch of phonon impedance or number of modes between the two dissimilar materials constituting the interface. Therefore, the temperature distribution is expected to show linear trends in two materials, with an abrupt jump at the boundary. However, Fig. 2 indicates that although the temperature profile in CNT is linear when away from the boundary, it shows a clear nonlinear region near the boundary. In particular, the slope is higher than in the linear region, implying higher local resistance when approaching the boundary. Such nonlinear regions have also been observed for planar interfaces

^{26}and have been attributed to a distinct interfacial influence. That is, the atoms in such interfacial regions interact with atoms on the other side of the boundary, resulting in deviations in the LDOS relative to the bulk DOS.

^{27}In Fig. 2, the “interface” is defined as the region composed of the CNT interfacial region, the boundary, and the graphene interfacial region.

The effect of the interfacial region on thermal interfacial resistance is two-fold. On one hand, the interfacial region introduces additional resistance inside the region relative to conventional mismatch concepts, and we term this additional resistance as the interfacial region resistance. On the other hand, the existence of interfacial regions will make the mismatch right across the boundary smaller than the conventional mismatch models, hence reducing the boundary resistance. In other words, in the conventional theory, the mismatch only exists at the boundary, but now it exists over a transition space that includes the boundary and two interfacial regions. Hence, as shown in Fig. 2, $\Delta Tj$ can be divided into three parts: (1) the temperature drop between CNT temperature linear fitting and the true CNT temperature at the boundary $\Delta TCNT,IR$; (2) the temperature drop between the true CNT and graphene temperatures at the boundary $\Delta TB$; and (3) the temperature drop between graphene temperature and graphene temperature linear fitting at the boundary $\Delta Tgra,IR$. Therefore, from the decomposition of $\Delta Tj$, the corresponding *R _{j}* can also be divided into three parts as shown in Fig. 2: (1) graphene interfacial region resistance $Rgra,IR$ due to the phonon DOS mismatch between bulk graphene and the graphene interfacial region; (2) the boundary resistance $RB$ due to the mismatch between the CNT interfacial region and graphene interfacial region; and (3) the CNT interfacial region resistance $RCNT,IR$ due to the mismatch between the CNT bulk and CNT interfacial region.

^{26}It is obvious that

The temperature jump $\Delta TB$ due to boundary resistance $RB$ is the largest among the three.

When the middle CNT pillar becomes short($1<PL<3$ nm), the temperature distribution in the CNT does not exhibit a linear portion any more, but only nonlinear profiles, as shown in Fig. 5(a). This phenomenon implies that the local phonon DOS never reaches the bulk DOS even at the center of the CNT due to its short length, and hence the two junctions are coupled. When the CNT pillar becomes extremely short (*PL* < 1 nm), the interfacial regions disappear entirely, and the two junctions are directly coupled as shown in Fig. 5(b). The CNT pillar now behaves essentially as one scatterer between the two graphene sheets. Hence, the decomposition of the junction resistance is not shown.

The decomposed boundary and CNT interfacial region resistances for different CNT pillar lengths described above have been calculated, and the results are shown in Fig. 3(a). Because the graphene interfacial region resistance is very low and negligible compared to others, it is not plotted in the figure. Importantly, the boundary resistance does not change much for different pillar lengths, which means the mismatch between the CNT interfacial region and graphene interfacial region is insensitive to *PL*. However, the CNT interfacial resistance exhibits a similar trend as the total junction resistance, i.e., decreasing as the pillar length decreases. This implies that changes in the overall junction resistance with respect to *PL* are largely dictated by the CNT interfacial region resistance.

#### 2. van der Waals junctions

Unlike *sp*^{2} covalent junctions, the temperature profile of the van der Waals junction shows no clear interfacial regions as shown in Fig. 4. The temperature jump at the junction nearly equals to the boundary temperature jump. Using the decomposition approach employed for *sp*^{2} covalent junctions, the boundary resistance $RB$ is found to dominate the junction resistance *R _{j}*, while the interfacial region resistances are negligible. As mentioned above, the interfacial region by definition possesses a modified DOS. For van der Waals junctions, the atoms near the boundary have very weak interactions with atoms on the other side of the junction, and therefore, their local DOS is little modified, in contrast to the

*sp*

^{2}covalent junction.

### C. Phonon density of states

As mentioned above, the three decomposed resistances correspond to different regions of the structure. Here, the LDOS^{28} is calculated at different regions of *sp*^{2} covalent junction structure for three different *PLs*. The phonon density of states is calculated by taking the Fourier transform of the velocity autocorrelation function (VAC)

where

and *t* is the correlation time. The *i*-sum is over the number of atoms in the chosen domain, while the *j*-sum is over different starting time of the correlation. The Fourier transform of the VAC is

and is computed by the MATLAB fast Fourier transform (FFT) function. After the FFT results are normalized, the LDOS of different regions can be assessed.

In our results, it is found that the LDOS of the graphene interfacial region is not sensitive to *PL*. Hence, the graphene interfacial region resistance will be little affected by changes in *PL*.

Figures 6(a) and 6(b) show the LDOS of different regions. When the CNT pillar is long, the LDOS of CNT center region is like the LDOS of CNT bulk, and the LDOS of the CNT interfacial region and that of the CNT center region can be distinguished. However, as *PL* decreases, the phonon spectrum does not fully reach the bulk spectrum at the center of CNT as shown in Fig. 6(a). As a result, the mismatch between the CNT interfacial region and the center region is less than the mismatch between CNT interfacial region and CNT bulk, and the mismatch becomes smaller with decreased *PL*. When the CNT pillar becomes extremely short, the LDOS of the center of CNT pillar develops similar features to these of the CNT interfacial region as shown in Fig. 6(b). As a result the CNT interfacial region resistance $RCNT,IR$ diminishes. Therefore, the junction thermal resistance *R _{j}* will be lower than that in other cases.

To show these arguments more quantitatively, an overlap factor *η* is defined between the LDOS of two different regions named region 1 and region 2 as

where and *D*_{1} and *D*_{2} are the normalized LDOS of region 1 and region 2. The numerator of Eq. (6) means the overlap area between region 1 and region 2, while the denominator means the normalized LDOS area of a region. Hence, a high overlap factor *η* means a small mismatch between two regions.

We have calculated the LDOS overlap factor *η* between CNT interfacial region and center region, as well as *η* between graphene interfacial region and CNT interfacial region as a function of *PL*. *η* between CNT center and CNT interfacial regions determines the CNT interfacial region resistance $RCNT,IR$, and *η* increases with the decrease of *PL* as shown in Fig. 6(c). It indicates that $RCNT,IR$ will decrease when *PL* decreases. Conversely, *η* between the graphene interfacial and CNT interfacial regions, which dictates the boundary resistance $RB$, changes little with the decrease of *PL* as shown in Fig. 6(c). These results are consistent with our decomposition results, and elucidate why as *PL* decreases, $RCNT,IR$ decreases, while $RB$ remains nearly the same. Here, we show the overlap factor *η* for $PL\u22645.2$ nm, because *η* does not change noticeably for longer *PL*.

### D. Network model

After the junction thermal resistance is obtained, the thermal resistance of the entire pillared graphene structure ($R\Sigma $) can be calculated with a network model. The cross-plane thermal conductivity of the structure can be calculated based on Fourier's law

where *A _{n}* is the cross-sectional area of the network, and

*L*is the length of the system in

*z*direction.

Figure 7 provides a thermal conductivity comparison between our network model prediction and prior direct MD simulation results.^{13} Both results with and without the CNT resistance are shown. Our prediction based on the network without CNT resistance agrees well with the MD results, and that with CNT resistance is lower than the direct MD result. This is a bit surprising since the former is supposed to over-predict the thermal conductivity, and the latter to agree with direct MD data. However, our simulation domain is only one cell of pillared graphene, and the junction resistance would be slightly overvalued because of the simulation domain size effect. This factor contributes to a reduction of thermal conductivity, making the results without the CNT resistance matching with direct MD data, while those with CNT resistance are under-estimated. Nevertheless, the network prediction results are reasonable even when the CNT resistance is neglected. However, when *PL* or *MIPD* becomes large, our assumption that the CNT and graphene resistances can be neglected will not be appropriate anymore. Instead these resistances should be included in the network model.

## IV. SUMMARY

Because of the high thermal conductivity of both CNT and graphene, we show that the thermal resistance of pillared graphene structure is dominated by the CNT-graphene junction resistance, which can be confirmed by the large temperature jump at the junction and small temperature fall in CNT and graphene. Therefore, we have proposed to use the network model with junction thermal resistance as building blocks to find the total resistance of the structure. The non-equilibrium molecular dynamics is used to predict the junction thermal resistance, and the resistance of *sp*^{2} covalent junction is obtained to be about $6\xd710\u221211$ m^{2}K/W, which is much lower than typical thermal boundary resistance across other dissimilar materials. In contrast, the resistance across an individual CNT-graphene junction with van der Waals force is about $4\xd710\u22128$ m^{2}K/W. The *sp*^{2} covalent junction resistance is found to be primarily related to the pillar length *PL* of CNT and less to the inter-pillar distance *MIPD* or the reservoir distance *RD*, and it decreases as CNT pillar length *PL* decreases. To explore the underlying mechanism of this behavior, the junction resistance *R _{j}* is decomposed into one boundary resistance and two interfacial region resistances. We observed that the CNT interfacial region resistance $RCNT,IR$ decreases as PL decreases, while the boundary resistance $RB$ barely changes. This implies that change in the overall junction resistance with respect to PL is largely dictated by the CNT interfacial region resistance. We then used the LDOS approach to further investigate this phenomenon and found that the phonon spectra mismatch between the CNT interfacial region and the CNT center region decreases when the PL decreases, while the mismatch between graphene interfacial region and CNT interfacial region changes little, which is consistent with our decomposition result. The junction resistances are then put in series and parallel networks to easily interpret the thermal conductivity of such complicated and hierarchical structure. The thermal conductivity of the pillared graphene structure predicted by the network model agrees well with the direct simulation results. This indicates that our simple network model is effective in estimating thermal transport in hierarchical structures. Our work has revealed the critical role of junction resistance to thermal transport in pillared graphene structures, and the developed network model offers several orders of magnitude higher computation speed.

## ACKNOWLEDGMENTS

The authors would like to acknowledge the support by the Air Force Office of Scientific Research (AFOSR) MURI grant (FA9550-12-1-0037).