Length and torsion dependence of thermal conductivity in twisted graphene nanoribbons

The twist-to-writhe conversion (TWC) phenomenon, mentioned in the Introduction, is well known in twisted filamentary structures [49, 50]. It consists of releasing a rod’s torsional stress by spontaneous bending and folding, after the twist density reaches a critical value. This TWC has been shown to satisfy the Călugăreanu-White-Fuller linking number ($Lk$) theorem [50, 51, 52]:

where $Tw$ and $Wr$ are the total (real) amount of twist and a quantity called writhe of a curve which measures its non-planarity, respectively. The linking number, $Lk$, is a geometric parameter of a pair of closed curves and although it is well defined in terms of a double integral along them, it has been shown to be an integer equal to half the number of times one curve crosses the other [53].

The total twist, $Tw$, of a pair of curves and the writhe, $Wr$, of one space curve, are given by the following integrals along the corresponding curves [53]:

where $s$ is the arc-length of the curve x, t is its unitary tangent vector and u is a unitary vector orthogonal to t and pointing from the curve x to its parallel curve. All these vector quantities are functions of $s$. The total length of the curve x is simply given by $L=\int_{0}^{L}ds$.

Suppose we initially prepare two space closed curves such as to present a certain amount of $Lk$. The Călugăreanu-White-Fuller theorem guarantees that $Lk$ is always conserved no matter how the curves change along the time, provided they remain closed. Changes to the curves mean changes to their values of $Tw$ and $Wr$ through equations (2) and (3), respectively. According to the theorem, these changes are such that $Tw+Wr$ remains constant as along as the curves remain closed. An interesting feature is that the theorem has also been shown to hold for a pair of non-closed curves if their extremities are flat and belong to the same plane [53]. There exists a pair of parallel open curves with $Tw=Wr=0$ that connect the first two ones at infinity [53]. A similar argument can be made for a pair of open curves having a semi-integer value of $Lk$. As long as the ends of the curves lay on parallel planes, there exists a pair of twisted parallel curves with $Wr=0$ and $Tw=0.5$ that connect the first pair at infinity.

The two space curves required to calculate the twist and writhe of our TGRNs can be defined prior to applying the initial turns to the straight GNR. The first curve can be the GNR axis, and the second curve can be a line parallel to the first. Numerically, both curves can be defined by sets of co-linear carbon atom positions along the main length of the straight untwisted GNR. Once the curves are defined, if one fixes one end of a GNR and applies $n$ turns to the other end while keeping the GNR straight, the twist will be $Tw=n$ and the writhe $Wr=0$. Thus, the initial linking number applied to the, now twisted, GNR is $Lk=n$. If the ends of this TGNR are placed such that they belong to the same plane and are not allowed to rotate back to release the initially applied torsional stress, the value of $Lk$ will remain constant.

As shown by Fonseca [48], if the extremities of a TGNR are laid down on two different planar substrates, and their planes coincide, the $Lk$ theorem can be applied to the TGNR and eqs. (2) and (3) can be used to infer the values of $Tw$ and $Wr$ of the suspended part. Additionally, it is possible to investigate how $Tw$ and $Wr$ vary with several other parameters and physical conditions, such as the distance between the substrates, temperature, etc. The interesting thing is that as long as the extremities of the TGNR are kept on the substrates (and van der Waals forces guarantee that), no matter how $Tw$ and $Wr$ change with other physical conditions, $Lk$ will remain the same. For instance, Fonseca [48] showed that it holds true for changing the distance between substrates and the temperature of the system.

Here, we will investigate how the TC of TGNRs depends on $Lk$, its length, $d$, as well as its $Tw$ and $Wr$ taking into account that these last two quantities change with both $Lk$ and $d$. We will analyze the twist and writhe of non-closed TGNRs to which an integer or semi-integer number of turns was initially applied.

The TC calculation for the TGNRs will be performed by non-equilibrium molecular dynamics (NEMD) simulations using the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) potential [54, 55] as implemented in LAMMPS [56]. AIREBO is an extension of the REBO potential originally developed by Brenner et al. [54], which includes Lennard Jones and torsional potential terms [55]. After more than two decades, the AIREBO potential is still being successfully used to simulate structural [57, 58, 59, 60] and thermal properties [61, 62, 63, 64, 65, 66, 67] of carbon nanostructures, including heat transport simulations [35, 36, 37, 38, 68, 69]. One fundamental aspect of our choice is the computational time involved.

Nonetheless, we must keep in mind that AIREBO does not quantitatively reproduce absolute TC values for carbon nanostructures. In order to overcome this limitation, we will focus on how the TC of TGNRs depends on their geometric features, and not on its absolute value. Zhang and collaborators, for example, have performed a similar study using the original REBO/AIREBO to investigate TC trends in graphene with the number of isotopes [70], and in graphene oxide with the percentage of oxygen coverage [71].

The TC simulation protocol can be described as follows. TGNRs having $Lk=0$ (non-twisted and straight), $0.5,1.0,1.5$ and $2.0$ are generated by fixing one of their ends and applying $2\pi Lk$ rads with respect to the ribbon axis to the opposite end. With the extremities fixed, the TGNRs are optimized with the conjugate gradient energy minimization algorithm as implemented in LAMMPS (with energy and force tolerances of $10^{-8}$ eV and $10^{-8}$ eV/Å, respectively). Then, the TGNR extremities are placed at $\sim 3.4$ Å of distance to two different substrates modeled as large area square shape graphene single layers of $\sim 287$ Å of side, distanced by $d$. The amount of area of the TGNR extremities laid on each substrate is such that its suspended part has one of the following sizes: $d=100,200,300,400$ and $500$ Å. Each TGNR is further equilibrated at 300 K for about 1 ns using a Langevin thermostat, with 0.5 fs as time step and 1 ps as thermostat damping factor. Long time simulations are required in order to guarantee that the suspended part of the TGNR reaches equilibrium. During these simulations, the substrates are kept fixed and the objective of this part of the protocol is to get the equilibrium shape of the TGNR at the chosen 300 K temperature. From the equilibrium shapes of TGNRs, $Tw$ and $Wr$ can be calculated using the algorithm described in Ref. [48].

Armchair GNRs of about 600 Å length and 33 Å width are considered in the present study. They are fully hydrogen-passivated. As classical MD simulations show no special dependence of TC with the direction of thermal conduction in pristine graphene we have not repeated the simulations with zigzag GNRs [35].

The TC was calculated as follows. In a real situation, the substrates play the role of thermal baths. However, the simulations to determine the TC of TGNRs will be performed in the nanoribbons alone, without the substrates, to save time in the thermal equilibration of the substrates and to avoid the unknown heat transmission and thermal resistance at the interface between the substrate and the TGNR. Thermostats at $T_{\mbox{\scriptsize{HOT}}}=350$ K and $T_{\mbox{\scriptsize{COLD}}}=250$ K are, then, applied to the carbon atoms that, in the simulations with substrates, laid on each of them. In the absence of the substrates, a free TGNR would rotate and release its torsional stress. In order to avoid that, the hydrogen atoms that also laid on the substrates are kept fixed during the simulations to determine the TC. The carbon and hydrogen atoms that are suspended in the simulations with the substrates, are allowed to freely evolve, i. e., no thermostat or constraints are applied to them. The simulations to determine the TC of each TGNR were performed for, at least, 40 ns.

Figure 1 depicts three TGNRs with different values of $Lk$ and $d$. There, red and blue atoms are thermostated at $T_{\mbox{\scriptsize{HOT}}}$ and $T_{\mbox{\scriptsize{COLD}}}$ temperatures, respectively, while black atoms are kept fixed. Cyan, white, red and pink atoms at the suspended part of the TGNR are allowed to evolve freely. Red and pink atoms in the suspended part of the TGNR are those whose coordinates will be used to obtain the space curves needed to calculate $Tw$ and $Wr$ using equations (2) (3), respectively. For every TGNR, the TC, $Tw$ and $Wr$ of the suspended part were calculated.

The TC, $\kappa$, of a system along a direction $\mathbf{x}$, can be obtained from Fourier law:

where $J_{\mathbf{x}}$ is the heat flux along $\mathbf{x}$ direction and $\nabla_{\mathbf{x}}\equiv\partial/\partial x$. The heat flux is calculated by the energy per time per cross-sectional area that the thermostats provide to the system. The temperature gradient is calculated by dividing the TGNR in several slabs of length about 10 Å, and determining the local temperature of each slab through the average kinetic energy of the moving atoms over 10000 timesteps every 10000 timesteps. Figure 2 shows a typical temperature profile along the TGNR after the system reaches the steady state.

The functional form the the fitting function by itself is not so important at the moment. Different $d^{n}Lk^{m}$ terms could have been added to the fitting equation with no significant difference in the final result. The point is that it is possible to obtain an empirical analytical function for TC = TC$(d,Lk)$ from computational and/or experimental data and, then, use it for future predictive purposes. Here, Figure 4 serves to reinforce the conclusion that the TC of TGNRs cannot be described in a simplistic manner, solely in terms of the number of initially applied turns or $Lk$.

The TC of the TGNRs can be correlated to their geometric features twist, $Tw$, and writhe, $Wr$, instead of $d$ and $Lk$, because the present simulations were conducted in such a way that the linking number (or the initial number of turns applied to the straight GNR) remained fixed. Then, the Călugăreanu-White-Fuller theorem, eq. (1), can be used to distinguish groups of TC surfaces in a $Tw\times Wr$ space, each one corresponding to a value of $Lk$.

Figure 5 displays four non-zero constant-$Lk$ TC surfaces as function of both $Tw$ and $Wr$ corresponding to the values of $d$ and $Lk$ considered in the present study. The area of the surfaces increases with $Lk$ which reflects the ability of the TGNRs to convert twist to writhe to release at least part of the torsional stress. As the surfaces do not intersect one another, each pair of geometric coordinates, $(Tw,Wr)$ univocally characterizes the TC of a TGNR.

Separated plots of TC versus $Tw$ and versus $Wr$, for different values of suspended length, $d$, are shown in Figure 6. They allow for a better observation of the complexity of the dependence of TC of the TGNRs on their geometric parameters than Figure 5.

Figure 6(a) shows that TGNRs with smaller suspended lengths, $d$, present larger variations of the TC with the real twist, $Tw$. In other words, as the value of $d$ increases, the TC becomes less dependent on $Tw$. This shows that the TGNR TC could not be described uniquely even by the real twist, $Tw$. In fact, the above result is not unexpected. As can be seen in the examples of TGNRs shown in Figure 1, the suspended part of the structure becomes less curved as $d$ increases. This is a consequence of the decrease of the linear twist density with increasing $d$ of the TGNR having the same $Lk$. Literature has shown that rods and ribbons become unstable when the applied twist is such that the twist density becomes larger than a critical value [2, 3, 4]. Because of this, we decided to also investigate the dependence of TC on the twist density, as we will discuss in subsection III.3.

The curves in Figure 6(b) look different from those of panel (a) but they are consistent and reflect the fact that $Tw$ and $Wr$ are, in fact, connected by Eq. (1). In fact, Figure 6(b) shows that for TGNRs with larger suspended length, $d$, there is a larger number of TC points corresponding to $Wr<0.1$. This particular observation is coherent with the fact that increasing $d$ decreases the twist density of the TGNRs. If the twist density is small, the twisted but straight ribbon is a stable spatial conformation. In other words, the structure becomes less curved when the twist density is low. The points corresponding to values of $Wr>0.2$ are those obtained for the largest values of $Lk$ considered in the present study, or $Lk\geq 1.5$. Finally, Figure 6(b) also shows that the TGNR TC cannot be described by only its writhe, $Wr$.

Figures 6(a) and (b) can be re-arranged if we group data points by their $Lk$ values rather than $d$. Indeed, Figure 7 shows TC as a function of $Tw$ (top panel) and $Wr$ (bottom panel), for each value of $Lk$. Now, the curves display additional results about the complex dependence of the TC of TGNRs on their geometric features. It can be seen that the rate of change of TC with either $Tw$ or $Wr$, depends on $Lk$. However, as $Lk=Tw+Wr$ is constant, $dTw=-dWr$ for the curves corresponding to the same values of $Lk$. Thus, for a given value of $Lk$, the rate of change of TC with $Tw$ is equal in modulus to the rate of change of TC with $Wr$. For $Lk=0$, TC does not depend on $Tw$ or $Wr$, since there is no twist and the different values of TC correspond only to different values of $d$. It is equivalent to say that $\Delta\mbox{TC}/\Delta Tw\rightarrow\infty$ if $Lk,Tw,Wr\rightarrow 0$.

However, for $Lk\neq 0$, we observe that TC varies with either $Tw$ or $Wr$. The difference amongst the curves is the rate of change of TC with $Tw$ or $Wr$, $\Delta\mbox{TC}/\Delta Tw$, whose average values are given in Table 1. Although we have only a few values of $\Delta\mbox{TC}/\Delta Tw$, and their values might have large uncertainties, we can infer that they roughly decrease from a large amount and converge with increasing $Lk$.

The above analysis confirms that the TC of TGNRs in realistic conditions (large size GNRs, suspended and at different temperatures) is much more complex than the simple consideration of its dependence on the number of initially applied turns can deal with. The literature has only presented predictions for the dependence of TC on twist, for zero K, non-writhed, small width TGNRs.

The study of the dependence of TC on both $Tw$ and $Wr$ is not practical from the experimental point of view. It is easier to record the number of times the straight GNR was initially rotated and measure the suspended length of the produced TGNR than defining two space curves along the GNR and develop computational tools to extract its points to compute the corresponding $Tw$ and $Wr$. Therefore, although the above results show the TC can be well characterized by the pair $(Tw,Wr)$, they are not unique as shown in subsections III.1 and III.3.

The analyses in the previous sections suggest one more attempt to characterize the TC of TGNRs on just one quantity: the twist density. In fact, there are two possible twist densities to consider, $Lk/d$ and $Tw/d$. We will stick to the first one since, as mentioned in Section III.2, it is easier for an experiment to measure $Lk$ then to measure $Tw$ in a TGNR.

Figure 8 shows the TC of TGNRs as a function of $Lk/d$ for different values of $Lk$. It can be seen that each set of data points corresponding to one value of $Lk$ seems to follow one particular decaying curve. We, then, fitted each one to the following equation:

where $C$ and $\alpha$ are constants that depend on $Lk$, and $x_{0}=Lk/d_{0}$ is the smallest twist density of the set of data points corresponding to the same $Lk$, with $d_{0}$ being the largest suspended size of the TGNRs that, in our case, is 500 Å. The TC of TGNRs is, then, described by the following equation:

Table 2 shows the values of $C$ and $\alpha$ obtained from the fitting of the data points in the main panel of Figure 8, for each value of $Lk$. While $C$ is shown to be weakly dependent on $Lk$, $\alpha$ seems to be a significant function of $Lk$. In fact, the inset of Figure 8 shows that $\alpha(Lk)$ displays an approximately linear decreasing behavior with $Lk$. The meaning of this result is quite interesting. As $C(Lk)\sim const.$, one can conclude that larger values of $Lk$ lead to a weaker dependence of the TC on the twist density $Lk/d$. It is a remarkable, and apparently contradictory, result because larger values of $Lk$ imply larger values of the twist density and, therefore, a stronger dependence of TC on twist density. Thus, one would expect that $\alpha(Lk)$ should increase and not decrease with $Lk$. But the richness of the twist-to-writhe phenomenon can help understand this feature. As shown by Fonseca [48], a TGNR on top of two separated substrates can present a twist-to-writhe transition, where part of its twist is converted to writhe through curving and curling in the space. As this phenomenon decreases the torsional stress on the nanoribbon, it decreases the torsional stress/strain contribution to its TC.

While the dependence of the TC of TGNRs on $(Lk,d)$ is relatively simple, equations similar to (5) are difficult to interpret in physical terms. However, although the dependence on $(Lk,LK/d)$ is a bit more complicated than that on $(Lk,d)$, equation (7) carries a simple and physically meaningful form of describing the TC of TGNRs.