First-principles study of the dynamic Jahn-Teller distortion of the neutral vacancy in diamond

0.77 MB·

by  Joseph C. A. Prentice

First principles study of the dynamic Jahn-Teller distortion of the neutral vacancy in
diamond
Joseph C. A. Prentice,1 Bartomeu Monserrat,1,2 and R. J. Needs1
1TCM Group, Cavendish Laboratory, University of Cambridge,
J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom
2Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854-8019, USA
(Dated: January 6, 2017)
First-principles density functional theory methods are used to investigate the structure, ener-
getics, and vibrational motions of the neutral vacancy defect in diamond. The measured optical
absorption spectrum demonstrates that the tetrahedral T point group symmetry of pristine dia-
7 d
1 mond is maintained when a vacancy defect is present. This is shown to arise from the presence of
0 adynamicJahn-Tellerdistortionthatisstabilisedbylargevibrationalanharmonicity. Ourcalcula-
2 tionsfurtherdemonstratethatthedynamicJahn-Teller-distortedstructureofTd symmetryislower
n in energy than the static Jahn-Teller distorted tetragonal D2d vacancy defect, in agreement with
experimentalobservations. Thetetrahedralvacancystructurebecomesmorestablewithrespectto
a
the tetragonal structure by increasing temperature. The large anharmonicity arises mainly from
J
quartic vibrations, and is associated with a saddle point of the Born-Oppenheimer surface and a
4
minimuminthefreeenergy. ThisstudydemonstratesthatthebehaviourofJahn-Tellerdistortions

ofpointdefectscanbecalculatedaccuratelyusinganharmonicvibrationalmethods. Ourworkwill
]
l open the way for first-principles treatments of dynamic Jahn-Teller systems in condensed matter.
l
a
h
- I. INTRODUCTION imentally observed in the silicon vacancy. When ap-
s
plied to diamond, this approach also predicts a D dis-
e 2d
m Point defects in crystals introduce electron energy lev- tortion, but in this case the experimental observations
. els within the electronic band gap, trap charge carriers, showthattheneutraldiamondvacancyhasTd symmetry
at emit and absorb light, and phonon scattering from them instead.23,24 (Vacancy structures with both symmetries
m limits thermal and electrical conductivities.1,2 Such ef- are shown in Figs. 1(b) and 1(c).) The experimental ob-
fects influence many of the most desirable properties of servationscanberationalizedbytheappearanceofady-
-
d diamond, including its optical properties and high ther- namic Jahn-Teller distortion,22,23,25,26 due to strong an-
n mal conductivity, and point defects in diamond have harmonic vibrational motion. The dynamic Jahn-Teller
o been the subject of much previous theoretical work.3,4 effectisobservedinavarietyofsystems,includingdoped
c Some defects such as the Si-V5,6 and N-V centres7 manganites,27,28 fullerides,29,30 octahedral complexes of
[ −
have been identified as potential “qubits” in quantum d9 ions,31 and the excited states of the N-V− centre in

1 computers.2,8–13 Both these defects include a lattice va- diamond.32
v
cancy,whichisanimportantdefectinitsownright. The InadynamicJahn-Tellersystem,therearetwoormore
8
neutralvacancyindiamondisparticularlysignificantbe- minima in the Born-Oppenheimer (BO) energy surface,
1
1 causeitisknowntobestableoverawiderangeofdoping which are separated by energy barriers. In the vacancy
1 levels,14 and plays a central role in defect diffusion.15 It in diamond, there are three minima, corresponding to
0 is also important for its optical absorption and lumines- tetragonal distortions in each Cartesian direction. If the
. cence properties – it is associated with a series of lines energy barriers between the minima are low enough, the
1
0 in the absorption spectrum of diamond, including the wavefunctionofthesystemissharedbetweentheminima
7 strongandsharpGR1lineat1.673eV16–anditsapplica- instead of localising in a single minimum, resulting in a
1 tions in quantum information17 and precision sensing.18 dynamic Jahn-Teller effect. This means that the vibra-
v: Carbon and silicon are isoelectronic, and their pris- tionalwavefunctionpossessesthesymmetryoftheundis-
i tine lattices have tetrahedral Td point group symme- tortedstructure.26Structuresatsaddlepointsormaxima
X try, as depicted in Fig. 1(a). Watkins modelled the of the BO surface can be stabilised by anharmonic vi-
r neutral vacancy in silicon using a linear combination of brations, even at zero temperature, through zero-point
a
atomic orbitals including only the four atoms surround- motion. Previous estimates, from experimental data, of
ing the vacancy,19 showing that two electrons should oc- the energy barriers between minima and an Einstein-like
cupy a triply-degenerate state. The energy of these elec- frequencyforthetetragonaldefectmodesofthevacancy
trons can be lowered by splitting this degenerate energy in diamond show that the vibrational energy quantum is
state via a static Jahn-Teller distortion19,20 of tetrago- largerthanorcomparabletothebarrierenergy,implying
nal D symmetry.14,21,22 Theoretically, the presence of a dynamic Jahn-Teller effect close to 0 K.33
2d
a Jahn-Teller distortion is revealed by the emergence in The dynamic Jahn-Teller effect in the neutral vacancy
the undistorted structure of harmonic vibrational soft indiamondiswell-establishedexperimentally,butfroma
modes, i.e., modes with imaginary frequencies; this ap- first-principlesstandpointtheviewislessclear. Thecom-
proach successfully predicts the D distortion exper- monly used harmonic approximation for lattice dynam-
2d
2
(a) (b) (c)
FIG.1. Structureofthepristinediamondlatticeandpossibledistortionsofthevacancy. (a)showsasiteinthepristinelattice
and its four nearest neighbours. (b) shows the nearest neighbours of the vacancy with a distortion of T symmetry, and (c)
d
shows the vacancy structure with D symmetry. The lengths indicate the distances between atoms in the relaxed structures.
2d
The four atoms in (c) form two pairs.
ics cannot account for the presence of a dynamic Jahn- eachsymmetrystate. InSec.IVwegiveabriefsummary
Teller distortion, as this is an intrinsically anharmonic of our results and some concluding remarks.
effect. First-principlescalculationsoftheanharmonicvi-
brationalwavefunctionofthegroundstateoftheneutral
vacancy in diamond have not been reported previously,
II. CALCULATIONAL METHODS
meaning that a proper description of the dynamic Jahn-
Teller effect in this system is lacking. In this study, we
determineforthefirsttimeananharmonicwavefunction A. Electronic calculations
which accurately describes the dynamic Jahn-Teller dis-
tortion, providing a theoretical description of the exper- DFT calculations were performed with version 7.0.3
imentally observed T symmetry of the neutral vacancy of the castep code36 and the corresponding ‘on-the-
d
in diamond. fly’ ultrasoft carbon pseudopotential.37 We have used
Wehaveusedfirst-principlesDFTmethods34 tocalcu- the local density approximation (LDA) as parametrised
late the electronic and vibrational properties of the neu- by Perdew and Zunger,38, which has been widely used
tral vacancy. The electronic ground state is expected to in previous calculations involving diamond and similar
be well-described within DFT, although some of its elec- materials.14,21,39–41 LDA-DFTcalculationsprovidealat-
tronic excited states are not because they have many- tice constant for diamond34 of 3.529 ˚A, compared to the
body (multideterminant) character. Our study is con- experimental value of 3.567 ˚A.42 In our calculations, the
fined to the electronic ground state. The vibrational cal- LDA-DFT lattice constant is generally used, although
culations were carried out using a recently proposed vi- some calculations were also carried out using the exper-
brational self-consistent field (VSCF) method,35 which imental lattice constant to investigate the effect on the
neglects non-adiabatic effects but includes anharmonic- dynamic stability of the tetrahedral symmetry state. In-
ity. creasing the lattice constant was found to increase the
We find that formation of the dynamic Jahn-Teller T stability.
d
vacancystructureleadstoasmallreductioninthevibra- To ensure the existence of a locally stable tetragonal
tional energy compared with the static Jahn-Teller D state(i.e.,onethatdoesnotrelaxbacktothetetrahedral
2d
structure when anharmonic effects are included. Our re- symmetry when geometry optimised), a 256-atom (255
sultsshowthatanharmonicnuclearmotionleadstoady- when the vacancy is present) supercell is used through-
namic Jahn-Teller distortion of the vacancy in diamond out this work. In previous DFT work on the neutral
very close to zero temperature, and we find that this re- vacancy in silicon, it was observed that supercells of at
mains true at least up to 400 K. least 256 atoms were required to obtain a stable tetrag-
The rest of this work is organised as follows: in Sec. II onal distortion,14,21 constructed as a 2 2 2 supercell
× ×
weoutlinethecomputationalmethodsemployedandgive of a 32-atom bcc unit cell. Calculations with supercells
technicaldetailsofthecalculations. InSec.IIIwepresent with less than 256 atoms show that the same is true for
the main results of our study, detailing the effects of in- diamond. The 256-atom supercells also allow us to ef-
cluding harmonic and anharmonic contributions to the fectively isolate the periodic defects from one another so
total energy, and to the vibrational density of states in that accurate structures for the tetrahedral and tetrago-
3
nal vacancies can be obtained. rections arising from the 2-dimensional subspaces of the
The geometry was optimised such that the root mean BO energy surface that they span. The terms V and V
1 2
square of the forces on all of the atoms was below in this expansion are found by mapping the BO energy
0.001 eV/˚A. A simple 2-atom unit cell was relaxed to surface as a function of the amplitude u of each nor-
obtaintherelaxedLDA-DFTlatticeconstant,anda256- mal mode, using a series of DFT calculations and then
atomsupercellwasconstructedfromthis. Avacancywas fitting cubic splines to the results. In this work, 19 dif-
then created and the internal co-ordinates relaxed. The ferent amplitudes per mode are used in mapping the BO
tetrahedral state was found by imposing the pristine lat- surface. MappingtheBOsurfaceisbyfarthemostcom-
tice symmetry before relaxing, and the tetragonal state putationally expensive part of the entire calculation –
wasfoundbyimposingatetragonaldistortiononthefour for example, mapping the V terms scales as the number
1
nearest neighbours of the vacancy, before relaxing with of modes multiplied by the number of mapping points
no symmetry constraints. In both vacancy structures, n multiplied by the cost of a DFT calculation, which
the nearest neighbours relaxed away from the vacancy scales asymptotically as nN4 for a simulation cell with
(by 0.11 ˚A in the tetrahedral case) in order to increase N atoms. TheresultinganharmonicnuclearSchr¨odinger
their sp2-like bonding , as observed in previous work.22 equation is solved using the VSCF method to give the
A plane-wave cut-off energy of 650 eV was used for anharmonic vibrational energy and free energy. The an-
relaxing the structures and for the harmonic vibrational harmonicvibrationalwavefunction Φ(u) iswrittenasa
| i
calculations, with a 5×5×5 Monkhorst-Pack k-point Hartree product of the normal modes, i|φi(ui)i. The
grid43, as the energy differences between the structures states φi(ui) arerepresentedintermsofbasisfunctions
| i Q
were very well converged for these parameters. A larger which are the eigenstates of one-dimensional harmonic
value of the cut-off energy was used for some calcula- oscillators. A total of 40 basis functions were used for
tions, corresponding to even stricter convergence crite- each degree of freedom.
ria. TheSCFcyclethresholdfortheenergywastakento
In most of this work only the V terms are included
be 10 10 eV per atom to ensure a very accurate charge 1
− in the expansion of E(u). The number of calculations
density, and thus accurate forces.
required to map the BO surface accurately as a function
of two normal co-ordinates is approximately the square
of the number required for one co-ordinate, making this
B. Vibrational calculations computationally unfeasible if all pairs of modes were to
be included. Our tests show that the V terms are small
2
A finite displacement method was used to obtain the for most modes, as the harmonic approximation works
matrixofforceconstants,whichwasFouriertransformed well for vibrations in diamond.35 The only large V2 term
to obtain the dynamical matrix.44 This was diagonalised corresponds to the 2-dimensional BO subspace spanned
toobtaintheharmonicvibrationalfrequenciesandeigen- bytwosoftmodes(u1,u2)thatarepresentinthetetrahe-
vectors. Atomic displacements of 0.00529 ˚A were used. dralstructureofthediamondvacancy. ThisV2 termhas
Theanharmonicvibrationalcalculationswereconducted three equivalent minima, corresponding to a tetragonal
using the VSCF method described in Ref. 35 and used distortionalongeachCartesiandirection. Nomatterhow
successfully several times since.45–47 the orthogonal normal co-ordinates (u1,u2) are chosen,
In this method, the BO energy surface is described in it is not possible for the axes to pass precisely through
a basis of harmonic normal co-ordinates of amplitude u. more than one of these minima. This means that it is
Inthesecoordinates,theharmonicpotentialisseparable, impossible to capture the behaviour of the system with
and each degree of freedom contributes 1ω2u2, where ω 3 minima using only the V1 terms for these soft modes –
2
is a harmonic frequency. Starting from the harmonic the associated V2 term must therefore be included. Fur-
approximation, weimprovetherepresentationoftheBO thermore, the vibrational wavefunction for this V2 sub-
surfacebyusingaprincipalaxesapproximation,48 which space of the BO surface cannot be described correctly
takestheformofasumovermany-bodyterms,truncated as a product of two states labelled by the two corre-
at second order: sponding normal modes (u1,u2). Separating the wave-
function in this way breaks the rotational symmetry of
1
the problem. Instead, for this subspace, we separate the
E(u)=E(0)+ V (u )+ V (u ,u ). (1)
1 i 2 2 i i0 wavefunction using polar co-ordinates, r = u2+u2
Xi iX06=i and θ = arctan(u /u ), where the u uare writt1en in2
u 2 1 i
p
In this expression, u is a vector containing the normal unitsof1/ 2ω ,withω theimaginaryharmonicsoft
i,s i,s
| |
mode amplitudes u , Cartesian indices are collective la- phononfrequencies. Thispreservesthecorrectsymmetry
i
p
belsforthequantumnumbers(q,ν)whereqisaphonon of the problem and allows an accurate wavefunction to
wavevectorandν aphononbranch,andE(u)isthevalue bedetermined. Theusualwavefunction φ (u ) φ (u )
1 1 2 2
| i| i
of the BO energy surface when the atomic nuclei are in is re-expressed as R(r ) T(θ ) , where R(r ) is writ-
u u u
| i| i | i
configuration u. Anharmonicity is already included in ten as a sum of isotropic harmonic oscillator radial basis
theV terms,astheyarenotconstrainedtotheharmonic states, and T(θ ) is a sum of sinusoids. A total of 80
1 u
| i
form. The V terms provide additional anharmonic cor- angular basis functions and 20 radial basis functions are
2
4
used for these calculations. As all of the other modes
Pristine
are almost harmonic and therefore well-described by the
Tetrahedral
Vsu1btsepramcse,ctahneesnimerpglyycboentardibduetdiotnoftrhoemetnheirsgpyaortfictuhlearreVs2t nits) Tεetragonal
u
of the V1 terms of the tetrahedral defect. b. F
It would be extremely computationally expensive to ar
(
map all 762 anharmonic modes for both the tetrahe- s
e
dral and the tetragonal symmetry states of the vacancy, at
even if only the V terms were included. However, we st
can take advantage1 of the small anharmonicity in pris- of -2 -1 0 1 2 3
tine diamond.35 The fact that pristine diamond is well- sity
described by the harmonic approximation implies that a n
e
necessary condition for modes to have significant anhar- d
c
monic character is that they are strongly affected by the ni
presence of the vacancy. We expect that these modes ro
will generally have short wavelengths and high energies, ect
as over longer length scales the effect of the vacancy El
will be averaged out. The effect of the vacancy on each
mode can be calculated as the difference between the
-30 -20 -10 0
vacancy harmonic vibrational density of states (vDoS)
Electronic energy (eV)
and the pristine vDoS at the frequency of each mode,
∆g(ω) = gvac(ω) gpris(ω). We can then obtain an ac-
vib − vib
curate value for the anharmonic correction to the energy
FIG. 2. Electronic density of states of diamond. The dotted
without having to map all 762 modes by simply map-
linemarkstheFermienergy,ε . TheinsetcontainstheeDoS
F
ping the modes in descending order of ∆g(ω), using the
of the vacancy gap state, showing that the tetragonal distor-
harmonic approximation for all unmapped modes, and tion splits the state into a singlet at around 0.3 eV and a
convergingtheanharmoniccorrectionwithrespecttothe doublet at around 0.3 eV. −
number of modes mapped. Converging the correction to
within 0.1 meV requires 32 modes to be mapped for the
tetrahedralstate,but350forthetetragonalstate,imply-
as predicted by the Watkins model. The inset to Fig. 2
ing that the distortion away from tetrahedral symmetry
shows that the peak splits into two upon introduction
leads to a large increase in anharmonicity.
of the tetragonal distortion, as would be expected for a
To further reduce the computational expense of the
static Jahn-Teller distortion. These peaks correspond to
anharmoniccalculations,anenergycut-offof350eVwas
the singlet and doublet states predicted by the Watkins
used, with a 5 5 5 Monkhorst-Pack k-point grid43, as model.19 Harmonic calculations show that the tetrahe-
× ×
the shape of the BO surface is well converged for these
dral state is dynamically unstable, that is, it exhibits
parameters. In some calculations larger values of these
two soft modes, with frequencies ω and ω , that cor-
1,s 2,s
parameterswereused,correspondingtoevenstrictercon-
respondtotetragonaldistortions. Unlikethetetrahedral
vergence criteria. The energy differences between struc-
symmetry state, the tetragonal configuration is dynami-
tures were converged to within 0.5 meV per atom, with
cally stable, leading to the conclusion that, at this level
the self-consistency energy threshold for the DFT calcu-
of theory, a static Jahn-Teller distortion of tetragonal
lationssetto10 6 eVperatom. Theanharmoniccorrec-
− symmetry is favoured.
tion to the energy in the pristine 256-atom supercell was
Includinganharmonicityinthetreatmentofthetetra-
calculatedfromtheanharmoniccorrectionfora16-atom
hedral configuration provides a very different picture.
fcc supercell.
Figure 3 shows the BO surface mapped along one of the
soft mode directions, split into symmetric and antisym-
metric contributions to the anharmonic energy surface.
III. RESULTS
Theantisymmetricpartmostlyarisesduetothefactthat
anygivendirectioncannotpassthroughtwominimaand
A. Jahn-Teller effect and dynamical stability the origin, meaning that both minima cannot be well
mapped by a one-dimensional slice, as discussed above.
Figure 2 shows the calculated electronic density of However, the symmetric contribution clearly shows that
states (eDoS) for the pristine, tetrahedral and tetrago- quartic anharmonicity is present. Due to the impossi-
nal vacancy structures. There is littledifference between bility of describing all three minima correctly using one-
the eDoS of the pristine and vacancy states, apart from dimensional subspaces spanned by either soft mode indi-
theappearanceofapeakinthebandgapalmostexactly vidually, it is necessary to include the full 2-dimensional
at the Fermi level in the vacancy structures. This state subspacespannedbybothsoftmodes. Figure4(a)shows
arisesfromthefour‘danglingbonds’aroundthevacancy, the V term of the BO surface mapped on the plane
2
5
FIG. 3. Anharmonic Born-Oppenheimer energy surface
mapped along the direction of one of the two soft modes of
the tetrahedral structure, as well as its decomposition into
symmetric and antisymmetric parts. The symmetric part of
the anharmonic energy surface is clearly quartic, giving two
minima.
spanned by the two soft modes, and Fig. 4(b) shows
theanharmonicgroundstatenucleardensityinthissub-
space. In addition to the slice shown in Fig. 3, fur-
ther slices through the BO surface in Fig. 4(a) can be
found in the Supplemental Material. The anharmonic
groundstatevibrationaldensityofthetetrahedralstruc-
turehaspeaksineachofthethreeminimaoftheBOsur- FIG. 4. (a) Born-Oppenheimer energy surface in the plane
face, which lowers the overall energy of the system. The spanned by the two soft modes of the tetrahedral structure.
The static tetrahedral structure lies on the local maximum
fact that the wavefunction is shared between the min-
at the origin. Three equivalent minima corresponding to the
ima shows that, when anharmonic vibrational effects are
three possible tetragonal distortions are arranged symmet-
included, the Jahn-Teller effect in this system becomes
rically around the tetrahedral structure. Each contour line
dynamicratherthanstatic,withthesystemmaintaining
represents an energy increase of 0.0615 eV.
the full Td point symmetry of the pristine lattice. (b) Anharmonic vibrational ground state probability density
The dynamical stability of the tetrahedral state is asafunctionoftheamplitudesofthesoftmodesofthetetra-
somewhat sensitive to the exact form of the BO sur- hedralstructure. Thedensityhasthreepeaksthatcorrespond
face found in the DFT calculations. If the relaxed LDA to the minima of the BO surface in (a).
lattice constant (3.529 ˚A) is used when mapping the 2-
dimensional BO surface of the pair of soft modes, the
tetrahedralstateisstilldynamicallyunstableat0K,be-
comingstableat16.9K.However,usingtheexperimental
lattice constant of 3.567 ˚A42 reduces the size of the dy- (ω +ω )/2 ω ω . The values of θ depend
1,s 2,s 1,s 2,s u
namical instability significantly, decreasing the absolute on the precise c’hoice of’the axes defined by the two soft
valueofthealreadysmallgroundstateenergyassociated modes. The displacement patterns corresponding to the
with the V2 subspace by an order of magnitude. Using modesu1andu2inthisworkarepresentedintheSupple-
the experimental lattice constant the tetrahedral state is mentalMaterial,allowingtheseminimatobeunambigu-
calculated to become dynamically stable at 8.6 K. Given ouslyidentified. Attheseminima,thefournearestneigh-
the errors inherent in DFT calculations, our results are boursofthevacancyaredisplacedfromtheirtetrahedral
consistent with the tetrahedral state being dynamically equilibriumpositions;theyaredisplacedby0.074˚Aaway
stabledowntoliquidheliumtemperatures,asimpliedby from the vacancy along the distortion direction, but by
experiment,25 and even to absolute zero. half this distance towards the vacancy in the other two
The minima in the BO surface are, in the polar co- directions. The tetrahedral structure is at a maximum
ordinates defined above, at r = 1.64/ 2ω , θ = of the BO surface along the direction of the soft modes,
u s u
| |
39 ,159 ,279 , which correspond to tetragonal distor- but a minimum along all of the other modes, placing the
◦ ◦ ◦
p
tions along the x, z and y directions, respectively. ω = tetrahedralstructureatasaddlepointoftheBOsurface.
s
6
three structures – tetrahedral, tetragonal and pristine
) – are dynamically stable at low temperature at the an-
s
nit harmonic level, we turn to their thermodynamics. The
u
b. staticlattice,vibrational,andformationenergiesat20K
ar are reported in Table I for all three structures, as at
(
s this temperature all three are dynamically stable. The
e
at harmonic energy Ehar and the anharmonic correction
y of st ∆thEeavnahca=ncEyafnohrm−aEtiohanr,enpeerrgya,towmh,ichariesgciavlecnu,laatesdwaesl:l14as
sit 0 50 100 150
en N 1
al d PTeritsrtainheedral Ef =Evac− N− Epris,
n Tetragonal
o
ati where Evac and Epris are the total energies of the system
br with and without the vacancy, respectively, and N is the
Vi number of atoms in the pristine supercell. The values of
E for the two different symmetry states of the vacancy
f
100 125 150 175 are presented in the third part of Table I, at three levels
Vibrational energy (meV)
of theory – static (electronic), harmonic vibrational, and
(a) anharmonic vibrational.
Because the tetrahedral state is dynamically unstable
)
s
nit at the harmonic level (although not at the anharmonic
u level), due to the presence of the two soft modes, a har-
b. monic vibrational energy cannot strictly be defined for
ar Pristine
( this structure. Despite this, we have included an esti-
e Tetrahedral
c Tetragonal mated value for the tetrahedral harmonic energy, calcu-
n
e lated by simply cutting out the contribution of the two
r
e softmodes, toenablecomparisonsbetweenthetwosym-
diff metry states. The fact that this value involves the un-
S
o physical removal of two modes is noted in Table I.
D
With this caveat in mind, we can look at the thermo-
v
e dynamic stability of the two symmetry states at 20 K.
v
ati When only electronic and harmonic effects are included
ul intheformationenergy,thestatewithtetragonalsymme-
m
u tryisthemoststable,althoughtheinclusionofharmonic
C
0 25 50 75 100 125 150 175 vibrational effects reduces the Jahn-Teller relaxation en-
Vibrational energy (meV) ergy – the energy difference between the tetrahedral and
tetragonal structures – from 0.275 eV to 0.158 eV. Upon
(b)
inclusion of anharmonic effects, the tetrahedral state be-
comes the most stable, as observed experimentally, by
3meV.Thepredictedfinalformationenergyfortheneu-
FIG. 5. (a) shows the harmonic vibrational density of states
tral vacancy, including anharmonic effects, is therefore
for the pristine, tetrahedral vacancy and tetragonal vacancy
8.373 eV, which is close to the estimates from experi-
structuresindiamond,shownathighenergiesabove0.1eVin
ments of 9–15 eV.49 The formation energy of the un-
the main plot, and in full in the inset. Note the main differ-
ences between the vacancy states and the pristine structure relaxed vacancy at the static level is calculated to be
occur at high energies. 8.166 eV, implying a total relaxation energy of 0.989 eV
(b) shows the difference between the harmonic and anhar- at this level of theory. (This result is not included in
moniccumulativevDoS,∆G(ε)= 0εdε0ghar(ε0)−ganh(ε0),of Table I).
the vacancy and pristine structures, constructed using Gaus- ComparingthevDoSofthevacancystructurestothat
R
sian smearing of the frequencies.
of pristine diamond gives further insight into the effect
of the vacancy on the vibrational properties. Figure
5(a) shows the harmonic vDoS for all three structures
B. Thermodynamics at high energies, with the full vDoS as an inset, and Fig.
5(b) shows the difference between the harmonic and an-
ω
The tetragonal and pristine structures are dynami- harmonic cumulative vDoS, ∆G(ω) = dω g (ω )
0 0 har 0 −
cally stable at the harmonic level, and the previous sec- g (ω ),forbothsymmetrystatesofthevacancyaswell
anh 0
R
tion shows that the tetrahedral structure is also dynam- asthepristinelattice. Thecumulativedensitiesofstates
ically stable at low temperatures when anharmonicity were formed by broadening the mode frequencies with
is accounted for. Having therefore established that all Gaussians (of width 8.163 10 4 eV for the vacancy
−
×
7
Static energies Vibrational energies Vacancy formation energies
Structure E (eV/atom) E (eV/atom) ∆E (meV/atom) Estatic (eV) Ehar (eV) Eanh (eV)
static har anh f f f
Tetrahedral 0.0292 (0.1826) 0.190 7.451 (8.343) 8.373
Tetragonal 0.0281 0.1831 0.821 7.176 8.185 8.376
Pristine 0.0000 0.1791 0.071 0.000 0.000 0.000
TABLE I. DFT static lattice energies, vibrational energies and formation energies for each structure at 20 K. The second
columnshowstheelectronicstaticlatticeenergyE peratomforthepristine,tetrahedralandtetragonalstructuresrelative
static
tothepristinestructure. ThenexttwocolumnsshowtheharmonicvibrationalenergyE andanharmonicenergycorrection
har
∆E = E E per atom for the three structures. The last three columns show the formation energy at the static,
anh anh har
harmonic, and a−nharmonic levels of theory, Estatic, Ehar and Eanh, respectively. The tetrahedral structure is dynamically
f f f
unstable at the harmonic level, as marked by (italicised brackets).
states and 5.442 10 3 eV for the pristine lattice) and Figure6showstheanharmonicformationenergyofeach
−
×
cumulatively summing them. This allows us to see how symmetry state of the vacancy for a range of tempera-
the presence of anharmonicity changes the frequencies tures up to 400 K. It is clear that the tetrahedral struc-
themselves. Figure5(a)confirmsthatthepresenceofthe tureremainsthemoststableoverthistemperaturerange
vacancyonlyhasasignificanteffectonhighenergyvibra- – indeed, the difference in the formation energies of the
tions. This justifies our approach of including only the two symmetry states increases from 0.003 to 0.177 eV
highest energy vibrational modes in the anharmonic cal- at 400 K. The calculated value of the vacancy formation
culations. The atoms neighbouring the vacancy tend to energy at room temperature (300 K) is 8.172 eV, again
have larger vibrational amplitudes than the other atoms in reasonable agreement with experimental estimates of
intheveryhighestenergymodesforbothsymmetrycon- 9–15 eV.49
figurations, while for lower energy modes the amplitudes
8.40
are comparable.
8.35
Figure 5(b) shows that the effect of anharmonicity is
much more pronounced in the tetragonal configuration V)8.30
ttvheDtaronaShiaendrertaholfeaatnseditmrapihlraiesrdtsirniazele,osatrrnupdcrtaisurterienmse,ustcthhreusmccthaualrlneegsr.etshIainnnttthhheee ergy (e88..2250
n
e8.15
changes seen in the tetragonal case. This demonstrates n Tetrahedral
that distortions away from the tetrahedral symmetry of atio8.10 Tetragonal
thepristinelatticestronglyincreasetheanharmonicityof m8.05
the phonon modes, with the tetrahedral vacancy retain- For8.00
ing the weak anharmonic character of pristine diamond.
7.95
The optical modes at high energies are clearly more af-
fected by the inclusion of anharmonicity than the low 7.90
0 100 200 300 400
energy acoustic modes. In the pristine and tetrahedral Temperature (K)
structures, anharmonicity raises the frequency of some
modes whilst lowering those of others, leading to both FIG.6. Temperaturedependenceoftheformationenergiesof
positive and negative values of ∆G. In the tetragonal thetetrahedralandtetragonalsymmetryconfigurationsofthe
configuration, however, it generally raises the frequency neutral vacancy, including anharmonic effects but neglecting
of the modes by a small amount, showing that the lead- the small thermal expansion.
ing anharmonic term is quartic in character, as cubic
anharmonicity always acts to lower the energy in one
dimension.35
Given the above results at 20 K, we briefly examine IV. CONCLUSIONS
the temperature dependence of the formation energies
of the tetrahedral and tetragonal structures. Neglecting Ourresultsshowthatthetetrahedralsymmetrystruc-
the effect of thermal expansion, which is very small for tureoftheneutralvacancyindiamondisstabiliseddown
diamond over the range of temperatures considered,50 to almost zero temperature by anharmonic vibrations.
we calculate the anharmonic vibrational contribution to The vacancy undergoes a dynamic Jahn-Teller distor-
the free energy at a set of finite temperatures, using the tion which has the full T point group symmetry of the
d
excited states of the VSCF Hamiltonian to construct a pristine system, as observed experimentally. The an-
partition function , from which we can calculate a free harmonicvibrationalwavefunctionofthetetrahedralde-
Z
energy F = k T ln .35 For these calculations, 80 ba- fect has been calculated, and shown to be shared evenly
B
− Z
sis functions were used to obtain accurate excited states. among the three minima in the Born-Oppenheimer sur-
8
face, which correspond to the three tetragonal distor- defects.
tions. We have also calculated the temperature depen-
dence of the vacancy formation energy up to 400 K. Our
value for the formation energy of the neutral vacancy V. ACKNOWLEDGEMENTS
agrees well with experimental estimates of 9–15 eV.49
J.C.A.P. and R.J.N. thank the Engineering and Phys-
ical Sciences Research Council (EPSRC) of the UK
Ourresultsfortheisolatedneutralvacancyalsoimply for financial support [EP/J017639/1]. B.M. acknowl-
that our method could be used to calculate anharmonic edges Robinson College, Cambridge, and the Cambridge
properties of other point defects, including those in di- PhilosophicalSocietyforaHenslowResearchFellowship.
amond. Two examples of such defects are the Si-V or ComputationalresourceswereprovidedbytheHighPer-
N-V centres, which are especially of interest due to their formance Computing Service at the University of Cam-
possibleuseasqubitsinquantumcomputing.7–13Studies bridge and the Archer facility of the UK’s national high-
of their vibrational properties, including anharmonicity, performancecomputingservice, forwhichaccesswasob-
would lead to a fuller understanding of these important tained via the UKCP consortium [EP/K014560/1].
1 A. M. Stoneham, Theory of Defects in Solids: Electronic L.Childress, andR.Hanson,“Heraldedentanglementbe-
Structure of Defects in Insulators and Semiconductors tween solid-state qubits separated by three metres,” Na-
(Oxford University Press, 1975). ture 497, 86 (2013).
2 T. Schro¨der, S. L. Mouradian, J. Zheng, M. E. Trusheim, 12 H. S. Knowles, D. M. Kara, and M. Atatu¨re, “Observ-
M. Walsh, E. H. Chen, L. Li, I. Bayn, and D. Englund, ing bulk diamond spin coherence in high-purity nanodia-
“Quantum nanophotonics in diamond,” J. Opt. Soc. Am. monds,” Nat. Mater. 13, 21 (2014).
B 33, B65 (2016). 13 F. Dolde, I. Jakobi, B. Naydenov, N. Zhao, S. Pezzagna,
3 C. Freysoldt, B. Grabowski, H. Tilmann, J. Neugebauer, C. Trautmann, J. Meijer, P. Neumann, F. Jelezko, and
G. Kresse, A. Janotti, and C. G. Van de Walle, “First- J. Wrachtrup, “Room-temperature entanglement between
principles calculations for point defects in solids,” Rev. single defect spins in diamond,” Nat. Phys. 9, 139 (2013).
Mod. Phys. 86, 253 (2014). 14 F.CorsettiandA.A.Mostofi,“System-sizeconvergenceof
4 D. A. Drabold and S. Estreicher, eds., Theory of Defects point defect properties: The case of the silicon vacancy,”
in Semiconductors, Topics in Applied Physics, Vol. 104 Phys. Rev. B 84, 035209 (2011).
(Springer-Verlag, 2007). 15 R.Q.Hood,P.R.C.Kent,R.J.Needs, andP.R.Briddon,
5 L. J. Rogers, K. D. Jahnke, M. W. Doherty, A. Diet- “Quantum Monte Carlo study of the optical and diffusive
rich, L. P. McGuinness, C. Mu¨ller, T. Teraji, H. Sumiya, properties of the vacancy defect in diamond,” Phys. Rev.
J.Isoya,N.B.Manson, andF.Jelezko,“Electronicstruc- Lett. 91, 076403 (2003).
tureofthenegativelychargedsilicon-vacancycenterindi- 16 J. Walker, “Optical absorption and luminescence in dia-
amond,” Phys. Rev. B 89, 235101 (2014). mond,” Rep. Prog. Phys. 42, 1605 (1979).
6 A.GaliandJ.R.Maze,“Abinitiostudyofthesplitsilicon- 17 J. Wrachtrup and F. Jelezko, “Processing quantum infor-
vacancy defect in diamond: Electronic structure and re- mation in diamond,” J. Phys.: Condens. Matter 18, S807
lated properties,” Phys. Rev. B 88, 235205 (2013). (2006).
7 P. C. Maurer, G. Kucsko, C. Latta, L. Jiang, N. Y. Yao, 18 R. Schirhagl, K. Chang, M. Loretz, and C. L. Degen,
S. D. Bennett, F. Pastawski, D. Hunger, N. Chisholm, “Nitrogen-vacancy centers in diamond: nanoscale sensors
M. Markham, D. J. Twitchen, J. I. Cirac, and M. D. for physics and biology,” Ann. Rev. Phys. Chem. 65, 83
Lukin, “Room-temperature quantum bit memory exceed- (2014).
ing one second,” Science 336, 1283 (2012). 19 G. D. Watkins, “The Lattice Vacancy in Silicon,” in Deep
8 L. J. Rogers, K. D. Jahnke, M. H. Metsch, A. Sipahigil, Centers in Semiconductors, edited by Sokrates T. Pan-
J.M.Binder,T.Teraji,H.Sumiya,J.Isoya,M.D.Lukin, telides (Gordon and Breach, 1986).
P. Hemmer, and F. Jelezko, “All-optical initialization, 20 H. A. Jahn and E. Teller, “Stability of polyatomic
readout,andcoherentpreparationofsinglesilicon-vacancy moleculesindegenerateelectronicstates.I.Orbitaldegen-
spins in diamond,” Phys. Rev. Lett. 113, 263602 (2014). eracy,” Proc. R. Soc. A 161, 220 (1937).
9 A.Sipahigil,K.D.Jahnke,L.J.Rogers,T.Teraji,J.Isoya, 21 M. I. J. Probert and M. C. Payne, “Improving the con-
A.S.Zibrov,F.Jelezko, andM.D.Lukin,“Indistinguish- vergence of defect calculations in supercells: An ab initio
able photons from separated silicon-vacancy centers in di- study of the neutral silicon vacancy,” Phys. Rev. B 67,
amond,” Phys. Rev. Lett. 113, 113602 (2014). 075204 (2003).
10 G. Balasubramanian, P. Neumann, D. Twitchen, 22 S. J. Breuer and P. R. Briddon, “Ab initio investigation
M.Markham,R.Kolesov,N.Mizuochi,J.Isoya,J.Achard, of the native defects in diamond and self-diffusion,” Phys.
J. Beck, J. Tissler, V. Jacques, P. R. Hemmer, F. Jelezko, Rev. B 51, 6984 (1995).
and J. Wrachtrup, “Ultralong spin coherence time in iso- 23 G.Davies,“TheJahn-Tellereffectandvibroniccouplingat
topicallyengineereddiamond,”Nat.Mater.8,383(2009). deeplevelsindiamond,”Rep.Prog.Phys.44,787(1981).
11 H.Bernien,B.Hensen,W.Pfaff,G.Koolstra,M.S.Blok, 24 C. D. Clark and J. Walker, “The neutral vacancy in dia-
L.Robledo,T.H.Taminiau,M.Markham,D.J.Twitchen, mond,” Proc. R. Soc. A 334, 241 (1973).
9
25 M. Lannoo and A. M. Stoneham, “The optical absorption 38 J.P.PerdewandA.Zunger,“Self-interactioncorrectionto
oftheneutralvacancyindiamond,”J.Phys.Chem.Solids density-functional approximations for many-electron sys-
29, 1987–2000 (1968). tems,” Phys. Rev. B 23, 5048 (1981).
26 F. S. Ham, “The Jahn-Teller effect,” Int. J. Quantum 39 O. K. Al-Mushadani and R. J. Needs, “Free-energy calcu-
Chem. 5, 191 (1971). lations of intrinsic point defects in silicon,” Phys. Rev. B
27 A. J. Millis, B. I. Shraiman, and R. Mueller, “Dy- 68, 235205 (2003).
namicJahn-Tellereffectandcolossalmagnetoresistancein 40 K. Kunc, I. Loa, and K. Syassen, “Equation of state and
La Sr MnO ,” Phys. Rev. Lett. 77, 175 (1996). phonon frequency calculations of diamond at high pres-
1 x x 3
28 V. D−ediu, C. Ferdeghini, F. C. Matacotta, P. Nozar, and sures,” Phys. Rev. B 68, 094107 (2003).
G. Ruani, “Jahn-Teller dynamics in charge-ordered man- 41 R. Maezono, A. Ma, M. D. Towler, and R. J. Needs,
ganites from Raman spectroscopy,” Phys. Rev. Lett. 84, “EquationofstateandRamanfrequencyofdiamondfrom
4489 (2000). quantum Monte Carlo simulations,” Phys. Rev. Lett. 98,
29 V. Brouet, H. Alloul, T. Le, S. Garaj, and L. Forro, 025701 (2007).
“Role of dynamic Jahn-Teller distortions in Na C and 42 O.Madelung,ed.,Semiconductors-BasicData (Springer,
2 60
Na CsC studied by NMR,” Phys. Rev. Lett. 86, 4680 1996).
2 60
(2001). 43 H. J. Monkhorst and J. D. Pack, “Special points for
30 S.E.Canton,A.J.Yencha,E.Kukk,J.D.Bozek,M.C.A. Brillouin-zoneintegrations,”Phys.Rev.B13,5188(1976).
Lopes, G. Snell, and N. Berrah, “Experimental evidence 44 K. Kunc and R. M. Martin, “Ab initio force constants of
of a dynamic Jahn-Teller effect in C+,” Phys. Rev. Lett. GaAs: A new approach to calculation of phonons and di-
60
89, 045502 (2002). electric properties,” Phys. Rev. Lett. 48, 406 (1982).
31 M.C.M.O’Brien,“ThedynamicJahn-Tellereffectinoc- 45 S. Azadi, B. Monserrat, W. M. C. Foulkes, and R. J.
tahedrally co-ordinated d9 ions,” Proc. Roy. Soc. A 281, Needs, “Dissociation of high-pressure solid molecular hy-
323 (1964). drogen: A quantum Monte Carlo and anharmonic vibra-
32 K. C. Fu, C. Santori, P. E. Barclay, L. J. Rogers, N. B. tional study,” Phys. Rev. Lett. 112, 165501 (2014).
Manson, and R. G. Beausoleil, “Observation of the dy- 46 B. Monserrat, N. D. Drummond, Chris J. Pickard, and
namic Jahn-Teller effect in the excited states of nitrogen- R.J.Needs,“Electron-phononcouplingandthemetalliza-
vacancycentersindiamond,”Phys.Rev.Lett.103,256404 tion of solid helium at terapascal pressures,” Phys. Rev.
(2009). Lett. 112, 055504 (2014).
33 A.M.Stoneham,“Thelow-lyinglevelsoftheGR1centre 47 E.A.Engel,B.Monserrat, andR.J.Needs,“Anharmonic
in diamond,” Solid State Commun. 21, 339 (1977). nuclearmotionandtherelativestabilityofhexagonaland
34 R.O.Jones,“Densityfunctionaltheory: Itsorigins,riseto cubic ice,” Phys. Rev. X 5, 021033 (2015).
prominence,andfuture,”Rev.Mod.Phys.87,897(2015). 48 J.O.JungandR.B.Gerber,“Vibrationalwavefunctions
35 B.Monserrat,N.D.Drummond, andR.J.Needs,“Anhar- and spectroscopy of (H O) , n = 2,3,4,5: Vibrational
2 n
monic vibrational properties in periodic systems: energy, self-consistentfieldwithcorrelationcorrections,”J.Chem.
electron-phonon coupling, and stress,” Phys. Rev. B 87, Phys. 105, 10332 (1996).
144302 (2013). 49 J. C. Bourgoin, “An experimental estimation of the va-
36 S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, cancy formation energy in diamond,” Radiat. Eff. Defects
M. I. J. Probert, K. Refson, and M. C. Payne, “First Solids 79, 235 (1983).
principles methods using CASTEP,” Z. Kristallogr. 220, 50 S.StoupinandY.V.Shvyd’ko,“Thermalexpansionofdi-
567 (2005). amondatlowtemperatures,”Phys.Rev.Lett.104,085901
37 D. Vanderbilt, “Soft self-consistent pseudopotentials in a (2010).
generalized eigenvalue formalism,” Phys. Rev. B 41, 7892
(1990).
Supplemental Material for “First principles study of the dynamic Jahn-Teller
distortion of the neutral vacancy in diamond”
Joseph C. A. Prentice,1 Bartomeu Monserrat,1,2 and R. J. Needs1
1TCM Group, Cavendish Laboratory, University of Cambridge,
J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom
2Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854-8019, USA
(Dated: January 6, 2017)
7 I. PSEUDOPOTENTIALS
1
0 All density functional theory calculations were performed using CASTEP version 7.0.3, and its own “on-the-fly”
2
ultrasoft pseudopotential for carbon. The definition string for the pseudopotential was:

n
a 2|1.4|1.4|1.3|6|10|12|20:21(qc=6)
J

4
II. EQUILIBRIUM POSITIONS
]
l
l
a The equilibrium positions of the 255 atoms in the tetrahedral symmetry structure with the experimental lattice
h constant are given in the file TdSymmetry.cif. The 4 nearest neighbours of the vacancy are the 1st, 37th, 41st and
-
s 69th atoms listed.
e
m
. III. SOFT MODE DISPLACEMENT PATTERNS
t
a
m
The displacement patterns of the two soft modes, labelled by u and u in the main text, are given in the files
1 2
- u1Displacements.cif and u2Displacements.cif. The 4 nearest neighbours of the vacancy are the 1st, 37th, 41st
d
and 69th atoms listed.
n
o
c
[ IV. BORN-OPPENHEIMER ENERGY SURFACE SLICES

1
v Figure 1 shows two cross-sections of the BO surface in the plane spanned by the two soft modes of the tetrahedral
8 structure. Figure 1(a) shows a slice running from one of the three minima of the surface to the opposite high point
1 passing through the origin of coordinates (u ,u )=(0,0). Figure 1(b) shows a slice that runs from one minimum to
1 2
1 an adjacent minimum, and this line does not cross the origin of coordinates.
1
0
.
1
0
7
1
:
v
i
X
r
a