# A substitutional quantum defect in WS2 discovered by high-throughput
computational screening and fabricated by site-selective STM manipulation

John C. Thomas  Molecular Foundry, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, United States of America  Materials Sciences Division,
Lawrence Berkeley National Laboratory, Berkeley, CA, United States of America
Wei Chen  Yihuang Xiong  Bradford A. Barker  Department of Physics, University
of California, Merced, Merced, CA 95343, USA  Junze Zhou  Molecular Foundry,
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, United States of
America  Weiru Chen  Thayer School of Engineering, Dartmouth College, Hanover,
NH 03755, USA  Antonio Rossi  Molecular Foundry, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, United States of America  Materials Sciences
Division, Lawrence Berkeley National Laboratory, Berkeley, CA, United States
of America  Advanced Light Source, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, United States of America  Nolan Kelly  Department of
Physics, University of California, Merced, Merced, CA 95343, USA  Zhuohang Yu
Department of Materials Science and Engineering, The Pennsylvania State
University, University Park, PA 16082 United States of America  Center for
Two-Dimensional and Layered Materials, The Pennsylvania State University,
University Park, PA, 16802 United States of America  Da Zhou  Department of
Physics, The Pennsylvania State University, University Park, PA, 16802 United
States of America  Shalini Kumari  Department of Materials Science and
Engineering, The Pennsylvania State University, University Park, PA 16082
United States of America  Center for Two-Dimensional and Layered Materials,
The Pennsylvania State University, University Park, PA, 16802 United States of
America  Edward S. Barnard  Molecular Foundry, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, United States of America  Joshua A. Robinson
Department of Materials Science and Engineering, The Pennsylvania State
University, University Park, PA 16082 United States of America  Center for
Two-Dimensional and Layered Materials, The Pennsylvania State University,
University Park, PA, 16802 United States of America  Department of Physics,
The Pennsylvania State University, University Park, PA, 16802 United States of
America  Department of Chemistry, The Pennsylvania State University,
University Park, PA, 16802 United States of America  Mauricio Terrones
Department of Materials Science and Engineering, The Pennsylvania State
University, University Park, PA 16082 United States of America  Center for
Two-Dimensional and Layered Materials, The Pennsylvania State University,
University Park, PA, 16802 United States of America  Department of Physics,
The Pennsylvania State University, University Park, PA, 16802 United States of
America  Department of Chemistry, The Pennsylvania State University,
University Park, PA, 16802 United States of America  Adam Schwartzberg
Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
United States of America  D. Frank Ogletree  Molecular Foundry, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, United States of America
Eli Rotenberg  Advanced Light Source, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, United States of America  Marcus M. Noack  Applied
Mathematics and Computational Research Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, United States of America  Sinéad Griffin
Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
United States of America  Materials Sciences Division, Lawrence Berkeley
National Laboratory, Berkeley, CA, United States of America  Archana Raja
Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
United States of America  Materials Sciences Division, Lawrence Berkeley
National Laboratory, Berkeley, CA, United States of America  David A. Strubbe
Department of Physics, University of California, Merced, Merced, CA 95343, USA
Gian-Marco Rignanese  Institute of Condensed Matter and Nanoscicence,
Université catholique de Louvain, Louvain-la-Neuve 1348, Belgium  Alexander
Weber-Bargioni  Molecular Foundry, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, United States of America  Geoffroy Hautier

## Abstract

Point defects in two-dimensional materials are of key interest for quantum
information science. However, the space of possible defects is immense, making
the identification of high-performance quantum defects extremely challenging.
Here, we perform high-throughput (HT) first-principles computational screening
to search for promising quantum defects within WS2, which present localized
levels in the band gap that can lead to bright optical transitions in the
visible or telecom regime. Our computed database spans more than 700 charged
defects formed through substitution on the tungsten or sulfur site. We found
that sulfur substitutions enable the most promising quantum defects. We
computationally identify the neutral cobalt substitution to sulfur (Co${}_{\rm
S}^{0}$) as very promising and fabricate it with scanning tunneling microscopy
(STM). The Co${}_{\rm S}^{0}$ electronic structure measured by STM agrees with
first principles and showcases an attractive new quantum defect. Our work
shows how HT computational screening and novel defect synthesis routes can be
combined to design new quantum defects.

## Introduction

Point defects in semiconductors are considered as building blocks for quantum
information science (QIS) applications. Optically-active quantum defects
(OQDs) can be used in quantum sensing, memory, and networks1, 2, 3, 4. The
performance of an OQD depends on its fundamental properties and limitations
that can vary across defects5, 6. Certain defects, such as the silicon-
divacancy center in diamond, show robust optical coherence but low spin
coherence, while the NV- center in diamond shows high spin coherence but lower
optical coherence7, 8. The identification of new OQDs in a specific host with
optimal spin, optical, and electronic properties is essential to the
development of QIS applications.

Two-dimensional (2D) materials, particularly transition metal dichalcogenides
(TMDs), provide an enormous phase space of functionality with tunable and
exceptional spin, optical and electronic properties9, 10, 11, 12, 13, 14, 15,
16, 17. Additionally, as materials are reduced from bulk to lower
dimensionality, the spin-coherence lifetime of an OQD is expected to increase
18. A decisive factor for an OQD is the appearance of in-gap localized states
making it important to understand and measure the electronic levels induced by
a defect in a given 2D host. While a number of techniques can routinely
resolve the atomic lattice, the electronic levels introduced by the defect in
the host are not easily accessible by most experimental techniques. However,
scanning tunneling microscopy (STM) and scanning tunneling spectroscopy (STS)
can probe atomic-scale defects at the required length scale19, 20. This has
been used to characterize many defects in 2D materials, e.g., carbon radical
dopants, chalcogen vacancies, oxygen substitutions, and a variety of metal
substitutions21, 22, 16, 15, 23, 19. Next to these experimental developments,
first-principles approaches have been successfully used to compute and
understand the properties of quantum defects in bulk semiconductors and 2D
hosts24, 25, 26, 27. First principles techniques have even been used to
suggest OQDs in 2D materials, but these studies have remained targeted on a
few defects and have not browsed the large elemental space of possible
defects28, 29, 30, 31, 32.

Here, we use first principles high-throughput (HT) computing to build a
database of point defects in WS2 considering all possible substitutional
defects from 57 elements. We use this database to identify a handful of
promising defects and show that the substitution of cobalt on sulfur (CoS) in
WS2 is especially appealing. First principles computations indicate that the
neutral CoS shows several localized levels in the band gap, spin multiplicity,
and a potential for bright telecom emission. This defect is then synthesized
_in situ_ , examined with STM/STS, and the measured energy levels confirm and
benchmark the theoretical predictions, which highlights a unique two-level
quantum system.

## High-Throughput Search

A greatly sought-after electronic structure for an OQD involves two localized
defect levels (one occupied, the other unoccupied) well within the band gap33.
This requires a precise matching of defect and band edge levels. Additionally,
the optical transition between these defect levels should be bright and
exhibit large transition dipole moments (TDMs). While having localized defect
levels within the band gap is not in itself necessary for developing OQDs,
this electronic structure has advantages in terms of brightness and robustness
versus temperature34. With the 2.4 eV electronic band gap for WS2, finding
defect levels that are at the same time isolated within the band gap and with
transitions in the telecom or visible range (from 750 meV up to 2 eV) should
be achievable. However, identifying defects that could act as an OQD within
WS2 is challenging.

To search for such a defect, we have built a database with the computed
electronic structure of 757 charged point defects in WS2 considering either
the tungsten (MW) or sulfur (MS) substitution site (see Fig. 1a). All the
elements from the periodic table are used with the exception of rare-earths
and transuranides. We start our screening by computing the relaxed structure
and formation energies of the defect in multiple charge states within Density
Functional Theory (DFT) in the generalized gradient approximation (GGA).
Single-particle energies and band gaps are notoriously underestimated within
DFT and one of the gold standards in defect computation is to use hybrid
functionals such as PBE0 which adds a fraction of Fock exchange to the GGA
functional35, 36.

Recently, we have shown that for 2D materials, using a unique fraction of Fock
exchange for the defect and the host is not adequate and we use here an
approach combining a different amount of Fock exchange for defect levels and
band edges (see Methods)36. The use of hybrid functionals leads to a
significantly higher computational cost and can preclude broad screening.
Here, we accelerate the hybrid computation by fixing the wave function from
DFT and applying the hybrid functional Hamiltonian from PBE0 37. This single-
shot PBE0 approach (or PBE00) is similar to the single-shot $GW$
($G_{0}W_{0}$) approach and enables single-particle energy predictions that
are much improved compared to DFT at a minimal computational overhead, which
we have used for defects in silicon 37.

Our computational database includes formation energies, spin state, and
single-particle electronic energy levels for all the possible charged defects.
It also contains the TDMs between these single-particle levels indicating
optical transition brightness. We use this database to search for attractive
OQD candidates. We only consider defects with a charge state that is stable
within a certain Fermi level (EF) and with transitions between defect levels
localized within the band gap. For these defects, we evaluate their single-
particle excitation energies and TDMs. Fig. 1b shows the TDM versus excitation
energy for all defects. We differentiate between MS and MW defects as well as
singlets and multiplets. Few defects show high brightness (with a TDM of at
least 2.5 D) and an excitation energy within the telecom or visible range
(single-particle excitation energy > 750 meV) (see Supplementary Table 1 and
Supplementary Figs. 1 and 2 for a full list with their single-particle
levels). We identify a series of potential singlet OQDs that could act as
single-photon emitters and are formed through the substitution of W with a
main group element: Sb${}_{\rm W}^{-1}$, P${}_{\rm W}^{-1}$, Pb${}_{\rm
W}^{-2}$, N${}_{\rm W}^{-1}$, and C${}_{\rm W}^{-2}$. Only two transition
metal defects appear as promising singlets: OsW and TiS. Spin multiplets that
are of greater technological interest only appear through sulfur substitution:
Co${}_{\rm S}^{0}$, Fe${}_{\rm S}^{0}$, Zn${}_{\rm S}^{0}$, Si${}_{\rm
S}^{-1}$, and W${}_{\rm S}^{+1}$ except for Ru${}_{\rm W}^{0}$. The WS defect
has been suggested as an OQD by Tsai et al. as well, but in the zero charge
state28. Notably, common substitutional defects to tungsten in WS2: Re, V, Nb,
Mo, and Cr, do not show an adequate electronic structure (see Supplementary
Fig. 3) 38, 39, 40, 41, 19. They all have at most one level in the band gap of
the substitutional $d$ orbital character that is slightly above the valence
band edge (V) or below the conduction band (Cr and Re). They are only
excitable optically through a transition between a localized defect state and
a delocalized band level forming a bound exciton34. Our findings agree with
experimental results from photoluminescence or STS on MW defects40, 41, 19,
38.

Fig. 1: Two-Level Quantum Defect Screening in WS2. a Two defect configurations
that are considered in this work: substitution on W site (MW) and on S site
(MS). b Transition dipole moment vs. single-particle excitation energy at the
single-shot PBE0. The marker and color scheme stand for the defect structure
and whether the ground state is singlet or not. Each point stands for a charge
defect that is thermodynamically stable within a certain EF range in the band
gap, and with electronic structures that possess two localized defect levels
within the band gap, as shown in the inset.

While substitutional transition metals on W sites are easy to synthesize, our
screening results show that this is not the most promising approach for OQD
discovery. All our candidate transition metal OQDs except Ru${}_{\rm W}^{0}$
and Os${}_{\rm W}^{0}$ show up instead as MS. Fig. 2a shows the different
electronic structure for MW and MS in a molecular orbital diagram picture when
M is a transition metal42. For both substitutions, the $d$ orbitals of the
defect mix with either sulfur (MW) or tungsten (MS) forming bonding and anti-
bonding states separated by $\Delta_{AB}$. Additionally the different $d$
orbitals are split in three groups: ($d_{xz}$,$d_{yz}$),
($d_{xy}$,$d_{y^{2}-x^{2}}$) and $d_{z^{2}}$ with an energy $\Delta d$
according to crystal field splitting theory. For the sake of simplicity, we
assume here a C3v and D3h point group respectively for the MS and MW defects,
where lower symmetry through Jahn-Teller distortions are also possible. We
performed bonding analysis, and determined density of states, for all 3$d$
transition metal defects (Supplementary Fig. 4) and we observed a smaller
splitting between bonding and anti-bonding states for MS versus MW
($\Delta_{AB}$). This can be rationalized by the different atomic positions
for sulfur and tungsten orbitals. Additionally, the splitting between $d$
orbitals ($\Delta d$) is higher for MS versus MW. Fig. 2b shows the positions
of bonding and anti-bonding molecular orbitals across the 3$d$ series (MW
(blue) and MS (yellow) in the neutral charge state), where 3$d$ atomic
orbitals shift to lower energy from Ti to Cu. MS substitutions show clear
advantage in terms of a smaller $\Delta_{AB}$ and larger $\Delta_{d}$, which
leads to $d$-$d$ transitions in the telecom or visible range and enables more
potential for OQDs with two levels localized in the band gap.

Fig. 2: Molecular orbital trend within the 3$d$ transition metal series for MS
and MW defects. a The molecular orbital diagram shows the splitting between
anti-bonding and bonding state ($\Delta$AB) as well as the splitting with $d$
orbitals ($\Delta$d) for a typical MW and MS defect. b A schematic of the
bonding and anti-bonding state for different 3$d$ transition metals in MW
(blue) and MS (yellow) positions.

## Candidates

While our analysis shows that within gap $d$-$d$ transitions are more likely
in MS and rationalizes why there are still differences between MS defects.
Fig. 1 shows that Co${}_{\rm S}^{0}$ is by far the most attractive OQD
considering its non-singlet spin multiplicity, its large excitation energy,
and transition dipole moment. We compute the electronic structure and
formation energy for this Co${}_{\rm S}^{0}$ defect within full-fledged PBE0
computations including structural relaxation and self-consistency. We plot the
defect formation energy for different charge states of CoS versus EF in Fig.
3a. The defect is stable in its zero charge state spanning a large EF range.
The two thermodynamic charge transition levels correspond to ($+1$/$0$) at 0.4
eV above the valence band maximum (VBM) and ($0$/$1-$) at 0.4 eV below the
conduction band minimum (CBM).

Fig. 3: Thermodynamic charge transition levels and electronic structure of
CoS. a Formation energy of CoS as a function of Fermi level. The charge
transition levels are referenced to the band-edge positions of pristine WS2 as
obtained with PBE0 incorporating 22% of Fock exchange PBE0(0.22). b Orbital
diagram of the localized defect states for neutral CoS. Resonant states within
the valence band and conduction band manifolds are not depicted. The occupied
(unoccupied) states are shown by the filled (empty) rectangles, the height of
which indicates the degree of dispersion. The band-edge positions refer to
those of the pristine WS2 obtained with PBE0(0.22). Energies are referenced to
the vacuum level. SOC is not taken into account for the localized defect
states. c Top view of the charge density (in blue) for the three Co${}_{\rm
S}^{0}$ defect states as indicated in b. The isovalue is 0.001 $e$/Å3.

The electronic structure of the neutral Co${}_{\rm S}^{0}$ is shown in Fig.
3b. A full description of the electronic structure for all three charge states
is given in Supplementary Fig. 5. The neutral defect undergoes a Jahn-Teller
distortion towards the Cs symmetry. While there is significant mixing with the
host, the projection on the Co-3$d$ orbitals is provided in Fig. 3b with the
wavefunctions illustrated in Fig. 3c. The defect shows occupied $d_{xy}$ and
$d_{z^{2}}$ states well within the band gap that can be excited to the
unoccupied $d_{x^{2}-y^{2}}$, $d_{xz}$ or $d_{yz}$ levels that are also below
the conduction band. The lowest energy transition is between the $d_{z^{2}}$
and $d_{x^{2}-y^{2}}$ states and sits at a 1.5 eV difference in single-
particle energies and shows a TDM of 1.2 D. All these values are obtained from
full PBE0 but confirm the prediction from our screening at the single-shot
level. The zero-phonon lines (ZPL) associated with this transition can be
computed in constrained DFT by imposing the occupation of the unoccupied
$d_{x^{2}-y^{2}}$ state and relaxing the structure. We compute a ZPL of 0.96
eV, well within the telecom region. Transition from the $d_{z^{2}}$ to the
next orbital ($d_{xy}$) is significantly higher with a ZPL of 1.18 eV (and a
TDM of 3.0 D). All these results confirm the interest of the neutral Co
substitutional defect as it combines emission in the telecom and spin
multiplicity.

## CoS Fabrication and Characterization

In order to benchmark this novel screening approach, we create and
characterize the CoS defect. Comparisons between the specific energy levels
and effective orbital symmetries enable a direct comparison with the HT
screening approach and first-principles computations in general. In order to
fabricate the CoS defect in WS2, we make use of a unique experimental workflow
inside a low temperature and ultrahigh vacuum (UHV) scanning probe microscope
(SPM) that is shown in Fig. 4a-c. Sulfur vacancies (VS) within otherwise as-
grown WS2 are created by resistively heating the sample and, in tandem,
exposing it to a low incidence angle Ar+ sputtering beam (Fig. 4a)43. This
technique produces a high density of VS available for functionalization and
subsequent reactivity. As adsorbed cobalt has been shown to be unstable on
pristine TMD systems, such as MoS2 and WS2, we are able to make use of
adsorbed instability near the VBM of WS2 (below $-1.3$ V) to systemically
induce diffusion and/or evaporation events with the SPM tip44, 45. A Co
physical vapor deposition apparatus (Fig. 4b) in UHV deposits randomly
adsorbed Co to a defective WS2/MLG/SiC(0001) sample, which is held at liquid
helium temperatures, at submonolayer coverage. The bias over an adsorbed Co
atom can then be ramped towards the tip-induced diffusion energy range to
effectively excite the Co adatom into a VS for CoS defect creation. Fig. 4d-g
shows resulting data at each fabrication workflow step with scanning tunneling
micrographs. Linear defects are identified as one dimensional inversion
domains, which are a result of the VS creation process and has been described
in detail elsewhere43. We then focus on the realization of CoS. STM images,
taken in constant-current mode, over a single Co defect are acquired before a
Co diffusion event and after CoS formation, where the apparent height is
reduced by 0.14 nm (Fig. 4h).

Fig. 4: CoS Defect Formation and Characterization. a The process of forming a
high density of VS, b low-temperature deposition of Co atoms _in situ_ , and c
subsequent placement into a VS with the assistance of the STM probe that is
used to selectively manipulate atoms at voltage ranges below $-1.3$ V is shown
schematically. Corresponding scanning tunneling micrographs that capture
WS2/Gr/SiC(0001) d after defect introduction via Ar+ bombardment and e post Co
deposition are plotted ($I_{tunnel}$ = 30 pA, $V_{sample}$ = $1.2$ V). Scale
bars, 2 nm. STM images f before a voltage excitation and g after Co
substitution within an identified VS are also shown ($I_{tunnel}$ = 30 pA,
$V_{sample}$ = $1.2$ V, Vexcitation = $-2.1$ V). Scale bars, 2 nm. h The
apparent height difference of CoS compared to adsorption atop as-grown WS2 is
measured to be 0.15 nm, taken from linescans across both f and g red
highlighted regions.

In order to investigate the evolved electronic structure with SPM, we make use
of STS and differential conductance mapping, which are representative of the
local density of states (LDOS) over a given defect. Point STS over CoS is
shown in Fig. 5a-b, where in-gap states near $0.36$ eV and $0.47$ eV are
measured. We attribute peak broadening to electronic-phonon coupling, where
effective electron-phonon coupling strength is estimated with a single-mode
Franck-Condon model16. We include multiple phonon modes and additional quanta
of each mode (available for co-excitation) in the description detailed in
Supplementary Fig. 6 to explain dI/dV signal strength and broadening observed
beyond the model approximation. Additionally, a resonance peak is identified
at negative voltages ($-0.84$ $\pm$ 0.06 eV) that is attributed to electronic
charging from the underlying substrate to CoS, which shifts the defect to an
anion state, where an electron is, on average, donated to available CoS defect
levels. Spatially resolved DOS below the charging onset is comparable to that
of the occupied orbitals in the anionic state and to the charge neutral state
above this onset. Fig. 5c-e shows high-resolution differential conductance
image maps that detail electronic orbital densities measured at $-0.9$ eV,
$0.36$ eV, and $0.47$ eV. The LDOS at these energies are further benchmarked
against calculations at the PBE0 level of theory and shown in Fig. 5f-h for
each energy range experimentally measured, where CoS unoccupied orbitals are
hybridized with bonded W atoms and are $\sim$1.5 nm in diameter (see
Supplementary Fig. 7 for simulated STS for charge states presented). We find
strong agreement between experiment and theoretically obtained energy levels
and orbital symmetries, where we can then assign the dI/dV peak at 0.36 eV to
predominately $d_{x^{2}-y^{2}}$ orbital density, and the peak measured at 0.47
eV to a mixing of $d_{yz}$ and $d_{xz}$ orbitals at the CoS charge neutral
state. The peak at $-0.84$ eV is attributed to the CoS charging (to Co${}_{\rm
S}^{-1}$), and is discussed in further detail below. Quantitatively, the
$d_{x^{2}-y^{2}}$ state is experimentally 0.64 eV below the CBM while theory
predicts a level $0.5$ eV below the CBM, indicating a good agreement.

Fig. 5: Experimental and Simulated CoS Scanning Tunneling Spectroscopy. a STS
spectra recorded on a CoS defect and the as-grown WS2 monolayer on graphene
($V_{modulation}$ = 5 mV). b In-gap states identified are located at peak
maxima of $0.36$ eV and $0.47$ eV, each with a full-width half maximum near
0.045 eV. Differential conductance (dI/dV) imaging maps over the defect are
depicted at c $-0.9$ eV, d $0.373$ eV, and e $0.486$ eV ($V_{modulation}$ = 5
mV), showing CoS orbital geometries. Scale bars, 0.25 nm. f-h Simulated STS
maps using PBE0 over CoS orbitals identifying energy range densities near
experimentally measured values. Scale bars, 0.25 nm. A charging peak is
identified in a, where the i lowest unoccupied Co${}_{\rm S}^{0}$ state
becomes j resonant with the EF of the substrate and an electron is donated to
produce the Co${}_{\rm S}^{-1}$ defect. Both c and f are representative of the
Co${}_{\rm S}^{-1}$ orbital densities collected at the specified energy (the
charging ring onset in c is removed for clarity).

We attribute the sharp peak at $-0.84$ eV to a charging process of the neutral
cobalt to the anionic Co${}_{\rm S}^{-1}$ state. This charging is due to the
localized tip-induced band bending process has been described in the
literature on similar systems46, 16. The CoS lowest unoccupied state is
occupied at adequate negative voltages and alters the CoS charge state making
it anionic, detailed in Fig. 5i-j. The neutral/anionic charge transition level
is computed to be around 2.1 eV above the VBM (see Fig. 3) which is close to
the charge transition level for VS (see Supplementary Fig. 13) for which a
charging peak at a similar position is observed for the same type of sample16.
The charging peak near $-0.84$ eV varies spatially as the bias is ramped to
more negative values: the radius of the ring around the defect center
increases. In order to increase STS statistics, we perform an autonomous
hyperspectral experiment over CoS (see Supplementary Notes 1 and 2 in addition
to Supplementary Figs. 8-11)20. The charging peak is found to energetically
shift between a minimum of $-0.924$ eV and a maximum of $-0.627$ eV during
point STS measurements, which amounts to a 0.3 eV tip-induced bending range of
available states. This is near the onset of the lowest unoccupied state at
$\sim$0.3 eV above the EF, enabling CoS to behave as an electron acceptor.
Spatially-resolved charging ring formation as a function of applied bias is
shown in Supplementary Fig. 12, where linescans taken across differential
conductance maps from the defect center to outside the charging region
highlight a shift to larger distances at more negative voltages. Outside the
defect charging region, the substrate remains in a neutral state, which is
verified with STS around pristine WS2 regions. While the charging process
makes the identification of states closer to the VBM less straightforward, we
note that, inside the charging ring, a state around the cobalt is observed.
This state has the form of a $d_{z^{2}}$ orbital as expected from the computed
LDOS of the neutral cobalt defect in that energy range (Fig. 5c and
Supplementary Fig. 7b). The better comparison is with the LDOS of the
Co${}_{\rm S}^{-1}$ as the defect should be charged within the ring. Theory
predicts a reorganization of orbitals, an upward shift of the $d_{z^{2}}$ and
a change of symmetry going from Cs to C3v when CoS becomes negatively charged
(see Supplementary Fig. 5). From this picture, we expect the $d_{z^{2}}$ state
for Co${}_{\rm S}^{-1}$ to be 1 eV lower than the $d_{x^{2}-y^{2}}$ state from
Co${}_{\rm S}^{0}$. We found experimentally a value of 1.3 eV. If there is an
upward shift of $d_{z^{2}}$ when charged, it is smaller in experiment than in
theory. This discrepancy could come from the influence of the dielectric
environment of the graphene/SiC contacts that is not modeled in our WS2 system
in vacuum. In any case, next to the $d_{x^{2}-y^{2}}$, $d_{yz}$ and $d_{xz}$
Co state within the band gap, an additional Co $d_{z^{2}}$ state is observed
within the band gap (and 1.3 eV lower than the $d_{x^{2}-y^{2}}$ state)
confirming the theoretical results that Co in WS2 can lead to a two-level
system of great interest as a OQD.

## Conlcusion

We use HT computational screening to search for promising quantum defects in
WS2. Based on a database gathering computed properties for 757 charged defects
in WS2, we identify a handful of promising quantum defects with high
brightness and in-gap defect states compatible with optical emission in the
telecom or visible range. We fabricate the Co${}_{\rm S}^{0}$ defect, which
shows brightness, a spin-doublet ground state, and a computed ZPL in the
telecom at 0.966 eV, through metal deposition and subsequent sulfur vacancy
substitution by cobalt with an STM tip. STM and STS analysis indicates cobalt-
related defect states within the band gap confirming the computational
prediction and the interest of CoS as an OQD.

Our HT data indicates that fundamental electronic structure reasons make
transition metal substitution on sulfur sites more likely to lead to a OQD
with in-gap defect states that could emit in the telecom or visible than for
the tungsten substitution. This motivates more efforts in the community along
that direction. The fabrication process and HT computational screening used to
identify CoS highlight the capability of combining HT screening and advanced
synthesis techniques to identify and realize new OQDs. This can be performed
across a wide range of atomic species within 2D materials and other hosts with
many yet to be experimentally realized, which can be executed for a number of
different desired material properties, e.g., from catalysis to QIS.

## Materials and Methods

### Scanning probe microscopy (SPM) measurements

All measurements were performed with a Createc GmbH scanning probe microscope
operating under ultrahigh vacuum (pressure < 2$\times$10-10 mbar) at liquid
helium temperatures ($T<6$ K). Either etched tungsten or platinum iridium tips
were used during acquisition. Tip apexes were further shaped by indentations
onto a gold substrate for subsequent measurements taken over a defective
substrate. STM images are taken in constant-current mode with a bias applied
to the sample. STS measurements were recorded using a lock-in amplifier with a
resonance frequency of 683 Hz and a modulation amplitude of 5 mV.

### Sample preparation

Monolayer islands of WS2 were grown on graphene/SiC substrates with an ambient
pressure CVD approach. A graphene/SiC substrate with 10 mg of WO3 powder on
top was placed at the center of a quartz tube, and 400 mg of sulfur powder was
placed upstream. The furnace was heated to 900 ∘C and the sulfur powder was
heated to 250 ∘C using a heating belt during synthesis. A carrier gas for
process throughput was used (Ar gas at 100 sccm) and the growth time was 60
min. The CVD grown WS2/MLG/SiC was further annealed _in vacuo_ at 400 ∘C for 2
hours. Cobalt was deposited at a pressure of $1\times 10^{-9}$ mbar for 60
seconds with the sample held at 5 K.

### Neural network and Gaussian process implementation

The acquisition software used for autonomous experimentation was gpSTS, which
is a library for autonomous experimentation for scanning probe microscopy47,
20. An Intel Xeon E5-2623 v3 CPU with 8 cores and 64 GB of memory combined
with a Tesla K80 with 4992 CUDA cores was used for training the neural
network. Training data for WS2 and VS was combined with CoS spectra obtained
from an extended autonomous run.

### First-principles calculations

We considered 57 elements that could substitute for W and S in the
construction of a WS2 quantum defect database, as highlighted in the periodic
table in Supplementary Fig. 1. This collection covers the majority of the
elements except the rare-earth elements and noble gases. All defect
computations were performed at DFT level using automatic defect workflows that
are implemented in atomate software package48, 49, 50. The defect structure
generations and the formation energy computations are performed using PyCDT.
The DFT calculations were performed using Vienna Ab-initio Simulation Package
(VASP)51, 52 and the projector-augmented wave (PAW) method53 with the Perdew-
Burke-Ernzerhof (PBE) functional54. Each charged defect is simulated in a
144-atom orthorhombic supercell and with a vacuum of approximately 14 Å. A
plane-wave basis energy cutoff of 520 eV was used and the Brillouin zone is
sampled using $\Gamma$ point only. The defect structures were optimized at a
fixed volume until the forces on the ions are smaller than 0.01 eV/Å. The
charge states of each defect are determined by considering all the oxidation
states of the elements documented in the ICSD database49 and taking into
account the formal charges in WS2 (W4+ and S2-). The total energy of the
charged defects were further corrected to overcome the finite-size effect
using the method of Komsa et al.55, 56 as implemented in slabcc57.

The above procedures generated overall 757 substitutional charged defects in
monolayer WS2. Based on the defect formation energy, we first identified 260
charged defects that are thermodynamically stable, meaning their charge states
are accessible in a certain EF range. Among these stable defects, we further
search for the ones that possess two in-gap, localized levels that would
enable the optical intra-defect transition. The localization is defined using
inverse participation ratio (IPR) as detailed below. We considered levels with
IPR larger than 0.05 as localized states (bulk-like states in general have IPR
smaller than 0.01 in WS2). This trimmed down the list to 143 candidates, among
which 112 have non-singlet ground states. The classification of singlets and
multiplets is based on the electronic structure of the defect. In this case,
the singlets and multiplets refer to the total magnetic quantum number of the
unpaired electrons. Thus, defects with all electrons paired are classified as
singlet, while those with one or two paired electrons are classified as
doublet, triplet, etc. We note that due to limitations of Kohn-Sham (KS) DFT
and possibility of spin contamination for spin-polarized systems, more
powerful methods such as spin-flip Bethe-Salpeter are required in general to
rigorously determine the total spin $S$58. Finally, we screened out the ones
that would emit at telecom wavelength with reasonable brightness. The emission
wavelength is approximated using the single-particle KS energy difference
using the single-shot PBE0 incorporating 7% of Fock exchange. We refrained
from applying potential corrections at this stage as the KS energy difference
is largely unaffected by the electrostatic finite-size effect. The brightness
of the optical transition is approximated by the transition dipole moment
(TDM) as detailed below. To search for the most relevant transitions, we
consider the transitions that give the smallest energy difference, while also
allowing an energy window of up to 100 meV to take into account the errors and
band degeneracy. We then identified the transition with the largest TDM as the
most relevant transition. The above procedures recommend 17 non-singlet
candidates that emit at least 750 meV with a TDM of 3 D, as shown in
Supplementary Table 1.

The localization of an orbital is described using the IPR. For a given KS
State, the IPR is evaluated based on the probabilities of finding an electron
with an energy $E_{i}$ close to an atomic site $\alpha$59, 60, 61:

$\chi\left(E_{i}\right)=\frac{\sum_{\alpha}\rho_{\alpha}^{2}\left(E_{i}\right)}{\left[\sum_{\alpha}\rho_{\alpha}\left(E_{i}\right)\right]^{2}},$
(1)

where the summation runs over all atomic sites $\alpha$. The participation
ratio $\chi^{-1}$ stands for the number of atomic sites that confine the wave
function. Thus, a larger (smaller) IPR indicate a localized (delocalized)
state. IPR is unitless ranging between 0 to 1. We computed IPR using VASP
PROCAR. The optical transition dipole moment was evaluated by the Pyvaspwfc
code based on the single-particle wavefunction calculated at the PBE level62.
The transition dipole moment is written as:

$\boldsymbol{\mu}_{k}=\frac{\mathrm{i}\hslash}{\left(\epsilon_{\mathrm{f},k}-\epsilon_{\mathrm{i},k}\right)m}\left\langle\psi_{\mathrm{f},k}|\mathbf{p}|\psi_{\mathrm{i},k}\right\rangle,$
(2)

where $\hslash$ is the Planck constant, $\epsilon_{\mathrm{i},k}$ and
$\epsilon_{\mathrm{f},k}$ are the eigenvalues of the initial and final states,
$m$ is the electron mass, $\psi_{\mathrm{i}}$ and $\psi_{\mathrm{f}}$ are the
initial and final wavefunctions, and $\mathbf{p}$ is the momentum operator.

For selected substitutional defects, we carried out the fully self-consistent
hybrid functional (PBE0) calculations including structural relaxations. In
line with the single-shot PBE0 calculations and following previous work36, we
described the defect levels using the mixing parameter $\alpha=0.07$ for the
Fock exchange, which generally satisfies the Koopmans’ condition for localized
defects in monolayer WS2. On the other hand, we used $\alpha=0.22$ for the
pristine WS2 to determine the band-edge position. The alignment of defect
levels with respect to the band edges was then achieved through the vacuum
level which serves a common reference level. Spin-orbit coupling (SOC) is
taken into account unless otherwise specified. We used a planewave cutoff
energy of 400 eV and a $2\times 2\times 1$ $\mathbf{k}$-point mesh for ground-
state calculations. The zero-phonon line was assessed using a single $\Gamma$
point by imposing occupation constraints (constrained DFT 27). For charged
defects, the total energies are subject to finite-size effects and were
corrected by the method of Komsa et al.55, 56 as implemented in slabcc57,
whereas the single-particle KS levels were corrected by the potential
correction scheme (scpc) of Chagas da Silva et al.63. The simulated STM images
were plotted at a constant height of 3.5 Å above the surface using the
STM-2DScan package64 based on the Tersoff-Hamann theory65.

### Data availability

The data that support the findings of this study are available from the
corresponding authors on reasonable request.

### Code availability

The code used for the findings of this study are available from the
corresponding authors on reasonable request.

### Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences in Quantum Information Science under Award Number DE-
SC0022289. This work was supported as part of the Center for Novel Pathways to
Quantum Coherence in Materials, an Energy Frontier Research Center funded by
the U.S. Department of Energy, Office of Science, Basic Energy Sciences. Work
was performed at the Molecular Foundry and at the Advanced Light Source
supported by the Office of Science, Office of Basic Energy Sciences, of the
U.S. Department of Energy under contract no. DE-AC02-05CH11231. S.K and J.A.R.
acknowledge support from the National Science Foundation Division of Materials
Research (NSF-DMR) under awards 2002651 and 2011839. N.K. and D.A.S
acknowledge support from the National Science Foundation award DMR-2144317,
and the Merced nAnomaterials Center for Energy and Sensing (MACES), a NASA-
funded research and education center, under award NNH18ZHA008CMIROG6R. B.A.B.
was supported by the U.S. Department of Energy, Office of Science, Basic
Energy Sciences, CTC and CPIMS Programs, under Award DE-SC0019053. This
research used resources of the National Energy Research Scientific Computing
Center, a DOE Office of Science User Facility supported by the Office of
Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231
using NERSC award BES-ERCAP0020966. Additional computational resources were
provided by the Multi-Environment Computer for Exploration and Discovery
(MERCED) cluster at UC Merced, funded by National Science Foundation Grant No.
ACI-1429783.

### Author Contributions

G.H., A.W.-B, J.C.T., S.G., A. R. conceived the overall project. W.C., Y.X.,
B.A.B., W. C., S.G., D.A.S., G.-M.R., and G.H. performed the theoretical
simulations along with compiling the database for quantum defect search.
J.C.T, J.Z., and A.R. performed STM/STS experiments and subsequent analysis
with support from E.S.B., A.S., D.F.O., E.R., A.R., and A.W.-B. J.C.T. and
M.M.N. implemented autonomous experimentation. Z.Y., D.Z., S.K., J.A.R., and
M.T. carried out sample growth. All authors discussed the results and
contributed towards the manuscript.

### Competing interests

The authors declare that they have no competing interests.

## References

  * 1 Maletinsky, P. et al. A robust scanning diamond sensor for nanoscale imaging with single nitrogen-vacancy centres. Nat. Nanotechnol. 7, 320 (2012).
  * 2 Bradley, C. E. et al. A Ten-Qubit Solid-State spin register with quantum memory up to one minute. Phys. Rev. X 9, 031045 (2019).
  * 3 Pompili, M. et al. Realization of a multinode quantum network of remote solid-state qubits. Science 372, 259 (2021).
  * 4 Wolfowicz, G. et al. Quantum guidelines for solid-state spin defects. Nat. Rev. Mater. 6, 906 (2021).
  * 5 Bassett, L. C., Alkauskas, A., Exarhos, A. L., & Fu, K.-M. C. Quantum defects by design. Nanophotonics 8, 1867 (2019).
  * 6 Atatüre, M., Englund, D., Vamivakas, N., Lee, S.-Y., & Wrachtrup, J. Material platforms for spin-based photonic quantum technologies. Nat. Rev. Mater. 3, 38 (2018).
  * 7 Sukachev, D. D. et al. Silicon-Vacancy spin qubit in diamond: A quantum memory exceeding 10 ms with Single-Shot state readout. Phys. Rev. Lett. 119, 223602 (2017).
  * 8 Bourgeois, E., Gulka, M., & Nesladek, M. Photoelectric detection and quantum readout of nitrogen-vacancy center spin states in diamond. Adv. Opt. Mater. 8, 1902132 (2020).
  * 9 Stern, H. L. et al. Room-temperature optically detected magnetic resonance of single defects in hexagonal boron nitride. Nat. Commun. 13, 618 (2022).
  * 10 Kianinia, M., Xu, Z.-Q., Toth, M., & Aharonovich, I. Quantum emitters in 2D materials: emitter engineering, photophysics, and integration in photonic nanostructures. Appl. Phys. Rev. 9, 011306 (2022).
  * 11 Lin, Z. et al. Defect engineering of two-dimensional transition metal dichalcogenides. 2D Mater. 3, 022002 (2016).
  * 12 Manzeli, S., Ovchinnikov, D., Pasquier, D., Yazyev, O. V., & Kis, A. 2D transition metal dichalcogenides. Nat. Rev. Mater. 2, 1 (2017).
  * 13 Li, C. et al. Engineering graphene and TMDs based van der waals heterostructures for photovoltaic and photoelectrochemical solar energy conversion. Chem. Soc. Rev. 47, 4981 (2018).
  * 14 Ugeda, M. M. et al. Observation of topologically protected states at crystalline phase boundaries in single-layer WSe2. Nature Communications, vol. 9, p. 3401, Aug 2018.
  * 15 Schuler, B. et al. Electrically driven photon emission from individual atomic defects in monolayer WS2. Sci. Adv. 6, eabb5988 (2020).
  * 16 Schuler, B. et al. Large spin-orbit splitting of deep in-gap defect states of engineered sulfur vacancies in monolayer WS2. Phys. Rev. Lett. 123, 076801 (2019).
  * 17 Montblanch, A. R. P., Barbone, M., Aharonovich, I., Atatüre, M., & Ferrari, A. C. Layered materials as a platform for quantum technologies. Nat. Nanotechnol. 18, 555 (2023).
  * 18 Ye, M., Seo, H., & Galli, G. Spin coherence in two-dimensional materials. npj Comput. Mater. 5, 44 (2019).
  * 19 Schuler, B. et al. How substitutional point defects in two-dimensional WS2 induce charge localization, spin-orbit splitting, and strain. ACS Nano 13, 10520 (2019).
  * 20 Thomas, J. C. et al. Autonomous scanning probe microscopy investigations over WS2 and Au{111}. npj Comput. Mater. 8, 99 (2022).
  * 21 Barja, S. et al. Identifying substitutional oxygen as a prolific point defect in monolayer transition metal dichalcogenides. Nat. Commun. 10, 3382 (2019).
  * 22 Cochrane, K. A. et al. Spin-dependent vibronic response of a carbon radical ion in two-dimensional WS2. Nat. Commun. 12, 7287 (2021).
  * 23 Stolz, S. et al. Layer-dependent schottky contact at van der waals interfaces: V-doped WSe2 on graphene,” npj 2D Mater. Appl. 6, 66 (2022).
  * 24 Gali, A. Ab initio theory of the nitrogen-vacancy center in diamond. Nanophotonics 8, 1907 (2019).
  * 25 Freysoldt, C. et al. First-principles calculations for point defects in solids. Rev. Mod. Phys. 86, 253 (2014).
  * 26 Dreyer, C. E., Alkauskas, A., Lyons, J. L., Janotti, A. , & Van de Walle, C. G. First-principles calculations of point defects for quantum technologies. Annu. Rev. Mater. Res. 48, 1 (2018).
  * 27 Ivády, V., Abrikosov, I. A., & Gali, A. First principles calculation of spin-related quantities for point defect qubit research. npj Comput. Mater.4, 76 (2018).
  * 28 Tsai, J.-Y., Pan, J., Lin, H., Bansil, A., & Yan, Q. Antisite defect qubits in monolayer transition metal dichalcogenides. Nat. Commun. 13, 492 (2022).
  * 29 Frey, N. C., Akinwande, D., Jariwala, D., & Shenoy, V. B. Machine learning-enabled design of point defects in 2D materials for quantum and neuromorphic information processing. ACS Nano 14, 13406 (2020).
  * 30 Ping, Y. & Smart, T. J. Computational design of quantum defects in two-dimensional materials. Nat. Comput. Sci. 1, 646 (2021).
  * 31 Narang, P., Ciccarino, C. J., Flick, J., & Englund, D. Quantum materials with atomic precision: artificial atoms in solids: ab initio design, control, and integration of single photon emitters in artificial quantum materials. Adv. Func. Mater. 29 1904557 (2019).
  * 32 Li, S., Thiering, G., Udvarhelyi, P., Ivády, V., & Gali, A. Carbon defect qubit in two-dimensional WS2. Nat. Commun. 13, 1 (2022).
  * 33 Gupta, S., Yang, J.-H., & Yakobson, B. I. Two-level quantum systems in two-dimensional materials for single photon emission. Nano Lett. 19, 408 (2019).
  * 34 Xiong, Y., Mathew, M., Griffin, S. M., Sipahigil, A., & Hautier, G. Midgap state requirements for optically active quantum defects. Preprint available at https://arxiv.org/abs/2302.10767 (2023).
  * 35 Freysoldt, C. et al. First-principles calculations for point defects in solids. Rev. Mod. Phys. 86, 253 (2014).
  * 36 Chen, W., Griffin, S. M., Rignanese, G.-M., & Hautier, G. Nonunique fraction of Fock exchange for defects in two-dimensional materials. Phys. Rev. B106, L161107 (2022).
  * 37 Xiong, Y. et al. High-throughput identification of spin-photon interfaces in silicon. Preprint available at https://arxiv.org/abs/2303.01594 (2023).
  * 38 Loh, L. et al. Impurity-induced emission in Re-doped WS2 monolayers. Nano Lett. 21, 5293 (2021).
  * 39 Qin, Z. et al. Growth of nb-doped monolayer WS2 by liquid-phase precursor mixing. ACS Nano 13, 10768 (2019).
  * 40 Zhang, F. et al. Monolayer vanadium-doped tungsten disulfide: A room-temperature dilute magnetic semiconductor. Adv. Sci. 7, 2001174 (2020).
  * 41 Han, A. et al. One-step synthesis of single-site vanadium substitution in 1T-WS2 monolayers for enhanced hydrogen evolution catalysis. Nat. Commun. 12, 709 (2021).
  * 42 Pike, N. A. et al. Origin of the counterintuitive dynamic charge in the transition metal dichalcogenides. Phys. Rev. B 95, 201106 (2017).
  * 43 Rossi, A. et al. WS2 Band Gap Renormalization Induced by Tomonaga Luttinger Liquid Formation in Mirror Twin Boundaries. Preprint available at https://arxiv.org/abs/2301.02721 (2023).
  * 44 Kandel S. A. & Weiss, P. S. Binding and mobility of atomically resolved cobalt clusters on molybdenum disulfide. J. Phys. Chem. B 105, 8102 (2001).
  * 45 Tang, W. et al. Identically sized Co quantum dots on monolayer WS2 Featuring Ohmic Contact. Phys. Rev. Applied 13, 024003 (2020).
  * 46 Teichmann, K. et al. Controlled charge switching on a single donor with a scanning tunneling microscope. Phys. Rev. Lett. 101, 076103 (2008).
  * 47 Noack, M. M. et al. gpCAM. https://github.com/lbl-camera/gpCAM (2022).
  * 48 Mathew, K. et al. Atomate: A high-level interface to generate, execute, and analyze computational materials science workflows. Comput. Mater. Sci. 139, 140 (2017).
  * 49 Ong, S. P. et al. Python materials genomics (pymatgen): A robust, open-source python library for materials analysis,” Comput. Mater. Sci. 68, 314 (2013).
  * 50 Jain, A. et al. Commentary: The Materials Project: A materials genome approach to accelerating materials innovation. APL Mater. 1, 11002 (2013).
  * 51 Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 6, 15 (1996).
  * 52 Kresse G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169 (1996).
  * 53 Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 50, 17953 (1994).
  * 54 Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865 (1996).
  * 55 Komsa, H.-P. & Pasquarello, A. Finite-size supercell correction for charged defects at surfaces and interfaces. Phys. Rev. Lett.110, 095505 (2013).
  * 56 Komsa, H.-P., Berseneva, N., Krasheninnikov, A. V., & Nieminen, R. M. Charged point defects in the flatland: Accurate formation energy calculations in two-dimensional materials. Phys. Rev. X 4, 031044 (2014).
  * 57 Farzalipour Tabriz, M., Aradi, B., Frauenheim, T., & Deák, P. SLABCC: Total energy correction code for charged periodic slab models. Comput. Phys. Commun. 240, 101 (2019).
  * 58 Barker, B. A. & Strubbe, D. A. Spin-flip Bethe-Salpeter equation approach for ground and excited states of open-shell molecules and defects in solids. Preprint available at http://arxiv.org/abs/2207.04549. (2022).
  * 59 Wegner, F. Inverse participation ratio in 2+$\epsilon$ dimensions. Zeitschrift für Phys. B Condens. Matter 36, 209 (1980).
  * 60 Pashartis C. & Rubel, O. Localization of electronic states in III-V semiconductor alloys: A comparative study. Phys. Rev. Applied 7, 064011 (2017).
  * 61 Konstantinou, K., Mocanu, F. C., Lee, T.-H., & Elliott, S. R. Revealing the intrinsic nature of the mid-gap defects in amorphous Ge2Sb2Te5. Nat. Commun. 10, 3065 (2019).
  * 62 Zheng, Q. Vasp band unfolding. https://github.com/QijingZheng/VaspBandUnfolding (2018).
  * 63 Chagas da Silva, M. et al. Self-consistent potential correction for charged periodic systems. Phys. Rev. Lett. 126, 076401 (2021).
  * 64 Leung, S. Stm 2d scan. https://https://github.com/ShuangLeung/STM_2DScan (2020).
  * 65 Tersoff, J. & Hamann, D. R. Theory of the scanning tunneling microscope. Phys. Rev. B 31, 805 (1985).


††thanks: These two authors contributed equally.††thanks: These two authors
contributed equally.

# Urban green space and happiness in developed countries

Oh-Hyun Kwon  Department of Physics, Pohang University of Science and
Technology, Pohang 37673, Republic of Korea.  Inho Hong  Center for Humans and
Machines, Max Planck Institute for Human Development, Berlin 14195, Germany.
Jeasurk Yang  Department of Geography, National University of Singapore,
Singapore 119260, Singapore.  Donghee Yvette Wohn  Department of Informatics,
New Jersey Institute of Technology, Newark, NJ 07103, USA.  Woo-Sung Jung
<EMAIL_ADDRESS>Department of Physics, Pohang University of Science and
Technology, Pohang 37673, Republic of Korea.  Department of Industrial and
Management Engineering, Pohang University of Science and Technology, Pohang
37673, Republic of Korea.  Asia Pacific Center for Theoretical Physics, Pohang
37673, Republic of Korea.  Meeyoung Cha<EMAIL_ADDRESS>Data Science Group,
Institute for Basic Science, Daejeon 34126, Republic of Korea.  School of
Computing, Korea Advanced Institute of Science and Technology, Daejeon 34141,
Republic of Korea.

###### Abstract

Urban green space has been regarded as contributing to citizen happiness by
promoting physical and mental health. However, how urban green space and
happiness are related across many countries of different socioeconomic
conditions has not been explained well. By measuring urban green space score
(UGS) from high-resolution Sentinel-2 satellite imagery of 90 global cities
that in total cover 179,168 km2 and include 230 million people in 60 developed
countries, we reveal that the amount of urban green space and the GDP can
explain the happiness level of the country. More precisely, urban green space
and GDP are each individually associated with happiness; happiness in the 30
wealthiest countries is explained only by urban green space, whereas GDP alone
explains happiness in the 30 other countries in this study. Lastly, we further
show that the relationship between urban green space and happiness is mediated
by social support and that GDP moderates the relationship between social
support and happiness, which underlines the importance of maintaining urban
green space as a place for social cohesion in promoting people’s happiness.

##  I Introduction

The advantages of urban green space for public health and urban planning have
been of great interest in recent years. Green spaces such as parks, gardens,
street trees, riversides, and even private backyards facilitate physical
activity, social events, mental relaxation, and relief from stress and heat,
thereby leading to direct and indirect benefits for mental and physical health
de Vries et al. (2003); Dadvand et al. (2016). Thus, worldwide policy changes
and efforts have been made to build more urban green space to create
sustainable and comfortable living environments UN (2015).

Urban green space and happiness are known to have an implicit positive
correlation. Although this association is still unclear, five pathways through
which greenery might have beneficial effects have been reported: relieving
stress, stimulating physical activity, facilitating social interactions,
generating aesthetic enjoyment, and facilitating a sense of shelter from and
adjustment to environmental stressors de Vries et al. (2013); Dadvand et al.
(2016); Liu et al. (2019). Studies have suggested that the same pathways exist
in numerous countries Dzhambov et al. (2018). Among them, social interaction
facilitation has been confirmed with strong evidence. Studies Maas et al.
(2009); Jennings and Bamkole (2019) have shown that open green space promotes
social cohesion by providing places for social contact; people can naturally
encounter neighbors in local green spaces while walking dogs, gardening, and
having outdoor parties, which enhances community engagement. Moreover, larger
green areas such as parks can hold larger events and activities, enabling
social mixing between communities.

The amount of urban green space can be captured mainly by three kinds of
measurements: qualitative ratings of observers Kweon et al. (1998); de Vries
et al. (2013), national land-use and land-cover database Maas et al. (2006);
MacKerron and Mourato (2013); Alcock et al. (2014), and geographic information
system (GIS) techniques. Among these measurements, GIS techniques are the most
recently developed method. One example is utilizing the normalized difference
vegetation index (NDVI), a vegetation index computed from Landsat series
satellite images (30 m resolution) Beyer et al. (2014); Dzhambov et al.
(2018); Liu et al. (2019). Studies such as by Tsai et al. Tsai et al. (2018)
introduced multiple landscape metrics based on GIS and showed a strong
association between green space and mental health in U.S. metropolitan areas.
These studies assume the distance from an individual’s residence to the
nearest green space has associations with health data Stigsdotter et al.
(2010); Dadvand et al. (2016). The green space level was then measured as the
fraction of areas with NDVI values above a certain threshold (e.g., 0.2 to 0.4
for sparse vegetation and 0.6 for highly dense vegetation) EOS (2019).
However, this method raises the question of how to set an appropriate NDVI
threshold for global cities.

Despite the rich literature on green space’s mental benefits, they still have
limitations as global-scale comparative research. First, the analytical
settings are based on a limited number of Western countries Liu et al. (2019);
most of these studies have been conducted in the United States Beyer et al.
(2014); Tsai et al. (2018) and Europe Dadvand et al. (2016); Dzhambov et al.
(2018). Moreover, only a few are based on multi-country settings that enable
comparative analysis van den Berg et al. (2016). As a result, it is unclear
whether the association between green space and mental health is robust in
developing countries or only in developed countries. The main limitation
arises because there is no global medical dataset providing reliable and
standardized mental health surveys from different countries. Moreover, no
studies have established which green space measurement is appropriate for
analysis across countries. Various methods of measuring green space –
questionnaires, qualitative interviews, satellite images, Google Street View
images, and even smartphone technology Markevych et al. (2017) – still rely on
individual-level measurements (e.g., calculating the greenery level around
residential buildings) and hence are not scalable to the global level.

This paper presents a new way to analyze the effects of green space on
happiness at the planetary scale, incorporate the different countries’
different contexts, and achieve robust results. First, we measure the amount
of urban green space from high-resolution satellite images for different
countries by developing a globally comparable green space metric. Our metric
based on the total NDVI of built-up areas enables this comparison as it does
not require an arbitrary threshold that varies for different regions. It also
overcomes the limitations of official statistics based on national land-use
land cover data that tend to have different criteria by countries and often
include only official parks and open space. Our analysis on high-resolution
(10 m) Sentinel-2 satellite images provides more accurate information of urban
green space than the previous studies on the Landsat series images (i.e., the
resolution of 30 m) Beyer et al. (2014); Dadvand et al. (2016); Dzhambov et
al. (2018); Liu et al. (2019).

Next, this study uses selected happiness scores from the World Happiness
Report Helliwell et al. (2018), which provides reliable and standardized data
on multiple countries’ mental health and allows comparisons among nations. As
happiness is a criterion of emotional well-being, it is interconnected with
mental health. From the perspective that economic studies distinguish between
emotional well-being (happiness) and life satisfaction (life
evaluation)Kahneman and Deaton (2010), we focus on the impact of green space
on emotional happiness. Specifically, we study this relationship in the
developed countries of the highest Human Development Index (HDI), where green
environments in cities are considered more important for well-being.

Using these datasets from satellite imagery, we explore the relationship
between urban green space and happiness globally. Additionally, we identify
conditional indirect effects by national wealth and social support by
employing a moderated mediation regression model on socioeconomic indicators.

##  II Urban green space and happiness in countries

Figure 1: The distributions of urban green space and happiness over the world.
(a) The map of urban green space and happiness in 60 developed countries. The
size and color of circles represent the level of happiness and urban green
space in a country, respectively. The markers are placed on the most populated
cities of each country. (b-d) The histograms of (b) happiness, (c) urban green
space (UGS) and (d) logarithmic GDP per capita (log-GDP). We use the logarithm
of the total NDVI per capita as an indicator of urban green space, and the
logarithm of GDP per capita as a measure of wealth.

We examine the global relationship between urban green space and happiness in
60 developed countries ranked by the Human Development Index. Using the
Sentinel-2 satellite imagery dataset, we define each country’s urban green
space score (UGS) as a logarithmic total vegetation index per capita in the
most populated cities (i.e., those that include at least 10% of the national
population). Among the various vegetation indices available, NDVI Miura et al.
(2019) is used based on the robustness of the results for different tested
indices. The happiness score and the gross domestic product based on
purchasing power parity (GDP (PPP)) per capita of each country are from the
World Happiness Report Helliwell et al. (2018) and the International Monetary
Fund (IMF) estimation IMF (2018), respectively (see the Methods section and
the Supplementary Information for details).

Figure 1(a) shows an overall view of urban green space and the happiness of
countries around the world. This map highlights regional differences in the
green space distribution due to climate; countries near the equator in
tropical climates have relatively high UGS values, while countries located in
the 20-30∘ latitude range have exceptionally low UGS values due to the dry
climate. UGS increases with latitude in higher-latitude regions. On the other
hand, Northern and Western European and North American countries display
relatively large happiness. Western Asian countries also show relatively high
happiness with low UGS value, indicating that the relationship between
happiness and green space is not trivial.

Figure 1(b-d) shows the distribution of happiness, UGS, and log-GDP, and they
all show unimodal distributions with low skewness, which is appropriate for
linear regression analyses. Note that the probability distributions of NDVI
per capita and GDP per capita converge to a normal distribution after
logarithmic scaling. Our comparison of several green space measures shows that
the logarithmic NDVI per capita is most suitable for the following analysis in
terms of its distribution and explanatory power. We hence choose the
logarithmic NDVI per capita as the primary green space indicator in this
research. (see Supplementary Information). We also use the logarithmic GDP per
capita (PPP) (hereinafter referred to as the log-GDP) as a measure of the
wealth of the country, as noted in the Happiness Report Helliwell et al.
(2018).

Countries |  | All |  |  | Lower 30 |  |  | Top 30 |
---|---|---|---|---|---|---|---|---|---
Model | (1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9)
log-GDP | 1.0120*** | - | 1.1319*** | 0.9034** | - | 0.8517* | -0.0809 | - | 0.2581
| (0.6603) |  | (0.6234) | (1.6305) |  | (1.7493) | (1.3559) |  | (1.0314)
UGS | - | 0.1165 | 0.2249*** | - | 0.1497 | 0.0567 | - | 0.2785*** | 0.2946***
|  | (0.3545) | (0.2643) |  | (0.6042) | (0.6051) |  | (0.2313) | (0.2403)
Const | -4.2945** | 5.9007*** | -6.4709*** | -3.3428 | 5.1767*** | -3.0629 | 7.7712** | 5.8110*** | 2.9312
| (6.9672) | (1.4910) | (6.8998) | (16.5490) | (2.6490) | (17.1094) | (14.8065) | (0.9455) | (11.5463)
Adjusted $R^{2}$ | 0.3832 | 0.00123 | 0.4786 | 0.1296 | 0.0012 | 0.1013 | -0.0335 | 0.4457 | 0.4468
Observations | 60 | 60 | 60 | 30 | 30 | 30 | 30 | 30 | 30
Table 1: The regression analysis for happiness, UGS, and log-GDP. The values
denote the regression coefficients and the confidence intervals of each
independent variable with its significance (i.e., ***p¡0.01; **p¡0.05;
*p¡0.1). The regression model (1-3), model (4-6) and model (7-9) are examined
for the data of all countries, the lower 30 countries and the top 30 countries
ranked by GDP, respectively.

As per-country wealth is an important indicator of its citizens’ quality of
life, wealth (i.e., log-GDP) should be considered in analyzing urban green
space and happiness. Our regression analysis finds that UGS, together with
log-GDP, explains happiness. We make new observations from Table 1. Although
UGS is not substantially correlated with happiness in the simple linear
regression (i.e., model (2)), the multilinear model with log-GDP (i.e., model
(3)) has a substantial increase in prediction ability compared to the simple
regression analysis on log-GDP (i.e., model (1)). Therefore, urban green space
adds explanatory power to the correlation between wealth and happiness across
countries. The regression analyses with other green space-variant measures
further confirm this result’s robustness, confirming a substantial increase in
the adjusted R-squared value when including UGS in the regression.
Specifically, UGS based on the logarithmic NDVI per capita shows the best
regression performance (see the Supplementary Information for the results for
the different measures).

##  III Urban green space is effective in rich countries

Our results show that happiness is correlated with urban green space and the
GDP of a country. But, is this green space-happiness effect uniform across all
countries? Previous studies on the marginal effect of income on happiness
suggest that happiness may have a nonlinear relationship with GDP, presumably
showing saturation after a specific GDP — a concept known as the Easterlin
paradox Easterlin et al. (2010). This paradox tells us that increases in
happiness through GDP reach a saturation point, yet what factors promote
happiness beyond the saturation point is unknown.

Figure 2: The effect of GDP on the green-happiness relation. (a, b) The
relations of (a) log-GDP and happiness, and (b) urban green space (i.e., UGS)
and happiness across 60 developed countries. The top 30 and the lower 30
countries ranked by GDP are sized by the population size and colored by red
and black. The dotted lines are the linear fit for each GDP group. (c) Changes
of coefficients between urban green space and happiness for different sets of
GDP rank with increasing window size from top 10 to 60. (d) The rank
correlations between UGS and happiness for the groups of increasing countries
in the GDP rank order.

To test the Easterlin paradox, we repeated the analysis over clusters of
countries grouped by GDP. Figure 2(a) shows a high correlation between GDP and
happiness in the 30 lower-GDP countries (i.e., $\rho=0.40$), whereas the
correlation is no longer evident in the 30 higher-GDP countries (i.e.,
$\rho=-0.04$). These results suggest that economic prosperity (as measured by
GDP) is crucial for people’s happiness but fails to further promote happiness
in rich countries. The GDP appears to reach a happiness-correlation threshold
around the 30th wealthiest country, which corresponds to a GDP of 38,518
dollars. Previous research on the Easterlin paradox has stated that the GDP
per capita can increase happiness until it reaches a certain threshold but
cannot further increase happiness above that threshold. We observe a similar
pattern for wealth and happiness across countries. On the other hand,
happiness in the 30 wealthiest countries is well explained by urban green
space. As shown in Figure 2(b), urban green space is positively correlated
with happiness in the richest countries (i.e., $\rho=0.66,p<0.01$), but this
correlation is not significant in the 30 lower-GDP countries (i.e.,
$\rho=0.19,p=0.32$). Thus, urban green space is a factor that further
increases the happiness of a country after its GDP reaches a certain level.

The regressions for each of the 60 countries ranked by GDP in Table 1 confirm
the individual effects of urban green space and GDP on happiness. GDP is the
only substantial factor explaining happiness in the 30 lower-GDP countries
(models 4-6). In contrast, for the 30 higher-GDP countries, happiness is
explained only by the UGS (7-9). These findings suggest that GDP is critical
for happiness until it reaches a certain GDP threshold (i.e., the Easterlin
paradox), after which urban green space explains happiness better.

The correlation between UGS and happiness also corroborates the effect of UGS
in rich countries. The correlation in Figure 2(d) decreases as more countries
in the decreasing order of GDP are added. The correlation is substantial
(i.e., $\rho$ is approximately 0.8) among the countries excluding the top 30.
Figure 2(c) summarizes the effects of urban green space and GDP that cross
over each other around the 30th wealthiest country. For the top 30 countries,
urban green space has positive coefficients, but the GDP effect is not
significant. These relationships are reversed for less affluent countries.

In summary, economic support seems to promote happiness until the essential
requirements and living standards are met. However, economic success alone
fails to add persistent promotion of happiness. After some level, urban green
space appears to be related to other social factors that can further promote
happiness.

##  IV Urban green space for social cohesion

Our findings highlight urban green space as an indicator that might be
correlated with social factors promoting happiness beyond the achievement of
economic success. The question then arises, which social factors connect urban
green space with happiness? To identify this connection, we first examine the
correlation between UGS and socioeconomic variables reported in the World
Happiness Report: GDP per capita, social support, life expectancy, freedom,
generosity, and perceptions of corruption. Of these six variables, only
“social support” has a significantly positive correlation ($\rho=0.43,p<0.01$)
with UGS as we can see from Fig. 3(a), implying that social support could
mediate between urban green space and happiness. This relationship is
consistent with several existing studies that suggested urban green space as a
place of social cohesion Maas et al. (2009); Jennings and Bamkole (2019). On
the other hand, as indicated by life expectancy, physical health does not
display a significant relationship with green space ( $\rho=0.19,p=0.15$),
contradicting common sense. The regression analysis on happiness with UGS and
six socioeconomic variables also captures the interchangeability of urban
green space and social support (see the Supplementary Information for
details).

Figure 3: The moderated mediation model for UGS, happiness and socioeconomic
indicators. (a) Scatter plot of social support and UGS across countries. (b)
Diagram for the moderated mediation model. The boxes denote the model
variables. Solid black arrows denote a statistically significant relationship
between a pair of variables with the regression coefficient and the p-value
(i.e., ***p¡0.01). The gray dashed arrow represents a non-significant
relationship. Note that the coefficients are calculated with z-scores of the
variables to compare the effect size directly.

Here, we employ a moderated mediation model Preacher et al. (2007) to
characterize the complicated relationships among urban green space, social
support, GDP, and happiness. In moderated mediation models, the moderator
describes a variable’s conditional effect through the interaction term, and
the mediator describes a variable’s indirect effect connecting the other two
variables. Accordingly, moderated mediation models determine the pathway of
directed interactions between multiple variables.

First, we examine the mediation effect of urban green space and social support
as independent variables. The mediation regression model shows that social
support mediates the relationship between urban green space and happiness such
that (1) urban green space improves social support and (2) social support
promotes happiness. The mediation effect is significant (3) only when GDP is
considered in the model. Consequently, our moderated mediation model combining
these three effects, shown in Fig. 3(b), presents the pathway by which green
space affects happiness through social support, given that GDP moderates the
effect of social support. If we describe this relationship in equations,

$H=\beta_{0}+\beta_{1}M+\beta_{2}S+\beta_{3}SM,$ (1)
$S=\beta_{4}+\beta_{5}\ln{G}.$ (2)

where $H$, $M$, $S$, and $G$ represent happiness, log-GDP, social support, and
NDVI per capita, respectively, and the $\beta$ values denote the coefficients
of the regression models (see the Supplementary Information for details).

Our moderated mediation model can be used to estimate the amount of urban
green space required to increase happiness by a certain amount according to

$\Delta
H=\left(\beta_{1}^{\prime}+\beta_{2}^{\prime}M\right)\ln{\frac{G_{f}}{G_{i}}}.$
(3)

In the equation 3, the required ratio of urban green spaces in a country
decreases as its log-GDP increases. The required increase in urban green space
per capita can be estimated for each country based on its current GDP value.
For example, the United States needs an additional 36.1908 NDVI of urban green
space per capita to increase its happiness score by 0.0546. In contrast, 3,416
USD per capita is required to achieve the same increment in happiness. Here,
we used a 0.0546 happiness score as a reference value of $\Delta H$, which is
the average value between happiness ranks. Note that the NDVI per capita is
interpreted as a weighted area of green space, with a unit of $m^{2}$.
Similarly, Qatar needs 0.4981 NDVI per capita or 7,556 dollars per capita, and
South Korea needs 4.1332 NDVI per capita or 2,315 dollars per capita to
achieve the reference happiness score increase.

##  V Discussion

This paper revealed a global relationship between urban green space and
happiness in over 60 countries using high-resolution satellite imagery. Urban
green space has a higher impact in developed countries (i.e., countries with
higher GDPs), which suggests urban green space as a key to promoting happiness
beyond economic success. Our moderated mediation model further elucidates this
relationship as social support mediates the green-happiness relation, and GDP
moderates social support and happiness. This sophisticated model could
estimate additional green space needed to promote happiness for each country.

The current study newly defined the concept of UGS (urban green space score),
which can be used to calculate the amount of green space at any spatial scale
accounting for population density. We compared several green space measures
and proposed to use the logarithmic NDVI per capita as a preferred measure of
UGS. This index was validated through experiments and it makes it possible to
investigate green space at a global level, allowing us to perform cross-
sectional research on green space. Furthermore, the method obtaining UGS can
be utilized to investigate any spatial areas such as blue space (i.e., aquatic
environments such as lake and shore) Foley and Kistemann (2015); Raymond et
al. (2016).

Our findings have multiple policy-level implications. First, public green
space should be made accessible to urban dwellers to enhance social support.
In doing so, one critical aspect is public safety. If public safety in urban
parks is not guaranteed Groff and McCord (2012); Han et al. (2018), its
positive role in social support and happiness may diminish. The meaning of
public safety may change; for example, ensuring biological safety will be a
priority in keeping the urban parks accessible during the COVID-19 pandemic
Ugolini et al. (2020).In fact, the high indoor transmission rate of the virus
Lolli et al. (2020) will increase awareness and importance of open spaces like
urban parks. While some urban parks may be closed during lockdowns, some
reports suggest that viewing them from home could also help relax stress
during the pandemic Hedblom et al. (2019). Second, urban planning of public
green space is needed for both developed and developing countries. While our
findings confirmed a strong impact of urban green space on happiness in
developed countries, the same positive effect holds for developing countries,
albeit to a smaller degree. Furthermore, it is challenging or nearly
impossible to secure land for green space after built-up areas are developed
in cities. Therefore, urban planning for parks and green recovery (new
greening in built-up areas) should be considered in developing economies where
new cities and suburban areas rapidly expand Ewing (2008); Liu et al. (2020).

In addition to the above, recent climate changes can create substantial
volatility in sustaining urban green space. Extreme events such as wildfires,
floods, droughts, and cold waves could endanger urban forests around the world
Allen et al. (2010). On the other hand, global warming could also accelerate
tree growth in cities more than in rural areas due to the urban heat island
effect Pretzsch et al. (2017). In the end, the environmental influence is
bidirectional; urban green spaces affect local climates by reducing carbon
dioxide levels Nowak et al. (2013) and providing a cooling effect inside the
city that indirectly affects people’s well-being. Thus, we need more attention
to predicting climate changes and discovering their impact on public places
since the extreme changes could hamper the benefits of urban green space.

As an exciting future direction, satellite images of higher spatiotemporal
resolutions can be used to compute urban green space scores. This paper
focused on the correlation across countries fixed in time, given the short
span of the Sentinel-2 dataset launched in 2015. A causal analysis could be
done with satellite imagery data for a longer span. Also, our dataset does not
cover all the countries in the world. Fortunately, our observations from the
30 lower-income countries anticipate the substantial effect of GDP in other
developing countries excluded in our analysis. We have analyzed the highest-
resolution public dataset of satellite imagery in this study. However, our
method still has room for application to higher-resolution non-public datasets
such as the household level (less than 10m resolution) available in the
national-scale health dataset Houlden et al. (2017). Since satellite imagery
cannot account for green space inside buildings (such as green walls), future
research could quantify the effect of these mini-scale green spaces using
computer vision Seiferling et al. (2017).

##  VI Methods

###  VI.1 Collecting happiness and remote sensing data

Figure 4: Measuring urban green spaces. (a) Measurement methods to compute the
size of urban green space in each country. First, we find cities occupying
more than 10 percent of the total population in a country. Then, we extract
the built-up area of the cities with Copernicus global landcover data.
Finally, we calculate vegetation indices (e.g., NDVI) within the area using
Sentinel-2 satellite images.(b) Urban green space measured by the UGS (upper
row) and the vegetation ratio (lower row) in four world cities.The red areas
in the upper row indicate vegetation for the NDVI threshold of 0.4. The lower
row images show the adjusted NDVI per capita (i.e., UGS) for every 10m by 10m
pixel.

To identify the relationship between happiness and green space, we use
happiness scores from the World Happiness Report Helliwell et al. (2018) and
the NDVI scores from Sentinel-2 satellite imagery as remote sensing data. The
World Happiness Report from 2018 covered 156 countries. The report provides an
annual survey of how happy citizens perceive themselves to be and ranks the
countries by happiness score. The score is the average of the participants’
responses asked to rate how happy they are on a scale from 0 and 10. While
many socioeconomic indicators (e.g., unemployment and inequality) may affect
happiness, not all of these factors are measured annually across 156
countries. The report instead describes happiness with six primary
socioeconomic indicators: GDP per capita, social support, life expectancy,
freedom to make life choices, generosity, and perceptions of corruption. For
example, the social support variable is based on binary responses (yes/no) on
a Gallup World Poll question: ”If you were in trouble, do you have relatives
or friends you can count on to help you whenever you need them, or not?”

To quantify urban green space in global cities, we use the Sentinel-2 dataset
that provides the highest spatial resolution (10 m) among the publicly
available satellite imagery datasets (e.g., 30 m resolution in Landsat series)
Beyer et al. (2014); Dadvand et al. (2016); Dzhambov et al. (2018); Liu et al.
(2019). With this high resolution, we can identify granular green space,
including street vegetation and home gardens that could not be detected in
other public datasets. When using satellite imagery to detect small
vegetation, it is critical to consider the season in which the images were
obtained Beyer et al. (2014); Dadvand et al. (2016); Dzhambov et al. (2018);
Liu et al. (2019). We use the images from summer: June to September 2018 for
the Northern Hemisphere and December 2017 to February 2018 for the Southern
Hemisphere. Satellite images with below 10% cloud cover were used; when such
images could not be obtained for the study period, data from 2019 were used
instead.

Normalized difference vegetation index (NDVI) is a well-known remote sensing
indicator of green vegetation areas in satellite images Miura et al. (2019).
It detects vegetation as the difference between near-infrared and red light,
in the value range from -1 to +1. In general, high NDVI scores include urban
green spaces such as official parks, backyards, street trees, mountains,
riverbanks, golf courses, and urban farmlands. There are a few well-known
variants of NDVI Markevych et al. (2017), such as the soil-adjusted vegetation
index (SAVI) Huete (1988), which is corrected for soil brightness, and the
enhanced vegetation index (EVI) Jiang et al. (2008), which is corrected for
atmospheric effects. All NDVI, SAVI, and EVI2 scores can be calculated from
the two spectral bands of Sentinel-2, red (band 4) and near-infrared (NIR,
band 8), as follows:

$NDVI=\frac{NIR-RED}{NIR+RED},$ (4) $SAVI=\frac{(1+L)(NIR-RED)}{NIR+RED+L},$
(5) $EVI2=2.5\frac{NIR-RED}{NIR+2.4RED+1}.$ (6)

The robustness of the results for the three green space measures was verified
using NDVI as the primary metric.

###  VI.2 Measuring the amount of green space

The vegetation indices are measured in three steps, as illustrated in Fig.
4(a). The first step is to identify target cities containing at least 10% of
each country’s total population and represent the country’s overall happiness.
The second step is to extract only the built-up areas within the identified
cities’ administrative boundaries. As cities’ boundaries are historically and
culturally constructed and often arbitrary, the cities’ size needs to be
standardized; some cities include vast suburban areas (e.g., Istanbul) or
natural areas (e.g., deserts in Dubai). Thus, referring to the global land
cover data from the EU’s Copernicus Programme Buchhorn et al. (2020), we focus
on urban built-up areas to quantify the urban green space. Finally, the
vegetation indices (NDVI, EVI2, and SAVI) are calculated for the extracted
urban areas.

The final step is to compute the amount of green space in each country,
determined from the measured vegetation indices. Here, we define the amount of
green space as the logarithm of the total NDVI of built-up areas in the target
cities divided by the cities’ total population, called UGS, as a metric for
urban green space. UGS is calculated as follows:

$UGS=\log{\left(\frac{\sum_{c}\sum_{b(c)}NDVI(b)}{\sum_{c}N_{c}}\right)},$ (7)

where $NDVI(b)$ is the value of NDVI of pixel $b$ within built-up areas $b(c)$
in city $c$ and $N_{c}$ is the population of city $c$. In this calculation, we
adjusted negative NDVI values to zero Markevych et al. (2017) to prevent
errors caused by the accumulation of negative values in areas next to bodies
of water (see the Supplementary Information for the entire dataset).

## Acknowledgements

The authors thank to Farnoosh Hashemi, Ali Behrouz, and Taekho You for useful
comments. M.Cha work was supported by the Institute for Basic Science
(IBS-R029-C2).

## Author contributions statement

M.C. and D.Y.W. conceived the research, I.H., W.-S.J. and M.C. designed the
research, O.-H.K. and J.Y. collected the data, O.-H.K. performed the research,
O.-H.K., I.H. and M.C. analysed the data, O.-H.K., I.H. and J.Y. wrote the
manuscript. All authors reviewed the manuscript.

## References

  * de Vries et al. (2003) S. de Vries, R. A. Verheij, P. P. Groenewegen, and P. Spreeuwenberg, Environment and Planning A: Economy and Space 35, 1717 (2003).
  * Dadvand et al. (2016) P. Dadvand, X. Bartoll, X. Basagaña, A. Dalmau-Bueno, D. Martinez, A. Ambros, M. Cirach, M. Triguero-Mas, M. Gascon, C. Borrell, et al., Environmental International 91, 161 (2016).
  * UN (2015) UN, _Sustainable Development Goals_ (2015), available at https://sdgs.un.org/goals. Date accessed 4 November 2020.
  * de Vries et al. (2013) S. de Vries, S. M. E. van Dillen, P. P. Groenewegen, and P. Spreeuwenberg, Social Science & Medicine 94, 26 (2013).
  * Liu et al. (2019) Y. Liu, R. Wang, G. Grekousis, Y. Liu, Y. Yuan, and Z. Li, Landsacpe and Urban Planning 190 (2019).
  * Dzhambov et al. (2018) A. Dzhambov, T. Hartig, I. Markevych, B. Tilov, and D. Dimitrova, Environmental Research 160, 47 (2018).
  * Maas et al. (2009) J. Maas, S. M. E. van Dillen, R. A. Verheij, and P. P. Groenewegen, Health & Place 15, 586 (2009).
  * Jennings and Bamkole (2019) V. Jennings and O. Bamkole, International journal of environmental research and public health 16, 452 (2019).
  * Kweon et al. (1998) B. S. Kweon, W. C. Sullivan, and A. R. Wiley, Environment and Behavior 30, 832 (1998).
  * Maas et al. (2006) J. Maas, R. A. Verheij, P. P. Groenewegen, S. de Vries, and P. Spreeuwenberg, Journal of Epidemiology & Community Health 60, 587 (2006).
  * MacKerron and Mourato (2013) G. MacKerron and S. Mourato, Global Environmental change 23, 992 (2013).
  * Alcock et al. (2014) I. Alcock, M. P. White, B. W. Wheeler, L. E. Fleming, and M. H. Depledge, Environmental Science & Technology 48, 1247 (2014).
  * Beyer et al. (2014) K. M. M. Beyer, A. Kaltenbach, A. Szabo, S. Bogar, F. J. Nieto, and K. M. Malecki, International Journal of Environmental Research and Public Health 11, 3453 (2014).
  * Tsai et al. (2018) W.-L. Tsai, M. R. Mchale, V. Jennings, O. Marquet, J. A. Hipp, Y.-F. Leung, and M. F. Floyd, International Journal of Environmental Research and Public Health 15 (2018).
  * Stigsdotter et al. (2010) U. K. Stigsdotter, O. Ekholm, J. Schipperijn, M. Toftager, F. Kamper-Jørgensen, and T. B. Randrup, Scandinavian Journal of Public Health 38, 411 (2010).
  * EOS (2019) EOS, _NDVI FAQ: ALL YOU NEED TO KNOW ABOUT NDVI_ (2019), available at https://eos.com/blog/ndvi-faq-all-you-need-to-know-about-ndvi/. Date accessed 22 June 2020.
  * van den Berg et al. (2016) M. van den Berg, M. van Poppel, I. van Kamp, S. Andrusaityte, B. Balseviciene, M. Cirach, A. Danileviciute, N. Ellis, G. Hurst, D. Masterson, et al., Health & Place 38, 8 (2016).
  * Markevych et al. (2017) I. Markevych, J. Schoierer, T. Hartig, A. Chudnovsky, P. Hystad, A. M. Dzhambov, S. de Vries, M. Triguero-Mas, M. Brauer, M. J. Nieuwenhuijsen, et al., Environmental Research 158, 301 (2017).
  * Helliwell et al. (2018) J. F. Helliwell, R. Layard, and J. D. Sachs, New York: UN Sustainable Development Solutions Network (2018).
  * Kahneman and Deaton (2010) D. Kahneman and A. Deaton, Proceedings of the National Academy of Sciences 107, 16489 (2010).
  * Miura et al. (2019) T. Miura, S. Nagai, M. Takeuchi, K. Ichii, and H. Yoshioka, Scientific Reports 9 (2019).
  * IMF (2018) IMF, _World Economic Outlook Database_ (2018), available at https://www.imf.org/en/Publications/SPROLLS/world-economic-outlook-databases. Date accessed 16 November 2020.
  * Easterlin et al. (2010) R. A. Easterlin, L. A. McVey, M. Switek, O. Sawangfa, and J. S. Zweig, Proceedings of the National Academy of Sciences 107, 22463 (2010).
  * Preacher et al. (2007) K. J. Preacher, D. D. Rucker, and A. F. Hayes, Multivariate Behavioral Research 42, 185–227 (2007).
  * Foley and Kistemann (2015) R. Foley and T. Kistemann, Health & Place 35, 157–165 (2015).
  * Raymond et al. (2016) C. M. Raymond, S. Gottwald, and M. Kytta, Landscape and Urban Planning 153, 198–208 (2016).
  * Groff and McCord (2012) E. Groff and E. S. McCord, Security journal 25, 1 (2012).
  * Han et al. (2018) B. Han, D. A. Cohen, K. P. Derose, J. Li, and S. Williamson, American journal of preventive medicine 54, 352 (2018).
  * Ugolini et al. (2020) F. Ugolini, L. Massetti, P. Calaza-Martínez, P. Cariñanos, C. Dobbs, S. K. Ostoić, A. M. Marin, D. Pearlmutter, H. Saaroni, I. Šaulienė, et al., Urban Forestry & Urban Greening 56, 126888 (2020).
  * Lolli et al. (2020) S. Lolli, Y.-C. Chen, S.-H. Wang, and G. Vivone, Scientific Reports 10, 16213 (2020).
  * Hedblom et al. (2019) M. Hedblom, B. Gunnarsson, B. Iravani, I. Knez, M. Schaefer, P. Thorsson, and J. N. Lundström, Scientific Reports 9, 10113 (2019).
  * Ewing (2008) R. H. Ewing, _Characteristics, Causes, and Effects of Sprawl: A Literature Review_ (Springer US, Boston, MA, 2008).
  * Liu et al. (2020) X. Liu, Y. Huang, X. Xu, X. Li, X. Li, P. Ciais, P. Lin, K. Gong, A. D. Ziegler, A. Chen, et al., Nature Sustainability 3, 564 (2020).
  * Allen et al. (2010) C. D. Allen, A. K. Macalady, H. Chenchouni, D. Bachelet, N. McDowell, M. Vennetier, T. Kitzberger, A. Rigling, D. D. Breshears, E. H. Hogg, et al., Forest Ecology and Management 259, 660 (2010).
  * Pretzsch et al. (2017) H. Pretzsch, P. Biber, E. Uhl, J. Dahlhausen, G. Schütze, D. Perkins, T. Rötzer, J. Caldentey, T. Koike, T. van Con, et al., Scientific Reports 7 (2017).
  * Nowak et al. (2013) D. J. Nowak, E. J. Greenfield, R. E. Hoehn, and E. Lapoint, Environmental Pollution 178, 229 (2013).
  * Houlden et al. (2017) V. Houlden, S. Weich, and S. Jarvis, BMC public health 17, 460 (2017).
  * Seiferling et al. (2017) I. Seiferling, N. Naik, C. Ratti, and R. Proulx, Landscape and Urban Planning 165, 93 (2017).
  * Huete (1988) A. R. Huete, Remote Sensing of Environment 25, 295 (1988).
  * Jiang et al. (2008) Z. Jiang, A. R. Huete, K. Didan, and T. Miura, Remote Sensing of Environment 112, 3833 (2008).
  * Buchhorn et al. (2020) M. Buchhorn, B. Smets, L. Bertels, B. De Roo, M. Lesiv, N.-E. Tsendbazar, M. Herold, and S. Fritz, _Copernicus Global Land Service: Land Cover 100m: collection 3 epoch 2015, Globe_ (2020), available at https://lcviewer.vito.be/download.