# Condensate fraction in a 2D Bose gas measured across the Mott-insulator
transition

I. B. Spielman<EMAIL_ADDRESS>W. D. Phillips  J. V. Porto  Joint
Quantum Institute, National Institute of Standards and Technology, and
University of Maryland, Gaithersburg, Maryland, 20899, USA

###### Abstract

We realize a single-band 2D Bose-Hubbard system with Rb atoms in an optical
lattice and measure the condensate fraction as a function of lattice depth,
crossing from the superfluid to the Mott-insulating phase. We quantitatively
identify the location of the superfluid to normal transition by observing when
the condensed fraction vanishes. Our measurement agrees with recent quantum
Monte Carlo calculations for a finite-sized 2D system to within experimental
uncertainty.

Measurements of condensed matter systems realized by cold atoms in optical
lattices are now performed with sufficient accuracy to compare with ab-initio
calculations Gerbier et al. (2005a); Spielman et al. (2007); Mun et al.
(2007). Bosonic atoms in a lattice nearly perfectly realize the iconic Bose-
Hubbard (BH) Hamiltonian. Here, we study the system’s momentum distribution,
measure the condensate fraction, and accurately identify the transition point
from the the low temperature superfluid (SF) phase by identifying when the
condensate fraction vanishes.

The SF to Mott insulator (MI) transition can be accessed by changing the depth
of the optical potential Jaksch et al. (1998), and has been observed in 1D
Köhl et al. (2005), 2D Spielman et al. (2007) and 3D Greiner et al. (2002). A
range of studies have verified a detailed understanding of the MI phase in 2D
and 3D Gerbier et al. (2005a, b); Fölling et al. (2006); Spielman et al.
(2007). In contrast, the SF phase and the details of the transition to MI have
gone largely unstudied. Indeed, the only quantitative measurement locating the
transition is in 3D and is not in agreement with calculations Mun et al.
(2007). Here we focus specifically on the superfluid phase of a 2D system and
its transition to a normal state: we observe the expected increasing momentum
spread and vanishing condensate fraction as the system leaves the SF phase.
Our measured transition point agrees with the best available calculations
Wessel et al. (2004), thereby locating a point on the non-zero temperature 2D
BH phase diagram. Interestingly, the condensate fraction in our non-zero
temperature system vanishes more sharply than expected for a zero temperature
inhomogenous system, confirming that the superfluid regions are rapidly driven
normal as soon as an insulator appears Ho and Zhou (2007); Gerbier (2007).

The physics of interacting systems frequently depends spectacularly on
dimensionality: in 3D the SF is a conventional Bose-Einstein condensate (BEC);
in 2D, a Berezinskii-Kosterlitz-Thouless (BKT) SF; finally, in 1D there is no
true SF. In contrast, only the detailed properties of the MI phase depend on
dimensionality. In the $T>0$, 2D case studied here, the very existence of Bose
condensation is a consequence of the finite size of our trapped system. We
associate the presence of a bimodal momentum distribution with the SF phase,
and use fits to the distribution to identify the Bose-condensed fraction and
thereby measure the transition point between SF and normal.

At low temperature the transition from SF is to a normal state which crosses
over to a MI phase as the lattice depth increases Gerbier (2007); Ho and Zhou
(2007). As a result any $T>0$ measurement based on condensate fraction will
identify the SF to normal transition but be largely insensitive to the
subsequent crossover into the MI phase.

We study samples of ultra-cold rubidium atoms in a combined sinusoidal plus
harmonic potential. For atom occupancy per lattice site larger than unity
Greiner et al. (2002); Campbell et al. (2006), the low temperature SF phase
(shallow lattice) is expected to evolve into a structure composed of
alternating shells of SF and integer-occupied MI (deep lattice). As the
lattice deepens, each successive MI region appears and grows, as probed in
Ref. Fölling et al. (2006). At $T=0$ the amount of SF varies smoothly with
lattice depth giving no abrupt changes in the momentum distribution to
indicate a phase transition. In this work, we simplify the situation by
working near unit filling, where the only insulating phase is unit occupied
MI; thus, any observed signature can only be the transition from SF to normal.
Absent the lattice, recent experiments have shown that weak contact
interactions lead to a decrease in the 2D BEC transition temperature Kruger et
al. (2007). Lattice potentials increase the relative importance of
interactions; indeed, the onset of the MI phase corresponds to driving the
critical temperature to zero.

The BH model describes lattice-bosons that have a hopping matrix element $t$,
and an on-site interaction energy $U$. The physics of the BH model depends
only on $U/t$ Fisher et al. (1989). In an infinite, homogenous $T=0$ 2D
system, the transition from SF to MI occurs at the critical value $(U/t)_{\rm
c}\approx 16.5$ Krauth and Trivedi (1991); Elstner and Monien (1999); Wessel
et al. (2004); Kato et al. (2007). Remarkably, we observe a sharp transition
at $U/t=15.8(20)$ 111All uncertainties herein reflect the uncorrelated
combination of single-sigma statistical and systematic uncertainties. in our
$T>0$, finite-sized, harmonically trapped system.

Our data consists of images of atom density after sudden release and time-of-
flight (TOF), approximating the in-situ momentum distribution. Figure 1 shows
2D momentum distribution (right), and cross-sections through each distribution
(left). As evidenced by Fig. 1-a and -b, each diffraction order in the
momentum distributions consists of a narrow peak on a broad pedestal. Fitting
to a bimodal distribution (see below), we determine $f$, the fractional
contribution of the narrow component, and identify $f$ as the “condensate”
fraction. We associate images with non-negligible $f$ as being in the SF phase
Yi et al. (2007). We emphasize, however, that superfluidity is a transport
phenomena and cannot unambiguously be associated with features in the momentum
distribution Diener et al. (2007); Gerbier (2007); Ho and Zhou (2007). This
association is also imperfect at $T>0$ because in our 2D trapped system we
expect to observe a discernible “condensate” fraction even after the vortex-
pairs of a BKT SF unbind Hadzibabic et al. (2006), destroying the 2D SF. $f$
vanishes only when the resulting phase-fluctuating quasi-condensate vanishes
Kruger et al. (2007); Cladé and Helmerson (2008).

To characterize the transition from SF to normal, we extract two independent
quantities from an analysis of TOF images: $f$, and an “energy scale”
$\sigma$. We also measure a related quantity, the full width at half maximum
(FWHM) $\Gamma$ of the quasi-momentum distribution, which we compare to
theory. As the lattice depth is increased we find that $f$ vanishes
concurrently with a sudden increase in $\Gamma$, abrupt signatures that we
associate with the transition.

Figure 1: Momentum distributions and cross-sections at $U/t$ = 4(1), 8(1), and
20(2). Each row shows a single momentum distribution normalized by the total
atom number; the lines in the top-right panel indicate trajectories along
which four cross-sections were taken. The left panel shows the average of
these four sections (black solid line); the red-dashed lines denote the fit to
the bimodal distribution.

We produce nearly pure 3D ${}^{87}\rm{Rb}$ BECs with $N_{T}=1.2(4)\times
10^{5}$ atoms in the $\left|F=1,m_{F}=-1\right>$ state Spielman et al. (2007).
A pair of linearly polarized, $\lambda=820{\ {\rm nm}}$ laser beams forms a
$30(2){{E_{R}}}$ deep vertical optical lattice along $\hat{z}$ that divides
the 3D BEC into about $70$ 2D systems (turn-on time = $200{\ {\rm ms}}$). The
single photon recoil wave-vector and energy are ${{k_{R}}}=2\pi/\lambda$ and
${{E_{R}}}=\hbar^{2}{{k_{R}}}^{2}/2m=h\times 3.4{\ {\rm kHz}}$; $m$ is the
atomic mass and $h$ is Planck’s constant. The largest 2D system, containing
$\approx 3000$ atoms, has a chemical potential $\mu_{\rm 2D}=h\times
600(100){\ {\rm Hz}}$ and we measure a temperature $k_{B}T=k_{B}\times 33(4){\
{\rm nK}}=h\times 700(70){\ {\rm Hz}}$. Since the first vibrational spacing
$h\times 33(1){\ {\rm kHz}}\gg\mu_{\rm 2D},\ k_{B}T$, this system is well into
the 2D regime. In addition, a weaker, square 2D lattice in the
$\hat{x}$-$\hat{y}$ plane is produced by a second beam arranged in a folded-
retroreflected configuration Sebby-Strabley et al. (2006), linearly polarized
in the $\hat{x}$-$\hat{y}$ plane (turn-on-time = $100{\ {\rm ms}}$ 222The
beams for the two lattices originate from the same Titanium-Sapphire laser but
differ in frequency by about $160{\ {\rm MHz}}$.). The intensities of both
lattices follow exponentially increasing ramps, with $50{\ {\rm ms}}$ and
$25{\ {\rm ms}}$ time constants respectively, and reach their peak values
concurrently. These time-scales are chosen to be adiabatic with respect to
mean-field interactions, vibrational excitations, and tunneling within each 2D
system. The final depth of the $\hat{x}$-$\hat{y}$ lattice determines $U/t$
and ranges from $V=0$ to $25(2){{E_{R}}}$ 333The depth of the lattice along
$\hat{x}$ and $\hat{y}$ differ by 6%, and $V$ is the average.. The lattice
depths are calibrated by pulsing the lattice for $3{\ \mu{\rm s}}$ and
observing the resulting atom diffraction Ovchinnikov et al. (1998).

We calculate $U/t$ using a 2D band-stucture model and the $s$-wave scattering
length van Kempen et al. (2002). The uncertainty in $U/t$ stems from the
uncertainty in lattice depth 444The uncertainty in the $\hat{x}$-$\hat{y}$
lattice depth affects both $t$ and $U$, while vertical lattice uncertainties
affect only $U$. The $\pm 0.2\%$ uncertainty van Kempen et al. (2002) in the
${}^{87}\rm{Rb}$ $s-$wave scattering length is a negligible contribution to
the overall uncertainty.. The resulting uncertainty in $U/t$ is $\pm 10\%$.

Once both lattices are at their final intensity, the atomic system consists of
an array of 2D gasses each in a square lattice of depth $V$ and with a typical
density of 1 atom per lattice site. The atoms are held for $30{\ {\rm ms}}$,
and all confining potentials are abruptly removed (the lattice and magnetic
potentials turn off in $\lesssim 1{\ \mu{\rm s}}$ and $\simeq 300{\ \mu{\rm
s}}$, respectively). As a result, the initially confined states are projected
onto free particle states which expand for a $20.1{\ {\rm ms}}$ TOF 555Some of
the data were not taken under exactly these conditions: some had a $29.1{\
{\rm ms}}$ TOF and in others the lattice depth was rapidly (in $50{\ \mu{\rm
s}}$) increased to $30{{E_{R}}}$ before imaging (changing the single-site wave
functions, not the correlations which structure to the momentum distribution);
these differences do not affect the measurement, and are included in the
figures. , when they are detected by resonant absorption imaging. Apart from
effects of atomic interactions during expansion and the initial size of the
sample, initial momentum maps into final position, so each image approximates
the $\hat{x}$-$\hat{y}$ projection of the momentum distribution. We fit the
each momentum distribution to a simple function which describes the
distributions over the full range of $U/t$ studied here, with just three free
parameters.

First, we model the broad background as a thermal distribution of non-
interacting classical particles in a single 2D sinusoidal band where states
are labeled by their quasi-momentum $q_{x}$ and $q_{y}$,
$n(q_{x},q_{y})\propto\exp\left[2(\cos\pi q_{x}/{{k_{R}}}+\cos\pi
q_{x}/{{k_{R}}})/\sigma\right]$; this contributes two fitting parameters:
$\sigma$ and the non-condensed atom-number. In the shallow lattice limit,
$\sigma$ gives the temperature, $\sigma=k_{B}T/t$. This fit does not
distinguish atoms thermally occupying higher momentum states from atoms
occupying these states in the ground state wavefunction, i.e., from the
quantum depletion of the SF. In fact, $n(q_{x},q_{y})$ multiplied by a
suitable Wannier function, correctly describes the momentum distribution of
atoms in the MI phase to first order in $t/U$ where $\sigma$ is unconnected to
temperature, and is given by $\sigma=U/4t$. Our function fits the random phase
approximation (RPA) momentum distribution fairly well even as higher order
terms become important Sengupta and Dupuis (2005); Spielman et al. (2007).

The second portion of the momentum distribution consists of a narrow peak,
which we interpret as Bose-condensed atoms. We therefore take the narrow peak
to be the inverted parabola of a Thomas-Fermi profile (of fixed width for all
comparable data 666We do not allow the width of the condensate-peak to vary
with each fit; instead we first fit all of the SF data with the condensate-
width as a free parameter, and then repeat the fits with it held constant at
the average value: for $20.1{\ {\rm ms}}$ TOF we found $R_{\rm TF}=19(2){\
\mu{\rm m}}$, and for $29.1{\ {\rm ms}}$ TOF we found $R_{\rm TF}=26(2){\
\mu{\rm m}}$.), characterized by a single fitting parameter, condensed number.

The observed condensate peak width after TOF stems largely from initial system
size, not interaction effects during TOF or the initial momentum spread. Here
interactions during TOF are reduced due to rapid expansion along $\hat{z}$
after release from the tightly confining vertical lattice. Our analysis
further reduces these interaction effects by excluding data inside the $1^{\rm
st}$ Brillouin zone, with the highest density. This decreases the measured
FWHM of the peak from $30(1){\ \mu{\rm m}}$ to $22(1){\ \mu{\rm m}}$ and the
inferred momentum width from $0.26{{k_{R}}}$ to $0.21{{k_{R}}}$. Changing the
TOF from $20.1{\ {\rm ms}}$ to $29.1{\ {\rm ms}}$ only increased the FWHM from
$22(1){\ \mu{\rm m}}$ to $28(1){\ \mu{\rm m}}$ (decreasing the observed
momentum width from $0.21{{k_{R}}}$ to $0.17{{k_{R}}}$).

Figure 2: Condensed fraction $f$ and $\sigma$ vs. $V$ (bottom axis) or $U/t$
(top axis). The dots denote values determined from 2D fits to the full
momentum distribution: small dots result from one image and the large dots
indicate data averaged over about 20 separate images. The uncertainties are
their RMS variation, and are indicative of the single-image uncertainties. (a)
Condensate fraction. The red-dashed line is computed from our MFT model. (b)
fit parameter $\sigma$. There are two distinct regimes: at low $U/t$ it is
nearly constant (blue dashed line), from which we infer an initial temperature
$k_{B}T\approx 2t$; and at large $U/t$ it monotonically increases, consistent
with predictions of perturbation theory in the MI phase (red-dashed line).

Figure 2-a shows that as $V$ increases, $f$ vanishes at a critical value
$V_{\rm crit}$, while the total atom number remains constant. We verified that
this disappearance does not result from excessive irreversible heating of the
system by exceeding $V_{\rm crit}$, then lowering the lattice and observing a
condensed fraction Greiner et al. (2002).

To gain a qualitative understanding of the vanishing condensate fraction, we
performed a non-zero temperature mean-field theory (MFT) simulation of an
array of 2D BH systems in a 3D harmonic trap Sheshadri et al. (1993). To model
the non-zero temperature experimental system, we determine the entropy at
small $U/t$ that gives the observed $\approx 45\%$ condensate fraction, and
assume this entropy is unchanged as $V$ increases. The red-dashed line in Fig.
2-a shows the MFT condensate fraction vs. $V$ at constant entropy. Given that
$T=0$ MFT overestimates the transition ($(U/t)_{\rm MFT}=23.3$, compared to
$(U/t)_{\rm c}=16.5$ from more accurate calculations), the curve unexpectedly
lies on the data. MFT also gives $f$ as function of $U/t$ in units of
$(U/t)_{\rm c}$. We identify the transition point by fitting this function to
the data allowing $(U/t)_{\rm c}$ to vary, yielding $(U/t)_{\rm c}=15.8(20)$
(a lattice depth $V_{\rm crit}=9.0(5){{E_{R}}}$).

Figure 2-b displays $\sigma$ from the uncondensed background portion of the
distribution. At large $V$ we recover the behavior expected in the MI phase;
this measurement is equivalent to observations of the modulated momentum
distribution in the MI phase Gerbier et al. (2005a); Spielman et al. (2007).
At higher total atom number, our system would develop doubly and triply
occupied MI shells, expected to manifest as kinks in this curve. $\sigma$ is
monotonic with $V$, varying smoothly across $V_{\rm crit}$. This is in
agreement with RPA theory where the onset of superfluidity affects only states
near zero quasi-momentum. Figure 2-b shows that when $V\lesssim 4{{E_{R}}}$
($U/t\lesssim 3$), $k_{B}T/t\approx 2.0(3)$. Extrapolating to $V=0$ gives
$k_{B}T={{E_{R}}}\sigma/\pi^{2}\approx k_{B}\times 33{\ {\rm nK}}$ (valid when
$T\ll{{E_{R}}}$). This temperature is well below the $k_{B}\times 45{\ {\rm
nK}}$ expected for non-interacting particles in our 2D harmonic trap with
$f=0.45$, this reduction is similar to that observed in Ref. Kruger et al.
(2007), which focused on the critical temperature in interacting 2D atomic
systems with no 2D lattice.

Figure 3: Quasi-momentum width vs. $V$ (bottom axis) or $U/t$ (top axis). The
symbols denote the average FWHM of the quasi-momentum distribution along the
axes of highest symmetry (Top: averaged along $\hat{x}$+$\hat{y}$ and
$\hat{x}$-$\hat{y}$; Bottom: averaged along $\hat{x}$ and $\hat{y}$). The
small and large dots and uncertainties are as explained at Fig. 2. The red-
dashed line is the horizontally displaced RPA momentum width, as discussed in
the text, and the vertical grey line denotes the location of the SF-normal
transition identified from the sudden increase in $\Gamma$.

A related characterization of the system is the FWHM $\Gamma$ of the full
quasi-momentum distribution Greiner et al. (2002); Wessel et al. (2004); Kato
et al. (2007); Pollet and Troyer ; Yi et al. (2007). Figure 3 shows the width
of the 2D distributions (see Ref. Spielman et al. (2007)) as a function of
$V$. In the SF phase $\Gamma$ hardly depends on $V$ since the dominant feature
of the distribution is the condensate peak. $\Gamma$ only begins to change
very close to the SF to normal phase transition when the heights of the two
components of the bimodal distribution become comparable when the condensate
disappears, consistent with calculations in homogenous and trapped systems
Wessel et al. (2004). We calculate $\Gamma$ in the MI phase from the RPA
Sengupta and Dupuis (2005) quasi-momentum distribution which accurately
describes the large $U/t$ limit ($\sqrt{2}{{k_{R}}}$ along $\hat{x}+\hat{y}$
and ${{k_{R}}}$ along $\hat{x}$). In the RPA, $\Gamma$ is a function of $U/t$
in units of $(U/t)_{\rm c}$. The red dashed lines are fits to the measured
widths using two free parameters (joint between both panels): in the MI region
we use the RPA functional form with $(U/t)_{\rm c}$ as the first fit
parameter, and the constant width in the SF phase is the second. We obtain
$(U/t)_{\rm c}=16.7(20)$, in accord with the $(U/t)_{\rm c}=15.8(20)$ from our
fit to the condensate fraction.

We identify the point when the condensate fraction vanishes (and $\Gamma$
abruptly increases) with the onset of the SF to normal transition, i.e., when
a normal region begins to rapidly expand in our inhomogeneous system. (Our
measured visibility, computed as in Refs. Gerbier et al. (2005a, b), abruptly
drops from near unity at $U/t\approx 16$.) Increasingly accurate numerical
calculations give values of $(U/t)_{\rm c}$: $16.25(10)$ Wessel et al. (2004)
and $16.77$ Kato et al. (2007). Perhaps most relevant are QMC calculations
which include the effects of harmonic confinement; in this case, Wessel et al.
Wessel et al. (2004) find that a MI region first forms at $(U/t)_{\rm c}=17.2$
(the exact value of $(U/t)_{\rm c}$ depends on the details of the harmonic
potential). Both values lie within our experimental uncertainty.

The calculations Wessel et al. (2004); Kato et al. (2007) are at zero
temperature, and while they agree with our observed $(U/t)_{\rm c}$, they do
not predict a sudden increase in peak width or a vanishing condensate fraction
at $(U/t)_{\rm c}$. At $T=0$ and as $U/t$ increases past $(U/t)_{\rm c}$,
where an inhomogeneous system first develops a unit-occupied Mott core, the
shell of SF persists to large $U/t$. Thus, at $T=0$, $f$ drops rapidly at
$(U/t)_{\rm c}$, but does not vanish. Our system, however, is at small but
non-zero temperature, with a reduced condensate fraction of $\approx 45\%$ for
small $V$. Our MFT model shows that this temperatures quickly drives the SF
shells to the normal phase as $U/t$ increases past $(U/t)_{\rm c}$. As a
result the SF shells rapidly disappear as the normal, and then Mott regions
form (see Refs. Rey et al. (2006); Ho and Zhou (2007); Gerbier (2007)). That
this feature is seen in preliminary non-zero temperature QMC calculations
Pollet and Troyer  underscores the need for further non-zero temperature
calculations to compare with experiment.

This experiment constitutes the measurement of a single point of the non-zero
temperature 2D BH phase diagram. We expect future experiments will expand on
this result at different temperatures, densities, in different dimensions, and
in traps with more homogenous density distributions; new theory should aide in
the interpretation of these experiments.

We appreciate enlightening conversations with C. A. R. Sa de Melo, N. Trivedi,
L. Pollet, and C. J. Williams. We acknowledge the financial support of
ODNI/IARPA, and ONR; and I.B.S. thanks the NIST/NRC program.

## References

  * Gerbier et al. (2005a) F. Gerbier, A. Widera, S. Fölling, O. Mandel, T. Gericke, and I. Bloch, Physical Review Letters 95, 050404 (2005a).
  * Spielman et al. (2007) I. B. Spielman, W. D. Phillips, and J. V. Porto, Physical Review Letters 98, 080404 (2007).
  * Mun et al. (2007) J. Mun, P. Medley, G. K. Campbell, L. G. Marcassa, D. E. Pritchard, and W. Ketterle, Physical Review Letters 99, 150604 (2007).
  * Jaksch et al. (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Physical Review Letters 81, 3108 (1998).
  * Köhl et al. (2005) M. Köhl, H. Moritz, T. Stöferle, C. Schori, and T. Esslinger, Journal of Low Temperature Physics 138, 635 (2005).
  * Greiner et al. (2002) M. Greiner, O. Mandel, T. Esslinger, T. Hänsch, and I. Bloch, Nature 415, 39 (2002).
  * Gerbier et al. (2005b) F. Gerbier, A. Widera, S. Fölling, O. Mandel, T. Gericke, and I. Bloch, Physical Review A 72, 053606 (2005b).
  * Fölling et al. (2006) S. Fölling, A. Widera, T. Müller, F. Gerbier, and I. Bloch, Physical Review Letters 97, 060403 (2006).
  * Wessel et al. (2004) S. Wessel, F. Alet, M. Troyer, and G. G. Batrouni, Physical Review A 70, 053615 (2004).
  * Campbell et al. (2006) G. K. Campbell, J. Mun, M. Boyd, P. Medley, A. E. Leanhardt, L. Marcassa, D. E. Pritchard, and W. Ketterle, Science 313, 649 (2006).
  * Hadzibabic et al. (2006) Z. Hadzibabic, P. Krüger, M. Cheneau, B. Battelier, and J. Dalibard, Nature 441 (2006).
  * Cladé and Helmerson (2008) P. Cladé and K. Helmerson (2008), private communication.
  * Kruger et al. (2007) P. Kruger, Z. Hadzibabic, and J. Dalibard, Physical Review Letters 99, 040402 (2007).
  * van Kempen et al. (2002) E. G. M. van Kempen, S. J. J. M. F. Kokkelmans, D. J. Heinzen, and B. J. Verhaar, Phys. Rev. Lett. 88, 093201 (2002).
  * Fisher et al. (1989) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Physical Review B 40, 546 (1989).
  * Krauth and Trivedi (1991) W. Krauth and N. Trivedi, Europhysics Letters 14, 627 (1991).
  * Elstner and Monien (1999) N. Elstner and H. Monien, arXiv:cond-mat p. 9905367 (1999).
  * Kato et al. (2007) Y. Kato, N. Kawashima, and N. Trivedi (2007).
  * Yi et al. (2007) W. Yi, G.-D. Lin, and L.-M. Duan, arXiv p. 0705.4352v1 (2007).
  * Diener et al. (2007) R. B. Diener, Q. Zhou, H. Zhai, and T.-L. Ho, Physical Review Letters 98, 180404 (2007).
  * Gerbier (2007) F. Gerbier, Physical Review Letters 99, 120405 (2007).
  * Ho and Zhou (2007) T.-L. Ho and Q. Zhou, Physical Review Letters 99, 120404 (2007).
  * Sebby-Strabley et al. (2006) J. Sebby-Strabley, M. Anderlini, P. S. Jessen, and J. V. Porto, Phys. Rev. A 73, 033605 (2006).
  * Ovchinnikov et al. (1998) Y. B. Ovchinnikov, J. H. Müller, M. R. Doery, E. J. D. Vredenbregt, K. Helmerson, S. L. Rolston, and W. D. Phillips, Phys. Rev. Lett. 83, 284 (1998).
  * (25) L. Pollet and M. Troyer, private communication.
  * Sengupta and Dupuis (2005) K. Sengupta and N. Dupuis, Physical Review A 71, 033629 (2005).
  * Sheshadri et al. (1993) K. Sheshadri, H. R. Krishnamurthy, R. Pandit, and T. V. Ramakrishnan, Europhysics Letters 22, 257 (1993).
  * Rey et al. (2006) A. M. Rey, G. Pupillo, and J. V. Porto, Physical Review A 73, 023608 (2006).


IITM]Department of Metallurgical and Materials Engineering, Indian Institute
of Technology (IIT) Madras, Chennai 600036, India. Center for Atomistic
Modelling and Materials Design, Indian Institute of Technology (IIT) Madras,
Chennai 600036, India.

# 3-D Ti/TiN Interface

Prince Gollapalli  Varalakshmi Jalligampala  Kishor Peddapuvvala  Prajeet Oza
Satyesh Kumar Yadav<EMAIL_ADDRESS>[

## Keywords

3-D metal/ceramic interface, mechanical properties, density functional theory

## Abstract

Interface by definition is two-dimensional (2-D) as it separates 2 phases with
an abrupt change in structure and chemistry across the interface. The
interface between a metal and its nitride is expected to be atomically sharp,
as chemical gradation would require the creation of N vacancies in nitrides
and N interstitials in metal. Contrary to this belief, using first-principles
density functional theory (DFT), we establish that the chemically graded
Ti/TiN interface is thermodynamically preferred over the sharp interface. DFT
calculated N vacancy formation energy in TiN is 2.4 eV, and N interstitial in
Ti is -3.8 eV. Thus, diffusion of N from TiN to Ti by the formation of N
vacancy in TiN and N interstitial in Ti would reduce the internal energy of
the Ti-TiN heterostructure. We show that diffusion of N is thermodynamically
favorable till $\sim$23% of N has diffused from TiN to Ti, resulting in an
atomically chemically graded interface, which we refer to as a 3-D interface.
Experiments’ inability to identify a 3-D interface in Ti/TiN could be
attributed to limitations in identifying chemical composition and structure
with atomic-level resolution at interfaces. We define the sum of N vacancy
formation energy and N interstitial formation energy as driving-force, which
could be used as a convenient way to assess the possibility of forming a 3-D
interface in metal/ceramic heterostructures. We also show gradual variation in
lattice parameters and mechanical properties (like bulk modulus, shear
modulus, Young’s modulus, and hardness) across the Ti/TiN interface. 3-D
interfaces open a new way to control properties of metal/ceramic
heterostructures, in line with the already established advantage of gradation
at interfaces in micrometer length scale. For widely explored Ti/TiN
multilayer nano-heterostructures, the possibility of forming 3-D interface
could lead to enhanced wear and erosion resistance.

##  1 Introduction

#####

Interfaces play a critical role in the properties of nano-heterostructures
because of the large volume fraction 1. The structure of metal/ceramic
interface, like misfit dislocation network, and crystallographic orientation,
can significantly affect many properties 2, 3, 4, 5 such as dislocation
annihilation and nucleation, twinnability, recovery from radiation-induced
defects, thermal, optical, and electrical conductivity 6. These features can
represent interface structure in flat and sharp interfaces. Atomically
chemically graded interfaces present many more interface features like
continually varying mechanical properties, chemical, and atomic structure.
Chemical and microstructural gradation has been shown to improve wear and
erosion properties 7, 8, suppress crack formation in ceramic heterostructures,
and improve resistance to contact deformation and damage 9, 10. Similar
effects could be expected due to atomically chemically graded interface for
heterostructure.

#####

There have been several efforts to create atomically graded interfaces in
metal/ceramic heterostructures. All of them have focused on changing the ratio
of metal and anions that form ceramic. For example, Raveh et al. intended to
create a graded interface between Al and AlN by controlling the partial
pressure of N 11. Although it seems (based on AES depth profile) that they
have achieved a graded interface, there are no reports of stable compounds
with stoichiometry AlNx(x$<$1). The reason for such observation could be a
gradual change in the phase fraction of Al and AlN across the interface.
Limitations of experimental tools used to measure the atomic fraction of
materials across the interfaces invariably conclude that interfaces are
graded. Fe/MgO interface is also shown to be graded despite MgO being well
known to be a line compound; hence no intermixing of Fe in MgO should be
expected 12, 13. The depth profile of the sharp and rough interface between
metal and ceramic, as shown in Figure 1(a), would be interpreted as chemically
graded. Similarly, a sharp and flat interface, as shown in Figure 1(b), could
be construed as a chemically graded interface if the interface is tilted with
respect to the line of sight.

Figure 1: Schematic of the depth profile of metal/ceramic interfaces a) sharp
and rough b) sharp and flat. (c) Schematic of sharp Ti/TiN interface, with the
tendency of N to migrate from TiN to Ti. (d) Atomically chemically graded
interface formed either due to migration of N from TiN to Ti or formation of
compounds with a range of stoichiometry.

#####

From a purely thermodynamic consideration, atomic-level chemical gradation can
be achieved if 1) there are stable compounds with a range of stoichiometry or
2) mixing of anions or metal across the interface is thermodynamically
favorable. There are only a handful of ceramics with a range of stoichiometry,
for example, Ti-O system, Ti-C, Ti-N, and Fe-O, which we could find from
equilibrium phase diagrams 14, 15, 16, 17. Out of these, there has been
extensive effort to form a chemically graded Ti/TiN interface by varying
ratios of Ti and N during synthesis 18, 19. Ti/TiN heterostructure finds a
range of applications. Multilayer coating showed improved wear and corrosion
protection on metals compared to bare TiN 20, 21.

#####

We would like to particularly mention Watkins et al. 18 work which showed a
very smooth gradation in a mere distance of approximately 10 nm by controlling
N partial pressure while growing the Ti/TiN heterostructure. This smooth
gradation of N points to the possibility of a range of phases of Ti-N present
at the interface. Although from the equilibrium phase diagram, only a few
phases like solid solution of Ti and N with up to 20% atomic solubility, Ti2N,
and TiN at room temperature are known 22. Using DFT, we assess the relative
stability of phases reported in the equilibrium phase diagram along with a few
recently reported phases like Ti3N2, Ti4N3, and Ti6N5 23, 24. If these phases
with increasing N concentration are indeed stable, they can form part of the
graded interface.

#####

Starting from a sharp interface, forming a chemically graded interface between
metal and ceramics (oxide, nitride, or carbide) by diffusion of anions (O, N,
or C) across the interface is not expected as it involves the formation of
defects in both metal and ceramics. To our surprise, using first-principles
density functional theory (DFT), we show that Ti favors assimilation of N
atoms as interstitials within Ti over the formation of N vacancy in
stoichiometric TiN. This assimilation of N atoms in Ti suggests that a graded
interface is thermodynamically preferred over a sharp interface. The tendency
of N to migrate from TiN to Ti in the initially sharp Ti/TiN is in line with
the experimental observation by Susmita et al 25. Figure 1 (c) schematically
shows N migration from TiN to Ti by forming N vacancies in TiN and N
interstitials in Ti, starting from a sharp interface. If N migration continues
until it becomes unfavorable, the resulting interface would be chemically
graded, as shown in Figure 1 (d). We refer to this atomically chemically
graded interface as a 3-D interface.

#####

Using first-principles DFT calculations, we show that a range of Ti-N systems
is stable, explaining the formation of a chemically graded interface by
controlling the growth parameters. Then we show that it is thermodynamically
favorable for N to migrate from TiN to Ti, and approximately 23% of N can
diffuse from TiN to Ti. We also calculate misfit strains among various stable
phases likely to be found in the 3-D interface and their mechanical
properties.

##  2 Methodology

###  2.1 First-principles Calculations

#####

The first-principles electronic structure calculations were carried out within
the density functional theory (DFT) 26, 27 as implemented in Vienna Ab initio
Simulation Package (VASP) code 28, 29. The DFT calculations employed the
Perdew, Burke, and Ernzerhof (PBE) 30 generalized gradient approximation (GGA)
exchange-correlation functional and the projector-augmented wave (PAW) method
31. For all calculations, a plane wave kinetic energy cutoff of 520 eV was
used for the plane wave expansion of the wave functions to obtain highly
accurate forces. The Brillouin zone sampling was performed using Monkhorst-
Pack meshes 32 with a resolution of less than $2\pi\times 0.03~{}\AA^{-1}$.
All the structures were relaxed using conjugate gradient scheme, and the total
energies were converged to less than 1 meV per atom. The forces acting on each
atom are less than 0.01 eV/$\AA$. For Ti and N, 3p63d24s2 and 2s22p3 electrons
were considered as valence electrons, respectively. All the crystal structures
were drawn using VESTA 33 software.

##  3 Results and Discussion

###  3.1 Growth controlled formation of graded interface

#####

If graded interface forms by controlling N concentration, Ti-N phases with a
range of stoichiometries should exist. Gibbs free energy plot at the
temperature of interest would give accurate relative stability of various
phases. Assuming room temperature as the temperature of interest, we assess
the thermodynamic stability of various phases. Convex hull plot at zero kelvin
can give reliable results, as shown in Figure 2. The PV term contribution to
the Gibbs free energy is six orders of magnitude lower than internal energy,
and the contribution of entropy is also negligible. So we calculated
enthalpies of formation per atom for compounds and solid solutions as shown in
equation 1 and equation 2. The most stable Ti (hcp Ti) and N (N2 gas) phases
were considered as reference states in the hull plot.

$\Delta
H^{f}\left(Ti_{1-x}N_{x}\right)=E\left(Ti_{1-x}N_{x}\right)-\left(1-x\right)\
E\left(Ti\right)-\frac{1}{2}x\ E\left(N_{2}\right)$ (1)

$\Delta H^{f}\left(Ti-xN\right)=E\left(Ti-
xN\right)-E\left(Ti\right)-\frac{1}{2}x\ E\left(N_{2}\right)$ (2)

#####

Where $E(Ti_{1-x}N_{x})$ and $E(Ti-xN)$ are the DFT total energies of the
considered structures per atom in compounds, and solid solutions,
respectively, and $x$ is the atomic fraction of N. $E(Ti)$ and $E(N_{2})$ are
DFT total energies of Ti and N2.

Figure 2: Convex hull diagram for the Ti-N system. Solid circles represent
stable phases, and open circles represent phases above the hull.

#####

From the literature, we found the following thermodynamically stable Ti-N
phases: 1) solid solution of N in Ti (around 17-23 at.% N) 34, 35, 36, it is
represented as Ti-xN, where x is the atomic fraction of N, 2) compounds TiN,
Ti6N5, Ti4N3, Ti3N2, and Ti2N 23, 24, 34. Hence, we have considered the
following phases in our study: Ti, Ti-0.01N, Ti-0.14N, Ti-0.20N, Ti-0.25N,
Ti2N, Ti3N2, Ti4N3, Ti6N5, and TiN. Their detailed structural information like
space group, lattice parameters, and Wyckoff positions are given in Table 1.

Table 1: The space group, lattice parameters, and Wyckoff positions for various Ti-N phases. In Ti-xN, x is the percentage of N in Ti as an interstitial solid solution, where 0$\leq$ x$\leq$0.25. Structure | Lattice parameters | Wyckoff positions
---|---|---
a ($\AA$) | b ($\AA$) | c ($\AA$) | $\alpha$ | $\beta$ | $\gamma$ | Atom | Site | x | y | z
TiN (Fm$\bar{3}$m) | 4.252 | 4.252 | 4.252 | 90 | 90 | 90 | Ti | 4a | 0 | 0 | 0
cubic |  |  |  |  |  |  | N | 4b | 0.5 | 0.5 | 0.5
Ti6N5 (C2/m) | 5.214 | 9.015 | 8.517 | 90 | 144.8 | 90 | Ti | 8j | 0.487 | 0.824 | 0.746
monoclinic |  |  |  |  |  |  | Ti | 4i | 0.001 | 1 | 0.738
|  |  |  |  |  |  | N | 4h | 1 | 0.335 | 0.5
|  |  |  |  |  |  | N | 4g | 0.5 | 0.166 | 0
|  |  |  |  |  |  | N | 2c | 1 | 0 | 0.5
Ti4N3 (C2/m) | 9.934 | 3.018 | 10.302 | 90 | 150.3 | 90 | Ti | 4i | 0.266 | 0.5 | 0.633
monoclinic |  |  |  |  |  |  | Ti | 4i | 0.258 | 0 | 0.885
|  |  |  |  |  |  | N | 4i | 0.001 | 0.5 | 0.749
|  |  |  |  |  |  | N | 2a | 0 | 1 | 1
Ti3N2 (Immm) | 4.156 | 3.035 | 9.16 | 90 | 90 | 90 | Ti | 2d | 0.5 | 0 | 0.5
orthorhombic |  |  |  |  |  |  | Ti | 4j | 0 | 0.5 | 0.316
|  |  |  |  |  |  | N | 4i | 0 | 0 | 0.162
Ti2N (P42/mnm) | 4.956 | 4.956 | 3.036 | 90 | 90 | 90 | Ti | 4f | 0.203 | 0.797 | 0.5
tetragonal |  |  |  |  |  |  | N | 4f | 0 | 0 | 0
Ti-xN (P63/mmc) | 2.935 | 2.935 | 4.647 | 90 | 90 | 120 | Ti | 2c | 0.333 | 0.666 | 0.25
hexagonal |  |  |  |  |  |  | N(x) | 2a | 0 | 0 | 0
Ti (P63/mmc) | 2.935 | 2.935 | 4.647 | 90 | 90 | 120 | Ti | 2c | 0.333 | 0.666 | 0.25
hexagonal |  |  |  |  |  |  |  |  |  |  |

#####

To represent a solid solution of N in Ti within the periodic boundary
conditions, we took a primitive cell of hcp Ti having two atoms and repeated
4x4x3 along X, Y, and Z directions, forming a supercell with 96 Ti atoms. N is
more stable in octahedral voids than tetrahedral voids; the number of
octahedral voids is equal to the number of Ti atoms in the supercell. For
example, to create Ti-0.14N supercell, we took 96 Ti atoms supercell, and 16 N
atoms were placed in octahedral sites such that the average distance between N
atoms is maximum. A N atom was placed in a randomly chosen octahedral void in
Ti as the starting point, and the second void was chosen such that it is the
farthest from the first void. The third void was chosen such that it is the
farthest from the first two voids. The subsequent voids are filled similarly.
The same approach is used to construct other supercells representing Ti-xN
solid solutions.

#####

Figure 2 shows the convex hull plot of the Ti-N system. TiN, Ti6N5, Ti4N3,
Ti3N2, and Ti2N fall on the convex hull hence are stable, while solid
solutions of N in Ti (Ti-xN) are above the convex hull. Although phases above
the hull plot are considered unstable, the enthalpies of formation per atom
for Ti-0.01N, Ti-0.14N, Ti-0.20N, and Ti-0.25N phases are 0.002 eV, 0.024 eV,
0.033 eV, and 0.043 eV, respectively above the hull. These slight energy
differences suggest that they can be stable under suitable conditions. Being
part of the graded interface is one such condition that solid solutions could
be stable. If the solid solution splits into Ti and Ti2N, this will result in
an abrupt change in lattice strain, while the presence of Ti-N solid solutions
would provide a gradual change in lattice strain.

#####

The reason to obtain a smoothly varying chemically graded interface by
controlling N concentration during the growth could be two-fold: 1) phases
that are shown in the hull plot have smoothly varying N concentration, and 2)
all phases have a close-packed hexagonal arrangement of Ti atoms (Ti3N2,
Ti4N3, and Ti6N5 phases are variants of TiN with ordered N vacancies). While
forming a graded interface, hexagonal close-packed planes of various phases
can epitaxially align.

###  3.2 Diffusion controlled formation of graded interface

#####

If we assume that graded interface can indeed form by just migration of N from
TiN to Ti, Gibbs free energy of the Ti/TiN heterostructure should decrease on
the migration of N. Assuming N migration occurs in two steps: N vacancy
formation in TiN and N interstitial formation in Ti; an overall change in
Gibbs free energy, which we refer to as driving force ($\Delta G^{drive}$)
shown in equation 3 is the sum of Gibbs free energy of formation of N
vacancies ($\Delta G^{vac}$) given by equation 4 and Gibbs free energy of
formation of N interstitials ($\Delta G^{int}$) given by equation 5. During
the migration process, it is expected that more N vacancies and interstitials
are created near the interface; hence, vacancy and interstitial concentration
are not uniform. This non-uniformity of N migration makes it very difficult to
estimate the driving force. Thus we assume N vacancies are created uniformly
on the TiN side, and corresponding N interstitials are created uniformly in
Ti. N migration would stop only when the sum of the N vacancy formation energy
in Ti1-xNx (Ti in fcc lattice) and interstitial formation energy in Ti-xN (Ti
in hcp lattice) becomes positive. We also assume the interface structure would
not significantly alter the driving force ($\Delta G^{drive}$).

$\Delta G^{drive}=\Delta G^{vac}+\Delta G^{int}$ (3)

$\begin{split}\Delta
G^{vac}=E^{\left[Ti_{1-x}N_{x}\right]-y}-E^{Ti_{1-x}N_{x}}+\frac{1}{2}\mu_{N_{2}}+P\left(V^{\left[Ti_{1-x}N_{x}\right]-y}-V^{Ti_{1-x}N_{x}}\right)\\\
-T\left(S_{vib}^{\left[Ti_{1-x}N_{x}\right]-y}-S_{vib}^{Ti_{1-x}N_{x}}+S_{con}^{\left[Ti_{1-x}N_{x}\right]-y}-S_{con}^{Ti_{1-x}N_{x}}\right)\end{split}$
(4)

$\begin{split}\Delta G^{int}=E^{\left[Ti-xN\right]+y}-E^{Ti-
xN}-\frac{1}{2}\mu_{N_{2}}+P\left(V^{\left[Ti-xN\right]+y}-V^{Ti-xN}\right)\\\
-T\left(S_{vib}^{\left[Ti-xN\right]+y}-S_{vib}^{Ti-xN}+S_{con}^{\left[Ti-
xN\right]+y}-S_{con}^{Ti-xN}\right)\end{split}$ (5)

#####

Where $x$ is the atomic fraction of N in the Ti-N system, and $y$ is the
number of N atoms such that the percentage of N vacancies or interstitials is
in the order of ppm. $Ti_{1-x}N_{x}$ and $Ti-xN$ represent pristine compounds
and solid solutions, respectively. $[Ti_{1-x}N_{x}]-y$ and $[Ti-xN]+y$
represent systems with vacancies and interstitials, respectively. $V$,
$S_{vib}$, $S_{con}$ are the volume, vibrational entropy, and configurational
entropy, respectively. $E^{Ti_{1-x}N_{x}}$ and $E^{Ti-xN}$ are the internal
energies of compounds (TiN, Ti6N5, Ti4N3, Ti3N2, and Ti2N) and solid solutions
of N in Ti. $E^{\left[Ti_{1-x}N_{x}\right]-y}$ and $E^{\left[Ti-xN\right]+y}$
are the internal energies of bulk with N vacancies and N interstitials.
$\mu_{N_{2}}$ is the chemical potential of nitrogen, and $P$ is pressure.

#####

Adding up equations 4 and 5, the term $\frac{1}{2}\mu_{N_{2}}$ will cancel
out. At ambient pressure, the terms,
$P\left(V^{\left[Ti_{1-x}N_{x}\right]-y}-V^{Ti_{1-x}N_{x}}\right)$ from
equation 4 and $P\left(V^{\left[Ti-xN\right]+y}-V^{Ti-xN}\right)$ from
equation 5, are insignificant as the difference in the volume of condensed
matter due to creation of vacancies or interstitials are negligible.
Therefore, we can safely ignore this contribution. Vibrational entropies of
$Ti_{1-x}N_{x}$ with vacancy ($S_{vib}^{\left[Ti_{1-x}N_{x}\right]-y}$) and
without vacancy ($S_{vib}^{Ti_{1-x}N_{x}}$) are expected to be similar.
Similarly, the difference in vibrational entropies between system $Ti-xN$ with
interstitial ($S_{vib}^{\left[Ti-xN\right]+y}$) and without interstitial
($S_{vib}^{Ti-xN}$) should be infinitesimal. Hence, vibrational entropy
contribution could be ignored.

#####

The configurational entropy, $S_{con}^{Ti_{1-x}N_{x}}$ is zero as they are
ordered compounds. Contribution of configurational entropy to Gibbs free
energy $S_{con}^{\left[Ti_{1-x}N_{x}\right]-y}$ is positive, owing to the
creation of vacancy in an ordered compound. $S_{con}^{\left[Ti-
xN\right]+y}-S_{con}^{Ti-xN}$ is positive because there are more N
interstitials in $[Ti-xN]+y$ than $Ti-xN$. Hence the overall contribution of
configurational entropy to Gibbs free energy is negative, thus lowering the
$\Delta G^{drive}$. The contribution of configurational entropy is expected to
be small, and it is not easy to accurately account for a range of systems
considered, so we ignore its contribution. Thus, the driving force $\Delta
G^{drive}$ to form a chemically graded interface can be written as shown in
equation 6.

$\Delta G^{drive}=E^{\left[Ti_{1-x}N_{x}\right]-y}\ -E^{Ti_{1-x}N_{x}}+E^{[Ti-
xN]+y}-E^{Ti-xN}$ (6)

#####

To calculate $\Delta G^{drive}$, only the internal energies of defected and
pristine structures would suffice and can be reliably calculated using DFT. At
0 K, vacancy formation energy and interstitial formation energy can be
calculated according to equations 7 and 8, respectively.

$\Delta
E^{vac}=E^{\left[Ti_{1-x}N_{x}\right]-1}-E^{Ti_{1-x}N_{x}}+\frac{1}{2}E\left(N_{2}\right)$
(7)

${\Delta E}^{int}=E^{\left[Ti-xN\right]+1}-E^{Ti-
xN}-\frac{1}{2}E\left(N_{2}\right)$ (8)

#####

Where $E^{\left[Ti_{1-x}N_{x}\right]-1}$ and $E^{[Ti-xN]+1}$ are the DFT-
calculated total energies of the defect supercell containing one N vacancy and
one N interstitial, respectively. $E^{Ti_{1-x}N_{x}}$ and $E^{Ti-xN}$ are the
DFT-calculated total energies of the pristine supercells with compounds and
solid solutions. $E\left(N_{2}\right)$ is the DFT total energy of N2 molecule.
The sum of vacancy formation energy and interstitial formation energy is equal
to $\Delta G^{drive}$.

#####

Vacancy and interstitial formation energies are defined for a small vacancy or
interstitial concentration in a few ppm. Thus supercell containing a few
million atoms is required. Since it is not practical to model millions of
atoms using density functional theory, we have modeled at least 1.04% N
vacancies or interstitials and confirmed that further reducing the percentage
of vacancies and interstitials does not significantly affect their formation
energies. Even if the concentration of defects is reduced by half, formation
energies change by a maximum of 100 meV. We would also like to point out that
supercells considered are large, and the images of vacancies or interstitials
are at least 10 $\AA$ apart.

#####

In principle, one should calculate vacancy formation energy in all Ti1-xNx and
interstitial formation energy in all Ti-xN. Given computational limitation, we
only calculate for a few stable compounds identified from hull plot: TiN,
Ti6N5, Ti4N3, and Ti3N2. We correspondingly chose Ti, Ti-0.14N, Ti-0.20N, and
Ti-0.25N such that the fraction of N vacant/filled sites in TiN/Ti,
Ti6N5/Ti-0.14N, Ti4N3/Ti-0.20N, Ti3N2/Ti-0.25N are the same. The fraction of N
vacant sites is the ratio of the unfilled N sites to the available N sites and
the fraction of N-filled sites is the ratio of the filled N sites to the
available N sites. All systems were modeled using supercells containing 96 Ti
atoms; thus, 96 octahedral sites are available for N to occupy.

Table 2: Octahedral N vacant and filled sites in the Ti-N systems. Compound | Pristine | Defect | Vacant | % Vacant | Solid | Pristine | Defect | Filled | % Filled
---|---|---|---|---|---|---|---|---|---
Ti | N | Ti | N | N sites | N sites | Solution | Ti | N | Ti | N | N sites | N sites
TiN | 96 | 96 | 96 | 95 | 1 | 1.04 | Ti | 96 | 0 | 96 | 1 | 1 | 1.04
Ti6N5 | 96 | 80 | 96 | 79 | 17 | 17.71 | Ti-0.14N | 96 | 16 | 96 | 17 | 17 | 17.71
Ti4N3 | 96 | 72 | 96 | 71 | 25 | 26.04 | Ti-0.20N | 96 | 24 | 96 | 25 | 25 | 26.04
Ti3N2 | 96 | 64 | 96 | 63 | 33 | 34.38 | Ti-0.25N | 96 | 32 | 96 | 33 | 33 | 34.38
Figure 3: Vacancy/interstitial formation energy as a function of atomic % N
vacant/filled sites. Red squares represent vacancy formation energies ($\Delta
E^{vac}$), and blue circles represent interstitial formation energies ($\Delta
E^{int}$) of Ti-xN phases, where x is the concentration of N. The sum of
vacancy and interstitial formation energies is the driving force ($\Delta
G^{drive}$) as indicated in black upper triangles.

#####

Table 2 shows the number of Ti and N atoms in the pristine and defect
supercells. Defect supercells were constructed by removing/adding one N atom
from pristine supercells. The vacancy formation energy of Ti-N compounds could
depend on the choice of Wyckoff positions for removal of N. So we removed N
from all available Wyckoff positions, and vacancy formation energy shown in
figure 3 corresponds to the lowest formation energy because the first N to
come out will be the one which is easiest to remove. In Ti6N5, N occupies 2c,
4g, and 4h Wyckoff positions, and the N removed from the 4g position has the
lowest energy value of 3.22 eV (2c has 3.25 eV and 4h has 3.29 eV). Similarly,
in Ti4N3, N occupies 2a and 4i Wyckoff positions, and the N removed from the
2a position has the lowest energy value of 3.43 eV (4i has 3.65 eV).

#####

Figure 3 shows vacancy and interstitial formation energies and driving force.
Vacancy formation energy increases as more and more N is removed from FCC Ti.
In contrast, interstitial formation energy does not change significantly as N
is added to HCP Ti, at least for the range considered. The driving force for N
diffusion from FCC Ti to HCP Ti is most negative (-1.39 eV/atom) when 1.04% of
N atoms in TiN is removed, and the corresponding percentage Ti octahedral
sites are filled by N. Ti6N5 and Ti-0.14N correspond to 17.71% N vacant and
filled sites in FCC and HCP Ti, respectively. The driving force for N atom to
diffuse from FCC lattice to HCP lattice now reduces to -0.58 eV/atom. As more
N diffuses to HCP Ti, the driving force further reduces to -0.35 eV/atom and
becomes positive (0.21 eV/atom) for 26.04% and 34.38% N atoms changing sides,
respectively. As the diffusion of N would be thermodynamically unfavorable
when the driving force becomes positive, it can be inferred from Figure 3 that
a maximum of 23 at.% of N can migrate.

#####

Our calculations suggest that starting from a sharp interface by creating N
vacancies in TiN and their diffusion into Ti is thermodynamically favorable,
and the driving force decreases as more N diffuse. A sufficient amount of N
vacancies in TiN near the interface transforms the region to Ti6N5 and then
progressively to Ti4N3; on the other side, more and more N interstitials lead
to the formation of Ti-xN (Ti-0.20N and Ti-0.14N) phases. Although our results
suggest that migration should happen until all TiN transforms to the Ti4N3
phase, the kinetics of migration of N would hinder in completion of the
transformation, and such transformation is expected to be only near the
interface. This should result in: 1) on TiN side phases with the most amount
of N vacancies would be closest to the interface like Ti4N3, followed by Ti6N5
and TiN; 2) on Ti side N rich solid solutions would be closest to the
interface and N concentration would decrease away from the interface. The
interplay of driving force and kinetics thus leads to forming a chemically
graded interface, as shown in Figure 4. We want to point out that we only
discuss phases that we calculated to estimate driving force, but we believe
there will be other phases that have N concentration between the phases
examined. They can form by forming vacancies in ordered structures like TiN,
Ti6N5, and Ti4N3 or varying N interstitial concentrations in solid solutions.

Figure 4: Schematic representation of diffused interface and atomic
representation of the diffused Ti/TiN interface composing of Ti-N phases
considered in the calculation of driving force. Ti (0001) plane is parallel to
TiN (111) plane, and both are parallel to close-packed planes of all other
phases. The average Ti-Ti bond length in the close-packed plane of all phases
is assumed to be the same in the atomic representation.

#####

As Ti4N3 and Ti6N5 have rocksalt structures with ordered N vacancy, the
proposed diffused interface can occur by diffusion of N from TiN side to Ti,
without changing the atomic position of Ti. As the diffusion of N controls the
formation of a diffused interface, temperature and exposure time will dictate
the width of the diffused interface. Ti generally grows along $<$0001$>$
direction and TiN grows along $<$111$>$ direction, thus
${\left(111\right)_{TiN}\ ||\ \left(0001\right)}_{Ti}$. Since the atomic
arrangement of Ti does not change on the formation of the diffused interface,
close-packed planes of other phases will also be parallel to
$\left(111\right)_{TiN}$ and $\left(0001\right)_{Ti}$.

####  3.2.1 Misfit among Ti-N phases in 3-D Ti/TiN interface

Figure 5: Atomic structure of close-packed (0001) plane of HCP Ti and (111)
plane of rocksalt TiN. Blue spheres represent Ti atoms. AB and AC are the Ti-
Ti bond lengths along
$\left[10\bar{1}0\right]_{Ti}/\left[11\bar{2}\right]_{TiN}$ and
$\left[1\bar{2}10\right]_{Ti}/\left[1\bar{1}0\right]_{TiN}$ directions,
respectively. (b) Variation of Ti-Ti bond lengths along AB and AC directions.
(c)-(f) represents the atomic structure of close-packed planes and shows two
orthonormal lattice vectors used to calculate average AB and AC distances.

#####

In a sharp Ti/TiN interface, misfit is accommodated in one plane near the
interface if not at the interface by forming misfit dislocations. The
atomistic structure of close-packed planes of hcp Ti and rocksalt TiN is shown
in Figure 5a. In close-packed planes of Ti and TiN,
$\left[1\bar{1}0\right]_{TiN}\ ||\ \left[11\bar{2}0\right]_{Ti}$ and
$\left[11\bar{2}\right]_{TiN}\ ||\ \left[10\bar{1}0\right]_{Ti}$. AB and AC
are Ti-Ti bond lengths along
$\left[1\bar{1}0\right]_{TiN}/\left[11\bar{2}0\right]_{Ti}$ and
$\left[11\bar{2}\right]_{TiN}/\left[10\bar{1}0\right]_{Ti}$, respectively. The
calculated percentage difference in the Ti-Ti bond length in Ti and TiN along
the two orthonormal directions is 2.34%, and percentage difference along these
directions has been used to assess misfit strain across the Ti/TiN interface
37.

#####

All the phases mentioned in Figure 4 are expected to be a part of the
chemically graded interface with the hexagonal arrangement (not necessarily
regular hexagon) of Ti atoms, so equivalent Ti-Ti bond lengths AB and AC can
be calculated for all phases. These lengths can serve as good indicators of
misfit strain among various phases present in the Ti/TiN graded interface. The
phases other than Ti and TiN do not have a regular hexagonal arrangement (as
the Ti-Ti bond length depends on the distribution of N atoms around the chosen
bond). The atomic structure of the close-packed planes of Ti-0.14N, Ti-0.20N,
Ti4N3, and Ti6N5 are shown in Figure 5 (c)-(f). We identified orthonormal
lattice vectors in the close-packed plane along the directions similar to
$\left[1\bar{1}0\right]_{TiN}/\left[11\bar{2}0\right]_{Ti}$ and
$\left[11\bar{2}\right]_{TiN}/\left[10\bar{1}0\right]_{Ti}$. These lattice
vectors were used to calculate average AB and AC distances.

#####

The measured AB and AC lengths for various phases are plotted in Figure 5b.
There is a gradual increase in length of both AB and AC as N concentration
increases. This increase points to the possibility of misfit dislocations
being distributed over several planes, between various phases. Thus
dislocation density in any given plane would be significantly lower than the
sharp interface plane of Ti/TiN, which could result in unusual interaction of
dislocations and cracks with the interface.

####  3.2.2 Mechanical Properties of Ti-N phases at Ti/TiN interface

#####

Heterostructures of Ti/TiN, in the form of a single-layer TiN coating on Ti
alloys or Ti/TiN nano-multilayer coating, have been extensively explored for
various structural applications, like wear 38 and erosion resistance 39.
Mechanical properties of the heterostructure would depend on the mechanical
properties of the individual phases in the heterostructure and its variation
across the interface. We calculate mechanical properties such as bulk modulus
(B), shear modulus (G), Young’s modulus (E), Poisson’s ratio ($\nu$), Pugh’s
(G/B) ratio, and Vicker’s hardness (HV). The calculated elastic constants Cij
are reported in Table 3. The second-order elastic constants, Cij, were
calculated using the energy-strain method implemented in vaspkit 40 code by
straining the lattice from -1.5% to +1.5% with a 0.5% step size for each
independent elastic constant. Voigt-Reuss-Hill (VRH) approximation 41, 42, 43
was employed to obtain bulk modulus, shear modulus, Young’s modulus, and
Poisson’s ratio. The theoretical Vickers hardness Hv was estimated by using
Chen’s model 44, according to equation 9.

$H_{V}=2\left(\kappa^{2}G\right)^{0.585}-3$ (9)

Where $\kappa$ is the Pugh’s ratio, $\kappa$= G/B.

Figure 6: Variation of Poisson’s ratio, G/B ratio, Vickers hardness Hv, bulk
modulus B, shear modulus G, and Young’s modulus E, of Ti-N phases expected to
be present in chemically graded interface, as a function of N concentration.

#####

Poisson’s ratio decreased from Ti to TiN, indicating the material becoming
more covalent as the concentration of N increases in Ti. The value of
$\sim$0.2 and $\sim$0.4 indicates strong directional covalent and metallic
bonding, respectively 23, 45. A value of Pugh’s (G/B) ratio of $>$0.57
indicates brittleness suggesting Ti4N3, Ti6N5, and TiN are brittle and solid
solutions of N in Ti, including Ti, are ductile. In line with Poisson’s and
Pugh’s ratios, we observed a gradual increase in Vickers hardness from Ti to
TiN.

#####

Figure 6 (d)-(f) shows a gradual variation of bulk modulus, shear modulus, and
Young’s modulus as a function of N concentration of phases in the graded
interface. All of them increase smoothly as the concentration of N increases.
It has been shown that smooth variation of modulus at micrometer length scale
across the interface improves heterostructure properties like wear resistance
9 and crack suppression 46. Hence a similar gradation at atomic scale could
improve properties of heterostructures containing the graded interface. In
fact, we suspect that some of the properties of Ti/TiN heterostructures could
be due to the 3-D interface.

Table 3: The calculated elastic constants $C_{ij}$ (GPa) of the Ti–N phases. Phase | C11 | C12 | C13 | C15 | C22 | C23 | C25 | C33 | C35 | C44 | C46 | C55 | C66
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Ti | 182 | 81 | 76 |  |  |  |  | 191 |  | 45 |  |  |
Ti-0.14N | 204 | 101 | 84 |  |  |  |  | 230 |  | 70 |  |  |
Ti-0.20N | 230 | 105 | 86 |  |  |  |  | 262 |  | 75 |  |  |
Ti4N3 | 370 | 131 | 149 | -17 | 400 | 112 | 1 | 410 | -23 | 132 | 7 | 116 | 137
Ti6N5 | 423 | 131 | 135 | -16 | 425 | 134 | 15 | 428 | 2 | 156 | 15 | 169 | 150
TiN | 592 | 123 |  |  |  |  |  |  |  | 164 |  |  |

##  4 Conclusions

Using first-principles DFT calculations, we show that the atomically
chemically graded interface (which we refer to as 3-D interface) between
Ti/TiN is thermodynamically stable over a sharp interface. To establish the
existence of a 3-D interface, we define the term driving-force, which is the
sum of N vacancy formation energy in TiN and N interstitial formation energy
in Ti. The concept of driving-force can be used to assess other metal/ceramic
interfaces for the possibility of forming an atomically chemically graded
interface. We show that diffusion of N from TiN to Ti side alone is required
to form a 3-D interface. We also show that a range of stoichiometry of the
Ti-N system is stable; this explains the experimental effort to make a
chemically graded interface by controlling N partial pressure during the
growth of Ti/TiN to form a graded interface.

We show that there is smooth variation in lattice parameter and mechanical
properties like (bulk modulus, shear modulus, Young’s modulus, and hardness)
across the interface of Ti/TiN. Such structure and properties can result in
novel mechanical properties of Ti/TiN heterostructures. Gradation of structure
and properties in heterostructures at micrometer length scale has been shown
to improve their properties; a similar improvement could be expected in
heterostructures having such graded interface. Our findings suggest that
several promising properties of Ti/TiN heterostructure reported in literature
could be due to a chemically graded interface. Any attempt to understand
Ti/TiN heterostructure should involve an atomically chemically graded
interface, otherwise assumed to be atomically sharp.

### AUTHOR INFORMATION

### Corresponding Author

*Author to whom correspondence should be addressed. Electronic mail<EMAIL_ADDRESS>
### Author Contributions

The manuscript was written through contributions of all authors. All authors
have given approval to the final version of the manuscript.

### Funding

This work was supported by Science and Engineering Research Board (SERB), New
Delhi (Grant number: SRG/2019/000455) and Ministry of Human Resource and
Development (Grant number: SB20210844MMMHRD008277)

### Notes

The authors declare no competing financial interest.

We acknowledge the use of the computing resources at High Performance
Computing Environment (HPCE), IIT Madras.

## References

  * Beyerlein et al. 2015 Beyerlein, I.; Demkowicz, M.; Misra, A.; Uberuaga, B. Defect-interface interactions. _Progress in Materials Science_ 2015, _74_ , 125–210
  * Rühle and Evans 1989 Rühle, M.; Evans, A. Structure and Chemistry of Metal/Ceramic Interfaces. _Materials Science and Engineering: A_ 1989, _107_ , 187–197, Poceedings of the Symposium on Interfacial Phenomena in Composites: Processing Characterization and Mechanical Properties
  * Lackner et al. 2013 Lackner, J. M.; Waldhauser, W.; Major, B.; Major, L.; Kot, M. Plastic deformation in nano-scale multilayer materials — A biomimetic approach based on nacre. _Thin Solid Films_ 2013, _534_ , 417–425
  * Finnis 1996 Finnis, M. W. The theory of metal - ceramic interfaces. _Journal of Physics: Condensed Matter_ 1996, _8_ , 5811–5836
  * Gutiérrez Moreno and Nolan 2017 Gutiérrez Moreno, J. J.; Nolan, M. Ab Initio Study of the Atomic Level Structure of the Rutile TiO2(110)–Titanium Nitride (TiN) Interface. _ACS Applied Materials & Interfaces_ 2017, _9_ , 38089–38100, PMID: 28937740
  * Xiao and Watanabe 2015 Xiao, B.; Watanabe, S. Interface Structure in Cu/Ta2O5/Pt Resistance Switch: A First-Principles Study. _ACS Applied Materials & Interfaces_ 2015, _7_ , 519–525, PMID: 25532539
  * Bonu et al. 2019 Bonu, V.; Jeevitha, M.; Praveen Kumar, V.; Bysakh, S.; Barshilia, H. C. Ultra-thin multilayered erosion resistant Ti/TiN coatings with stress absorbing layers. _Applied Surface Science_ 2019, _478_ , 872–881
  * Swaminathan et al. 2010 Swaminathan, V. P.; Wei, R.; Gandy, D. W. Nanotechnology Coatings for Erosion Protection of Turbine Components. _Journal of Engineering for Gas Turbines and Power_ 2010, _132_ , 082104
  * Suresh 2001 Suresh, S. Graded Materials for Resistance to Contact Deformation and Damage. _Science_ 2001, _292_ , 2447–2451
  * Zhang et al. 2021 Zhang, H.; Li, Z.; He, W.; Ma, C.; Chen, J.; Liao, B.; Li, Y. Damage mechanisms evolution of TiN/Ti multilayer films with different modulation periods in cyclic impact conditions. _Applied Surface Science_ 2021, _540_ , 148366
  * Raveh et al. 1999 Raveh, A.; Weiss, M.; Pinkas, M.; Rosen, D.; Kimmel, G. Graded Al-AlN, TiN, and TiAlN multilayers deposited by radio-frequency reactive magnetron sputtering. _Surface and Coatings Technology_ 1999, _114_ , 269–277
  * Li et al. 2017 Li, N.; Yadav, S. K.; Xu, Y.; Aguiar, J. A.; Baldwin, J. K.; Wang, Y. Q.; Luo, H. M.; Misra, A.; Uberuaga, B. P. Cr incorporated phase transformation in Y2O3 under ion irradiation. _Scientific Reports_ 2017, _7_ , 1–9
  * Ljungcrantz et al. 1996 Ljungcrantz, H.; Odén, M.; Hultman, L.; Greene, J. E.; Sundgren, J. Nanoindentation studies of single‐crystal (001)‐, (011)‐, and (111)‐oriented TiN layers on MgO. _Journal of Applied Physics_ 1996, _80_ , 6725–6733
  * 14 Sigrid Wagner Ingo Arpshofen, H. J. S.; Materials Science International Team, MSIT®, The binary system Ti-O [99Wal]: Datasheet from MSI Eureka in SpringerMaterials (https://materials.springer.com/msi/phase-diagram/docs/sm_msi_r_10_010801_01_full_LnkDia0). https://materials.springer.com/msi/phase-diagram/docs/sm_msi_r_10_010801_01_full_LnkDia0
  * 15 Ivanchenko, V.; Materials Science International Team, MSIT®, C-Ti binary phase diagram: Datasheet from MSI Eureka in SpringerMaterials (https://materials.springer.com/msi/phase-diagram/docs/sm_msi_r_10_010408_01_full_LnkDia0). https://materials.springer.com/msi/phase-diagram/docs/sm_msi_r_10_010408_01_full_LnkDia0
  * 16 Zhaohui Yuan Yong Du, H. X. W. X.; Materials Science International Team, MSIT®, N-Ti phase diagram at 1 bar [1998Zen]: Datasheet from MSI Eureka in SpringerMaterials (https://materials.springer.com/msi/phase-diagram/docs/sm_msi_r_10_012892_01_full_LnkDia0). https://materials.springer.com/msi/phase-diagram/docs/sm_msi_r_10_012892_01_full_LnkDia0
  * 17 Ortrud Kubaschewski Rainer Schmid-Fetzer, u. b. L. R. L. C. O. F.; Materials Science International Team, MSIT®, Temperature - composition phase diagram of the Fe-O system: Datasheet from MSI Eureka in SpringerMaterials (https://materials.springer.com/msi/phase-diagram/docs/sm_msi_r_10_011481_02_full_LnkDia2). https://materials.springer.com/msi/phase-diagram/docs/sm_msi_r_10_011481_02_full_LnkDia2
  * Watkins et al. 2016 Watkins, E. B.; Majewski, J.; Baldwin, J. K.; Chen, Y.; Li, N.; Hoagland, R. G.; Yadav, S. K.; Liu, X. Y.; Beyerlein, I. J.; Mara, N. A. Neutron reflectometry investigations of interfacial structures of Ti/TiN layers deposited by magnetron sputtering. _Thin Solid Films_ 2016, _616_ , 399–407
  * Kupczyk 1993 Kupczyk, M. Influence of coating thickness on tool life and wear. _Surface and Coatings Technology_ 1993, _60_ , 446–449
  * Subramanian et al. 2011 Subramanian, B.; Ananthakumar, R.; Jayachandran, M. Structural and tribological properties of DC reactive magnetron sputtered titanium/titanium nitride (Ti/TiN) multilayered coatings. _Surface and Coatings Technology_ 2011, _205_ , 3485–3492
  * Jiinjen December 4, 1996 Jiinjen, A. S. Multilayer coating of a nitride-containing compound and method for producing it. European Patent EP0745699A1. December 4, 1996
  * Wolff Lodewijk, Bastin Giel and Hans 1985 Wolff Lodewijk, Bastin Giel,; Hans, H. The Titanium-Nitrogen System. _Solid State Ionics_ 1985, _16_ , 105–112
  * Yu et al. 2015 Yu, S.; Zeng, Q.; Oganov, A. R.; Frapper, G.; Zhang, L. Phase stability, chemical bonding and mechanical properties of titanium nitrides: A first-principles study. _Physical Chemistry Chemical Physics_ 2015, _17_ , 11763–11769
  * Okamoto 1993 Okamoto, H. N-Ti (Nitrogen-Titanium). _Journal of Phase Equilibria_ 1993, _14_ , 536
  * Chowdhury et al. 2021 Chowdhury, S.; Gupta, R.; Prakash, S.; Behera, L.; Gupta, M. Role of interlayer thickness on interdiffusion in Ti/TiN multilayers. _Applied Surface Science_ 2021, _564_ , 150430
  * Hohenberg and Kohn 1964 Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. _Phys. Rev._ 1964, _136_ , B864–B871
  * Kohn and Sham 1965 Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. _Phys. Rev._ 1965, _140_ , A1133–A1138
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Physical Review B_ 1996, _54_ , 11169–11186
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. _Computational Materials Science_ 1996, _6_ , 15–50
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Physical Review Letters_ 1996, _77_ , 3865–3868
  * Blöchl 1994 Blöchl, P. E. Projector augmented-wave method. _Physical Review B_ 1994, _50_ , 17953–17979
  * Monkhorst and Pack 1976 Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. _Physical Review B_ 1976, _13_ , 5188–5192
  * Momma and Izumi 2011 Momma, K.; Izumi, F. VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data. _Journal of Applied Crystallography_ 2011, _44_ , 1272–1276
  * Xue et al. 2014 Xue, J. X.; Zhang, G. J.; Guo, L. P.; Zhang, H. B.; Wang, X. G.; Zou, J.; Peng, S. M.; Long, X. G. Improved radiation damage tolerance of titanium nitride ceramics by introduction of vacancy defects. _Journal of the European Ceramic Society_ 2014, _34_ , 633–639
  * Dupressoire et al. 2021 Dupressoire, C.; Descoins, M.; Vande Put, A.; Epifano, E.; Mangelinck, D.; Emile, P.; Monceau, D. The role of nitrogen in the oxidation behaviour of a Ti6242S alloy: a nanoscale investigation by atom probe tomography. _Acta Materialia_ 2021, _216_ , 117134
  * Bars et al. 1983 Bars, J. P.; David, D.; Etchessahar, E.; Debuigne, J. Titanium $\alpha$\- Nitrogen Solid Solution Formed by High Temperature Nitriding: Diffusion of Nitrogen, Hardness, and Crystallographic Parameters. _Metallurgical Transactions A_ 1983, _14A_ , 1537–1543
  * Yang et al. 2017 Yang, W.; Ayoub, G.; Salehinia, I.; Mansoor, B.; Zbib, H. Deformation mechanisms in Ti/TiN multilayer under compressive loading. _Acta Materialia_ 2017, _122_ , 99–108
  * Cheng et al. 2010 Cheng, Y. H.; Browne, T.; Heckerman, B.; Bowman, C.; Gorokhovsky, V.; Meletis, E. I. Mechanical and tribological properties of TiN/Ti multilayer coating. _Surface and Coatings Technology_ 2010, _205_ , 146–151
  * Cao et al. 2019 Cao, X.; He, W.; Liao, B.; Zhou, H.; Zhang, H.; Tan, C.; Yang, Z. Sand particle erosion resistance of the multilayer gradient TiN/Ti coatings on Ti6Al4V alloy. _Surface and Coatings Technology_ 2019, _365_ , 214–221
  * Wang et al. 2021 Wang, V.; Xu, N.; Liu, J. C.; Tang, G.; Geng, W. T. VASPKIT: A user-friendly interface facilitating high-throughput computing and analysis using VASP code. _Computer Physics Communications_ 2021, _267_ , 108033
  * Voigt, W 1928 Voigt, W, _Lehrbuch der Kristallphysik (Leipzig: Springer-Verlag)_ 1928,
  * Reuss 1929 Reuss, A. Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle . _ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik_ 1929, _9_ , 49–58
  * Hill 1952 Hill, R. The Elastic Behaviour of a Crystalline Aggregate. _Proceedings of the Physical Society. Section A_ 1952, _65_ , 349–354
  * Chen et al. 2011 Chen, X. Q.; Niu, H.; Li, D.; Li, Y. Modeling hardness of polycrystalline materials and bulk metallic glasses. _Intermetallics_ 2011, _19_ , 1275–1281
  * Pugh 1954 Pugh, S. XCII. Relations between the Elastic Moduli and the Plastic Properties of Polycrystalline Pure Metals. _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ 1954, _45_ , 823–843
  * Jitcharoen et al. 1998 Jitcharoen, J.; Padture, N. P.; Giannakopoulos, A. E.; Suresh, S. Hertzian-Crack Suppression in Ceramics with Elastic-Modulus-Graded Surface. _Journal of the American Ceramic Society_ 1998, _81_ , 2301–2308