Atomistic simulations of phonon behaviors in isotopically doped graphene with Sierpinski carpet fractal structure

Two-dimensional (2D) graphene monolayer has been attached importance because of the fantastic physical properties. In this work, we conduct the atomistic simulations to evaluate the phonon behaviors in isotopically doped graphene with Sierpinski Carpet (SC) fractal structure. The thermal conductivities (k) with different fractal numbers are calculated by molecular dynamics simulation. The relationship between the k and the fractal number from 0 to 8 shows a first decreasing and then stable trend. The maximum reduction ratio of the k in SC fractal structures is 52.37%. Afterwards, we utilize the molecular dynamics simulation, phonon wave packet simulation and lattice dynamics simulation to investigate the phonon density of states (PDOS), energy transmission coefficient (ETC), phonon group velocity and participation ratio (PR) in SC fractal structures. In SC fractal structures, the PDOS increases in the low frequency region and the G-band will soften with the enhanced fractal number. We also observe that the isotopic doping atoms can lead to continuous reflected waves in SC fractal structure regions. Moreover, phonon modes in SC fractal structures possess the lower ETCs, phonon group velocities and PRs in comparison with the pristine graphene monolayer. Therefore, we attribute the lower k in SC fractal structures to the stronger phonon-impurity scattering and the increasing localized phonon modes.


Introduction
Graphene, consisting of covalently bonded carbon (C) atoms, is a novel two-dimensional (2D) material with a honeycomb lattice, which has triggered gigantic investigations because of the high charge carrier mobility [1], strong fracture strength [2], and ultrahigh thermal conductivity [3]. Due to the absence of direct bandgap in pristine graphene, the modified graphene structures have been investigated to gain the desired electrical properties, such as unfolded band gap [4], excellent Seebeck coefficient [5][6][7], and adjustable carrier mobility [8]. Moreover, the thermal properties of graphene play a crucial role in the thermal management of nanoelectronic devices to obtain an outstanding electrical performance [9]. In addition, the thermal-to-electric energy conversion capability relies partly on the thermal properties of materials [10]. Therefore, tuning the thermal properties has emerged as an indispensable channel to obtain fascinating physical properties we need, which can be applied to thermal management and thermoelectric applications.
For the past few years, the major methods to tune the thermal conductivity (k) of graphene are introducing the defects or isotopic doping atoms [11], conducting superlattice structures [12], adding different strains [13] and building some specific structures [14][15][16]. Hu et al investigate the effects of the defects and isotopic doping atoms on the thermal conductivity of graphene using the molecular dynamics (MD) simulation [11]. They find that both defects and isotopic doping atoms can decrease the k monotonically with the increasing defect ratio and isotopic doping ratio. By changing the period length of the graphene and boron nitride superlattice from 0.992 to 14.88 nm, Chen et al achieve the minimum thermal conductivity with different temperatures [12]. It is reported that the k of graphene nanoribbon (GNR) is insensitive to compressive strain but drops sharply under Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. tensile strain [13]. The maximum reduction ratios of the thermal conductivity in zigzag GNR and armchair GNR are ∼60% and ∼40% by applying different strains, respectively. Moreover, conducting some specific structures is an effective approach to manipulate the thermal conductivity of the graphene monolayer [14][15][16][17]. Cui et al investigate the k of a novel graphene wrinkle structure using the MD method, which is dropped by 36% and 52% in the parallel and perpendicular directions at 300 K [14]. They attribute this phenomenon to the increasing phonon scatterings and the localized phonon modes in graphene wrinkle structure. Another especial structure is the phononic crystal, which can significantly reduce the thermal conductivity in graphene by changing the factors (period length from 10.4 to 62.5 nm and porosity from 0 to 66%) [15]. Furthermore, the k of bilayer graphene phononic crystal with an 8 nm neck width is experimentally measured by Xiong et al [16], which is Recently, as the nanofabrication techniques make great progress, the self-similar structures formed by the nanocrystal self-assembly have been accomplished [18]. The Sierpinski Carpet (SC) fractal structure is one of the self-similar sets, which is fabricated by the technique of dividing a square into some smaller squares, removing one square, and continuing recursively [19]. It is reported that the bandgap in phononic crystals dominating the elastic wave propagation can be changed by conducting the SC fractal structures, which is widely used in the fractal antennas design [20][21][22][23]. Furthermore, the effect of SC fractal structure on the k of the Si/Ge nanocomposites is evaluated by using MD method [24]. It can be observed that the enhanced phonon scattering at the Si/Ge interfaces makes the larger reduction of thermal conductivity in SC fractal structures in comparison with the regular carpet structures. Moreover, Kang et al evaluate the thermal conductivity of graphene with SC fractal defects by employing the MD method [25]. They observe that when the fractal number increases from 0 to 3, the thermal conductivity drops from 164.40 to 19.60 W m −1 K −1 . However, they only investigate the effect of small fractal number (from 0 to 3) on thermal conductivity in SC graphene monolayer. In addition, the phonon behaviors in graphene monolayer, such as the phonon density of states (PDOS), phonon group velocity, and participation ratio (PR), are not analyzed in that work. Accordingly, we presume that the SC fractal structure is a promising structure to tune the thermal conductivity of 2D materials, as we illustrate above. Here, we investigate the combined effects of the SC fractal structures and the isotopic doping atoms in graphene monolayer on the phonon behaviors to enrich the adjustable channel of the thermal properties in 2D materials.
In this work, we investigate the phonon behaviors in isotopically doped graphene with SC fractal structure using the typical atomistic simulations. First, the thermal conductivity of graphene monolayer with the isotopic doping and SC fractal structure is calculated by using the MD simulation. Furthermore, we analyze the phonon behaviors, such as PDOS, the propagation of phonon wave packet (PWP), energy transmission coefficient (ETC), phonon group velocity and PR using the MD simulations, PWP simulations and lattice dynamics (LD) simulations, respectively, which are conducive to uncover the underlying phonon transport mechanisms in graphene monolayer.

Model
The Sierpinski Carpet fractal structure is regarded as the self-similar sets, and the establishment process is as follows: We first divide a square to 9 identical sub-squares by trisecting each side of this square. Then the central sub-square can be removed to create the SC fractal structure with the fractal number ( f n ) of 1. As shown in figure 1(a), this structure is named as SC1. Afterwards, the same trisection procedure is continuously applied to divide and remove the remaining 8 sub-squares. The SC2 and SC3 structures are displayed in figures 1(b), (c). The self-similar fractal dimension (D) of the SC fractal structure is calculated by D=lgb/lga=lg8/lg3=1.89, where a is the number of square side lengths evenly divided and b is the ratio of the remaining area to the removing area of the square at the f n th SC fractal structure. In our work, the SC fractal structures are modeled with the fractal number ranging from 0 to 8, and the fractal number of the pristine graphene monolayer is 0. It should be noticed that all removed central sub-squares are substituted by isotopic doping atoms. We model the unit cell with a size of 9.37×9.35 nm 2 (x×y nm 2 ) and the total length obtains 9 unit cells in the x direction, which includes 30 096 atoms with a size of 84.33×9.35 nm 2 . Schematic of the SC2 graphene monolayer with 3 unit cells is displayed in figure 1(d). The armchair and zigzag directions in graphene monolayer are along the x and y directions, respectively. In addition, 13 C can be selected as the isotopic doping atoms [26], which substitutes the removing atoms in SC fractal structures, and the other atoms can be modeled by natural atoms ( 12 C).

Molecular dynamics simulation
In this work, all MD simulations are performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [27]. The optimized Tersoff potential [28] is adopted to model the C-C interactions in the graphene monolayer, which accurately reproduces the phonon dispersion curves and predicts the thermal properties of graphene monolayer [29][30][31][32]. Moreover, it is reported that the optimized Tersoff potential is a better choice to investigate the thermal properties of graphene monolayer in comparison with the original Tersoff, the REBO, and the AIREBO potentials, in which the optimized Tersoff has some advantages as follows: (1) the thermal conductivity value closer to the experimental results, (2) the accurate phonon dispersion curve, (3) the accurate frequency value of G-band in PDOS, (4) the correct phonon group velocity and relaxation time, and (5) the larger contribution of ZA phonon mode to the total thermal conductivity [31,33]. The periodic boundary conditions are applied to the y and z directions, while the free boundary is used in the x direction. To prevent the interaction with the system image in the z direction, we add a large vacuum distance with 6 nm. The outmost atoms on both sides are fixed with 0.3 nm. We set the time step as 0.5 fs. The thermal conductivity of graphene monolayer is calculated by using a non-equilibrium molecular dynamics (NEMD) method, and the simulation process is as follows: we first place the initial structure in the canonical ensemble (NVT) for 500 ps at 300 K and then switch it to the microcanonical ensemble (NVE) for 500 ps to relax the structure. The NEMD simulation is conducted under the NVE ensemble for 6 ns. The atoms next to the fixed atoms on both sides serve as the heat source and heat sink with 1 nm as shown in figure 2(a), which are set as 350 K and 250 K by using a Langevin thermostat, respectively, thus the heat flux direction is along the x direction. To calculate the temperature distribution, we divide some slices with 0.2 nm along the x direction and record the kinetic energies of atoms in each slice for the last 1 ns. Afterwards, the temperature distribution based on the energy equipartition theorem is calculated by the following formula where 〈E〉 represents the mean kinetic energy, N is the number of atoms in each slice, m is the atomic mass, v i represents the velocity of atom i, and k B is the Boltzmann constant. Moreover, the thermal conductivity based on the Fourier's law of heat transfer is given by where Q is the heat flux, A represents the cross-sectional area, k is the thermal conductivity, T is the temperature, and x represents the heat flux direction. The exchanged energies in the heat source and heat sink are recorded to calculate the heat flux, which are depicted in figure 2(b). We select the graphite interlayer distance (0.34 nm) as the thickness of the graphene monolayer to calculate the cross-sectional area [34]. For each simulation, we conduct five simulations with different velocity seeds, and the average final values are calculated with the standard deviations shown as error bars.

Lattice dynamics simulation
In LD simulations, we model the SC0, SC2 and SC6 graphene monolayers with a size of 5.54×5.66 nm 2 (x×y nm 2 ) using the same trisection approach mentioned in the previous section, which consists of 1196 carbon atoms. After performing the LD simulations using General Utility Lattice Program (GULP) [35], the phonon dispersion curves are calculated along the Γ to Χ direction. The phonon group velocities are obtained based on the phonon dispersion curves. The phonon participation ratio (P λ ) can be investigated by [36][37][38] å å e e = l a a l a l - where N is the number of the total atoms for system, ε iα,λ is the eigenvector component of mode λ at the α direction of the ith atom and the superscript ' * ' represents the complex conjugate. The eigenvector can be calculated by LD simulations using GULP.

Phonon wave packet simulation
Periodic boundary conditions are applied in all three directions, and the vacuum distance is still 6 nm. We choose the SC0 structure as the initial material to excite the phonon wave packet. For SC0+SC0, SC0+SC2 and SC0+SC6 systems, the length ratio between SC0 and SC0/SC2/SC6 region is 1:1. To avoid the interaction with the phonon wave packet image along the x direction, the length of the system is set as 1107.78 nm. Before adding the phonon wave packet, all systems are relaxed using the Polak-Ribiere version of the conjugate gradient (CG) algorithm at 0 K. The phonon wave packet is formed by the following formula [39,40] where u il,α is the displacement component at the α direction (x, y and z) of ith atom in lth unit cell, M represents the amplitude of the phonon wave packet and k 0 means the wavevector. The phonon wave packet is centered at x 0 in real space. The parameter σ represents the width of the phonon wave packet. The eigenvector and wavevector in equation (4) can be calculated by LD simulations using GULP [41]. To initialize phonon wave packet velocities, we calculate the initial velocity by

Thermal conductivity
We first evaluate the thermal conductivity of the isotopically doped graphene monolayer with the increasing fractal number from 0 to 8 to reveal the effects of isotopic doping and fractal structure. The MD results are shown in figure 3. The SC0 monolayer has the thermal conductivity of 1228.59 W m −1 K −1 , which is similar to our previous work [42]. It can be observed that when the fractal number increases, the k of SC fractal structures first decreases and then tends to be stable. The minimum thermal conductivity of SC fractal structures is 585.20 W m −1 K −1 with a fractal number of 6, and the corresponding reduction ratio is 52.37%. A similar decline trend in low fractal number region is observed in other work [25]. We attribute this phenomenon to the phonon mode localization. The increasing fractal number indicates the elevated isotopic doping percentage. The 13 C atoms in SC fractal patterns can generate localize phonon modes and enhance the phonon-impurity scattering, which leads to the rapid decline of thermal conductivity. These localized phonon modes will constantly produce with the increasing fractal number, which further strengthen the phonon-impurity scattering, thus the thermal conductivity can be continually reduced until the fractal number reaches 6. For the fractal number >6, the percentage of the 13 C atoms is over 50%, these isotopic doping atoms become dominant and the natural atoms become the dopants. The effect of phonon-impurity scattering produced by isotopic doping atoms will be weakened. Therefore, the thermal conductivity will tend to be stable and exhibits a slight ascent. As the fractal number increases continuously (>8), the SC graphene monolayer is trending to a pure graphene monolayer with 13 C atoms. The phonon-impurity scattering will be weakened continuously, which leads to a higher thermal conductivity.

PDOS
The analysis of the phonon behaviors plays an essential role in explaining our results obtained by MD simulation and expounds on the physical mechanisms of heat conduction. One of the effective methods of analyzing phonon behavior differences is to calculate the PDOS, which can reflect the variation of the intrinsic vibration mode after the system is affected by different modulators. In this work, the SC0, SC2 and SC6 graphene monolayers are selected to investigate the PDOS. We first relax the structures in the NVT ensemble for 0.5 ns at 300 K and then switch them in the NVE ensemble for 0.5 ns. The velocities of atoms are recorded for the next 0.1 ns, and the PDOS can be calculated by taking the Fourier transform of the velocity autocorrelation function [43] where v is atomic velocity, and the angle brackets denote ensemble averaging. The PDOS of SC0, SC2 and SC6 graphene monolayers are depicted in figure 4. There are two dominating reasons that can affect the thermal conductivity in the system. On one hand, the flattened PDOS curves (for low frequency and medium frequency regions) and the softened G-band (for high frequency region) in SC2 and SC6 structures can broaden the phonon modes, which can decrease the mean free path of these phonon modes [44]. In addition, the PDOS in SC6 is more flattened and softened than that in SC2. Therefore, the contribution of these variational phonon modes to thermal conductivity is diminished. On the other hand, the phonon relaxation time due to the isotopic doping atoms (τ d ) can be defined as [11] where f is the isotopic doping ratio and G(ω) represents the PDOS value. The phonon relaxation time is inversely proportional to the isotopic doping ratio and PDOS. It is reported that the low frequency phonon modes have more contribution to the thermal conductivity in graphene monolayer [45]. According to figure 4 and other work [44], we define the region smaller than 15 THz as the low frequency region. We calculate the average PDOS in low frequency region are 6.76, 6.88 and 7.64 for SC0, SC2, and SC6 graphene monolayers, respectively. The increase of PDOS and isotopic doping ratio in low frequency region can lead to a further reduction in the phonon relaxation time based on equation (7). Apart from the above factors, a G-band redshift in SC fractal structures is observed from figure 4, which reveals that the flattened phonon dispersion curves decrease the phonon group velocity. Hence, the phonon transport will be impeded with the elevated fractal number.

PWP simulation and ETC
Furthermore, to understand the phonon transmission mechanism, we perform the PWP simulations. As we mention in the previous section, in SC0+SC0, SC0+SC2 and SC0+SC6 systems (figure 5), the SC0 structure is selected as the initial material to excite the PWP. We choose the wavevector k 0 as 0.35×(2π/a) to generate the TA phonon wave packet, where a is the lattice constant, the corresponding frequency is 11.64 THz. The propagation of phonon wave packets in SC0+SC0, SC0+SC2, and SC0+SC6 systems are displayed in figure 6. The moving directions of phonon wave packets in three systems are along the x direction and the real space position of launching phonon wave packet is 250 nm. In SC0+SC0 system ( figure 6(a)), the phonon wave packet is not deforming and continuously moving along the x direction although the wave packets cross the interface at 24 ps. We can attribute this phenomenon to the ballistic phonon transport due to the negligible anharmonic phonon scattering. However, for SC0+SC2 and SC0+SC6 systems, when the time is 24 ps, the phonon wave packets are separated into the transmitted wave and the reflected wave as shown in figures 6(b), (c).
After the wave packets cross the interface, the transmitted waves continue to move along the positive x direction, nevertheless, the reflected waves propagate along the negative x directions. We can observe that there are continuous reflected waves in both SC0+SC2 and SC0+SC6 systems. Moreover, the phonon wave packet can separate more reflected waves in SC0+SC6 system. This phenomenon can be ascribed to the consecutive phonon-impurity scatterings generated by the isotopic doping atoms in SC fractal patterns.
To intuitively display the amount of transmitted energy when the phonon wave packets impact different interfaces, we calculate the ETC with different phonon frequencies. We choose the wavevector k 0 from 0.04×(2π/a) to 0.83×(2π/a) to generate the frequency-resolved PWP, the corresponding frequency is from 1.42 to 22.85 THz. The moving directions and real space positions of phonon wave packets in SC0+SC0, SC0 +SC2, and SC0+SC6 systems remain unchanged. It should be noticed that due to the periodic boundary condition, the wave packet crossing the right boundary represents entering the left boundary, thus we only calculate the energy when the phonon wave packet reaches the right boundary to avoid the effect of the excess energy in the left region. The frequency-resolved ETCs are depicted in figure 7. In the SC0+SC0 system, the ETC is almost equal to 1 in overall frequencies because of the ballistic phonon transport, the corresponding wave packet is not split and maintain the same shape until it gets to the right boundary as shown in figure 6(a). However, in SC0+SC2 and SC0+SC6 systems, the phonon modes with small frequency process the longer wavelength and can ignore some isotopic doping atoms in SC2/SC6 region, in which the ETCs are close to 1. We can notice that at the edge of the Brillouin zone (k 0 →2π/a), the ETCs drop sharply, where the phonon modes have the shorter wavelength and lower phonon group velocity with high frequency. Therefore, these phonon modes have a small contribution to thermal transport in SC fractal structures. In order to make a more direct comparison, we calculate the average ETC in different systems. The average ETCs of SC0+SC0, SC0+SC2 and SC0+SC6 systems are 1.00, 0.87 and 0.79, respectively. It can be qualitatively deduced that the phonon modes carry fewer thermal energies in SC fractal structures than that in pristine graphene monolayer, which can occupy lower thermal conductivity.

Phonon group velocity and PR
Afterwards, the phonon group velocities and PRs of SC0, SC2, and SC6 graphene monolayers are evaluated by using LD simulations to gain insight into the phonon behaviors. We can observe from figure 8 that the SC fractal  structures will shrink the phonon group velocity. At Γ point, the highest phonon group velocity of SC0 graphene monolayer is 21.92 km s −1 , which is similar to the first principles result (∼22 km s −1 ) [46]. In whole frequency region, the SC0 and SC6 graphene monolayers have the highest and lowest phonon group velocities, respectively. The phonon group velocities of SC2 graphene monolayer fall in between them. To legibly compare the phonon group velocity, the average values in SC0, SC2 and SC6 graphene monolayers are 4.32 km s −1 , 3.06 km s −1 , and 2.61 km s −1 , separately. Therefore, the declined phonon group velocities result in the decrease of the thermal conductivities of SC fractal structures. This phenomenon coincides with the MD and PWP results displayed in figures 3 and 6.
We further analyze the PR to evaluate the phonon localization in SC0, SC2, and SC6 graphene monolayers, which is an effective approach to interpret the restraint of the thermal conductivity. The PRs calculated by equation (3) are illustrated in figure 9. All atoms participate in the eigenvibration in a bulk material without defects under the harmonic approximation, which leads to the non-localized phonon modes. These phonon modes can naturally diffuse in the whole space and the corresponding PRs is 1. Nonetheless, the phonon modes cannot spread to the whole space and are frequently localized in imperfect materials due to the phonon-impurity scatterings, phonon-defect scatterings and phonon-boundary scatterings, which can reduce the PRs [47]. In other words, the PR is equal to 1/N when only one atom takes part in the motion. The phonon modes can predictably transport the thermal energy with larger PRs. The phenomenon we can discern from figure 9 is that the isotopic doping atoms in SC fractal patterns promote the reduction of the PRs. Moreover, the SC6 graphene monolayer can further decline the PR of phonon modes. In whole frequency region, we calculate the average PRs of SC0, SC2, and SC6 graphene monolayers, which are equals to 0.787, 0.561 and 0.559, respectively. These results can prove that the SC fractal structures possess more localized phonon modes. When we conduct an SC fractal pattern in the graphene monolayer, the isotopic doping atoms can give rise to the different atom edges, which can generate the different phonon edge-modes. The localized phonon edge-modes can strengthen the localization of the phonon modes in SC fractal structures [48]. In addition, the 13 C atoms can act as the phonon  scattering cores and can decrease the phonon relaxation time and mean free path as we mention in figure 4, which can increase the number of localized phonon modes in SC fractal structures. Therefore, we can conclude that the thermal conductivities in SC fractal structures are suppressed by the localized phonon modes.
In short, the reduction of the thermal conductivities in SC fractal structures is attributed to three aspects: (1) the introduction of the isotopic doping atoms can increase the average PDOS in low frequency region and flatten the G-band in high frequency region, which will drop the phonon relaxation time and mean free path; (2) the phonon modes with lower ETCs in SC fractal structures transmit smaller thermal energy in comparison with the pristine graphene monolayer; (3) the lower phonon group velocities and PRs due to the isotopic doping atoms in SC fractal patterns will weaken the contribution of phonon modes to the thermal conduction and obstruct the thermal transmission.

Conclusions
In this work, we utilize the atomistic simulation techniques to evaluate the effects of the isotopic doping and fractal structure on phonon behaviors in graphene monolayers. First, the thermal conductivities of different SC fractal structures are calculated by using MD simulation. It can be observed that the k of SC fractal structures initial declines and then tends to be stable as the fractal number increases. The maximum reduction ratio of thermal conductivity is 52.37% when the fractal number is 6. Second, we investigate the PDOS, ETC, phonon group velocity and PR to get insight into the phonon behaviors. The increasing average PDOS in low frequency region and the flattened G-band in high frequency region can decrease the phonon relaxation time and mean free path with the elevated fractal number. Moreover, by performing PWP simulations, we observe that there are continuous reflected waves in both SC0+SC2 and SC0+SC6 systems in comparison with the SC0+SC0 system. The lower ETCs in SC0+SC6 system can prove the less contribution of phonon modes to thermal transport in the SC6 graphene monolayer. Finally, the average phonon group velocities and PRs of SC0, SC2 and SC6 graphene monolayers are 4.32 km s −1 and 0.787, 3.06 km s −1 and 0.561, 2.61 km s −1 and 0.559, respectively. The lower phonon group velocity and PR in SC fractal structures indicate that the stronger phonon-impurity scattering and localized phonon modes will impede the phonon thermal transport. This work is expected to provide a new approach to tune the thermal properties of 2D materials.