Luc Lesoil

Mathieu Acher

Arnaud Blouin

Deep Software Variability: Towards Handling Cross-Layer Configuration

Configuring software is a powerful means to reach functional and performance goals of a system. However, many layers (hardware, operating system, input data, etc.), themselves subject to variability, can alter performances of software configurations. For instance, configurations' options of the x264 video encoder may have very different effects on x264's encoding time when used with different input videos, depending on the hardware on which it is executed. In this vision paper, we coin the term deep software variability to refer to the interaction of all external layers modifying the behavior or non-functional properties of a software. Deep software variability challenges practitioners and researchers: the combinatorial explosion of possible executing environments complicates the understanding, the configuration, the maintenance, the debug, and the test of configurable systems. There are also opportunities: harnessing all variability layers (and not only the software layer) can lead to more efficient systems and configuration knowledge that truly generalizes to any usage and context.

INTRODUCTION

Software systems can often be configured to reach specific functional and performance goals, either statically at compile time or through the choice of command line options at runtime. Beyond having a strong impact on non-functional properties [START_REF] Han | An Empirical Study on Performance Bugs for Highly Configurable Software Systems[END_REF], this configuration stage may introduce bugs [START_REF] Abal | 42 Variability Bugs in the Linux Kernel: A Qualitative Analysis[END_REF][START_REF] Attariyan | Using Causality to Diagnose Configuration Bugs[END_REF][START_REF] Zhang | Which Configuration Option Should I Change?[END_REF], or even prevent the software from starting altogether [START_REF] Yin | An Empirical Study on Configuration Errors in Commercial and Open Source Systems[END_REF]. Several researchers have been investigating the use of testing or statistical learning to determine whether a specific configuration would raise an error [START_REF] Thüm | A Classification and Survey of Analysis Strategies for Software Product Lines[END_REF][START_REF] Zhang | Automated diagnosis of software configuration errors[END_REF], or to predict performances [START_REF] Pereira | Learning Software Configuration Spaces: A Systematic Literature Review[END_REF][START_REF] Siegmund | Performance-Influence Models for Highly Configurable Systems[END_REF], typically in a controlled environment. However, it has been shown that the software environment (in a broad sense, see Figure 1) can interact with its configuration, thus changing its validity [START_REF] Cohen | Interaction Testing of Highly-Configurable Systems in the Presence of Constraints[END_REF] or its performances [START_REF] Fichte | A Time Leap Challenge for SAT-Solving[END_REF][START_REF] Jamshidi | Learning to sample: exploiting similarities across environments to learn performance models for configurable systems[END_REF]. That is, only considering the software layer might hide important bugs or provide non-optimal values for configurations' options of the software. Environment properties (hardware, operating system, etc.) can be organized into variability layers, as those depicted on Figure 1. We call deep software variability the interaction of all variability layers that could modify the behavior or non-functional properties of a software.

We illustrate this concept of deep software variability with concrete examples in Section 2. We emphasize the challenges and opportunities of configuring a software w.r.t. this deep software variability in Section 3. We highlight existing related work in Section 4.

Let us illustrate the concept of deep software variability with x264, a popular, highly-configurable video encoder. In Figure 2, we show the duration of the video compression (in seconds) and the size of the encoded video (in MB), for two configurations 1 under two different environments 2 . Performances were measured with two different input videos 3 and the same version of x264 4 . The remainder of this section describes five challenges related to x264's deep variability.

Impacts of variability layers

According to Figure 2, changing environment A for environment B would make the compression approximately 15 times longer. However, the compression durations of configurations 1 and 2 remain stable, and roughly proportional between A and B for both videos. There are more complicated cases, where the variability layers change the relative importance of configuration options. For an example, switching from configuration 1 to configuration 2 (i.e., deactivating the feature mbtree) decreases the size of the encoded animation video (33MB↘ 21MB), while it increases the size of the encoded vertical video (28MB↗ 34MB). It becomes difficult to determine the effect that mbtree may have on the size of the output video. Similar results have been reported in the literature [START_REF] Pereira | Sampling Effect on Performance Prediction of Configurable Systems: A Case Study[END_REF][START_REF] Jamshidi | Learning to sample: exploiting similarities across environments to learn performance models for configurable systems[END_REF][START_REF] Valov | Empirical comparison of regression methods for variability-aware performance prediction[END_REF]: hardware, input data, software version (or a combination thereof) can impact performance properties of x264.

Deep testing and benchmarking

To go further than the example of Figure 2, let us consider a general x264's performance model that has to work in environments A and B, but for 10 input videos and 20 boolean features of x264; assuming that each video is compressed in ten seconds, testing exhaustively the 21 millions of possible configurations would take about 8 months of computation. Regardless of the time needed to compute these measurements, it would require both a lot of human and computational resources to test multiple environments and a framework allowing to measure performances on an higher level than the software layer (e.g., automatically testing multiple workloads and input videos).

Improve documentation

Occasionally, requirements between layers can break the build of x264. For instance, manually compiling x264(version 0.152.2854, git commit e9a5903) on Ubuntu 18.04 LTS requires the upgrade of the Linux package nasm (version > 2.13) 5 , unless one deactivates asm during the compilation (i.e., specifying the configuration option -disable-asm). To the best of our knowledge, this information is neither centralized nor included in the official documentation ; we have to search in forums, where developers gave the answer to this problem. Add: (1) software dependencies (2) the effects of a variability layer on software performances would certainly improve x264's documentation.

Generalize the configuration knowledge

Between environments A and B, both the hardware and its operating system are changed. What would be x264's performances on the environment C, composed of A's hardware and B's operating system? Can we estimate them without measuring configurations in environment C (i.e., just with measurements of environments A and B)? As long as we are unable to quantify the marginal effect of each layer, changing a single property of the environment makes us incapable of predicting x264's performance.

Cross-layer tuning

By generalizing the measures of Figure 2 on other environments (cf. Section 2.2), a streaming website could then automatically choose an environment, and configure a server fully dedicated to the encoding of the animation video. In order to achieve this specialization, experts would have to fix the layers' variability, selecting the optimal configuration corresponding to this setting, and hopefully increasing the overall performances of x264.

CHALLENGES AND OPPORTUNITIES

We propose to break the deep software configuration challenge down to five steps:

(1) Identify influential variability layers that impact software's non-functional properties (2) Test software's performances in multiple environments (3) Improve knowledge about how configurable systems interact with their environments (4) Transfer performances across environments and build performance models robust to environment changes (5) Specialize the environment for the software, i.e., configure each layer to improve software performances

For each step, we describe the current problems encountered by users or developers, the challenges that researchers could tackle, and the potential benefits of facing them.

Impacts of variability layers

We identify at least four potential variability layers:

• Hardware: material part of the computer on which the software is executed. Examples of properties: the computer or server model and brand, the amount and type of Random-Access Memory available or the Graphics Processing Unit frequency and memory, etc. • Operating system: the operating system and dependencies that could be used by the software. Examples of properties : the distribution of an operating system, its version, the linux kernel version and modules, the packages that are installed on the operating system, etc. • Input data: the raw data fed and processed by the software. Examples of properties: a video for x264, a word or a sentence for an online translator, a C program for a compiler, etc. • Software: how the software was configured and built.

Examples of properties: any configuration option (e.g., -no-cabac for x264), how the software was compiled (compiler and version), the version of the software, etc.

Figure 1 shows an example of these layers and possible impacts on two performances (encoding time and size).

3.1.1 Problem. Software practitioners may not be aware of the impacts of some layers over some performance properties. In fact, it is hard to know a priori which and to what extent configurations' layers do have an impact on software. For software variability researchers e.g., seeking to train predictive performance models, these variability layers can disrupt the model's training and weaken the results, preventing their generalization and making them useless to operate under different conditions.

3.1.2 Challenges. Can we predict how much (combination of) variability layers will change non-functional properties of a software? Owing to the cost of executing cross-layer configurations', identifying the influential layers remains a challenge; their importances might vary depending on the subject system, the targeted non-functional property as well as the software configurations considered. While characterizing the impact of variability layers individually is necessary, it is not sufficient. As observed for software options, we conjecture that cross-layer configurations interact and can alter software performances and configuration option's effects.

3.1.3 Opportunities. Once their effects have been identified and quantified, the influential variability factors can be leveraged to improve software performance, while the others remain fixed and can be forgotten (e.g., when benchmarking or transfer learning, see hereafter).

Deep testing and benchmarking

3.2.1 Problem. Since each user applies the software on his input data, with his dependencies on his laptop, there are almost as many possible environments as users. Because of the combinatorial explosion of possible environments, it becomes time-consuming to test all possible environment performances (cf. Section 2.2). Moreover, testing all variability layers requires to test variability on the system level, which is not supported by existing tools.

Challenges.

How to sample a good benchmark of environments? A good benchmark should be representative of the real usage of the software. It has to contain not only a wide variety of possibilities within a layer (e.g., different brands or models for the hardware layer), but also crosslayer diversity (e.g., several variants of software running in different versions of operating systems). These two rules should be followed, in order to avoid missing corner case environments, in which combination of layers interact with software features, causing performance bugs.

Can we build a framework that tests configurable systems in multiple environments? In general, testing software configurations requires to automatically observe their properties, being functional or related to performance. Deep software variability further challenges the automation: (1) each layer (e.g., operating system) may come with its own specific tooling; (2) integrating all layers together to test the software is not straightforward: for instance, instrumenting the tests of a software configuration in variants of hardware and operating systems. Another open challenge is: Can we detect corner cases configurations of layers that eventually cause software's performance bugs? As exemplified in Figure 2, a change in the input data or hardware can cause performance differences. There are two hypotheses worth investigating.

First, such differences are expected and can be explained. For example, it is normal that a given software option has more or less effect when processing over a specific input. A second hypothesis is that performance differences among variability layers are in fact a manifestation of a software bug that causes an accidental and unexpected decrease in performance. For validating or refuting the "bug or feature" hypothesis, software developers and domain experts of the different layers should be involved e.g., by reviewing performance results or formalizing expected impacts of layers.

3.2.3

Opportunities. Some software bugs, especially related to performance, never manifest in a fixed environment. The diversification of the different layers is an opportunity to test the robustness and resilience of the software layer in multiple environments. Another observation is that an absolute number for e.g., execution time has little meaning in itself and in a fixed environment. On the contrary, some performance issues are noticeable under the conditions the distributions are put in perspective and compared to others. That is, another opportunity to detect unexpected performance differences relatively to different configuration layers. Overall, developers and operationals can exploit deep software variability to detect more (performance) bugs.

Improve documentation

Challenges. Can developers include deep software variability constraints and effects in the documentation?

The challenge is to develop user-friendly documentation that alerts users to the effects of deep software variability, its requirements, and gives them an overview of how the software should be configured for their environment. This documentation could be used as a pretext to centralize measurements from volunteers, benefiting to the community members ; users could compare their usage of the software to others running it in similar environments.

3.3.3

Opportunities. Building such a documentation would encourage the fine-grained, informed customization of configurable systems. When the effects of software options are captured, users may try to optimize and conveniently replace the default software configuration. Advanced users could even build their own variant of the software, working in their environment and fulfilling their own needs. Furthermore, different experts can participate to the customization process, each focusing on their domain of expertise (e.g., hardware, operating system) with the help of a precise documentation.

Generalize the configuration knowledge

3.4.1 Problem. Changing a simple property of its executing environment can alter the performances of a software. For instance, operating systems evolve frequently, introducing numerous changes in software environments; performance models has to evolve and adapt their prediction following the same rule. Unfortunately, practitioners cannot afford to train one new model per environment.

Challenges. Is it possible to accurately predict software performance for any environment?

The goal is to train a general model that aims at predicting software performances for any user's environment (i.e., robust to major changes, like changing simultaneously the hardware and the input data). Another achievement is to identify the limits of the model, and include the tests (cf.Section 3.2) as a part of the learning process when transferring performances leads to poor predictions.

Opportunities.

Being able to transfer performance from one environment to another avoids the need to build a performance model for each environment, thus saving time and energy resources. From the user's point of view, it allows to estimate the performance of a configuration without testing it, facilitating the configuration of the software for its environment.

Cross-layer tuning

3.5.1 Problem. Some operating systems and packages include options that may be over-adapted to the environment in which they run. In a sense, these layers are too generalpurpose. It can lead to performance bugs when the configuration options misidentified the software environment. In the same vein, a software always configured with the default values will not benefit from the optimization coming from its environment (e.g., parallelization if the default sticks to sequential tasks), wasting energy and computing time.

3.5.2 Challenges. How can we tune the software for one specific executing environment? Instead of enduring deep software variability, the goal is to pre-select the right environment for the software, tuning each layer separately in such a way it improves the overall software performances. The tuning process may have a cost since it requires to know which variability layer has an effect on software, what will be its interactions with other layers, and how to configure it in symbiosis with the rest of the environment. This cost should be confronted to the underlying benefits.

Opportunities.

Since cross-layer tuning works on multiple variability layers, in comparison to a simple software tuning (i.e., mono-layer optimization), we expect it to largely outperform default values configurations in terms of performances (e.g., energy consumption, execution time, etc.) while keeping the same level of service. Cross-layer tuning should be more efficient at accomplishing a unique task; however, it requires to restrict some flexibility (e.g., at runtime) and forget some variability of each layer.

RELATED WORK

For each step defined in Section 3, we refer to existing research papers that have identified and addressed parts of the deep software variability.

Impacts of variability layers

Numerous research papers study the effects of variability layers on software's performances or on configuration options: hardware [START_REF] Fichte | A Time Leap Challenge for SAT-Solving[END_REF][START_REF] Teng | Reliability Modeling of Hardware and Software Interactions, and Its Applications[END_REF][START_REF] Zhang | Maximizing Performance Under a Power Cap: A Comparison of Hardware, Software, and Hybrid Techniques[END_REF], workloads [START_REF] Kuenzer | Unikernels Everywhere: The Case for Elastic CDNs[END_REF][START_REF] Madhavapeddy | Unikernels: Library Operating Systems for the Cloud[END_REF][START_REF] Olivier | A Binary-Compatible Unikernel[END_REF][START_REF] Sayagh | On Cross-Stack Configuration Errors[END_REF][START_REF] Zhang | Maximizing Performance Under a Power Cap: A Comparison of Hardware, Software, and Hybrid Techniques[END_REF], variants [START_REF] Muehlbauer | Identifying Software Performance Changes Across Variants and Versions[END_REF][START_REF] Thüm | Performance Analysis Strategies for Software Variants and Versions[END_REF], versions [START_REF] Hasreiter | Evolution of Performance Influences in Configurable Systems[END_REF][START_REF] Sandoval Alcocer | Performance evolution blueprint: Understanding the impact of software evolution on performance[END_REF][START_REF] Thüm | Performance Analysis Strategies for Software Variants and Versions[END_REF][START_REF] Werner | Energy and Performance Evolution of Configurable Systems: Case Studies and Experiments[END_REF], compilation options [START_REF] Hoste | Cole: Compiler Optimization Level Exploration[END_REF][START_REF] Plotnikov | Automatic Tuning of Compiler Optimizations and Analysis of their Impact[END_REF] and input data [START_REF] Pereira | Sampling Effect on Performance Prediction of Configurable Systems: A Case Study[END_REF][START_REF] Ding | Autotuning Algorithmic Choice for Input Sensitivity[END_REF]. Such studies provide evidence that some layers have a noticeable impact on the software (configuration) layer. Deep software variability provides a conceptual, unified framework to systematically investigate the significance of the problem. To the best of our knowledge, there are two open research directions: (1) a systematic and comprehensive assessment of the influence of individual layers onto the distribution of software configurations: only a few hardware, input data, or software configurations are usually considered; (2) empirical knowledge about how layers composed and interact each other, for example how hardware and input data both influence the performance properties of a configurable system (e.g., x264).

Deep testing and benchmarking

Related work proposes several methods to test the software layer, typically to improve coverage of the software configuration space: combinatorial interaction testing [START_REF] Cohen | Interaction Testing of Highly-Configurable Systems in the Presence of Constraints[END_REF], twise [START_REF] Perrouin | Automated and Scalable T-wise Test Case Generation Strategies for Software Product Lines[END_REF][START_REF] Thüm | A Classification and Survey of Analysis Strategies for Software Product Lines[END_REF], genetic algorithms [START_REF] Henard | Multi-Objective Test Generation for Software Product Lines[END_REF], adversarial machine learning [START_REF] Temple | Towards Adversarial Configurations for Software Product Lines[END_REF], to name a few. These approaches are usually applied at the software layer, but need to be adapted at the system level (e.g., operating system) and integrated into a comprehensive testing or benchmarking pipeline. The effectiveness of existing testing methods (e.g., differential testing [START_REF] William M Mckeeman | Differential testing for software[END_REF], metamorphic testing [START_REF] Segura | A Survey on Metamorphic Testing[END_REF], multimorphic testing [START_REF] Temple | Empirical Assessment of Multimorphic Testing[END_REF]) in the context of deep software variability is another research direction.

Improve documentation

State of the art provides evidences that documentation sometimes lacks of details [START_REF] Aghajani | Software Documentation Issues Unveiled[END_REF][START_REF] Forward | The Relevance of Software Documentation, Tools and Technologies: A Survey[END_REF], and give insights on how to build a configurable system's documentation [START_REF] Bakman | Automatic documentation of configurable systems by outputting explanatory information of configuration parameters in a narrative format and configuration parameters differences[END_REF][START_REF] Bakman | Method and system for automatic documentation of configurable systems[END_REF][START_REF] Robillard | On-demand Developer Documentation[END_REF]. Deep software variability further challenges developers of software configurable systems: the configuration documentation should mention the potential effects of other layers on the tuning of configuration options.

Generalize the configuration knowledge

Several research papers have shown that it is possible to transfer software performances from a source environment to a target environment, using either the similarities between environments [START_REF] Dorn | Generating Attributed Variability Models for Transfer Learning[END_REF][START_REF] Jamshidi | Learning to sample: exploiting similarities across environments to learn performance models for configurable systems[END_REF][START_REF] Valov | Transferring Performance Prediction Models Across Different Hardware Platforms[END_REF], or selecting the best source environment [START_REF] Nair | Transfer Learning with Bellwethers to find Good Configurations[END_REF]. Including causalities in the performance model to improve the transfer [START_REF] Ali Javidian | Transfer Learning for Performance Modeling of Configurable Systems: A Causal Analysis[END_REF][START_REF] Krishna | CADET: A Systematic Method For Debugging Misconfigurations using Counterfactual Reasoning[END_REF] is another promising research direction to explore. Deep software variability provides a concrete case for transfer learning: a change in a specific layer (e.g., hardware) is a new target environment and calls to adapt the prediction model obtained out of a source (i.e., original environment).

Cross-layer tuning

Various research papers [START_REF] Hoste | Cole: Compiler Optimization Level Exploration[END_REF][START_REF] Kuenzer | Unikernels Everywhere: The Case for Elastic CDNs[END_REF][START_REF] Madhavapeddy | Unikernels: Library Operating Systems for the Cloud[END_REF][START_REF] Olivier | A Binary-Compatible Unikernel[END_REF][START_REF] Plotnikov | Automatic Tuning of Compiler Optimizations and Analysis of their Impact[END_REF][START_REF] Zhang | Maximizing Performance Under a Power Cap: A Comparison of Hardware, Software, and Hybrid Techniques[END_REF] configure their variability layer to optimize non-functional properties of software, while others, as [START_REF] Ding | Autotuning Algorithmic Choice for Input Sensitivity[END_REF][START_REF] Sayagh | On Cross-Stack Configuration Errors[END_REF], adapt configuration options of software to variability layers. Such works are in line with the vision exposed in Section 3.5 though only one layer at a time is usually considered as part of the specialization process.

CONCLUSION

In this vision paper, we proposed deep software variability and described its many challenges and opportunities. As illustrated with x264 and partly observed in the literature, many layers (hardware, operating system, input data, etc.), themselves subject to variability, can alter the configurable software layer. Deep software variability calls to investigate how to systematically handle cross-layer configuration. At this step of the research, we ignore to what extent deep software variability impacts performances of software. May configurable systems be more or less sensitive to deep variability, depending on quantitative properties of interest (e.g., energy consumption, execution time)?

In fact, deep software variability concerns all scientists that aim at publishing reproducible research works; if your results are not consistent from one environment to another, your conclusions could be either weakened or invalidated. Many scientific domains are potentially impacted. Applied to ethical decisions, deep software variability could also have dramatic consequences; what would you do if a machine learning model deciding on the release of a prisoner changes its prediction when modifying its executing environment?

To make the vision of this paper more explicit, we recall the general idea of this paper: when considering software variability, we encourage you to broaden the scope of your research, and include parts of software environments in the study. Lastly, we illustrate how we could use the different steps detailed throughout the paper to concretely handle deep variability for a given software:

(1) First, we have to identify influential variability layers for this software. For instance, to estimate how the input data layer interacts with its configuration options, we could consider a large amount of inputs and measure their effects on performances. If the effects of features remain the same across all inputs, we could argue that input data is a negligible layer of variability, and could remain fixed when benchmarking the software. If not, this layer is influential, and a representative set of inputs should be extracted. This process should be applied to each layer. (2) Each influential layer of variability is now represented by elements of reference condensing its inner variability (e.g., for the hardware layer, at least a server, a laptop and an Arduino). Then, we propose to measure software performances in different environments, composed of different layer's elements (e.g., a laptop extracted of the hardware layer, an operating system coming from the operating system's layer, etc.). (3) Computing these measurements will allow us to detect some bugs and reference them in the documentation. (4) With this sample of measurements, we could build a performance model working for all environments.

A novel environment would be associated to a set of reference elements. With transfer learning and sampling techniques, we would be able to build a dedicated model designed for this new environment. (5) Being able to predict accurately the performance of each environment would allow us to improve the default configuration of software w.r.t. their environment, i.e., tune them for their environment.

Figure 1 :

 1 Figure 1: x264's deep variability

Figure 2 :

 2 Figure 2: Combinatorial explosion of cross-layer configurations for measuring x264's deep variability

1 1 :

 1 mbtree activated, 2 : mbtree deactivated 2 A : Ubuntu 20.04 LTS on a Dell Latitude 7400, B : Debian 10.4 on a Raspberri Pi 4 model B 3 animation for the 𝐴𝑛𝑖𝑚𝑎𝑡𝑖𝑜𝑛_1080𝑃 -3𝑑67.𝑚𝑘𝑣 video, vertical for the 𝑉 𝑒𝑟𝑡𝑖𝑐𝑎𝑙𝑉 𝑖𝑑𝑒𝑜_1080𝑃 -2195.𝑚𝑘𝑣 video, both extracted from the Youtube User General Content Dataset [40] 4 x264version 0.155.2917, git commit 0a84d98

3. 3 . 1

 31 Problem. Assisting users in charge of configuring software is a still an open issue per se; deep software variability exacerbates this problem. As described in Section 2.3, impacts of different layers are only partially reported in the documentation. Due to the lack of documentation on the relative effects of configuration options, practitioners have to ask directly to experts in forums or open topics to make their software work, or just to configure it properly.

We cross-checked the information provided in this webpage

This research was funded by the ANR-17-CE25-0010-01 Vary-Vary project.


Baltej Singh

Ziliang Wang

# Sunkyu Park

Sai Gopalakrishnan

Gautam

Jean- Noël Chotard

Laurence Croguennec

Dany Carlier

Anthony K Cheetham

Christian Masquelier
email: ^christian.masquelier@u-picardie.fr

Pieremanuele Canepa
email: pcanepa@nus.edu.sg

A Chemical Map of NaSiCON Electrode Materials for Sodium-ion Batteries

Na-ion batteries are promising devices for smart grids and electric vehicles due to cost effectiveness arising from the overall abundance of sodium (Na) and its even geographical distribution. Among other factors, the energy density of Na-ion batteries is limited by the positive electrode chemistry. NaSICON-based positive electrode materials are known for their wide range of electrochemical potentials, high ionic conductivity, and most importantly their structural and thermal stabilities. Using first-principles calculations, we chart the chemical space of 3d transition metal-based NaSICON phosphates of formula Na x MM'(PO 4 ) 3 (with M and M'= Ti, V, Cr, Mn, Fe, Co and Ni), to analyze their thermodynamic stabilities and the intercalation voltages for Na + ions. Specifically, we computed the Na insertion voltages and related properties of 28 distinct NaSICON compositions. We investigated the thermodynamic stability of Na-intercalation in previously unreported Na x Mn 2 (PO 4 ) 3 and Na x VCo(PO 4 ) 3 . The calculated quaternary phase diagrams of the Na-P-O-Co and Na-P-O-Ni chemical systems explain the origin of the suspected instability of Ni and Co-based NaSICON compositions. From our analysis, we are also able to rationalize anomalies in previously reported experimental data in this diverse and important chemical space.

Introduction

Developing new battery systems capable of storing increasing quantities of energy poses extraordinary scientific and economic challenges. Lithium (Li)-ion batteries supplies the world's portable devices, such as mobile phones, cameras, laptops, etc., but the supply chains of Li and the required transition metals may soon be at risk due to geopolitical considerations. [1], [2] Sodium (Na) appears as a viable alternative to Li for battery applications due to its earth abundance, as well as the possibility of harvesting it directly from sea water.

With Na being ~50 times more affordable compared to Li, Na-ion batteries (NIBs) are being explored by researchers worldwide, including a number of commercialization attempts. [3][4][5][6][7][START_REF] Mukherjee | Batter[END_REF] In addition, inexpensive stainless-steel current collectors are typically used in NIBs instead of the expensive copper ones found in Li-ion cells.

Cathode materials factor prominently in the overall energy density stored by NIBs, and the optimization of electrode chemistries to provide high intercalation voltages and gravimetric/volumetric capacities remains a crucial aspect in the design of competitive NIBs.

In theory, transition metal layered oxide-based cathode materials for NIBs would offer the largest theoretical energy densities. [5][6][7]9,10] However, the dominant two dimensional character of their structures affects strongly the longevity of the cathode material, with the electrochemical cell exhibiting a shorter-than-desired cycle life. [6,7,[10][11][12][13][14][15] Promising alternatives to layered oxides are polyanion-based cathode materials. [16] Polyanion positive electrode materials leverage a combination of multivalent cations (e.g., P 5+ , Si 4+ ), which are thermodynamically stable by virtue of the strong binding energy of their polyanionic covalent bonds. [17][18][19] An important class of phosphate electrodes discovered by Hong and Goodenough is the Natrium Super Ionic CONductors (NaSICONs), [20], [21] with formula Na x MM'(XO 4 ) 3 , where M and M' are metals and X = Si, P and/or S. NaSICON electrode materials and electrolytes typically display significant Na + -mobility. [22] The NaSICON framework is important for the development of new NIB materials, [17,[22][23][24][25][26][27][28][29] since its structural tunability enables the exploration of a vast chemical space, which greatly diversifies the potential, the electrochemical properties and related structural features. In theory, a charged MM'(PO 4 ) 3

NASICON framework with redox-active transition metals M and M' would allow the intercalation of up to 4 Na + ions, [9] which make these materials promising in terms of energy density. [30] Notwithstanding the challenges of working with carbon-based anodes in NIBS, [9,[31][32][33] one could benefit from the presence of 4 Na + ions in NaSICONs, e.g., Na 4 Mn II V III (PO 4 ) 3 . [34,35] However, the reversible extraction/intercalation of 4 Na + ions is still to be reported in NaSICON frameworks.

As an example of a NaSICON positive electrode, Na 3 V III V III (PO 4 ) 3 (NVP) can reversibly exchange two electrons (via the activation of the V IV /V III redox couple) delivering Na 1 V IV V IV (PO 4 ) 3 at an average voltage of ~3.4 V vs. Na/Na + with a capacity of ~110 mAh/g with an energy density of ~370 Wh/kg. [27,36,37] One additional Na + can be inserted (electro)chemically into NVP to achieve Na 4 V II V III (PO 4 ) 3 with a voltage of ~1.63 V vs. Na/Na + .

The reversible extraction of 3 Na + ions (i.e., starting from Na 4 V II V III (PO 4 ) 3 as the pristine material to Na 1 V IV V IV (PO 4 ) 3 ) is possible but not practical due to the 1.8 V difference between the two processes (V IV /V III and V III /V II ) and the difficult synthesis of Na 4 V 2 (PO 4 ) 3 . Uebou et al. [38] have demonstrated the extraction of 2.5 Na + ions with partial reversibility.

Symmetric Na electrochemical cells entirely made of NaSICON frameworks adopt NVP as both electrodes; [39] in the charged state of the battery the anode is Na 4 V III V II (PO 4 ) 3 and the cathode Na 1 V IV V IV (PO 4 ) 3 with Na 3 Zr 2 Si 2 PO 12 as the electrolyte. This could deliver an average voltage of ~1.8 V. [39][40][41] NVP suffers from poor intrinsic electronic conductivity due to the isolated VO 6 octahedra, which are never face-, edge-or corner-sharing with each other, [25,42] in contrast to their arrangement in layered-oxide cathodes where edge-sharing exists between MO 6 moieties. [6,9] The Ti analogue of NVP is the Ti IV -containing NaTi IV Ti IV (PO 4 ) 3 composition, which is stable in air and can reversibly exchange two electrons by benefitting from the Ti IV /Ti III redox couple, but at a lower voltage of ~2.1 V vs. Na/Na + compared to that of ~3.4 V for V IV /V III redox in NVP. [25,26] NaTi IV Ti IV (PO 4 ) 3 is also commonly used as a negative electrode material. [28] The NaSICON Na 3 Fe III Fe III (PO 4 ) 3 exhibits a voltage-composition plateau at ~2.5 V vs. Na/Na + with a discharge capacity of ~61 mAh/g and associated to the redox couple Fe III /Fe II . [43] The Fe IV /Fe III redox was claimed to be observed in Na 3 Fe III Fe III (PO 4 ) 3 at ~3.4V vs. Na/Na + , with a limited specific discharge capacity of ~22 mAh/g. [44] Nevertheless, the formation of tetravalent iron ions has not yet been reproduced.

Yamada et al. [45] reported a voltage of ~4.5 V vs. Na/Na + (enabled by the reversible Cr IV /Cr III redox couple) when Na + is extracted from Na 3 Cr III Cr III (PO 4 ) 3 , which corresponds to the highest voltage ever reported among NaSICON materials with a single TM species.

Importantly, the Mn, Co and Ni, Na x M 2 (PO 4 ) 3 NaSICON analogues have not yet been reported.

The high voltage delivered by the Cr IV /Cr III redox couple offers scope for exploring new highvoltage/capacity materials using more than one TM, e.g., Na 4 Cr III Mn II (PO 4 ) 3 . [30,46] Other NaSICON materials, such as Na 3 V III Cr III (PO 4 ) 3 , [47] Na 2 Ti IV V III (PO 4 ) 3 , [48,49] Na 3 Ti IV Mn II (PO 4 ) 3 , [50,51] Na 2 Ti IV Fe III (PO 4 ) 3 , [52] Na 4 V III Mn II (PO 4 ) 3 , [34,35,[53][54][55][56] Na 3 V III Fe III (PO 4 ) 3 [29,35] and Na 4 V III Ni II (PO 4 ) 3 , [35] have also been studied to enable high reversible capacity, longevity and as anode and cathode electrodes.

To access three electrons (i.e., 3 Na + ions per 2 TMs per formula unit) in any mixed TM NaSICON, two criteria must be met: i) the TM must be in the +2 and +3 oxidation states in the fully discharged state (i.e., 4 Na + ions), and ii) at least one of the TMs must accommodate multiple redox active oxidation states. Compounds that meet these criteria are mainly either Mn or V-based NaSICONs, e.g., for example Na 2 Ti IV V III (PO 4 ) 3 , Na 3 Ti IV Mn II (PO 4 ) 3 , Na 4 Mn II V III (PO 4 ) 3 , and Na 4 Cr III Mn II (PO 4 ) 3 , which undergo reversible electrochemical reactions up to three electrons with average voltages of ~2.4, ~3.2, ~3.4, and ~4.1 V vs. Na/Na + . While Na 4 Cr III Mn II (PO 4 ) 3 exhibits a high theoretical energy density (566 Wh kg -1 ) for a gravimetric capacity of 160 mAh g -1 , [30] only ~40 mAh g -1 was reported at voltages below 1.6 V vs. Na/Na + .

However, there remains a wide chemical space yet to be explored, which could yield new

Results:

Structural features of Na x MM'(PO 4 ) 3 electrodes

Depending on the TM species, temperature and/or Na content (x) per formula unit (f.u.)

NaSICON materials adopt typically a rhombohedral ( ), [45,57] an ordered monoclinic ( or ), [20,21,47] or, in the specific case of -Na 3 Ti III Ti III (PO 4 ) 3 , a triclinic structure ( ). [58] For example, Na 4 V II V III (PO 4 ) 3 and Na 4 Fe II Fe III (PO 4 ) 3 [59] are rhombohedral, while Na 3 V III V III (PO 4 ) 3 , [57] Na 3 Fe III Fe III (PO 4 ) 3 , [43] Na 3 Ti III Ti III (PO 4 ) 3 and Na 3 Cr III Cr III (PO 4 ) 3 form monoclinic [58] structures at room temperature due to Na/vacancy orderings. Figure 1 shows the archetypal structure of the rhombohedral Na x MM'(PO 4 ) 3 . the "lantern units" --a recurrent motif of NaSICON structures. The lantern units assemble into a 3-D MM'(PO 4 ) 3 framework producing two types of available sites for Na in the rhombohedral phase (Figure 1b): Na1 (one per f.u.) and Na2 (three per f.u.). The sixcoordinated Na1 sites are located between two (M,M')O 6 octahedra, while the Na2 sites remain eight-coordinated by the O atoms of the PO 4 tetrahedra. At x = 4 in Na x MM'(PO 4 ) 3 , all the Na positions (i.e., 1 x Na1 and 3 x Na2 per f.u.) are fully occupied and the structure is ordered and rhombohedral. In several NaSICONs, the Na + ions are fully ordered at x = 3 and the structure adopts monoclinic ( ) symmetry (Figure 1c), where the Na1 site splits into two distinct sites (i.e., fully occupied Na1a and Na1b). The Na2 site splits into five sites (3 fully occupied, Na2a, Na2b and Na2c, and 2 vacancies).

Na intercalation in Na x M 2 (PO 4 ) 3 NaSICONs

The energetics of Na intercalation and the corresponding average voltages (see system. The convex hull highlights the thermodynamic phase behavior at 0 K of Na (de)intercalation from/into the Na x M 2 (PO 4 ) 3 frameworks. Since the convex hull at 0 K does not include any entropic and pV effects, it informs on the propensity of Na to bond with the M 2 (PO 4 ) 3 frameworks. Among the M = Ti, V, Cr, Mn, Fe, Co and Ni compounds, NaSICONs based on Ti, V, Cr and Fe have been reported experimentally. [28], [36], [45], [43], [44] From previous experimental work, the redox activities of Ti and V provide accessible Na concentrations ranging in 1 ≤ x ≤ 4, while Na x Cr 2 (PO 4 ) 3 spans a narrower range of Na concentration (1 ≤ x ≤

methodology
3), being limited by the Cr IV /Cr III redox couple. Notably, the Cr III /Cr II redox couple has not been cycled reversibly in NaSICON frameworks so far. ) is found to be slightly metastable (~1.8 meV/atom above the stability line) and will phase separate into Na 1 V IV V IV (PO 4 ) 3 and Na 3 V III V III (PO 4 ) 3 .

Consequently, the stable structures lying on the convex hull (e.g., Na 3 V III V III (PO 4 ) 3 ) give rise to a "step" in the calculated Na voltage curve (see Figure S2). The convexity, that is the depth of the convex hull, is an indication of the thermodynamic stabilities of the specific Na/vacancy configurations. The calculated ground states on the convex hull do not necessarily correspond to topotactic structures since we allow changes to the symmetry (e.g., rhombohedral  monoclinic) of the host during deintercalation of Na from Na 4 MM'(PO 4 ) 3 in our calculations.

Analyzing Figures 2a and2b, the following general trends can be deduced:

i)

The end member compositions in Na x M 2 (PO 4 ) 3 , i.e., x = 1 and 4, typically adopt the rhombohedral ( or ) space group. Exceptions to this trend are Na 4 Cr II Cr III (PO 4 ) 3

and Na 4 Mn II Mn III (PO 4 ) 3 , which are predicted to have triclinic ( ) symmetry due to Jahn-Teller distortions of high-spin Cr 2+ and Mn 3+ d 4 ions. These have not yet been observed experimentally.

ii) A global minimum in the formation energy curves, i.e., the highest magnitude of formation energy, is present across all transition-metal Na x M 2 (PO 4 ) 3 systems at x = 3 (Figure 2a). [41] iii) Among all the Na 3 M III M III (PO 4 ) 3 compounds investigated, Na 3 Cr III Cr III (PO 4 ) 3 displays the deepest or lowest (i.e., most favorable) formation energy, followed by To validate our methodology, we have benchmarked our results against experimental observations in Na x V 2 (PO 4 ) 3 , which has been extensively studied. [57], [60] Experimentally, Na 3 V III V III (PO 4 ) 3 adopts monoclinic symmetry ( ) at room temperature, which is in qualitative agreement with our DFT calculations identifying the monoclinic ordering ( ) as the stable structure at x = 3. From Figure 2a, Na 3 V III V III (PO 4 ) 3 is on the convex hull with Na 1 V IV V IV (PO 4 ) and Na 4 V II V III (PO 4 ) 3 , which adopt the rhombohedral space groups and , respectively, in our calculations; this is consistent with X-ray diffraction experiments. [61], [36], [57] In Figure 2b, the extraction of 2 Na atoms from Na 3 V III V III (PO 4 ) 3 gives rise to an average voltage of ~2.96 V vs. Na/Na + and corresponds to the V IV /V III redox, while the insertion of Na into Na 3 V III V III (PO 4 ) 3 results in a voltage of 1.54 V (associated with the V III /V II redox). Experimentally, both the V redox couples (V IV /V III and V III /V II ) are involved in two-phase reaction mechanisms, yielding average voltages of ~3.40 V and ~1.63 V vs.

Na/Na + , respectively, [60] which are in reasonable agreement with our computed values.

Notably, both GGA+U and hybrid functionals underestimate the computed voltages of Na x V 2 (PO 4 ) 3 (see Section S2 in SI). Another important assessment to verify the validity of our simulations is ensuring that changes in oxidation states of the TMs upon Na extraction/insertion are captured. We verify that the correct redox processes do occur in our calculations by tracking the on-site magnetic moments on the TM atoms (Table S1), as well as the transition metal coordination environments (Table S2), which are discussed in Section S3 of the SI.

Whenever we could not assign the expected oxidation state of the TMs directly from the computed magnetic moments (e.g., in cases of Fe, Co and Ni based NaSICONs), we inspected the appropriate density of states (DOS). The DOS are reported in Section S4 of the SI (Figures S9-S12). As an example, we show in Figure 3, the total and atom-projected DOS in Na x V 2 (PO 4 ) 3 , for x = 1, 3 and 4, which are the ground states identified in Figure 2.

The DOS are displayed for spin up and spin down states in the range of interest (-2.5 --2 eV) around the Fermi energy. Unless a given structure is metallic, the top of the Fermi energy is arbitrarily set to the top of the valance band. Panels a and b in Figure 3 show the results for Na 1 V IV V IV (PO 4 ) 3 and Na 3 V III V III (PO 4 ) 3 , respectively. From the projected DOS, the valence bands near the fermi energy are occupied by 3d electrons of vanadium (red line), while the 2p O states (orange line) lie at lower energies. In general, the vanadium 3d states tend to dominate the valence band and shift at higher energies as more Na is inserted, closing the band gap (> 1eV for both compounds). However, the monoclinic distortion of Na 3 V III V III (PO 4 ) 3 is responsible for an increase in band gap (~1.8 eV) contrary to the Na = 1 and Na = 4 trends. As expected, the intercalation of more Na + ions as in Na 4 V III V II (PO 4 ) 3 (Figure 3c) increases (destabilizes)

further the valance band, with this composition becoming gapless and a Fermi energy dominated by V(3d) states when GGA+U is used. For the Na 4 V III V II (PO 4 ) 3 , the band gap opens to ~0.3 eV when the more accurate HSE06 hybrid functional is used at the relaxed HSE06 structure. Very recent hybrid functional simulations on Na x Ti 2 (PO 4 ) 3 have verified these orders of magnitudes, with Na 4 Ti III Ti II (PO 4 ) 3 showing a small band gap of ~0.59 eV. [62] In the case of Na x Ni 2 (PO 4 ) 3 the projected DOS of Figure S12 suggests metallic behavior at all the Na concentrations explored (i.e., x = 1, 3 and 4). Na x Mn 2 (PO 4 ) 3 and Na x Co 2 (PO 4 ) 3

show band gaps >1.0 eV at x = 3 and 4, but become metallic at x =1, as shown by Figures S9 andS11. Na 1 Fe IV Fe IV (PO 4 ) 3 also appears, surprisingly, metallic (Figure S10), which introduces difficulties in localizing the Fe III /Fe IV redox couple. Therefore, due to the apparent metallic behavior of some of these systems (Fe, Co and Ni-based NaSICONs), we were unable to assign unique oxidation states to their TMs. Notably, DFT is a ground state theory and not adequate in the prediction of accurate band gaps. [63]

Na + ion intercalation in Na x MM'(PO 4 ) 3 NaSICON structures

We extend our analysis to the reversible extraction of Na from mixed Na 4 MM ' (PO 4 ) 3 compounds where, for the sake of simplicity, the ratio M:M' is fixed as 1:1. While all the M:M' is also reported in Figures 4a and4c as a reference.

Among the Na x TiM'(PO 4 ) 3 systems, only Na x Ti 2 (PO 4 ) 3 , Na x TiV(PO 4 ) 3 and Na x TiMn(PO 4 ) 3 have been experimentally investigated over the entire Na composition range (1 ≤ x ≤ 4), [28], [49], [48], [50], [51] while Na x TiCr(PO 4 ) 3 [52], [64] and Na x TiFe(PO 4 ) 3 [52] have been studied in limited ranges of 1 ≤ x ≤ 3 and 2 ≤ x ≤ 4, respectively. Na x TiCo(PO 4 ) 3 and Na x TiNi(PO 4 ) 3 are theoretical compounds and have not yet been reported experimentally. In general, the end member compositions at x = 1 and 4 (Figures 4a and4b) of the Ti-M' NaSICON electrodes adopt rhombohedral symmetry ( or , or ) except for Na 4 Ti III Cr II (PO 4 ) 3 ( ) and Na 1 Ti IV Fe IV (PO 4 ) 3 ( ), which are theoretically derived. This could be due to Jahn-Teller effects caused by the high-spin d 4 electronic configurations of Cr 2+ and Fe 4+ , respectively, as well as the low stability of Fe 4+ . For Na x TiM'(PO 4 ) 3 with M'= V, Mn, Co and Ni, a global minimum in the convex hull is located at x = 3, and the magnitudes of the formation energy follow the trend Ni > Co > Mn > V. On the other hand, the global minima for the M' = Ti, Cr and Fe NaSICONs are located at x = 2 (Figure 4a), which is driven by the high stability of their stable oxidation states Ti 4+ , Cr 3+ and Fe 3+ , respectively. Overall, the depth of the convex hulls (i.e., maximum magnitude of formation energies) for mixed TM NaSICONs follows the trend Fe > Cr ~ Ni > Co > Ti > Mn ~ V. and4d present the computed voltages for Na extraction from Na x TiM'(PO 4 ) 3 in the composition range 1 ≤ x ≤ 4. Expectedly, the compound exhibiting the lowest intercalation voltage (~0.96 V) is Na x Ti 2 (PO 4 ) 3 for 3 ≤ x ≤ 4, which overestimates the experimental voltage (~0.4 V) reported by Senguttuvan et al. [28] The DFT data of Figure 4a suggest that an unreported stable phase appears at Na 1.5 TiFe(PO 4 ) 3 , but with an oxidation state for Fe of 3.5 (Fe IV + Fe III ). Although, the highest Na intercalation voltage (~4.99 V) is computed for Na x TiFe(PO 4 ) 3 between x = 1.5 and x = 1, the oxidation state of Fe would have to be +4, which may not be experimentally accessible in a reversible manner. We find that Na 1 Fe IV Fe IV (PO 4 ) 3, where iron is expected to be Fe 4+ appears metallic from the DOS in

Figures 4c

To assess the quality of our computed data, we analyzed in detail the experimentally reported data for Na x TiV(PO 4 ) 3 . At x = 1, 2 and 4, Na x TiV(PO 4 ) 3 has rhombohedral symmetry ( ), while the structure of Na 3 TiV(PO 4 ) 3 is monoclinic ( ). [49], [48] Our calculations reproduce the rhombohedral symmetry ( ) at x = 1 and 4, and the monoclinic symmetry ( ) at x = 3, whereas we find a triclinic structure ( ) for Na 2 TiV(PO 4 ) 3 . [48] The Na x TiV(PO 4 ) 3 convex hull (Figure 4a) displays a minimum at x = 3, analogous to the observed minima in the Na x M 2 (PO 4 ) 3 systems (Figure 2a); this also provides the largest step in the corresponding voltage profile of Figure 4c. In particular, the Na extraction from Na 4 Ti III V II (PO 4 ) 3 takes place through the redox couples V III /V II 1.54 V (experimentally ~1.6 V [49], [48] ), Ti IV /Ti III 2.50 V (~2.1 V), and V IV /V III 2.90 V (~3.4 V) vs. Na/Na + , respectively. The mechanism of Na extraction from Na 4 Ti III V II (PO 4 ) 3 has been shown to occur via a solidsolution mechanism for the V III /V II redox couple, followed by a two-phase reaction for both the Ti IV /Ti III and V IV /V III redox couples. [48] Figure 5 The calculated average voltage per extracted Na + ion vs. Na/Na + over the x range [1-4] in Na x MM'(PO 4 ) 3 compounds, where M = M' or M = Ti and M' = V, Cr, Mn, Fe, Co and Ni.

Due to the existence of multiple voltage steps in Na x MM'(PO 4 ) 3 (where M' = Ti, V, Cr, Mn, Fe, Co and Ni and M = Ti), we compare directly the computed average voltage per Na + exchanged over the entire composition range (i.e., 1 ≤ x ≤ 4), as shown in Figure 5.

Importantly, the overall average voltage monotonically increases from Ti to Ni in cases

where M = M' and M = Ti, which is consistent with trends in standard reduction potentials of the transition metals. We note that the M 4+ ions, in particular, become more oxidizing as one crosses the transition series from Ti to Ni, with the lack of voltage increase from Mn to Fe reflecting the stability of their 3d 5 configurations. Specifically, the calculated voltages increase from ~1.96 V (for Ti) to ~4.40 V vs. Na/Na + (for Ni) in Na x M 2 (PO 4 ) 3 (green line in exceptions to this trend are for V and Nb-based NaSICONs, which can also operate on the V V /V IV couple (see mixed Na x VM'(PO 4 ) 3 systems in Figure 6) and the Nb V /Nb IV couple in Na x NbTi(PO 4 ) 3 not discussed here. [65] 2.5 Thermodynamic stabilities of selected Na x M 2 (PO 4 ) 3 compounds

The calculated quaternary phase diagrams for the systems Na-P-O-Mn, Na-P-O-Co, and Na-P-O-Ni (Figures 7 andS34 in SI) reflect the stabilities of the unreported NaSICONs with these transition metals compared with the stable elemental, binary, ternary and quaternary phases. While it remains extremely challenging to characterize quaternary (and beyond) phase diagrams solely based on sparse experimental data, theory provides a viable alternative for investigating such complex systems. We calculated all the binary/ternary/quaternary compounds available in the inorganic crystal structure database (ICSD) [66] within the Na-P-O-Mn, -Co and -Ni quaternary systems. shows that Na 3 Mn III Mn III (PO 4 ) 3 and Na 1 Mn IV Mn IV (PO 4 ) 3 are stable compounds.

Na 3 Mn III Mn III (PO 4 ) 3 has the structure of the monoclinic ) polymorph, whereas Na 1 Mn IV Mn IV (PO 4 ) 3 has a rhombohedral ( ) structure. From our computations, we observed a Jahn-Teller distortion driven by the Mn 3+ 3d 4 ions and highlighted by four longer equatorial bonds ~2.10 Å and two shorter axial bonds of ~1.96 Å . No Ni-containing NaSICON phases appear stable in the Na-Ni-P-O system (Figure S34). In the Na-Co-P-O phase diagram (Figure 7b), only Na 3 Co III Co III (PO 4 ) 3 appears stable and exhibits rhombohedral symmetry ( ). From this analysis, it is predicted that Na 3 Mn III Mn III (PO 4 ) 3 , Na 1 Mn IV Mn IV (PO 4 ) 3 , and Na 3 Co III Co III (PO 4 ) 3 have potential for successful experimental synthesis, given their thermodynamic stabilities.

Discussion

Using first-principles calculations, we investigated the electrochemical properties of 28 transition-metal-based NaSICONs as cathode materials for NIBs. To guide our discussion, we summarize the reported experimental voltages, the experimental gravimetric, and the theoretical gravimetric capacities of some Na x MM'(PO 4 ) 3 systems in Table I. For a transparent and fair comparison between theoretical and experimentally reported capacities in Table I, the reported capacities are renormalized from the reported composition (pristine)

to the molecular weight of the Na 4 MM'(PO 3 ) 4 . The theoretical capacities are also computed with respect to the Na 4 MM'(PO 3 ) 4 compounds. processes. ? is used whenever the space group of a specific phase is unknown or not reported. Intercalation voltages account for the number of Na + ions exchanged. Whenever available, the space group of each phase is reported.

For a transparent comparison between theoretical and experimentally reported capacities, the reported capacities are renormalized to the molecular weight of Na 4 MM'(PO 3 ) 4 . The theoretical capacities (Theo.) are also computed with respect to the Na 4 MM'(PO 3 ) 4 compound. oxidation state, with Cr 2+ (high-spin d 4 ) being typically an unstable oxidation state for Cr and easily prone to Jahn-Teller distortions in octahedral environments, e.g., NaSICON. [START_REF] Cox | Transition Metal Oxides: An Introduction to Their Electronic Structure and Properties[END_REF] Na x V 2 (PO 4 ) 3 (NVP) exhibits the 2 nd highest measured voltage of the single TM NaSICONs, with a computed average voltage of ~2.48 V vs. Na/Na + (Figure 5). With the same number of maximum electrons exchanged in Na x M 2 (PO 4 ) 3 (where M= Ti or V), and Ti being lighter than V, Na x Ti 2 (PO 4 ) 3 is expected to show the largest gravimetric capacity. Na x Ti 2 (PO 4 ) 3 and NVP provide the largest theoretical gravimetric capacities (~170.1 mAh/g and ~167.9 mAh/g for Na 4 M 2 (PO 4 ) 3 ) in single TM NaSICONs. In theory, additional capacity and an increased voltage could be achieved by exploiting the V V /V IV redox couple in the reaction Na 1 V IV V IV (PO 4 ) 3  Na + + 1e -V V V IV (PO 4 ) 3 . While the seminal report by Gopalakrishnan et al. suggested the possibility of chemically extracting the last Na + ion to form V V V IV (PO 4 ) 3 , [START_REF] Gopalakrishnan | [END_REF] subsequent attempts to do so have proven unsuccessful. Note that, V V V IV (PO 4 ) 3 is not a thermodynamically stable compound [69] and is predicted to decompose into VPO 5 + VP 2 O 7 .

We assessed the existence of uncharted single-TM NaSICONs, which are Ni, Co and Mn.

Except for Na 1 Mn IV Mn IV (PO 4 ) 3 , Na 3 Mn III Mn III (PO 4 ) 3 and Na 3 Co III Co III (PO 4 ) 3 that appear stable according to our convex hulls but have not been reported experimentally, our analysis of the phase diagrams indicates that Ni-, and Co-NaSICONs are generally unstable. However, Mn has been mixed effectively with Ti, V and Cr, while Ni-V mixed NaSICONs have also been made (Table I). [34] In particular, Zhou et al. [35] and later Chen and collaborators [34] have synthesized Mn II -containing Na 4 Mn II V III (PO 4 ) 3 , which upon Na extraction exploits the Mn III /Mn II redox couple and Mn 3+ 3d 4 is Jahn-Teller active.

In the case of mixed Na 4 MM ' (PO 4 ) 3 NaSICONS, we kept the M:M' ratio to 1:1. From an extensive assessment of polyanion electrodes, Masquelier and Croguennec [17] observed that the redox potentials of specific TMs remain nearly independent [70] from those of other TMs that may be present. This general rule is also well supported by the simulations that are summarized in Figure 6. Based on this general principle, one can envision selected combinations of TMs delivering reversible high voltages vs. Na/Na + . For example, the Ti IV /Ti III redox couple showed a similar voltage [17] (∼2.1 V vs. Na/Na + ) in three related NaSICONs systems: Na x Ti 2 (PO 4 ) 3 , Na x TiNb(PO 4 ) 3 and Na x TiFe(PO 4 ) 3 , respectively.

Notably, in experiments, there is a consistent shift of ~0.3 V between the voltages measured vs. Na/Na + compared to those measured vs. Li/Li + ).

The Na x TiV(PO 4 ) 3 system has been shown to provide the largest gravimetric capacity (~187.1 mAh/g, Table I) among NaSICONs as a result of the accessibility of V IV /V III /V II and Ti IV /Ti III redox couples. Nevertheless, the computed average voltage (~2.3 V, Figure 6) of Na x TiV(PO 4 ) 3 highlights the low energy density of the system. [49], [48] More promising in terms of overall capacities and voltages, and thus energy density, are Na x TiCr(PO 4 ) 3 , [52], [64] Na x TiMn(PO 4 ) 3 , [50], [51] Na x CrMn(PO 4 ) 3 , [30], [46] Na x VCr(PO 4 ) 3 , [47] and Na x VMn(PO 4 ) 3 , [34,54] but the reversible extraction/insertion of Na from some of these systems remains to be verified. Excluding some specific exceptions, especially Fe, Co and Ni based NaSICONs since they appear metallic, our simulations are also able to capture the mechanism of Na intercalation in mixed NaSICON electrodes. For example, Na 3 VFe(PO 4 ) 3 has been reported to reversibly insert Na with the activation of the Fe III /Fe II and V III /V II redox couples, respectively (see Table I), [35] and this behavior is reproduced by our calculations (Figure 6 and Figure S21 of SI). Similar conclusions can be advanced for mixed TiV, TiCr, VCr and CrMn NaSICON materials.

From our analysis (Figure 6), specific unexplored combinations of TMs appear to be worthy of future study. Our calculations suggest promising NaSICON compositions in terms of average voltages, such as TiCo, VCo, CrFe and CrCo NaSICONs as reported in Figure 5 and Table 1. Both CrFe and CrCo NaSICONs are predicted to exploit the high voltage Cr IV /Cr III redox couple, which has been proven effective into other mixed NaSICONs, e.g., TiCr, VCr and MnCr. [64], [47] Although the addition of Co into Ti, V and Cr appears interesting, our analysis of the Na x Co 2 (PO 4 ) 3 system suggests that only Na 3 Co III Co III (PO 4 ) 3 with Co 3+ 3d 6

(low spin) is thermodynamically stable, but this compound has not been yet synthesized.

While one would target compounds with Co 2+ to ensure high voltages , it appears possible to introduce Co in smaller proportions [71] than a 1:1 ratio, but more experimental and theoretical work is required.

We do not anticipate Ni to be of importance in mixed NaSICON systems due to the high instability of quaternary Na x Ni 2 (PO 4 ) 3 . Manoun et al. [72] have briefly reported the synthesis of Na 4 Cr III Ni II (PO 4 ) 3 . A recent report by Zhou et al. [35] claims the extraction of Na from Na x VNi(PO 4 ) 3 (Table I) involving the redox couples V IV /V III , V V /V IV , and Ni III /Ni II .

Nevertheless, we speculate that the Ni content of the materials is not in a 1:1 ratio but appears to be significantly lower. Furthermore, the highest voltage accessed (<4.3 V vs.

Na/Na + ) during the electrochemical cycling is lower than that of the expected Ni redox couple. While our DFT calculations seem to verify the activity of the V IV /V III and V V /V IV redox couples in Na x VNi(PO 4 ) 3 (Figure 6), we are unable to verify the oxidation states involved in the Ni III /Ni II reaction from the computed magnetic moments as the material becomes metallic in our simulations (Figure S12). Our findings also cast doubt on the activity of Ni in this compound as reflected by the limited gravimetric capacity reported experimentally (~ 80 mAh/g). [35]

Conclusions

Using ab initio density functional theory and thermodynamics, we explored the full chemical map of 3d TM-based NaSICON phosphates Na x MM'(PO 4 ) 3 (M, M'= Ti, V, Cr, Mn, Fe, Co and Ni) cathode materials for high-energy density and safe sodium-ion batteries. We identified the ground state structures at various Na contents as well as the corresponding Na-intercalation voltages and redox processes for 28 distinct NaSICON compositions of which only 13 have been reported experimentally.

Further, we investigated the Na-intercalation properties of previously unreported We performed a complete analysis of the wide-ranging chemical space of NaSICON phosphate cathode materials for Na-ion batteries, and our work can be used to guide further experimental synthesis of the new and promising compositions identified here. We plan to synthesize promising electrode compositions identified in this work.

Methodology

We used the Vienna ab initio simulation package(VASP) [73], [74] for ab initio density functional theory (DFT) total energy calculations. The projector augmented wave (PAW) potentials were used for the core wave-functions. The Perdew-Burke-Ernzerhof (PBE) parameterized spin-polarized generalized gradient approximation (GGA) was used for the exchange and correlation energy. [75] The strong on-site coulomb correlation of 3d electrons of the TM (Ti, V, Cr, Mn, Fe, Co and Ni) is addressed with Hubbard's U correction [76] according to the Dudarev method. [77] The effective U parameters used are 4.0 eV for Ti, [28] 3.1 for V, 3.5 for Cr, 3.9 for Mn, 4.0 for Fe, 3.4 for Co, and 6.0 for Ni, [78] and were fitted to match the experimental formation energies of binary oxides. We also introduced an empirical correction of 1.36 eV/O 2 [79] to remedy the spurious error originating from the well-known over-binding of O 2 while using GGA.

The periodic wave functions were expanded in terms of plane waves up to a kinetic energy cut-off of 520 eV. The PAW potentials used to describe the core electrons were Na 08Apr2002 3s 1 , P 17Jan2003 2s 2 3p 3 , O 08Apr2002 2s 2 2p 4 , Ti 08Apr2002 3d 3 4s 1 , V_pv 07Sep2000 3p 6 3d 4 4s 1 , Cr 06Sep2000 3d 5 4s 1 , Mn 06Sep2000 3d 6 4s 1 , Fe 06Sep2000 3d 7 4s 1 , Co 06Sep2000 3d 8 4s 1 , and Ni 06Sep2000 3d 9 4s 1 . Additionally, a -centered Monkhorst-Pack [80] k-point mesh with 25 subdivisions along each reciprocal lattice vector was applied to all structures. Using these settings, the total energy of each structure was converged to within 10 -5 eV/cell, atomic forces within 10 -2 eV and the stress within 0.29 GPa.

An intercalation battery based on the NaSICON cathode electrodes implies the reversible insertion/extraction of Na + -ions into/from the Na y MM'(PO 4 ) 3 framework according to the redox reaction of Equation 1.

(

where y and x are the initial and final Na concentration in NaSICON framework and ΔG 0 is the change of Gibbs energy at 0K for the reaction of Equation 1. Here, we approximated the Gibbs energy of each component by the DFT total energies (i.e., ), thus neglecting the pV and entropic contributions. The average voltage across an intercalation extent ( ) can be calculated from the ΔG 0 , as in Equation 2.

(

) 2
where is the Na chemical potential (set to bulk Na metal) and is the Faraday constant.

To establish the general phase behavior of Na (de)intercalation into the Na y MM'(PO 4 ) 3 structure, we monitored, using Equation 3, the formation energies ( ) of various orderings at different Na concentrations, x (1 ≤ x ≤ 4), with respect to the DFT energies of the fully discharged (i.e., E[Na 4 MM'(PO 4 ) 3 ]), and fully charged (E[NaMM'(PO 4 ) 3 ]) configurations.

(3) Different NaSICON structures were studied according to specific Na concentrations in Na x MM'(PO 4 ) 3 with M and M'= Ti, V, Cr, Mn, Fe, Co, and Ni. In choosing the amount of mixing of TM in these NaSICON structures, we have considered only two distinct situations: i) M= M', which leads to Na x M 2 (PO 4 ) 3 , and ii) M and M' in the ratio 1:1.

When fully sodiated, the fully ordered rhombohedral symmetry ( ) of the high temperature NaSICON structure is typically observed (e.g., Na 4 Fe III Fe II (PO 4 ) 3 [59] ) and is therefore our starting model to study Na removal and TM mixing. Na vacancies are created in the fully sodiated structure (Na 4 MM'(PO 4 ) 3 ), resulting in Na x MM'(PO 4 ) 3 , where the Na content varies in the range 1 ≤ x ≤ 4 in steps of Δx = 0.5. The possible configurations originating from the various orderings of Na and vacancies (Va) in the NaSICON are obtained using the pymatgen library. [81] A ranking according to the classical Ewald energy [82] based on integer charges (i.e., Na=+1, P=+5, O=-2) and variable charge on the TM (+2, +3 and +4) is applied to limit the number of possible structures to a computationally tractable level. For mixed TM NaSICONs, the ordering of M and M' is simultaneously performed with that of Na and Va.

DFT calculations are performed on the primitive cell and supercells (2 1 1) of these orderings. A formula unit of the fully discharged NaSICON (Na 4 MM'(PO 4 ) 3 ) contains 21 atoms.

Figure 1 a

 1 Figure 1 a The rhombohedral ( ) structure of Na 4 MM'(PO 4 ) 3 with M and M' = Ti, V, Cr, Mn, Fe, Co and Ni. The "lantern unit" is made up of two MO 6 (or M'O 6 ) octahedra (red) sharing corners with PO 4 3tetrahedra (blue). The gray box in (a) highlights the arrangement of lantern units. Two distinct Na ions, i.e., Na1 (black) and Na2 (green), are present in the rhombohedral NaSICONs. Panel b lantern representation of the Na disordering (with partial Na occupancies on Na1 and Na2 sites) in the rhombohedral Na 3 MM'(PO 4 ) 3 polymorph. Panel c shows a local view (lantern unit) of the complete ordering of Na + ions in the corresponding Na 3 MM'(PO 4 ) 3 monoclinic ( ) phase. Panel d shows the local environment of sites Na1b (violet) and Na2a (yellow) in the same monoclinic arrangement. Vacancies are represented by open circles.

in Section 5 )and b of Figure 2 .

 52 for single TM (M=M') NaSICON electrodes are shown in panels a The computed formation energies for the most stable orderings in Na x M 2 (PO 4 ) 3 are displayed in Figure 2a as a function of Na concentration, where the lower energy envelopes form the so-called convex hull (solid lines) for a given Na x M 2 (PO 4 ) 3

Figure 2 Figure S2 ,

 2S2 Figure 2 Panel a shows the computed formation energies and respective convex hulls for Na vacancy orderings as a function of Na concentration (x) in Na x M 2 (PO 4 ) 3 where M = Ti, V, Cr, Mn, Fe, Co and Ni. Panels b and c show the intercalation voltages vs. Na/Na + for the experimentally known (M = Ti, V, Cr and Fe) and predicted (M = Mn, Co and Ni) Na x M 2 (PO 4 ) 3 compounds, respectively. Only the stable orderings forming the convex hull are displayed in panel a. The energies of unstable configurations for each convex hull are shown in Figures S1-S7 of the SI.

Na 3 Fe

 3 III Fe III (PO 4 ) 3 > Na 3 Ti III Ti III (PO 4 ) 3 ≈ Na 3 V III V III (PO 4 ) 3 respectively, which reflects the stability of the M 3+ oxidation states of these TMs. The stability of Na 3 Cr III Cr III (PO 4 ) 3 , in particular, is enhanced by the large ligand field stabilization energy of the Cr 3+ (3d 3 ) cation. The other compounds, Na 3 Co III Co III (PO 4 ) 3 , Na 3 Mn III Mn III (PO 4 ) 3 and Na 3 Ni III Ni III (PO 4 ) 3 , display shallower (i.e., less favorable) formation energies (< 25 meV/atom, Figure2a). The shallow convex hulls for Na X Mn 2 (PO 4 )3 and Na x Ni 2 (PO 4 ) 3 are caused by the Mn 3+ and Ni 3+ ions which are both Jahn-Teller active. The magnetic moment of Co 3+ Na 3 Co III Co III (PO 4 ) 3 is ~3.0 μ B , indicating an intermediate spin state and possible Jahn-Teller activity. iv) Unsurprisingly, the deep minima in the convex hull plot of Figure 2a lead to a large voltage step as shown in Figures 2b and 2d, following the sequence Cr > Fe > Ti > V > Co > Mn > Ni. v) The M(III/II) redox couples that correspond to Na contents from x = 3 to x = 4, follow the voltage trend Ni > Co > Mn > Fe > Cr > V > Ti, while the M(IV/III) couples (x = 1 to x = 3) follow the sequence Ni > Co > Fe > Cr > Mn > V > Ti. Thus, Ni and Co (V and Ti) display the highest (lowest) average voltages for both M redox couples.

Figure 3

 3 Figure 3The total (gray) and atom projected (V red, O orange and Na yellow) DOS of Na x V 2 (PO 4 ) 3 NaSICON (NVP), where for Na concentrations x = 1, 3 and 4. The vertical line denotes the Fermi energy level and E g is the calculated band gap at the GGA+U level of theory.

  combinations have been considered, to simplify our discussion we concentrate on the cases where M = Ti, while M' = Ti, V, Cr, Mn, Fe, Co and Ni. Other mixed Na 4 MM ' (PO 4 ) 3 are reported in Section S5 of the SI (FiguresS13-S33). The main reason behind choosing Tibased mixed TM NaSICONs is their maximum experimental electrochemical data availability as compared to other combinations. Figures4a and 4bshow the calculated convex hulls for experimentally reported and theoretical Na x TiM ' (PO 4 ) 3 compounds, respectively, with Na contents in the range 1 ≤ x ≤ 4. Figure4cdisplays the corresponding voltages. Na x Ti 2 (PO 4 ) 3

Figure 4

 4 Figure 4 Panels a and b show the computed convex hulls as a function of Na concentrations (x) for Na x TiM'(PO 4 ) 3 , where M' = Ti, V, Cr, Mn, Fe in panel a and M' = Co and Ni in panel b. Panels c and d show the corresponding Na (de-)intercalation voltages over changes in Na-content. The energies of unstable configurations for each convex hull, as well as complete voltage curves for each compound, are given in the SI (Figure S13 -S33). Symbols in panels a and b, indicate stable ordered phases. To facilitate the visualization of the voltage curves near ~2.8 V vs. Na/Na + the Ti-Cr curve in panel c has been slightly lifted, but it is expected to overlap the Ti-V and Ti-Mn voltage curves.

Figure 5 ) 2 . 4 3 Figure 6 A

 52436 Figure 5), while the maximum voltage is ~3.73 V for Ti-Ni in the case of Na x MM'(PO 4 ) 3 (M = Ti) systems (orange line). Notably, the voltage values for Na x TiM'(PO 4 ) 3 are consistently

Figure 6 2 .Figure 6 Figure 5 .

 6265 Figure 6 consolidates the computed voltages of 28 plausible combinations of 3d TM in Na x MM'(PO 4 ) 3 with M and M' in a 1:1 ratio. In Figure 6, boxes outlined in red represent the

Figure 7

 7 Figure 7 Panels a and b show portions of the computed phase diagrams of Na-P-O-Mn and Na-P-O-Co at 0 K. Two stable NaSICON structures, i.e., NaMn 2 (PO 4 ) 3 and Na 3 Mn 2 (PO 4 ) 3 are observed in the phase diagram of Na-P-O-Mn, which are in equilibrium with the highly stable binary (P 2 O 5 and MnO 2 ) and ternary (NaMnO 2 and NaPO 3 ) compounds. In panel b, a stable NaSICON structure, i.e., Na 3 Co 2 (PO 4 ) 3 is observed in the phase diagram of Na-P-O-Co, which is in equilibrium with the highly stable binary (P 2 O 5 and Co 3 O 4 ) and ternary (NaCoO 2 and NaPO 3 ) compounds.

Figure 7a shows the

  Figure 7a shows the computed phase diagram of the Na-P-O-Mn system, where the open

Na x Mn 2 (

 2 PO 4 ) 3 and Na x VCo(PO 4 ) 3 (1 ≤ x ≤ 4), among others. The calculated quaternary phase diagrams of the Na-P-O-M (M= Mn, Co and Ni) chemical systems indicate the instability of Ni and most Co-based NaSICONs, though the Na 1 Mn IV Mn IV (PO 4 ) 3 , Na 3 Mn III Mn III (PO 4 ) 3 and Na 3 Co III Co III (PO 4 ) 3 compounds are identified as stable compositions.

  and S 6+ ) and anions (mostly O 2-) arranging into phosphate, silicate and sulfate (PO 4

	3-, SiO 4 4-and SO 4 2-

Table I

 I Experimentally observed Na x MM'(PO 4 ) 3 phases along with the corresponding Na concentrations, structures, Na

+ reversible intercalation voltages, redox couples, mechanism of Na + intercalation, and gravimetric capacities. The intercalation voltages (in V) and theoretical (Theo.) capacities (in mAh/g) are shown. Rev. and Irrev. indicate experimentally reversible and irreversible

Acknowledgements P. C., C. M., A. K. C., and J.-N. C. are grateful to the ANR-NRF NRF2019-NRF-ANR073 Na-MASTER. P.C. and B. S acknowledge funding from the National Research Foundation under his NRF Fellowship NRFF12-2020-0012. L. C., D. C. and C. M. acknowledge the ANRT and TIAMAT for the funding of S. P.'s PhD thesis as well as the financial support of Région Nouvelle Aquitaine and of the French National Research Agency (STORE-EX Labex Project ANR-10-LABX-76-01). The computational work was performed on resources of the National Supercomputing Centre, Singapore (https://www.nscc.sg).

M-M' (M=M')

x = 1 2 3 4 Reported Capacity

Theo. Capacity Ti-Ti [28] ? P-1 R-3c 142.7 170.1 | 2.1 V, Ti IV /TI III , Rev. | 0.4 V, Ti III /TI II , Irrev. |

V-V [36,37,39,60] R-3c

Cr-Cr [45] R-3c

Theo. Capacity

Ti-V [48,49] R-3c [52] R-3c [50,51] R-3c R-3c R-3c ? 152.3 167.6 | 4. V, Mn IV /Mn III , Rev.| 3.5 V, Mn III /Mn II , Rev. | 2.1 V, Ti IV /TI III , Rev. | Ti-Fe [52] R-3c R-3c R-3c -111.5 | 2.4 V, Fe III /Fe II , Rev. | 2.1 V, Ti IV /TI III , Rev. | V-Cr [47] ? ? R-3c 85.7 111.7 | 4.1 V, V V /V IV , Rev. | 3.4 V, V IV /V III , Rev. |

Ti-Cr

V-Mn [34,54,56] R-3c R-3c R-3c R-3c 156.0 166.6 | 3.9 V, V V /V IV , Irrev. | 3.6 V, Mn III /Mn II , Rev. | 3.4 V, VI V /V III , Rev. | V-Fe [29,35] ? C2/c ? 98.1 110.8 | 3.3 V, V IV /V III , Rev.

| 2.5 V, Fe III /Fe II , Rev. | V-Ni [35] ? ? R-3c 80.0 110.2 | 3.9 V,V V /V IV :Ni III /Ni II , Rev. | 3.5 V, V IV /V III , Rev. | Cr-Mn [30,46] R-3c R-3c R-3c R-3c 160.5 166.2 | 4.4 V, Cr IV /Cr III , Rev. | 4.2 V, Mn IV /Mn III , Rev. | 3.6 V, Mn III /Mn II , Rev. | Among the Na x M 2 (PO 4 ) 3 compounds, Na x Cr 2 (PO 4 ) 3 delivers the highest measured voltage of ~4.5 V (~4 V calculated theoretically) vs. Na/Na + , [45] albeit with a limited capacity of just ~98 mAh/g, and corresponding to the extraction of 2 Na from Na 3 Cr III Cr III (PO 4 ) 3 to form Na 1 Cr IV Cr III (PO 4 ) 3 . Na 4 Cr 2 (PO 4 ) 3 has not yet been obtained by either a chemical or an electrochemical process. Na 4 Cr III Cr II (PO 4 ) 3 would require Cr to exist in a mixed, Cr III /Cr II