A method for calculating the thermal conductivity of black phosphorene based on external electric field regulation

Abstract

The present invention provides a method for calculating the thermal conductivity of black phosphorene based on the regulation of an external electric field. The method comprises: constructing a computational model of two-dimensional black phosphorene; adjusting various parameters of the computational model to optimize the structure of the computational model, and calculating various parameter values ​​of the optimized computational model; applying several external electric fields of different strengths to the optimized computational model, and returning to the step of constructing the computational model to obtain the computational model of two-dimensional black phosphorene under different working conditions; obtaining the harmonic force constant and the non-harmonic force constant of two-dimensional black phosphorene by compressed sensing lattice dynamics; solving the phonon Boltzmann transport equation by an iterative method, obtaining the regulation of the anisotropy of thermal conductivity of black phosphorene by several external electric fields of different strengths and calculating the corresponding value of thermal conductivity of black phosphorene. The present invention realizes the reversible regulation of the anisotropy of thermal conductivity of black phosphorene by an external electric field, and reveals the mechanism of the anisotropy of thermal conductivity of black phosphorene by electric field from the phonon scale.

Description

Technical Field

[0001] The present invention relates to the technical field of solid material thermal conductivity regulation, and in particular to a method for calculating the thermal conductivity of black phosphorene based on external electric field regulation. Background Art

[0002] As a new type of two-dimensional semiconductor material, black phosphorene has high carrier mobility, adjustable band gap, excellent electrical and optical properties, and high thermal stability, and can achieve sensitive thermal conductivity regulation under the action of an electric field. At the same time, the layered structure of black phosphorene gives it significant thermal conductivity anisotropy in the plane, which shows very important application potential in the fields of nanoelectronics, optoelectronics, and thermoelectrics.

[0003] In the past few decades, people have tried countless methods to regulate the thermal conductivity of solid materials. However, most methods, such as doping, defects, and deformation, are irreversible. Once the material is synthesized, the thermal conductivity is almost impossible to change, and they have the disadvantages of high cost and long time. Therefore, we urgently need a reversible and active method to control thermal conductivity. Summary of the invention

[0004] The purpose of the present invention is to provide a method for calculating the thermal conductivity of black phosphorene based on external electric field regulation in view of the deficiencies of the prior art. The method does not require changing the atomic structure and has high stability and reversibility. The method can achieve reversible regulation of the anisotropy of the thermal conductivity of black phosphorene by an external electric field, and reveal the mechanism of the electric field on the anisotropy of the thermal conductivity of two-dimensional black phosphorene from the phonon scale.

[0005] In order to achieve the above-mentioned purpose, the technical solution of the present invention provides a method for calculating the thermal conductivity of black phosphorene based on external electric field regulation, the method comprising: step S1: constructing a calculation model of two-dimensional black phosphorene; step S2: adjusting various parameters of the calculation model of two-dimensional black phosphorene to optimize the structure of the calculation model of two-dimensional black phosphorene, and calculating the values ​​of various parameters of the optimized calculation model of two-dimensional black phosphorene; step S3: applying several external electric fields of different intensities to the optimized calculation model of two-dimensional black phosphorene, and returning to step S2 to obtain the calculation model of two-dimensional black phosphorene under different working conditions; step S4: obtaining the harmonic force constant and the non-harmonic force constant of two-dimensional black phosphorene by compressed sensing lattice dynamics; step S5: using an iterative method to solve the phonon Boltzmann transport equation, obtain the regulation of the anisotropy of the thermal conductivity of black phosphorene by the external electric fields of several different intensities, and calculate the corresponding value of the thermal conductivity of black phosphorene.

[0006] Preferably, the step S1 specifically comprises: constructing a computational model of two-dimensional black phosphorene by using crystal structure visualization software VESTA.

[0007] Preferably, the various parameters of the two-dimensional black phosphorene calculation model include: k-space grid, energy convergence error, force convergence error, and electron convergence error; and the values ​​of various parameters of the optimized two-dimensional black phosphorene calculation model are calculated by Quantum Espresso software.

[0008] Preferably, the step S4 specifically comprises: step S4.1: obtaining a Taylor expansion of the total energy in terms of atomic displacement power; the expression of the Taylor expansion is:

[0009]

[0010] In formula (1), V is the total energy; V0 is the energy in the reference state; Φ a is the displacement u from atom a a The associated first-order force constant; Φ ab is the displacement u from atoms a and b a and u b The associated second-order force constant; Φ abc is the displacement u from atoms a, b, and c a 、u b and u c The relevant third-order force constant; u a is the displacement of atom a in Cartesian direction i; u b is the displacement of atom b in Cartesian direction i; u c is the displacement of atom c in Cartesian direction i;

[0011] Step S4.2: Given the training data, the expansion coefficient is automatically selected by compressed sensing lattice dynamics, and the value of the expansion coefficient is determined by applying mathematical techniques. Then, formula (1) can be written as a force-displacement relationship expression; the force-displacement relationship expression is:

[0012] F a =-Φ a -Φ ab u b -Φ abc u b u c / 2-…(2)

[0013] In formula (2), F a is the force acting on atom a; the F a obtained by calculation from first principles using any common density functional theory code, and the F a A set of L atomic configurations used in a supercell to achieve a linear relationship for the unknown force constant tensor Φ; the expression of the linear relationship is:

[0014]

[0015] In formula (3), F is the force vector; Φ is the unknown force constant vector, Φ includes all high-order force constants; is the sensing matrix, The elements of are the products of the atomic displacements, The elements of correspond to the different terms in the force expansion, The specific expression is:

[0016]

[0017] In formula (4), is the displacement of atom b in the first configuration, where the subscript b denotes a specific atom and the superscript 1 denotes the first configuration; is the product of the displacements of atoms b and c in the first configuration; Each row in corresponds to a force component, The expression for the total number of rows M is:

[0018] M=3LN at (5)

[0019] In formula (5), L is the number of atoms in the supercell; N at is the number of atoms in the supercell;

[0020] Step S4.3: Minimizing coefficients based on lattice dynamics via compressive sensing The norm of ||Φ||1≡∑ I |Φ I |To solve formula (3); where, The norm of The approximate value of The norm represents the number of non-zero elements in a vector, Φ I | is the absolute value of the force constant tensor element, and I is the index of the force constant tensor; solving formula (3) yields:

[0021]

[0022] In formula (6), Φ CS represents the force constant tensor optimized by the compressed sensing method; argmin Φ ||Φ||1 means finding the force constant value that minimizes the function ||Φ||1, where ||Φ||1 represents the value of Φ. norm; represents the sum of squares of error terms, μ is a weight parameter; represents the sum of the absolute values ​​of all force constant tensor elements; represents the sum of squares of all error terms.

[0023] Preferably, the step S5 specifically comprises: step S5.1: when the only scattering source is a two-phonon and a three-phonon process, a linearized Boltzmann transport equation is obtained; the expression of the linearized Boltzmann transport equation is:

[0024]

[0025] In formula (7), F λ represents the deviation of the corrected phonon distribution function; is the relaxation time of the λ mode, Obtained from perturbation theory; ν λ represents the group velocity of the phonon mode λ; Δ λ is a measure of how much a particular phonon mode deviates from the relaxation time approximation prediction, Δ λ Take speed as the dimension;

[0026] Among them, the Δ λ The calculation expression is as follows:

[0027]

[0028] Among them, the The calculation expression is as follows:

[0029]

[0030] In formula (8) and formula (9), N is the total number of discrete grid points in the Brillouin zone; ξ λλ” is the frequency ratio; F λ” is the driving function of the phonon mode λ"; Γ λλ' is the two-phonon scattering rate; and are the three-phonon scattering rates in the absorption process and the emission process, respectively; is the sum of the three-phonon scattering rates in the absorption process; is the sum of the three-phonon scattering rates during the emission process; is the sum of the two-phonon scattering rates during the emission process;

[0031] Step S5.2: Introduce the Brillouin zone discretization into a Γ-centered regular grid of q points N = N1 × N2 × N3, and set ξ λλ' =ω λ' / ω λ ; Among them, ω λ' is the angular frequency of the phonon mode λ', ω λ is the angular frequency of the phonon mode λ; N1, N2 and N3 are the numbers of discrete points along three different directions of the lattice respectively;

[0032] Step S5.3: Normalized characteristic function e based on three phonons p,q and the anharmonic force constant tensor The three-phonon scattering rates of the absorption process are obtained respectively Three-phonon scattering rate in the emission process and the scattering matrix elements for the three-phonon scattering process Expression of

[0033] Among them, the anharmonic force constant tensor The specific expression is:

[0034]

[0035] In formula (10), E is the total energy of the system; is the displacement component of the i-th atom in the α direction; is the displacement component of the jth atom in the β direction; is the displacement component of the kth atom in the γ direction;

[0036] Among them, the The specific expression is:

[0037]

[0038] Among them, the The specific expression is:

[0039]

[0040] In formula (11) and formula (12), is Planck's constant; f0' and f0" represent the equilibrium distribution functions of phonon modes λ and λ", respectively; ω λ ,ω λ' and ω λ" denote the angular frequencies of the phonon modes λ, λ' and λ", respectively; and are the three-phonon scattering matrix elements in the absorption and emission processes respectively; δ is the Dirac function;

[0041] Among them, the The specific expression is:

[0042]

[0043] In formula (13), M i 、M j and M k denote the masses of atoms i, j, and k respectively; represents the α component of the eigenfunction of mode λ at atom i; represents the β component of the eigenfunction of mode p' at atom j, where ±q' depends on the specific absorption or emission process; represents the γ component of the eigenfunction of mode p" at atom k, where -q" represents the inverse wave vector;

[0044] Step S5.4: Using an iterative method to solve formula (7) to obtain an expression for the thermal conductivity of black phosphorene; wherein the convergence criterion used in the iterative method is to make the relative change of the conductivity tensor calculated by formula (13) less than a preset value; the expression for the thermal conductivity of black phosphorene is:

[0045]

[0046] In formula (14), represents the component of the thermal conductivity tensor, describing the contribution of the temperature gradient along the α direction to the heat flow along the β direction; Ω is the volume of the unit cell; k B is the Boltzmann constant; T is the temperature; is the group velocity of the phonon mode λ in the α direction; is the driving function related to the λ and β directions of the phonon modes; f0 is the equilibrium distribution function of the phonon;

[0047] Step S5.5: Based on formula (14), the regulation of the anisotropy of thermal conductivity of black phosphorene by the external electric fields of different intensities is obtained, and the corresponding value of thermal conductivity of black phosphorene is calculated.

[0048] Preferably, after step S5, the method further includes: step S6: automatically extracting phonon parameter information; and step S7: based on the phonon parameter information, drawing a conclusion on the mechanism of the effect of an external electric field on the anisotropy of thermal conductivity of two-dimensional black phosphorene.

[0049] Preferably, the step S6 specifically includes: step S6.1: automatically extracting phonon lifetime information and frequency information; step S6.2: automatically extracting phonon group velocity information and frequency information; step S6.3: automatically extracting Greeneson parameter information and frequency information.

[0050] Preferably, the step S7 specifically includes: step S7.1: based on the phonon lifetime information and the frequency information, plotting the curve of the phonon lifetime of black phosphorus ene changing with frequency under different electric field strengths, and drawing the conclusion that applying an external electric field of a certain intensity can prolong the phonon lifetime; step S7.2: based on the phonon group velocity information and the frequency information, plotting the curve of the group velocity of black phosphorus ene changing with frequency under different electric field strengths, and drawing the conclusion that applying an external electric field of a certain intensity can reduce the phonon group velocity; step S7.3: based on the Greenesen parameter information and the frequency information, plotting the curve of the Greenesen parameter of black phosphorus ene changing with frequency under different electric field strengths, and drawing the conclusion that applying an external electric field of a certain intensity can reduce the Greenesen parameter.

[0051] Beneficial effects of the present invention:

[0052] (1) The method provided by the present invention for calculating the thermal conductivity of black phosphorene based on external electric field regulation realizes the reversible and active control of thermal energy transfer by regulating the anisotropic thermal conductivity of black phosphorene by electric field. Specifically, black phosphorene exhibits significant changes in thermal conductivity under the action of the electric field, and this change can be adjusted according to the strength and direction of the electric field, thereby achieving precise control of thermal energy transfer. This method not only overcomes the problem of irreversible thermal conductivity in the prior art, but also provides the ability to actively control the temperature of the device under different heat load conditions.

[0053] (2) The method provided by the present invention for calculating the thermal conductivity of black phosphorene based on the regulation of an external electric field can better understand and optimize the thermal conductivity characteristics of the material by studying the behavior of phonons in black phosphorene and its influence on heat conduction. This not only helps to improve the performance of nanoelectronic devices and thermoelectric materials, but also provides theoretical guidance for the design of new materials with customized heat transfer properties. It helps to understand the microscopic mechanism of phonons on heat transfer, reveals the potential for regulating thermal conductivity under an external electric field, and also demonstrates the efficiency and reversibility of the method of the present invention for calculating the thermal conductivity of black phosphorene based on the regulation of an external electric field, opening up a new path for future technological applications.

[0054] Additional aspects and advantages of the invention will become apparent from the following description, or may be learned by practice of the invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 A schematic flow chart of a method for calculating thermal conductivity of black phosphorene based on applied electric field regulation according to an embodiment of the present invention is shown;

[0056] Figure 2 A schematic diagram of the atomic structure of black phosphorene according to an embodiment of the present invention is shown;

[0057] Figure 3A schematic diagram showing phonon spectrum curves of black phosphorene at different electric field intensities according to an embodiment of the present invention is shown;

[0058] Figure 4 A schematic diagram showing a thermal conductivity variation curve of black phosphorene under different temperature conditions according to an embodiment of the present invention;

[0059] Figure 5 A schematic diagram showing a thermal conductivity variation curve of black phosphorene under different electric field intensities according to an embodiment of the present invention;

[0060] Figure 6 A schematic diagram showing the dependence of the phonon lifetime of black phosphorene in two directions in a plane of an embodiment of the present invention on temperature;

[0061] Figure 7 A schematic diagram showing the variation of the group velocity of black phosphorene with frequency under different electric field strengths according to an embodiment of the present invention;

[0062] Figure 8 A schematic diagram showing the variation of Grün-Eisen parameters of black phosphorene with frequency under different electric field intensities according to an embodiment of the present invention. DETAILED DESCRIPTION

[0063] In order to more clearly understand the above-mentioned objectives, features and advantages of the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be noted that the embodiments of the present application and the features in the embodiments can be combined with each other without conflict.

[0064] In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention may also be implemented in other ways different from those described herein. Therefore, the protection scope of the present invention is not limited to the specific embodiments disclosed below.

[0065] Figure 1 FIG. 1 is a schematic flow chart of a method for calculating the thermal conductivity of black phosphorene based on external electric field regulation according to an embodiment of the present invention. Figure 1 As shown, the method for calculating the thermal conductivity of black phosphorene based on external electric field regulation includes:

[0066] Step S1: constructing a computational model of two-dimensional black phosphorene;

[0067] Step S2: adjusting various parameters of the calculation model of the two-dimensional black phosphorene to optimize the structure of the calculation model of the two-dimensional black phosphorene, and calculating the values ​​of various parameters of the optimized calculation model of the two-dimensional black phosphorene;

[0068] Step S3: applying a plurality of external electric fields of different strengths to the optimized two-dimensional black phosphorene calculation model, and returning to step S2 to obtain the two-dimensional black phosphorene calculation model under different working conditions;

[0069] Step S4: obtaining the harmonic force constant and the anharmonic force constant of the two-dimensional black phosphorene by compressive sensing lattice dynamics;

[0070] Step S5: using an iterative method to solve the phonon Boltzmann transport equation, obtain the regulation of the anisotropy of the thermal conductivity of black phosphorene by external electric fields of several different intensities, and calculate the corresponding value of the thermal conductivity of black phosphorene.

[0071] In this embodiment, the heat transfer characteristics are precisely controlled by regulating the direction of the electric field. The external electric field changes the charge distribution and polarization, affects the interaction between atoms, regulates the phonon behavior and anharmonicity, and enhances the shielding effect of the charge, thereby significantly affecting the thermal conductivity. The method of regulating thermal conductivity by electric field provided by the present invention is efficient and reversible, does not require changing the atomic structure, and has the advantages of being environmentally friendly and responsive.

[0072] In this embodiment, the method provided by the present invention for calculating the thermal conductivity of black phosphorene based on the regulation of an external electric field realizes the reversible and active control of thermal energy transfer by regulating the anisotropic property of thermal conductivity of black phosphorene by electric field. Specifically, black phosphorene exhibits a significant change in thermal conductivity under the action of the electric field, and this change can be adjusted according to the strength and direction of the electric field, thereby achieving precise control of thermal energy transfer. This method not only overcomes the problem of irreversible thermal conductivity in the prior art, but also provides the ability to actively regulate the temperature of the device under different heat load conditions.

[0073] Furthermore, the method provided by the present invention for calculating the thermal conductivity of black phosphorene based on the regulation of an applied electric field can better understand and optimize the thermal conductivity characteristics of the material by studying the behavior of phonons in black phosphorene and its effect on heat conduction. This not only helps to improve the performance of nanoelectronic devices and thermoelectric materials, but also provides theoretical guidance for the design of new materials with customized heat transfer properties. It helps to understand the microscopic mechanism of phonons on heat transfer, reveals the potential for regulating thermal conductivity under an external electric field, and also demonstrates the efficiency and reversibility of the method of the present invention for calculating the thermal conductivity of black phosphorene based on the regulation of an applied electric field, opening up a new path for future technological applications.

[0074] In one embodiment of the present invention, step S1 specifically includes: constructing a computational model of two-dimensional black phosphorene by using crystal structure visualization software VESTA. In a specific embodiment, using crystal structure visualization software VESTA, a Cmca spatial structure is selected to construct a computational model of two-dimensional black phosphorene; wherein the expansion cell of the Cmca spatial structure in the x, y, and z directions is 6×6×1 times.

[0075] In one embodiment of the present invention, the various parameters of the computational model of the two-dimensional black phosphorene include: k-space grid, energy convergence error, force convergence error, and electron convergence error; and the values ​​of the various parameters of the optimized computational model of the two-dimensional black phosphorene are calculated by Quantum Espresso software.

[0076] In one embodiment of the present invention, the k-space grid is set to 15×11×1, and the energy convergence error is set to 1×10 -8 Ry, the force convergence error is set to 1×10 -7 Ry, and the electronic convergence error is 1×10 -10 Ry.

[0077] In one embodiment of the present invention, 9 external electric fields of different strengths are applied to the optimized two-dimensional black phosphorene calculation model, and the process returns to step S2 to obtain the two-dimensional black phosphorene calculation model under different working conditions; the electric field strengths of the 9 external electric fields of different strengths are respectively and

[0078] In one embodiment of the present invention, the step S4 specifically includes: step S4.1: obtaining the Taylor expansion of the total energy in terms of atomic displacement power; the expression of the Taylor expansion is:

[0079]

[0080] In formula (1), V is the total energy; V0 is the energy in the reference state; Φ a is the displacement u from atom a a The associated first-order force constant; Φ ab is the displacement u from atoms a and b a and u b The associated second-order force constant; Φ abc is the displacement u from atoms a, b, and c a 、u b and u c The relevant third-order force constant; u a is the displacement of atom a in Cartesian direction i; u b is the displacement of atom b in Cartesian direction i; u c is the displacement of atom c in Cartesian direction i;

[0081] Step S4.2: Given the training data, the expansion coefficient is automatically selected by compressed sensing lattice dynamics, and the value of the expansion coefficient is determined by applying mathematical techniques. Then, formula (1) can be written as a force-displacement relationship expression; the force-displacement relationship expression is:

[0082] F a =-Φ a -Φ ab u b -Φ abc u b u c / 2-…(2)

[0083] In formula (2), F a is the force acting on atom a; the F a obtained by calculation from first principles using any common density functional theory code, and the F a A set of L atomic configurations used in a supercell to achieve a linear relationship for the unknown force constant tensor Φ; the expression of the linear relationship is:

[0084]

[0085] In formula (3), F is the force vector; Φ is the unknown force constant vector, Φ includes all high-order force constants; is the sensing matrix, The elements of are the products of the atomic displacements, The elements of correspond to the different terms in the force expansion, The specific expression is:

[0086]

[0087] In formula (4), is the displacement of atom b in the first configuration, where the subscript b denotes a specific atom and the superscript 1 denotes the first configuration; is the product of the displacements of atoms b and c in the first configuration; Each row in corresponds to a force component, The expression for the total number of rows M is:

[0088] M=3LN at (5)

[0089] In formula (5), L is the number of atoms in the supercell; N at is the number of atoms in the supercell;

[0090] Step S4.3: Minimizing coefficients based on lattice dynamics via compressive sensing The norm of ||Φ||1≡∑ I |Φ I |To solve formula (3); where, The norm of The approximate value of The norm represents the number of non-zero elements in a vector, ΦI | is the absolute value of the force constant tensor element, and I is the index of the force constant tensor; solving formula (3) yields:

[0091]

[0092] In formula (6), Φ CS represents the force constant tensor optimized by the compressed sensing method; argmin Φ ||Φ||1 means finding the force constant value that minimizes the function ||Φ||1, where ||Φ||1 represents the value of Φ. norm; represents the sum of squares of error terms, μ is a weight parameter; represents the sum of the absolute values ​​of all force constant tensor elements; represents the sum of squares of all error terms.

[0093] In this embodiment, formula (3) is equivalent to the general expression of formula (2). a It is a component of the vector F in formula (3). determines the phonon dispersion in the harmonic approximation, is the third-order anharmonic force constant tensor.

[0094] In this embodiment, The norm of approximation of , and yields a computationally tractable convex optimization problem.

[0095] In one embodiment of the present invention, the step S5 specifically includes: step S5.1: when the only scattering source is a two-phonon and a three-phonon process, a linearized Boltzmann transport equation is obtained; the expression of the linearized Boltzmann transport equation is:

[0096]

[0097] In formula (7), F λ represents the deviation of the corrected phonon distribution function; is the relaxation time of the λ mode, Obtained from perturbation theory; ν λ represents the group velocity of the phonon mode λ; Δ λ is a measure of how much a particular phonon mode deviates from the relaxation time approximation prediction, Δ λ Take speed as the dimension;

[0098] Among them, the Δ λ The calculation expression is as follows:

[0099]

[0100] Among them, the The calculation expression is as follows:

[0101]

[0102] In formula (8) and formula (9), N is the total number of discrete grid points in the Brillouin zone; ξ λλ” is the frequency ratio; F λ” is the driving function of the phonon mode λ"; Γ λλ' is the two-phonon scattering rate; and are the three-phonon scattering rates in the absorption process and the emission process, respectively; is the sum of the three-phonon scattering rates in the absorption process; is the sum of the three-phonon scattering rates during the emission process; is the sum of the two-phonon scattering rates during the emission process;

[0103] Step S5.2: Introduce the Brillouin zone discretization into a Γ-centered regular grid of q points N = N1 × N2 × N3, and set ξ λλ' =ω λ' / ω λ ; Among them, ω λ' is the angular frequency of the phonon mode λ', ω λ is the angular frequency of the phonon mode λ; N1, N2 and N3 are the numbers of discrete points along three different directions of the lattice respectively;

[0104] Step S5.3: Normalized characteristic function e based on three phonons p,q and the anharmonic force constant tensor The three-phonon scattering rates of the absorption process are obtained respectively Three-phonon scattering rate in the emission process and the scattering matrix elements for the three-phonon scattering process Expression of

[0105] Among them, the anharmonic force constant tensor The specific expression is:

[0106]

[0107] In formula (10), E is the total energy of the system; is the displacement component of the i-th atom in the α direction; is the displacement component of the jth atom in the β direction; is the displacement component of the kth atom in the γ direction;

[0108] Among them, the The specific expression is:

[0109]

[0110] Among them, the The specific expression is:

[0111]

[0112] In formula (11) and formula (12), is Planck's constant; f0' and f0" represent the equilibrium distribution functions of phonon modes λ and λ", respectively; ω λ ,ω λ' and ω λ" denote the angular frequencies of the phonon modes λ, λ' and λ", respectively; and are the three-phonon scattering matrix elements in the absorption and emission processes respectively; δ is the Dirac function;

[0113] Among them, the The specific expression is:

[0114]

[0115] In formula (13), M i 、M j and M k denote the masses of atoms i, j, and k respectively; represents the α component of the eigenfunction of mode λ at atom i; represents the β component of the eigenfunction of mode p' at atom j, where ±q' depends on the specific absorption or emission process; represents the γ component of the eigenfunction of mode p" at atom k, where -q" represents the inverse wave vector;

[0116] Step S5.4: Using an iterative method to solve formula (7) to obtain an expression for the thermal conductivity of black phosphorene; wherein the convergence criterion used in the iterative method is to make the relative change of the conductivity tensor calculated by formula (13) less than a preset value; the expression for the thermal conductivity of black phosphorene is:

[0117]

[0118] In formula (14), represents the component of the thermal conductivity tensor, describing the contribution of the temperature gradient along the α direction to the heat flow along the β direction; Ω is the volume of the unit cell; k B is the Boltzmann constant; T is the temperature; is the group velocity of the phonon mode λ in the α direction; is the driving function related to the λ and β directions of the phonon modes; f0 is the equilibrium distribution function of the phonon;

[0119] Step S5.5: Based on formula (14), the regulation of the anisotropy of thermal conductivity of black phosphorene by the external electric fields of different intensities is obtained, and the corresponding value of thermal conductivity of black phosphorene is calculated.

[0120] In this embodiment, N1, N2 and N3 are the numbers of discrete points along three different directions of the lattice (usually the three main axis directions of the lattice).

[0121] In this embodiment, an iterative method is used to solve formula (7) to obtain an expression for the thermal conductivity of black phosphorene. In a specific embodiment, formula (7) is iteratively solved starting from the zero-order approximation: The convergence criterion is to make the relative change of the conductivity tensor calculated by formula (13) less than a default value of 10 -5 .

[0122] In one embodiment of the present invention, after step S5, the method further includes: step S6: automatically extracting phonon parameter information; and step S7: based on the phonon parameter information, drawing a conclusion on the mechanism of the effect of an applied electric field on the anisotropy of thermal conductivity of two-dimensional black phosphorene.

[0123] In one embodiment of the present invention, the step S6 specifically includes: step S6.1: automatically extracting phonon lifetime information and frequency information; step S6.2: automatically extracting phonon group velocity information and frequency information; step S6.3: automatically extracting Grüneisen parameter information and frequency information.

[0124] In this embodiment, the steps of automatically extracting phonon parameter information are specifically described.

[0125] In one embodiment of the present invention, the step S7 specifically includes: step S7.1: based on the phonon lifetime information and the frequency information, plotting the phonon lifetime of black phosphorus as a function of frequency under different electric field strengths, and drawing the conclusion that applying an external electric field of a certain intensity can prolong the phonon lifetime; step S7.2: based on the phonon group velocity information and the frequency information, plotting the group velocity of black phosphorus as a function of frequency under different electric field strengths, and drawing the conclusion that applying an external electric field of a certain intensity can reduce the phonon group velocity; step S7.3: based on the Greenesen parameter information and the frequency information, plotting the Greenesen parameter of black phosphorus as a function of frequency under different electric field strengths, and drawing the conclusion that applying an external electric field of a certain intensity can reduce the Greenesen parameter.

[0126] The technical solution of the present invention will be demonstrated with a specific embodiment below.

[0127] (1) Construct a computational model of two-dimensional black phosphorene.

[0128] The crystal structure visualization software VESTA was used to select the Cmca space group structure to construct a two-dimensional black phosphorene computational model; the structure was expanded to 6×6×1 in the x, y, and z directions. Figure 2 Figure 2 shows the atomic model of black phosphorene. (a) is the top view and (b) is the side view, showing the typical arrangement of wrinkles. The applied electric field is along the zigzag direction, armchair direction and vertical out-of-plane direction, as shown in Figure 2. Figure 2 As shown by the arrow in (b), its influence on the thermal conductivity of the material is studied.

[0129] (2) Optimize the structure of the computational model of two-dimensional black phosphorene.

[0130] The parameters of the two-dimensional black phosphorene calculation model include: k-space grid, energy convergence error, force convergence error, and electron convergence error; the k-space grid is set to 15×11×1, the energy convergence error is set to 1×10 -8 Ry, the force convergence error is set to 1×10 -7 Ry, and the electronic convergence error is 1×10 -10 Ry. The calculation software is Quantum Espresso software package.

[0131] (3) Applying several external electric fields of different strengths to the optimized two-dimensional black phosphorene calculation model, and returning to step (2) to obtain the two-dimensional black phosphorene calculation model under different working conditions.

[0132] Specifically, nine external electric fields of different strengths are applied to the optimized two-dimensional black phosphorene calculation model. The electric field strengths of these nine external electric fields of different strengths are respectively and The external electric field Figure 2 As shown by the arrows in , the pressure is applied in the Zigzag direction, the Armchair direction and the direction perpendicular to the plane to study its influence on the thermal properties of the material. Step (2) is repeated to obtain the calculation model of two-dimensional black phosphorene under different working conditions.

[0133] (4) The harmonic and anharmonic force constants of two-dimensional black phosphorene were obtained by compressive sensing of lattice dynamics.

[0134] The step (4) specifically comprises: step S4.1: obtaining the Taylor expansion of the total energy in terms of atomic displacement power; the expression of the Taylor expansion is:

[0135]

[0136] In formula (1), V is the total energy; V0 is the energy in the reference state; Φ a is the displacement u from atom a a The associated first-order force constant; Φ ab is the displacement u from atoms a and b a and u b The associated second-order force constant; Φ abc is the displacement u from atoms a, b, and c a 、u b and u c The relevant third-order force constant; u a is the displacement of atom a in Cartesian direction i; u b is the displacement of atom b in Cartesian direction i; u c is the displacement of atom c in Cartesian direction i;

[0137] Step S4.2: Given the training data, the expansion coefficient is automatically selected by compressed sensing lattice dynamics, and the value of the expansion coefficient is determined by applying mathematical techniques. Then, formula (1) can be written as a force-displacement relationship expression; the force-displacement relationship expression is:

[0138] F a =-Φ a -Φ ab u b -Φ abc u b u c / 2-…(2)

[0139] In formula (2), F a is the force acting on atom a; the F a obtained by calculation from first principles using any common density functional theory code, and the F a A set of L atomic configurations used in a supercell to achieve a linear relationship for the unknown force constant tensor Φ; the expression of the linear relationship is:

[0140]

[0141] In formula (3), F is the force vector; Φ is the unknown force constant vector, Φ includes all high-order force constants; is the sensing matrix, The elements of are the products of the atomic displacements, The elements of correspond to the different terms in the force expansion, The specific expression is:

[0142]

[0143] In formula (4), is the displacement of atom b in the first configuration, where the subscript b denotes a specific atom and the superscript 1 denotes the first configuration; is the product of the displacements of atoms b and c in the first configuration; Each row in corresponds to a force component, The expression for the total number of rows M is:

[0144] M=3LN at (5)

[0145] In formula (5), L is the number of atoms in the supercell; N at is the number of atoms in the supercell;

[0146] Step S4.3: Minimizing coefficients based on lattice dynamics via compressive sensing The norm of ||Φ||1≡Σ I |Φ I |To solve formula (3); where, The norm of The approximate value of The norm represents the number of non-zero elements in a vector, Φ I | is the absolute value of the force constant tensor element, and I is the index of the force constant tensor; solving formula (3) yields:

[0147]

[0148] In formula (6), Φ CS represents the force constant tensor optimized by the compressed sensing method; argmin Φ ||Φ||1 means finding the force constant value that minimizes the function ||Φ||1, where ||Φ||1 represents the value of Φ. norm; represents the sum of squares of error terms, μ is a weight parameter; represents the sum of the absolute values ​​of all force constant tensor elements; represents the sum of squares of all error terms.

[0149] After the calculation is completed, the simple harmonic force constant of black phosphorene is used to draw a schematic diagram of the phonon spectrum curve of black phosphorene under different electric field strengths. Figure 3 As shown in Figure 2, when the electric field is first applied in the zigzag direction, the phonon dispersion in the high symmetry path (Γ–X–S–Y–Γ) of the first Brillouin zone of black phosphorene is and Figure 3The longitudinal acoustic (LA) and transverse acoustic (TA) phonon branches corresponding to in-plane vibrations, and the ZA branch corresponding to out-of-plane vibrations are depicted in Figure 2. Mainly due to its intrinsic wrinkled geometry, the three acoustic branches of black phosphorene have linear dispersion relations with the wave vector close to the Γ point.

[0150] It is worth noting that the phonon band gap of black phosphorene is 2.6THz. In terms of the phonon group velocity near the Γ point, the velocity along the Γ-X (Zigzag) direction is higher than the velocity along the Γ-Y ​​(Armchair) direction, highlighting the anisotropic nature of phonon transport in black phosphorene. When the phonon dispersion is examined under different electric field strengths, it can be observed that the dispersion of acoustic phonons changes very little, indicating that the group velocity changes very little with the change of electric field strength.

[0151] (5) The phonon Boltzmann transport equation is solved by an iterative method to obtain the regulation of the anisotropy of thermal conductivity of black phosphorene by the external electric fields of different intensities, and the corresponding value of thermal conductivity of black phosphorene is calculated.

[0152] The step (5) specifically comprises:

[0153] Step S5.1: When the only scattering sources are two-phonon and three-phonon processes, a linearized Boltzmann transport equation is obtained; the expression of the linearized Boltzmann transport equation is:

[0154]

[0155] In formula (7), F λ represents the deviation of the corrected phonon distribution function; is the relaxation time of the λ mode, Obtained from perturbation theory; ν λ represents the group velocity of the phonon mode λ; Δ λ is a measure of how much a particular phonon mode deviates from the relaxation time approximation prediction, Δ λ Take speed as the dimension;

[0156] Among them, the Δ λ The calculation expression is as follows:

[0157]

[0158] Among them, the The calculation expression is as follows:

[0159]

[0160] In formula (8) and formula (9), N is the total number of discrete grid points in the Brillouin zone; ξ λλ” is the frequency ratio; F λ”is the driving function of the phonon mode λ"; Γ λλ' is the two-phonon scattering rate; and are the three-phonon scattering rates in the absorption process and the emission process, respectively; is the sum of the three-phonon scattering rates in the absorption process; is the sum of the three-phonon scattering rates during the emission process; is the sum of the two-phonon scattering rates during the emission process;

[0161] Step S5.2: Introduce the Brillouin zone discretization into a Γ-centered regular grid of q points N = N1 × N2 × N3, and set ξ λλ' =ω λ' / ω λ ; Among them, ω λ' is the angular frequency of the phonon mode λ', ω λ is the angular frequency of the phonon mode λ; N1, N2 and N3 are the numbers of discrete points along three different directions of the lattice respectively;

[0162] Step S5.3: Normalized characteristic function e based on three phonons p,q and the anharmonic force constant tensor The three-phonon scattering rates of the absorption process are obtained respectively Three-phonon scattering rate in the emission process and the scattering matrix elements for the three-phonon scattering process Expression of

[0163] Among them, the anharmonic force constant tensor The specific expression is:

[0164]

[0165] In formula (10), E is the total energy of the system; is the displacement component of the i-th atom in the α direction; is the displacement component of the jth atom in the β direction; is the displacement component of the kth atom in the γ direction;

[0166] Among them, the The specific expression is:

[0167]

[0168] Among them, the The specific expression is:

[0169]

[0170] In formula (11) and formula (12), is Planck's constant; f0' and f0" represent the equilibrium distribution functions of phonon modes λ and λ", respectively; ω λ ,ω λ' and ω λ" denote the angular frequencies of the phonon modes λ, λ' and λ", respectively; and are the three-phonon scattering matrix elements in the absorption and emission processes respectively; δ is the Dirac function;

[0171] Among them, the The specific expression is:

[0172]

[0173] In formula (13), M i 、M j and M k denote the masses of atoms i, j, and k respectively; represents the α component of the eigenfunction of mode λ at atom i; represents the β component of the eigenfunction of mode p' at atom j, where ±q' depends on the specific absorption or emission process; represents the γ component of the eigenfunction of mode p" at atom k, where -q" represents the inverse wave vector;

[0174] Step S5.4: Use the iterative method to solve formula (7) to obtain the expression of thermal conductivity of black phosphorene; Formula (7) is iteratively solved starting from the zero-order approximation The convergence criterion is to make the relative change of the conductivity tensor calculated by formula (13) less than a default value of 10 -5 .

[0175] The expression of thermal conductivity of black phosphorene is:

[0176]

[0177] In formula (14), represents the component of the thermal conductivity tensor, describing the contribution of the temperature gradient along the α direction to the heat flow along the β direction; Ω is the volume of the unit cell; k B is the Boltzmann constant; T is the temperature; is the group velocity of the phonon mode λ in the α direction; is the driving function related to the λ and β directions of the phonon modes; f0 is the equilibrium distribution function of the phonon;

[0178] Step S5.5: Based on formula (14), the regulation of the anisotropy of thermal conductivity of black phosphorene by the external electric fields of different intensities is obtained, and the corresponding value of thermal conductivity of black phosphorene is calculated.

[0179] After the calculation is completed, a schematic diagram of the thermal conductivity change curve of black phosphorene under different temperature conditions is drawn. Figure 4 As shown in Figure 1, at 300K, the thermal conductivity of black phosphorene in the zigzag direction is 20.81W / mK, and the thermal conductivity in the armchair direction is 6.52W / mK. Furthermore, a schematic diagram of the thermal conductivity change curve of black phosphorene under different electric field strengths is drawn. Figure 5 As shown, with the increase of the electric field intensity along the Zigzag direction, the thermal conductivity of black phosphorene gradually decreases, showing that the anisotropy is reduced by up to 28%.

[0180] (6) Use self-written programs to automatically extract phonon parameter information.

[0181] The step (6) specifically comprises the following steps:

[0182] 1) Automatically extract phonon lifetime information and frequency information:

[0183] ① Read the total number of Q points from the 'BTE.qpoints_full' file and the number of irreducible Q points from the 'BTE.qpoints' file.

[0184] ②Initialize the arrays of frequency and life data.

[0185] ③ Read the life and frequency data from the 'BTE.w_anharmonic' file and process it.

[0186] ④ Write the processed data to two output files 'lifetime_freq.dat', one containing all the frequency and lifetime data, ⑤ Calculate and write the average lifetime of each mode to the 'lifetime_ilmd_ave.dat' file.

[0187] After extracting the phonon lifetime, a schematic diagram of the dependence of the phonon lifetime of black phosphorene in two directions in the plane on temperature is drawn. Figure 6 As shown in Figure 2, we elucidate the temperature dependence of the phonon lifetime of black phosphorene in two in-plane directions (zigzag and armchair). Compared with the case without electric field, the phonon lifetime is significantly prolonged under the electric field, especially at a frequency of about 3 THz.

[0188] 2) Automatically extract phonon group velocity information and frequency information:

[0189] ① Read the total number of Q points from the 'BTE.qpoints_full' file and the number of irreducible Q points from the 'BTE.qpoints' file.

[0190] ②Initialize the frequency and speed arrays.

[0191] ③ Read frequency data from the 'BTE.omega' file, read speed data from the 'BTE.v' file, and process them.

[0192] ④Write the frequency and speed data into the 'v_freq.dat' file

[0193] After extracting the phonon group velocity, a schematic diagram of the change of the group velocity of black phosphorene with frequency under different electric field strengths was drawn. Figure 7 As shown in Figure 2, the phonon group velocity near the Γ point is higher along the Γ-X (sawtooth) direction than along the Γ-Y ​​(armchair) direction, highlighting the anisotropic nature of phonon transport in black phosphorene. When the electric field (indicated by the blue dots) is applied, the group velocity decreases slightly. In addition, the group velocity changes in the zigzag and armchair directions are 3.5% and 2.2%, respectively.

[0194] Two key conclusions are drawn from this: (1) The differences in the rate of change of group velocities in these directions lead to changes in the anisotropic thermal conductivity under the action of the electric field. (2) However, these changes in group velocity are relatively small and do not significantly affect the overall change in lattice thermal conductivity.

[0195] 3) Automatically extract GreenEisen parameter information and frequency information:

[0196] ① Read the total number of Q points from the 'BTE.qpoints_full' file and the number of irreducible Q points from the 'BTE.qpoints' file.

[0197] ②Initialize the frequency and Grueneisen parameter arrays.

[0198] ③Read frequency data from the 'BTE.omega' file.

[0199] ④ Read the Grueneisen parameter data from the 'BTE.gruneisen' file.

[0200] ⑤Write the frequency and Grueneisen parameter data into the 'gruneisen_freq.dat' file.

[0201] After extracting the Grüneisen parameters, a schematic diagram of the Grüneisen parameters of black phosphorene with frequency under different electric field strengths was drawn. Figure 8 As shown, for B 1uThe Grüneisen parameter analysis of the branch shows that The Grüneisen parameter decreases significantly when an electric field is applied. and The average Grüneisen parameters under θ are -1.42 and -5.8, respectively. This observation emphasizes the principle that the larger the absolute value of the Grüneisen parameter, the stronger the anharmonicity. Under electric field, B 1u The absolute value of the branched Grüneisen parameter increases, resulting in a decrease in thermal conductivity. Therefore, this change in phonon anharmonicity is the basis for the change in phonon lifetime and thermal conductivity of black phosphorene under the action of an electric field.

[0202] In summary, the present invention achieves reversible and active control of thermal energy transfer by regulating the anisotropy of thermal conductivity of black phosphorene through electric field, and reveals the mechanism of electric field on the anisotropy of thermal conductivity of two-dimensional black phosphorene from the phonon scale. Some useful conclusions are obtained:

[0203] 1. At room temperature, the thermal conductivity of black phosphorene in the zigzag direction exceeds that in the armchair direction, showing an anisotropy ratio of 3.13.

[0204] 2. When an electric field in the zigzag direction is applied, the anisotropy of the thermal conductivity of black phosphorene oscillates slightly, eventually leading to a significant decrease of up to 28%.

[0205] 3. Phonon correlation analysis shows that the anisotropy of the intrinsic thermal conductivity of black phosphorene is mainly due to the difference in phonon group velocities in the zigzag and armchair directions.

[0206] 4. When a moderate external electric field is applied, the phonon lifetime becomes a key regulating factor, significantly reducing the anisotropy of thermal conductivity of black phosphorene. Under the electric field strength of , the phonon lifetime is significantly shortened, which significantly reduces the thermal conduction anisotropy.

[0207] 5. This invention provides valuable insights and guidance for designing nanomaterials with customized thermal conductivity properties, which is critical to optimizing performance in specific energy-related applications.

[0208] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for calculating the thermal conductivity of black phosphorene based on external electric field regulation, characterized in that: include: Step S1: constructing a computational model of two-dimensional black phosphorene; Step S2: adjusting various parameters of the calculation model of the two-dimensional black phosphorene to optimize the structure of the calculation model of the two-dimensional black phosphorene, and calculating the values ​​of various parameters of the optimized calculation model of the two-dimensional black phosphorene; Step S3: applying a plurality of external electric fields of different strengths to the optimized two-dimensional black phosphorene calculation model, and returning to step S2 to obtain the two-dimensional black phosphorene calculation model under different working conditions; Step S4: obtaining the harmonic force constant and the anharmonic force constant of the two-dimensional black phosphorene by compressive sensing lattice dynamics; Step S5: solving the phonon Boltzmann transport equation by an iterative method, obtaining the regulation of the anisotropy of thermal conductivity of black phosphorene by the external electric fields of different intensities, and calculating the corresponding value of thermal conductivity of black phosphorene; The step S4 specifically includes: Step S4.1: Obtain the Taylor expansion of the total energy in terms of atomic displacement power; the expression of the Taylor expansion is: In formula (1), V is the total energy; V0 is the energy in the reference state; Φ a is the displacement u from atom a a The associated first-order force constant; Φ ab is the displacement u from atoms a and b a and u b The associated second-order force constant; Φ abc is the displacement u from atoms a, b, and c a 、u b and u c The relevant third-order force constant; u a is the displacement of atom a in Cartesian direction i; u b is the displacement of atom b in Cartesian direction i; u c is the displacement of atom c in Cartesian direction i; Step S4.2: Given the training data, the expansion coefficient is automatically selected by compressed sensing lattice dynamics, and the value of the expansion coefficient is determined by applying mathematical techniques. Then, formula (1) can be written as a force-displacement relationship expression; the force-displacement relationship expression is: F a =-Φ a -Φ ab in b -Φ abc in b in c / 2-…(2) In formula (2), F a is the force acting on atom a; a obtained by calculation from first principles using any general density functional theory code, and the F a A set of L atomic configurations used in a supercell to achieve a linear relationship for the unknown force constant tensor Φ; the expression of the linear relationship is: In formula (3), F is the force vector; Φ is the unknown force constant vector, Φ includes all high-order force constants; is the sensing matrix, The elements of are the products of the atomic displacements, The elements of correspond to the different terms in the force expansion, The specific expression is: In formula (4), is the displacement of atom b in the first configuration, where the subscript b denotes a specific atom and the superscript 1 denotes the first configuration; is the product of the displacements of atoms b and c in the first configuration; Each row in corresponds to a force component, The expression for the total number of rows M is: <h2 style=";text-align:left;direction:ltr">M=3LN<h2 style=";text-align:left;direction:ltr"> at <h2 style=";text-align:left;direction:ltr"> (5) In formula (5), L is the number of atoms in the supercell; N at is the number of atoms in the supercell; Step S4.3: Minimizing coefficients based on lattice dynamics via compressive sensing The norm of ||Φ||1≡∑ I |Φ I |To solve formula (3); where, The norm of The approximate value of The norm represents the number of non-zero elements in a vector, |Φ I | is the absolute value of the force constant tensor element, and I is the index of the force constant tensor; solving formula (3) yields: In formula (6), Φ CS represents the force constant tensor optimized by the compressed sensing method; argmin Φ ||Φ||1 means finding the force constant value that minimizes the function ||Φ||1, where ||Φ||1 represents the value of Φ. norm; represents the sum of squares of error terms, μ is a weight parameter; represents the sum of the absolute values ​​of all force constant tensor elements; represents the sum of squares of all error terms.

2. The method for calculating the thermal conductivity of black phosphorene based on external electric field regulation according to claim 1 is characterized in that: The step S1 specifically includes: A computational model of two-dimensional black phosphorene was constructed by using the crystal structure visualization software VESTA.

3. The method for calculating the thermal conductivity of black phosphorene based on external electric field regulation according to claim 1, characterized in that: The parameters of the two-dimensional black phosphorene calculation model include: k-space grid, energy convergence error, force convergence error, and electron convergence error; and The values ​​of various parameters of the optimized two-dimensional black phosphorene calculation model are calculated using Quantum Espresso software.

4. The method for calculating the thermal conductivity of black phosphorene based on external electric field regulation according to claim 1, characterized in that: The step S5 specifically includes: Step S5.1: When the only scattering sources are two-phonon and three-phonon processes, a linearized Boltzmann transport equation is obtained; the expression of the linearized Boltzmann transport equation is: In formula (7), F λ represents the deviation of the corrected phonon distribution function; is the relaxation time of the λ mode, Obtained from perturbation theory; ν λ represents the group velocity of the phonon mode λ; Δ λ is a measure of how much a particular phonon mode deviates from the relaxation time approximation prediction, Δ λ Take speed as the dimension; Among them, the Δ λ The calculation expression is as follows: Among them, the The calculation expression is as follows: In formula (8) and formula (9), N is the total number of discrete grid points in the Brillouin zone; ξ λλ” is the frequency ratio of λ and λ”; F λ” is the driving function of the phonon mode λ"; Γ λλ' is the two-phonon scattering rate; and are the three-phonon scattering rates in the absorption process and the emission process, respectively; is the sum of the three-phonon scattering rates in the absorption process; is the sum of the three-phonon scattering rates during the emission process; is the sum of the two-phonon scattering rates during the emission process; Step S5.2: Introduce the Brillouin zone discretization into a Γ-centered regular grid of q points N = N1 × N2 × N3, and set ξ λλ' =ω λ' / ω λ ; Among them, ω λ' is the angular frequency of the phonon mode λ', ω λ is the angular frequency of the phonon mode λ; N1, N2 and N3 are the numbers of discrete points along three different directions of the lattice respectively; Step S5.3: Normalized characteristic function e based on three phonons p,q and the anharmonic force constant tensor The three-phonon scattering rates of the absorption process are obtained respectively Three-phonon scattering rate in the emission process and the scattering matrix elements for the three-phonon scattering process Expression of Among them, the anharmonic force constant tensor The specific expression is: In formula (10), E is the total energy of the system; r i α is the displacement component of the i-th atom in the α direction; is the displacement component of the jth atom in the β direction; is the displacement component of the kth atom in the γ direction; Among them, the The specific expression is: Among them, the The specific expression is: In formula (11) and formula (12), is Planck's constant; f0' and f0" represent the equilibrium distribution functions of phonon modes λ and λ", respectively; ω λ ,ω λ' and ω λ" denote the angular frequencies of the phonon modes λ, λ' and λ", respectively; and are the three-phonon scattering matrix elements in the absorption and emission processes respectively; δ is the Dirac function; Among them, the The specific expression is: In formula (13), M i 、M j and M k denote the masses of atoms i, j, and k respectively; represents the α component of the eigenfunction of mode λ at atom i; represents the β component of the eigenfunction of mode p' at atom j, where ±q' depends on the specific absorption or emission process; represents the γ component of the eigenfunction of mode p" at atom k, where -q" represents the inverse wave vector; Step S5.4: Using an iterative method to solve formula (7) to obtain an expression for the thermal conductivity of black phosphorene; wherein the convergence criterion used in the iterative method is to make the relative change of the conductivity tensor calculated by formula (13) less than a preset value; the expression for the thermal conductivity of black phosphorene is: In formula (14), represents the component of the thermal conductivity tensor, describing the contribution of the temperature gradient along the α direction to the heat flow along the β direction; Ω is the volume of the unit cell; k B is the Boltzmann constant; T is the temperature; is the group velocity of phonon mode λ in the α direction; is the driving function related to the λ and β directions of the phonon modes; f0 is the equilibrium distribution function of the phonon; Step S5.5: Based on formula (14), the regulation of the anisotropy of thermal conductivity of black phosphorene by the external electric fields of different intensities is obtained, and the corresponding value of thermal conductivity of black phosphorene is calculated.

5. The method for calculating the thermal conductivity of black phosphorene based on external electric field regulation according to any one of claims 1 to 4, characterized in that: After step S5, the method further includes: Step S6: Automatically extract phonon parameter information; Step S7: Based on the phonon parameter information, a conclusion is drawn about the mechanism of the effect of an external electric field on the anisotropy of thermal conductivity of two-dimensional black phosphorene.

6. The method for calculating the thermal conductivity of black phosphorene based on external electric field regulation according to claim 5, characterized in that: The step S6 specifically includes: Step S6.1: Automatically extract phonon lifetime information and frequency information; Step S6.2: Automatically extract phonon group velocity information and frequency information; Step S6.3: Automatically extract Greenisen parameter information and frequency information.

7. The method for calculating the thermal conductivity of black phosphorene based on external electric field regulation according to claim 6, characterized in that: The step S7 specifically includes: Step S7.1: based on the phonon lifetime information and the frequency information, plotting the curve of the phonon lifetime of black phosphorene as a function of frequency under different electric field intensities, and drawing a conclusion that applying an external electric field of a certain intensity can prolong the phonon lifetime; Step S7.2: based on the phonon group velocity information and the frequency information, a curve of the group velocity of black phosphorene changing with frequency under different electric field intensities is plotted, and a conclusion is drawn that applying an external electric field of a certain intensity can reduce the phonon group velocity; Step S7.3: Based on the Grün-Eisen parameter information and the frequency information, a Grün-Eisen parameter variation curve of black phosphorene at different electric field strengths is plotted, and a conclusion is drawn that the Grün-Eisen parameter can be reduced by applying an external electric field of a certain strength.

Images