Effect of the in-plane aspect ratio of a graphene filler on anisotropic heat conduction in paraffin/graphene composites

Enhancement of polymer thermal conductivity using nanographene fillers and clarification of its molecularscale mechanisms are of great concern in the development of advanced thermal management materials. In the present study, molecular dynamics simulation was employed to theoretically show that the in-plane aspect ratio of a graphene filler can have a significant impact on the effective thermal conductivity of paraffin/graphene composites. Our simulation included multiple graphene fillers aggregated in a paraffin matrix. The effective thermal conductivity of a paraffin/graphene composite, described as a second-rank tensor in the framework of equilibrium molecular dynamics simulation, was calculated for two types of graphene fillers with the same surface area but in-plane aspect ratios of 1 and 10. The filler with the higher aspect ratio was found to exhibit a much higher thermal conductivity enhancement than the one with the lower aspect ratio. This is because a high in-plane aspect ratio strongly restricts the orientation of fillers when they aggregate and, consequently, highly ordered agglomerates are formed. On decomposing the effective thermal conductivity tensor into various molecular-scale contributions, it was identified that the thermal conductivity enhancement is due to the increased amount of heat transfer inside the graphene filler, particularly along the longer in-plane axis. The present result indicates a possibility of designing the heat conduction characteristics of a nanocomposite by customizing the filler shapes so as to control the aggregation structure of the fillers.


Introduction
Composites consisting of a paraffin matrix and graphene fillers have attracted interest as thermal materials, such as phase change materials (PCMs) for thermal energy storage 1 and thermal interface materials (TIMs) for reducing thermal resistance in semiconductor devices. 2 While paraffin has many advantages including light weight, high flexibility, and low cost, it has the disadvantage of a low thermal conductivity, of the order of 10 À1 W (m K) À1 . Graphene, whose thermal conductivity has been reported to be in the range of 1000-5000 W (m K) À1 , 3,4 is a suitable filler to enhance the thermal conductivity of paraffin. If the effect of the filler is proportional to the contact area between the matrix and the filler, a smaller filler would be more advantageous, as a larger contact area could be achieved at the same filler concentration. Therefore, attempts have been made to use various types of nanographene fillers. 5,6 The observed values for the effective thermal conductivity of actual paraffin/graphene composites are rather scattered in the range of the order of 10 À1 -10 1 W (m K) À1 for filler concentrations less than 20 wt%, [7][8][9][10] and these values are much lower than expected when compared with the high thermal conductivity of graphene. To understand this result, many molecular-scale factors must be considered in addition to thermodynamic conditions including temperature, pressure, and filler concentration.
The geometric shape of graphene fillers is one such factor. The thermal conductivity of a single-layer graphene (SLG) with nanoscale in-plane dimensions is typically lower by one or two orders of magnitude than that of macroscopic graphene 4,11 in particular by enhanced scattering of ballistic phonons at the filler edges. In addition, the effect of interfacial thermal resistance between paraffin and graphene, estimated to be of the order of 10 À9 -10 À7 m 2 ÁK W À1 , 7,10,12 increases with decreasing filler size. These effects of the filler size may be a reason for the insufficient thermal conductivity enhancement.
The aspect ratio of the graphene filler is also an important factor. Some authors have reported that filler percolation dramatically enhances effective thermal conductivity. 13,14 In general, filler particles with a higher aspect ratio can achieve the percolation threshold at lower filler concentration. 15 The thermal percolation, however, may not always be visible because the contrast in thermal conductivity between matrix and filler is several orders of magnitude lower than the contrast in electric conductivity in the case of electric percolation. 16,17 Some of the recent effective medium theories take into account the aspect ratio of the in-plane width to the cross-plane thickness of multi-layer graphene (MLG) fillers, 17,18 but the effect of the in-plane aspect ratio (the ratio of one in-plane side to the other in-plane side) has rarely been discussed.
Effective thermal conductivity also depends strongly on the detailed arrangement of matrix and filler molecules in a composite. It is known that pristine graphene fillers aggregate easily in a polymer matrix, especially by stacking in the cross-plane direction, resulting in a poor filler dispersion. 19,20 Although it is generally believed that better filler dispersion leads to higher effective thermal conductivity, the effect of filler aggregation in reality is more complicated as summarized in the review paper of Gu et al. 11 The thermal conductivity of an MLG suspended across a trench decreases with increasing number of graphene layers, N layer , because of the increased scattering of the out-ofplane phonon modes. 21 Molecular dynamics (MD) simulation showed that an MLG in a vacuum behaves similarly. 22 Since the interaction with a substrate also scatters the out-of-plane phonon modes of graphene, the thermal conductivity of an SLG supported on a hard substrate reduces to the order of 10 2 W (m K) À1 . 23,24 In contrast to the suspended case, the thermal conductivity of the supported graphene increases with N layer , as the graphene-substrate interactions per single sheet become weaker. Both MD simulation 18 and experiments 7 indicate that a sufficiently stiff matrix can play the same role as the substrate. Thus, graphene stacking (poor filler dispersion) either increases or decreases the effective thermal conductivity of a composite, depending on whether the graphene filler in the composite is close to the suspended or supported graphene, respectively.
As for molecular orientation, the effective thermal conductivity of a polymer/graphene composite is usually higher when graphene fillers are more aligned, 2,25,26 although the opposite behavior has been observed in some cases. 27 In addition, both MD simulations 28,29 and experiments 8 have shown that paraffin molecules in the vicinity of a graphene or carbon nanotube filler are forced to align along the fillers and thus contribute to the enhancement of thermal conductivity. Consequently, heat conduction in a paraffin/graphene composite is anisotropic.
MD simulation has been utilized to investigate the effect of molecular-scale factors on the thermal conductivity enhancement in paraffin/graphene composites, 10,12,25,[28][29][30] as well as other polymer/graphene composites. 5,18 Most of these MD studies dealt with a system where a single filler of SLG or MLG was embedded in a matrix or multiple SLGs were uniformly dispersed, 25 and therefore the aggregation of fillers except for simple stacking has not been fully discussed. In addition, the effective thermal conductivity of composites has been calculated mostly using non-equilibrium MD (NEMD) simulation. In order to examine the anisotropy in thermal conductivity using the NEMD method, one must perform multiple simulations with a temperature gradient imposed in different directions. [28][29][30][31] Ideally, it is desirable to express thermal conductivity as a second-rank tensor and examine heat conduction in the principal directions. To this aim, the equilibrium MD (EMD) simulation based on the Green-Kubo formula is more suitable than the NEMD method because all components of the thermal conductivity tensor can be obtained in a single simulation run.
In the present study, we used such an EMD evaluation of the thermal conductivity tensor to calculate the effective thermal conductivity of a paraffin/graphene composite as a function of filler concentration. The results of two nanographene fillers with the same surface area but different in-plane aspect ratios were compared to show the significant effect of the shape of the graphene filler on thermal conductivity enhancement. In our simulation, the filler aggregation and the resultant anisotropy in heat conduction were explicitly considered by using multiple graphene nanofillers and by representing effective thermal conductivity as a tensor. In addition, heat conduction in each principal direction was further decomposed into various molecular-scale contributions 32 to obtain a detailed picture of heat conduction in the composite.
To provide reference information, we also computed the thermal conductivity of an SLG and the melting point of paraffin for the models used in the present study. The computational methods for these values are briefly described in the ESI. †

Thermal conductivity tensor for a binary system
Let us consider a two-component system consisting of molecular species I = 1 and 2 in a volume V in a non-equilibrium steady state under a constant temperature gradient rT. The volume-averaged mass flux vector of species I is defined as where m i and _ r i are the mass and velocity vector of atom i, respectively, and r I ¼ P i2I m i =V and u I ¼ P i2I m i _ r i = P i2I m i are the mass density and streaming velocity of species I, respectively. We assume that there is no net mass flow at any moment of time as and that mass flux for each species on average is zero as where h. . .i ne denotes the time average in the non-equilibrium steady state. The conductive part of the heat flux, which essentially contributes to thermal conductivity, can be obtained by removing, from the total heat flux, the contribution from the convection of each species: 33,34 where f i and r i are the potential energy and the position vector of atom i, respectively, r Ã ij is the fragment of vector r ij = r i À r j contained in the volume V, and the operator # indicates the tensor product. If V denotes the entire MD box under the threedimensional periodic boundary conditions, r Ã ij can be replaced with r ij . The i-j force vector f ij = Àqf j /qr i , in general, depends on the coordinates of i, j, and many other atoms. In the right-hand side of eqn (4), the first term is called the transport term since it expresses the thermal energy transfer associated with the transport of atoms excluding the macroscopic mass flow. The second term is named the interaction term since it explains the thermal energy transfer between different atoms via inter-and intramolecular interactions. Within the framework of linear irreducible thermodynamics, the steady state averages of heat and mass fluxes can be linearly expanded using the gradients of temperature and chemical potentials as 33,35 where m I is the chemical potential of the species I and r T indicates the gradient measured at a constant temperature. In addition, J 2 was removed using the condition of eqn (2). The coefficient L AB with A, B = Q or 1 is the Onsager coefficient in the form of a second-rank tensor. Using eqn (3), (5) can be solved as hJ Q i ne = À(L QQ À L Q1 L À1 11 L 1Q )/T 2 ÁrT, which defines the thermal conductivity tensor l as l = l Q + l Mass (6) with l Q = L QQ /T 2 and l Mass = ÀL Q1 L À1 11 L 1Q /T 2 . Here, l Q describes the main part of the thermal conductivity that is associated with the heat flow directly induced by the temperature gradient, while the heat-mass coupling term l Mass explains an indirect process in which the temperature gradient first changes the density distributions of the two species and then the change in density distribution causes the heat flow. Since the tensor L Q1 L À1 11 L 1Q is positive definite, the eigen values of l Mass are negative, i.e., the heat-mass coupling effect always decreases thermal conductivity.
The Green-Kubo formula relates the running value of the Onsager coefficient, L AB (t), to a correlation function in the reference equilibrium state as 36,37 where k B is the Boltzmann constant and J A ðt 0 Þ J B ð0Þ h i eq is the time correlation function between J A and J B at a time difference t 0 . The reference equilibrium state signifies the equilibrium state which is spontaneously reached when rT applied in the non-equilibrium steady state is turned off. After a sufficiently large time t, L AB (t) becomes constant, and this value is equal to L AB appeared in eqn (5). Using eqn (6) and (7), one can derive the thermal conductivity tensor l as a 3 Â 3 matrix from EMD simulation. The three eigenvectors of l, e 1 , e 2 , and e 3 , indicate the principal axes of heat conduction, and the corresponding eigenvalues l 1 , l 2 , and l 3 represent the thermal conductivity with respect to these axes. We refer to the largest and lowest eigen values as the first and third principal thermal conductivities, respectively.
The heat flux in eqn (4) can be expressed by the sum of the partial heat flux J X due to a specific interaction or a specific group of atoms, J Q ¼ P X J X . In the present study, molecular species I was either paraffin (P) or graphene (G), and the total heat flux in eqn (4) was decomposed as follows: where the first and the second parentheses in the right-hand side correspond to the transport and interaction terms in eqn (4), respectively. The interaction term was decomposed according to whether interacting atoms i and j belong to paraffin or graphene and whether they belong to the same molecule (Intra) or not (Inter). The partial thermal conductivity tensor l X corresponding to J X can be calculated as 32 and, taking eqn (6) into account, the total thermal conductivity tensor can be decomposed as follows: The contribution to the principal thermal conductivity l a (a = 1, 2, or 3) from a partial thermal conductivity can be obtained by the transformation P À1 l X P, where P is the matrix that diagonalizes the total thermal conductivity tensor l, and not l X .

Equilibrium MD simulations of a paraffin/graphene composite
Nonacosane (C 29 H 60 ) was considered as the paraffin matrix in the present study. This molecule was also used in the experimental study on the thermal conductivity enhancement of paraffin/graphene PCM by Shi et al. 9 The nonacosane molecule was modeled using the NERD potential. 38 This potential employs a united atom model, where a CH 2 unit and a CH 3 unit are considered as the single interaction sites and the sites interact via van der Waals (vdW) interactions and intramolecular bond stretching, angle bending, and dihedral interactions. The graphene filler was modeled using the AIREBO potential, 39 which satisfactorily describes the covalent bonds between C and H atoms in the sp, sp 2 , and sp 3 hybridization states as well as vdW interactions among the non-bonded atoms. For both NERD and AIREBO potentials, vdW interactions are described using the Lennard-Jones (LJ) potential, where the distance parameter s ab and the energy parameter e ab for heterogeneous atom types a and b are given using the Lorentz-Berthelot combining rules. In the present study, the vdW interactions between nonacosane and graphene atoms were modeled using the LJ potential in the same manner.
The cutoff distance for the LJ interaction between graphene atoms was set to 3s ab according to the original parametrization of the AIREBO potential, 39 while the cutoff distance of 15 Å was used for other atom pairs.
To investigate the effect of the in-plane aspect ratio of the graphene filler, two filler shapes with a similar surface area were used. One had an in-plane aspect ratio of B1 as shown in Fig. 1(a) and is referred to as the square filler hereafter, while the other had an in-plane aspect ratio of B10 as shown in Fig. 1(d) and is called the ribbon-like filler. For both fillers, the shorter and the longer edges had armchair and zigzag structures, respectively, and all edges were terminated with hydrogen atoms. The square filler consisted of 284 C atoms and 78 H atoms, and the ribbon-like filler consisted of 276 C atoms and 46 H atoms.
All MD simulations in the present study were performed using the LAMMPS software package. 40 The algorithm of Martyna et al. 41 was used with a timestep of 0.5 fs for the numerical integration of the equation of motion to generate various statistical ensembles such as the NVE (micro canonical), NVT (canonical), and NpT (isobaric-isothermal) ensembles. We modified some source codes so that heat flux was computed according to eqn (4) with the decomposition as in eqn (8). The detailed expression for the i-j force f ij in eqn (4) that we employed here for NERD and AIREBO potentials can be found in ref. 42 and 43, respectively. It should be noted that these expressions satisfy the energy conservation law, which is not guaranteed by the original heat flux formula in LAMMPS. 44,45 We constructed several MD systems of the composites at 360 K using either the square or ribbon-like filler with different filler mass fractions w G up to relatively high values of filler concentration. The simulation settings for all cases examined are summarized in Table 1. The system with w G = 0 wt% corresponds to the pure paraffin matrix, containing 1440 nonacosane molecules. The maximum filler concentration was w G = 60.5 wt%, containing 140 sheets of square graphene.
The MD box was initially a cube with side length of 170 Å, in which graphene fillers and paraffin molecules were arranged regularly in a lattice-like pattern with fixed molecular orientations. Graphene sheets were sufficiently far from each other (an example of this molecular arrangement is shown in Fig. S3 in the ESI †). The initial velocities of atoms were randomly assigned from the Maxwell velocity distribution at 10 K. The system temperature was then increased to 600 K at a constant rate by a 0.5 ns NVE run with velocity scaling at 50 fs interval. Next, the system was agitated by a 9.5 ns NpT run at 600 K and 1 atm to erase the memory of the initial molecular arrangement. After this, the system was cooled down to 360 K over 3 ns and equilibrated at 360 K for 6 ns. In these NpT runs, a pressure of 1 atm was , and (f), the axes l 1 , l 2 , and l 3 indicate the first, second, and third principal axes of the thermal conductivity tensor, respectively. Table 1 Conditions for simulations of a paraffin/graphene composite, where w G is the mass fraction of the graphene fillers, N P is the number of nonacosane molecules, N G is the number of graphene sheets, L is the side of the cubic simulation box, and t is the length of the production run. The last column shows the simulation results of effective thermal conductivity l (isotropic average), where the figure in parentheses indicates uncertainty in the last digit applied isotropically. At this point, the volume of the MD box was fixed, and the system was further relaxed by a 3 ns NVT run followed by a 3 ns NVE run. It was considered that an equilibrium state was reached at this point, since no significant drift in the potential energy was recognized during the last NVE run. Finally, the production NVE run was then performed for t = 30.0-42.5 ns. The total simulation time t was equally partitioned into five segments to estimate the statistical error of a physical quantity A as sðAÞ= ffiffi ffi 5 p , where s(A) is the standard deviation of the five average values computed for the five segments. We note that the temperature chosen here, 360 K, is higher than the melting temperature of pure nonacosane, which was calculated to be 348 K as described in the ESI. † Thus, we focused on the liquid-like states of the composite. While the self-diffusion coefficient of the nonacosane matrix in the composite was D B 2.8-6.2 Â 10 À10 m 2 s À1 , that of a graphene filler was much lower and the filler arrangement did not significantly change within the simulation time t. Here, selfdiffusion coefficient D for each component was estimated using the Einstein relation 6D ¼ lim t!1 dxðtÞ=dt, 46 where x(t) is the mean square displacement of the center-of-mass coordinates of a paraffin molecule or a graphene sheet at time difference t, and its slope for t -N was approximately obtained from the linear fit to the x-t curve for 1 r t r 1.5 ns.

Effective thermal conductivity
The running values of the principal thermal conductivities were calculated for each time point t by diagonalizing the temporal thermal conductivity tensor constituted of L AB (t) in eqn (7). As an example, the results of the sample with 20 wt% ribbon-like fillers are plotted in Fig. 2. The three curves show different behaviors, indicating the anisotropy of heat conduction in the composite. For all cases examined, the curve was approximately constant after t = 8 ps, and the steady-state value of the principal thermal conductivity was determined as the time average over t = 8-10 ps.
The isotropic thermal conductivity was calculated as the average of the three principal thermal conductivities, l = (l 1 + l 2 + l 3 )/3. The results are shown in Fig. 3(a), as a function of graphene mass fraction w G . The right axis shows the enhancement factor Z = ( l À l mat )/l mat Â 100 relative to the matrix thermal conductivity l mat . As shown in Table 1, the thermal conductivity of pure nonacosane was calculated to be l mat = 0.130 AE 0.002 W (m K) À1 . The isotropic thermal conductivity increased with graphene mass fraction. In the case of square filler, the enhancement was almost invisible for small graphene mass fractions until the enhancement factor increased to Z = 9.8% at w G = 10 wt%.
The enhancement factor for the liquid-like states is usually lower than that for the solid-like states. 9,28,29 Shi et al. measured the effective thermal conductivity of liquid-like nonacosane/ graphene composites for w G r 10 wt%. They observed the maximum enhancement of Z B 100% at 2 wt%, but the Fig. 2 First, second, and third principal thermal conductivities, l 1 , l 2 , and l 3 , respectively, as a function of time t for the paraffin/graphene composite with 20 wt% ribbon-like fillers. l of the principal thermal conductivities, l 1 , l 2 , and l 3 , for the ribbon-like and square graphene fillers; (b) principal thermal conductivities for the square filler; (c) principal thermal conductivities for the ribbon-like filler. enhancement was not significant at some other filler concentrations. 9 Our MD results are qualitatively consistent with their observation, but in our case, another reason for the low enhancement factor comes from the small size of the graphene filler. As described in the ESI, † the thermal conductivity of a single square filler was 54 AE 5 W (m K) À1 and that of a ribbonlike filler was 66 AE 4 along the longest axis, and these values are two orders of magnitude smaller than those of macroscopic graphene. We note that in many experiments, Z reaches several hundred percent at filler concentrations less than 10 wt%, 10,27 but these values were obtained for the solid-like states using graphene fillers of micrometer sizes. The ribbon-like filler shows a much higher enhancement than the square filler. When compared at w G = 40 wt%, the enhancement factor was Z = 79% for the ribbon-like filler while Z = 14% for the square filler, demonstrating clearly the effect of the in-plane aspect ratio of the graphene filler.
The principal thermal conductivities l 1 , l 2 , and l 3 for the square and ribbon-like fillers are shown in Fig. 3(b) and (c), respectively. For both the filler types, the anisotropy in heat conduction, indicated by the difference among the three principal components, increased with increasing graphene mass fraction. We note that the pure paraffin case also showed a small anisotropy as indicated by l 1 = 0.142 AE 0.006 W (m K) À1 , l 2 = 0.128 AE 0.002 W (m K) À1 , and l 3 = 0.122 AE 0.006 W (m K) À1 , since the molecular alignment is not perfectly isotropic at a certain moment of time owing to the limited size of the system. Therefore, if the anisotropy in the composite is as small as at this level, we cannot judge whether or not it is due to the addition of graphene fillers. As w G increases, l 3 remains close to the pure matrix value, indicating that the enhancement in the isotropic thermal conductivity l is mostly due to l 1 and l 2 . This also means that the different degrees of enhancement between the square and ribbon-like fillers seen in Fig. 3(a) originate from l 1 and l 2 . The effect of filler shape observed in Fig. 3 cannot be explained by the 22% difference in thermal conductivity of a single layer between the square and ribbonlike fillers. A more important factor is that fillers with different shapes aggregate differently in a composite, as will be discussed in the next section.

Molecular structure of the composite
In our MD simulations, the graphene fillers aggregated into several clusters during the process of reaching the equilibrium state at 360 K, which is consistent with the poor dispersion of graphene fillers in actual composites. 20 When the filler concentration was lower than 5 wt%, each cluster was an MLG, which is an aggregate of several graphene sheets stacked in the crossplane direction and has an almost rectangular-parallelepiped shape. At higher filler concentrations, agglomerates of multiple MLGs were formed. These structures in the composites with 20 wt% fillers are demonstrated in Fig. 1. It is well-known 47 that the adsorption layers of paraffin molecules are formed at the filler-matrix interface. During the agglomeration process, some of these layers failed to escape from the gap between the MLGs and became part of the agglomerate as shown in Fig. 1(b) and (c).
As shown in Fig. 1(b), MLGs in an aggregate are oriented in various directions and there is no apparent correlation between the agglomerate structure and the principal axes of the thermal conductivity tensor. By contrast, in the case of the ribbon-like filler, multiple MLGs orient in the same direction. That is to say, the longest axes of all constituent graphene sheets point in the same direction as shown in Fig. 1(e) and (f), and this direction is nearly equal to the first principal axis of thermal conductivity tensor. This highly ordered alignment of the ribbon-like fillers explains the large anisotropy in thermal conductivity exhibited in Fig. 3(c).
The l 2 curve of the ribbon-like filler in Fig. 3(c) makes a sudden dip at 20 wt%. This dip is likely because the filler orientation in the directions other than the longest-axis direction was less uniform in the 20 wt% case than those in other cases with ribbon-like fillers (see Fig. S4 in the ESI † for the simulation snapshots of all cases examined in the present study). The dip in l 2 explains the relatively small inclination of the curve of isotropic average l of the ribbon-like filler at 10-20 wt% in Fig. 3(a). The size of our MD system may appear to be small in comparison to the longest axis of the ribbon-like filler. To check the effects of the system size and the periodic boundary conditions, another simulation was performed for 20 wt% ribbon-like fillers using a simulation box with a side of L = 156.2 Å, which was 1.45 times larger than the original one, and contained 114 graphene sheets and 3873 nonacosane molecules. In this system, the graphene fillers aggregated into more than one cluster having different orientations (see Fig. S4 in the ESI †). However, the effective thermal conductivity and its components were calculated to be l = 0.172 AE 0.004, l 1 = 0.23 AE 0.01, l 2 = 0.159 AE 0.005, and l 3 = 0.128 AE 0.007 in the units of W (m K) À1 . These values are close to those obtained with the original system size (Run 14 in Table 1), l = 0.174 AE 0.002, l 1 = 0.25 AE 0.01, l 2 = 0.146 AE 0.007, and l 3 = 0.13 AE 0.01 W (m K) À1 . Therefore, although the system size can affect the details of filler arrangement, the highly ordered filler orientation and the resultant high effective thermal conductivity of the ribbon-like filler are not artifacts of the limited size of the simulation box. This result is also consistent with the fact that the size effect in EMD is generally much smaller than that in NEMD simulation. 48 The structural difference between the agglomerates of the square and ribbon-like fillers can be understood by the fact that MLGs tend to stick together in such a way that the area of the joint surface is maximized as much as possible. This tendency restricts the possible patterns of aggregation for the ribbon-like fillers, for example, Fig. 4(a) and (b) illustrate the agglomerates of MLG1, MLG2, and MLG3, for the square and ribbon-like fillers, respectively. In the case of square filler, MLG1 in Fig. 4(a) has a sufficiently large surface area on any face and can form an agglomerate in which the in-plane directions of the constituent MLGs are not the same as shown in the figure. A similar aggregation pattern for the ribbon-like filler is shown in Fig. 4(b). In this case, the sticking pattern like that of MLG3 is unlikely to occur because the joint surface area with other MLGs is too small. Consequently, MLGs tend to stick together so that their longest axes point in the same direction, as MLG1 and MLG2 do in Fig. 4(b).
Although rectangular fillers only are dealt with here, it would be possible to consider the effect of in-plane aspect ratio for other filler shapes in the same way. For example, in the case of a triangular filler, the five faces of an MLG are all flat, then the effect will be quite similar to that of a rectangular filler. In contrast, for an elliptical filler, since the side faces of an MLG are curved, the sticking of two MLGs using these faces may not be stable regardless of the in-plane aspect ratio. In such a case, the effect of in-plane aspect ratio is considered to be small.
In order to further examine the structural order of filler and matrix molecules, we computed the orientation of the C-C bonds in terms of the orientational order parameter S used for the isotropic-nematic transition of liquid crystals. 49 Here, S was obtained as the largest eigenvalue of the second-rank ordering tensor S ab ¼ ð3s a s b À d ab Þ=2 , where s a is the a component (a = x, y, or z) of the unit vector along a C-C bond, d ab is the Kronecker delta, and h. . .i denotes the average over time and different bonds. The value of S varies from 0 to 1 according to the degree to which the C-C bonds are aligned in the same direction. The results for the square and ribbon-like fillers are compared in Fig. 5(a), where the order parameters computed for the C-C bonds in graphene, in paraffin, and in both molecules are separately shown.
The results confirm that the ribbon-like filler gives a higher orientational order of molecules than the square filler. In the case of the square filler, S of graphene decreases with graphene mass fraction w G as expected from the random orientation of the MLGs. In contrast, S of paraffin always increases with filler concentration, indicating the effect of graphene filler in inducing the alignment of paraffin molecules. 28,29 At w G = 40 wt%, a relatively high value of S is observed because the in-plane direction of the graphene fillers is aligned by chance as shown in Fig. 5(b). This irregularity is, however, not reflected in the curve of the effective thermal conductivity for the square filler in Fig. 3(a) and (b).
For the ribbon-like filler, S of graphene remained high up to a filler concentration of 20 wt%, which is consistent with a highly ordered filler alignment. In the case of 40 wt% ribbon-like fillers, two domains of different filler orientations are formed as shown in Fig. 5(c), reducing the S value. This reduction in S explains the relatively small slope of l 1 at 20-40 wt% in Fig. 3(c), since the S of the ribbon-like filler approximately represents how much the long axis of the filler is aligned in the first principal direction of heat conduction. The value of isotropic average l is therefore considered lower than when the two domains were oriented in the same direction. As a result of the reduction in S, on the other hand, the square and ribbon-like fillers have similar values of S at 40 wt%. Nevertheless, effective thermal conductivity is much higher for the ribbon-like filler because the size of a single orientational domain is still larger for the ribbon-like filler as can be understood by comparing Fig. 5(b) and (c). In addition to the graphene fillers, the ordering of paraffin molecules is also more significant in the case of the ribbon-like filler than in the case of the square filler.
As we saw above, effective thermal conductivity depends on the detailed filler arrangement in a composite. In the case of ribbon-like filler, effective thermal conductivity is linearly approximated as l = 0.129 + 0.00328w G for w G = 0.591-10.2 wt% (Run 11 to 13). In this range of w G , most of the fillers in a composite had the same orientation with respect to all molecular axes (see Fig. S4 in the ESI †). Therefore, this linear equation is considered to predict the upper bound of l for each value of w G that is reached when all fillers have the same orientation. The linear equation predicts that the effective thermal conductivity values for the perfect filler alignment will be l = 0.196 and 0.260 W (m K) À1 at w G = 20.5 and 40 wt%, respectively. These values are about 10% higher than our results, l = 0.174 and  0.233 W (m K) À1 , respectively. Thus, about 10% difference in effective thermal conductivity is typically expected with different filler arrangements in a composite. This level of difference does not change the superiority of the ribbon-like filler shown in Fig. 3(a). In the present study, the high computational cost restricted us to consider only one molecular configuration for each filler concentration. Therefore, the curves in Fig. 3 are not the ideal ones obtained by averaging over different initial molecular configurations. However, based on the above discussion, the main conclusion of the present study would not be affected by a specific filler arrangement, as long as the relaxation simulation is sufficiently long for graphene sheets to move in a paraffin matrix and form agglomerates of MLGs.

Molecular scale heat transfer
To examine the molecular mechanism of thermal conductivity enhancement in Fig. 3, the decomposition of the isotropic thermal conductivity l into the molecular contributions according to eqn (10) is shown in Fig. 6. For all cases, the heat-mass coupling term is a negative value as expected. Its magnitude was at most 2% of the total thermal conductivity of the composite, meaning that the heat-mass coupling is insignificant in the present study. As is typically the case with liquids in general, 50-52 the contribution of the transport term (l TransP + l TransG ) in the present case was also small with a maximum of 16% of the total thermal conductivity. The effective thermal conductivity is therefore determined by heat transfer via intra-and intermolecular interactions.
At low filler concentrations, effective thermal conductivity can be mainly explained by heat transfer among paraffin molecules through their intra-and intermolecular interactions (l IntraP + l interP ), and l IntraP is the largest contributor as is typically the case with liquid alkanes longer than decane. 51,52 As the graphene mass concentration increases, these paraffin contributions decrease their percentage and, instead, the graphene counterparts (l IntraG + l interG ) increase. In particular, increase in l IntraG is exceptional, demonstrating that the enhancement of effective thermal conductivity is mainly achieved by the increased amount of heat transfer inside graphene. Comparing Fig. 6(a) and (b) reveals that the curves of the partial thermal conductivities other than that of l IntraG are quite similar for the two filler shapes. These results indicate that it is the amount of heat transfer inside the graphene filler that is most affected by both filler concentration and filler shape. As discussed above, the contributions relevant to graphene filler and paraffin matrix increase and decrease, respectively, with increasing filler concentration. This trend is mainly due to the change in the number density of graphene sheets, n G , and that of paraffin molecules, n P . To understand this, the heat transfer efficiency per graphene sheet and that per paraffin molecule are defined as L G = tr(l TransG + l IntraG + l InterGG )/(3n G ) and and L P = tr(l TransP + l IntraP + l InterPP )/(3n P ), respectively, and plotted in Fig. 7, where tr(A) means the trace of a tensor A. As shown in Fig. 7(a), the dependence of L G on filler concentration is not significant. The relatively large value at 10 wt% for the ribbon-like filler is presumably due to a specific filler arrangement in the composite rather than to filler concentration, because in this case all fillers had the same orientation.
The change in L G is considered to correspond roughly, if not exactly, to the change in the intrinsic thermal conductivity of graphene in the composite. As described in the Introduction, thermal conductivity per single graphene sheet in an MLG depends on the number of stacking layers, and the dependence is inverted depending on whether the MLG is suspended or supported. 11 The weak dependence of L G on filler concentration in Fig. 7(a) indicates that the nature of the graphene fillers in our paraffin matrix is between that of the suspended and supported graphenes, considering that the average number of layers in an MLG is proportional to filler concentration. This behavior is reasonable as our composite is liquid-like but exhibits very low fluidity. As for paraffin, the dependence of L P on filler concentration, as shown in Fig. 7(b), is also insignificant but has a slightly increasing trend. This result implies that the thermal conductivity of the paraffin matrix is Fig. 6 Decomposition of the isotropic thermal conductivity of the paraffin/ graphene composite with (a) square filler and (b) ribbon-like filler. The total thermal conductivity was decomposed into contributions from the thermal energy transfer associated with intermolecular interaction (Inter), intramolecular interaction (Intra), molecular transport (Trans), and heat-mass coupling (Mass), and P and G indicate paraffin and graphene, respectively. The partial contribution l X is shown as the percentage to the total thermal conductivity l (see Fig. S5 in the ESI † for the plot of the absolute value of l X ).
enhanced in the composite because of the filler-induced ordering of the paraffin molecules shown in Fig. 5(a). A similar phenomenon of ordering of paraffin molecules induced by graphene fillers has also been reported in other computational 28,29 and experimental 8 studies.
Next, the same decomposition as that for the isotropic thermal conductivity was examined for each principal thermal conductivity to investigate the origin of anisotropy. As an example, Fig. 8 compares the results for the cases with 20 wt% fillers, whose simulation snapshots are displayed in Fig. 1. For the square filler, Fig. 8(a) indicates that the moderate anisotropy in the effective thermal conductivity (Fig. 3(b)) originates from that in the intramolecular heat transfers in paraffin and graphene, i.e., l IntraP and l IntraG . For the ribbon-like filler, as shown in Fig. 8(b), the intragraphene contribution in the first principal direction is prominent, and this direction is approximately parallel to the longer in-plane axis of the ribbon-like fillers as discussed before in relation to Fig. 1

(e) and (f).
This result clearly shows that the enhancement in the intragraphene heat conduction along the longest axis is the reason for the high enhancement in l 1 . As shown in Fig. 3(c), the second principal thermal conductivity l 2 also increases to some extent with filler concentration, but this enhancement is due not only to the intra-graphene heat transfer, but also to other contributions, although we do not explicitly show it here. In contrast to the intramolecular contributions, the partial thermal conductivities associated with intermolecular heat transfer were mostly isotropic for all graphene-graphene, graphene-paraffin, and paraffin-paraffin pairs, regardless of the filler shape.
Our results clearly indicate that an increased amount of heat transfer inside the graphene fillers is the dominant factor for thermal conductivity enhancement upon adding graphene fillers. In contrast, the NEMD study of octadecane/graphene composite by Babaei et al. 28,29 reached the conclusion that the thermal conductivity enhancement can be explained by increase in the thermal conductivity of the paraffin matrix. They employed a 4 Â 4 nm 2 graphene filler and observed 28% thermal conductivity enhancement with a filler concentration as low as 1.55 wt% above the melting temperature of paraffin. In their case, the orientational order parameter of paraffin, which is based on the end-to-end vector of a paraffin chain, dramatically changed by one order of magnitude upon adding one sheet of graphene. Since, as shown in Fig. 6, effective thermal conductivity at several wt% values of filler concentration is mostly explained by heat conduction in the paraffin matrix, it is reasonable that this sudden ordering of paraffin dramatically enhances heat conduction in paraffin, thereby leading to a large thermal conductivity enhancement. In our case, however, such a phase transition was not observed at least for filler concentrations of a few wt%. Thus, the detailed mechanism of thermal conductivity enhancement in their simulations and ours at low filler concentrations is different. The identification of the condition under which the strong ordering of paraffin is induced by graphene filler is a subject for future study.

Conclusions
In the present study, the effect of the in-plane aspect ratio of the graphene filler on the effective thermal conductivity of a paraffin/ graphene composite was investigated using MD simulation, where filler aggregation, anisotropy in heat conduction, and the detailed picture of molecular-scale heat transfer were explicitly considered. At relatively high filler concentrations, thermal conductivity  enhancement was explained to be mostly due to the increased amount of heat transfer inside the graphene filler. It was shown that a higher in-plane aspect ratio restricts the variation in filler orientation more when they aggregate and therefore the agglomerates formed are more ordered. Since this filler alignment further improves the intra-graphene heat transfer, fillers with a higher in-plane aspect ratio can more efficiently enhance the effective thermal conductivity. Here, the superiority of high aspect ratio fillers was demonstrated from the viewpoint of filler aggregation rather than thermal percolation. We used graphene fillers with rather small in-plane width and length, thereby leading to modest values of thermal conductivity enhancement. It is expected that, for larger fillers, thermal conductivity enhancement is higher, and accordingly the effect of the in-plane aspect ratio is more significant. The present result indicates the possibility of controlling the aggregation structure of fillers by customizing filler shapes, which is an intriguing strategy for the rational design of the heat transfer characteristics of composites. 6

Conflicts of interest
The authors declare that there are no conflicts of interest.