On the Mechanical Properties and Fracture Patterns of Biphenylene-Based Nanotubes: A Reactive Molecular Dynamics Study

: Carbon-based materials have garnered significant attention since the groundbreaking synthesis of graphene. The exploration of novel 2D carbon allotropes has led to the discovery of materials with intrinsic properties distinct from graphene. Within this context, the biphenylene network (BPN) was successfully synthesized. In this study, we used molecular dynamics (MD) simulations with the Reactive Force Field (ReaxFF) to delve into the thermomechanical properties and fracture patterns of biphenylene-based nanotubes (BPN-NTs) exhibiting armchair (AC-BPN-NT) and zigzag (ZZ-BPN-NT) chiralities. Throughout the longitudinal deformation process, we observed significant morphological transformations preceding the structural fracture of the system. These transformations unfolded in distinct inelastic phases. In both cases, AC-and ZZ-BPN-NT, stress accumulation in four-membered rings led to the creation of octagonal structures; however, in AC, this occurs in the fracture region, subsequently causing the presence of nanopores. On the other hand, for ZZ-BPN-NT, stress accumulation in the rectangular rings occurred in bonds parallel to the deformation, with elongated octagonal structures. The Young’s modulus of these nanotubes ranged from 746 to 1259 GPa, with a melting point of around 4000 K. Our results also explore the influence of diameter and curvature, drawing comparisons with BPN monolayers.


Introduction
Carbon exhibits excellent versatility due to its tendency to be present in various hybridizations [1,2].Its remarkable characteristics have sparked massive interest in the scientific community, which shows a growing and conspicuous effort in exploring the rich allotropy of this element to make it the basis for an unprecedented technological revolution [3].In addition to well-known carbon allotropes such as diamond [4], graphite [5], and graphene [6,7], several research groups have successfully reported the synthesis route of new two-dimensional carbon-based materials, biphenylene [8], γ-graphyne [9], and the fullerene network [10].
Carbon-based nanomaterials, such as those mentioned above, have broad applications in various technological fields, including energy storage [11] and conversion [12] devices.They offer a cost-effective alternative that may soon replace conventional silicon-based technology [7].Thus, among the industrial sectors benefiting from the implementation of these materials, it is worth highlighting the organic photovoltaic industry [13], lightemitting diodes [14], and thin-film transistors [15].
In addition to carbon-based materials, other 2D materials exhibit diverse properties with advantages and disadvantages depending on the application [16], formed from different chemical species.For instance, transition metal dichalcogenides (TMDs) can comprise transition metals such as tungsten and molybdenum and chalcogens such as sulfur, selenium, and tellurium [17].These materials may function as semiconductors and share similarities with carbon-based materials, albeit with a relatively larger thickness [18].Additionally, materials based on silicon, nitrogen, and boron display intriguing two-dimensional characteristics of interest [19].
The development of new materials is a process of great theoretical [20] and experimental [21] complexity, capable of posing challenges that test the methods and knowledge considered canonical by the scientific community.Regarding the theoretical aspect of this endeavor, the current scientific literature has a plethora of methodologies that describe atomic and molecular movements with excellent efficiency [22], enabling a detailed analysis of the physicochemical properties of the material under study even before resorting to its laboratory synthesis [23].Among various experimentally obtained characteristics, thermomechanical properties present a considerable range of interest within the scientific community working on carbon-based nanomaterials [24][25][26][27].
Recently, Fan et al. provided a comprehensive account of the synthetic pathway for a novel two-dimensional non-benzenoid carbon allotrope as thin as graphene, known as the Biphenylene Network (BPN) in their publication [8,28].This material is characterized by rings comprising 4, 6, and 8 carbon atoms with sp 2 hybridization, forming a regular lattice.Even before its synthesis, numerous publications explored theoretical predictions concerning its electronic and structural properties [29][30][31].By the first quarter of 2024, nearly three years after the authors reported the synthesis of BPN, the paper has garnered over 400 citations on the Scholar platform.
For instance, Ferguson et al. delved into the potential use of BPN in lithium-based batteries [31], while research by Hosseini, based on density functional theory (DFT) calculations, examined the impact of aluminum doping on the electronic structure and gas detection properties of BPN nanosheets [32].The first theoretical proposal of BPN dates back to 1968 [23].However, post-synthesis, there has been a substantial surge in publications exploring biphenylene's physical and chemical properties when combined with various materials [33][34][35][36][37].
Among these works, studies employing quantum chemistry methods to investigate the adsorption of aromatic organic pollutants on graphene doped with BPN have emerged [38].Additionally, a study utilizing DFT calculations focused on exploring the application of biphenylene nanotubes (BPN-NT) as a CO 2 gas sensor [39].These collective efforts signify the growing interest and diverse applications of BPN in 2D nanomaterials.The mechanical and thermal properties of BPN have already been reported in the literature in the form of monolayers [40].Nevertheless, the existing literature lacks in-depth studies on the thermomechanical properties and fracture patterns of BPN-NT, as well as an exploration of the impact of parameters associated with this one-dimensional topology (nanotubes).The present study seeks to address and bridge this gap in our understanding.
In this work, we conducted fully atomistic reactive (ReaxFF) MD simulations to investigate fracture dynamics and thermomechanical properties of a given set of BPN nanotubes with different chiralities and diameters.Our simulations considered different tube diameters at room temperature.Young's modulus values of BPN-NTs were estimated to range from 746 to 1259 GPa.Upon a longitudinal deformation process, we observed significant morphological transformations preceding the structural fracture of BPN-NTs.These transformations unfolded in distinct inelastic phases.

Methodology 2.1. Biphenylene-Based Nanotube Modeling
The unit cell of BPN comprises six carbon atoms.It is a rectangle measuring 4.26 Å by 3.88 Å, as shown in Figure 1.To form BPN-NTs, we replicate this unit cell vertically and horizontally (see Figure 2).To generate a BPN-NT, a methodology analogous to creating single-walled carbon nanotubes (SWCNTs) [41] is employed.The chiral vector C h is defined as √ 3) = 3.88 Å.These values were obtained using a bond length (a c ) equal to that found for C-C bonds in graphene (1.42 Å) [42], for model simplification, as it has been reported that BPN possesses bond lengths of 1.41, 1.45, and 1.46 Å [43,44].
Creating a nanotube also involves defining the translational vector T = pa 1 + qa 2 (with p and q integers), corresponding to the smallest vector orthogonal to C h (to nanotube unit cell).Given that C h • T = 0, then (na 1 + ma 2 ) • (pa 1 + qa 2 ) = 0, i.e., 3np + (1 + √ 3)mq = 0.The only way to obtain integers p and q is if n = 0 or m = 0. Chiral nanotubes can be achieved with reasonable approximations for the ratio p/q ≈ −(m/n) • 0.91, which is only possible for considerably high integers p, q, m, n.In our approach, with periodic boundary conditions, these situations occur in systems with many atoms.Therefore, this discussion focuses on the (n, 0)-(armchair) and (0, m)-BPN-NT (zigzag) systems.
The BPN-NTs investigated in this study were constructed using the VMD software [45].This process involves rolling BPN sheets along the chiral vectors (n, 0) and (0, m), resulting in nanotubes as illustrated in Figure 1.
The length of the nanotubes is obtained from the translational vector L = |T|.In the case of (n, 0)-BPN-NT, we have T = qa 2 , while the nanotube with the edge (0, m) concerns the translation vector T = pa 1 .Values of p = 24 and q = 26 were chosen, representing the smallest integers that yield a nanotube with a length of at least 100 Å.Through this procedure, we constructed five BPN-NTs of each chirality, with diameters ranging between D ≈ 5 Å and D ≈ 25 Å.Table 1 lists the geometric and structural parameters of the nanotubes analyzed in this study.

Molecular Dynamics Simulations Settings
We employed the LAMMPS software [46,47] to conduct fully atomistic molecular dynamics (MD) simulations, utilizing the reactive force field ReaxFF [48] to delve into the fracture dynamics of BPN-NTs outlined in Table 1.The ReaxFF force field considers bond order, providing an accurate representation of bond breaking and formation [49].Notably, it rectifies the incorrect prediction that two carbon atoms in the system form a triple bond in the rectangular ring, destabilized by terminal radical electrons.As a result, the carbon-carbon bond is no more robust than a double bond.The force field incorporates a partial energy contribution (E C 2 ), presenting the BPN in a stable form [50] to account for the stability of this bond type.Among the potentials tested [51][52][53], the one presented in reference [48], which incorporates the EC2 correction, best reproduced the structural characteristics of BPN.This study has utilized it in the discussion of thermal and mechanical properties.The ReaxFF has been utilized in other studies involving BPN, encompassing mechanical properties [54], wettability [55], and thermal transport [56].
We numerically integrated the equations of motion using the velocity Verlet algorithm [57] with a time step of 0.1 fs.A uniaxial tensile deformation was applied along the length of the nanotubes (axis z), with an engineering strain rate of 1.0 × 10 −6 fs −1 .Before initiating the stretching process of the nanostructures, they were equilibrated in an NPT ensemble at a constant temperature of 300 K, using a Nosé-Hoover thermostat [58] for 200 ps, to prevent the occurrence of residual stresses from system modeling.Before this, all systems were optimized at 0 K.The maximum deformation applied to the structures was 50%, and the values of von Mises (VM) stress [59] were calculated for each atom k through the expression where σ k xx , σ k yy , and σ k zz are the normal stress components, and σ k xy , σ k yz , and σ k zx are the shear stress components.
The computed values of σ k V M aid in locating the point or region of fracture within the structure, providing insights into the fracture processes.From the obtained stress-strain curves, we gather information about the elastic properties of the nanotubes, such as Young's modulus (Y M ), fracture strain (FS), and ultimate strength (US).In addition to Y M , which represents the material's stiffness by relating the stress-strain relationship, we consider the US, the maximum stress value reached by the system.FS represents the percentage of strain relative to the US.Snapshots and molecular dynamics trajectories are also obtained using the visualization and analysis software VMD [45].
Regarding the thermal properties of BPN-NTs, we subjected the systems to gradual heating within an NVT ensemble spanning from 300 to 10,000 K.The temperature increase rate was adjusted to around 200 K/ms.Before initiating the heating process, the nanotubes underwent the equilibration and thermalization procedures mentioned earlier.All nanotubes following this equilibration are available in the Supplementary Material in Protein Data Bank (PDB) format.To measure the phase changes of BPN-NTs as the temperature increased, we recorded the total energy, temperature, and thermal capacity [60].

Mechanical Stability
We explore the correlation between stress (σ) and tensile strain (ε), applying strain along the unit vectors a 2 and a 1 corresponding to chiralities (n, 0) and (0, m), respectively.Figure 3 illustrates the stress-strain relationship of BPN-NTs.Our simulations reveal distinct structural phases as deformation progresses, initiating with a linear (elastic) region and transitioning to atomic rearrangements in the plastic regime.Irrespective of diameter, a new nanotube configuration (phase transition) occurs at around 8% strain for (n, 0) and 12% strain for (0, m), particularly evident in the (n, 0) systems, where a complete bond breakage in one of the four-sided rings is observed.Subsequently, in all cases, the final rupture process forms Linear Atomic Chains (LACs).Complete system rupture is observed in all (n, 0)-BPN-NTs, while in (0, m)-BPN-NTs, separation is not yet evident even at 50% strain.
In Figure 3a, the systems (4,0), (7,0), (11,0), (15,0), and (19,0)-BPN-NT are shown in black, red, blue, purple, and yellow, respectively.These systems, with an edge identical to that seen in armchair carbon nanotubes, exhibit an abrupt transition from the elastic to the plastic regime.The Young's modulus of the systems is higher for nanotubes with a larger diameter.Table 2 presents the numerical values of the elastic properties.Young's modulus (Y M ) was calculated through a linear fit of stress during the first 3% of deformation (within the linear regime).The ultimate strength (US) and fracture strain (FS) were also obtained from the stress-strain curve.
The (4, 0)-BPN-NT system, presented in Figure 3a, has two aspects that make it more susceptible to deformation: the relatively small diameter (≈ 5 Å), which imparts a more significant curvature, and van der Waals interactions between its internal walls, making this system less stable than the others.Despite the relatively different approach with periodic systems, the Young's modulus values presented are consistent with those reported in the literature for 2D BPN, around 1020 GPa, in the equivalent direction [40].In the same perspective, a similar behavior is presented for the systems (0,5), (0,9), (0,14), (0,18), and (0,22)-BPN-NT in Figure 3b, with curves in black, red, blue, purple, and yellow, respectively.As in the armchair case, here, the (0, m)-BPN nanotubes (zigzags) show a behavior not significantly affected by the diameter but with a change from elastic to plastic regimes occurring at considerably smaller deformations, around 5%.The system appears more malleable in this type of edge, with significantly lower maximum stress than the armchair case.One can note that the US of the systems in this chirality revolves around 60 GPa (see Table 2), while the armchair systems are close to 100 GPa.
Still in Figure 3b, the (0, m)-BPN-NT systems undergo a phase transition to another system composed of 8-membered rings (as discussed below in Figure 4).These rings are arranged horizontally in the strain direction, generating a considerable stress accumulation, where later, the system breaks (around 13%) with pores and a constant formation of LACs until the 50% applied deformation.The Young's modulus of the systems is considerably close, around 790 GPa, with only the most minor diameter case, (0,5)-BPN-NT, having Y M = 746 GPa, given its structural characteristics.The equivalent 2D case reports a Young's modulus on the order of 740 GPa [40].
The difference in the US between the (n, 0) and (m, 0) cases lies in the arrangement of the involved rings.BPN comprises systems of 4, 6, and 8-membered rings, with the first two being close to geometrically regular rings.However, the 8-membered ring exhibits considerable asymmetry, which aligns with the deformation direction in the armchair case, making the system stiffer in this direction.Conversely, in the zigzag chirality, this 8membered ring is oriented with its shorter direction aligned with the deformation direction, making it more flexible in this direction, as observed in the fracture pattern discussed in Figure 4.
For a more comprehensive discussion on the dynamics of BPN-NT under tensile loading, we illustrate in Figure 4  In Figure 4a,b, distinctive behaviors of BPN-NTs emerge even in the early stages of the deformation process.For instance, the (11,0)-BPN-NT maintains its structural integrity throughout deformation, enduring up to its ultimate strain percentage (7.6%)before experiencing bond breakages.In Panel Figure 4c, a snapshot at 8.2% reveals the presence of pores in the system.Panel (e) of this figure further illustrates the enlargement of these pores, accompanied by small LACs.This particular feature leads to a temporary drop in stress within the system, which rises again with increasing deformation, elucidating the stress oscillations depicted in Figure 3.In the case of (0,14)-BPN-NT (Figure 4b,d,f), at 13%, after the FS of the system, the snapshot shows the system with several structural failures, including bond breakages in the fracture region.Panels (d) and (f) of this figure depict the system with deformations of 46% and 50%, respectively.The fracture process in this direction is analogous to the two-dimensional system, where atomic reorganization allows the system to undergo a high level of plastic deformation without complete rupture, maintaining stress around 60 GPa throughout the deformation process.
It is worth discussing further that in both cases, some bonds involved in the fourmembered rings are always the first to dissociate.This is due to the sp 2 configuration with a 90 • angle.In all cases, the dissociation of a bond involved in the rectangular ring leads to the formation of octagonal fragments in the system.In the armchair case, this occurs near the fracture region (see Figure 4c,e), while in the zigzag case, this leads to the appearance of these elongated fragments parallel to the deformation direction (see Figure 4d,f).
Characteristics such as fracture patterns and response to heating have already been partially addressed in our current results, as the heating of the intermediate system obtained around 1000 K is not expected to undergo significant changes from what has been discussed here, while the system with aligned octagonal rings in the strain direction, in the case of (0, m) nanotubes, has not fully ruptured even up to the 50% strain applied here, so its characteristics are not yet fully known.Ultimately, despite further studies being necessary to characterize these intermediate structures, they are likely recoverable, as they are part of the plastic deformation region of BPN.
In summary, concerning the fracture pattern and stress-strain response observed for BPN nanotubes, the armchair case exhibits a stress drop with increasing deformation after the FS.This trend occurs due to the breaking of C-C bonds in the system and the formation of nanopores, which lead to a decrease in stress but not to zero immediately.The curve for the zigzag case follows a typical stress-strain curve since the system is not entirely fractured under the total deformation applied.

Thermal Stability
Finally, we analyze the thermal stability of BPN-NT.The system temperature and total energy as a function of time are standard output quantities provided by LAMMPS.Using a post-processing code, we calculate the derivative of total energy concerning the temperature of normalized data, thus obtaining the C V .We conducted simulations of gradual heating (melting process) in a temperature range from 300 to 10,000 K over 200 ps. Figure 5 illustrates the relationship between total energy (represented in red) and heat capacity (C V , represented in cyan), which in this study was obtained by the energy variation concerning temperature in a constant-volume system [60], during the temperature increase process from a canonical ensemble.This figure shows that the total energy increases almost linearly concerning temperature, with three clearly defined regimes indicated by the observed slopes: 300-1000 K, 1100-3000 K, and 5000-10,000 K.
The first heating stage (300-1000 K) is characterized by the linear increase of the total power curve with temperature, starting from room temperature until the first peak in the C V curve.During this phase, the BPN-NT's structural integrity remains intact.However, as the temperature surpasses this heating range, thermal vibrations cause significant deformations in the original morphology of the nanotube, resulting in a structure similar to that presented in Figure 4f.This trend indicates that this new system is more stable than the pristine nanotube under extreme conditions, such as increased temperature and deformation.
The second heating stage (1100-3000 K) corresponds to the continuous heating of the newly formed structure, leading to the transition from the solid to the gas phase.This phenomenon is represented by the steep slope in the total energy curve and the formation of a prominent peak in the C V curve, indicating the melting point of BPN-NT at 4000 K.It is noteworthy that this value is comparable to the melting point of a graphene monolayer (4095 K), amorphous graphene monolayer (3626 K) [61], and BPN monolayer (4024 K) [40].Previous studies also predicted a melting point for graphene from 4000 K to 6000 K [62][63][64].Like BPN monolayers, the melting of BPN-NT occurs at an approximate temperature of 4000 K [40].Before reaching this critical temperature, the total energy curve does not show significant changes in its slope.However, in the range between 4100 K and 5000 K, we observe a drastic change in the slope of the curve, representing the increase in kinetic energy due to higher particle velocities in the gas phase.Additionally, the vibrational, oscillatory, and torsional energies in the solid phase convert into kinetic energy throughout the melting process, contributing to the increase in total power.Temperatures above 5000 K characterize the third heating regime where the gas phase predominates.
In Figure 6, we present a series of snapshots representing stages of a heating ramp simulation, with a temperature variation ranging from 300 K to 5000 K.We used a color scheme to represent the temperature per atom, with white representing 300 K and red representing 5000 K, based on the equipartition theorem.At temperatures of 1000 K and 2000 K (Figures 6b and 6c, respectively), thermal fluctuations cause changes in morphology, resulting in structures similar to those presented in Figures 4b,c.Thus, we can see that temperature and deformation similarly affect structural changes in low-temperature regimes.At temperatures near 3000 K, the melting process of the structure begins (Figure 6d), and above 4000 K (Figure 6e), thermal vibrations have acquired sufficient intensity to allow the formation of only a few LACs in an amorphous structure, with no resemblance to pristine BPN-NT.Above this temperature, a complete structural transition to a gaseous phase occurs, with continuous breaking of the LACs in the system.

Summary and Conclusions
In summary, we employed reactive molecular dynamics (MD) simulations to explore the fracture dynamics and thermomechanical stability of biphenylene nanotubes with chiralities (n, 0) and (0, m).Our findings uncovered a sequence of morphological transformations in BPN-NTs during uniaxial tensile loading, culminating in a complete fracture.As the deformation progresses, distinct inelastic phases emerge.Upon reaching the critical deformation point initiating the fracture process, nanopores form in the system, followed by the development of Linear Atomic Chains (LACs) until the structure experiences complete rupture.
Similar to BPN sheets, BPN-NTs also exhibit different stages of deformation, which are related to the formation of other morphologies.During the structural deformation process, we observe that in (n,0)-BPN-NTs, stress accumulation in the early stages of stretching leads to sudden failure of the structures, inducing a structural conformation conducive to the formation of LACs that resist high degrees of stretching until reaching complete rupture.Regarding (0,m)-BPN-NTs, stress accumulation during the stretching process induces four and six-atom rings to convert into elongated octagonal structures in the direction of the applied deformation; these systems do not withstand a high tensile value but show high deformation.
Through BPN-NT stretching simulations, we obtained Young's modulus values ranging between 746.26 and 1259.92GPa.The values for (n, 0)-BPN-NTs are considerably higher than those for zigzag chirality (0, m); these values are in the same order of magnitude as Young's modulus of graphene nanotubes, approximately 1 TPa [65], but less resistant in terms of their elastic regime.Regarding the heating simulations, we observed that the melting point of BPN-NT occurs at 4000 K, a value very close to the melting point of the BPN sheet (4024 K).
Theoretically, using quantum techniques, the new crystalline systems formed upon lattice stretching can be characterized by their structural and electronic properties.As biphenylene is already synthesized (in a 2D form [8]), deformation techniques can be applied to verify these phase transitions experimentally.Once this initial challenge is overcome, various experimental methods such as X-ray diffraction (XRD), transmission electron microscopy (TEM), and scanning electron microscopy (SEM) help us to understand the atomic arrangement and crystallographic structure of these new material phases.

Figure 1 .
Figure 1.(a) Schematic representation of the BPN sheet and its unit cell (a), (b) zoomed-in region of the structure containing the unit cell with lengths of the unit vectors, (c) and side view of the structure.

Figure 2 .
Figure 2. (a) Side and front views of a (n,0)-BPN-NT and (b) (0,m)-BPN-NT, where the region in red highlights the chirality of the structure.

where a 1
and a 2 are orthogonal vectors, and m and n are integers.In this context, |a 1 | = 3a c = 4.26 Å and |a 2 | = a c (1 +

Figure 3 .
Figure 3. Stress-strain curves for uniaxial deformations applied in the longitudinal direction of nanotubes in armchair (a) and zigzag (b) chiralities, for different diameters.
a sequence of representative MD snapshots for uniaxial deformation applied along the longitudinal axis of (11,0)-BPN-NTs (panels (a), (c), and (e)) and (0,14)-BPN-NTs (panels (b), (d), and (f)), displaying the von Mises stress values (σ k V M ) [33,58].The σ k V M values provide relevant information about the local structure during fracture, making it possible to identify regions of high-stress concentration and predict where structural failures are more likely to occur.

Figure 4 .
Figure 4. Representative MD snapshots with von Mises stress values σ k V M per atom for (11,0)-(illustrated in panels (a,c,e)), and (0,14)-BPN-NT (illustrated in panels (b,d,f)) under the effect of longitudinal stress.The upper panels (a,b) show the transition from the elastic to the plastic regime.The middle panels (c,d) represent the deformation undergone by the structures immediately before the first fracture.Finally, panels (e,f) illustrate the moment of the first structural failure.

Figure 5 .
Figure 5.Total energy (E) and heat capacity (C V ) for BPN-NT as a function of temperature during the heating process.

Figure 6 .
Figure 6.Representative MD snapshots for the heating ramp simulations (melting process) for 300 K at 0 ps (a), 1000 K at 20 ps (b), 2000 K at 40 ps (c), 3000 K at 60 ps (d), 4000 K at 80 ps (e), and 5000 K at 100 ps (f).The temperature per atom varies from 300 to 5000 K, represented by the colors in the BWR scale.

Table 1 .
Geometrical and structural parameters for the modeled BPN-NTs.