scaletikzpicturetowidth[1]

# Impact of Dynamical Fermions on the Centre Vortex Gluon Propagator

James C. Biddle  Waseem Kamleh  Derek B. Leinweber  Centre for the Subatomic
Structure of Matter, Department of Physics, The University of Adelaide, SA
5005, Australia

###### Abstract

The impact of $SU(3)$ centre vortices on the Landau-gauge gluon propagator is
calculated in the presence of dynamical fermions and compared to the pure
Yang-Mills case. The presence of dynamical fermions is found to alter the
behaviour of the centre vortex propagator when compared to the established
pure-gauge result. The gluon spectral representation is also explored from the
centre vortex perspective, where centre vortices are shown to exhibit clear
signs of positivity violation, which is an indicator of confinement. Vortex
removal subsequently restores positivity, demonstrating the crucial role
centre vortices play in the confinement of gluons.

††preprint: ADP-22-23/T1194

##  I Introduction

The first step in any lattice calculation is to simulate the ground state QCD
fields that permeate the vacuum. These fields give rise to the fundamental QCD
features of confinement and dynamical chiral symmetry breaking resulting in
dynamical generation of mass. Naturally, there has been substantial effort
from the community to deduce what aspect of these fields gives rise to these
features. Techniques have been developed to identify topological defects such
as Abelian monopoles [1, 2], instantons [3, 4, 5, 6] and centre vortices [7,
8, 9, 10] within these ground state fields that present possible answers to
this question. Centre vortices are of particular interest as they offer access
to the most fundamental mechanism that could underpin these phenomena.

The centre vortex picture of QCD has experienced a renewed interest in recent
years due to a host of lattice results reinforcing their importance in an
understanding of the non-perturbative properties of QCD. The majority of these
studies have been performed in the context of pure-gauge QCD, where the
effects of fermion loops are omitted from the lattice Monte-Carlo procedure.
These results have shown that vortex removal results in a loss of dynamical
mass generation [11, 12, 13], loss of string tension [14, 15] and the
suppression of the infra-red landau gauge gluon propagator [16, 15]. Vortices
alone have been capable of reproducing the qualitative picture of non-
perturbative QCD through reproduction of the linear static quark potential
[14, 17, 18], an infrared dominated gluon propagator [16], and the re-
introduction of dynamical mass generation in the low-lying hadron spectrum
[13].

Despite the success of these results, there has been a persistent quantitative
discrepancy between centre-vortex results and those from the original gauge
fields on which the vortices are identified. This manifests in a variety of
ways, most notably in a lower string tension obtained from the static quark
potential [14, 17, 18], lower hadron masses in the low-lying hadron spectrum
[13], and in residual infrared strength being retained by the vortex-removed
gluon propagator [16]. The origin of these discrepancies is as-yet unknown,
but they suggest a common property present in existing vortex studies.

Recently, work has been done to examine the effect of dynamical fermions on
the structure of centre vortices and the corresponding impact on the static
quark potential [19, 20]. These results showed for the first time a
quantitative agreement between vortex-only results and those obtained from
unmodified gauge fields. These results motivate further exploration of the
relationship between centre vortices and dynamical fermions. Here we continue
this line of investigation by calculating the Landau gauge gluon propagator on
vortex-modified configurations in the presence of dynamical fermions.

We will also examine the gluon spectral density by calculating the Euclidean
correlator to determine the presence or absence of positivity violation, as
positivity violation serves as an indicator of gluon confinement [21]. It is
well understood that positivity violation in the gluon and quark propagators
is a necessary condition for light-quark confinement [22]. As such, positivity
violation arising from centre vortices serves as a strong indication that the
centre vortex mechanism underpins the confinement of physical particles.

This paper is structured as follows. Section II outlines the vortex
identification procedure on the lattice and describes the ensembles used in
this work. Section III introduces the gluon propagator and presents the
results from our centre-vortex studies. Section IV discusses the notion of
positivity violation and presents results based on an examination of the
Euclidean correlator. Section V will summarise the findings of this paper.

##  II Vortex Identification

The procedure of centre vortex identification on the lattice is now well-
established [16, 23]. To identify the ‘thin’ lattice vortices that correspond
to the centre of physical ‘thick’ vortices [24, 25], it is necessary to bring
all configurations in an ensemble into maximal centre gauge (MCG). This is
done by applying a gauge transformation $\Omega(x)$ such that the functional

$R=\frac{1}{V\,N_{\text{dim}}\,n_{c}^{2}}\sum_{x,\mu}\left|\operatorname{Tr}U^{\Omega}_{\mu}(x)\right|^{2}\,,$
(1)

is maximised. This serves to bring every gauge link $U_{\mu}(x)$ as close as
possible to a centre element of $SU(3)$, where the group centre is comprised
of the three elements proportional to the identity

$\mathbb{Z}_{3}=\exp\left(m\frac{2\pi i}{3}\right)I,~{}m\in\\{-1,0,+1\\}.$ (2)

A more detailed description of the algorithm used to perform this updating
procedure can be found in Ref. [19].

Once fixed to MCG, a new ensemble is constructed by projecting each link onto
its nearest centre element. This centre-projected configuration is known as
the vortex-only configuration, $Z_{\mu}(x)$. From this construction, we also
define the vortex-removed ensemble as
$R_{\mu}(x)=Z_{\mu}^{\dagger}(x)\,U_{\mu}(x)$. The result of this procedure is
three ensembles:

  1. 1.

Original, untouched (UT) fields, $U_{\mu}(x)$,

  2. 2.

Vortex-only (VO) fields, $Z_{\mu}(x)$,

  3. 3.

Vortex-removed (VR) fields, $R_{\mu}(x)$,

We refer to the latter two collectively as the vortex-modified ensembles. It
is these three ensembles that we study to determine the impact of centre
vortices.

For this work, we make use of three original (UT) ensembles of 200
$32^{3}\times 64$ lattice gauge fields. Two of these are $(2+1)$ flavour
dynamical ensembles from the PACS-CS collaboration [26]. The heaviest pion
mass of $701~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ and lightest pion mass of
$156~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ are chosen to provide the greatest
range of masses to best see the impact of centre vortices as the physical
point is approached. The third ensemble is pure Yang-Mills, and provides a
reference point to previous studies. This pure-gauge ensemble has been tuned
to have a similar lattice spacing as the dynamical ensembles so that finite
volume effects should be similar across all ensembles used in this work. A
summary of the ensemble parameters is provided in Table 1.

Table 1: A summary of the lattice ensembles used in this work [26]. Type | $a\,($\mathrm{f}\mathrm{m}$)$ | $\beta$ | $\kappa_{\rm u,d}$ | $m_{\pi}\,($\mathrm{M}\mathrm{e}\mathrm{V}$)$
---|---|---|---|---
Pure gauge | 0.100 | 2.58 | - | -
Dynamical | 0.102 | 1.9 | 0.13700 | 701
Dynamical | 0.093 | 1.9 | 0.13781 | 156

##  III Gluon Propagator

###  III.1 Definition

In the continuum, the momentum-space Landau gauge gluon propagator is of the
form

$D^{ab}_{\mu\nu}(q)=\left(\delta_{\mu\nu}-\frac{q_{\mu}q_{\nu}}{q^{2}}\right)\,\delta^{ab}\,D(q^{2})\,,$
(3)

where where $D(q^{2})$ is the scalar gluon propagator. On the lattice, the
scalar propagator for $p^{2}\neq 0$ is calculated by considering [16]

$D(p^{2})=\frac{2}{3\,(n_{c}^{2}-1)\,V}\big{\langle}{\rm
Tr}\,A_{\mu}(p)\,A_{\mu}(-p)\big{\rangle}\,.$ (4)

where $n_{c}=3$ is the number of colours, $V$ is the lattice volume and
$A_{\mu}(p)$ is calculated via the discrete Fourier transform of the midpoint
definition of the gauge potential [27],

$A_{\mu}(x+\hat{\mu}/2)=\frac{1}{2i}\,\left(U_{\mu}(x)-U^{\dagger}_{\mu}(x)\right)\\\
-\frac{1}{6i}\,{\rm Tr}\left(U_{\mu}(x)-U^{\dagger}_{\mu}(x)\right)+{\cal
O}(a^{2})\,.$ (5)

As the gauge fields used in this analysis are generated using the
$\mathcal{O}(a^{2})$-improved Iwasaki action [28], the tree-level behaviour of
the gluon propagator is improved by making the substitution [29, 30, 31]

$p_{\mu}\rightarrow
q_{\mu}=\frac{2}{a}\sqrt{\sin^{2}\left(\frac{p_{\mu}a}{2}\right)+\frac{1}{3}\sin^{4}\left(\frac{p_{\mu}a}{2}\right)}\,,$
(6)

where $p_{\mu}$ are the usual lattice momentum variables

$p_{\mu}=\frac{2\pi\,n_{\mu}}{aN_{\mu}},~{}~{}n_{\mu}\in\left(-\frac{N_{\mu}}{2},\frac{N_{\mu}}{2}\right]\,,$
(7)

and $N_{\mu}$ is the lattice extent in the $\mu$ direction. The tree-level
continuum scalar propagator is then given by $D(q^{2})=\frac{1}{q^{2}}$. This
choice of momentum variables reduces the sensitivity of the gluon propagator
to finite lattice spacing effects at large momenta [32].

The perturbative scalar propagator is defined as
$D(q^{2})=\frac{Z(q^{2})}{q^{2}}$. For the remainder of this section we will
focus on the renormalisation function $Z(q^{2})=q^{2}\,D(q^{2})$. We then
renormalise $Z(q^{2})$ in the momentum space subtraction (MOM) scheme [33, 34]
on the untouched configurations by enforcing the condition that $Z^{\rm
UT}(\mu^{2})=1$ at the largest available momentum on all ensembles,
$\mu=5.5~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$. This is performed via
determination of a constant $Z_{3}^{\rm UT}$ satisfying

$\frac{Z^{\rm UT}_{\rm bare}(\mu^{2})}{Z_{3}^{\rm UT}}=Z^{\rm
UT}(\mu^{2})=1\,.$ (8)

Renormalising the vortex-modified results requires more careful consideration,
as there is no a priori method by which it should be performed. Specifically,
the problem arises from the absence of a perturbative expectation for the
vortex-only propagator. The vortex-removed results are expected to encapsulate
the high-momentum behaviour, and as such one can reasonably expect that the
MOM scheme method would apply to these ensembles. However, the vortex-only
results are dominated by infrared strength and a decay to $0$ at high
momentum. Hence, a multiplicative renormalisation based on a perturbative
expectation does not apply.

To approach this renormalisation issue, we present two sets of results. The
first set will display all propagators from an ensemble divided by $Z_{3}^{\rm
UT}$ as determined via the MOM scheme described in Eq. (8). This allows us to
readily compare the vortex-modified propagators across all ensembles.

Based on the findings of Ref. [16], we also consider renormalising the vortex-
modified propagators via a best-fit approach. To do this, we consider taking a
linear combination of the vortex-only and vortex-removed bare renormalisation
functions, $Z^{\rm VO}_{\rm bare}(q^{2})$ and $Z^{\rm VR}_{\rm bare}(q^{2})$
respectively, such that the “reconstructed” propagator

$Z^{\rm recon}(q^{2})=\frac{\zeta^{\rm VO}\,Z_{\rm bare}^{\rm
VO}(q^{2})+\zeta^{\rm VR}\,Z_{\rm bare}^{\rm VR}(q^{2})}{Z_{3}^{\rm UT}}$ (9)

is fit to $Z^{\rm UT}(q^{2})$ via a linear least-squares fit. Here,
$\zeta^{\rm VO}$ and $\zeta^{\rm VR}$ are fit parameters defined such that the
renormalised vortex-modified propagators are

$\displaystyle Z^{\rm VO}(q^{2})$ $\displaystyle=\frac{\zeta^{\rm
VO}}{Z_{3}^{\rm UT}}\,Z_{\rm bare}^{\rm VO}(q^{2})\,,$ (10) $\displaystyle
Z^{\rm VR}(q^{2})$ $\displaystyle=\frac{\zeta^{\rm VR}}{Z_{3}^{\rm
UT}}\,Z_{\rm bare}^{\rm VR}(q^{2})\,.$ (11)

Fitting the reconstructed propagator is subject to the constraint

$Z^{\rm recon}(\mu^{2})=Z^{\rm UT}(\mu^{2})=1\,,$ (12)

so that the MOM scheme is replicated in the fit. This reduces the fit to a
single parameter, as we can constrain e.g. $\zeta^{\rm VR}$ to be

$\zeta^{\rm VR}=\frac{Z_{3}^{\rm UT}-\zeta^{\rm VO}\,Z^{\rm VO}_{\rm
bare}(\mu^{2})}{Z^{\rm VR}_{\rm bare}(\mu^{2})}\,.$ (13)

Once a fit is found, the renormalisation defined in Eqs. (10, 11) is applied
such that the reconstructed propagator is simply given by the sum

$Z^{\rm recon}(q^{2})=Z^{\rm VO}(q^{2})+Z^{\rm VR}(q^{2})\,.$ (14)

As we shall see, this fitting approach is appealing as it produces excellent
agreement between the pure-gauge untouched propagator and the reconstructed
propagator as defined in Eq. (9).

###  III.2 Results

We first present the pure-gauge calculation of the scalar propagator, with all
results renormalised using the untouched renormalisation constant, $Z_{3}^{\rm
UT}$. The results from the three ensembles, UT, VO and VR are shown in Fig. 1.
As expected, these results agree with those of Ref. [16], with the untouched
propagator defined by an infrared peak and an ultraviolet plateau to tree-
level. The vortex-modified counterparts qualitatively capture these two
features, with the vortex-only propagator featuring an infrared peak, whereas
the vortex-removed results retain the ultraviolet plateau. However, there is
still significant infrared strength present in the vortex-removed propagator,
which indicates that some long-range physics remains in the vortex-removed
ensemble.

Figure 1: Pure-gauge gluon propagator as calculated on the untouched (green),
vortex-removed (orange) and vortex-only (purple) ensembles. All propagators a
renormalised by applying the renormalisation constant found by applying the
MOM scheme to the untouched propagator. A black line at $Z(q^{2})=1$ is used
to show the asymptotic behaviour.

We now consider the dynamical ensemble with the heaviest pion mass. We plot
the gluon propagator calculated on this ensemble in Fig. 2. We observe that
even at this unphysically large pion mass, the impact on the propagator is
significant. Qualitatively, the propagators retain the same features as
described for the pure-gauge sector, however the untouched propagator is
noticeably screened, as is to be expected from the introduction of dynamical
fermions [21]. The vortex-only propagator also exhibits screening, which is a
heretofore unseen effect. Furthermore, the infrared enhancement of the vortex-
removed propagator is significantly reduced when compared to the pure-gauge
results shown in Fig. 1, and now displays behaviour completely consistent with
the perturbative expectation. These two changes indicate a noticeable shift in
the behaviour of the centre vortices under the introduction of dynamical
fermions.

Figure 2: $m_{\pi}=701~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ gluon propagator.
The data scheme is as described in Fig. 1.

The story is similar for the results of the lightest pion mass, shown in Fig.
3. Screening effects are further enhanced in the untouched propagator as the
pion mass is reduced, although it is difficult to observe any change in
screening in the vortex-only propagator. To aid in this, we plot a comparison
of the vortex-only propagators across all three ensembles in Fig. 4. Here we
can clearly see the presence of screening upon introduction of dynamical
fermions. Between the two dynamical ensembles, screening effects are slightly
enhanced as the pion mass decreases, however the effect is very subtle. The
vortex-removed propagator also retains the suppression of infrared enhancement
found at $m_{\pi}=701~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$. Given that the
behaviour of the vortex-modified propagators is so similar between the two
pion masses, it appears that the mere presence of dynamical fermions plays a
substantial role in altering centre vortex structure and the manner in which
they generate the gluon propagator.

Figure 3: $m_{\pi}=156~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ gluon propagator.
The data scheme is as described in Fig. 1. Figure 4: The vortex-only
propagators from all three ensembles. Screening is distinctly visible as we
transition from pure-gauge to dynamical gauge fields.

An interesting trend in the results presented in Figs. 1, 2 and 3 is the fact
that the vortex-removed results exceed the untouched results at high momentum
with the same renormalisation constant applied. It is well understood that a
larger renormalisation constant is necessary to account for increased
roughness in an ensemble [35]. Given that the vortex removal process
represents a significant change in the texture of the gauge field, it appears
that such roughness has been induced in the vortex-removed fields. This
finding supports the need for more detailed consideration of the
renormalisation of the vortex-modified propagators.

We now repeat the above presentation but with the second renormalisation
method applied, as defined at the end of Sec. III.1. The pure-gauge results
are presented in Fig. 5. The shape of the propagators is naturally the same as
before, with the interesting addition from the renormalisation method being
the reconstructed propagator. Here we observe good agreement between the
untouched and reconstructed propagators. This indicates that the additional
degree of freedom in the renormalisation method is to some extent
encapsulating the manner in which the untouched propagator is partitioned into
its vortex-modified components.

Figure 5: Pure-gauge gluon propagator. All propagators are renormalised by
applying the renormalisation method described at the end of Sec. III.1. The
“recon” data (red) is an attempted reconstruction of the original propagator
by summing the vortex-only and vortex-removed propagators.

The dynamical ensembles with this renormalisation method applied show a
reduced agreement between the untouched and reconstructed propagators relative
to the pure gauge results. The significance of this disagreement is unknown,
and represents another interesting shift in behaviour when transitioning from
pure-gauge to dynamical QCD. The fit constants as described in Eqns. (8) and
(9) are presented in Table 2.

Figure 6: $m_{\pi}=701~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ gluon propagator. The data scheme is as described in Fig. 5. Figure 7: $m_{\pi}=156~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ gluon propagator. The data scheme is as described in Fig. 5. Table 2: The MOM scheme renormalisation constants, $Z_{3}^{\rm UT}$, as well as the fitted renormalisation constants defined in Eq. (9). Ensemble | $Z_{3}^{\rm UT}$ | $\zeta^{\rm VO}$ | $\zeta^{\rm VR}$
---|---|---|---
Pure gauge | 7.112 | 0.5543 | 0.8985
$m_{\pi}=701~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ | 9.316 | 0.6916 | 0.8251
$m_{\pi}=156~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ | 11.51 | 0.5834 | 0.8780

When comparing the vortex-only propagators with this new renormalisation
scheme we observe that screening effects remain apparent, as is evident from
Fig. 8. Furthermore, it is also possible to see a distinct increase in
screening behaviour as we transition from the heavy to light pion mass. This
suggests that perhaps this renormalisation method is more representative of
the relative contributions of the vortex-only and vortex-removed propagators
to the untouched propagator.

Figure 8: The vortex-only propagators from all three ensembles. The best fit
renormalisation method produces a greater distinction between vortex-only
propagators compared to Fig. 4.

In summary, the vortex-modified propagators undergo significant changes in
behaviour upon the introduction of dynamical fermions. Residual infrared
strength present in the pure-gauge vortex-removed propagator is suppressed in
full QCD. The vortex-only propagators effectively capture screening effects
manifesting as suppressed infrared enhancement, indicating that the long-range
behaviour of the vortex-only fields mirrors their untouched counterparts.
Best-fit renormalisation provides further insight into the structure of these
vortex fields, where we find that the sum of vortex components reconstructs
the original propagator to a fair degree. This further supports the idea that
the vortex-only and vortex-removed propagators embody a splitting of the
vacuum into long- and short-range strength respectively.

##  IV Positivity Violation

###  IV.1 Discussion

For an arbitrary two-point function $D(x-y)$ to represent correlations between
physical particles in the sense of a Wightman quantum field theory [36], it is
necessary by the Osterwalder-Schrader axioms [37] for $D(x-y)$ to satisfy

$\int
d^{4}x\,d^{4}y\,f^{*}(-x_{0},\mathbf{x})\,D(x-y)\,f(y_{0},\mathbf{y})\geq
0\,,$ (15)

for a suitable complex test function $f$. If this axiom is satisfied, then the
scalar propagator defined in Eq. (4) has spectral representation

$D(p^{2})=\int_{0}^{\infty}\,dm^{2}\,\frac{\rho(m^{2})}{p^{2}+m^{2}}\,,$ (16)

with $\rho(m^{2})\geq 0$, known as the Källen-Lehmann representation.

To investigate the behaviour of the spectral representation, we consider the
Euclidean correlator, $C(t)$, obtained by taking the Fourier transform of
$D(p_{0},\mathbf{0})$ as defined in Eq. (16) such that,

$C(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}dp_{0}\int_{0}^{\infty}dm^{2}\frac{\rho(m^{2})}{p_{0}^{2}+m^{2}}\,e^{-ip_{0}\,t}\,.$
(17)

Extending the $p_{0}$ integral to the complex plane and employing the residue
theorem, one arrives at

$C(t)=\int_{0}^{\infty}dm\,e^{-mt}\,\rho(m^{2})\,.$ (18)

Clearly if $C(t)<0$ for any $t$ then $\rho(m^{2})$ is not positive definite,
and we say that positivity has been violated. This implies that there is no
Källen-Lehmann representation as defined in Eq. (16), and as such the
propagator does not represent a correlation between physical states. Hence,
the states do not appear in the physical spectrum. In the context of the gluon
propagator, this can be taken as an indication that gluons are confined.

On the lattice [33], the Euclidean correlator, $C(t)$, is given by the
discrete Fourier transform of the temporal component of Eq. (4),

$C_{\text{lat}}(t)=\frac{1}{N_{t}}\sum_{n_{t}=0}^{N_{t}-1}e^{-2\pi
in\,t/N_{t}}D(q_{4}(n_{t})^{2})\,,$ (19)

where $N_{t}$ is the lattice extent in the temporal direction and $q_{4}$ is
the lattice momentum described in Eq. (6) and Eq. (7). As $D(0)$ is associated
with the lowest frequency mode of the propagator, it is a dominant term in Eq.
(19). As such, it is essential to ensure that finite volume effects are
accounted for.

On the lattice, finite volume effects alter the tensor structure of the
propagator given in Eq. (3) such that it has the general form [33]

$D_{\mu\nu}^{ab}(q)=\left(\delta_{\mu\nu}-\frac{h_{\mu\nu}(q)}{f\left(q^{2}\right)}\right)\,\delta^{ab}\,D\left(q^{2}\right)\,,$
(20)

where $f(q^{2})\rightarrow q^{2}$ and $h_{\mu\nu}\rightarrow q_{\mu}q_{\nu}$
for large $q_{\mu}$, but $f^{-1}(q^{2})$ is finite at $q=0$. We define

$\tilde{h}_{\mu\nu}(q)=\frac{h_{\mu\nu}(q)}{f(q^{2})}\,,$ (21)

and note that in the infinite volume limit,

$\displaystyle\tilde{h}_{\mu\mu}(q)$
$\displaystyle=f^{-1}(q^{2})\,h_{\mu\mu}(q)$
$\displaystyle=\frac{q_{\mu}q_{\mu}}{q^{2}}$
$\displaystyle=\frac{q^{2}}{q^{2}}\rightarrow 1\text{ as }q^{2}\rightarrow
0\,.$

However, on a finite volume lattice, $f^{-1}(q^{2})$ cannot approach infinity.
Since $q_{\mu}$ can take the value of $0$ and $f^{-1}(q^{2})|_{q^{2}=0}$ is
finite, $\tilde{h}_{\mu\mu}=0$ for $q_{\mu}=0$ in a finite volume. Thus, the
extraction of the scalar propagator $D(0)$ from the lattice propagator
requires a normalisation different from that of Eq. (4).

This change in normalisation can be implemented by noting that the quantity

$\left(\delta_{\mu\nu}-\frac{q_{\mu}q_{\nu}}{q^{2}}\right)\,,$ (22)

changes in the finite volume of the lattice to

$\left(\delta_{\mu\nu}-\tilde{h}_{\mu\nu}(q)\right)\,.$ (23)

For $q_{\mu}\neq 0$, setting $\mu=\nu$ and summing provides

$\sum_{\mu}\left(\delta_{\mu\mu}-\frac{q_{\mu}q_{\mu}}{q^{2}}\right)=4-1=3\,.$
(24)

But for $q_{\mu}=0$ on the lattice, $\tilde{h}_{\mu\nu}(0)=0$ and

$\sum_{\mu}\left(\delta_{\mu\mu}-\tilde{h}_{\mu\mu}(q)\right)=4\,.$ (25)

This results in

$D_{\mu\mu}^{aa}(0)=4\,(n_{c}^{2}-1)D(0)\,,$ (26)

as opposed to

$D_{\mu\mu}^{aa}(q)=3\,(n_{c}^{2}-1)D(q),~{}q\neq 0\,.$ (27)

To verify the validity of this factor we explore the behaviour of the ratio of
off-diagonal to diagonal propagator components for $q_{\mu}=0$, i.e. ratios of
the form

$\frac{D_{\mu\nu}(0)}{D_{\rho\rho}(0)}=\frac{\tilde{h}_{\mu\nu}(0)}{1-\tilde{h}_{\rho\rho}(0)},~{}\mu\neq\nu\,,$
(28)

where $\rho$ is not summed. As $\tilde{h}_{\mu\nu}\approx 0$, this ratio
provides direct access to $\tilde{h}_{\mu\nu}$ relative to the Kronecker delta
of $1$.

The values of these ratios calculated on the pure-gauge untouched
configurations are shown in Fig. 9. It is clear that these ratios are
consistent with 0 at $1\,\sigma$, indicating that both the diagonal and off-
diagonal components of $\tilde{h}_{\mu\nu}$ are small relative to $1$. These
results are corroborated by the other ensembles used in this work. This
determination justifies the use of a factor of 4 instead of 3 in calculating
the scalar propagator at zero momentum to address the impact of the finite
volume.

Figure 9: A plot of the 0-momentum ratio of the off-diagonal to diagonal
tensor gluon propagator as described in Eq. (28). We observe that the majority
of values are consistent with zero, indicating that the lattice correction
function $\tilde{h}_{\mu\nu}\rightarrow 0$ as $q\rightarrow 0$.

###  IV.2 Results

With this understanding developed, it is now possible to calculate $C(t)$ as
defined in Eq. (19). The results for the pure-gauge ensembles are shown in
Fig. 10. As expected [21], the untouched correlator shows clear signs of
positivity violation. Interestingly, the vortex-only correlators also exhibit
robust positivity violation. The positivity violation present in the vortex-
removed result at large distances is consistent with the observations made in
Fig. 1, where residual infrared strength in the vortex-removed gluon
propagator is apparent. Thus, the separation of perturbative and non-
perturbative physics through vortex modification is imperfect in the pure
gauge sector.

Figure 10: Pure-gauge Euclidean correlator. Shown are the results from the
untouched (green), vortex-removed (orange) and vortex-only (purple) ensembles.
A dashed line at $C(t)=0$ is provided to aid in observing positivity
violation.

The results from the dynamical ensembles, shown in Figs. 11 and 12 demonstrate
an interesting change in behaviour. Here we observe a similar robust violation
of positivity in the vortex-only results as observed on the pure-gauge
ensemble. However, the untouched results show a lesser degree of positivity
violation, especially on the lightest pion mass ensemble shown in Fig. 12.
Note however that violation is still present at large times.

As with the gluon propagator results in the previous section, the most
striking change is in the vortex-removed correlator. In this sector we now
observe consistency with positivity. This supports the interpretation of the
positivity violation in the vortex-removed pure-gauge results as being related
to the residual non-perturbative infrared strength in the gluon propagator. As
this residual strength is significantly diminished on the dynamical ensembles,
we now see that the residual $q^{2}$ dependence in the VR renormalisation
function may be purely perturbative in origin. In this case, vortex
modification has been successful in separating perturbative and non-
perturbative physics.

Figure 11: $m_{\pi}=701~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ Euclidean
correlator. Data is as described in Fig. 10. Figure 12:
$m_{\pi}=156~{}$\mathrm{M}\mathrm{e}\mathrm{V}$$ Euclidean correlator. Data is
as described in Fig. 10.

In summary, vortex-only configurations exhibit significant positivity
violation, as would be expected of a confining infrared-dominated theory.
Conversely, the vortex removed configurations show a loss of this positivity
violation, admitting the possibility that they do support a spectral
representation of the propagator constructed from perturbative gluon
interactions. These results provide additional support for the fact that
centre vortices encapsulate the confining aspects of QCD.

##  V Conclusion

Calculations of the behaviour of centre vortices in the presence of dynamical
fermions are new, and each calculation provides further insight into the
fascinating shift that centre vortices appear to undergo upon the introduction
of dynamical fermions. Here we have found that centre vortices in the presence
of dynamical fermions are effective in capturing the non-perturbative physics
of QCD. Moreover, vortex removal appears to also be far more effective at
removing the infrared strength of the propagator.

In regards to positivity violation, we establish the known result that
unmodified lattice ensembles give rise to positivity violation in the
Euclidean correlator [21]. We then determined that both with and without the
presence of dynamical fermions there is clear evidence that vortex-only
ensembles exhibit significant positivity violation. On our pure-gauge
ensemble, the vortex-removed correlator showed slight positivity violation at
long distances, but on both dynamical ensembles this effect vanished. In full
QCD, centre-vortex modification of the ground-state vacuum fields appears to
provide an effective separation of perturbative and non-perturbative physics.
These results present evidence that centre vortices in the QCD ground-state
vacuum fields provide the origin of confinement.

The results presented here add to the growing number of investigations into
the impact of dynamical fermions on centre vortices. We find marked
differences between the pure-gauge and dynamical case, which implores further
study. The origin of these disparities is currently unknown, and further work
into the geometric structure of centre vortices in full QCD is planned. It is
clear that there is an intimate relationship between dynamical fermions and
centre vortices, and that this relationship has significant implications for
the QCD vacuum.

###### Acknowledgements.

We thank the PACS-CS Collaboration for making their 2 +1 flavour
configurations available via the International Lattice Data Grid (ILDG). This
research was undertaken with the assistance of resources from the National
Computational Infrastructure (NCI), provided through the National
Computational Merit Allocation Scheme and supported by the Australian
Government through Grant No. LE190100021 via the University of Adelaide
Partner Share. This research is supported by Australian Research Council
through Grants No. DP190102215 and DP210103706. WK is supported by the Pawsey
Supercomputing Centre through the Pawsey Centre for Extreme Scale Readiness
(PaCER) program.

## References

  * ’t Hooft [1981] G. ’t Hooft, Topology of the Gauge Condition and New Confinement Phases in Nonabelian Gauge Theories, Nucl. Phys. B 190, 455 (1981).
  * Mandelstam [1976] S. Mandelstam, Vortices and Quark Confinement in Nonabelian Gauge Theories, Phys. Rept. 23, 245 (1976).
  * Belavin  _et al._ [1975] A. A. Belavin, A. M. Polyakov, A. S. Schwartz, and Yu. S. Tyupkin, Pseudoparticle Solutions of the Yang-Mills Equations, Phys. Lett. B 59, 85 (1975).
  * Witten [1979] E. Witten, Instantons, the Quark Model, and the 1/n Expansion, Nucl. Phys. B 149, 285 (1979).
  * Callan  _et al._ [1978] C. G. Callan, Jr., R. F. Dashen, and D. J. Gross, Toward a Theory of the Strong Interactions, Phys. Rev. D 17, 2717 (1978).
  * Aharonov  _et al._ [1978] Y. Aharonov, A. Casher, and S. Yankielowicz, Instantons and Confinement, Nucl. Phys. B 146, 256 (1978).
  * ’t Hooft [1978] G. ’t Hooft, On the Phase Transition Towards Permanent Quark Confinement, Nucl. Phys. B 138, 1 (1978).
  * ’t Hooft [1979] G. ’t Hooft, A Property of Electric and Magnetic Flux in Nonabelian Gauge Theories, Nucl. Phys. B 153, 141 (1979).
  * Cornwall [1979] J. M. Cornwall, Quark Confinement and Vortices in Massive Gauge Invariant QCD, Nucl. Phys. B 157, 392 (1979).
  * Nielsen and Olesen [1979] H. B. Nielsen and P. Olesen, A Quantum Liquid Model for the QCD Vacuum: Gauge and Rotational Invariance of Domained and Quantized Homogeneous Color Fields, Nucl. Phys. B 160, 380 (1979).
  * Trewartha  _et al._ [2015a] D. Trewartha, W. Kamleh, and D. Leinweber, Evidence that centre vortices underpin dynamical chiral symmetry breaking in SU(3) gauge theory, Phys. Lett. B 747, 373 (2015a), arXiv:1502.06753 [hep-lat] .
  * O’Malley  _et al._ [2012] E.-A. O’Malley, W. Kamleh, D. Leinweber, and P. Moran, SU(3) centre vortices underpin confinement and dynamical chiral symmetry breaking, Phys. Rev. D 86, 054503 (2012), arXiv:1112.2490 [hep-lat] .
  * Trewartha  _et al._ [2017] D. Trewartha, W. Kamleh, and D. B. Leinweber, Centre vortex removal restores chiral symmetry, J. Phys. G 44, 125002 (2017), arXiv:1708.06789 [hep-lat] .
  * Langfeld [2004] K. Langfeld, Vortex structures in pure SU(3) lattice gauge theory, Phys. Rev. D 69, 014503 (2004), arXiv:hep-lat/0307030 .
  * Bowman  _et al._ [2011] P. O. Bowman, K. Langfeld, D. B. Leinweber, A. Sternbeck, L. von Smekal, and A. G. Williams, Role of center vortices in chiral symmetry breaking in SU(3) gauge theory, Phys. Rev. D 84, 034501 (2011), arXiv:1010.4624 [hep-lat] .
  * Biddle  _et al._ [2018] J. C. Biddle, W. Kamleh, and D. B. Leinweber, Gluon propagator on a center-vortex background, Phys. Rev. D 98, 094504 (2018), arXiv:1806.04305 [hep-lat] .
  * O’Cais  _et al._ [2010] A. O’Cais, W. Kamleh, K. Langfeld, B. Lasscock, D. Leinweber, P. Moran, A. Sternbeck, and L. von Smekal, Preconditioning Maximal Center Gauge with Stout Link Smearing in SU(3), Phys. Rev. D 82, 114512 (2010), arXiv:0807.0264 [hep-lat] .
  * Trewartha  _et al._ [2015b] D. Trewartha, W. Kamleh, and D. Leinweber, Connection between center vortices and instantons through gauge-field smoothing, Phys. Rev. D 92, 074507 (2015b), arXiv:1509.05518 [hep-lat] .
  * Biddle  _et al._ [2022] J. Biddle, W. Kamleh, and D. Leinweber, Static quark potential from centre vortices in the presence of dynamical fermions, (2022), arXiv:2206.00844 [hep-lat] .
  * Virgili  _et al._ [2022] A. Virgili, W. Kamleh, and D. Leinweber, Smoothing algorithms for projected center-vortex gauge fields, (2022), arXiv:2203.09764 [hep-lat] .
  * Bowman  _et al._ [2007] P. O. Bowman, U. M. Heller, D. B. Leinweber, M. B. Parappilly, A. Sternbeck, L. von Smekal, A. G. Williams, and J.-b. Zhang, Scaling behavior and positivity violation of the gluon propagator in full QCD, Phys. Rev. D 76, 094505 (2007), arXiv:hep-lat/0703022 .
  * Roberts  _et al._ [2007] C. D. Roberts, M. S. Bhagwat, A. Holl, and S. V. Wright, Aspects of hadron physics, Eur. Phys. J. ST 140, 53 (2007), arXiv:0802.0217 [nucl-th] .
  * Biddle  _et al._ [2020] J. C. Biddle, W. Kamleh, and D. B. Leinweber, Visualization of center vortex structure, Phys. Rev. D 102, 034504 (2020), arXiv:1912.09531 [hep-lat] .
  * Engelhardt and Reinhardt [2000] M. Engelhardt and H. Reinhardt, Center projection vortices in continuum Yang-Mills theory, Nucl. Phys. B 567, 249 (2000), arXiv:hep-th/9907139 [hep-th] .
  * Bertle  _et al._ [2000] R. Bertle, M. Faber, J. Greensite, and S. Olejnik, Center vortices and color confinement in lattice QCD (2000) arXiv:hep-lat/0009017 [hep-lat] .
  * Aoki  _et al._ [2009] S. Aoki _et al._ (PACS-CS), 2+1 Flavor Lattice QCD toward the Physical Point, Phys. Rev. D 79, 034503 (2009), arXiv:0807.1661 [hep-lat] .
  * Alles  _et al._ [1997] B. Alles, D. Henty, H. Panagopoulos, C. Parrinello, C. Pittori, and D. G. Richards, $\alpha_{s}$ from the nonperturbatively renormalised lattice three gluon vertex, Nucl. Phys. B 502, 325 (1997), arXiv:hep-lat/9605033 [hep-lat] .
  * Iwasaki [1983] Y. Iwasaki, Renormalization Group Analysis of Lattice Theories and Improved Lattice Action. II. Four-dimensional non-Abelian SU(N) gauge model, (1983), arXiv:1111.7054 [hep-lat] .
  * Symanzik [1983] K. Symanzik, Continuum Limit and Improved Action in Lattice Theories. 1. Principles and phi**4 Theory, Nucl. Phys. B 226, 187 (1983).
  * Weisz and Wohlert [1984] P. Weisz and R. Wohlert, Continuum Limit Improved Lattice Action for Pure Yang-Mills Theory. 2., Nucl. Phys. B 236, 397 (1984), [Erratum: Nucl.Phys.B 247, 544 (1984)].
  * Weisz [1983] P. Weisz, Continuum Limit Improved Lattice Action for Pure Yang-Mills Theory. 1., Nucl. Phys. B 212, 1 (1983).
  * Bonnet  _et al._ [2001] F. D. R. Bonnet, P. O. Bowman, D. B. Leinweber, A. G. Williams, and J. M. Zanotti, Infinite volume and continuum limits of the Landau gauge gluon propagator, Phys. Rev. D 64, 034501 (2001), arXiv:hep-lat/0101013 .
  * Leinweber  _et al._ [1999] D. B. Leinweber, J. I. Skullerud, A. G. Williams, and C. Parrinello (UKQCD), Asymptotic scaling and infrared behavior of the gluon propagator, Phys. Rev. D 60, 094507 (1999), [Erratum: Phys.Rev.D 61, 079901 (2000)], arXiv:hep-lat/9811027 .
  * Bowman  _et al._ [2004] P. O. Bowman, U. M. Heller, D. B. Leinweber, M. B. Parappilly, and A. G. Williams, Unquenched gluon propagator in Landau gauge, Phys. Rev. D 70, 034509 (2004), arXiv:hep-lat/0402032 .
  * Bowman  _et al._ [2002] P. O. Bowman, U. M. Heller, D. B. Leinweber, and A. G. Williams, Gluon propagator on coarse lattices in Laplacian gauges, Phys. Rev. D 66, 074505 (2002), arXiv:hep-lat/0206010 .
  * Wightman [1956] A. S. Wightman, Quantum Field Theory in Terms of Vacuum Expectation Values, Phys. Rev. 101, 860 (1956).
  * Osterwalder and Schrader [1973] K. Osterwalder and R. Schrader, Axioms for Euclidean Green’s Functions, Commun. Math. Phys. 31, 83 (1973).


# Role of chalcogen vacancies and hydrogen in the optical and electrical
properties of bulk transition-metal dichalcogenides

Shoaib Khalid Princeton Plasma Physics Laboratory, P.O. Box 451, Princeton,
New Jersey 08543, USA  Anderson Janotti<EMAIL_ADDRESS>Department of
Materials Science and Engineering, University of Delaware, Newark, DE 19716,
USA  Bharat Medasani Princeton Plasma Physics Laboratory, P.O. Box 451,
Princeton, New Jersey 08543, USA

###### Abstract

Like in any other semiconductor, point defects in transition-metal
dichalcogenides (TMDs) are expected to strongly impact their electronic and
optical properties. However, identifying defects in these layered two-
dimensional materials has been quite challenging with controversial
conclusions despite the extensive literature in the past decade. Using first-
principles calculations, we revisit the role of chalcogen vacancies and
hydrogen impurity in bulk TMDs, reporting formation energies and thermodynamic
and optical transition levels. We show that the S vacancy can explain recently
observed cathodoluminescence spectra of MoS2 flakes and predict similar
optical levels in the other TMDs. In the case of the H impurity, we find it
more stable sitting on an interstitial site in the Mo plane, acting as a
shallow donor, and possibly explaining the often observed n-type conductivity
in some TMDs. We also predict the frequencies of the local vibration modes for
the H impurity, aiding its identification through Raman or infrared
spectroscopy.

There is great interest in the development of transition-metal dichalcogenide
(TMD) semiconductors for electronics [1, 2, 3, 4, 5, 6, 7, 8]. These materials
have band gaps in the near infrared range, varying with the combination of
metal(M) and chalcogen(X) atoms. Due to their layered structure, they can
potentially be combined with other 2D materials such as hexagonal boron
nitride or graphene forming van der Waals heterostructures, creating novel
materials with properties suitable for applications in flexible electronics
[9, 10, 11, 12]. As semiconductors, understanding the behavior of native
defects that can incorporated during growth or processing is an important step
toward controlling the electrical conductivity in TMDs. Defects have been
shown to affect carrier concentrations, scattering, and carrier mobilities
[13, 14, 15]. More specifically, defect engineering has become prominent in
altering the electrical [1], optical [16, 17, 18], and mechanical [19, 20]
behavior TMDs.

Identification of defects is the most important step in defect engineering.
There are various techniques used to detect the defects including microscopy
techniques such as scanning tunneling microscopy (STM), [21, 22] atomic force
microscopy (AFM), [22, 23] and optical spectroscopy techniques such as
cathodoluminescence (CL), and photoluminescence (PL) [24, 25, 26, 27]. There
is an added advantage in using optical techniques since they do not destroy
the sample. Optical techniques involve creating electron-hole pairs and
measuring the broad emission peak in the optical spectrum created by electrons
in the defect states recombining with holes. First-principles calculations
have been used previously to identify the defects by computing the optical
emission from the defects and the sub-band levels created by the defects in
the band gap for different materials [28, 29, 30]. In this work, using first
principle calculations we discuss the possible mechanism behind the near-
infrared broad emission peak related to bulk MoS2 recently measured by Fabri
et al. [24]. We also calculate the emission energy corresponding to chalcogen
vacancy for other members of the bulk TMD family. Combining the first
principle calculations with optical spectroscopy techniques provides a
powerful tool to identify individual point defects with great accuracy.

Figure 1: (a) Primitive cell of bulk TMD’s. The 180 atoms supercell used for
defect calculations with (b) H-X defect(side view) (c) Hi defect (side view)
(d) Hi defect (top view) and (e) chalcogen vacancy in the middle layer of the
defect supercell (view from just above the middle layer).

TMDs have also attracted attention as candidate materials for hydrogen storage
applications [31, 32] because of their weak interlayer interaction and large
interlayer spacing [33]. MoS2 nanotubes have already shown promising signs as
a hydrogen storage material [34, 35]. Hydrogen is also a commonly occurring
impurity during the growth of TMDs, so understanding the possible
configurations of hydrogen and its electronic behavior is important. To fully
understand the behavior of H in TMDs, we studied H in several different
configurations. We show that hydrogen interstitial, Hi, is the most stable
configuration in bulk TMDs except for MoS2. Vibrational frequencies related to
Hi in bulk TMDs are discussed and can provide direct evidence of the presence
of H in these materials through Raman and IR spectroscopy measurements.

Our density functional calculations [36, 37] are based on meta-GGA strongly
constrained and appropriately normed (SCAN) functional [38]. The interaction
between valence electrons and ionic cores is treated by projector augmented
wave (PAW) potentials [39] as implemented in VASP code [40, 41]. We used a
plane wave cutoff of 400 eV with an energy convergence of 1 meV/atom. The
Brillouin zone is sampled with 6$\times$6$\times$4 $\Gamma$-centered mesh. The
van der Waals interaction between the layers was treated by empirical
correction proposed by Grimme [42]. For calculating the Huang-Rhys (HR)
factors and total mass-weighted distortions, $\Delta$Q, we used the Nonrad
code [43]. For the defect calculations, orthogonal lattice vectors are used to
make an 180 atoms supercell, as shown in Fig. 1. Spin polarization is included
in cases with unpaired electrons. For the charge defect calculations, we used
the recently developed self-consistent charge state correction (SCPC) [44] as
implemented in the VASP code. The effect of spin-orbit coupling is also
included.

Table 1: Lattice parameter and band gap of bulk TMDs calculated using SCAN. The experimental results are taken from [45, 46, 47, 48, 49, 50]. Material | Lattice parameter (Å) | Band gap (eV)
---|---|---
a=b, c | Exp. | Eg | Exp.
MoS2 | 3.17, 12.45 | 3.15, 12.29 | 1.07 | 1.23,1.2,1.17
MoSe2 | 3.29, 13.12 | 3.29, 12.93 | 1.03 | 1.09
MoTe2 | 3.50, 14.26 | 3.51, 13.97 | 0.84 | 0.85
WS2 | 3.16, 12.51 | 3.17, 12.36 | 1.17 | 1.35
WSe2 | 3.28, 13.19 | 3.28, 12.95 | 1.13 | 1.20

For bulk TMDs including MoS2, 2$H$ is the most stable phase. It is composed of
$AB$ stacked atomic layers, where each layer forms a hexagonal lattice in the
primitive cell as shown in Fig. 1. Each layer interacts with the neighboring
layer by weak van der Waals forces. The equilibrium lattice parameters are
obtained with six atoms primitive cell using SCAN functional with Grimme Van
der Waals dispersion correction [42]. In Table 1, we list the equilibrium
lattice parameters and band gap, which are in good agreement with the
experimental data [50, 47]. In general, the in-plane and out-of-plane lattice
parameters depend strongly on the chalcogen atom and very weakly on the metal
atom. The band gap decreases as we move from the lighter S atom to the heavier
Te atom.

Figure 2: Orbital-resolved electronic band structures of (a) MoS2, (b) MoSe2,
(c) MoTe2, (d) WS2 and (e) WSe2 using the SCAN functional and including
spin–orbit coupling. The valence band maximum (VBM) is set to zero in each
case. Note that the band gap is indirect for all the studied bulk TMDs with
VBM at $\Gamma$ and CBM lies in direction between K and $\Gamma$, for MoTe2
the VBM is at K.

In Fig. 2, we show the orbital-resolved electronic band structure of bulk
TMDs. All the studied TMDs are semiconductors with indirect band gaps; the
valence-band maxima (VBM) occur at $\Gamma$ (except for MoTe2 which has VBM at
K) and the conduction-band minima (CBM) are along the high symmetry $\Gamma$-K
direction. The spin-orbit splitting increases from the lighter S atom to the
heavier Te atom for the same the metal atom. It is interesting to note that
the band dispersion decreases from S to Te, such that the VBM shifts from
$\Gamma$ to K in MoTe2, whereas the position of CBM remains unchanged. The VBM
and CBM have most of the contribution from metal $d$ states with some
contribution from the chalcogen $p$ states. The calculated band gap using SCAN
is in good agreement with the previous reports [51, 52, 53, 49] and
experimental data [45, 46, 47, 48].

The chalcogen vacancies have been observed by direct techniques such as high-
resolution transmission electron microscopy [54, 55, 56, 57, 58, 59] and also
by indirect methods based on chemical reactions and x-ray photoelectron
spectroscopy [60, 61, 62, 58, 63]. Komsa et al. [53] studied different native
defects including vacancy complexes in MoS2 and showed that the sulfur vacancy
has the lowest formation energy. In this work, we first focus on the chalcogen
vacancy. The formation energy of the chalcogen vacancy as a function of Fermi
level is shown in Fig. 3 for all the studied bulk TMDs. The chalcogen vacancy
in bulk TMDs is a deep acceptor with the (0/-1) transition level deep inside
the band gap in agreement with previous studies [53, 64, 65]. The relaxation
around the vacancy is very minimal. The formation energy of neutral chalcogen
vacancy in the chalcogen poor limit and the (0/-1) transition levels are shown
in Table 2. In the neutral charge state, the chalcogen vacancy leads to an
occupied state just below the VBM and doubly degenerate empty states within
the band gap. The degenerate states in the band gap are highly localized on
neighboring metal atoms and can be occupied by four electrons. In the -1
charge state, one of the states is half occupied. The -2 charge state is
unstable within the band gap for all the studied TMDs, in agreement with
previous studies [53, 64, 65].

Table 2: Formation energy of chalcogen vacancy in neutral charge state and (0/-1) transition level with respect to VBM for all the studied TMDs. Material | Formation energy (eV) | Transition level (eV)
---|---|---
neutral | (0/-1)
MoS2 | 1.50 | 0.68
MoSe2 | 1.77 | 0.68
MoTe2 | 2.13 | 0.58
WS2 | 1.78 | 0.73
WSe2 | 2.03 | 0.79
Figure 3: Formation energies of chalcogen vacancy(dashed blues lines for chalcogen rich and solid blue lines for metal rich), H$i$(green lines) and H-X (X=chalcogen atom)(red lines) as a function of Fermi level for (a) MoS2, (b) MoSe2, (c) MoTe2, (d) WS2 and (e) WSe2. The VBM is used as a reference and the position of CBM is also shown in each case.  Table 3: Calculated Huang-Rhys (HR) factors for initial state $S$$i$ and final state $S$$f$, total mass-weighted distortions $\Delta$Q and emission energy due to chalcogen vacancy for bulk TMD’s. Material | HR factor |  $\Delta$Q |  $E$PL |  $E$FC
---|---|---|---|---
$S$$i$ |  $S$$f$ | ($\sqrt{amu}$ Å) | (eV) | (meV)
MoS2 | 2.73 | 1.42 | 1.05 | 0.62 | 68
MoSe2 | 4.33 | 1.94 | 1.44 | 0.72 | 98
MoTe2 | 6.14 | 3.33 | 2.07 | 0.61 | 103
WS2 | 5.30 | 5.23 | 1.62 | 0.81 | 90
WSe2 | 7.54 | 8.34 | 2.17 | 0.97 | 108
Figure 4: (a) Schematic of optical emission due to chalcogen vacancy in bulk
TMDs (b) Configuration coordinate diagram for the chalcogen vacancy in TMDs.
The chalcogen vacancy in the neutral charge state is considered as the ground
state, and the energy variation as a function of atomic displacement around
the stable configuration is shown. The photoluminescence emission energy
related to chalcogen vacancy is labeled as $E$PL.

Fabri et al. [24] in their cathodoluminescence experiment observed a near-
infrared emission peak around 0.75 eV for exfoliated bulk MoS2 flakes. Gupta
et al. [66] studied the Franck-Condon shift associated with monolayer MoS2 and
related the emission energy to the S vacancy in monolayer MoS2 with the one
measured by Fabri et al.. It is important to mention here that the exfoliated
MoS2 flakes measured by Fabri et al. [24] lie in the bulk regime, and that
exciton binding energies in the bulk are much smaller than the exciton binding
energy in the monolayer MoS2. In this work, we studied the optical emission
due to chalcogen vacancy in bulk MoS2 and other members of the bulk TMD
family. The schematic of optical emission due to chalcogen vacancy in bulk
TMDs is shown in Fig. 4(a). An electron-hole pair is generated in a
photoluminescence (PL)/cathodoluminescence (CL) excitation. The generated
electron in the conduction band does not recombine with the hole in the
valence band, instead, the electron in the -1 defect state of chalcogen
vacancy recombines with a hole in the valence band emitting a photon with
energy $E$PL. In Table 3, we list the emission energy, $E$PL, for all the
studied TMDs. Our calculated emission energy for MoS2 bulk is around 0.62 eV,
in reasonable agreement with experimental observation [24]. In Fig. 4(b), we
show the schematic of the configuration coordinate diagram for the sulfur
vacancy in neutral and -1 charge states; the electron in the -1 charge state
combines with the hole in the neutral charge state emitting a photon with
energy $E$PL. The Franck-Condon shift, $E$FC, which is the energy difference
of $V$S0 in ground state configuration and $V$S0 in $V$S-1 configuration, the
total mass-weighted distortions, $\Delta$Q, and Huang-Rhys (HR) factors for
initial state $S$$i$ and final state $S$$f$ are calculated and listed in Table
3. $\Delta$Q corresponds to the difference in geometry between the initial
state and final state configuration, whereas the HR factors quantify the
electron-phonon coupling strength.

Table 4: H-X (X=Chalcogen atom) bond length for H-X configuration and H-M (M=Metal atom) bond length for H$i$ configuration for all the studied TMDs. All the values are in Å. Material | H-X bond length | H-M bond length
---|---|---
neutral | +1 | +1
MoS2 | 1.40 | 1.36 | 1.85
MoSe2 | 1.57 | 1.49 | 1.90
MoTe2 | 1.84 | 1.70 | 1.98
WS2 | 1.46 | 1.36 | 1.86
WSe2 | 1.59 | 1.49 | 1.92

Hydrogen is one of the most important impurities commonly found during the
growth of TMDs. Zhen et al. [67] showed that H2 molecule does not form a
covalent bond with MoS2 and is electrically inactive with a migration barrier
of 0.53 eV, which makes it highly unstable at high temperatures. In this work,
we studied Hi and H-X (X = chalcogen atom) as two possible configurations of H
in TMDs. We find that for Hi, the hydrogen atom sits at the interstitial site
in the plane of the metal atom inside the hexagonal lattice. In H-X
configuration, the position of H depends upon the charge state. In the neutral
charge state, the H-S bond is aligned with Mo-S bond pointing toward the empty
space between the layers, whereas the H atom lies on top of the chalcogen for
+1 charge state. Both configurations cause very minimal local relaxation. The
H-X bond length in neutral and +1 charge states are 1.40 Å and 1.36 Å,
respectively, in MoS2, which is very close to the calculated bond length of
1.40 Å in the H2S molecule. The H-X bond length for all the other TMDs is
reported in Table 4. For the Hi configuration, the neighboring metal atoms
which are in the same plane as H, are displaced outward away as a general
trend for all studied TMDs. The H-M(M=metal atom) in-plane bond length for
MoS2, MoSe2, MoTe2, WS2 and WSe2 for +1 charge state of the Hi configuration
is shown in Table 4. Singh et al. [68] discussed the interactions of hydrogen
and native point defects (hydrogenated defects) in monolayer MoS2 and found
them to act as compensating centers in the presence of hydrogen.

In Fig. 3, we show the formation energy of the two H configurations for all
the studied TMDs. Both configurations have the lowest formation energies in
chalcogen-rich limit. Hi configuration is more stable than H-X for the entire
range of Fermi level for all the TMDs except for MoS2, where H-X is slightly
more stable in agreement with previous studies for MoS2 [64]. Hi acts as a
shallow donor for all the TMDs and remains ionized in +1 charge state
throughout the range of band gap except for MoS2 and WS2. The (+/0) transition
levels are around 133 meV and 64 meV below the CBM for MoS2 and WS2,
respectively. The H-X acts as a deep donor for all the studied TMDs with (+/0)
transition level listed in Table 5. In the metal-rich limit, the Fermi level
will be pinned around 0.37 eV and 0.31 eV below the CBM for MoS2 and WS2,
respectively. At these energies, the formation energy of H-S (Hi) in +1 charge
and sulfur vacancy in -1 charge are equal for MoS2 (WS2). For all the other
TMDs, Hi acts as a shallow donor even in metal-rich limit and can lead to
n-type conductivity.

Table 5: (+/0) transition level for H-X configuration with respect to CBM for all the studied TMDs. Material | Transition level (eV)
---|---
(+/0)
MoS2 | 0.24
MoSe2 | 0.34
MoTe2 | 0.22
WS2 | 0.08
WSe2 | 0.42
Figure 5: Kohn-Sham defect states for Hi+1 configuration in bulk MoS2. The fully occupied defect state is below the VBM and localized on the H atom, whereas the two empty states are in the band gap and localized on the neighboring Mo atom. The VB and CB correspond to the valence band and conduction band respectively. Table 6: Vibrational frequencies of Hi+1 for all the TMDs. The in-plane modes are slightly different because of the different atomic environments around Hi+1. Material | Vibrational Frequencies (cm-1)
---|---
In-Plane modes | Out-of-Plane mode
MoS2 | 1550, 1526 | 603
MoSe2 | 1331, 1311 | 665
MoTe2 | 1028, 1021 | 691
WS2 | 1580, 1545 | 744
WSe2 | 1347, 1317 | 779

In Fig. 5, we show the electronic configuration of Hi in MoS2. When an H atom
is placed at the center of the hexagon in the Mo plane, it interacts with
neighboring three Mo atoms pushing them slightly outward. This interaction
leads to the creation of three defect states. In the neutral charge state, one
of the bonding-like fully occupied defect states lies around 8.28 eV below the
VBM, whereas the other two anti-bonding defect states lie in the band gap
around 0.13 eV below the CBM, one of these states is partially occupied with
one electron. In the +1 charge state, the occupied bonding state lies around
8.6 eV below the VBM, and the two empty anti-bonding states are located around
0.2 eV below the CBM. The charge densities show that these defect states are
highly localized. The bonding state has most of the contribution from the H
$s$-orbital whereas the anti-bonding empty defect states are highly localized
on neighboring Mo atoms. We also calculated the vibrational frequencies
associated with Hi in all the TMDs. There are three vibrational modes for the
Hi configuration; two of them are in-plane modes, and one mode is in the out-
of-plane direction. The calculated frequencies listed in Table. 6 will help
their identification. We also evaluated the stability of Hi by calculating the
migration barrier height of Hi+1 using nudged elastic band (NEB) method [69].
The energy barriers of Hi+1 for migration either in the plane or out-of-plane
lie in the range of 1.5 eV to 2.3 eV for all studied TMDs, making this defect
quite stable at room temperature.

In summary, the chalcogen vacancy is a deep acceptor in all the studied bulk
TMDs. The chalcogen vacancy can also lead to emission peaks near-infrared
region as measured by a recent cathodoluminescence experiment for bulk MoS2.
The calculated emission peaks for all the studied bulk TMDs provide a
guideline for probing chalcogen vacancy in future photoluminescence
(PL)/cathodoluminescence (CL) experiments. We also investigated the possible
configuration of H in bulk TMDs. Our results show that Hi is a shallow donor
in all the studied TMDs except MoS2 and WS2, and can lead to $n$-type
conductivity in these materials. The Fermi level is pinned well below the CBM
for MoS2 and WS2 in the metal-rich limit. The formation energy of Hi is lower
than H-X for most of the TMDs except MoS2. We also determined the two in-plane
and one of-plane frequency modes associated with Hi+1 for all the TMDs. Our
results provide a detailed understanding of the chalcogen vacancy and the
hydrogen impurity, providing a guideline for experiments to probe these
defects in bulk TMDs.

This work was supported by the Laboratory Directed Research and Development
(LDRD) Program (Grant No. PPPL-132) at Princeton Plasma Physics Laboratory
under U.S. Department of Energy Prime Contract No. DE-AC02-09CH11466. A.J.
acknowledges funding from the National Science Foundation (NSF) award
#OIA-2217786. The calculations were carried out at the National Energy
Research Scientific Computing Center (NERSC), a U.S. Department of Energy
Office of Science User Facility located at Lawrence Berkeley National
Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award
BES-ERCAP27253, the Stellar Cluster at Princeton University, and the DARWIN
computing system at the University of Delaware, using the NSF grant no.
1919839.

## References

  * Radisavljevic  _et al._ [2011] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Nat. Nanotechnol 6, 147 (2011).
  * Lopez-Sanchez  _et al._ [2013] O. Lopez-Sanchez, D. Lembke, M. Kayci, A. Radenovic, and A. Kis, Nat. Nanotechnol 8, 497 (2013).
  * Late  _et al._ [2013] D. J. Late, Y.-K. Huang, B. Liu, J. Acharya, S. N. Shirodkar, J. Luo, A. Yan, D. Charles, U. V. Waghmare, V. P. Dravid, _et al._ , ACS Nano 7, 4879 (2013).
  * Georgiou  _et al._ [2013] T. Georgiou, R. Jalil, B. D. Belle, L. Britnell, R. V. Gorbachev, S. V. Morozov, Y.-J. Kim, A. Gholinia, S. J. Haigh, O. Makarovsky, _et al._ , Nat. Nanotechnol 8, 100 (2013).
  * Bernardi  _et al._ [2013] M. Bernardi, M. Palummo, and J. C. Grossman, Nano Lett. 13, 3664 (2013).
  * Choi  _et al._ [2012] W. Choi, M. Y. Cho, A. Konar, J. H. Lee, G.-B. Cha, S. C. Hong, S. Kim, J. Kim, D. Jena, J. Joo, _et al._ , Adv Mater 24, 5832 (2012).
  * Koperski  _et al._ [2015] M. Koperski, K. Nogajewski, A. Arora, V. Cherkez, P. Mallet, J.-Y. Veuillen, J. Marcus, P. Kossacki, and M. Potemski, Nat. Nanotechnol. 10, 503 (2015).
  * Kern  _et al._ [2016] J. Kern, I. Niehues, P. Tonndorf, R. Schmidt, D. Wigger, R. Schneider, T. Stiehm, S. Michaelis de Vasconcellos, D. E. Reiter, T. Kuhn, _et al._ , Adv Mater 28, 7101 (2016).
  * Lee  _et al._ [2013] G.-H. Lee, Y.-J. Yu, X. Cui, N. Petrone, C.-H. Lee, M. S. Choi, D.-Y. Lee, C. Lee, W. J. Yoo, K. Watanabe, _et al._ , ACS Nano 7, 7931 (2013).
  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Li  _et al._ [2016] M.-Y. Li, C.-H. Chen, Y. Shi, and L.-J. Li, Mater Today 19, 322 (2016).
  * Velickỳ and Toth [2017] M. Velickỳ and P. S. Toth, Appl. Mater. Today 8, 68 (2017).
  * Najmaei  _et al._ [2015] S. Najmaei, J. Yuan, J. Zhang, P. Ajayan, and J. Lou, Acc. Chem. Res. 48, 31 (2015).
  * Zou and Yakobson [2015] X. Zou and B. I. Yakobson, Acc. Chem. Res. 48, 73 (2015).
  * Nan  _et al._ [2014] H. Nan, Z. Wang, W. Wang, Z. Liang, Y. Lu, Q. Chen, D. He, P. Tan, F. Miao, X. Wang, _et al._ , ACS Nano 8, 5738 (2014).
  * Mak  _et al._ [2010] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 105, 136805 (2010).
  * Xiao  _et al._ [2012] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev. Lett. 108, 196802 (2012).
  * Mak  _et al._ [2013] K. F. Mak, K. He, C. Lee, G. H. Lee, J. Hone, T. F. Heinz, and J. Shan, Nat. Mater 12, 207 (2013).
  * Bertolazzi  _et al._ [2011] S. Bertolazzi, J. Brivio, and A. Kis, ACS Nano 5, 9703 (2011).
  * Castellanos-Gomez  _et al._ [2012] A. Castellanos-Gomez, M. Poot, G. A. Steele, H. S. Van Der Zant, N. Agraït, and G. Rubio-Bollinger, Adv Mater 24, 772 (2012).
  * Liu  _et al._ [2015] H. Liu, H. Zheng, F. Yang, L. Jiao, J. Chen, W. Ho, C. Gao, J. Jia, and M. Xie, ACS Nano 9, 6619 (2015).
  * Barja  _et al._ [2019] S. Barja, S. Refaely-Abramson, B. Schuler, D. Y. Qiu, A. Pulkin, S. Wickenburg, H. Ryu, M. M. Ugeda, C. Kastl, C. Chen, _et al._ , Nat. Commun. 10, 3382 (2019).
  * Rosenberger  _et al._ [2018] M. R. Rosenberger, H.-J. Chuang, K. M. McCreary, C. H. Li, and B. T. Jonker, ACS Nano 12, 1793 (2018).
  * Fabbri  _et al._ [2016] F. Fabbri, E. Rotunno, E. Cinquanta, D. Campi, E. Bonnini, D. Kaplan, L. Lazzarini, M. Bernasconi, C. Ferrari, M. Longo, _et al._ , Nat Commun 7, 13044 (2016).
  * Yan  _et al._ [2019] J. Yan, C. Ma, Y. Huang, and G. Yang, Adv. Sci. 6, 1802092 (2019).
  * Verhagen  _et al._ [2020] T. Verhagen, V. L. Guerra, G. Haider, M. Kalbac, and J. Vejpravova, Nanoscale 12, 3019 (2020).
  * Wu  _et al._ [2017] Z. Wu, W. Zhao, J. Jiang, T. Zheng, Y. You, J. Lu, and Z. Ni, J. Phys. Chem. C 121, 12294 (2017).
  * Lyons  _et al._ [2010] J. Lyons, A. Janotti, and C. Van de Walle, Appl. Phys. Lett. 97, 152108 (2010).
  * Frodason  _et al._ [2020] Y. K. Frodason, K. Johansen, L. Vines, and J. Varley, J. Appl. Phys. 127, 075701 (2020).
  * Lyons  _et al._ [2017] J. Lyons, J. Varley, D. Steiauf, A. Janotti, and C. Van de Walle, J. Appl. Phys. 122, 035704 (2017).
  * Clauss [1972] F. J. Clauss, Solid lubricants and self-lubricating solids, Academic Press, New York, NY (1972).
  * Nikishenko  _et al._ [1982] S. Nikishenko, A. Slinkin, G. Antoshin, K. M. Minachev, and B. Nefedov, Kinet. Catal. 23, 695 (1982).
  * Makara  _et al._ [1997] V. Makara, N. Babich, N. Zakharenko, V. Pasechny, O. Rudenko, V. Surzhko, L. Kulikov, A. Semenov-Kobzar, M. Antonova, A. Chekhovsky, _et al._ , Int. J. Hydrogen Energy 22, 233 (1997).
  * Chen  _et al._ [2001] J. Chen, N. Kuriyama, H. Yuan, H. T. Takeshita, and T. Sakai, J. Am. Chem. Soc 123, 11813 (2001).
  * Chen  _et al._ [2003] J. Chen, S. Li, and Z. Tao, J. Alloys Compd. 356, 413 (2003).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev 140, A1133 (1965).
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Turiansky  _et al._ [2021] M. E. Turiansky, A. Alkauskas, M. Engel, G. Kresse, D. Wickramaratne, J.-X. Shen, C. E. Dreyer, and C. G. Van de Walle, Comput. Phys. Commun. 267, 108056 (2021).
  * Chagas da Silva  _et al._ [2021] M. Chagas da Silva, M. Lorke, B. Aradi, M. Farzalipour Tabriz, T. Frauenheim, A. Rubio, D. Rocca, and P. Deák, Phys. Rev. Lett. 126, 076401 (2021).
  * Kam and Parkinson [1982] K. Kam and B. Parkinson, J. Phys. Chem. 86, 463 (1982).
  * Lezama  _et al._ [2014] I. G. Lezama, A. Ubaldini, M. Longobardi, E. Giannini, C. Renner, A. B. Kuzmenko, and A. F. Morpurgo, 2D Mater. 1, 021002 (2014).
  * Murray and Evans [1979] R. Murray and B. Evans, J. Appl. Cryst. 12, 312 (1979).
  * Fortin and Sears [1982] E. Fortin and W. Sears, J. Phys. Chem. Solids 43, 881 (1982).
  * Tan  _et al._ [2020] A. M. Z. Tan, C. Freysoldt, and R. G. Hennig, Phys. Rev. Mater. 4, 064004 (2020).
  * James and Lavik [1963] P. B. James and M. Lavik, Acta Cryst. 16, 1183 (1963).
  * Ivanovskaya  _et al._ [2008] V. V. Ivanovskaya, A. Zobelli, A. Gloter, N. Brun, V. Serin, and C. Colliex, Phys. Rev. B 78, 134104 (2008).
  * Padilha  _et al._ [2014] J. E. Padilha, H. Peelaers, A. Janotti, and C. G. Van de Walle, Phys. Rev. B 90, 205420 (2014).
  * Komsa and Krasheninnikov [2015] H.-P. Komsa and A. V. Krasheninnikov, Phys. Rev. B 91, 125304 (2015).
  * Zhou  _et al._ [2013] W. Zhou, X. Zou, S. Najmaei, Z. Liu, Y. Shi, J. Kong, J. Lou, P. M. Ajayan, B. I. Yakobson, and J.-C. Idrobo, Nano Lett. 13, 2615 (2013).
  * Lehtinen  _et al._ [2015] O. Lehtinen, H.-P. Komsa, A. Pulkin, M. B. Whitwick, M.-W. Chen, T. Lehnert, M. J. Mohn, O. V. Yazyev, A. Kis, U. Kaiser, _et al._ , ACS Nano 9, 3274 (2015).
  * Zhu  _et al._ [2017] H. Zhu, Q. Wang, L. Cheng, R. Addou, J. Kim, M. J. Kim, and R. M. Wallace, ACS Nano 11, 11005 (2017).
  * Komsa  _et al._ [2012] H.-P. Komsa, J. Kotakoski, S. Kurasch, O. Lehtinen, U. Kaiser, and A. V. Krasheninnikov, Phys. Rev. Lett. 109, 035503 (2012).
  * Jeong  _et al._ [2017] H. Y. Jeong, Y. Jin, S. J. Yun, J. Zhao, J. Baik, D. H. Keum, H. S. Lee, and Y. H. Lee, Adv Mater 29, 1605043 (2017).
  * Leiter  _et al._ [2020] R. Leiter, Y. Li, and U. Kaiser, Nanotechnology 31, 495704 (2020).
  * Qiu  _et al._ [2013] H. Qiu, T. Xu, Z. Wang, W. Ren, H. Nan, Z. Ni, Q. Chen, S. Yuan, F. Miao, F. Song, _et al._ , Nat. Commun 4, 2642 (2013).
  * Han  _et al._ [2016] H.-V. Han, A.-Y. Lu, L.-S. Lu, J.-K. Huang, H. Li, C.-L. Hsu, Y.-C. Lin, M.-H. Chiu, K. Suenaga, C.-W. Chu, _et al._ , ACS Nano 10, 1454 (2016).
  * Chen  _et al._ [2015] B. Chen, H. Sahin, A. Suslu, L. Ding, M. I. Bertoni, F. Peeters, and S. Tongay, Acs Nano 9, 5326 (2015).
  * Tosun  _et al._ [2016] M. Tosun, L. Chan, M. Amani, T. Roy, G. H. Ahn, P. Taheri, C. Carraro, J. W. Ager, R. Maboudian, and A. Javey, ACS Nano 10, 6853 (2016).
  * Singh  _et al._ [2022] A. Singh, M. Dey, and A. K. Singh, Phys. Rev. B 105, 165430 (2022).
  * Kim  _et al._ [2022] J. Y. Kim, Ł. Gelczuk, M. P. Polak, D. Hlushchenko, D. Morgan, R. Kudrawiec, and I. Szlufarska, NPJ 2D Mater. Appl. 6, 75 (2022).
  * Gupta  _et al._ [2018] S. Gupta, S. N. Shirodkar, D. Kaplan, V. Swaminathan, and B. I. Yakobson, J. Phys. Condens. Matter 30, 095501 (2018).
  * Zhu  _et al._ [2016] Z. Zhu, H. Peelaers, and C. G. Van de Walle, Phys. Rev. B 94, 085426 (2016).
  * Singh and Singh [2019] A. Singh and A. K. Singh, Phys. Rev. B 99, 121201 (2019).
  * Henkelman  _et al._ [2000] G. Henkelman, B. P. Uberuaga, and H. Jónsson, J. Chem. Phys. 113, 9901 (2000).