Mostafa Sadeghi
email: mostafa.sadeghi@inria.fr

Paul Magron
email: paul.magron@inria.fr

A Sparsity-promoting Dictionary Model for Variational Autoencoders

Keywords: generative models, variational autoencoders, sparsity, dictionary model, speech spectrogram modeling

Structuring the latent space in probabilistic deep generative models, e.g., variational autoencoders (VAEs), is important to yield more expressive models and interpretable representations, and to avoid overfitting. One way to achieve this objective is to impose a sparsity constraint on the latent variables, e.g., via a Laplace prior. However, such approaches usually complicate the training phase, and they sacrifice the reconstruction quality to promote sparsity. In this paper, we propose a simple yet effective methodology to structure the latent space via a sparsity-promoting dictionary model, which assumes that each latent code can be written as a sparse linear combination of a dictionary's columns. In particular, we leverage a computationally efficient and tuning-free method, which relies on a zeromean Gaussian latent prior with learnable variances. We derive a variational inference scheme to train the model. Experiments on speech generative modeling demonstrate the advantage of the proposed approach over competing techniques, since it promotes sparsity while not deteriorating the output speech quality.

Introduction

Unsupervised representation learning [START_REF] Bengio | Representation learning: A review and new perspectives[END_REF] is defined as the task of automatically extracting useful information from unlabeled data, in the form of a feature or representation vector, which can be subsequently used in downstream tasks. To that end, one successful approach, which has gained much attention in the past years, is based on variational autoencoders (VAEs) [START_REF] Kingma | Auto-encoding variational bayes[END_REF]. These models explicitly consider a latent vector which encapsulates some information about the data. A VAE model consists of a stochastic encoder (a recognition network), which transforms the input data into a latent space whose dimension is usually much lower than the original data space, and a stochastic decoder (a generative network), which produces data back in the original space from the encoded latent representation. VAEs have been successfully exploited in a variety of speech processing related tasks, such as speech enhancement [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF], source separation [START_REF] Nguyen | Deep variational generative models for audio-visual speech separation[END_REF], speech recognition [START_REF] Tan | Learning utterance-level normalisation using variational autoencoders for robust automatic speech recognition[END_REF], and speech generation and transformation [START_REF] Hsu | Learning latent representations for speech generation and transformation[END_REF].

A common practice in a VAE framework consists in assuming a standard normal distribution as the prior distribution of the latent space. However, this is not effective, as it might not efficiently capture the underlying distribution of the data. Structuring the latent space by imposing some meaningful constraints and distributions is thus of paramount importance to obtain more expressive and interpretable representations. A promising approach is to favor parsimonious or sparse latent representations, that is, feature vectors with mostly zeroes or small magnitude entries. Indeed, sparsity is known to promote a variety of benefits, such as increasing interpretability of the latent features and reducing the risk of overfitting [START_REF] Asperti | Sparsity in variational autoencoders[END_REF]. It has proven very successful for feature learning in many fields such as computer vision [START_REF] Gregor | Draw: A recurrent neural network for image generation[END_REF] or natural language processing [START_REF] Prokhorov | Learning sparse sentence encoding without supervision: An exploration of sparsity in variational autoencoders[END_REF].

In this regard, a sparse prior has been proposed in [START_REF] Mathieu | Disentangling disentanglement in variational autoencoders[END_REF] in the form of a mixture of two Gaussian distributions with fixed variances, where one of them is very small in order to encourage sparsity. A mixing parameter balances the contribution of the two distributions, and thus the amount of sparsity. Tonolini et al. [START_REF] Tonolini | Variational sparse coding[END_REF] proposed a variational sparse coding (VSC) framework, where they consider a Spike-and-Slab distribution [START_REF] Andersen | Bayesian inference for structured spike and slab priors[END_REF] for the encoder, which is a mixture of a Gaussian distribution, as in a standard VAE, and a sparsity-promoting Delta function. The prior distribution has the same form as the encoder but instead of directly processing the original data as input, some (unknown) pseudo-input data are defined and learned jointly with the VAE parameters. The sparsity level is tuned via a userdefined parameter, which governs the prior distribution of the mixing variable. Following a different strategy, [START_REF] Miao | On incorporating inductive biases into VAEs[END_REF] proposed a deterministic dimension selector function, so that some dimensions of the latent vector are deactivated via a point-wise multiplication. The amount of sparsity is then monitored via an entropy-based regularization term. Similarly, in [START_REF] Moran | Identifiable variational autoencoders via sparse decoding[END_REF] a stochastic per-feature masking variable is applied to the latent representation via point-wise multiplication to zero out some dimensions of the latent vector. A hierarchical Spike-and-Slab Lasso prior is assumed for the masking variable, consisting of two Laplace distributions with different parameters. A set of hyper parameters controls the desired amount of sparsity. Even though these approaches have shown promising performance, they often involve several user-defined hyperparameters and non-Gaussian distributions or discrete latent variables, which complicates the learning process.

In this paper, we adopt a different standpoint. We propose to structure the latent space via the usage of a sparsitypromoting dictionary model: each latent vector is assumed to have a sparse representation over a dictionary, i.e., it can be written as a weighted sum of only few columns of the dictionary. As a sparsity regularizer, we adopt a zero-mean Gaussian prior distribution with learnable variances: this is inspired by the relevance vector machines [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF], which has proven effective for promoting sparsity. Furthermore, in contrast to non-Gaussian priors, it does not complicate the learning procedure, and does not introduce any extra sparsity-adjusting parameter. We derive an efficient variational inference scheme to learn the parameters, with a negligible extra computational burden compared to a plain VAE. Our experimental results on speech generative modeling demonstrate the effectiveness of the proposed model in promoting sparsity while preserving the reconstructed speech quality.

The rest of the paper is organized as follows. Section 2 reviews the standard VAE model and inference. Section 3 presents the proposed approach and its application to speech spectrogram modeling. Experiments are conducted in Section 4. Finally, Section 5 concludes the paper.

Background on VAE

Generative modeling

Let s = {s1, . . . , sN } denote a set of training data with si ∈ R n .

The core idea in a generative modeling context is to encode the data via some latent variables, denoted z = {z1, . . . , zN } with zi ∈ R m , where usually m ≪ n. Then, the goal is to model the joint distribution p(s, z) = p(s|z) • p(z). To this end, the variational autoencoding framework [START_REF] Kingma | Auto-encoding variational bayes[END_REF] assumes some parametric forms for the generative distribution p(s|z), which is also referred to as the decoder, and for the prior distribution p(z), which are typically expressed as Gaussian distributions:

p θ (s|z) = N (µ θ (z), diag(σ 2 θ (z)), p(z) = N (0, I), (1)
where N (µ, Σ) denotes a Gaussian distribution with mean µ and covariance matrix Σ, I denotes the identity matrix of appropriate size, and the power is applied element-wise. Furthermore, µ θ (.) and σ θ (.) are some non-linear functions denoting the mean and standard deviation, which are implemented via some deep neural networks (DNNs). To learn the model parameters, i.e., θ, one would need to compute the posterior distribution of the latent codes, that is p θ (z|s), which is, unfortunately, intractable to compute. The VAE framework proposes to approximate this distribution with a parametric Gaussian distribution, called the encoder, as follows:

q ψ (z|s) = N (µ ψ (s), diag(σ 2 ψ (s)), (2)
where µ ψ and σ ψ are also implemented using some DNNs with parameters ψ.

Parameter estimation

The whole set of parameters Φ = {θ, ψ} is learned using the following variational procedure. Instead of directly optimizing the data log-likelihood log p θ (s), which is intractable, a lower bound L is targeted such that L(Φ; s) ≤ log p θ (s). To that end, the evidence lower bound (ELBO) of the data log-likelihood is considered, which is defined as:

L(Φ; s) = E q ψ (z|s) [log p θ (s|z)] -DKL(q ψ (z|s)∥p(z)), (3)
where DKL(q∥p) represents the Kullback-Leibler (KL) divergence between q and p. The KL divergence acts as a regularizer in the above formula, whereas the first term measures the reconstruction quality of the model. Then, the overall training phase consists of optimizing L(Φ; s) over Φ using a gradient-based optimizer together with the application of the reparametrization trick, which enables backpropagation through the encoder parameters [START_REF] Kingma | Auto-encoding variational bayes[END_REF].

Proposed framework 3.1. Model

The main idea underlying our proposal consists in structuring the latent space by considering a sparse dictionary model for each latent code zi as follows:

zi = Dai, ∀i, (4)
where D ∈ R m×k is a dictionary with unit-norm columns (to avoid scale ambiguity), and ai ∈ R k denotes a (sparse) representation. In this work, we consider that D is a fixed dictionary, even though it could also be learned jointly with the code vector ai. Note that one way to achieve a sparser representation is to consider an overcomplete dictionary, that is, such that k > m. The proposed model is illustrated in Fig. 1.

The motivation behind the model in ( 4) is to keep the reconstruction quality of VAEs intact, while at the same time promoting interpretability of the latent space by enforcing sparsity of the code vector ai. To impose such a constraint, one possible approach is to leverage some sparsity-promoting prior, such as based on the Laplace distribution [START_REF] Mohamed | Bayesian and L1 approaches for sparse unsupervised learning[END_REF]. However, such distributions involve non-smooth or complicated forms, which makes parameter inference and training more cumbersome. Instead, we propose to use the following zero-mean Gaussian prior:

p(ai; γi) = N (0, diag(γi)), (5)
in which γi is a vector of learnable variances. This prior has already proven successful in obtaining sparse solutions for linear regression and classification tasks [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF][START_REF] Wipf | Sparse bayesian learning for basis selection[END_REF]. Indeed, as discussed in [START_REF] Wipf | Sparse bayesian learning for basis selection[END_REF], if one considers an inverse Gamma hyperprior on γi in a Bayesian setting, then the prior distribution, which is given by:

p(ai) = p(ai|γi) • p(γi)dγi, (6)
becomes a Student's t-distribution, which has been shown to promote sparsity due to its sharp peak at zero [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF]. Even though we do not consider such a hyperprior explicitly, we still observe experimentally (see Section 4) that the model given by ( 4) and ( 6) promotes sparsity. We refer the reader to [START_REF] Wipf | Sparse bayesian learning for basis selection[END_REF] for more detailed discussions.

Inference

The inference phase follows the standard procedure of VAEs described in Section 2.2, with the main difference being that the encoder is now used to encode each si into ai, then compute zi via the dictionary model in [START_REF] Nguyen | Deep variational generative models for audio-visual speech separation[END_REF], and finally pass it to the decoder to reconstruct si. Therefore, zi only implicitly participates in the inference. Then, the cost function to optimize is:

L(Φ, γ; s) = E q ψ (a|s) [log p θ (s|z)] -DKL(q ψ (a|s)∥p(a; γ)), with zi = Dai ∀i. (7)
It should be noted that Φ are called amortized parameters, as they are shared among all the training data, while γ = {γ1, . . . , γN } are parameters specific to each data point.

We use an alternating minimization approach to update Φ and γ by alternately fixing one variable and optimizing the cost function over the other one. This corresponds to solving the following two sub-problems iteratively:

Update γ : γ ← argmin γ DKL(q ψ (a|s)∥p(a; γ)) (8a) Update Φ : Φ ← argmax Φ L(Φ, γ; s). (8b)
Since the KL divergence between two Gaussian distributions can be obtained in closed-form, it is straightforward to

▷ γ (b) -update: Using (9) ▷ Reparametrization: ▷ ϵ (b) ∼ N (0, I) ▷ a (b) = µ ψ (s (b) ) + σ ψ (s (b) ) ⊙ ϵ (b) ▷ z (b) -update: Using z (b) = Da (b)
▷ Φ -update: Using one-step gradient ascent on [START_REF] Mathieu | Disentangling disentanglement in variational autoencoders[END_REF] 5: End for 6: End while solve (8a), which yields the following update rule:

γi = E q ψ (a i |s) [a 2 i ] = µ 2 ψ (si) + σ 2 ψ (si), ∀i. (9)
Then, we resort to the reparametrization trick [START_REF] Kingma | Auto-encoding variational bayes[END_REF] to approximate the expectation in [START_REF] Asperti | Sparsity in variational autoencoders[END_REF], and thus reformulate subproblem (8b). Specifically, we approximate it using a single sample given by ai = µ ψ (si) + σ ψ (si) ⊙ ϵi, where ⊙ denotes the element-wise multiplication, and ϵi ∼ N (0, I). We then compute zi = Dai for each data sample i in order to finalize computation of the following approximated loss:

L(Φ, γ; s) = log p θ (s|z) -DKL(q ψ (a|s)∥p(a; γ)). (10)
Then, we optimize (10) by a single gradient ascent step, as in the standard VAE learning framework, which yields the update on Φ. The proposed procedure, called the sparsity-promoting dictionary model VAE (SDM-VAE) is summarized in Algorithm 1. Note that for practical implementation, the training data is split into mini-batches, and a stochastic gradient algorithm is applied sequentially to each mini-batch.

Application to speech modeling

In this paper, we evaluate the effectiveness of the proposed approach for the modeling of speech signals, that is, speech analysis-resynthesis [START_REF] Bie | A benchmark of dynamical variational autoencoders applied to speech spectrogram modeling[END_REF]. To that end, we first compute the short-time Fourier transform (STFT) of the time-domain signals, which yields the complex-valued data x = {x1, . . . , xN }, where N is the total number of time frames. In practice, the VAE processes and retrieves power spectrograms s = |x|2 . In terms of probabilistic modeling, this boils down to assuming that the data follows a circularly-symmetric zero-mean complex-valued Gaussian distribution, which is a common assumption when modeling speech signals with VAEs [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF].

Once the VAE model is trained on the power spectrograms, we perform speech analysis-resynthesis on the test data, similarly to [START_REF] Bie | A benchmark of dynamical variational autoencoders applied to speech spectrogram modeling[END_REF]. More precisely, once the input power spectrogram is encoded to obtain the posterior parameters, the latent's posterior mean is fed to the decoder to reconstruct the input spectrogram. Then, the magnitude spectrogram estimate ŝ1/2 is combined with the input STFT's phase ̸ x in order to obtain a complex-valued STFT: x = ŝ1/2 ⊙ exp(j ̸ x). Finally, timedomain estimates are retrieved by applying the inverse STFT to x. Note that, as such, the VAE ignores the phase information, which is only exploited for re-synthesizing time-domain signals. We leave modeling the phase in VAEs [START_REF] Nugraha | A deep generative model of speech complex spectrograms[END_REF] or extending the proposed method to complex-valued autoencoders, e.g., by overcoming the circularly-symmetric assumption [START_REF] Magron | Phase-dependent anisotropic gaussian model for audio source separation[END_REF], to future work.

As for the dictionary, we consider a discrete cosine transform (DCT) matrix with k atoms, that is, the r-th column of the DCT dictionary D is defined as follows [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]:

dr = [cos ((ℓ -1)πr/k)] m ℓ=1 , (11)
followed by mean subtraction and normalization to ensure unitnorm columns.

Experiments

In this section, we assess the potential of our proposed VAE model in terms of speech modeling. The code implementing the proposed VAE model is available online for reproducibility. 14.1. Protocol

Data

We use the speech data in the TCD-TIMIT corpus [START_REF] Harte | TCD-TIMIT: An audio-visual corpus of continuous speech[END_REF] for training and evaluating the model. It includes speech utterances from 56 English speakers with an Irish accent, uttering 98 different sentences, each with an approximate length of 5 seconds, and sampled at 16 kHz. The total speech data duration is about 8 hours. 39 speakers are used for training, 8 for validation, and the remaining 9 speakers for testing. The STFT parameters are as follows. The STFT is computed with a 1024 samples-long (64 ms) sine window, 75% overlap and no zero-padding, which yields STFT frames of length n = 513.

Methods & model architecture

As baseline methods, we compare the performance of a standard VAE 2 [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF], the VSC model3 [START_REF] Tonolini | Variational sparse coding[END_REF], and the proposed SDM-VAE model. All the VAE models follow the same simple encoder-decoder architecture as the one proposed in [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF] in order to focus our experiments on the latent space structure. More precisely, both the encoder and decoder are fully-connected networks using a single hidden layer with 128 nodes and hyperbolic tangent activation functions. The dimension of the latent space is m ∈ {32, 64}. The dictionary in SDM-VAE is a DCT matrix, as defined in [START_REF] Tonolini | Variational sparse coding[END_REF], which contains k ∈ {32, 64} atoms. These parameters values are chosen according to prior studies [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF][START_REF] Tonolini | Variational sparse coding[END_REF], where they have shown good performance. For comparison, we also consider the identity matrix as the dictionary D = I. In this setting, the model becomes similar to a classical VAE, except the latent prior now is zi ∼ N (0, diag(γi)) instead of the standard normal prior zi ∼ N (0, I).

Parameter settings & training

All the VAE variants are implemented in PyTorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF], and trained with mini-batch stochastic gradient descent using the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with a learning rate equal to 0.0001, and a batch size of 128. As stopping criterion we use early stopping on the validation set with a patience of 20 epochs, meaning that the training stops if the validation loss does not improve after 20 consecutive epochs. For the sparsity parameter of the VSC model, denoted α, we have experimented three different values: {0.05, 0.5, 0.9}, where a lower value corresponds to a sparser representation. To evaluate the compared methods in terms of reconstruction quality, we compute the short-term objective intelligibility (STOI) measure [START_REF] Taal | An algorithm for intelligibility prediction of time-frequency weighted noisy speech[END_REF] and the perceptual evaluation of speech quality (PESQ) score [START_REF] Rix | Perceptual evaluation of speech quality (PESQ)-a new method for speech quality assessment of telephone networks and codecs[END_REF]. These metrics respectively belong to the [-0.5, 4.5] and [0, 1] range, and higher is better. We also compute the Hoyer metric [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF] in order to evaluate the sparsity of the latent codes. For a vector z ∈ R m , the Hoyer metric is defined as follows:

Hoyer(z) = √ m -∥z∥1/∥z∥2 √ m -1 , (12)
where ∥.∥1 and ∥.∥2 respectively denote the ℓ1 and ℓ2 norms. The Hoyer metric ranges in [0, 1], and a higher score corresponds to a more sparse latent feature. It should be emphasized that for the proposed SDM-VAE model, the sparsity is measured on a, since it is the interpretable latent code in this case.

Results

From the results, averaged over all the test samples, and displayed in Table 1, we can draw several conclusions. First, we observe that for the classical VAE, the PESQ decreases when going from m = 32 to 64. Overall, the opposite behavior is observed for the sparse VAE models, for which increasing the latent space dimension does not sacrifice the reconstruction quality (except for SDM-VAE when using a DCT dictionary with k = 32 atoms, which is expected as the dictionary becomes undercomplete for m = 64). These results confirm a known advantage of sparsity as a regularizer for the model to avoid overfitting.

For VSC, we see that increasing sparsity (i.e., decreasing α) results in sacrifying the reconstruction quality in terms of PESQ when m = 32, although a different trend is observed for m = 64. Our proposed SDM-VAE model using the DCT dictionary is less sensitive to the sparsity level, since it exhibits more stable STOI values, and an increasing PESQ when the sparsity score increases.

We also observe that increasing the number of atoms in the dictionary results in an increased sparsity score. In particular, the highest sparsity score is obtained for m = 32 and when using an overcomplete dictionary (k > m). Nonetheless, this setting does not improve reconstruction quality compared to using a complete dictionary. As a result, using k = m would be an appropriate and simple guideline, since it maximizes reconstruction quality and avoid fine-tuning the dictionary size.

Interestingly, we observe that the best SDM-VAE results are obtained when using an identity dictionary rather than the DCT. This shows that a simple change in the prior structure compared to the baseline VAE (as detailed in Section 4.1.2) is an efficient way to both improve performance and promote sparsity, rather than resorting to more involved models such as in VSC. Nevertheless, the DCT dictionary remains the best performing option when m = 32. This motivates a future investigation on the design of more optimal dictionaries, or the joint learning of the dictionary along with other parameters based on dictionary learning algorithms from the literature [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF][START_REF] Sadeghi | Dictionary learning with low mutual coherence constraint[END_REF].

Finally, we remark that SDM-VAE yields the best results in terms of both reconstruction quality (PESQ and STOI) and sparsity (Hoyer score). This outlines the potential of our proposal, which does not rely on a complicated inference scheme or extra-parameter tuning, and exhibits a stable performance while enabling interpretability and avoiding overfitting.

Conclusions

In this work, we have proposed a novel approach for promoting sparsity in VAEs based on structuring the latent space using a sparse dictionary model. We derived a simple inference scheme, which does not require any fine-tuning of hyperparameters. Experiments conducted on speech signals modeling have demonstrated the potential of this technique compared to other existing sparse approaches. In particular, the proposed method improves speech reconstruction quality, and efficiently exploits sparsity to improve interpretability of the latent representation while reducing the risk of overfitting. Future work will focus on learning the dictionary along with other parameters, possibly with some mutual coherence constraint [START_REF] Sadeghi | Dictionary learning with low mutual coherence constraint[END_REF]. Extending the developed model to dynamical VAEs [START_REF] Bie | A benchmark of dynamical variational autoencoders applied to speech spectrogram modeling[END_REF], and exploiting it in speech enhancement [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF] and source separation applications [START_REF] Nguyen | Deep variational generative models for audio-visual speech separation[END_REF] are other future directions to pursue.

Figure 1 :

 1 Figure 1: Schematic diagram of the proposed VAE with a sparse dictionary model.

Table 1 :

 1 Reconstruction quality and sparsity measure for various VAE-based methods in terms of PESQ, STOI, and Hoyer scores. The results are averaged over the test set (the standard deviation for all methods is in the order of 0.01 for PESQ, and 0.001 for STOI and the Hoyer score).

	Dimension of z		m = 32			m = 64
			PESQ STOI Hoyer PESQ STOI Hoyer
	VAE [3]		3.29	0.85	0.40	3.26	0.85	0.56
		α = 0.05	3.00	0.81	0.57	3.25	0.84	0.51
	VSC [11]	α = 0.5	3.25	0.84	0.54	3.32	0.85	0.65
		α = 0.9	3.25	0.84	0.47	3.26	0.85	0.60
		I	3.33	0.86	0.64	3.45	0.87	0.73
	SDM-VAE	DCT (k = 32)	3.37	0.86	0.66	3.28	0.84	0.66
		DCT (k = 64)	3.32	0.86	0.87	3.33	0.86	0.76
	4.1.4. Evaluation

https://gitlab.inria.fr/smostafa/sdm-vae

https://gitlab.inria.fr/smostafa/avse-vae

https://github.com/Alfo5123/ Variational-Sparse-Coding

Acknowledgements

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).


Vidur Tuli

Antoine Claisse

Luca Messina

P A Burr

Solubility and vacancy-mediated inter-diffusion in the Zr-Nb-Cr system

To better understand the microstructural evolution of Cr-coated Zr alloy cladding, we investigate, using density functional theory calculations, the vacancy-mediated diffusion behaviour of Zr and Nb solutes in BCC-Cr as well as Zr and Cr solutes in BCC-Nb. The calculated vacancy formation and migration energies are in good agreement with available literature, as are the self-diffusivities. We show that in BCC-Cr both Nb and Zr are faster diffusers than Cr (vacancy-mediated self-diffusion). It was also found that both Zr and Nb segregate towards vacancy sinks in BCC-Cr at normal reactor operating temperatures, but at elevated temperatures their flux is expected to be in opposite directions. In BCC-Nb, we show that Cr is a slower diffuser than Nb self-diffusion, while Zr is faster; and both Zr and Cr are expected to decorate vacancy sinks in BCC-Nb at all temperatures relevant to reactor operation. A similar behaviour is likely to occur in β-Nb phase found in Zr-Nb alloys. The implications of these findings for Cr-coated Zr alloy cladding are discussed.

Introduction

Following the Fukushima Daiichi accident, there has been an interest in accident tolerant fuels (ATF) in the nuclear industry [START_REF] Burr | Accident Tolerant Fuel Concepts for Light Water Reactors tech[END_REF]. The goal of developing ATFs is to develop fuel that can tolerate loss of coolant accidents for considerably longer than existing fuels while at least maintaining current fuel performance during normal operation [START_REF] Terrani | Accident tolerant fuel cladding development: Promise, status, and challenges[END_REF][START_REF] Goldner | Overview of the accident-tolerant fuel development tech[END_REF][START_REF] Zinkle | Accident tolerant fuels for LWRs: A perspective[END_REF]. Cr-coated Zr alloys are a leading candidate for ATF cladding [START_REF] Terrani | Accident tolerant fuel cladding development: Promise, status, and challenges[END_REF][START_REF] Tang | Protective coatings on zirconium-based alloys as accident-Tolerant fuel (ATF) claddings[END_REF]. Cr is the coating of choice as it improves the cladding's oxidation and corrosion resistance for normal and accident conditions [START_REF] Bischoff | AREVA NP's enhanced accident-tolerant fuel developments: Focus on Cr-coated M5 cladding[END_REF][START_REF] Park | High temperature steam-oxidation behavior of arc ion plated Cr coatings for accident tolerant fuel claddings[END_REF]. Initial testing of first generation Cr-based coatings shows that it provides an exceptional balance between corrosion resistance, adhesion of the coating, neutron absorption, fretting resistance and resistance to steam oxidation [START_REF] Brachet | Early studies on Cr-Coated Zircaloy-4 as enhanced accident tolerant nuclear fuel claddings for light water reactors[END_REF].

The microstructural evolution at the interface between the Zr alloy substrate and the coating layer under reactor conditions is not yet fully understood. It is expected that a degree of intermixing will occur between Zr-alloy substrate and the Cr coating. This could be caused by either thermally activated equilibrium processes or by (ballistic) collision cascades. Following the collision cascades, the defects generated will also migrate, annihilate and cluster according to equilibrium thermodynamic and transport processes. In addition, the competition between dissolution and precipitation of secondary phases (e.g. ZrCr 2 [START_REF] Lumley | The stability of alloying additions in Zirconium[END_REF]) is also dictated by thermodynamics. As interdiffusion of these elements is a relevant part of such phenomena, we will look at the inter-diffusion of Zr and Cr. We also consider Nb as it is an important alloying addition in Zr alloy, and a potential coating material.

The solubility and diffusivity of alloying elements in α-Zr, including Cr and Nb, has been extensively studied. Previous work [START_REF] Lumley | The stability of alloying additions in Zirconium[END_REF][START_REF] Burr | From solid solution to cluster formation of Fe and Cr in α-Zr[END_REF][START_REF] Lu | Firstprinciples investigation on diffusion mechanism of alloying elements in dilute Zr alloys[END_REF] showed that Cr occupies both interstitial and substitutional sites in Zr. Both interstitial and vacancy mediated diffusivities of Cr in Zr have been calculated recently [START_REF] Jain | First-principles calculations of solute transport in zirconium: Vacancy-mediated diffusion with metastable states and interstitial diffusion[END_REF], showing that fast interstitial diffusion dominates under equilibrium conditions, in agreement with experimental measurements by Hood et al. showing that Cr is an extremely fast diffuser in Zr [START_REF] Hood | Chromium diffusion in α-zirconium, Zircaloy-2 and Zr-2.5Nb[END_REF][START_REF] Balart | 51Cr diffusion in α-Zr single crystals[END_REF]. Cr in α-Zr is trapped by excess vacancies and as a result it has also been suggested that the solubility of Cr in α-Zr increases with radiation damage [START_REF] Burr | From solid solution to cluster formation of Fe and Cr in α-Zr[END_REF]. In addition to being trapped by vacancies, Cr is also dragged towards vacancy sinks in Zr [START_REF] Jain | First-principles calculations of solute transport in zirconium: Vacancy-mediated diffusion with metastable states and interstitial diffusion[END_REF]. This suggests that in a radiation damaged matrix, the flux of Cr is predicted to be substantially different than in an undefective matrix [START_REF] Jain | First-principles calculations of solute transport in zirconium: Vacancy-mediated diffusion with metastable states and interstitial diffusion[END_REF].

Another common alloying element used in Zr-based claddings is Nb. This includes Zirlo [START_REF] Sabol | Development of a Cladding Alloy for High Burnup[END_REF], M5 [START_REF] Garner | LWR Fuel Performance/Top Fuel[END_REF], Zr-1%Nb and Zr-2.5%Nb. Nb additions are known to improve the cladding's corrosion properties [START_REF] Jeong | Influence of Nb concentration in the α-matrix on the corrosion behavior of Zr-xNb binary alloys[END_REF][START_REF] Jeong | Correlation between microstructure and corrosion behavior of Zr-Nb binary alloy[END_REF][START_REF] Kim | Corrosion and microstructural characteristics of Zr-Nb alloys with different Nb contents[END_REF]. When added above its solubility limit, Nb precipitates as β-Nb particles [START_REF] Park | Corrosion behavior of Zr alloys with a high Nb content[END_REF][START_REF] Banerjee | Precipitation in zirconiumniobium martensites[END_REF], and at higher concentrations still, it will form a two-phase β alloy (eg. Zr-2.5 wt% Nb alloys) [START_REF] Perovic | A study of the distribution of Nb and Fe in two-phase Zr-2.5 wt% Nb alloys[END_REF]; but in all Zr alloys Nb retains some solubility in the α phase. Even though previously reported to be a vacancy mediated diffuser [START_REF] Hood | Nb diffusion in single-crystal α-Zr[END_REF][START_REF] Tendler | Atomic size and fast diffusion of metallic impurities in zirconium[END_REF], it has recently been suggested that it is predominantly an interstitial diffuser with the octahedral site being the most stable interstitial site [START_REF] Lu | Firstprinciples investigation on diffusion mechanism of alloying elements in dilute Zr alloys[END_REF][START_REF] Zou | Hf diffusion in dilute Fe-free Zr(Nb) alloys[END_REF]. Interstitial diffusion mechanism of Nb in α-Zr, despite it preferentially occupying substitution sites, may be explained by the fact that the binding energy between solute and vacancy is not large enough to reduce the formation energy of the vacancy to levels competing with the low barrier for migration of interstitial Nb. Nb diffusion in Zr is reported to be 2-3 orders of magnitude faster than Zr self-diffusivity [START_REF] Hood | Nb diffusion in single-crystal α-Zr[END_REF]. Nb has also been observed to segregate heavily towards grain boundaries in Zirlo [START_REF] Hudson | Initial observation of grain boundary solute segregation in a zirconium alloy (ZIRLO) by three-dimensional atom probe[END_REF], suggesting diffusion towards vacancy sinks.

On the other hand, diffusion behaviour of Zr and Nb solutes in body-centered-cubic (BCC) Cr, and that of Zr and Cr solutes in BCC-Nb has hardly been studied. Here we use ab-initio atomic-scale simulations to shed light on the role of vacancies on the solubility, mobility and segregation of Zr and Nb solutes in BCC-Cr, and Cr and Zr solutes in BCC-Nb. We combine this information with existing knowledge of Cr and Nb solutes in hexagonal-close-packed (HCP) Zr, to paint a comprehensive picture of the equilibrium thermodynamic and kinetic processes that govern the inter-diffusion at the Zr-alloy/Cr-coating interface, thus shedding light on the microstructural evolution of this system.

Methodology

In the framework of thermodynamics of irreversible processes, the diffusion coefficients of each species in the system can be derived from the Onsager coefficients (L ss , L sv , L vv ). These coefficients are tensorial transport coefficients that relate solute and vacancy chemical potential gradients, µ s and µ v , to solute and vacancy fluxes, j s and j v [START_REF] Lidiard | Atomic transport in solids: Models and their parameterization[END_REF].

j s = -L ss ∇µ s -L sv ∇µ v
(1)

j v = -L vv ∇µ v -L sv ∇µ s (2)
In the dilute limit, the solute diffusivity D s is directly proportional to L ss [START_REF] Ramunni | Diffusion behavior in Nickel-Aluminum and Aluminum-Uranium diluted alloys[END_REF][START_REF] Mehrer | Diffusion in Solids: Fundamentals, Methods, Materials, Diffusion-Controlled Processes[END_REF] .

D s = k B T ΩL ss c s (3)
where k B is the Boltzmann constant, T is the temperature, Ω is the volume per atom and c s is the solute concentration. When only one diffusion mechanism is at play (or dominates), the solute diffusivity can also be expressed in Arrhenius form

D s = D s 0 e -Q k B T (4)
where D s 0 is the pre-exponential factor, containing all entropy terms, and Q is the activation energy required for diffusion.

We quantify the vacancy drag of solutes by calculating the drag ratio (DR)

DR = L sv /L ss (5)
A positive drag ratio means that the vacancy drags the solute along with it, on the other hand a negative ratio does not imply that the solute will segregate away from vacancy sinks. To ascertain this, one must consider the partial diffusion coefficient ratio (PDR), which compares solute to solvent diffusivity. In the dilute limit, the PDR is related to transport coefficients by the following expression [START_REF] Messina | Exact ab initio transport coefficients in bcc Fe-X (X=Cr, Cu, Mn, Ni, P, Si) dilute alloys[END_REF]:

d sv /d av = l as + l ss L 0 vv (6)
where, a refers to the solvent, s refers to the solute, v refers to the vacancy, d refers to the partial diffusion coefficient, and

l as = -l ss (1 + DR) (7)
l ss = L ss /c s (8)
Solute drag by vacancies, thus solute enrichment at vacancy sinks, happens when d sv is negative (d av is always positive), which in turn means the PDR is negative. In case of solute and vacancy diffusing in opposite directions (inverse Kirkendall effect), the PDR is always positive. However, solute enrichment at sinks may still occur when the solute and solvent diffuse in the same direction but the solute diffuses slower than the host matrix atoms or 0 < PDR < 1. Solute depletion at sinks happens when the PDR is positive and solute atoms diffuses faster than the matrix atoms, which means d sv >d av or PDR>1

The Onsager coefficients were calculated using both the Onsager code by Dallas Trinkle [START_REF] Trinkle | Diffusivity and derivatives for interstitial solutes: activation energy, volume, and elastodiffusion tensors[END_REF][START_REF] Trinkle | Automatic numerical evaluation of vacancy-mediated transport for arbitrary crystals: Onsager coefficients in the dilute limit using a Green function approach[END_REF] and the KineClue code by Schuler et al. [START_REF] Schuler | A kinetic cluster expansion code to compute transport coefficients beyond the dilute limit[END_REF]. Within numerical error, the two codes yielded consistent results despite differences in the way the cutoff for solute-vacancy correlations is calculated. The solutevacancy binding energies, the migration energy (E m ) for each jump, and the associated attempt frequencies were calculated using DFT simulations. We calculated the solute binding energies to vacancies up to the 10 th nearest neighbour position, and considered all interactions up to the 5 th nearest neighbour to build the network of transition rates. Fig 1 shows the vacancy mediated jump network for a solute in a BCC structure with interactions extending up to the 5 th nearest neighbour. In a dilute binary alloy, the vacancy jump in the dilute limit is denoted as ω 0 (not shown in the figure) and the solute-vacancy exchange jump is denoted as ω 2 .

Figure 1: Compendium of vacancy jumps in the vicinity of the solute atom in a BCC lattice. The numbers in the circles refer to the nearest neighbour position with respect to the solute s. Each ω represents a symmetrically unique jump, following the nomenclature from [START_REF] Messina | Exact ab initio transport coefficients in bcc Fe-X (X=Cr, Cu, Mn, Ni, P, Si) dilute alloys[END_REF] Cr at room temperatures and above is paramagnetic [START_REF] Bates | The magnetic properties of chromium[END_REF]. This disordered magnetic structure is computationally difficult to model [START_REF] Steneteg | Equation of state of paramagnetic CrN from ab initio molecular dynamics[END_REF][START_REF] Anisimov | Investigation of real materials with strong electronic correlations by the LDA+DMFT method[END_REF], thus we considered two bounding magnetic ordering of BCC-Cr: anti-ferromagnetic (AFM) and non-magnetic (no spin polarisation, NM). We neglected possible difference in binding and migration energy due to different spins. Considering that the binding energy and diffusivity results of the two bounding cases were not too different, we expect paramagnetic structure to exhibit similar results.

DFT simulations were carried out using the VASP code [START_REF] Kresse | Efficient iterative schemes for ab initio totalenergy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF] with the PBE exchange correlation functional [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] and a consistent plane wave cut off energy of 350 eV. BCC-Cr and BCC-Nb were modelled using 128 atom supercells with 4 × 4 × 4 k-point grid. Atoms were described with PAW pseudopotentials from the VASP 5.4 repository. Partial occupancy were treated with a first-order Methfessel-Paxton smearing function of width 0.1 eV. The equilibrium lattice parameters for BCC-Cr were found to be 2.84 Å(NM) and 2.86 Å(AFM). These are in agreement with values reported in literature -2.85 Å(NM) and 2.87 Å(AFM) [START_REF] Medvedeva | Structural properties, electronic structure, Fermi surface, and mechanical behavior of bcc Cr-Re alloys[END_REF]. The equilibrium lattice parameter for BCC-Nb was found to be 3.32 Å, which is in agreement with the value of 3.29 Å reported by Jiang et al. [START_REF] Jiang | First-principles study of binary bcc alloys using special quasirandom structures[END_REF]. All point defects were then modelled using the equilibrium lattice parameters.

Defect formation energies (E f ) were calculated using:

E f = E DFT d -E DFT p ± i µ(i) (9)
where E DFT d and E DFT p are the total energies from the defective and perfect DFT simulations and µ i is the chemical potential of all species i that are added or removed from the perfect crystal to form the defect. The chemical potential µ is calculated as the DFT energy per atom of the metallic elements in their ground state.

The solute-vacancy binding energy was calculated using the following equation

E b = E f sv -(E f s + E f v ) (10)
where E b is the binding energy, E f sv is the formation energy of the defect cluster containing a solute atom and a vacancy, E f s is the formation energy of the isolated solute defect and E f v is the formation energy of a vacancy.

The E m of each jump was calculated as the saddle point in the minimum energy paths between two stable defect sites. This was calculated using the climbing nudged elastic band method (cNEB) as implemented in the VTST-modified version of VASP [START_REF] Henkelman | Climbing image nudged elastic band method for finding saddle points and minimum energy paths[END_REF][START_REF] Henkelman | Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points[END_REF]. The cNEBs were simulated with a minimum of 5 images, and up to 9 images. For E m of jumps that do not involve the 1 st , 2 nd or 3 rd nearest neighbour we use the Final Initial State Energy (FISE) approximation [START_REF] Messina | Exact ab initio transport coefficients in bcc Fe-X (X=Cr, Cu, Mn, Ni, P, Si) dilute alloys[END_REF] as it yields sufficiently accurate results, see appendix for details.

The jump's attempt frequencies were calculated by dividing hopping atom's vibrational eigenfrequencies in the ground state configuration by the vibrational frequencies in the transition state, following Vineyard's expression [START_REF] Vineyard | Frequency factors and isotope effects in solid state rate processes[END_REF]

ν = 3 n=1 ν GS 2 n=1 ν T S (11)
The vibrational frequencies were computed with pertubation theory (phonon calculations) where only the hopping atom's degrees of freedom were considered. Vibrational frequencies of the entire supercell for each state are computationally demanding and result in high uncertainty due to the large number of frequencies in both the numerator and denominator [START_REF] Messina | Exact ab initio transport coefficients in bcc Fe-X (X=Cr, Cu, Mn, Ni, P, Si) dilute alloys[END_REF]. Therefore, we considered the vibrational frequencies of the hopping atom only, which has shown to yield accurate and reliable results in a range of systems [START_REF] Lu | Firstprinciples investigation on diffusion mechanism of alloying elements in dilute Zr alloys[END_REF][START_REF] Jain | First-principles calculations of solute transport in zirconium: Vacancy-mediated diffusion with metastable states and interstitial diffusion[END_REF][START_REF] Messina | Exact ab initio transport coefficients in bcc Fe-X (X=Cr, Cu, Mn, Ni, P, Si) dilute alloys[END_REF][START_REF] Burr | Formation and migration of point defects in tungsten carbide: Unveiling the sluggish bulk self-diffusivity of WC[END_REF].

Results and Discussion

The E f and E m in BCC-Cr and BCC-Nb are presented in Table 1, and are in excellent agreement with available literature. There is a significant discrepency between the vacancy migration energy we obtain in BCC-Nb and the results published by Nguyen-Manh [START_REF] Nguyen-Manh | Systematic groupspecific trends for point defects in bcc transition metals: An ab initio study[END_REF]. However, our results are in better agreement with more recent modelling results published by Ma et al. who reported a value of 0.64 eV [START_REF] Ma | Effect of stress on vacancy formation and migration in body-centered-cubic metals[END_REF] and experimental results published by Schultz who reported a value of 0.55 eV [START_REF] Schultz | 2.2.1 Special remarks: Datasheet from Landolt-Börnstein -Group III Condensed Matter • Volume 25[END_REF]. The self-diffusivity of Cr and Nb, calculated considering just the vacancy exchange mechanism, and using vacancy formation entropy values of 2.25 k B and 2.2 k B [START_REF] Burton | Vacancy-Formation Entropy in Cubic Metals[END_REF] respectively, are shown in Fig 2 . The self-diffusivity results are in excellent agreement with experimental bulk-diffusivity results [START_REF] Mundy | Self-diffusion in chromium[END_REF][START_REF] Einziger | Niobium self-diffusion[END_REF]. Literature available suggested that the self diffusion of Cr can be represented by mono-vacancy diffusion up to 1700-2000 K [START_REF] Park | Diffusion of tungsten in chromium: Experiments and atomistic modeling[END_REF][START_REF] Askill | Self-diffusion in chromium[END_REF][START_REF] Campbell | Positron trapping and self-diffusion activation energies in chromium[END_REF][START_REF] Peterson | Self-diffusion in pure metals[END_REF]. The experimental results reported for Nb self diffusion account for both monovacancy and divacancy diffusion. For monovacancy diffusion, Einzinger et al. [START_REF] Einziger | Niobium self-diffusion[END_REF] report the pre-exponential factor as 8 ± 3 ×10 -3 cm 2 /s and the activation energy required for diffusion as 3.64 eV. These are in reasonable agreement with the values of 0.3×10 -3 cm 2 /s and 3.14 eV respectively, obtained from our modelling work. The reported activation energy for divacancy diffusion is 4.54 eV, leading to different slope of the diffusivity results in Fig 2b .  Moreover, it is worth noting that our calculations are performed in the ground state (0K) and are extrapolated to the temperature range in question through harmonic approximation. The further we extrapolate the ground state results, the weaker the assumption of harmonic behaviour become. At high temperature the effect of thermal expansion, not accounted for here, becomes more important. Similarly anharmonic effects become more pronounced [START_REF] Lejaeghere | Error estimates for solid-state density-functional theory predictions: An overview by means of the ground-state elemental crystals[END_REF][START_REF] Moustafa | Accurate and precise ab initio anharmonic free-energy calculations for metallic crystals: Application to hcp Fe at high temperature and pressure[END_REF]. Furthermore, the existence of impurities in the samples could account for the remaining discrepancy [START_REF] Einziger | Diffusion of zirconium in niobium: The influence of fast diffusing impurities on the self-diffusion isotope effect[END_REF]. The formation energies of substitution and interstitial defects in BCC-Cr and BCC-Nb are summarised in Table 2. The table also contains defect formation energies in HCP-Zr to provide a wholistic picture of the Zr-Nb-Cr system. The Cr defect formation energies in HCP-Zr are from Burr et al. [START_REF] Burr | From solid solution to cluster formation of Fe and Cr in α-Zr[END_REF], whereas the Nb defect formation energies in HCP-Zr are from Xin et al. [START_REF] Xin | Point defect properties in hcp and bcc Zr with trace solute Nb revealed by ab initio calculations[END_REF]. It is evident that all solutes are preferentially accommodated as substitutional species in BCC-Cr and BCC-Nb. This accommodation as substitutional species is expected in BCC-Cr, given that the metallic radii of Zr (161 pm) and Nb (147 pm) are significantly larger than that of Cr (128 pm) [START_REF] Kim | Laves phase fields in Cr-Zr-Nb and Cr-Zr-Hf alloy systems[END_REF].

BCC-Cr

Considering the high formation energy for interstitial defects, modelling was not carried out for interstitial diffusion as it is not expected to contribute meaningfully to equilibrium diffusion (i.e. not driven by ballistic collisions). It was observed that dumbbells only form in the <110> configuration. This is also expected considering that the <110> direction is the most common dumbbell orientation in BCC metals [START_REF] Vainshtein | Structure of Crystals 330-399[END_REF][START_REF] Yip | Handbook of Materials Modeling[END_REF]. Conversely, crowdion defects with a longer range one-dimensional structure form along the <111> directions in BCC-Cr, while in BCC-Nb they form in the <100> and <111> direction. We expect the energy of these defects to be more sensitive to the choice of supercell size than other point defects [START_REF] Burr | Formation and migration of point defects in tungsten carbide: Unveiling the sluggish bulk self-diffusivity of WC[END_REF][START_REF] Derlet | Multiscale modeling of crowdion and vacancy defects in body-centered-cubic transition metals[END_REF].

Even though this study focuses on vacancy mediated diffusion, under irradiation damage interstitial diffusion becomes as important as vacancy-mediated diffusion. Thus, based on the interstitial sites observed, we can infer the possible diffusion pathways for interstitial diffusion of solutes in BCC-Cr and BCC-Nb, once the defects are formed as a consequence of collision cascades. Since the <111> crowdion site is the lowest energy interstitial site in both BCC-Cr and BCC-Nb, diffusion along the crowdion direction is likely to be fast. However, it is noted that as the crowdion moves forwards it effectively becomes a self-interstitial atom, with the solute atom moving only one site, until a further crowdion defect crosses the path of the solute. In BCC-Nb crowdion may form also along the <100> direction, but at considerably higher energy. Considering that the dumbbell configuration has only slightly higher energy compared to the <111> crowdion site in BCC-Cr (0.4-0.7 eV), we can expect this configuration to contribute significantly to the overall solute diffusion, possibly as an intermediate jump between two different <111> crowdions. Similar behaviour is expected for Zr in BCC-Nb. Cr in BCC-Nb, which is only stable on the octahedral site, is expected to have a simple (and isotropic) interstitial diffusion pathway from octahedral to octahedral site, until it is trapped by a vacancy or other defect. Conversely, in BCC-Cr, the octahedral site is not likely to contribute significantly to diffusion of Zr and Nb owing to its high formation energy. It is worth noting that the interstitial diffusion pathways discussed here are based on just the formation energies of the interstitial sites. However, the actual contribution of each site will depend on the migration barrier of interstitial defects as well (which we do not look at in this work).

Modelling was carried out for vacancy mediated diffusion of Zr and Nb in BCC-Cr and Zr and Cr in BCC-Nb. Fig 3 shows the binding energies of solutes to vacancies as a function of neighbour positions. Both Zr and Nb are found to be strongly bound to Cr vacancies in the 1 st nearest neighbour position, but exhibit near-zero binding at all other distances from the solute (Fig 3a). This distinct lack of attraction between vacancies and solutes past the first nearest neighbour suggests that the diffusion of Cr vacancies is not affected by the presence of Zr or Nb atoms as long as they stay further away than the 1 st nearest In BCC-Nb, Zr and Cr solutes exhibit markedly different interaction with Nb vacancies (Fig 3b). Zr exhibits the strongest binding to vacancies in the 1 st nearest neighbour position while Cr is most strongly bound to vacancies only in the 2 nd nearest neighbour position. Cr also exhibits repulsion to vacancies at most other distances. The binding energies of both the solutes exhibits stark symmetry and this may be due to the opposite difference in radii between the solutes and Nb (+14 pm for Zr and -19 pm for Cr [START_REF] Kim | Laves phase fields in Cr-Zr-Nb and Cr-Zr-Hf alloy systems[END_REF]).

Fig 4 and5 show the E m obtained using the cNEB method for diffusion in BCC-Cr and BCC-Nb respectively. The difference in E m along the migration pathways between different solutes decreases as the solute moves further away than the vacancy. This is expected as we get progressively further away from the solute, the vacancy jumps profiles approach those of the pure BCC metals.

The distance between the Zr solute initial and final position in the solute-vacancy jump (ω 2 ) is small in both BCC-Cr and BCC-Nb. This is expected as it is an solute and as a result it is placed closer to the vacancy at the start of the jump. The barrier for the (ω 2 ) jump is very low, which causes on the one-hand easier solute-vacancy exchange, but on the other hand a strong correlation effect that compensates for the low barrier. It is also worth noting that we do not observe any oversized solute atom mechanism as reported in [START_REF] Bocquet | Migration mechanism for oversized solutes in cubic lattices: The case of yttrium in iron[END_REF][START_REF] Lambert | Diffusion mechanisms of Mo contamination in Si[END_REF]. That is, in the present study, the solute-vacancy complex does not relax to a configuration where the solute sits in the middle of two halfvacancies.

Fig 6 shows the diffusivity of solutes in both BCC-Cr and BCC-Nb at equilibrium vacancy concentrations, and the D 0 and Q values for all solutes are provided in Table 3. Both Zr and Nb exhibit faster diffusion that Cr self-diffusion in BCC-Cr (Fig 6a). The effect of magnetic ordering of BCC-Cr on diffusivity is found to be small (∼ 10 cm 2 /s at high temperature ∼ 100 cm 2 /s at lower temperatures) with the NM structure providing the upper bound and AFM structure providing the lower bound.

It is also observed that Cr exhibits slower diffusion than Zr and Nb self-diffusion in BCC-Nb (Fig 6b). This can be explained by the repulsion of Cr with vacancies at the 1 st NN position and the high barrier for the solute-vacancy exchange jump.

Fig 6b also shows that the experimental diffusivity for Cr in BCC-Nb [START_REF] Pelleg | Diffusion of 51Cr into niobium single crystals[END_REF] is in excellent agreement with the current predictions (activation energies of 3.67 eV and 3.62 eV, respectively), while for Zr in BCC-Nb [START_REF] Einziger | Diffusion of zirconium in niobium: The influence of fast diffusing impurities on the self-diffusion isotope effect[END_REF] there is some discrepancy. The higher experimentally measured activation energy (3.93 eV, compared to 2.92 eV for the current study) could be due to a number of reasons. The first being, as the authors have acknowledged themselves, errors could have been introduced due to the challenge of observing radioactive impurity levels with the experimental techniques used. Another reason could be the possibility of divacancy mechanism playing a role in the diffusion like it did in Nb self-diffusion. Divacancy mechanism would have a higher activation energy as reported for Nb selfdiffusion. Finally, high temperature effects being ignored in modelling work presented here could lead to further discrepancy between the two set of results. However, experimentally Zr diffusivity is shown to be ∼ 1 order of magnitude higher than Nb self-diffusivity at high temperatures [START_REF] Einziger | Diffusion of zirconium in niobium: The influence of fast diffusing impurities on the self-diffusion isotope effect[END_REF][START_REF] Ablitzer | Diffusion of niobium, iron, cobalt, nickel and copper in niobium[END_REF], which is consistent with modelling results.

Table 4 shows a direct comparison between the Q values obtained from our work against the activation energy values obtained by the classical approach (Q classic ). In the classical approach:

Q classic = E f + E * m + E * b ( 12
)
where E f is the vacancy formation energy, E * m is the migration energy for the ω 2 jump and E * b is the solute binding energy to the vacancy in the 1 st NN position. It is clear that where ω 2 is the rate limiting step (for Cr in BCC-Nb), the agreement is excellent. However, when ω 2 is significantly faster than other jumps, correlation between reconfiguration, exchange and escape jumps becomes important. This highlights the importance of performing a thorough analysis of vacancy-solute interactions in cases where the classical approach would fail. Fig 7 shows the drag ratios and PDR of the solutes in BCC-Cr and BCC-Nb. The drag ratios of solutes in both BCC-Cr and BCC-Nb start off positive and large and reduce with increasing temperature. The decline in drag ratios is more dramatic in BCC-Nb at T > 800 K. In BCC-Cr the behaviour of Nb and that of Zr solutes are different and we observe a more rapid decline in the Nb drag ratio in the 600-1000 K temperature range. Since the Nb drag ratio becomes negative, it suggests that under elevated temperatures of 700-900 K (e.g. dry storage) Nb would diffuse away from vacancy sinks. The PDR of Nb in BCC-Cr, shown in Fig 7c, confirms that there would be depletion of Nb solute at vacancy sinks at these elevated temperatures [START_REF] Messina | Exact ab initio transport coefficients in bcc Fe-X (X=Cr, Cu, Mn, Ni, P, Si) dilute alloys[END_REF], as PDR>1. This means that the flux of Nb will be in opposite direction to that of Zr in BCC-Cr in the temperature range 700 K < T < 1200 K as far as the vacancy mechanism is concerned. The rate of change of the drag ratio of Nb as a function of temperature is also greater than that of Zr, suggesting that Nb redistribution is more sensitive to relatively small variations in temperature, which may then affect the performance of the coating after prolonged period of radiation damage and thermal cycling. Both solutes (Zr and Cr) also exhibit a positive drag ratio in BCC-Nb at all temperatures of concern. This suggests that, like Zr in BCC-Cr, both solutes are expected to decorate vacancy sinks in BCC-Nb (Fig 7b). The PDR, shown in Fig 7d, further confirms this finding as both Zr and Cr have a PDR < 1 in BCC-Nb. Since the PDR for Cr is really small at all temperatures shown (∼ -10 -07 to -10 -02 ), Cr drag by vacancies and enrichment at sinks, although active, will be negligible compared to Zr.

The behaviour of these solutes in BCC-Nb and BCC-Cr is in stark contrast to that of Cr and Nb in HCP-Zr. Both solutes are predominantly interstitial diffusers in HCP-Zr, and both are significantly faster than Zr self-diffusion. However, the diffusion behaviour of Cr in the HCP-Zr matrix is expected to change with increase in irradiation induced damage, leading to a con- siderably slower diffusion, and significant vacancy drag [START_REF] Jain | First-principles calculations of solute transport in zirconium: Vacancy-mediated diffusion with metastable states and interstitial diffusion[END_REF]. Considering the β-Nb phase in Zr-Nb alloys, or a potential Nb coating, we can infer the following behaviour from the current results. Firstly, Cr has a much lower solution energy in BCC-Nb (0.29 eV) compared to HCP-Zr (1.89 eV), see Table 1. This suggests that, in Zr-Nb alloys, Cr solutes would have a preference to dissolve in the β-Nb phase compared to the bulk of the alloy. However, we know that the ZrCr 2 is stable in both the C14 and C15 structures [START_REF] Lumley | The stability of alloying additions in Zirconium[END_REF][START_REF] Soubeyroux | Phase stability and neutron diffraction studies of Laves phases Zr(Cr1-xMx)2 with M = Mn, Fe, Co, Ni, Cu and 0 < x < 0.2 and their hydrides[END_REF] with a much lower formation enthalpy of -0.3 eV/formula unit [START_REF] Lumley | The stability of alloying additions in Zirconium[END_REF][START_REF] Sun | Ab initio calculation of the phase stability, mechanical properties and electronic structure of ZrCr2 Laves phase compounds[END_REF]. This means that the precipitation of ZrCr 2 phase, typically with significant Fe content in the form of Zr(Fe,Cr) 2 , will compete with the Cr dissolution in the β-Nb phase. We also expect both Zr and Cr solutes in the β-Nb phase to segregate towards vacancy sinks. However, looking at the Cr-vacancy repulsion at the 1 st nearest neighbour position, the higher ω 2 barrier, and the small PDR value, the transport of Cr in Nb will be limited. Radiation damage induced dissolution of β-Nb phase may lead to increased Considering the interface between Cr coating and the Zr based alloy, it is reasonable to assume that Zr or Nb atoms will be displaced in to the coating layer by a few Å due to interdiffusion or ballistic intermixing. However, the formation of ZrCr 2 phase, typically with Fe in the form of Zr(Fe,Cr) 2 has been extensively reported at the interface between Cr coating and Zr substrate, making it a more likely outcome of the interdiffusion and making it less likely for Zr and Nb solutes to displaced through the intermetallic layer of thickness ranging from 1-2 µm [START_REF] Brachet | Early studies on Cr-Coated Zircaloy-4 as enhanced accident tolerant nuclear fuel claddings for light water reactors[END_REF][START_REF] Zhong | Response of Cr and Cr-Al coatings on Zircaloy-2 to high temperature steam[END_REF][START_REF] Brachet | High temperature steam oxidation of chromium-coated zirconium-based alloys: Kinetics and process[END_REF]. We can still expect ballistic damage to cause redistribution of Zr and Nb solutes in to the Cr coating, either from the host cladding or the intermetallics formed at the interface between the Zr substrate and the Cr coating. At normal reactor operating temperatures (573-673 K), we expect both Nb and Zr solutes to segregate towards vacancy sinks in the Cr coating (e.g. grain boundaries and dislocation loops). However, under elevated temperatures (700 K < T < 1200 K) we expect the flux of Nb to be in the opposite direction to that of Zr, with Zr decorating vacancy sinks and Nb moving away from vacancy sinks. We suggest that under normal reactor operating temperature, vacancy sinks in the coating material will have an increased likelihood of intermetallic formation due to increased concentration of solute.

Conclusion

We use DFT calculations to better understand the solubility and inter-diffusion behaviour in the Zr-Nb-Cr system. The DFT calculation supports the following conclusions:

• In BCC-Cr, Zr and Nb solutes are preferentially accommodated as substitutional species. Same can be said for Zr and Cr solutes in BCC-Nb.

• Zr is strongly bound to vacancies in the 1 st nearest neighbour position in both BCC-Cr and BCC-Nb followed by much weaker binding at all other positions. Nb in BCC-Cr exhibits a similar behaviour. However, Cr in BCC-Nb exhibits strong binding to vacancies in the 2 nd nearest neighbour position and repulsion in most other positions.

• In BCC-Cr, both Zr and Nb are faster diffusers than Cr, with both solutes likely to decorate vacancy sinks up to temperatures of 700 K. Above 700 K, Zr and Nb solutes exhibit flux in opposite directions with Zr following vacancy flux, while Nb displays inverse Kirkendall effect.

• Cr is a slower diffuser than Zr and Nb self-diffusion in β-Nb. Both solutes are dragged by vacancies even at high temperature.

The current work lays the foundation of understanding the intermixing in Cr coated Zr alloys. The findings clarify the atomic scale transport behaviour of alloying elements in the alloy. Further work is required to look at the diffusivity of solutes in solid solutions found in these alloys. This is important as solid solutions, such as β-Zr and β-Nb phases, are quite commonly observed in Zr alloys.

Acknowledgements

We would like to thank Westinghouse Electric (Sweden) for providing financial support to enable work on this project. This research was supported by an Australian Government Research Training Program (RTP) Scholarship. This work was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government; the Multi-modal Australian Sci-enceS Imaging and Visualisation Environment (MASSIVE); the Pawsey Supercomputing Centre, which is supported by the Australian Government and the Government of Western Australia; and was enabled by Intersect Australia Limited. Learn more at www.intersect.org.au.

A. Appendix

Fig A .1 compares the energies obtained using the FISEestimated transition state and the NEB-simulated transition state in BCC-Cr with Zr as solute. It is evident that the difference becomes negligibly small for jumps further away than the third nearest neighbour from the solute. We expect other clusters to converge to dilute behaviour faster than [Zr Cr :V Cr ] BCC-Cr . This is due to the fact that Zr in Cr causes the largest compressive stress. This means that at long distances FISE for other clusters should be as good as, or better than, FISE for [Zr Cr :V Cr ] BCC-Cr

Figure 2 :

 2 Figure 2: Comparison between calculated and experimental self-diffusivity in BCC-Cr and BCC-Nb[START_REF] Mundy | Self-diffusion in chromium[END_REF][START_REF] Einziger | Niobium self-diffusion[END_REF]

Figure 3 :

 3 Figure 3: Binding energies of solutes to vacancies in BCC-Cr and BCC-Nb. Positive binding energy denotes repulsive interactions while negative binding energy denotes attractive interactions.

Figure 4 :Figure 5 :

 45 Figure 4: Vacancy migration energy with Nb and Zr solutes in BCC-Cr

Figure 6 :

 6 Figure 6: Diffusivity of Nb, Cr and Zr solutes in BCC-Cr and BCC-Nb

Figure 7 :Figure A. 1 :

 71 Figure 7: Drag ratios and partial diffusion coefficient ratio of solutes in BCC-Cr and BCC-Nb as a function of T

Table 1 :

 1 Vacancy formation and migration energy in BCC-Cr and BCC-Nb

			BCC-Cr	BCC-Nb
			NM	AFM
		Current	2.73 eV	2.96 eV	2.68 eV
	E f	Work
		Literature	2.64	3.02	2.6-3.0
		(mod-	eV [46]	eV [47]	eV [46,
		elling)			47, 49,
					50]
		Literature	-	-	2.70
		(experi-			eV [48]
		ment)
		Current	0.90 eV	1.10 eV	0.44 eV
	E m	Work
		Literature	0.91	1.11	0.90
		(mod-	eV [46]	eV [47]	eV [46],
		elling)			0.64
					eV [47]
		Literature	-	-	0.55
		(experi-			eV [48]
		ment)

Table 2 :

 2 Formation energies of extrinsic defects in BCC-Cr, BCC-Nb and HCP-Zr. Self-interstitial defects are not included as they are not relevant to the current study.

	Solute	Site	BCC-Cr	BCC-Nb	HCP-Zr
			NM	AFM
		Octahedral	10.13 eV 14.94 eV		2.56 eV [61]
	Nb	Dumbbell <110> Crowdion <111>	9.18 eV 8.47 eV	9.34 eV 8.86 eV	N/A	--
		Substitution	0.82 eV	0.88 eV		0.67 eV [61]
		Octahedral			3.57 eV	1.88 eV [10]
	Cr	Dumbbell <110> Crowdion	N/A	N/A	not stable not stable 2.06 eV [10] -
		Substitution			0.29 eV	1.89 eV [10]
		Octahedral	10.80 eV 10.84 eV	5.74 eV
	Zr	Dumbbell <110> Crowdion <111>	9.89 eV 9.31 eV	10.03 eV 9.68 eV	4.89 eV 4.61 eV	N/A
		Crowdion <100> not stable not stable	5.28 eV
		Substitution	1.634 eV	1.77 eV	0.29 eV
	neighbour configuration.

Table 3 :

 3 Diffusivity equation coefficients obtained by performing an Arrhenius fit to the diffusion coefficients shown in Fig 6

	Solute	Solvent	D 0 (cm 2 /s) Q (eV)
	Zr	Cr(AFM)	0.05	3.69
	Zr	Cr(NM)	0.04	3.34
	Nb	Cr(AFM)	0.07	3.89
	Nb	Cr(NM)	0.08	3.49
	Cr	Nb	0.06	3.67
	Zr	Nb	0.04	2.92

Table 4 :

 4 Comparison between reported Q and Q classic . Excellent agreement observed for Cr in BCC-Nb, where ω 2 is the rate limiting step.

	Solute Solvent Q	E f	E m	E b	Q classic
			(eV)	(eV)	(eV)	(eV)	(eV)
	Zr	Cr(AFM)3.69	2.96	0.05	-1.13	1.88
	Zr	Cr(NM) 3.34	2.73	0.05	-0.95	1.83
	Nb	Cr(AFM)3.89	2.96	0.19	-0.43	2.72
	Nb	Cr(NM) 3.49	2.73	0.28	-0.45	2.56
	Cr	Nb	3.67	2.69	0.84	0.14	3.67
	Zr	Nb	2.92	2.69	0.18	-0.45	2.42

concentration of Nb in the bulk of the alloy and consequently increased decoration of vacancy sinks in the α-Zr matrix.