 Feryal Arabaci-Colak

Selon l'article L223-5 du Code de l'action sociale et des familles, « Le service élabore au moins une fois par an, ou tous les six mois pour les enfants âgés de moins de deux ans, un rapport, établi après une évaluation pluridisciplinaire, sur la situation de tout enfant accueilli ou faisant l'objet d'une mesure éducative. Ce rapport porte sur la santé physique et psychique de l'enfant, son développement, sa scolarité, sa vie sociale et ses relations avec sa famille et les tiers intervenant dans sa vie. »

Feryal Arabaci-Colak 

Écrits des professionnels en protection de l'enfance : éléments d'histoire des enfants confiés à l'Aide sociale à l'enfance L'histoire des enfants pris en charge par l'Aide sociale à l'enfance (ASE) n'est pas seulement celle qui se vit et s'écrit avec les membres de sa famille, mais aussi celle qui se construit avec l'intervention de professionnels de la protection de l'enfance. Les paroles dites par l'éducateur, le psychologue, ou l'assistante familiale de l'enfant vont agir sur sa manière de se percevoir, sa manière de concevoir sa propre histoire personnelle ou familiale, et de considérer sa place dans le monde. Les écrits des professionnels sont essentiels durant son suivi : écrits des professionnels éducatifs, ou jugement écrit des magistrats, ceux-ci vont appuyer des décisions qui donneront des orientations sur le devenir réel de l'enfant. Outre sa dimension procédurière et protocolaire, l'écrit du professionnel de la protection de l'enfance peut revêtir une dimension clinique qui permette une décision prise au plus près de la situation de l'enfant et de sa famille. Nous verrons ainsi l'objet et la finalité des écrits en protection de l'enfance, leur rôle quant à l'implication des professionnels ainsi que leurs apports dans les prises en charge des enfants placés.

L'écrit du professionnel : élément essentiel des procédures propres au système de la protection de l'enfance

Lorsque la situation d'un enfant présente des éléments d'inquiétude, cela doit être retranscrit dans un écrit qui s'appelle « Information préoccupante » (IP). Une IP débouchera soit sur une évaluation en protection de l'enfance, soit sur une séparation physique de l'enfant de son milieu familial. Une mesure de placement provisoire, peut être décidée par un magistrat dans une ordonnance de placement provisoire (OPP), dont la limite sera courte avant qu'une audience mène soit vers un retour de l'enfant dans sa famille, soit vers un maintien de la mesure de placement, et ce dans le cadre d'un jugement en assistance éducative (JAE). À chaque étape de la procédure, les écrits des professionnels de la protection de l'enfance seront nécessaires car ils permettront au magistrat pour enfants d'avoir des éléments sur la situation de l'enfant. Le rapport éducatif, obligatoire selon la Loi, « porte sur la santé physique et psychique de l'enfant, son développement, sa scolarité, sa vie sociale et ses relations avec sa famille et les tiers intervenant dans sa vie 1 ».

Éléments centraux des procédures en protection de l'enfance, les écrits sont aussi des éléments des dossiers des enfants qui y sont pris en charge. Le dossier d'ASE est composé de tout document administratif, médical, socio-éducatif ou juridique relatif aux jeunes qui font l'objet d'un accompagnement par le service de l'ASE. Traces des éléments de leur vie, traces de leur vécu en tant que mineurs confiés à l'ASE, il leur est possible de consulter leur dossier administratif, qui sera archivé. Notre écrit, qui consiste à donner des éléments de compréhension sur la dimension éducative et psychique de la situation de l'enfant, peut également être considéré comme un moyen de communication avec lui une fois devenu adulte s'il décide de le consulter.

Depuis 19782 , les lois ouvrent à chacun le droit d'accès aux documents nous concernant. Plusieurs lois, depuis, ont précisé les droits des personnes accueillies dans le secteur social, médical et médico-social. La loi du 2 janvier 2002 confirme la responsabilité des institutions et des services dans l'accès, la tenue et la conservation du dossier unique des personnes. À l'ASE, toute personne suivie ou ayant été suivie peut consulter ou avoir copie de son dossier. Si la personne est mineure, elle a besoin d'une autorisation du représentant légal. Si elle est majeure, elle peut le faire seule ou être accompagnée par la personne de son choix.

Pour Alice Bénavente et Henri Mialocq (2003) : « Le dossier est porteur d'une forte charge symbolique. Il est la trace de l'existence passée, l'album qui fixe les différentes étapes de l'enfance et de l'adolescence de la personne. Il confirme au consultant le sentiment d'avoir existé et celui d'avoir été reconnu par la société. En "écrivant sur", le professionnel marque son "intérêt pour" et conditionne "la capacité d'existence de" celui dont il parle ». Ces auteurs insistent sur le soulagement repéré chez les consultants de ces dossiers à la lecture des écrits les concernant, et affirment que, souvent, ces découvertes n'ont pas d'effet cataclysmique, mais au contraire apaisant.

Concrètement, dans le cadre de mes observations, étant intervenue comme psychologue durant treize années en service d'ASE, j'ai constaté que ce sont les travailleurs sociaux référents de l'enfant suivi qui rédigent ces rapports ; les écrits des psychologues sont plus rares. Lorsque les psychologues établissent un écrit se joignant à celui du référent éducatif, l'écrit du psychologue peut mettre en avant les éléments de la dimension psychique, affective et relationnelle, alors que le rapport éducatif peut insister sur les dimensions plus factuelles de la vie de l'enfant.

Vignette clinique

La situation que j'ai choisi d'évoquer est celle d'une enfant placée à l'ASE depuis sa naissance. Je relate ici des éléments que j'ai lus dans le dossier de Sarah, puisqu'au moment de son placement je ne travaillais pas encore au service. La mère de Sarah, madame H, ellemême confiée à l'ASE durant sa minorité, se présente aux urgences d'un centre hospitalier pour des douleurs abdominales. Il s'avère qu'elle est sur le point d'accoucher, Sarah vient donc au monde suite à un déni de grossesse de sa mère. Après l'accouchement, sa mère dira plusieurs fois que ce n'est pas possible qu'elle puisse être sa fille. Ayant quitté le service de maternité plusieurs heures sans prévenir le personnel soignant, ces derniers adressent une IP au service de protection de l'enfance. Le juge des enfants décide d'un placement de Sarah à l'ASE. À son retour à l'hôpital, madame aura la possibilité de voir sa fille dans le cadre de visites médiatisées (c'est-à-dire exclusivement en la présence de professionnels de la protection de l'enfance). Seules deux rencontres pourront avoir lieu, madame étant très en retard pour l'une d'elles, et écourtant l'autre. Sarah sera accueillie dans une famille d'accueil du service lorsqu'elle est âgée de deux mois. Auparavant, elle est restée en maternité puis en pédiatrie, le temps que l'ASE puisse trouver un lieu d'accueil. Une pouponnière avait été sollicitée au préalable, démarche n'ayant pu aboutir, Sarah arrive chez madame K, assistante familiale. Elle est alors dans un état chétif, dormant beaucoup, et présentant des difficultés d'interaction avec son environnement, faute de premiers soins privilégiés avec un adulte.

Madme H n'a honoré, par la suite, aucun des rendez-vous proposés par le service pour rencontrer sa fille. Le service ASE a déposé une requête auprès du juge des affaires familiales (JAF), pour la reconnaissance de la situation de délaissement parental3 dans laquelle se trouve alors Sarah. Une fois considérée comme pupille, Sarah peut être adoptable. Madame K souhaite adopter Sarah4 , qu'elle considère comme sa fille depuis qu'elle l'accueille chez elle.

Comme il s'agit d'une étape importante concernant la situation de Sarah, nous organisons avec l'éducatrice une rencontre afin de lui en parler. Sarah a seize mois lorsque je la rencontre pour la première fois ; l'assistante familiale assiste également à l'entretien. Nous expliquons à madame K que, pour l'instant, la première étape est celle de l'audience. La mère est également convoquée. Si elle s'y présente, le JAF peut très bien donner un avis négatif à la requête, ou reporter sa décision à une date ultérieure. L'éducatrice précise, sur un ton insistant, que ce fut déjà le cas pour d'autres situations qu'elle suit. Son insistance vient du fait qu'elle a remarqué que, pour madame K, Sarah est déjà sa fille, et que cette audience n'est qu'une démarche anodine, le « scénario » semble déjà tout tracé, selon elle.

Il se trouve que l'éducatrice avait contacté madame H pour une dernière tentative de rencontre avant l'audience auprès du JAF. Madame H s'était bien présentée à cette rencontre, durant laquelle elle avait affirmé qu'elle ne renoncerait pas à sa fille, et qu'elle ferait en sorte de la garder sous sa filiation. Le JAF ne se prononcera pas sur une déclaration de délaissement, puisqu'au vu des fragilités psychologiques relatées de la situation de madame, une expertise psychiatrique est nécessaire afin de clarifier si la situation de délaissement est volontaire du côté de madame, ou involontaires car liée à des troubles d'ordre psychiatriques.

Sarah a presque 3 ans quand elle voit sa mère pour la première fois, en ma présence et celle de l'éducatrice. Nous avons déjà organisé trois entretiens avec la mère, puis une rencontre avec l'assistante familiale et Sarah. Alors que Sarah nommait madame K « maman », suite aux rencontres avec sa mère, elle la désignera tantôt par « maman », tantôt par « tata » (« tata » étant la nomination que l'éducatrice et moi suggérons pour désigner madame K), et nommera sa mère biologique « maman-Iéna5 ». Même si la procédure de délaissement sera laissée de côté, le placement de Sarah se poursuivra. Maintenue au sein de la famille d'accueil, elle verra sa mère lors de rencontres médiatisées au service, rencontres qui seront d'abord assumées par madame H de manière régulière, puis de moins en moins dans leur fréquence. Sarah est âgée de 5 ans, lorsque sa mère n'honore plus aucun rendezvous à l'ASE, et cesse de voir sa fille.

Madame K, elle-même, a perdu sa mère très jeune, et a été élevée par une tante maternelle. Pour elle, sa mère d'adoption prend la place de sa mère qu'elle a perdue. Là où pour nous « tata » la distinguait de « maman », pour elle, sa « tata » a été sa mère d'adoption. Madame H, quant à elle, expliqua que sa propre mère avait été adoptée. Au moment de la naissance de Sarah, elle avait la crainte qu'on lui rapte sa fille. Son absence de contacts durant deux ans aurait pu générer une telle situation : madame K aurait effectivement pu adopter sa fille. Madame H dit qu'elle avait choisi le prénom « Sarah » car il appartient aux trois religions monothéistes ; madame K nous disait, avant que madame H renoue des contacts avec sa fille que, si elle avait eu une fille, elle l'aurait prénommée « Sarah », puisque ce prénom appartient aux trois religions monothéistes. Nous notons ici comment plusieurs éléments concernant l'histoire de Sarah se prêtent à des télescopages, que les professionnels sont amenés à identifier et à nommer dans les différents écrits la concernant, pour que cela ne soit pas des éléments opaques de son histoire.

L'éducatrice référente a changé de service, elle était pourtant présente depuis le début du placement. Madame K s'est maintenue dans une posture de rivalité à l'égard de la figure de la mère biologique, de telle sorte qu'une fois la mère n'étant plus présente lors des rencontres, elle reprenait la place de la mère auprès de Sarah, qui d'ailleurs l'appelait à nouveau « maman », sans que celle-ci le rectifie. Les écrits professionnels dans cette situation sont une trame, un fil conducteur de l'histoire de Sarah. Ils évoquent la chronologie des faits la concernant, et permettent aux différents professionnels qui interviendront durant son placement de savoir ce qui s'est passé, même si elle ne s'en souvient plus et même si son assistante familiale ne l'évoque pas ; il s'agira, pour la professionnelle qui remplacera l'éducatrice de départ, d'avoir des traces du travail de compréhension et d'analyse de la situation depuis le début du placement.

Analyse du cas clinique

Le lien particulier entre Sarah et son assistante familiale est celui d'une « adoption silencieuse », que Myriam [START_REF] David | Le placement familial. De la pratique à la théorie[END_REF] désigne comme une adoption de fait de l'enfant par sa famille d'accueil, hors de tout cadre juridique. Adoption illusoire, nous insistons sur l'importance d'en repérer les enjeux, afin de donner à l'enfant une lecture assez claire et diversifiée de sa situation, notamment au travers des rencontres avec ce dernier mais aussi via les écrits qui prennent place dans son dossier ASE.

Ainsi, M. [START_REF] David | Le placement familial. De la pratique à la théorie[END_REF] a noté la situation particulière à laquelle est confronté chaque enfant confié en famille d'accueil : l'enfant grandit alors avec le sentiment d'être divisé, car il a le sentiment d'appartenir tantôt à sa famille naturelle et tantôt à sa famille d'accueil. Elle précise que l'enfant lutte compulsivement pour et contre son appartenance à ces deux familles. Elle désigne cela sous le terme de « mal de placement », l'enfant souhaite mais redoute aussi de quitter et de retourner chez lui, tout comme il refuse de se laisser adopter par la famille d'accueil, alors qu'il le désire également. M. David pense qu'il est nécessaire que l'accompagnement des équipes de placements familiaux se centre sur l'enfant (parce qu'il est l'objet du placement), et sur le lien entre l'enfant et ses parents.

Jean Claude Rouchy (1997) évoque, de même, la complexité de la situation des familles d'accueil lorsqu'elles prennent en charge un enfant placé. À la question de la place de l'enfant placé, Rouchy se demande si la famille d'accueil est, pour l'enfant, un groupe d'appartenance secondaire. Il répond que, sous l'angle du placement et de la décision judiciaire, certainement, alors que du côté de l'enfant, la famille d'accueil serait incluse dans son groupe d'appartenance primaire, qui alors est structuré dans des rapports de dualité, avec des éléments de contradiction et de conflictualité liées à la bipolarité représentée par la famille d'accueil et la famille d'origine. L'individuation, pour Rouchy, s'élabore dans le groupe d'appartenance primaire, en rapport avec les groupes d'appartenance secondaire où la maturation s'effectue dans la dialectique de processus primaire et secondaire. L'individuation aurait un effet positif, si le groupe est contenant, et qu'il présente un espace transitionnel, tout en proposant des images identificatoires enrichissantes pour la structuration du Moi. Rouchy précise que des carences dans l'un de ces facteurs auront des répercussions sur la façon dont l'enfant peut se constituer en tant que sujet, s'individualiser et se socialiser progressivement. L'aspect contenant étant celui qui surdétermine les deux autres, car dès lors que celui-ci est carencé, l'espace transitionnel est impacté, et le processus de subjectivation se trouve ainsi en difficulté.

En ce qui concerne les alliances familiales, R. [START_REF] Kaës | L'idéologie, l'idéal, l'idée[END_REF] explique que la réalité psychique qui s'y forme et l'identité familiale qui en résulte, reposent sur un ensemble d'alliances, de pactes et de contrats divers et variables, non seulement dans leurs formes, leurs contenus, et leurs fonctions, mais aussi dans la qualité et la quantité des partenaires qu'ils lient. Il s'agit alors de veiller à ce qu'un enfant placé puisse grandir avec l'appartenance à ses deux familles : famille d'accueil et famille naturelle.

En ce qui concerne le lien familial biologique, nous pouvons nous référer à Maurice [START_REF] Berger | Ces enfants qu'on sacrifie… au nom de la protection de l'enfance[END_REF], qui l'interroge lorsqu'il évoque la prévalence de l'idéologie du lien biologique familial propre au système de la protection de l'enfance, qu'il définit ainsi : « position de principe selon laquelle le maintien du lien physique réel entre l'enfant et ses parents a une valeur absolue et intouchable ». Dans ces cas, le risque serait, pour l'ASE, d'occulter l'importance des liens affectifs entre l'enfant et sa famille d'accueil.

Pour [START_REF] Kaës | L'idéologie, l'idéal, l'idée[END_REF], l'idéologie est une construction de bases requises pour toute action ou pour toute pensée collective. Elle est donc nécessaire pour agir et penser, a fortiori pour agir et penser ensemble. Dans ce qu'il appelle le pacte dénégatif, qui est une forme pathologique du contrat narcissique, se produit un assèchement narcissique de l'investissement de l'objet, lorsque l'institution et le sujet ne parviennent pas à nouer leurs intérêts narcissiques dans un contrat identificatoire porteur d'un processus de subjectivation. Il note que cette forme pathologique est présente dans la position idéologique. Alors que, dans le contrat narcissique, les liens vivants se fabriquent et évoluent, dans le pacte narcissique, les liens sont fixés sans transformation, ils sont polarisés par le narcissisme de mort. Ici, madame K se saisit du désinvestissement narcissique de l'enfant par sa mère, pour l'investir comme sa propre fille. Dérive pathologique des alliances inconscientes, le pacte dénégatif se met en place sur la base du déni, du rejet et du désaveu, ce qui a pour conséquence de créer dans le lien, et en chacun de ses sujets, de l'énigmatique, du non-signifiable, du nontransformable. Le pacte dénégatif exige l'alliance de plusieurs sujets et assigne à chacun la place de gardien et de garant du pacte. Kaës explique que l'idéologie est un pacte dénégatif : on ne peut croire seul ; il faut que chacun y croie (ou fasse semblant d'y croire), pour que chacun continue à y croire (ou à faire semblant d'y croire). Un pacte dénégatif (celui pour la famille d'accueil de croire que Sarah est leur fille) est-il écrasé par un autre (l'idéologie du lien familial portée par l'ASE qui consiste à croire que c'est madame H qui est la mère) ? Les alliances aliénantes, qui peuvent alors apparaître, sont fondées sur le déni ou le désaveu, le rejet ou la forclusion. Elles rendent les sujets qui s'y lient radicalement étrangers à euxmêmes, telles qu'elles se présentent dans le pacte dénégatif. [START_REF] Kaës | Filiation et affiliation. Quelques aspects de la réélaboration du roman familial dans les familles adoptives, les groupes et les institutions[END_REF] définit la filiation comme impliquant le rapport d'au moins trois générations successives reconnues comme telles, et la référence commune à un mythe originant. Les écrits des professionnels, qui inscrivent l'enfant dans une histoire qui précède le placement, peuvent ainsi être un moyen d'évoquer ce mythe, ils concourent alors à lever les dénis et les désaveux, et à mettre en évidence indirectement les enjeux des pactes dénégatifs qui sont parfois à l'oeuvre dans les situations de l'ASE.

Rôle des écrits professionnels en protection de l'enfance : élément central dans les processus d'historicisation de l'enfant placé

Dans les situations de placement, l'enfant est délogé du milieu familial d'origine dans l'attente de le réintégrer, dès lors que la situation familiale est jugée adéquate si la sécurité physique, psychique, morale et affective de l'enfant y est garantie. L'enfant, dans un entredeux, est au risque d'être suspendu à un temps du présent, dans la nostalgie d'un passé qu'il n'a plus, et d'un futur incertain auprès de la famille. Comment lui permettre alors de composer avec l'absence de sa famille dans son quotidien ? La mise à mal des repères fragilise la situation psychoaffective de l'enfant placé, les processus de séparation/individuation sont alors en difficulté, surtout dans les cas où l'enfant est confronté à une multitude de lieux d'accueil. Comment peut-il donc se reconnaître comme un maillon d'une chaîne généalogique, d'appartenance ou d'affiliation, lorsque son vécu est marqué par des lieux d'accueil différents et des professionnels nouveaux ? J'ai ainsi pu rencontrer un certain nombre d'enfants qui n'avaient plus de souvenirs de leur vécu. Cette amnésie peut s'expliquer par les changements de lieux qu'ils ont connus, faute de repères stables dans leur environnement. Nous pouvons constater chez eux de nombreuses difficultés de symbolisation rencontrées chez les enfants placés à l'ASE. La difficulté à se faire une histoire de leur vécu met à mal le contrat narcissique (Aulagnier, 1975), pourtant essentiel pour tout processus d'identification, et ainsi de subjectivation. Le contrat narcissique désigne l'espace psychique d'articulation entre anticipation parentale et anticipation sociale, comme pré-investissement où va venir s'inscrire le sujet. En tant que professionnel, à travers les écrits qui parlent des enfants placés, il sera important de faire la part entre ce qui est naturel et ce qui est institué, afin de permettre aux enfants placés de se constituer au travers de l'articulation de ces deux registres, sans que l'un annule l'autre, ou désavoue l'autre. Ce travail de différenciation nécessite un travail pluridisciplinaire au sein même de l'équipe dont il est question, parler d'un enfant ne peut être celui d'une parole unique. La pluralité des discours permet une discussion de la situation de l'enfant. Le travail auprès des enfants confiés est souvent relaté dans des écrits éducatifs, qui émanent des référents éducatifs. Or, les psychologues de l'institution, de par leur regard porté sur la dimension psychique de l'enfant, permettent un nouage avec l'éducatif, qui est essentiel à une compréhension complexifiée de l'enfant pris en charge par l'ASE.

Les professionnels éducateurs sont contraints à établir un rapport pour chaque suivi, alors que les psychologues sont libres d'adresser une note dès lors qu'ils le jugent nécessaire par rapport à la situation d'un mineur confié à l'ASE. Même si l'écrit professionnel est un écrit procédural, puisqu'il s'inscrit au sein d'une procédure de mesure éducative, la possibilité de le rédiger pour le psychologue, s'il le juge nécessaire dans le travail auprès de l'enfant, en fait un outil clinique. Pour Michel Chauvière (1999), dans le travail social, « clinique sociale » signifie métaphoriquement : « au chevet du client », c'est-à-dire des personnes, des groupes ou des quartiers difficiles, au plus près des gens. Les interventions sont aussi faites d'observations et de dialogues, qui ne peuvent être réduits à l'accompagnement social.

Les institutions qui travaillent auprès de ces « ratages » de la socialisation sont dès lors définissables comme les lieux de la « mésinscription », selon Alain Noël [START_REF] Henri | Le secret de famille et l'enfant improbable[END_REF]. Le service de l'ASE peut être désigné comme une des « institutions de la mésinscription » (Pinel, Gaillard, 2015). Selon J.-P. Pinel et G. Gaillard, il a alors fondamentalement une fonction de « remaillage de la trame symbolique », en générant des processus qui visent à lier, à différencier, à métaboliser et à symboliser. Tout professionnel confronté aux déliaisons pathologiques, comme cela est le cas dans les prises en charge de l'ASE, a besoin d'être écouté, compris et soutenu par ses collègues. Pour pouvoir contenir un enfant, et contenir en soi-même les affects que son contact suscite, il faut pouvoir bénéficier soi-même de contenance. La prise en charge d'un jeune procède d'un engagement groupal, avec le travail en équipe, principalement. L'enveloppe groupale vise à soutenir le professionnel et à protéger le lien professionnel-enfant. Elle exerce en cela une fonction de méta-garance [START_REF] Eiguer | Le tiers : psychanalyse de l'intersubjectivité[END_REF].

Avant de conclure cet article, je voudrais aussi souligner la fonction de l'écrit professionnel comme la possibilité d'inscrire le travail du professionnel dans un maillage institutionnel qui le surpasse. Élément de l'histoire des interventions professionnelles, le dossier ASE sera un élément du vécu de l'enfant durant son placement. Il s'agit d'un écrit qui s'inscrit dans un ensemble d'autres écrits précédents, et qui s'adresse à d'autres, ce qui lui confère une prise en compte de l'altérité et de la fonction tierce (exercée par le juge, par

La loi n° 78-753 du 17 juillet 1978 prévoit, dans son article 4, que « l'accès aux documents administratifs s'exerce : a) par consultation gratuite sur place… b)… par la délivrance d'une copie facilement intelligible… ». Seuls peuvent être communiqués les documents administratifs et médicaux :-Documents administratifs : ce sont tous ceux qui ont été élaborés par l'ASE en l'absence d'intervention d'un juge ou avant l'intervention d'un juge ; la Commission d'accès aux documents administratifs (CADA) considère que même les rapports transmis au juge mais élaborés avant l'intervention de celui-ci sont de nature administrative.-Les photographies sont remises de droit.

Article 350 du Code civil : l'enfant recueilli par un particulier, une oeuvre privée ou un service d'Aide sociale à l'enfance, dont les parents se sont manifestement désintéressés pendant l'année qui précède l'introduction de la demande en déclaration d'abandon, peut être déclaré abandonné par le tribunal de Grande instance.

Loi n°84-422 du 6 juin 1984 relative aux droits des familles dans leurs rapports avec les services chargés de la protection de la famille et de l'enfance, et au statut des pupilles de l'État : le deuxième alinéa de l'article 32 estime que l'assistant familial est prioritaire pour l'adoption d'un enfant qui est accueilli lorsque des liens affectifs ont été noués.

Prénom modifié pour des raisons de confidentialité.

Aide sociale à l'enfance, contrat narcissique, individuation, alliances inconscientes, pacte dénégatif, mésinscription, filiation.


Shuyin Yu

Qingfeng Zeng
email: qfzeng@nwpu.edu.cn

Artem R Oganov

Gilles Frapper

Litong Zhang

Phase stability, chemical bonding and mechanical properties of titanium nitrides: A first-principles study

Keywords: crystal structure prediction, Ti-N system, electronic properties, hardness

We have performed first-principles evolutionary searches for all stable titanium nitrides and have found, in addition to the well-known rocksalt-type TiN, new ground states Ti 3 N 2 , Ti 4 N 3 , Ti 6 N 5 at atmospheric pressure, and Ti 2 N and TiN 2 at higher pressures. The latest nitrogen-rich structure presents encapsulated N 2 dumbbells with a N-N distance of 1.348 Å at 60 GPa and TiN 2 is predicted to be mechanically stable (quenchable). Our calculations of the mechanical properties (bulk modulus, shear modulus, Young's modulus, Poisson's ratio, and hardness) are in excellent agreement with the available experimental data and show that the hardness of titanium nitrides increases with increasing nitrogen content. The hardness of titanium nitrides is enhanced by strengthening directional covalent bonds and disappearance of Ti-Ti metallic bonds. Among the predicted compounds, TiN 2 has the highest hardness of 27.2 GPa.

Introduction

The discovery of new ultra-incompressible, superhard materials with novel mechanical and electronic properties is of great fundamental interest and practical importance. Conventional superhard materials are often formed by light elements (B, C, N, O), such as diamond [START_REF] Occelli | Properties of diamond under hydrostatic pressures up to 140 GPa[END_REF], γ-B 28 [START_REF] Oganov | Boron: a hunt for superhard polymorphs[END_REF], c-BN [START_REF] Zheng | Superhard hexagonal transition metal and its carbide and nitride: Os, OsC, and OsN[END_REF], B 6 O [START_REF] He | Boron suboxide: As hard as cubic boron nitride[END_REF], and BC 5 [START_REF] Solozhenko | Ultimate metastable solubility of boron in diamond: synthesis of superhard diamondlike BC 5[END_REF]. An alternative is to combine these elements with transition metals [START_REF] Holleck | Material selection for hard coatings[END_REF][START_REF] Kaner | Designing superhard materials[END_REF], e.g., OsB 2 [START_REF] Cumberland | Osmium diboride, an ultra-incompressible, hard material[END_REF], ReB 2 [START_REF] Chung | Synthesis of ultra-incompressible superhard rhenium diboride at ambient pressure[END_REF], IrN 2 and OsN 2 [START_REF] Young | Synthesis of novel transition metal nitrides IrN 2 and OsN 2[END_REF]. Here we explore the Ti-N system, which gives promise for materials with superior properties: chemical stability, thermal stability, oxidative resistance, good adhesion to substrate, high fracture toughness and high hardness, and these properties make them suitable for many applications, as a thin film or coating for highspeed cutting tools.

TiN was first separated by Story-Maskelyne [START_REF] Bannister | Osbornite, meteoritic titanium nitride[END_REF] from a meteorite, followed by numerous studies on the film growth [START_REF] Dubois | Model studies of low temperature titanium nitride thin film growth[END_REF][START_REF] Dubois | Infrared studies of the surface and gas phase reactions leading to the growth of titanium nitride thin films from tetrakis (dimethylamido) titanium and ammonia[END_REF][START_REF] Sherman | Growth and properties of LPCVD titanium nitride as a diffusion barrier for silicon device technology[END_REF], elasticity and mechanical properties [START_REF] Liu | Structural and elastic properties of TiN under high pressure[END_REF][START_REF] Brik | First-principles studies of the electronic and elastic properties of metal nitrides XN (X= Sc, Ti, V, Cr, Zr, Nb)[END_REF][START_REF] Yang | First-principles calculations of mechanical properties of TiC and TiN[END_REF]. The crystal structure of TiN has been intensively investigated, however its phase transition under pressure is still a controversy [START_REF] Zhao | Isostructural phase transition of TiN under high pressure[END_REF][START_REF] Ojha | Pressure-induced structural phase transformation and elastic properties of transition metal mononitrides[END_REF][START_REF] Chauhan | Structural stability of TiO and TiN under high pressure[END_REF]. The other well-known compounds are Ti 2 N and Ti 3 N 4 . ε-Ti 2 N [START_REF] Holmberg | Structural studies on the titanium-nitrogen system[END_REF] (P4 2 /mnm) is the most stable structure at normal conditions, while the δ' phase (I4 1 /amd) can only exist at high temperatures [START_REF] De Novion | Order and disorder in transition metal carbides and nitrides: experimental and theoretical aspects[END_REF]. Ivashchenko et al. [START_REF] Ivashchenko | Firstprinciples study of phase stability of Ti 2 N under pressure[END_REF] predicted phase transformation sequence ε-Ti 2 N → Au 2 Te-type → Al 2 Cu-type with transitions occurred at 77.5 and 86.7 GPa, respectively. For the metastable Ti 3 N 4 , Kroll et al. [START_REF] Kroll | Assessment of the Hf-N, Zr-N and Ti-N phase diagrams at high pressures and temperatures: balancing between MN and M 3 N 4 (M= Hf, Zr, Ti)[END_REF] proposed that it adopts the CaTi 2 O 4 -type structure at normal conditions. With increasing pressure, it was proposed to transform first into the Zr 3 N 4 -type structure, and then into the cubic Th 3 P 4type structure.

Although titanium nitrides have been studied for decades, a detailed theoretical study of their phase equilibria, electronic and mechanical properties would help to reconcile controversies and possibly predict new technologically useful materials. This is exactly the purpose of the present paper. We hope this study will provide guidance for experimental groups aiming to synthesize those novel crystal structures.

Computational Methodology

To find all potential Ti x N y structures in the binary Ti-N phase diagram, we performed an extensive computational search for the most stable ground-state structures with the lowest enthalpy at selected pressures (P=1 atm, 20 and 60 GPa), using the evolutionary algorithm approach as implemented in the USPEX code [START_REF] Oganov | Crystal structure prediction using ab initio evolutionary techniques: Principles and applications[END_REF][START_REF] Oganov | How evolutionary crystal structure prediction works and why[END_REF] in its variable-composition mode [START_REF] Oganov | Evolutionary crystal structure prediction as a method for the discovery of minerals and materials[END_REF], interfaced with the VASP density functional package [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF]. This approach of crystal structure prediction is based on purely chemical composition with no experimental input required. In our searches, we allowed all possible compositions in the Ti-N system with structures containing up to 24 atoms in the unit cell. The first generation of structures was produced randomly, and the succeeding generations were obtained by applying heredity (50%), atom transmutation (20%), lattice mutation (15%) operations and produced randomly (15%).

The first-principles electronic structure calculations were carried out within the framework of density functional theory (DFT) with exchange-correlation functionals in the generalized gradient approximation (GGA) according to Perdew-Burke-Ernzerhof (PBE) parameterization [START_REF] Perdew | Generalized gradient approximation made simple[END_REF]. The interactions between the ions and the electrons are described by the projectoraugmented wave (PAW) method [START_REF] Blöchl | Projector augmented-wave method[END_REF] with a cutoff energy of 600 eV. Structures were optimized until the maximum residual force component was smaller than 1 meV/Å. All the enthalpy calculations are well converged to an accuracy of <10 -6 eV/atom. Brillouin zone was sampled by uniform Γ-centered Monkhorst-Pack meshes [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF] with the resolution 2π×0.06 Å -1 . The most stable structures were studied further with increased accuracy. k-points meshes of resolution better than 2π×0.03 Å -1 were adopted.

Theoretical phonon spectra were calculated based on the supercell method using the PHONOPY package [START_REF] Togo | First-principles calculations of the ferroelastic transition between rutile-type and CaCl 2 -type SiO 2 at high pressures[END_REF] in order to probe the dynamic stability of the predicted Ti x N y compounds at different pressures. The elastic constants were calculated from the stress-strain relations, and Voigt-Reuss-Hill (VRH) approximation [START_REF] Hill | The elastic behaviour of a crystalline aggregate, Pro-ceedings of the Physical Society[END_REF][START_REF] Voigt | Lehrbuch derkristallphysik[END_REF][START_REF] Reuss | A calculation of the bulk modulus of polycrystalline materials[END_REF] was employed to obtain the bulk modulus B, shear modulus G, Young's modulus E, and Poisson's ratio ν. The theoretical Vickers hardness H v was estimated by using the Chen's model [START_REF] Chen | Modeling hardness of polycrystalline materials and bulk metallic glasses[END_REF], according to the expression:

H v = 2(κ 2 G) 0.585 -3 ( 1
)
where κ is the Pugh ratio: κ=G/B. The directional dependence of the Young's modulus for crystals of different symmetries encountered in this study is given in Supplementary Materials.

Results and Discussion

Thermodynamics of titanium nitrides can be quantified by constructing the thermodynamic convex hull, which is a complete set of phases stable against transformation into any other phase and decomposition into any set of other phases. These thermodynamically stable phases can be synthesized in principle [START_REF] Ghosh | First-principles calculations of the structural and thermodynamic properties of bcc, fcc and hcp solid solutions in the Al-TM (TM=Ti, Zr and Hf) systems: A comparison of cluster expansion and supercell methods[END_REF]. The convex hull curves were reconstructed after phonon calculations for the selected thermodynamically stable phases, which are shown in Figure 1. Noticeably, our evolutionary simulations succeed in finding the well-known NaCl-type TiN and ε-Ti 2 N observed by experiments at ambient conditions. In addition, we uncovered three new stable compounds, the orthorhombic Ti 3 N 2 , and monoclinic Ti 4 N 3 and Ti 6 N 5 . At 60 GPa, ε-Ti 2 N transforms into the orthorhombic Cmcm phase, and a new stable nitrogen-rich compound is identified as TiN 2 . The detailed crystallographic data and enthalpies of formation are listed in Table 1s (see Supplementary Materials). Note that the enthalpies of formation without the ZPE correction are only ∼0.01 eV/atom different from the ZPE-corrected energies, and it will not change the relative phase stabilities without considering ZPE correction. Thus, for the computational intensity reason, the zero point energy is omitted for the systematic structural search.

Rocksalt TiN and related subnitrides Ti n+1 N n

While TiN has the ideal cubic rocksalt structure, Ti 3 N 2 , Ti 4 N 3 and Ti 6 N 5 are versions of this structure with ordered N-vacancies (see Figure 2) -in Ti 3 N 2 one third of nitrogen sites are vacant, in Ti 4 N 3 -one quarter, and in Ti 6 N 5 -one sixth. Similar vacancy-ordered phases were earlier reported to be stable for hafnium and titanium carbides M n+1 C n (M=Hf and n=2, 5 [START_REF] Zeng | Prediction of stable hafnium carbides: Stoichiometries, mechanical properties, and electronic structure[END_REF]; M=Ti and n=1, 2, 5 [START_REF] Jiang | Pressure-composition phase diagram of Ti-C from first principles[END_REF]). Generally, the decrease of vacancy concentration makes the structure denser, less compressible, harder, and more stable under pressure. Indeed, among defective rocksalt-type phases only Ti 6 N 5 (with the lowest vacancy concentration) survives as a (barely) stable phase at 60 GPa. Ti 3 N 2 has space group Immm, while Ti 4 N 3 and Ti ε-Ti 2 N has the anti-rutile structure with space group P4 2 /mnm. At pressures above 20.8 GPa, tetragonal ε-Ti 2 N transforms into an orthorhombic Cmcm structure (see Figure 2e). This structure contains a distorted close packing of Ti atoms, where one half of octahedral voids are occupied by nitrogen atoms, and edge-sharing NTi 6octahedra form double slabs alternating with nitrogen-free layers. It is noteworthy that the enthalpy of Cmcm-Ti 2 N is lower by more than 0.24 eV/f.u. than that of the Au 2 Te-type or Al 2 Cu-type structures proposed by Ivashchenko et al. [START_REF] Ivashchenko | Firstprinciples study of phase stability of Ti 2 N under pressure[END_REF]. Having investigated the crystal structures and energetics of the various titanium subnitrides Ti n+1 N n and TiN, we now consider their electronic properties. The total and partial densities of states (DOS) of the understudied compounds are shown in Figure 3(a). The crystal orbital Hamilton population (-COHP) curves for all Ti-N compounds were calculated by using the LOB-STER package [START_REF] Dronskowski | Crystal orbital hamilton populations (COHP): energy-resolved visualization of chemical bonding in solids based on density-functional calculations[END_REF] and are displayed in Figure 3(b). These orbital-pair interactions can provide a quantitative measure of bond strengths. The positive and negative energy regions in the -COHP curves correspond to bonding and antibonding states, respectively.

For Ti n+1 N n structures, the DOS are decomposed into three well separated regions energy: (1) a deeply lowest valence band, s N ; (2) a hybridized Ti-3d/N-2p band, d M p N ; (3) a partially filled higher lying Ti-3d band, d M . The s N band is dominated by the 2s orbitals of the nitrogen atoms. Therefore, the contribution of this band to the bonding is not so large. The following upper group of valence bands d M p N is a result of strong hybridization from the 3d states of Ti atoms and the 2p states of N atoms. These peaks correspond to the Ti-3d/N-2p bonding orbitals contribution. Their antibonding counterparts appear well above the Fermi level.

When n goes from 1 to 5 in Ti n+1 N n series, e.g., the number of N-filled Ti 6 octahedra increases, the energy region for the Ti-3d/N-2p peaks is expanded. This energy dispersion is caused by more extensive mixing between Ti-3d and N-2p orbitals, and suggests the enhancement of the covalency of titanium-nitrogen bonding network. Let turn our discussion on the upper energy region. From both DOS and COHP analyses, one may see that the bottom of the d M band is mainly dominated by the bonding orbitals of 3d Ti atoms and the corresponding DOS exhibits a feature of "near free electron" responsible for metallicity (see Figure 3 As N content increases, valence band broadens (d M p N band), valence electron concentration increases, and covalent interactions become stronger. This is reflected in shortening of Ti-N bonds: the shortest bond lengths are 2.082, 2.071, 2.069 and 2.022 Å for Ti 2 N, Ti 3 N 2 , Ti 4 N 3 and Ti 6 N 5 at 0 GPa, respectively. In the same sequence the N(E F ) gradually decreases, reaching the lowest value ∼0.098 states/eV/electron for TiN (see Figure 3(a)). To further explore the bonding nature in each Ti n+1 N n and TiN structures, we examined the electron localization function (ELF) [START_REF] Becke | A simple measure of electron localization in atomic and molecular systems[END_REF] of these Ti-N compounds (see Figure 4). We find that there are relatively large ELF values between Ti and N atoms, indicating the partially covalent Ti-N interactions. ELF=0.5 can be identified between titanium atoms, suggesting that the metallic Ti-Ti bonding. All of these results are consistent with the discussion about the DOS and COHP for each structure of Ti x N y . At 60 GPa, we uncovered a nitrogen-rich compound TiN 2 with the CuAl 2 -type structure (SG: I4/mcm, see Figure 5a). The energy of this I4/mcm structure is lower than that of the previously proposed thermodynamic ground state, the CaC 2 -V-type structure predicted by Kulkarni et al. [START_REF] Kulkarni | Structure prediction of binary pernitride MN 2 compounds (M= Ca, Sr, Ba, La, and Ti)[END_REF] at 60 GPa. This new structure is stable against the decomposition into the mixture of TiN and N 2 at pressures above 26.6 GPa. To examine the stability of the I4/mcm structure, the phonon frequencies were calculated at atmospheric pressure and at 60 GPa. Imaginary frequency mode was not observed in the entire Brillouin zone (see Figure 5b and Supplementary Materials), indicating that the predicted I4/mcm structure is dynamical stable and, if synthesized, might be quenchable at atmospheric pressure.

TiN 2 , a high-pressure structure with N 2 dumbbells

The structure of TiN 2 contains TiN 8 face-sharing tetragonal antiprisms stacked along the c-axis. Unlike in normal CuAl 2 -type structures, here we have N 2dumbbells encapsulated in cubic Ti 8 hexahedra. The N-N bond length is calculated as 1.348 Å and 1.378 Å at 60 GPa and 1 atm, respectively. Projected N-2p DOS shows that antibonding 1π g * levels are almost fully occupied at 60 GPa (see Figure 5c). Consequently, for electron counting purposes, the dinitrogen units should be formally considered as N 2 4-, a pernitride unit isoelectronic to fluorine F 2 molecules. Its electronic ground state configuration is 1σ g 2 , 1σ u 2 , 1π u 4 , 2σ g 2 , 1π g * 4 for 14 valence electrons and its bond order is one. Formally, this leaves the titanium atoms of TiN 2 in a d 0 configuration (Ti 4+ ). Note that a large variety of structures are known with the square antiprismatic geometry and a formally d 0 metal center such as TaF 8

3-in Na 3 TaF 8 [START_REF] Burdett | Eight-coordination[END_REF]. Stabilizing Ti-Ti interactions in binary Ti-N structures are weakening when nitrogen content increases, e.g., when formal valence d electrons d n decrease (Ti oxidation number increases). This behavior is structurally reflected by the elongation of the Ti-Ti distances upon the nitrogen content (shortest Ti-Ti bonds of 2.255 Å and 2.517 Å in Ti 2 N and TiN 2 at 60 GPa, respectively). Note that the pernitride N 2 4-units point directly towards the tetragonal faces of the eight-coordinated titanium atoms (Ti-N bonding length is 2.088 Å at 60 GPa) and are perpendicular to each other in or-der to minimize the steric clashes between the nitrogen σ-lone pairs (Pauli repulsions). Based on these considerations, the covalent networks in TiN 2 with the CuAl 2 -type structure are built on transition metal -main group elements and main-group-based directional covalent bonds, a prerequisite for superhard materials.

Mechanical properties of Ti-N stable structures

The dynamical stability of the predicted Ti-N compounds has been confirmed by the absence of any imaginary phonon modes (see Supplementary materials). The mechanical stability was examined using the calculated elastic constants (see Table 1), and we found that all these phases satisfy the Born-Huang stability criteria [START_REF] Wu | Crystal structures and elastic properties of superhard IrN 2 and IrN 3 from first principles[END_REF]. Our calculated elastic constants are in good agreement with the available data [START_REF] Liu | Structural and elastic properties of TiN under high pressure[END_REF][START_REF] Ahuja | Structural, elastic, and high-pressure properties of cubic TiC, TiN, and TiO[END_REF][START_REF] Kim | Elastic constants of single-crystal transition-metal nitride films measured by linefocus acoustic microscopy[END_REF][START_REF] Wang | Structural and elastic properties of cubic and hexagonal TiN and AlN from first-principles calculations[END_REF]. One can note the large C 33 values for Ti 2 N, Ti 4 N 3 , TiN and especially TiN 2 (691 GPa), which indicate very low compressibility along the c-axis. For TiN 2 , we also find a very large elastic constant C 44 (350 GPa).

The calculated bulk modulus B and shear modulus G are also listed in Table 1 and the trend of these mechanical properties as a function of N content is shown in Figure 6. Our calculated results for the well-known rocksalt TiN and ε-Ti 2 N are in good agreement with reported values [START_REF] Ahuja | Structural, elastic, and high-pressure properties of cubic TiC, TiN, and TiO[END_REF][START_REF] Liu | Structural and elastic properties of TiN under high pressure[END_REF][START_REF] Kim | Elastic constants of single-crystal transition-metal nitride films measured by linefocus acoustic microscopy[END_REF][START_REF] Wang | Structural and elastic properties of cubic and hexagonal TiN and AlN from first-principles calculations[END_REF], supporting the accuracy and reliability for other Ti-N compounds. TiN has the largest B of 317 GPa, slightly larger than TiN 2 of 310 GPa. Both values are very high and comparable to WB 4 with a bulk modulus of 328 GPa [START_REF] Wang | Origin of hardness in wb4 and its implications for ReB 4 , TaB 4 , MoB 4 , TcB 4 , and OsB 4[END_REF]. Bulk modulus is mainly related to valence electron density (VED) [START_REF] Gilman | Design of hard crystals[END_REF], and indeed the VED of TiN 2 (0.478 e -/Å 3 ) is very close Figure 7 shows the directional dependence of Young's moduli of the considered Ti-N compounds, and one can clearly see anisotropy from the deviations of its shape from sphere. The Young's modulus along the caxis direction is higher than that along the other directions for Ti 2 N, Ti 4 N 3 , and TiN 2 , which agrees well with our above analysis on elastic constants. The Young's modulus is more isotropic in Ti 4 N 3 , while TiN 2 shows more anisotropy than other Ti-N compounds, because of its high C 44 .

Poisson's ratio ν and G/B ratio [START_REF] Pugh | Relations between the elastic moduli and the plastic properties of polycrystalline pure metals[END_REF] are indicative of the degree of directionality of the covalent bonding. A typical value for ν is ∼0.2 for strong directional covalent materials and ∼0.4 for good metals. A low Poisson's ratio results from directional bonds, which increase the shear modulus and limit the motion of dislocations, thereby increasing a material's hardness. From Table 1, we find that ν drops down from 0.276 in Ti 2 N to 0.218 in TiN 2 with increasing N content, except Ti 4 N 3 , which shows a local maximum. For Ti 4 N 3 , we can find that there are large Ti-Ti metallic bonding in contrast with other Ti-N compounds (see Figure 3(b)). The large Ti-Ti metallic bonds impair the covalency of this phase, which lead to a local minimum of hardness (see Figure 6f). The Vickers hardness was estimated by Chen's empirical method and listed in Table 1. Among the studied Ti-N compounds, TiN 2 has the largest hardness of 27.2 GPa. Note that this is an orientationally averaged value (anisotropy is not considered in Equation 1). It is no wonder that increasing nitrogen content changes bonding to more directional, the more N is added, the more Ti-N bonds are formed and the less Ti-Ti bonding remains. Together with the analyses on electronic properties and chemical bonding in those stable Ti-N compounds, We can conclude that with increasing N content, the enhancement of directional covalent interactions and decline of metallicity lead to the increase of the hardness.

Conclusions

In summary, we have extensively explored the stable structures and possible stoichiometries in the Ti-N system by first-principles evolutionary crystal structure prediction. In addition to the well-known phases Ti 2 N and TiN, we have uncovered three new intriguing structures at ambient conditions (Ti 3 N 2 , Ti 4 N 3 and Ti 6 N 5 ). At high pressures, two new phases Cmcm-Ti 2 N and I4/mcm-TiN 2 were discovered. All those phases are mechanically and dynamically stable at ambient conditions. The calculated elastic constants C i j of Ti-N compounds are in good agreement with the available reported data. Other mechanical properties including bulk modulus, shear modulus, Young's modulus, and Poisson's ratio have been further computed from C i j . Among the studied Ti-N compounds, TiN 2 exhibits significant elastic anisotropy and possess the highest hard-ness of 27.2 GPa. We found a strong correlation between the mechanical properties and N content, e.g., the more N content, the more directional covalent Ti-N bonding and the less Ti-Ti metallic bonding, which lead to the enhancement of the hardness. The materials discovered here are attractive for technological applications because of a compromise between hardness and ductility, due to a peculiar interplay between metallicity and covalency.

Figure 1 :

 1 Figure 1: (Color online). Convex hull diagrams for Ti-N system at ambient pressure (a) and at 60 GPa (b). Circles denote stable (solid circles) and metastable (open circles) structures. The hexagonal ω-Ti [38] and α-N 2 [39] are adopted as reference states.

Figure 2 :

 2 Figure 2: (Color online). Crystal structures of subnitrides Ti-N compounds. (a) Immm-Ti 3 N 2 at 0 GPa, (b) C2/m-Ti 4 N 3 at 0 GPa, (c) C2/m-Ti 6 N 5 at 0 GPa, (d) Cmcm-Ti 2 N at 60 GPa. The Ti and N atoms are represented as big blue and small gray spheres, respectively, nitrogen vacancies are shown in the green octahedra.

  6 N 5 belong to space group C2/m. It is convenient to visualize these structures by N-centered octahedral (NTi 6 ), as shown in Figure 2. In the Ti n+1 N n series, Ti 5 N 4 is missing as a stable compound. Nevertheless, a structure (SG: I4/m, No. 87) is very close to the convex hull curve but lies above it (15 meV/atom), i.e., Ti 5 N 4 is a metastable phase.

Figure 3 :

 3 Figure 3: (Color online). (a) The total and partial density of states (DOS), (b) Crystal orbital Hamilton population (-COHP) for the Ti n+1 N n subnitrides and TiN at zero pressure. Vertical dashed line is the Fermi energy.

  (a)). Vacant antibonding Ti-3d/N-2p and antibonding Ti-3d/N-3d states compose the bottom of the conduction band. These results confirm the mixed covalent Ti-N and metallic Ti-Ti bonding nature in titanium subnitrides and TiN compounds.

Figure 4 :

 4 Figure 4: (Color online). Calculated electron localization function (ELF) maps for (a) P4 2 /mnm-Ti 2 N, (b) Immm-Ti 3 N 2 , (c) C2/m-Ti 4 N 3 , (d) C2/m-Ti 6 N 5 , (e) Fm 3m-TiN, and (f) I4/mcm-TiN 2 , blue spheres represent Ti atoms while gray spheres represent N atoms.

Figure 5 :

 5 Figure 5: (Color online). (a) Crystal structure of I4/mcm-TiN 2 , (b) phonon dispersion curve at 60 GPa, (c) projected and total DOS with symmetry labels of molecular N 2 4-filled valence orbitals assigned to each DOS peak at 60 GPa, and (d) -COHP of I4/mcm-TiN 2 at 60 GPa.

Figure 6 :

 6 Figure 6: (Color online). Calculated bulk modulus B, shear modulus G, Young's modulus E, Poisson's ratio ν, G/B ratio, and Vickers hardness H v of titamium nitrides as a function of N content. The structures from left to right are ε-Ti 2 N, Immm-Ti 3 N 2 , C2/m-Ti 4 N 3 , C2/m-Ti 6 N 5 , Fm 3m-TiN, and I4/mcm-TiN 2 .

Figure 7 :

 7 Figure 7: (Color online). Directional dependence of Young's moduli (in GPa) for (a) P4 2 /mnm-Ti 2 N, (b) Immm-Ti 3 N 2 , (c) C2/m-Ti 4 N 3 , (d) C2/m-Ti 6 N 5 , (e) Fm 3m-TiN, and (e) I4/mcm-TiN 2 .

Table 1 :

 1 The calculated elastic constants C i j , bulk modulus B (GPa), shear modulus G (GPa), Young's modulus E (GPa), Poisson's ratio ν, κ=G/B ratio, and Vickers hardness H v (GPa) of the Ti-N compounds at 0 GPa.

	Phase		C 11	C 12	C 13	C 16	C 22	C 23	C 26	C 33	C 36	C 44	C 45	C 55	C 66	B	G	E	ν	κ	H v
	Ti 2 N	GGA	300 214 115					442		157			137 214 112 287 0.277 0.523 11.87
	(P4 2 /mnm) LDA	353 235 138					489		184			147 246 130 332 0.275 0.528 13.38
		Ave	327 225 127					467		171			142 230 121 310 0.276 0.526 12.63
	Ti 3 N 2	GGA	458 99	125		400 115		332		119		145 86	206 124 309 0.250 0.602 15.53
	(Immm)	LDA	545 113 140		449 140		370		127		172 92	237 138 347 0.256 0.583 16.00
		Ave	502 106 133		425 128		351		123		159 89	222 131 328 0.253 0.593 15.77
	Ti 4 N 3	GGA	368 150 155 -17	393 119 3	421 -29	140 10	126 114 225 125 317 0.265 0.556 14.05
	(C2/m)	LDA	412 181 179 -16	444 135 8	484 -42	167 19	139 130 258 141 357 0.269 0.547 14.84
		Ave	390 166 167 -17	419 127 6	453 -36	154 15	133 122 242 133 337 0.267 0.552 14.45
	Ti 6 N 5	GGA	424 140 145 -20	429 146 20	429 1	147 21	158 171 238 151 373 0.239 0.634 18.97
	(C2/m)	LDA	484 161 175 -33	490 175 33	481 1	166 34	187 201 275 170 424 0.243 0.618 20.07
		Ave	454 151 160 -27	460 161 27	455 1	157 28	173 186 257 161 399 0.241 0.626 19.52
	TiN	GGA	590 145								169				294 189 466 0.235 0.643 22.55
	(Fm 3m)	LDA	704 157								183				339 215 533 0.238 0.634 24.22
		Ave	647 151								176				317 202 500 0.237 0.639 23.39
		GGA[15] 583 129								179				278					24.19[45]
		LDA[15] 698 139								198				321
		Exp[46] 625 165								163				320					23.00[47]
	TiN 2	GGA	535 279 71					653		336			148 284 197 481 0.218 0.693 25.64
	(I4/mcm) LDA	631 301 106					729		364			185 335 233 567 0.218 0.696 28.75
		Ave	583 290 89					691		350			167 310 215 524 0.218 0.695 27.20

Acknowledgements

We thank the Natural Science Foundation of China (Grants No. 51372203 and No. 51332004), the Basic Research Foundation of NWPU (Grant No. JCY20130114), the Foreign Talents Introduction and Academic Exchange Program (Grant No. B08040), the National Science Foundation (Grants No. EAR-1114313 and No. DMR-1231586), DARPA (Grants No. W31P4Q1310005 and No. W31P4Q1210008), NSW (grant DMR-1231586) and the Government of the Russian Federation (Grant No. 14.A12.31.0003) for financial support. The authors also acknowledge the High Performance Computing Center of NWPU for the allocation of computing time on their machines.


Somayeh Faraji
email: sfaraji86@gmail.com

Busheng Wang

Hubert Okadome Valencia

Gilles Frapper
email: gilles.frapper@univ-poitiers.fr

Computational Discovery of 2D Copper Chalcogenides CuX (S, Se, Te)

We investigate the energy landscape of 2D CuS, CuSe, and CuTe compounds using a global evolutionary algorithm in combination with density functional theory (DFT). Four low-lying energy Cu2X2 (X=S, Se, Te) slabs with P 3m1, P 4/mmm, P 4/nmm, and P mmn symmetries have been identified on their respective potential energy surfaces. Three structures present tetrahedral CuX4 motifs which pave the 2D slabs, while the latter is based on shared Cu4X2 octahedra. The viability of each phase was examined by looking at dynamical, thermodynamic, and thermal stability criteria. We find that 2D copper monosulfide crystallizes in the P 3m1 phase, reminiscent of the covalent slab found in the CaAl2Si2 prototype. Two polymorphs of 2D copper selenide with symmetries P 3m1 and P 4/mmm are close in energy. The latter Cu2Se2 slab is made of fused Cu4Se2 square bipyramids with electronegative Se atoms in apical positions. Finally, the ground-state 2D-CuTe has a P mmn space group containing distorted CuTe4 tetrahedra with some Te-Te bonding. The electronic and bond analyses show that all of these predicted 2D CuX phases are metallic with iono-covalent bonds.

I. INTRODUCTION

Two-dimensional (2D) materials present a new class of materials with novel physics which have attracted significant interest in recent years, particularly since the discovery of graphene as the first room-temperature 2D material. 1 Such 2D materials exhibit a wide range of extraordinary optical, mechanical, and electronic properties which are markedly different from their bulk counterparts. Owing to such properties, they have attracted the interest of many experimental and theoretical researchers to keep searching for novel 2D materials and understand their specific physical and chemical properties. While the focus of most of the initial researches has been concentrated on 2D semiconductors, 2 a growing interest has been shown to metallic 2D materials, especially transition metal mono-chalcogenides MX (M = transition metal like Cu, Fe, Pb, etc.; X=S, Se, Te). [3][4][5] Binary copper chalcogenides CuX (X = S, Se, and Te) are potential materials with applications in photothermal therapy 6 , batteries [7][8][9] , and room-temperature sensors 10 . Copper sulfide nanocrystals have received special attention due to their unique optoelectronic properties. 11 Among them, CuS, known as covellite, is one of the most studied copper sulfides due to its large ion diffusion and adsorption band gap. It is, therefore, a good candidate for technological applications such as ammonia gas sensor 10 , anticancer drugs 12 , catalysts 9,13 , glucose biosensors 14 , and electrode materials for rechargeable Li-ion 11,15 or Mg-based 16 batteries. The copperselenide system presents many phases with potential applications for photothermal therapy 6 , thermoelectric materials 17 , or electrodes for supercapacitors 18 . Finally, the copper-telluride system is viewed as the most complex copper-chalcogenide, owing to its several nonidentical polymorphs for different Cu:Te ratios [19][20][21] . Its large thermoelectric power has attracted the interest of many researchers. 22,23 It can also be viewed as a quasi-one-dimensional charge density wave system 24 which is an intriguing fundamental phenomenon. 25 In addition to the studies briefly introduced above on bulk CuX materials, their low-dimensional structures have also drawn a lot of attention. In 2012, a singlecrystal nanosheet of CuS with a thickness of 3.2 nm was fabricated 11 and used for rechargeable magnesium batteries 16 . Recently, a computational work based on DFT calculations has demonstrated that monolayer graphenelike structure CuS was unstable and the three-layered CuS (3L-CuS) with a thickness of 0.773 nm was predicted to be intrinsically stable. 26 Besides, a CuS structure with flower-like microspheres was synthesized 27 . It was shown that its loose shell was constructed from a graphene-like CuS nanosheet. To our knowledge, the 3L-CuS is therefore proposed to be the thinnest stable free-standing CuS nanosheet with covellite stoichiometry. The structural, optical, and electronic properties of thin films of CuSe were also investigated. 28 The fabrication of a monolayer of CuSe grown on a Cu substrate was first reported in 2017. 29 It presents periodic patterns of triangular nanopores with uniform size and was reported as a honeycomb lattice of CuSe layers. A similar 1D moiré pattern of planar CuSe monolayers was then synthesized. 30 Then, it was shown that there is a strong electrical anisotropy and current-voltage curves along the zigzag and armchair directions of this graphenelike structure. 31 Thin films of copper telluride have been experimentally synthesized by many groups 32,33 and are used as a back-contact for CdTe-based solar cells. 34 The latter structure was characterized by XRD as a polycrystalline structure, with preferred (101), (112), and (020) orientations. 33 Despite experimental breakthrough in copper chalcogenides chemistry, the discovery of novel 2D materials is usually time-consuming and expensive. In this respect, crystal structure prediction (CSP) searches emerge as an efficient approach to accelerate the discovery and de-sign of new 2D materials. 35,36 To extend the 2D CuX based structures family, we performed extensive structural searches of 2D CuX by employing an evolutionary algorithm in combination with DFT calculations.

In this paper, we report the discovery of six new phases of 2D CuX materials (X=S, Se, Te) from only the knowledge of their chemical composition. These phases are named 2D CuS-P 4/mmm, CuS-P 4/nmm, CuS-P mmn, CuSe-P 4/mmm, CuSe-P mmn, and CuTe-P 4/nmm. We investigated their viability 37 by the DFT study of their thermodynamic, dynamical, and thermal stabilities of twelve 2D CuX structures within P 3m1, P 4/nmm, P 4/mmm, and P mmn space groups as well as their bonding and electronic properties.

The present work is organized as follows. Section II provides information on the employed computational tools and the methods such as the configurational search approach, DFT calculations, and the main criteria for establishing the stability of the proposed 2D structures. The 2D evolutionary algorithm crystal structure search results are discussed in section III A. The energetic and thermal stabilities of the four selected 2D CuX phases are discussed in section III B. The chemical bonding and electronic properties are described in section III C. Our conclusion and final remarks are provided in the last section.

II. METHODS AND COMPUTATIONAL DETAILS

A. Crystal structure prediction searches

Configurational exploration methods based on computation enable the pre-screening for 2D materials by rapidly and accurately identifying 2D material candidates. In this study, we first examined the energy landscape of each 2D CuX (X=S, Se, Te) compounds by a series of searches using a fixed-composition evolutionary algorithm (EA) as implemented in the Universal Structure Predictor: Evolutionary Xtallography (USPEX) software. [38][39][40] We used an EA specifically designed in the USPEX code for 2D CSP. During the CSP search, the geometry optimization and total energy calculations of a 2D crystalline structure are carried out using DFT as implemented in the Vienna Ab initio Simulation Package (VASP version 5.4.4) 41,42 . The first-principles methodology is described in the next section (Sec.II B). Primitive cells containing up to 8 atoms are considered. In the first iteration of USPEX (i.e. the first EA generation), we chose a large number of 80 structures to be generated randomly. All subsequent generations contain at least 60 structures produced from the previous generation, using evolutionary variation operators defined in Refs 38-40. A ratio of 50%, 30%, and 20% is chosen for the generation of these 60 structures using heredity, random planar symmetry, and transmutation (of Cu and X atoms), respectively. In 3D space, layers of 2D materials are sepa-rated by a vacuum thickness of 15 Å to avoid interlayer interactions. In USPEX calculation, five successive steps of increasing accuracy are employed to improve the geometry convergence. The ab initio methodology corresponding to the latest step is described in Section II B. Then, all 2D CuX structures within 200 meV/atom of the convex-hull (formed by the structure with the most favorable total energy) are considered as possible thermodynamic (meta)stable phases. The geometrical and electronic properties of them are investigated in the following.

B. DFT methodology

The first-principles calculations were performed using a projected-augmented-wave (PAW) 43 method as implemented in the VASP package. The Perdew-Burke-Ernzerhof (PBE) 44 exchange-correlation functional was used within the generalized gradient approximation (GGA) 45 . A cutoff energy of 520 eV was used for all the systems. All of the 2D structures were optimized until the net forces on the atoms were below 1 meV/ Å and the total stress tensor deviated by ≤ 0.01 GPa, resulting in enthalpies that converged to less than 1 meV/atom (lower than a chemical accuracy of 1 kcal/mol, i.e., 0.04 eV/atom). The Brillouin zone was sampled using a Γ-centered, uniform k-point grid with a resolution of 2π × 0.04 Å-1 for all selected structures. All EA searches are done at the PBE level of theory: the enthalpy at zero Kelvin is computed at the same PBE-level methodology that was used to optimize the crystal structure. For the energetic aspect, the Strongly Constrained and Appropriately Normed (SCAN) meta-GGA functional 46 is also used. Finally, to determine the electronic properties, the Heyd-Scuseria-Ernzerhof (HSE06) functional 47 was used on the PBE-level optimized structure and referred to as HSE06//PBE. Because d-splitting is sensitive to the quantum mechanical methodologies, the electronic structures of all the predicted 2D CuX were also obtained by using Hubbard model, i.e. PBE+U 48,49 and employing U=4 eV. Furthermore, the magnetism was checked to determine the preferred magnetic ground for our 2D CuX, as presented in Table S2 of the Supporting Information (SI). All 2D CuX phases are nonmagnetic, therefore, spin polarization was not taken into account in this work.

C. Criteria of stability

In performing a CSP search for novel materials, an essential aspect is to judge the stability of the predicted phases. For this purpose, three main criteria have been used in our study to examine the dynamical, thermodynamical, and thermal stability of our 2D CuX candidates.

We assert the structure as dynamically stable by the absence of imaginary frequencies in the phonon dispersion diagram. The structure then corresponds to a local enthalpy minimum, i.e. a stable or metastable crystalline solid. Real-space force constants were calculated using the Density Functional Perturbation Theory (DFPT) 50 as implemented in the VASP package which relies on finite displacements. The calculations were done on supercells of the PBE-optimized structure candidates and phonon frequencies were calculated from the force constants using the Phonopy 51 code.

The thermodynamic stability of a 2D CuX structure was examined by calculating its enthalpy of formation and cohesive energy. The former was given at zero Kelvin as ∆ H f = H CuX -(H Cu + H X ). The enthalpy H was obtained by H = E e + PV, with E e , P, and V, the total electronic energy, external pressure, and structure cell volume, respectively. The bulk of Cu (F m 3m), S (P 2 1 /c), Se (P 2 1 /c), and Te (P 3 1 21) were taken as references for the calculation of H Cu , H S , H Se , and H Te , respectively.

The effect of the Zero Point Energy (ZPE) on the stability of bulk and 2D CuX compounds was also studied. Its inclusion only slightly modified the enthalpies of formation and the cohesive energies. Therefore, the energies presented in this work are not ZPE-corrected.

Besides the dynamical stability criterion, we investigated the thermal stability of each proposed CuX structure at different temperatures by performing ab initio molecular dynamic (AIMD) simulations. The chosen 2D supercells are of 2 × 2 × 1 or 3 × 3 × 1 sizes, with number of atoms up to 128. These AIMD calculations were carried out in the NVT canonical ensemble by employing the Nose-Hoover thermostat. The initial temperature was set at 1000 K, then if the 2D CuX structure decomposes, a second AIMD simulation was performed at 800 K, and so on at 600, 400 K, and finally at 300 K. The AIMD simulation time is 10 ps in total with a time step of 2 fs. The VASP software using the previously described PBE functional within GGA was used, while the Brillouin zone integration was reduced to the Γ point of the supercells, owing to the high computational cost of the AIMD simulations.

D. Bonding and electronic properties analysis

To perform chemical bonding analysis, we obtained band structures, total and projected density of states (DOS), and electron localization function from the optimized PBE geometries obtained by VASP. Additionally, crystal overlap Hamilton population, was calculated using the LOBSTER 52 package. The bond orders (BO) were computed by the Chargemol program 53 using the Manz' bond-order equation 54 with density derived electrostatic and chemical (DDEC6) partitioning 53,55 . The charge analyses were performed by employing the Bader charge analysis program of Henkelman et al. 56 Finally, the GDIS visualization software 57 which interfaced the USPEX and VASP codes was employed. In this article, images of the structures are produced using VESTA 58 software.

III. RESULTS AND DISCUSSION

A. In silico prediction of new copper chalcogenides phases

We performed extensive CSP searches of 2D CuX (X=S, Se, Te) as explained in Section II A. For each copper chalcogenide system, more than 800 structures were generated per USPEX run. These structures define a Potential Energy Surface (PES) on which we select all low-lying energetic 2D CuX structures within a 200 meV/atom energy window (see FIG. S1 of SI which illustrates some of these geometries). Following this criterion and phonon calculations, four stationary points are localized on the PES, that correspond to hexagonal P 3m1, square P 4/mmm, and P 4/nmm, and orthorhombic (centered rectangular) P mmn Bravais lattices. Their 2D crystal structures are displayed in Fig. 1 and their bonding is discussed in section III C. In the following, we discuss general structural aspects of these four geometries. Hexagonal P 3m1 slab presents the structural arrangement found in the anionic sublattice of CaAl 2 Si 2 structure type while square P 4/nmm structure is isostructural to the Cr 2 Si 2 layers found in ThCr 2 Si 2 structure type. In both P 3m1 and P 4/nmm structures, Cu atom resides at the center of an X 4 tetrahedron, and these tetrahedra pack differently to form the aforementioned types of 2D crystals. In P 3m1, each tetrahedron shares three of its six edges, while in P 4/nmm it shares four edges.

The third structure of 2D CuX, P 4/mmm, has also a square lattice. P 4/mmm and P 4/nmm slabs have both square copper and chalcogen lattices. These compounds differ from the packing of their two chalcogen monolayers, i.e. an AA packing is observed in P 4/mmm with Cu 4 X 2 square bipyramidal units sharing four Cu-Cu edges while an AB packing is present in P 4/nmm.

Finally, the fourth and last crystal structure type of 2D CuX has a centered rectangular lattice with a P mmn space group. The Cu 2 X 2 slab may be viewed as a "compressed" P 4/nmm structure along one lattice direction, leading to Cu 4 rhombi with square tilings. Here, a X 4 square net transforms into a X 4 rectangular network from P 4/nmm to P mmn 2D structures, and the copper lattice is slightly corrugated in P mmn phase.

All P 3m1, P 4/nmm and P mmn phases have CuX 4 tetrahedra (elongated in P mmn), while P 4/mmm presents CuX 4 square planar units. Finally, while P 3m1 presents three-coordinated Cu centers in its metal sublattice (nearest Cu atoms), the other 2D structures have four-coordinated Cu centers. All chalcogen atoms are connected to four nearest Cu neighbors. We will describe the calculated structural parameters in more detail further in Section III C.

We note that these four CuX slabs present a different energy ranking depending on the chalcogen element X (see Fig. 2), and a dedicated discussion will follow in section III B. We confirmed that these twelve stationary points are local minima, owing to their dynamical stability. This is shown by the absence of imaginary frequencies in their phonon dispersion curves (see FIG. S2 of the SI).

Our EA-DFT searches lead to novel 2D copper chalcogenide materials, namely P 4/mmm CuX phases (X=S, Se, Te), CuS P 4/nmm and CuSe P mmn. While other found structures have been previously proposed, i.e. CuS P 3m1 and P mmn, CuSe P 3m1 and P 4/nmm, and CuTe P 3m1, P 4/nmm, and P mmn, no detailed study of the dynamic, thermodynamic, and thermal stabilities was reported. It is therefore the topic of the next sections: the study of energetic, structural, thermal, bonding, and electronic aspects of the 12 phases of 2D copper chalcogenides discovered by our CSP search.

B. Energetic and thermal stabilities

The total PBE energies relative to the P 3m1 structure is plotted in FIG. 2 for each 2D CuX system. The order of stability of the four structures is roughly identical in both CuS and CuSe systems, P 3m1 being the ground-state 2D structure. The P 4/mmm and P 4/nmm phases are in a quasi-degenerate energy state (within our minimization accuracy), for CuS, while the difference is more pronounced but still low (8 meV/atom) for the corresponding phases of CuSe. Nevertheless, a rough estimation of activation energy from P 4/nmm to P 4/mmm which was obtained via step by step and concerted models (FIG. S8) shows that the associated activation barrier is high in energy (FIG. S9): 1.40, 1.31, and 2.18 eV for CuS, CuSe, and CuTe, respectively. This should prevent the transition between P 4/mmm and P 4/nmm at moderate temperature (details are given in section S6 of the SI). It also means that, while close in energy, the P 4/mmm and P 4/nmm polymorphs may be synthesized from different precursors and/or experimental procedures. Conversely, the 2D CuTe exhibits a ground-state structure in the P mmn geometry. It is noteworthy that the ground state phase of the 3D solid-state CuTe (bulk phase) presents a layered structure where our proposed 2D-P mmn CuTe slabs are stacked in a similar AA stacking configuration. Our ab initio CSP search thus leads to the experimentally characterized CuTe layer topology, assuming only the CuTe chemical formula.

The enthalpies of formation of the twelve 2D CuX phases are given in Table I. Excluding CuS P mmn and CuTe P 4/mmm, the formation enthalpies are negative, indicating the thermodynamic stability of the structures. These computed ∆H f at PBE level fall into the calculated formation enthalpies range of experimentally characterized known 2D materials 59 . This result allows us to invite the experimental materials community to investigate our predicted 2D CuX compounds. In addition, the formation enthalpies (∆H f ) of dynamically stable 2D CuX phases and reported 2D materials were calculated as summarized in Table S3 of SI. Notably, the synthesized P 6m2 CuSe is predicted in our USPEX searches, with formation enthalpy of 0.160 eV/atom, which is relatively higher than the CuSe phases reported in our work. Moreover, considering supported 2D compounds -CuSe (∆H f = 0.160 eV/atom), Cu 2 Si (∆H f = 0.413 eV/atom), and SiC (∆H f = 0.293 eV/atom) have been synthesized and characterized, we hope our theoretical findings will stress the interest of experimental colleagues and led to the experimental realization of 2D CuX on proper surfaces. Beside the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA), the SCAN meta-GGA and the hybrid HSE06 (HeydScuseriaErnzerhof) functionals have been employed to evaluate the relative stability of the 2D CuX structures. The 2D structures were fully optimized at the SCAN level of theory while single point energy HSE06 calculations are undertaken at PBE geometry (HSE06//PBE level). Table II compiled the results obtained from each functional. In order to simplify the comparison of values, the total energies are reported relative to the P 3m1 structure of each copper chalcogenide. a All the energies were obtained by doing a full geometry optimization, except for the HSE06 functional for which a single point calculation is performed on the optimized PBE structure, named HSE06//PBE.

For the 2D CuS structures, the stability order is not changed by the different functionals. For 2D copper selenide phases, the energy difference between the two lowest structures tends to be null when SCAN and HSE06/PBE levels are considered; P 4/mmm 2D CuSe becomes almost isoenergetic to P 3m1 2D CuSe. Therefore, one might conclude that both P 4/mmm and P 3m1 2D CuSe structures are potential synthesizable phases. P mmn CuTe phase is the lowest energetic structure at all levels of theory, but we note that P 4/nmm is located at only 2-5 meV/atom above P mmn at SCAN and HSE06 levels of theory. Finally, the order for the 3rd and 4th least stable 2D structures is not changed regardless of our level of theory.

The thermal stability of the structures was investigated via AIMD simulations at different temperatures. The calculated temperatures at which a structure remains stable, i.e. does not decompose, are summarized in Table III. Snapshots of the structures in the end of the 10 ps AIMD simulations are shown in figures S3-S5 of SI. According to AIMD results, the most thermally stable structure for 2D CuS has the P 3m1 geometry. This structure is found stable up to 800 K. The P 4/nmm and P 4/mmm structures are thermally stable up to 600 and 400 K, respectively. The simulations of P mmn showed that the structure is not stable even at 300 K. The results of AIMD simulations for 2D CuSe show that the structures with P 3m1 and P 4/mmm symmetry are the most stable structures with a stability temperature up to 600 K. The other two phases are not so stable: P mmn is found stable at 300 K only and P 4/nmm is not stable at this temperature.

Finally, for 2D CuTe structures, the P 3m1 geometry is found to be the most stable, with a stability temperature up to 800 K. Interestingly, the most energetically favorable geometry, P mmn, is not the most thermally stable one, which means that the activation barrier for decomposition is lower for P mmn. This result shows the importance of kinetic criteria in the prediction of viable compounds.

C. Structural, bonding, and electronic properties of Cu2X2 slabs

As discussed previously, 2D CuX compounds present four low-lying structure types. Different views of their associated Cu 2 X 2 slabs are depicted in figures 1 and 3.

In the following, we provide a more detailed structural description for each geometry, and discuss their chemical bonding mode. Both Manz' bond order (BO) and integrated crystal overlap Hamilton population (ICOHP) are computed. We found that both analysis give the same conclusion. Thus, only the ICOHP descriptor of bond strength will be discussed thereafter. It is presented in Table IV (see FIG. S10 of SI for BO results).

The first structure shown in Figure 3-(a) crystallizes in a trigonal system with the P 3m1 space group (No. 164, Z=2). This slab is constructed from two parallel layers of X chalcogen atoms arranged in monolayer triangular lattices. The Cu atoms similarly resemble two parallel triangular monolayers, sandwiched between the two X monolayers. In other words, the combination of each X monolayer with its adjacent Cu layer makes a buckled honeycomb Cu 2 X 2 lattice. Each Cu atom is then bonded to four equivalent X atoms to form a mixture of edge and corner-sharing CuX 4 tetrahedra (See Fig. 1-(a)). The interatomic distances between Cu-X and Cu-Cu atoms are also given in Table IV. The average bond length of Cu-X for this geometry is 2.31, 2.46, and 2.62 Å in CuS, CuSe, and CuTe, respectively. The size of the tetrahedra, therefore, increases with the atomic number. Conversely, Cu-Cu separation decreases significantly when selenium or tellurium atoms substitute sulfur ones in P 3m1 slab, a factor that can explain the relative energetic destabilisation of this phase in going from S to Se to Te in respect to the three other dynamically stable structures.

The second structure P 4/mmm Cu 2 X 2 slab, shown in Fig. 1-(b) and Fig. 3-(b), is formed by a square lattice of copper atoms with capping chalcogen atoms sitting in the apical positions of square bipyramids. Here, each Cu atom is then bonded in a square coplanar geometry to four X atoms. The Cu-X bond lengths are 2.34, 2.46, and 2.63 Å for CuS, CuSe, and CuTe, respectively, increasing with the atomic number. Contrary to the P 3m1 geometry, the same behavior is also observed for d Cu-Cu , an increase from 2.56 Å to 2.76 Å in going from S to Te (see Table IV). In 2D P 4/mmm CuTe, a short Te-Te separation of 3.51 Å is observed.

The third structure P 4/nmm also possesses a square lattice of Cu atoms. As can be seen in the Fig. 3-(c), the P 4/nmm phase (No. 129, Z=2) can be regarded as a structural modification of the P 4/mmm slab. From P 4/mmm to P 4/nmm, one chalcogen atom is migrating from its four-fold position to be located under a chalcogen atom making the square bipyramidal Cu 4 X 2 unit. Thus, P 4/nmm Cu 2 X 2 slab is then constructed from three nonequivalent square monolayers with a AB stacking. Cu 4 X 2 square bipyramids share their four basal Cu-Cu edges in this extended 2D Cu 2 X 2 net. The Cu-X bond lengths are similar to the ones found in P 4/mmm, showing the weak effect of the chalcogen position on the CuX bond lengths. Conversely, Cu-Cu distances decrease in Cu 2 X 2 slabs from S (2.74 Å) to Se (2.62 Å) then remain constant from Se to Te (2.62 Å). This structural behavior may explain the energetic stabilisation of P 4/nmm when going from S to Te. Finally, our fourth CuX candidate crystallizes in a rectangular 2D system and has the orthorhombic Pmmn space group (No. 59, Z = 2). Its topology resembles that of the P 4/nmm structure with rectangular Cu 4 units rather that square ones (see Fig. 1-(d)). There are two short and two long Cu-X bond lengths, and similarly, short and long X-X distances. This Cu 2 X 2 slab is encountered of distorted tetrahedral CuX 4 units is encountered as the lowest energy ground-state structure in 2D CuTe compound. This net allows short Te-Te separations of 3.26 Å and significant Te-Te bonding is computed (ICOHP(Te-Te)=0.89).

A question arises: what electronic factors govern the energy ranking of these four 2D geometries? To approach an answer to this question, we need to understand the structural details of these 2D nets. To do so, we have undertaken orbital interaction analysis with the assistance of total and projected DOS, ICOHP, and BO calculations. We refer to the seminal work of Zheng and Hoffmann 60 nomenon encountered in the four-electron two-orbital repulsion. In the extended CuSe net, the Fermi level crosses a high density of states, as expected for a d 9 metal. From the DOS and bandstructures calculations analysis at the PBE and PBE+U levels, all 2D compounds are metallic (see FIG. S6 of SI). So, all the structures are still metal and their band structures present globally the same trends. Let us focus now on P 3m1 and P 4/nmm structures, as representative nets with 3-and 4-coordinate copper atoms, respectively. One may see from Table IV that P 3m1 Cu 2 X 2 slab presents among the four proposed 2D structures with the strongest copper-chalcogen Cu-X ICOHP. This result validates the Zheng and Hoffmann orbital interaction model 60 : in a less dispersive metal sublattice, i.e. in P 3m1 3-coordinate Cu lattice, the metals interact better with the ligand sublattice. This criterion cannot however explain alone the relative stability of the four predicted 2D CuX nets; other subtle electronic effects must be involved.

The P 3m1 structure has a corrugated Cu net where each metal is coordinated to 3 others, while the bare metal layer in P 4/nmm (as well as P 4/mmm) is then a square lattice with four coordinated metal atoms. Thus, the copper d band dispersion is smaller in 3-coordinated Cu lattice P 3m1 than in 4-coordinated one P 4/nmm, as schematically illustrated in Fig. 5. This band is also Cu-X antibonding (more at the top of the band). Bonding and antibonding Cu-Cu orbital overlaps are weaker at the bottom and the top of the d band in P 3m1. Therefore, for a high band filling, d 9 Cu 2+ in CuX, 3-coordinated Cu P 3m1 should be then energetically favored while 4coordinated Cu structures are more destabilized. This is what we observed in 2D CuS. Sulfur is the most electronegative element in the chalcogen family, while Te is the less one (2.58, 2.55, and 2.10 on the Pauling scale for S, Se, and Te, respectively). Thus, in 2D CuS, one may safely assign a formal +2 charge for Cu, i.e. d 9 configuration, and the octet is obeyed for S 2-. As calculated at all levels of theory, P 3m1 Cu 2 S 2 slab is by far the most stable structure while the other structures with four-coordinated metals are located higher in energy. By replacing S by Se, the Cu-X electronegativity difference is lowered, thus the Cu-X charge transfer decreases. One may expect a formal d n configuration of copper higher than d 9 . When the d band is completely filled, 3-coordinated Cu and 4-coordinated Cu slabs must compete in energy. This crude model explains well why P 3m1 and P 4/mmm are closer in energy in 2D CuSe. Now, why is P 4/mmm with Cu 4 Se 2 square bipyramidal units lower in energy that P 4/nmm with Cu 4 Se square pyra-midal units? In P 4/mmm, the computed Se-Se ICOHP is larger in P 4/mmm than in P 4/nmm (0.52 vs 0.19, respectively) because of the through-space Se-Se stabilizing interaction present in Cu 4 Se 2 units. A short interlayer Se-Se separation of 3.21 Å is observed in P 4/mmm Cu 2 Se 2 slab, reflecting the delocalization of electrons over the whole cluster.

Finally, Te is the less electronegative element of this chalcogen trio, with the biggest element radius. This chemical composition adopts the P mmn geometry as its lowest energy structure at PBE level. A corrugated copper monolayer is observed in response to the packing of Te atoms. Each copper atom seats in an elongated Te 4 tetrahedron with small and large Te-Cu-Te bond angles (73.4 • and 101.4 • ). This structural arrangement leads to short (3.26 Å) and long (4.05 Å) in-plane Te-Te bond lengths. Thus, a significant Te-Te bonding is calculated along each 1D tellurium chain with a computed ICOHP(Te-Te) of 0.89 eV/pair (see Table IV). This bonding mode stabilizes this structure type, as in P 4/mmm which is almost isoenergetic to P mmn at SCAN and HSE06 levels of theory. In both P mmn and P 4/mmm structures, some delocalized Te-Te covalent bonding is observed (see Table IV).

Our ELF analysis shows that electron density is delocalized along the X-X chains in each X monolayer of P mmn CuX (see Fig 6). Meanwhile a strongly localized lone pair is observed on each chalcogen atom in the other structures, specifically on electronegative sulfur atoms. This is expected for a 6 valence electron atom in such coordination mode.

IV. CONCLUSION

This work presents CSP search for 2D copper chalcogenides CuX compounds. The extensive search was carried out by using evolutionary algorithm and ab initio calculations based on DFT. According to our results, four geometries with space groups P 3m1, P 4/mmm, P 4/nmm, and P mmn are predicted for 2D CuX (X=S, Se, and Te). And all the 2D CuX phases reported in this work are nonmagnetic. These geometries were confirmed to be dynamically stable. The P 3m1 geometry has the lowest formation energy for CuS and CuSe 2D materials. Conversely, the P mmn geometry is a global minimum for the 2D CuTe structure. Moreover, the thermal stability of the structures was investigated by AIMD simulations at different temperatures. The P 3m1 structure, which is determined by AIMD studies to be stable up to 800, 600, and 800 K for CuS, CuSe, and CuT e, respectively, is found to be the most thermally stable structure. For CuSe, a new structure with the P 4/mmm symmetry is also equivalently stable at the same temper- enlightens the chemical bonding in the 2D nets. The analysis of ICOHP gives insights into the structural stability of 2D CuX nets. All 2D CuX 2D structures have a metallic character with a iono-covalent nature of the Cu-X bonds. We expect that our findings will inspire the materials science community to synthesize our in silico proposed two-dimensional CuX phases.

V. ACKNOWLEDGMENTS

FIG. 1 .

 1 FIG. 1. (color online) Views of the four single Cu2X2 sheets of 2D CuX phases, (a) P 3m1, (b) P 4/mmm, (c) P 4/nmm and (d) P mmn. Their copper-chalcogen coordination polyhedra are enlightened. Cu and X (X= S, Se, Te) atoms are shown in blue and yellow, respectively.

FIG. 2 .

 2 FIG. 2. (color online)Total energies (in meV/atom, PBE level) of the four lowest structures of CuS, CuSe, and CuTe compounds, relative to the P 3m1 phase. (The dotted lines between each chalcogenide structure with the same geometry is just a guide for the eye.).

FIG. 3 .

 3 FIG. 3. (color online) Extended top and side views of the (a) P 3m1, (b) P 4/mmm, (c) P 4/nmm, and (d) P mmn structures. The primitive cell of each structure is shown in solid line. Cu and chalcogen (X= S, Se, Te) atoms are shown in blue and yellow, respectively. The c axis is perpendicular to the Cu2X2 slab.

FIG. 4 .

 4 FIG. 4. (color online) (a) Construction of the two-dimensional bands in 2D CuSe; (b) Calculated total and projected density of states for 2D CuSe P 3m1.

FIG. 5 .

 5 FIG. 5.Copper d band dispersion in corrugated and planer Cu net with 3-coordinated Cu lattice P 3m1 and 4coordinated one P 4/nmm, respectively.

  FIG. 6. (color online) 3D and 2D electron localization function (ELF) of CuX compounds in P 3m1, P 4/mmm, P 4/nmm, and P mmn structures. Both 3D and 2D plots are given for an ELF value of 0.8.

TABLE I .

 I Formation enthalpies (in meV/atom) of the predicted 2D CuS, CuSe, and CuTe structures at the PBE level of theory. Boldface entries emphasize the lowest energy phases.

	Structure	2D-CuS	2D-CuSe	2D-CuTe
	(space group)
	P 3m1	-82	-75	-46
	P 4/mmm	-28	-67	+17
	P 4/nmm	-29	-59	-50
	P mmn	+56	-43	-69

TABLE II .

 II Total energy per atom relative to the P 3m1 structure (in meV/atom) of the 2D CuX structures at PBE, SCAN, and HSE06 levels of theory. a Boldface entries emphasize the lowest energy structures.

	compound	structure	PBE	SCAN	HSE06//PBE
		P 3m1	0	0	0
	CuS	P 4/mmm	54	41	39
		P 4/nmm	53	93	76
		P mmn	138	98	160
		P 3m1	0	0	0
	CuSe	P 4/mmm	8	-1	-2
		P 4/nmm	16	52	40
		mmn	32	72	64
		P 3m1	0	0	0
	CuTe	P 4/mmm	63	94	63
		P 4/nmm	-4	-8	-9
		P mmn	-23	-13	-11

TABLE III .

 III Calculated Temperatures (in K) at which the predicted 2D CuX (X= S, Se, Te) structures were found thermally stable in the temperature range of 300-1000 K. Boldface entries emphasize the lowest energy phases.

		CuS	CuSe	CuTe
	P 3m1	800	600	800
	P 4/mmm	400	600	-a
	P 4/nmm	600	-a	300
	P mmn	-a	300	600

a At 300K, our AIMD simulation leads to the breaking of Cu-X bonds.

TABLE IV .

 IV Bond length l (in Å) and integrated crystal Hamilton population (ICOHP in eV/pair) in CuX (X= S,Se,Te) 2D structures. Cu'(X') denotes a second type of Cu (X) atom considered for BO calculation.

			CuS			CuSe		CuTe
	structure	bond	l	ICOHP	l	ICOHP	l	ICOHP
	P 3m1	Cu-Cu'	2.87	0.10	2.61	0.25	2.61	0.26
		Cu-X	2.26	1.69	2.41	1.57	2.61	1.31
		Cu'-X	2.36	1.19	2.50	1.40	2.63	1.42
		X-X	3.63	0.19	3.76	0.18	3.90	0.22
	P 4/mmm	Cu-Cu'	2.56	0.21	2.64	0.22	2.76	0.18
		Cu-X	2.34	1.39	2.46	1.42	2.63	1.30
		X-X a	2.95	0.40	3.21	0.52	3.51	0.46
		X-X b	3.62	0.03	3.73	0.08	3.91	0.11
	P 4/nmm	Cu-Cu'	2.74	0.11	2.62	0.20	2.62	0.23
		Cu-X	2.30	1.49	2.45	1.48	2.64	1.24
		X-X	3.70	0.17	3.71	0.19	3.70	0.33
	P mmn	Cu-Cu'	2.58	0.19	2.57	0.26	2.66	0.24
		Cu-X	2.29	1.57	2.42	1.57	2.62	1.37
		Cu'-X	2.45	0.93	2.59	0.88	2.73	0.88
		X-X c	2.97	0.53	3.09	0.69	3.26	0.89
		X-X d	3.99	0.15	3.90	0.15	4.05	0.16

a between monolayers, b in-plane, c short, d long.

We thank the ANR PRCI Predict 2D Nanomat (P.I. G.F.), the FEDER and Nouvelle Aquitaine Region (France).

We also acknowledge the High-Performance Computing Centers of Poitiers University (Mésocentre SPIN, France) and the Irene/TGCC, and Jean Zay/IDRIS GENCI (France) under projects no. X2016087539 and A0100807539 for allocation of computing time.