A First Principles Investigation of the Mechanical Properties of g-TlN

We investigate the structure and mechanical properties of proposed graphene-like hexagonal thallium nitride monolayer (g-TlN) using first-principles calculations based on density-functional theory. Compared to graphene-like hexagonal boron nitride monolayer (g-BN), g-TlN is much softer, with 12% in-plane stiffness, 25%, 22%, and 20% ultimate strengths in armchair, zigzag, and biaxial strains respectively. However, g-TlN has a larger Poisson’s ratio, 0.69, about 3.1 times that of g-BN. It was found that the g-TlN also sustains much smaller strains before rupture. We obtained the second, third, fourth, and fifth order elastic constants for a rigorous continuum description of the elastic response of g-TlN. The second order elastic constants, including in-plane stiffness, are predicted to monotonically increase with pressure while the Poisson’s ratio monotonically decreases with increasing pressure.

Hexagonal thallium nitride monolayer (g-TlN) is a proposed graphene-like 2D material, which is only monoatomically thick (Figure 1).Although g-TlN has not been fabricated currently, the theoretical study of the g-TlN could expand the range of possible applications of III-nitrides, and open new perspectives for miniaturization in engineering functional nano-devices and interconnects by a chemical modification.The bulk thallium nitride (TlN) was predicted to have a small energy gap, indicating a semi-metallic character [40,41].The combi- nation of thallium with other group III atoms which has wide gap in III-nitrides yields promising semiconductors for optical communication systems (laser diodes, detectors) with small band gap, down to the infrared energy region [42][43][44][45][46][47].Structural stability [48], defects [49], phonon and elastic instabilities under pressure [50] were well studied in bulk TlN.The structural and reactivity parameters of novel N 12 X 12 H 12 (X = B, Al, Ga, In, Tl) nitrides, in their coronene-like (C 24 H 12 ) structure was investigated [51], where Tl-N bond length is reported as 2.15 -2.20 Å.
Mechanical properties are critical in designing parts or structures with g-TlN regarding the practical applications.Strain engineering is a common and important approach to tailor the functional and structural properties of the nanomaterials [52].One can expect that the properties of g-TlN will be affected by applied strain too.In addition, g-TlN is vulnerable to be strained with or without intent because of the monoatomic thickness.For example, there are strains because of the mismatch of lattices constants or surface corrugation with substrate [53,54].Therefore, the knowledge of mechanical properties of g-TlN is highly desired.
Depending on the loading, the mechanical properties are divided into four strain domains: linear elastic, nonlinear elastic, plastic, and fracture.Materials in first two strain domains are reversible, i.e., they can restore to equilibrium status after the release of the loads.On the contrary, the last two domains are non-reversible.Defects are nucleated and accumulated with the increase of the strain, until rupture.As in graphene, the nonlinear mechanical properties are prominent since it remained elastic until the intrinsic strength was reached [55,56].Thus it is of great interest to examine the nonlinear elastic properties of g-TlN, which is necessary to understand the strength and reliability of structures and devices made of g-TlN.
Several previous studies have shown that 2D monolayers present a large nonlinear elastic deformation during the tensile strain up to the ultimate strength of the material, followed by a strain softening until fracture [33,[56][57][58][59].We expect that the g-TlN behaves in a similar manner.Under large deformation, the strain energy density needs to be expanded as a function of strain in a Taylor series to include quadratic and higher order terms.The higher order terms account for both nonlinearity and strain softening of the elastic deformation.They can also express other anharmonic properties of 2D nanostructures including phenomena such as thermal expansion, phonon-phonon interaction, etc. [55].
The goal of this paper is to study the mechanical behaviors of g-TlN at large strains and find an accurate continuum description of the elastic properties from abinitio density functional theory calculations.The total energies of the system, forces on each atom, and stresses on the simulation boxes are directly obtained from DFT calculations.The response of g-TlN under the nonlinear deformation and fracture is studied, including ultimate strength and ultimate strain.The high order elastic con-stants are obtained by fitting the stress-strain curves to analytical stress-strain relationships that belong to the continuum formulation [58].We compared this proposed new material with the well known 2D materials such as g-BN, graphene, and graphyne.Based on our result of the high order elastic constants, the pressure dependence properties, such as sound velocities and the second order elastic constants, including the in-plane stiffness, are predicted.Our results for the continuum formulation could also be useful in finite element modeling of the multiscale calculations for mechanical properties of g-TlN at the continuum level.The remainder of the paper is organized as follows.Section 2 presents the computational details of DFT calculations.The results and analysis are in Section 3, followed by conclusions in Section 4.

Density Functional Theory Calculations
We consider a conventional unit cell containing 6 atoms (3 thallium atoms and 3 nitrogen atoms) with periodic boundary conditions (Figure 2).The 6-atom conventional unit cell is chosen to capture the "soft mode", which is a particular normal mode exhibiting an anomalous reduction in its characteristic frequency and leading to mechanical instability.This soft mode is a key factor in limiting the strength of monolayer materials can only be captured in unit cells with hexagonal rings [60].
The total energies of the system, forces on each atoms, stresses, and stress-strain relationships of g-TlN under the desired deformation configurations are characterized via first-principles calculations with density-functional theory (DFT).DFT calculations were carried out with the Vienna Ab-initio Simulation Package (VASP) [61][62][63][64] which is based on the Kohn-Sham Density Functional Theory (KS-DFT) [65,66] with the generalized gradient approximations as parameterized by Perdew, Burke, and Ernzerhof (PBE) for exchange-correlation functions [67].The electrons explicitly included in the calculations are the 5d 10 6s 2 6p 1 electrons for thallium atoms and the 1s 2 2s 2 2p 1 electrons for nitrogen atoms.The core electrons of thallium are replaced by the projector augmented wave (PAW) and pseudo-potential approach [68,69].A plane-wave cutoff of 600eV is used in all the calculations.The calculations are performed at zero temperature.The criterion to stop the relaxation of the electronic degrees of freedom is set by total energy change to be smaller than 10 -6 eV.The optimized atomic geometry was achieved through minimizing Hellmann-Feynman forces acting on each atom until the maximum forces on the ions were smaller than 0.001 eV/Å.
The atomic structures of both deformed and un-deformed configurations are obtained by fully relaxing a 6 atom-unit cell with all atoms placed in one plane.The simulation invokes periodic boundary conditions for the two in-plane directions while the displacement to out ofplane direction is forbidden.
The irreducible Brillouin Zone was sampled with a Gamma-centered 17 × 17 × 1 k-mesh.Such a large kmesh was used to reduce the numerical errors caused by the strain of the systems.The initial charge densities were taken as a superposition of atomic charge densities.There was a 15 thick vacuum region to reduce the inter-layer interaction to model the single layer system.To eliminate the artificial effect of the out-of-plane thickness of the simulation box on the stress, we used the second Piola-Kirchhoff stress [33] to express the 2D forces per length with units of N/m.
For a general deformation state, the number of independent components of the second, third, fourth and fifth order elastic tensors are 21, 56, 126, and 252 respectively.However, there are only fourteen independent elastic constants need to be explicitly considered due to the symmetries of the atomic lattice point group D 6h which consists of a six-fold rotational axis and six mirror planes [56].
The fourteen independent elastic constants of g-TlN are determined by a least-squares fit to the stress-strain results from DFT based first-principles studies in two steps, detailed in our previous work [33].In the first step, we use a least-squares fit of five stress-strain responses.Five relationships between stress and strain are necessary because there are five independent fifth-order elastic constants (FFOEC).We obtain the stress-strain relationships by simulating the following deformation states: uni-axial strain in the zigzag direction (zigzag); uni-axial strain in the armchair direction (armchair); and equibiaxial strain (biaxial).From the first step, the components of the second-order elastic constants (SOEC), the thirdorder elastic constants (TOEC), and the fourth order elastic constants (FOEC) are over-determined (i.e., the number of linearly independent variables are greater than the number of constrains), and the fifth-order elastic constants are well-determined (the number of linearly independent variables are equal to the number of constrains).Under such circumstances, the second step is needed: least-square solution to these over-and well-determined linear Equations.

Atomic Structure
We first optimize the equilibrium lattice constant for g-TlN.The total energy as a function of lattice spacing is obtained by specifying nine lattice constants varying from 3.3 Å to 4.1 Å, with full relaxations of all the atoms.A least-square fit of the energies versus lattice constants with a fourth-order polynomial function yields the equilibrium lattice constant as a = 3.731 Å.The most energetically favorable structure is set as the strain-free structure in this study and the atomic structure, as well as the conventional cell is shown in Figure 2. Specifically, the bond length of Tl-N bond is 2.154 Å, which is 0.704 Å (or 49%) longer than the bond length of B-N bond in g-BN.The N-Tl-N and Tl-N-Tl angles are 120˚ and all atoms are within one plane.Our result of bond length is in good agreement with previous DFT calculations of N 12 Tl 12 H 12 .

Strain Energy
When the strains are applied, all the atoms are allowed full freedom of motion within their plane.A quasi-Newton algorithm is used to relax all atoms into equilibrium positions within the deformed unit cell that yields the minimum total energy for the imposed strain state of the super cell.
Both compression and tension are considered with Lagrangian strains ranging from −0.1 to 0.3 with an increment of 0.01 in each step for all three deformation modes.We define strain energy per atom E s = (E tot − E 0 )/n, where E tot is the total energy of the strained system, E 0 is the total energy of the strain-free system, and n = 6 is the number of atoms in the unit cell.This size-independent quantity is used for the comparison between different systems.Figure 3 shows the E s of g-TlN as a function of strain in uni-axial armchair, uni-axial zigzag, and biaxial deformation.E s is seen to be anisotropic with strain direction.E s is non-symmetrical for compression (η < 0) and tension (η > 0) for all three modes.This non-symmetry indicates the anharmonicity of the g-TlN structures.The harmonic region where the Es is a quadratic function of applied strain can be taken between −0.02 < η < 0.02.The stresses, derivatives of the strain energies, are linearly increasing with the increase of the applied strains in the harmonic region.The anharmonic region is the range of strain where the linear stress-strain relationship is invalid and higher order terms are not negligible.With even larger loading of strains, the systems will undergo irreversible structural changes, and the systems are in plastic region where they may fail.The maximum strain in the anharmonic region is the critical strain.The critical strain is 0.2 under armchair deformation.However, for other two directions, the critical strains are not spotted.The ultimate strains are determined as the corresponding strain of the ultimate stress, which is the maxima of the stress-strain curve, as discussed in the following section.

Stress-Strain Curves
The second P-K stress versus Lagrangian strain relationship for uni-axial strains along the armchair and zigzag directions, as well as biaxial strains are shown in Figure 4.The stresses are the derivatives of the strain energies with respect to the strains.The ultimate strength is the maximum stress that a material can withstand while being stretched, and the corresponding strain is the ultimate strain.Under ideal conditions, the critical strain is larger than the ultimate strain.The systems of perfect g-TlN under strains beyond the ultimate strains are in a metastable state, which can be easily destroyed by long wavelength perturbations, vacancy defects, as well as high temperature effects [70].The ultimate strain is determined by the intrinsic bonding strengths and acts as a lower limit of the critical strain.Thus it has a practical meaning in considering for its applications.
The ultimate strengths and strains corresponding to the different strain conditions are in Table 1, compared with that of g-BN, graphene, and graphyne.The material behaves in an asymmetric manner with respect to compressive and tensile strains.With increasing strains, the Tl-N bonds are stretched and eventually rupture.When the strain is applied in the armchair direction, the bonds of those parallel with this direction are more severely stretched than those in other directions.The ultimate strain in armchair deformation is 0.17, smaller than that of g-BN, graphene, and graphyne.The critical strain is 0.2 under armchair deformation, where there is big drop of the stresses (Figure 4 top panel), indicating the failure of the system.However, for other two directions, the critical strains are not spotted.Under the zigzag deformation, in which the strain is applied perpendicular to the armchair, there is no bond parallel to this direction.The bonds incline to the zigzag direction with an angle of 30˚ are more severely stretched than those in the armchair numerical methods such as the finite element technique.direction.The ultimate strain in this zigzag deformation is 0.21, smaller than that of g-BN and graphene, while the same as graphyne.At this ultimate strain, the bonds that are at an incline to the armchair direction appear to be ruptured (Figure 4 middle panel).Under the biaxial deformation, the ultimate strain is b m  = 0.16, which is the smallest among those of g-BN, graphene, and graphene.As such ultimate strain applied, all the Tl-N bonds are observed to be ruptured (Figure 4 bottom panel).
The second order elastic constants model the linear elastic response.The higher (>2) order elastic constants are important to characterize the nonlinear elastic response of g-TlN using a continuum description.These can be obtained using a least squares fit of the DFT data and are reported in Table 2. Corresponding values for graphene are also shown.
The in-plane Young's modulus Y s and Poison's ratio may be obtained from the following relationships: Y s = It should be noted that the softening of the perfect g-TlN under strains beyond the ultimate strains only occur for ideal conditions.The systems under this circumstance are in a meta-stable state, which can be easily destroyed by long wavelength perturbations and vacancy defects, as well as high temperature effects, and enter a plastic state [70].Thus only the data within the ultimate strain has physical meaning and was used in determining the high order elastic constants in the following subsection.

Elastic Constants
The elastic constants are critical parameters in finite element analysis models for mechanical properties of materials.Our results of these elastic constants provide an accurate continuum description of the elastic properties of g-TlN from ab-initio density functional theory calculations.They are suitable for incorporation into Y s a ry small to g-BN ( aphene and graphyne (21%).The reduction of in-plane stiffness from g-BN to g-TlN is a result of the weakened bond of Tl-N compared to the B-N bond in g-BN.While all other things being equal, bond length is inversely related to bond strength and the bond dissociation energy, as a stronger bond will be shorter.Considering the bond length, in g-TlN the bond length of Tl-N is 2.154 Å, about 49 percents larger than B-N bond length in g-BN (1.45 Å).The bonds can be viewed as being stretched in prior by the introducing of thallium atoms, in reference to g-BN.These stretched bonds are weaker than those un-stretched, leading to a reduction of the mechanical strength.
Knowledge of higher order elastic constants is very useful in understanding the anharmonicity.Using the higher order elastic continuum description, one can calculate the stress and deformation state under uni-axial stress, rather than uni-axial strain [56].Explicitly, when pressure is applied, the pressure dependent second-order elastic moduli can be obtained from the high order elastic continuum description [32,59,72,73].The third-order elastic constants are important in understanding the nonlinear elasticity of materials such as changes in acoustic velocities due to finite strain.As a consequence, the nano devices such as nano surface acoustic wave sensors and nano waveguides could be synthesized by introducing local strain [33,59].
Stress-strain curves in the previous section show that they will soften when the strain is larger than the ultimate strain.From the view of electron bonding, this is due to the bond weakening and breaking.This softening behavior is determined by the TOECs and FFOECs in the continuum aspect.The negative values of TOECs and FFOECs ensure the softening of g-TlN monolayer under large strain.
The hydrostatic terms (C , C , C , C , and so on) of g-TlN monolayers are smaller than those of g-BN and graphene, consistent with the conclusion that the g-TlN is "softer".The shear terms (C 12 , C 112 , C 1122 , etc.) in general are smaller than those of g-BN and graphene, which contributes to its high compressibility.Compared to graphene, graphyne, and g-BN, one can conclude that the mechanical behavior of g-TlN is similar to graphyne, and much softer than graphene and g-BN.

Pressure Effect on the Elastic Moduli
With third-order elastic moduli, we can study th of the second-order elastic moduli on the pressure P acting in the plane of g-TlN.Explicitly, C when pressure is applied, the pressure dependent second-order elastic moduli ( 11 The second-order elastic moduli of increase linearly with the applied pre H In summary, we studied the mechanical response of s using DFT based first princi-g-TlN are seen to ssure (Figure 5).owever, Poisson's ratio decreases monotonically with the increase of pressure.only occurs when the pressure is zero.This anisotropy could be the outcome of anharmonicity.

Conclusions
g-TlN under various strain ples calculations.It is observed that g-TlN exhibits a nonlinear elastic deformation up to an ultimate strain, which is 0.17, 0.21, and 0.16 for armchair, zigzag, and biaxial directions, respectively.The deformation and failure behavior and the ultimate strength are anisotropic.It has low in-plane stiffness (34.5 N/m) and large Poisson ratio compared to g-BN and graphene.Compared to g-BN, g-TlN has 12% in-plane stiffness, 25%, 22%, and 20% ultimate strengths in armchair, zigzag, and biaxial strains respectively, and a 3.1 times of Poisson's ratio.We also found that the g-TlN can sustain much smaller strains before the rupture.
The nonlinear elasticity of g-TlN was investigated.We found an accurate continuum description of the elastic properties of g-TlN by explicitly determining the fourteen independent components of high order (up to fifth order) elastic constants from the fitting of the stress- strain curves obtained from DFT calculations.This data is useful to develop a continuum description which is suitable for incorporation into a finite element analysis model for its applications in large scale.The second order elastic constants including in-plane stiffness are predicted to monotonically increase with pressure while Poisson's ratio monotonically decreases with increasing pressure.

Figure 2 .
Figure 2. Atomic structure of g-TlN in the conventional unit cell (6 atoms) in the un-deformed reference configuration.

Figure 3 .
Figure 3. Energy-strain responses for uni-axial strain in armchair and zigzag directions, and biaxial strains.
and v =12 11   C C .We have Y s = 34.5 (N/m) and v = 0.689.The in-plane stiffness of g-TlN is

Figure 5 .
Figure 5. Second-order elastic moduli and Poisson ratio as a function of the pressure.