# Categorizing models using Self-Organizing Maps: an application to modified
gravity theories probed by cosmic shear

Agnès Ferté1 Shoubaneh Hemmati2 Daniel Masters2 Brigitte Montminy1 Peter L.
Taylor1 Eric Huff1 Jason Rhodes1 1 Jet Propulsion Laboratory, California
Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, 91109, USA
2IPAC, California Institute of Technology, 1200 E California Blvd, Pasadena,
CA 91125, USA

###### Abstract

We propose to use Self-Organizing Maps (SOM) to map the impact of physical
models onto observables. Using this approach, we are be able to determine how
theories relate to each other given their signatures. In cosmology this will
be particularly useful to determine cosmological models (such as dark energy,
modified gravity or inflationary models) that should be tested by the new
generation of experiments. As a first example, we apply this approach to the
representation of a subset of the space of modified gravity theories probed by
cosmic shear. We therefore train a SOM on shear correlation functions in the
$f(R)$, dilaton and symmetron models. The results indicate these three
theories have similar signatures on shear for small values of their parameters
but the dilaton has different signature for higher values. We also show that
modified gravity (especially the dilaton model) has a different impact on
cosmic shear compared to a dynamical dark energy so both need to be tested by
galaxy surveys.

††©2023. All rights reserved.

##  1 Introduction

In many fields a large variety of theoretical models are proposed in the
literature and are to be tested with current or future data. For example, in
cosmology 74 inflation models were systematically tested in [Martin  _et
al._(2013)Martin, Ringeval, and Vennin] using cosmic microwave background 2013
data from the Planck satellite [Ade  _et al._(2014)Ade, Aghanim, Alves,
Armitage-Caplan, Arnaud, Ashdown, Atrio-Barandela, Aumont, Aussel, and et
al.], or in astrophysics, testing various models of exoplanets’ atmosphere
will be one of the goals of experiments like the JWST as reviewed by
[Madhusudhan(2019), Biller and Bonnefoy(2018)]. Selecting only a subset of
models to be tested given an observable is in some cases necessary as
obtaining robust constraints can be time consuming, needs modeling
developments specific to each theory or requires large computing power (to run
Monte Carlo Markov Chains for instance). However making a model selection
based on the inspection of observables can be inconvenient in the case where
models are characterized by several parameters or observables are described by
several components, such as tomographic real-space correlation functions in
cosmology. In this paper, we address this challenge by presenting a new
approach using unsupervised non-Gaussian dimensionality reduction algorithms
to facilitate model selection. We especially use Self-Organizing Maps (SOMs)
as they can produce 2-dimensional maps where models are assigned to cells
preserving topology. Such mapping indeed allows to easily determine models
that have different impacts on a given observable and the effect of
measurements on this mapping can also be added (but is not done in the present
study as described later on). This approach complements Principal Component
Analysis (PCA) as SOMs are adequate for non-Gaussian reduction and
vizualization. As a first step and for purpose of example, we apply our
approach to the case of modified gravity (MG) theories probed by cosmic shear.
A similar approach has been proposed in [Mancarella  _et al._(2022)Mancarella,
Kennedy, Bose, and Lombriser], using Bayesian Neural Networks as a classifier
trained on matter power spectra. In the present study, we instead use a
dimensionality reduction algorithm trained on cosmic shear, a probe of the
matter power spectrum.

Weak gravitational lensing is indeed a powerful probe of gravity as shown
early on in [Song(2005), Ishak  _et al._(2006)Ishak, Upadhye, and Spergel]
demonstrating that some MG theories can be distinguished from a dark energy
using cosmic shear, and later in [Song(2006), Tsujikawa and Tatekawa(2008),
Jain and Zhang(2008)]. Tomographic cosmic shear corresponds to the correlation
in redshift bins of the deformation of observed galaxy shapes due to weak
lensing. It directly probes the matter distribution and is sensitive to the
growth of large scale structures in the Universe, which depends on the laws of
gravity. Deviations from GR were thus tested using cosmic shear data from
CFHTLens111https://www.cfhtlens.org/ [Harnois-Déraps  _et al._(2015)Harnois-
Déraps, Munshi, Valageas, van Waerbeke, Brax, Coles, and Rizzo, Simpson  _et
al._(2013)Simpson, Heymans, Parkinson, Blake, Kilbinger, Benjamin, Erben,
Hildebrandt, Hoekstra, Kitching _et al._, Dossett  _et al._(2015)Dossett,
Ishak, Parkinson, and Davis, Ferté  _et al._(2019)Ferté, Kirk, Liddle, and
Zuntz], KIDS222http://kids.strw.leidenuniv.nl/ [Joudaki  _et
al._(2017)Joudaki, Mead, Blake, Choi, de Jong, Erben, Conti, Herbonnet,
Heymans, Hildebrandt _et al._, Tröster  _et al._(2021)Tröster, Asgari, Blake,
Cataneo, Heymans, Hildebrandt, Joachimi, Lin, Sá nchez, Wright _et al._] and
DES333https://www.darkenergysurvey.org/ [Abbott  _et al._(2019)Abbott,
Abdalla, Avila, Banerji, Baxter, Bechtol, Becker, Bertin, Blazek, Bridle _et
al._] (in the later cases galaxy-galaxy lensing and clustering data were also
used). Forecasts show these constraints will improve with future imaging
surveys including the Euclid satellite444https://sci.esa.int/web/euclid/
[Martinelli  _et al._(2011)Martinelli, Calabrese, De Bernardis, Melchiorri,
Pagano, and Scaramella, Euclid Collaboration(2020)],
LSST555https://www.lsst.org/ that will be produced by the Rubin Observatory
[The LSST Dark Energy Science Collaboration(2018), Hojjati  _et
al._(2016)Hojjati, Plahn, Zucca, Pogosian, Brax, Davis, and Zhao, Ferté  _et
al._(2019)Ferté, Kirk, Liddle, and Zuntz] and the Roman space
mission666https://roman.gsfc.nasa.gov/ [Doré  _et al._(2019)Doré, Hirata,
Wang, Weinberg, Eifler, Foley, Heinrich, Krause, Perlmutter, Pisani _et al._,
Eifler  _et al._(2021a)Eifler, Miyatake, Krause, Heinrich, Miranda, Hirata,
Xu, Hemmati, Simet, Capak, Choi, Doré , Doux, Fang, Hounsell, Huff, Huang,
Jarvis, Kruk, Masters, Rozo, Scolnic, Spergel, Troxel, von der Linden, Wang,
Weinberg, Wenzl, and Wu, Eifler  _et al._(2021b)Eifler, Simet, Krause, Hirata,
Huang, Fang, Miranda, Mandelbaum, Doux, Heinrich _et al._].

One of the challenges faced by these future surveys is to decide which
theories of MG to test amongst the large space of models proposed in the
literature. For instance [Ishak  _et al._(2019)Ishak, Baker, Bull, Pedersen,
Blazek, Ferreira, Leonard, Lin, Linder, Pardo, and Valogiannis] prioritizes
theories to be tested by LSST according to their theoretical maturity and
impact on LSST observables, however without comparing their signatures on
observables. Moreover, for surveys such as LSST, cosmic shear is expected to
be measured through correlation functions in 10 redshift bins, corresponding
to 55 cross- and auto-correlations, making a visual inspection of the impact
of MG models on cosmic shear difficult. The goal of the present approach is to
fill the gap between choosing models to test based on their theoretical
interest (such as in [Ishak  _et al._(2019)Ishak, Baker, Bull, Pedersen,
Blazek, Ferreira, Leonard, Lin, Linder, Pardo, and Valogiannis]) and testing a
wide variety of theories (e.g. [Planck Collaboration(2016)]) or using non-
parametric reconstruction [Raveri  _et al._(2021)Raveri, Pogosian, Koyama,
Martinelli, Silvestri, Zhao, Li, Peirone, and Zucca] or PCA [Linder and
Huterer(2005), Asaba  _et al._(2013)Asaba, Hikage, Koyama, Zhao, Hojjati, and
Pogosian]): we propose a way to categorize models in two dimensions according
to their impact on observables using SOMs. We do so by producing a training
set made of theoretical predictions of cosmic shear in MG models assuming a
fixed cosmology. We indeed want to know how MG models would compare to each
other if we had a perfect knowledge of cosmology. If two MG models have the
same impact at fixed cosmology then there is definitely no chance to
distinguish between them when doing a standard cosmological inference
analysis, studying one model is therefore sufficient. So in this paper, we set
$\Omega_{m}$ = 0.31, $H_{0}$ = 67, $\Omega_{b}$ = 0.04, $A_{s}=2.1\times
10^{-9}$, $n_{s}=0.96$, one massive neutrino with a mass of 0.06 eV. The
training sets are available on request and the script to reproduce results
presented in this paper is available
online777https://github.com/aferte/CosmoSOM-MG.

##  2 Self-Organizing Map of modified gravity probed by cosmic shear

Self-Organizing Maps – Self-Organizing Maps, introduced in [Kohonen(1982)],
are a class of unsupervised neural networks which learn and reduce dimensions
(usually to 2 dimensions) of a higher dimensional data set while preserving
the topology of the manifold. Dimensionality reduction makes SOMs perfect
tools for visualization. The training of a SOM is an iterative process where
weights of the neurons within the radius of the neighborhood function are
updated according to their distance to the input data, where the neuron
closest to the input data is called the Best Matching Unit (BMU). At the end
of the training, the neurons are weighted so that they match the distribution
of the input data. In astronomy, SOMs have been used from classifications of
stars/galaxies [Maehoenen and Hakala(1995)] to measurements of the physical
properties of galaxies as in [Hemmati  _et al._(2019)Hemmati, Capak,
Pourrahmani, Nayyeri, Stern, Mobasher, Darvish, Davidzon, Ilbert, Masters, and
Shahidi], to photometric redshift measurement improvements for weak lensing
cosmology as in [Masters  _et al._(2015)Masters, Capak, Stern, Ilbert,
Salvato, Schmidt, Longo, Rhodes, Paltani, Mobasher _et al._, Wright  _et
al._(2020)Wright, Hildebrandt, van den Busch, and Heymans, Myles  _et
al._(2021)Myles, Alarcon, Amon, Sánchez, Everett, DeRose, McCullough, Gruen,
Bernstein, Troxel _et al._].

In this work, we use SOMPY888https://github.com/sevamoo/SOMPY ([Moosavi  _et
al._(2014)Moosavi, Packmann, and Vallés] \- a python library for SOM) with a
rectangular grid. Each cell or neuron on the grid has a weight vector with the
size of the input data set dimensions. We assign the initial weights through
PCA and tune them in each epoch of the training phase by comparing them with
the training data using an Euclidean metric. Our choice of distance metrics is
currently limited by metrics available in the SOMPY package. Alternative
metrics taking into account the covariance on data vector measurements will be
explored in future work. The training is considered done when the weights of
the cells represent the training data space well. In our case, we use the
default limit radius of the neighborhood function in SOMPY as a stopping point
for the training.

The choice of grid size is not well-defined and depends mostly on the size and
variations of the training sample. A grid axis ratio of 1:1 has been shown in
[Davidzon  _et al._(2019)Davidzon, Laigle, Capak, Ilbert, Masters, Hemmati,
Apostolakos, Coupon, de la Torre, Devriendt, Dubois, Kashino, Paltani, and
Pichon] to slightly outperform the non-square choices, we hence choose a
$6\times 6$ cells grid here, which satisfies the needs of this study (we
tested that a larger grid led to a large fraction of empty cells and a smaller
grid didn’t allow to show the differences between theories).

Training set – We use theoretical predictions of the cosmic shear real-space
correlation functions $\xi_{\pm}^{ij}(\theta)$ at angular separation $\theta$
in redshift bins $i$ and $j$. The SOM is trained on theoretical predictions of
$\xi_{\pm}^{ij}(\theta)$ in various modified gravity (MG) theories, for
different values of their parameters. We model the correlation functions
$\xi_{\pm}^{ij}(\theta)$ as:

$\xi_{\pm}^{ij}(\theta)=\int\frac{d\ell\ell}{2\pi}J_{0/4}(\ell\theta)C_{\ell}^{\kappa_{i}\kappa_{j}},$
(1)

where $J_{0/4}$ are the Bessel functions of order 0 and 4, and
$C_{\ell}^{\kappa_{i}\kappa_{j}}$ is the convergence $\kappa$ power spectrum
between redshift bins $i$ and $j$:

$C_{\ell}^{\kappa_{i}\kappa_{j}}=\int
d\chi\frac{q_{\kappa}^{i}q_{\kappa}^{j}}{\chi^{2}}P(k=\frac{\ell+0.5}{\chi},z(\chi)),$
(2)

where $\chi$ is comoving distance, $P(k,z)$ the matter power spectrum and
$q_{\kappa}^{i}$ is the radial weight function:
$q_{\kappa}^{i}=\frac{3H_{0}^{2}\Omega_{m}\chi}{2a(\chi)}\int
d\chi^{\prime}\frac{\chi^{\prime}-\chi}{\chi}n^{i}(z(\chi^{\prime}))\frac{dz}{d\chi^{\prime}}$.

Figure 1: Redshift distribution of galaxies $n(z)$ used to compute the cosmic
shear correlation functions $\xi_{\pm}(\theta)$ forming the training set.

To determine the redshift distribution $n(z)$ entering the radial weight
function, we consider galaxies distributed in 10 redshift bins, similar to
what is expected for LSST Year 10 [The LSST Dark Energy Science
Collaboration(2018)]. The forecasted binned redshift distribution is shown on
figure 1, spanning from redshift 0 to 4, following
$z^{\alpha}\exp{-(z/z_{0})^{\beta}}$ with $\alpha=2$, $\beta=1$ and
$z_{0}=0.3$, convolved by a Gaussian photometric error of variance $\sigma$ =
0.02 at $z=0$. We use the cosmological code
CosmoSIS999https://bitbucket.org/joezuntz/cosmosis/wiki/Home [Zuntz  _et
al._(2015)Zuntz, Paterno, Jennings, Rudd, Manzotti, Dodelson, Bridle, Sehrish,
and Kowalkowski] to compute the theoretical predictions of cosmic shear, where
the impact of modified gravity on the matter power spectrum $P(k,z)$ is
computed using the latest version of MGCAMB101010https://github.com/sfu-
cosmo/MGCAMB [Zucca  _et al._(2019)Zucca, Pogosian, Silvestri, and Zhao,
Hojjati  _et al._(2011)Hojjati, Pogosian, and Zhao, Zhao  _et al._(2009)Zhao,
Pogosian, Silvestri, and Zylberberg], implemented in CosmoSIS. For the purpose
of this first study we assume a linear matter power spectrum. Indeed, while
various approaches have been proposed to model non-linearities in modified
gravity theories, an alternative is to remove measurements of cosmic shear on
non-linear scales for the cosmological analysis (as done in [Abbott  _et
al._(2019)Abbott, Abdalla, Avila, Banerji, Baxter, Bechtol, Becker, Bertin,
Blazek, Bridle _et al._] for instance). In a similar way, we use the cosmic
shear correlation functions in a range of $\theta$ roughly similar to the
conservative scale cuts used for the constraints on deviations to GR in
[Abbott  _et al._(2019)Abbott, Abdalla, Avila, Banerji, Baxter, Bechtol,
Becker, Bertin, Blazek, Bridle _et al._]111111Such scale cuts depend on the
analysis’ specifications (such as the number of redshift bins or binning in
$\theta$). We choose the scale cuts from [Abbott  _et al._(2019)Abbott,
Abdalla, Avila, Banerji, Baxter, Bechtol, Becker, Bertin, Blazek, Bridle _et
al._] for simplicity, the impact of non-linearities will be described in
future work.. We therefore use predictions of $\xi_{+}$ between 30 and 300
arcminutes and of $\xi_{-}$ between 250 and 300 arcminutes for all redshift
combinations. All the cosmic shear correlations are then turned into a one-
dimensional vector of 45,375 elements (corresponding to 765 values of $\theta$
for $\xi+$, 60 for $\xi_{-}$ in 55 redshift bin correlations.)

We consider modified gravity theories that are readily available in MGCAMB and
have an impact on cosmic probes, following [Hojjati  _et al._(2016)Hojjati,
Plahn, Zucca, Pogosian, Brax, Davis, and Zhao], and compute $\xi_{\pm}$ for
100 different values of the corresponding MG parameters. These models
reproduce a cosmological constant in the background evolution of the recent
Universe. These values are linearly spaced within a range that we choose to
not be informed by astrophysical constraints as cosmological analyses often
use such large priors. We use:
\- The Hu-Sawicki model of $f(R)$ gravity [Hu and Sawicki(2007)] which can
describe a cosmic acceleration while satisfying solar system tests of gravity.
In this model, the GR action is modified by a function of the Ricci scalar
$R$. This model is described by the Compton wavelength parameter $B_{0}$ [Song
_et al._(2007)Song, Hu, and Sawicki, Dossett  _et al._(2014)Dossett, Hu, and
Parkinson, Planck Collaboration(2016)]: we vary $\log_{10}{B_{0}}$ linearly
between -10 and -2. We focus on this parametrisation as one of the tightest
cosmological constraints on $f(\mathrm{R})$ is made on $B_{0}$ in [Planck
Collaboration(2016)]. An alternative $f(\mathrm{R})$ parametrisation uses
$f_{R_{0}}$, which corresponds to the value of $\mathrm{d}f(R)/\mathrm{d}R$
today and $n$, the index of the power law $f(R)$ follows. We produced a
training set for this parametrisation which is available on request, although
we focus on the $B_{0}$ parametrisation here.
\- The dilaton model [Damour and Polyakov(1994), Brax  _et al._(2010)Brax, van
de Bruck, Davis, and Shaw] parametrized by the value of the force range and
the coupling to matter today, respectively $\xi_{0}$ and $\beta_{0}$. We vary
$\log_{10}{\xi_{0}}$ linearly between -6 and -3.3 (above this value our
current pipeline gives a negative shear correlated functions) for $\beta_{0}$
set to $\Omega_{\Lambda,0}/\Omega_{m,0}$ = 2.2 as indicated in [Hojjati  _et
al._(2016)Hojjati, Plahn, Zucca, Pogosian, Brax, Davis, and Zhao] (for the
scalar field to explain dark energy).
\- The symmetron model [Hinterbichler and Khoury(2010)] parametrized by the
force range $\xi_{\star}$, the scale factor at the time of the force
activation $a_{\star}$ and the coupling to matter $\beta_{\star}$. We vary
$\log_{10}(\xi_{\star})$ between -6 and -2 for $\beta_{\star}$ set to 0.5, 1
and 1.5, while we fix $a_{\star}$ for simplicity, as is done in [Hojjati  _et
al._(2016)Hojjati, Plahn, Zucca, Pogosian, Brax, Davis, and Zhao] for
instance. We choose to set $a_{\star}$ to 0.5 as it is the standard choice,
but note that cosmic shear depends on the value of $a_{\star}$, a different
choice would thus lead to a different SOM.

Model | Parameters | Prior range | Current cosmological constraints
---|---|---|---
f(R) | $B_{0}$ | [$10^{-10}$,$10^{-2}$] | $B_{0}<0.86\times 10^{-4}$
Dilaton | $\xi_{0},\beta_{0}$ | [$10^{-6}$,$5\times 10^{-4}$], 2.2 |  $\xi_{0}<$ 3$\times 10^{-3}$
Symmetron | $\xi_{*},\beta_{*},a_{*}$ | [$10^{-6}$,$10^{-2}$], [0.5,1.5], 0.5 |  $\xi_{*}<$ 1.8$\times 10^{-3}$
Table 1: Parameters varied, their value range and current upper bounds (at
95$\%$ confidence level) for the modified gravity models used to compute the
SOM training data set (in the dilaton case, $\beta_{0}$ is fixed to
$\Omega_{\Lambda,0}/\Omega_{m,0}$).

The parameters’ range used to produce our training set mostly contain values
of modified gravity parameters still allowed by cosmological measurements.
Firstly, [Planck Collaboration(2016)] placed an upper bound on the $B_{0}$
parameter of $f(R)$ gravity. In this case the 95% confidence level upper limit
is $B_{0}<0.86\times 10^{-4}$ when combining temperature and polarisation CMB
measurements with cosmic shear, baryonic acoustic oscillation (BAO) and
redshift space distortion measurements. We note that cosmic shear measurements
were used to constrain $f_{R_{0}}$and $n$ in [Vazsonyi  _et
al._(2021)Vazsonyi, Taylor, Valogiannis, Ramachandra, Ferté, and Rhodes] using
HSC data and [Tröster  _et al._(2021)Tröster, Asgari, Blake, Cataneo, Heymans,
Hildebrandt, Joachimi, Lin, Sá nchez, Wright _et al._] using KiDS-1000 data,
and in both cosmological analyses, the results were prior dominated. Cosmic
shear data were also used in [Harnois-Déraps  _et al._(2015)Harnois-Déraps,
Munshi, Valageas, van Waerbeke, Brax, Coles, and Rizzo] leading to rejection
of $f(R)$ gravity with $f_{R_{0}}$$=10^{-4}$, however at fixed $n$.
Furthermore, using screening properties, tests of gravity on astrophysical
scales have also been developed, for instance in [Jain  _et al._(2013)Jain,
Vikram, and Sakstein] using local distance indicators and [Desmond and
Ferreira(2020)] using galaxy morphology. These analyses obtained much more
stringent upper bounds on $f_{R_{0}}$, ruling out models with $f_{R_{0}}$above
$\sim 10^{-7}$. Cosmological tests are however still relevant, at least to
reproduce such constraints with future large scale structures data. Secondly,
[Hojjati  _et al._(2016)Hojjati, Plahn, Zucca, Pogosian, Brax, Davis, and
Zhao] used CMB measurements in addition to BAO, cosmic shear and matter power
spectrum available at the time leading to 95% confidence limit upper bound of
$3\times 10^{-3}$ on $\xi_{0}$ for the dilaton and of $1.8\times 10^{-3}$ for
$\xi_{\star}$ in the symmetron model, when varying the neutrino mass. Table 1
summarizes the varied parameters, their corresponding ranges and current upper
bounds (at 95% confidence level) for each of the considered MG theories. We
therefore compute shear correlation functions for a total of 500 models.

##  3 Results and discussion

SOM of modified gravity theories – After training the $6\times 6$ square grid
SOM on the training set described above, we map the training set on the SOM.
Fig. 2 shows the resulting two dimensional SOM grid with cells colored by the
final number of models per neuron.

Figure 2: Self-Organizing Map trained on shear correlation functions
$\xi_{\pm}(\theta)$ in $f(R)$ gravity and the dilaton and symmetron models,
where cells are colored by the number of models per BMU of the SOM. Shear
correlation functions $\xi_{+}^{22}(\theta)$ corresponding to cells numbered
from 1 to 3 are represented on Fig. 3.

In this figure and the following ones, cells in white correspond to cells with
no models. The cell corresponding to the highest number of models lie in the
left top corner (numbered 1 in Fig. 2 with a total of 132 models in this
cell), which corresponds to General Relativity or small deviations to GR. This
indicates that shear correlation functions for small values of the modified
gravity parameters ($B_{0}$, $\xi_{0}$ and $\xi_{\star}$) are
indistinguishable from GR.

Figure 3: Shear auto-correlations $\xi^{22}_{+}(\theta)$ in the second
redshift bin associated to the BMUs numbered from 1 to 3 on Fig. 2,
respectively in black, blue and purple. The dashed and dotted lines correspond
to different models in the same BMU corresponding to the closest (in dashed
lines) and most distance (in dotted lines) model to the cosmic shear
prediction in GR within one cell. For purpose of illustration, the black error
bars show the expected error bars on $\xi_{+}^{22}$ for a survey like LSST.

Fig. 3 shows correlation functions $\xi^{22}_{+}(\theta)$ in redshift bin 2
for different cells (which are numbered from 1 to 3 in Fig. 2). The dotted and
dashed lines for cells 2 and 3 respectively show the closest and furthest
(according to the Euclidean distance metric) models associated to these two
cells to GR. The black line respresents the median of models in cell 1, these
models present negligible variations as the MG parameters are very small (e.g.
the median of $B_{0}$ in cell 1 is 3.1$\times 10^{-9}$). We also show for
purpose of illustration the expected error bars on $\xi_{+}^{22}(\theta)$ from
a survey like LSST (10 source redshift bin with 55 galaxies
$\mathrm{arcmin}^{-2}$ and a total shape noise of 0.3 over 50$\%$ of the sky),
which we computed using CosmoCov [Krause and Eifler(2017), Fang  _et
al._(2020a)Fang, Krause, Eifler, and MacCrann, Fang  _et al._(2020b)Fang,
Eifler, and Krause].

Other cells correspond to the various modified gravity theories used in the
training, so we show in Fig. 4 a different representation of this SOM, now
highlighting only the cells corresponding to each MG theory. In this case, the
cells are colored according to the median of the corresponding MG parameter in
each cell. The panels correspond to $f(R)$ gravity, the dilaton and symmetron
models from left to right, only showing the distribution of $\xi_{\star}$ for
$\beta_{\star}$ = 1 in the symmetron model.

Figure 4: SOM of modified gravity from Fig. 2 showing the median of the values
of the parameter for each cell for $f(R)$ gravity, the dilaton (for
$\beta_{0}$ = 2.2) and symmetron (for $\beta_{\star}$ = 1) models from left to
right.

The 3 theories have a large overlap in the top left part of the SOM for small
values of the MG parameters. However, for larger values, the dilaton model
correspond to cells different from those corresponding to $f(R)$ and the
symmetron model. This means the shear vectors for the dilaton were associated
to different BMUs than other models, which indicates that its signature is
different enough that it should be constrained by future analysis.

---
Figure 5: Number of models in each cell for ($w_{0}$,$w_{a}$) on the top left
panel and dilaton theories on the top right panel of the 6$\times$6 Self-
Organizing Map grid trained on shear correlation functions $\xi_{\pm}(\theta)$
computed in ($w_{0}$,$w_{a}$) model and $f(R)$ gravity, dilaton and symmetron
(for $\beta_{\star}=$ 0.5, 1 and 1.5) models. The bottom panel highlights the
$f(R)$ gravity model on the left and the symmetron model on the right panel.

Application: does a dynamical dark energy have the same signature as MG
theories – Current and future galaxy surveys aim at testing if the equation of
state $w$ of dark energy varies with time. The equation of state in this case
often follows [Linder(2003), Chevallier and Polarski(2001)] i.e. is
parametrized by $w_{0}$ (the value of $w$ today) and $w_{a}$ (its dependence
with time) as: $w(a)=w_{0}+w_{a}(1-a(t))$ where $a(t)$ is the scale factor. We
want to know if constraining a dynamical dark energy using cosmic shear is
sufficient to explain signatures from MG theories, using the present SOM
approach. To this end, we compute the shear correlation functions for
different combinations of $w_{0}$ and $w_{a}$ using the same pipeline as used
in modified gravity replacing MGCamb by CAMB [Lewis  _et al._(2000)Lewis,
Challinor, and Lasenby, Howlett  _et al._(2012)Howlett, Lewis, Hall, and
Challinor, Lewis(2014)]. We thus produce a set of shear correlation functions
corresponding to 96 combinations121212the pipeline failed for some
combinations of $w_{0}+w_{a}>$ 0 because of reionization properties leading to
convergence failure so we get 96 combinations instead of 100. of
($w_{0}$,$w_{a}$) randomly sampled from the chosen priors: [-1.5,-0.5] for
$w_{0}$, [-1,1] for $w_{a}$. To do so we use the apriori sampler that is
readily available in CosmoSIS. We then train a SOM on a combination of this
data set and the shear in modified gravity data set described above. The
resulting SOM is shown in Fig. 5 with the cells corresponding to dynamical
dark energy models highlighted on the left panel. We mapped the dilaton models
onto the SOM in the right panel to show that a dynamical dark energy does not
overlap with the dilaton model, except for parameters values close to GR (the
middle cell). We also show cells corresponding to cosmic shear in $f(R)$
gravity on the lower left panel and the symmetron model on the lower right
panel. Although the dilaton shows the smallest overlap with ($w_{0}$,$w_{a}$),
$f(R)$ gravity and the symmetron model also cover cells that are not covered
by a dynamical dark energy. This indicates these MG theories would need to be
constrained in addition to a dynamical dark energy in order to extract more
information from cosmic shear.

##  4 Conclusion

We presented a new approach to categorize models using a dimensionality
reduction algorithm such as Self-Organizing Map, as a promising way to help
guide analyses in the case of a large theory space, like in cosmology. To
summarize, one can follow the following steps to apply this approach: produce
the signal expected in the different models considered, train a SOM with this
training data set and inspect the resulting SOM to determine models of choice.
For this last step, we recommend imaging the number of models per cell in the
final SOM. Indeed, a large number of models in a cell or a group of cells
indicates models that could not be distinguished, however cells with a low
number of models correspond to one or a few models having different impact on
the observable and therefore models that are worth exploring further. More
generally, one can determine clusters of models on the SOM visually or through
clustering algorithms such as $k$-means which is included in SOMPY, this will
be especially useful for the case of a larger number of models (such as the 74
inflationary models tested in [Martin  _et al._(2013)Martin, Ringeval, and
Vennin]). Finally, a new data set can also be mapped on a trained SOM to
determine which models can explain it.

In this paper, we applied this approach to the case of a subset of modified
gravity theories, categorizing these theories through their impact on cosmic
shear. We showed for instance that the signatures left on cosmic shear from
$f(R)$ gravity, the dilaton and symmetron models are similar for small
deviations from GR but the dilaton can be distinguished from the other models
for larger parameters values, while signatures from $f(R)$ gravity are similar
to the ones left by the symmetron model for $\beta_{\star}=1$. Moreover we
showed MG theories have a different impact on cosmic shear compared to a
dynamical dark energy, indicating that MG models should indeed be explored by
future surveys such as LSST or Euclid in order to further exploit cosmic shear
data. However, this first analysis has several caveats: we assume a linear
matter power spectrum, consider a fixed cosmology and we use theoretical
predictions of the shear correlation functions $\xi_{\pm}(\theta)$ in finely
binned $\theta$. We will go beyond some of these assumptions in future work
and develop several applications. As such, this approach can be applied to
other cases in comsology such as cosmic inflation or any case of a large
theory space in physics, in order to easily identify interesting directions of
exploration.

## Acknowledgement

The research was carried out at the Jet Propulsion Laboratory, California
Institute of Technology, under a contract with the National Aeronautics and
Space Administration (80NM0018D0004). In particular, this work was initiated
with a JPL Data Science Pilot grant under a program run by Daniel Crichton and
Richard Doyle who provided early support and advice. PLT acknowledges support
for this work from a NASA Postdoctoral Program Fellowship. BM was supported by
the JPL/Caltech Summer Undergraduate Research Fellowship Program. AF would
like to thank Benjamin Giblin for indicating an issue in the implementation of
MGCAMB in CosmoSIS to compute $f(R)$ gravity cosmic shear predictions and Joe
Zuntz for his help to fix this issue; Xiao Fang for help with CosmoCov;
members of the JPL Dark Sector group and of the cosmology group at Caltech for
their support and inputs; and early audiences of this work for useful
discussions.

## References

  * [Martin  _et al._(2013)Martin, Ringeval, and Vennin] J. Martin, C. Ringeval, and V. Vennin, Encyclopaedia inflationaris (2013), arXiv:1303.3787 [astro-ph.CO] .
  * [Ade  _et al._(2014)Ade, Aghanim, Alves, Armitage-Caplan, Arnaud, Ashdown, Atrio-Barandela, Aumont, Aussel, and et al.] P. A. R. Ade, N. Aghanim, M. I. R. Alves, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, H. Aussel, and et al., Astronomy & Astrophysics 571, A1 (2014).
  * [Madhusudhan(2019)] N. Madhusudhan, Annual Review of Astronomy and Astrophysics 57, 617 (2019).
  * [Biller and Bonnefoy(2018)] B. A. Biller and M. Bonnefoy, Handbook of Exoplanets , 2107–2135 (2018).
  * [Mancarella  _et al._(2022)Mancarella, Kennedy, Bose, and Lombriser] M. Mancarella, J. Kennedy, B. Bose, and L. Lombriser, Phys. Rev. D 105, 023531 (2022).
  * [Song(2005)] Y.-S. Song, Phys. Rev. D 71, 024026 (2005).
  * [Ishak  _et al._(2006)Ishak, Upadhye, and Spergel] M. Ishak, A. Upadhye, and D. N. Spergel, Phys. Rev. D 74, 043513 (2006).
  * [Song(2006)] Song, arXiv:astro-ph/0602598 (2006).
  * [Tsujikawa and Tatekawa(2008)] S. Tsujikawa and T. Tatekawa, Physics Letters B 665, 325 (2008), arXiv: 0804.4343.
  * [Jain and Zhang(2008)] B. Jain and P. Zhang, Physical Review D 78, 063503 (2008), arXiv: 0709.2375.
  * [Harnois-Déraps  _et al._(2015)Harnois-Déraps, Munshi, Valageas, van Waerbeke, Brax, Coles, and Rizzo] J. Harnois-Déraps, D. Munshi, P. Valageas, L. van Waerbeke, P. Brax, P. Coles, and L. Rizzo, Monthly Notices of the Royal Astronomical Society 454, 2722 (2015).
  * [Simpson  _et al._(2013)Simpson, Heymans, Parkinson, Blake, Kilbinger, Benjamin, Erben, Hildebrandt, Hoekstra, Kitching _et al._] F. Simpson, C. Heymans, D. Parkinson, C. Blake, M. Kilbinger, J. Benjamin, T. Erben, H. Hildebrandt, H. Hoekstra, T. D. Kitching, _et al._ , Monthly Notices of the Royal Astronomical Society 429, 2249 (2013).
  * [Dossett  _et al._(2015)Dossett, Ishak, Parkinson, and Davis] J. N. Dossett, M. Ishak, D. Parkinson, and T. M. Davis, Phys. Rev. D 92, 023003 (2015).
  * [Ferté  _et al._(2019)Ferté, Kirk, Liddle, and Zuntz] A. Ferté, D. Kirk, A. R. Liddle, and J. Zuntz, Phys. Rev. D 99, 083512 (2019).
  * [Joudaki  _et al._(2017)Joudaki, Mead, Blake, Choi, de Jong, Erben, Conti, Herbonnet, Heymans, Hildebrandt _et al._] S. Joudaki, A. Mead, C. Blake, A. Choi, J. de Jong, T. Erben, I. F. Conti, R. Herbonnet, C. Heymans, H. Hildebrandt, _et al._ , Monthly Notices of the Royal Astronomical Society 471, 1259 (2017).
  * [Tröster  _et al._(2021)Tröster, Asgari, Blake, Cataneo, Heymans, Hildebrandt, Joachimi, Lin, Sá nchez, Wright _et al._] T. Tröster, M. Asgari, C. Blake, M. Cataneo, C. Heymans, H. Hildebrandt, B. Joachimi, C.-A. Lin, A. G. Sá nchez, A. H. Wright, _et al._ , Astronomy & Astrophysics 649, A88 (2021).
  * [Abbott  _et al._(2019)Abbott, Abdalla, Avila, Banerji, Baxter, Bechtol, Becker, Bertin, Blazek, Bridle _et al._] T. M. C. Abbott, F. B. Abdalla, S. Avila, M. Banerji, E. Baxter, K. Bechtol, M. R. Becker, E. Bertin, J. Blazek, S. L. Bridle, _et al._ (DES Collaboration), Phys. Rev. D 99, 123505 (2019).
  * [Martinelli  _et al._(2011)Martinelli, Calabrese, De Bernardis, Melchiorri, Pagano, and Scaramella] M. Martinelli, E. Calabrese, F. De Bernardis, A. Melchiorri, L. Pagano, and R. Scaramella, Physical Review D 83, 023012 (2011).
  * [Euclid Collaboration(2020)] Euclid Collaboration, Astronomy and Astrophysics 642, A191 (2020).
  * [The LSST Dark Energy Science Collaboration(2018)] The LSST Dark Energy Science Collaboration, The lsst dark energy science collaboration (desc) science requirements document (2018), arXiv:1809.01669 [astro-ph.CO] .
  * [Hojjati  _et al._(2016)Hojjati, Plahn, Zucca, Pogosian, Brax, Davis, and Zhao] A. Hojjati, A. Plahn, A. Zucca, L. Pogosian, P. Brax, A.-C. Davis, and G.-B. Zhao, Phys. Rev. D 93, 043531 (2016).
  * [Doré  _et al._(2019)Doré, Hirata, Wang, Weinberg, Eifler, Foley, Heinrich, Krause, Perlmutter, Pisani _et al._] O. Doré, C. Hirata, Y. Wang, D. Weinberg, T. Eifler, R. J. Foley, C. H. Heinrich, E. Krause, S. Perlmutter, A. Pisani, _et al._ , arXiv:1904.01174 [astro-ph] (2019).
  * [Eifler  _et al._(2021a)Eifler, Miyatake, Krause, Heinrich, Miranda, Hirata, Xu, Hemmati, Simet, Capak, Choi, Doré , Doux, Fang, Hounsell, Huff, Huang, Jarvis, Kruk, Masters, Rozo, Scolnic, Spergel, Troxel, von der Linden, Wang, Weinberg, Wenzl, and Wu] T. Eifler, H. Miyatake, E. Krause, C. Heinrich, V. Miranda, C. Hirata, J. Xu, S. Hemmati, M. Simet, P. Capak, A. Choi, O. Doré , C. Doux, X. Fang, R. Hounsell, E. Huff, H.-J. Huang, M. Jarvis, J. Kruk, D. Masters, E. Rozo, D. Scolnic, D. N. Spergel, M. Troxel, A. von der Linden, Y. Wang, D. H. Weinberg, L. Wenzl, and H.-Y. Wu, Monthly Notices of the Royal Astronomical Society 507, 1746 (2021a).
  * [Eifler  _et al._(2021b)Eifler, Simet, Krause, Hirata, Huang, Fang, Miranda, Mandelbaum, Doux, Heinrich _et al._] T. Eifler, M. Simet, E. Krause, C. Hirata, H.-J. Huang, X. Fang, V. Miranda, R. Mandelbaum, C. Doux, C. Heinrich, _et al._ , Monthly Notices of the Royal Astronomical Society , stab533 (2021b).
  * [Ishak  _et al._(2019)Ishak, Baker, Bull, Pedersen, Blazek, Ferreira, Leonard, Lin, Linder, Pardo, and Valogiannis] M. Ishak, T. Baker, P. Bull, E. M. Pedersen, J. Blazek, P. G. Ferreira, C. D. Leonard, W. Lin, E. Linder, K. Pardo, and G. Valogiannis, Modified gravity and dark energy models beyond $w(z)$cdm testable by lsst (2019), arXiv:1905.09687 [astro-ph.CO] .
  * [Planck Collaboration(2016)] Planck Collaboration, A&A 594, A14 (2016).
  * [Raveri  _et al._(2021)Raveri, Pogosian, Koyama, Martinelli, Silvestri, Zhao, Li, Peirone, and Zucca] M. Raveri, L. Pogosian, K. Koyama, M. Martinelli, A. Silvestri, G.-B. Zhao, J. Li, S. Peirone, and A. Zucca, A joint reconstruction of dark energy and modified growth evolution (2021), arXiv:2107.12990 [astro-ph.CO] .
  * [Linder and Huterer(2005)] E. V. Linder and D. Huterer, Phys. Rev. D 72, 043509 (2005).
  * [Asaba  _et al._(2013)Asaba, Hikage, Koyama, Zhao, Hojjati, and Pogosian] S. Asaba, C. Hikage, K. Koyama, G.-B. Zhao, A. Hojjati, and L. Pogosian, Journal of Cosmology and Astroparticle Physics 2013 (08), 029–029.
  * [Kohonen(1982)] T. Kohonen, Biological Cybernetics 43, 59 (1982).
  * [Maehoenen and Hakala(1995)] P. H. Maehoenen and P. J. Hakala, The Astrophysical Journal Letters 452, L77 (1995).
  * [Hemmati  _et al._(2019)Hemmati, Capak, Pourrahmani, Nayyeri, Stern, Mobasher, Darvish, Davidzon, Ilbert, Masters, and Shahidi] S. Hemmati, P. Capak, M. Pourrahmani, H. Nayyeri, D. Stern, B. Mobasher, B. Darvish, I. Davidzon, O. Ilbert, D. Masters, and A. Shahidi, The Astrophysical Journal 881, L14 (2019).
  * [Masters  _et al._(2015)Masters, Capak, Stern, Ilbert, Salvato, Schmidt, Longo, Rhodes, Paltani, Mobasher _et al._] D. Masters, P. Capak, D. Stern, O. Ilbert, M. Salvato, S. Schmidt, G. Longo, J. Rhodes, S. Paltani, B. Mobasher, _et al._ , The Astrophysical Journal 813, 53 (2015).
  * [Wright  _et al._(2020)Wright, Hildebrandt, van den Busch, and Heymans] A. H. Wright, H. Hildebrandt, J. L. van den Busch, and C. Heymans, Astronomy & Astrophysics 637, A100 (2020).
  * [Myles  _et al._(2021)Myles, Alarcon, Amon, Sánchez, Everett, DeRose, McCullough, Gruen, Bernstein, Troxel _et al._] J. Myles, A. Alarcon, A. Amon, C. Sánchez, S. Everett, J. DeRose, J. McCullough, D. Gruen, G. M. Bernstein, M. A. Troxel, _et al._ , Monthly Notices of the Royal Astronomical Society 505, 4249 (2021).
  * [Moosavi  _et al._(2014)Moosavi, Packmann, and Vallés] V. Moosavi, S. Packmann, and I. Vallés, Sompy: A python library for self organizing map (som) (2014), gitHub.[Online]. Available: https://github. com/sevamoo/SOMPY.
  * [Davidzon  _et al._(2019)Davidzon, Laigle, Capak, Ilbert, Masters, Hemmati, Apostolakos, Coupon, de la Torre, Devriendt, Dubois, Kashino, Paltani, and Pichon] I. Davidzon, C. Laigle, P. L. Capak, O. Ilbert, D. C. Masters, S. Hemmati, N. Apostolakos, J. Coupon, S. de la Torre, J. Devriendt, Y. Dubois, D. Kashino, S. Paltani, and C. Pichon, Monthly Notices of the Royal Astronomical Society 489, 4817 (2019).
  * [Zuntz  _et al._(2015)Zuntz, Paterno, Jennings, Rudd, Manzotti, Dodelson, Bridle, Sehrish, and Kowalkowski] J. Zuntz, M. Paterno, E. Jennings, D. Rudd, A. Manzotti, S. Dodelson, S. Bridle, S. Sehrish, and J. Kowalkowski, Astronomy and Computing 12, 45 (2015).
  * [Zucca  _et al._(2019)Zucca, Pogosian, Silvestri, and Zhao] A. Zucca, L. Pogosian, A. Silvestri, and G.-B. Zhao, Journal of Cosmology and Astroparticle Physics 2019 (05), 001.
  * [Hojjati  _et al._(2011)Hojjati, Pogosian, and Zhao] A. Hojjati, L. Pogosian, and G.-B. Zhao, Journal of Cosmology and Astroparticle Physics 2011 (08), 005.
  * [Zhao  _et al._(2009)Zhao, Pogosian, Silvestri, and Zylberberg] G.-B. Zhao, L. Pogosian, A. Silvestri, and J. Zylberberg, Physical Review D 79, 083513 (2009).
  * [Hu and Sawicki(2007)] W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007).
  * [Song  _et al._(2007)Song, Hu, and Sawicki] Y.-S. Song, W. Hu, and I. Sawicki, Phys. Rev. D 75, 044004 (2007).
  * [Dossett  _et al._(2014)Dossett, Hu, and Parkinson] J. Dossett, B. Hu, and D. Parkinson, Journal of Cosmology and Astroparticle Physics 2014 (03), 046.
  * [Damour and Polyakov(1994)] T. Damour and A. Polyakov, Nuclear Physics B 423, 532–558 (1994).
  * [Brax  _et al._(2010)Brax, van de Bruck, Davis, and Shaw] P. Brax, C. van de Bruck, A.-C. Davis, and D. Shaw, Phys. Rev. D 82, 063519 (2010).
  * [Hinterbichler and Khoury(2010)] K. Hinterbichler and J. Khoury, Phys. Rev. Lett. 104, 231301 (2010).
  * [Vazsonyi  _et al._(2021)Vazsonyi, Taylor, Valogiannis, Ramachandra, Ferté, and Rhodes] L. Vazsonyi, P. L. Taylor, G. Valogiannis, N. S. Ramachandra, A. Ferté, and J. Rhodes, Phys. Rev. D 104, 083527 (2021).
  * [Jain  _et al._(2013)Jain, Vikram, and Sakstein] B. Jain, V. Vikram, and J. Sakstein, The Astrophysical Journal 779, 39 (2013).
  * [Desmond and Ferreira(2020)] H. Desmond and P. G. Ferreira, Phys. Rev. D 102, 104060 (2020).
  * [Krause and Eifler(2017)] E. Krause and T. Eifler, Monthly Notices of the Royal Astronomical Society 470, 2100–2112 (2017).
  * [Fang  _et al._(2020a)Fang, Krause, Eifler, and MacCrann] X. Fang, E. Krause, T. Eifler, and N. MacCrann, Journal of Cosmology and Astroparticle Physics 2020 (05), 010–010.
  * [Fang  _et al._(2020b)Fang, Eifler, and Krause] X. Fang, T. Eifler, and E. Krause, Monthly Notices of the Royal Astronomical Society 497, 2699–2714 (2020b).
  * [Linder(2003)] E. V. Linder, Phys. Rev. Lett. 90, 091301 (2003).
  * [Chevallier and Polarski(2001)] M. Chevallier and D. Polarski, International Journal of Modern Physics D 10, 213–223 (2001).
  * [Lewis  _et al._(2000)Lewis, Challinor, and Lasenby] A. Lewis, A. Challinor, and A. Lasenby, ApJ 538, 473 (2000).
  * [Howlett  _et al._(2012)Howlett, Lewis, Hall, and Challinor] C. Howlett, A. Lewis, A. Hall, and A. Challinor, Journal of Cosmology and Astroparticle Physics 2012 (04), 027.
  * [Lewis(2014)] A. Lewis (2014), https://cosmologist.info/notes/CAMB.pdf.


# Opening band gaps of low-dimensional materials at the meta-GGA level of
density functional approximations

Bimal Neupane<EMAIL_ADDRESS>Hong Tang<EMAIL_ADDRESS>Niraj
K. Nepal  Santosh Adhikari  Adrienn Ruzsinszky  Department of Physics, Temple
University, Philadelphia, Pennsylvania 19122, United States

###### Abstract

The quasiparticle band structure can be properly described by the GW
approximation, at an increased computational cost. Semilocal density
functionals up to the generalized gradient approximation (GGA) level cannot
compete with the accuracy of hybrid-based approximations or GW. Meta-GGA
density functionals with a strong dependence on the kinetic energy density
ingredient can potentially give wider band gaps compared to GGA’s. The recent
TASK meta-GGA density functional [Phys. Rev. Research, 1, 033082 (2019)] is
constructed with an enhanced nonlocality in the generalized Kohn-Sham scheme,
and therefore harbors great opportunities for band gap prediction. Although
this approximation was found to yield excellent band gaps of bulk solids, this
accuracy cannot be straightforwardly transferred to low-dimensional materials.
The reduced screening of these materials results in larger band gaps compared
to their bulk counterparts, as an additional barrier to overcome. In this work
we demonstrate how the alteration of exact physical constraints in this
functional affects the band gaps of monolayers and nanoribbons, and present
accurate band gaps competing with the HSE06 approximation. In order to achieve
this goal, we have modified the TASK functional (a) by changing the tight
upper-bound for one or two-electron systems ($h_{X}^{0}$) from 1.174 to 1.29
(b) by changing the limit of interpolation function
$f_{X}(\alpha\rightarrow\infty$) of the TASK functional that interpolates the
exchange enhancement factor $F_{X}(s,\alpha)$ from $\alpha=$ 0 to 1. The
resulting modified TASK (mTASK) was tested for various materials from 3D to 2D
to 1D (nanoribbons), and was compared with the results of the higher-level
hybrid functional HSE06 or with the G0W0 approximation within many-body
perturbation theory. We find that mTASK greatly improves the band gaps and
band structures of 2D and 1D systems, without significantly affecting the
accuracy of the original TASK for the bulk 3D materials, when compared to the
PBE-GGA and SCAN meta-GGA. We further demonstrate the applicability of mTASK
by assessing the band structures of TMD nanoribbons with respect to various
bending curvatures.

##  I Introduction

Due to their reduced dimensionality, two-dimensional (2D) materials (e.g.,
transition metal dichalcogenides or TMDs) exhibit an extraordinary optical
response Ruppert  _et al._ (2014); Ugeda  _et al._ (2014) in comparison with
bulk counterparts, as has been shown early via the examples of graphene and 2D
MoS2. The spatial confinement and reduced dielectric screening of 2D materials
causes strong Coulomb interaction that allows more stable exciton formation
with large binding energy and oscillator strength compared to bulk crystals.
These attractive features can in principle be harvested for optoelectronics.
Fine tuning of the optical properties can be achieved by mechanical bending
that alters the electronic structure. Optical nanodevice functionality can
require the optimization of exciton binding energy, optical absorption and
most importantly the fundamental band gap with respect to the strains existing
in the bending space. Strain engineering Castellanos-Gomez  _et al._ (2013)
has been known as a tool to control the electronic properties of 2D materials,
but the impact of bending is not yet thoroughly explored.

Fundamental band gaps carry a great relevance from a theoretical and practical
view Perdew  _et al._ (2017). Since the fundamental band gap is the unbound
limit of an exciton series, it controls the optoelectronic response of
materials. Optical absorption spectra within time dependent DFT (TDDFT) are
usually evaluated on a scissor-operator-corrected DFT band gap Levine and
Allan (1989); Del Sole and Girlanda (1993). This approach is simpler than
applying the more expensive GW approximation, but the scissor-shift correction
is often evaluated from experiments, and the correction is very unlikely to be
available for strained or bent structures relevant for industrial
applications. In bulk crystals, GW yields a nearly constant shift to the
fundamental gap of semilocal DFT even for strained structures, but this is not
necessarily true any more for low-dimensional systems Spataru and Léonard
(2013); Wagner  _et al._ (2019). This difference is a consequence of enhanced
many-body effects present at low-dimensionality. The literature reports both
enhanced and weaker band gap change with GW compared to Kohn-Sham DFT, for the
fundamental and higher band gaps of stretched or compressed carbon nanotubes
Levine and Allan (1989). A similar effect can be expected for nanoribbons. As
an alternative to the missing scissor corrections for the bent structures, we
can use the hybrid functional HSE06 Heyd  _et al._ (2003, 2006) in a
generalized Kohn-Sham scheme to estimate the band gaps and bandstructures for
various bending curvatures. HSE06 is known to yield reasonably accurate band
gaps for bulk crystals. Less information is available about the accuracy of
band gaps of nanoribbons with HSE06, and the increased size of supercells can
make such calculations more demanding. The demand for an alternative
approximation with more computational feasibility is large, but semilocal
density functionals are far from being reliable for accurate band gap
prediction. The underestimation of band gap in semilocal functionals is a
consequence of the lack of derivative discontinuity within the Kohn-Sham
potential Perdew  _et al._ (1982); Perdew and Levy (1983); Sham and Schlüter
(1983); Perdew (1985). Accurate band gaps require an effective potential that
is either a discontinuous multiplicative operator or a continuous nonlocal
operator. Early band-structure calculations showed that the LSDA band gaps for
semiconductors were often about half the measured fundamental energy gaps.
Meta-GGA density functional approximatons are placed on the third rung of
density functional approximations Perdew  _et al._ (2005). Meta-GGA’s are
explicit functionals of the Kohn-Sham orbitals and implicit functionals of the
density. Some meta-GGA’s were already found to open band gaps more than GGA’s
Sun  _et al._ (2015). The successful SCAN Sun  _et al._ (2015) meta-GGA
usually yields a slight improvement for bulks solids, but it is still far from
approaching the accuracy of HSE06 or GW. Alternatively the recently developed
TASK meta-GGA Aschebrock and Kümmel (2019) for band gaps of bulk solids in a
generalized Perdew  _et al._ (2017); Seidl  _et al._ (1996) Kohn-Sham scheme
is a promising approximation. The TASK meta-GGA is not designed to yield
comparable accuracy for ground-state properties; capturing accurate ground-
state properties and accurate electronic structures simultaneously is not to
be expected from the first three rungs of Jacob’s ladder of density functional
approximations. Nevertheless, our focus in this work is on band gaps and not
on ground-state properties. TASK and similar meta-GGA’s, implemented within
the generalized Kohn-Sham scheme, could provide improved estimates of
quasiparticle energies. Further hints refer to the TASK meta-GGA as a
potential tool for excitonic peaks in the spirit of Nazarov and Vignale
Nazarov and Vignale (2011), taking the second functional derivative of the
exchange-correlation energy. With the consideration of meta-GGA’s for band
gaps we follow an approach like that of the Tran-Blaha potential functional
Tran and Blaha (2009), with the change that the potential is properly a
functional derivative of the exchange-correlation energy. While the TASK
functional has been tested on a set of bulk solids with great accuracy as a
result, its performance is completely unknown for low-dimensional materials.
By now, tremendous numbers of references point out the enhanced optical
response from the different screening in low-dimensional materials compared to
bulk solids Ugeda  _et al._ (2014); Qiu  _et al._ (2013); Ramasubramaniam
(2012); He  _et al._ (2014).

##  II Methods

Typical meta-GGA’s utilize the kinetic energy density as an additional
ingredient beyond the ones of GGA’s Perdew  _et al._ (2005); Sun  _et al._
(2015); Aschebrock and Kümmel (2019). The kinetic energy density makes the
meta-GGA form an implicit density functional, in principle, with potential
Perdew  _et al._ (2017) for band gap prediction and optical spectra. Accurate
band gap prediction is based on the derivative discontinuity Perdew  _et al._
(1982); a feature that semilocal density functional approximations are
missing. The derivative discontinuity manifests itself as a step structure in
the exact exchange-correlation potential as Kümmel and Perdew (2003)

$\Delta_{X}={v_{X}\left(\textbf{r}\right)|}_{+}-{v_{X}\left(\textbf{r}\right)|}_{-}{=\frac{{\delta
E}_{X}\left[n\right]}{\delta n\left(\textbf{r}\right)}|}_{+}-{\frac{{\delta
E}_{X}\left[n\right]}{\delta n\left(\textbf{r}\right)}|}_{-,}$ (1)

where the positive and negative signs refer to the left and right-hand sides
of the potential at integer electron numbers. Perdew and Levy showed that the
functional derivative discontinuity of the exchange-correlation density
functional plays a crucial role in the correct prediction of band gaps Perdew
and Levy (1983). The functional derivative of a meta-GGA within the Kohn-Sham
scheme contains the $\tau$-dependence in the third term below that increases
the derivative discontinuity Eich and Hellgren (2014)

$\frac{{\delta E}_{X}^{mGGA}\left[n\right]}{\delta
n\left(\textbf{r}\right)}=\frac{\delta}{\delta
n\left(\textbf{r}\right)}\int{e_{X}\left(n,\nabla
n,\tau\right)}d^{3}r^{\prime}=\frac{\partial e_{X}}{\partial
n}\left(\textbf{r}\right)-\nabla\cdot\left[\frac{\partial
e_{X}}{\partial\nabla n}\left(\textbf{r}\right)\right]+\int{\frac{\partial
e_{X}}{\partial\tau}\left(\textbf{r}^{\prime}\right)\frac{\partial\tau\left(\textbf{r}^{\prime}\right)}{\partial
n\left(\textbf{r}\right)}}d^{3}r^{\prime}$ (2)

To date, only a few meta-GGAs have been identified to display the necessary
dependence on the kinetic energy density Sun  _et al._ (2015). A recent effort
Aschebrock and Kümmel (2019) demonstrates that meta-GGA’s with an enhanced
dependence on the dimensionless kinetic energy density ingredient
$\alpha=\frac{\left(\tau-\tau^{w}\right)}{\tau^{unif}}$, give accurate band
gaps and can be applied to optical properties. Band gap calculations with the
TASK meta-GGA have been done within the generalized Kohn-Sham scheme, in which
the exchange-correlation potential is a differential operator with significant
nonlocality Arbuznikov and Kaupp (2003). The dimensionless iso-orbital
indicator $\alpha$ contains $\tau^{w}$, the kinetic energy of a one-electron
system and $\tau^{unif}$,the kinetic energy density of the homogeneous
electron gas. The condition that is responsible for an improved band gap as
well as for excitonic peaks is:

$\frac{\partial e_{X}}{\partial\tau}>0$ (3)

$\frac{\partial e_{X}}{\partial\tau}$ is a $\tau$-dependent factor that
contributes to a first or second functional derivative of the exchange-
correlation potential or the corresponding energy. For practical use it is
easier to show the $\tau$-dependence of the exchange enhancement factor
$F_{X}$. A sizeable nonlocality leads to the condition Nazarov and Vignale
(2011):

$\frac{\partial F_{X}}{\partial\alpha}\ll 0.$ (4)

The recent SCAN meta-GGA was found to exhibit an observable opening of band
gaps of various bulk solids. Still, SCAN is a ground-state density functional.
The TASK meta-GGA approximation is a thoughtful re-construction of the
exchange form of SCAN so that with the increased slope $\frac{\partial
F_{X}}{\partial\alpha}$ more nonlocality in the exchange potential can be
achieved. The TASK exchange basically keeps the exchange $F_{X}$ of SCAN

$F_{X}^{TASK}\left(s,\
\alpha\right)=h_{X}^{0}g_{X}\left(s\right)+\left[1-f_{X}(\alpha)\right]\left[h_{X}^{1}\left(s\right)-h_{X}^{0}\right]\left[g_{X}\left(s\right)\right]^{d}$
(5)

In the expression of $F_{X}$, $h_{X}^{0}$ is the tight upper-bound for one or
two-electron systems. The same $h_{X}^{0}$ was applied in the SCAN functional.
For $\alpha=$0, $F_{X}^{TASK}\leq 1.174$, which is also the conjectured bound
for this limit Perdew  _et al._ (2014). $f_{X}(\alpha)$ is a function that
satisfies the fourth-order gradient limit and interpolates between $\alpha$=0
and $\alpha=$1\. When $f_{X}\left(\alpha=0\right)=$1, the TASK functional
recovers SCAN by satisfying the two-dimensional scaling Perdew  _et al._
(2014); Levy (1991); Pollack and Perdew (2000) of the reduced gradient
$s=\frac{\left|\nabla n\right|}{2\left(3\pi^{2}\right)^{1/3}n^{4/3}}$ as,

${\ \ g}_{X}\left(s\right)=1-e^{\left(-cs^{\frac{1}{2}}\right)}.$ (6)

In equation (6), c is a parameter that was found by fitting to the exact
atomic energy of hydrogen. TASK deviates from SCAN in the construction of
$h_{X}^{1}\left(s\right)$. The construction of $h_{X}^{1}\left(s\right)$ and
$f_{X}(\alpha)$ aims to increase the slope defined by equation (4). The
fourth-order gradient expansion around $s=$0 and $\alpha=$1 recovered by both
$h_{X}^{1}\left(s\right)$ and $f_{X}(\alpha)$ together is

${F}_{X}^{GE4}\left(s,\alpha\right)\sim
1+\mu_{s}s^{2}+\mu_{\alpha}\left(\alpha-1\right)+C_{s}s^{4}+C_{s\alpha}s^{2}\left(\alpha-1\right)+C_{\alpha}\left(\alpha-1\right)^{2}+\mathcal{O}\left(\nabla\right)^{6}$
(7)

along with a fourth-order expansion for $\tau$ Brack  _et al._ (1976); Perdew
_et al._ (1986), $\mu_{s}$, $\mu_{\alpha}$, $C_{s}$, and $C_{s\alpha}$ are
chosen to correspond to the coefficients of the Taylor expansions for
$h_{X}^{1}\left(s\right)$ and $f_{X}\left(\alpha\right)$ to fourth-order,
respectively. All these coefficients and the exact conditions for
$f_{X}\left(\alpha\right)$ lead to eight equations for the eight coefficients
of Chebyshev expansions Boyd (1987); Press  _et al._ (1992) for
$h_{X}^{1}\left(s\right)$ and $f_{X}\left(\alpha\right)$ to fourth-order.
$\mu_{\alpha}$ is determined from a quadratic equation with the ingredient of
$h_{X}^{0}$. The more negative value of $\mu_{\alpha}$ yields an increased
nonlocality compared to any previous meta-GGA’s, so that TASK delivers
excellent fundamental band gaps for bulk solids with a large enough variety of
structures.

In the last decade, utilizing methodological and computational advances Gygi
(2009), hybrid functionals Bickelhaupt and Baerends (2000) have been
increasingly used to investigate a variety of periodic systems with plane wave
basis sets. Among hybrid functionals the HSE06 approximation Heyd  _et al._
(2003, 2006) is particularly popular, but some other hybrid approximations
Yanai  _et al._ (2004); Weintraub  _et al._ (2009); Adamo and Barone (1999);
Skone  _et al._ (2014); Refaely-Abramson  _et al._ (2015) with various ratios
of exact exchange admixture have also become beneficial for the condensed
matter community. The most successful hybridization schemes are based on
partitioning the Coulomb operator into short- and long-range components in a
two-parameter form in which the parameters $\alpha$ and $\beta$ control the
mixing ratio of long- and short-range exchange. This scheme is the basis of
various popular density functionals, such as CAM-B3LYP Yanai  _et al._ (2004),
LC-$\omega$PBE Weintraub  _et al._ (2009). HSE06 with $\beta=$ 0.25 and
$\omega=$ 0.11 Bohr-1 is a short-range screened hybrid that recovers PBE0 for
the short range therefore enabling computational efficiency for periodic
systems. The mixing of exact exchange in hybrid functionals is obviously
linked to spatial nonlocality. The hybrid parameter $\alpha$ is often treated
as an adjustable parameter to reproduce the experimental band gap of solids
Pozun and Henkelman (2011); Conesa (2012); Alkauskas  _et al._ (2008, 2011);
Broqvist  _et al._ (2010). Since in semiconductors and insulators the
screening of the long-range tail of the Coulomb interaction was found to be
proportional to the inverse of the static dielectric constant
${(\epsilon}_{\infty}^{-1})$, it has become natural to link $\alpha$ to
$\epsilon_{\infty}^{-1}$. It should be also noted that an analogy was
established between hybrid functionals and the static COulomb Hole plus
Screened EXchange (COHSEX) approximation Hedin (1965), in which the screened
Coulomb interaction is expressed as

$W(\mathbf{r},\mathbf{r}^{\prime})\approx\frac{1}{\epsilon_{\infty}}v(\mathbf{r},\mathbf{r}^{\prime})$
(8)

Our work, although in a completely different way, seeks the same spatial
nonlocality needed for band gaps, considering the reduced screening with
increasing low-dimensionality. Semilocal density functionals such as GGA’s do
not exhibit spatial nonlocality, but through $\alpha$ meta-GGA’s Becke and
Edgecombe (1990); Savin  _et al._ (1992, 1997) are implicit functionals of the
density, and can overcome this limitation Sun  _et al._ (2015, 2012, 2013) .
The dimensionless ingredient $\alpha$ in the TASK and SCAN functional is an
explicit functional of the Kohn-Sham orbitals, and therefore carries potential
spatial nonlocality. Note that the notational coincidence with the hybrid
functionals’ $\alpha$, traditionally called, is originally accidental, but now
gains some physical analogy. Semiconductors have a significantly lesser
screening than metals, as a consequence of their less-localized exchange-
correlation hole. Meta-GGA functionals based on the Laplacian of the electron
density and not on the kinetic energy density can better describe the more
localized exchange-correlation hole of metals than the one of semiconductors
Mejia-Rodriguez and Trickey (2018); Mejía-Rodríguez and Trickey (2020); Mejia-
Rodriguez and Trickey (2017); Kaplan and Perdew . Meta-GGA’s with gradient and
Laplacian-only ingredients were found more accurate for metals than for semi-
conductors, while SCAN and other meta-GGA’s with orbitals in their $\alpha$ or
$\tau$-dependence work better for semiconductors and insulators. With a proper
modification of the original TASK functional the slope $\frac{\partial
F_{X}}{\partial\alpha}$ can be further increased so that the modified TASK
(mTASK) has screening appropriate for low-dimensional materials.

In Section 4 we will discuss how the increased nonlocality with the
corresponding increased screening affects 3D, 2D and 1D materials, when the
fundamental band gaps of these materials with various dimensionality are
compared to HSE06. In our modification (modified TASK or mTASK) we focus on
the coefficient $\mu_{\alpha}$. An increased slope is equivalent to increased
$\alpha$-dependence of $F_{X}$. Based on the fourth-order gradient expansion
of $F_{X}$, $\mu_{\alpha}$ is determined from a quadratic equation with
$h_{X}^{0}$ as an ingredient

$\mu_{\alpha}^{\pm}=-\frac{97+3h_{X}^{0}\pm\sqrt{9\left(h_{X}^{0}\right)^{2}+74166h_{X}^{0}-64175}}{1200}.$
(9)

This TASK expression is able to deliver a large-enough nonlocality
$\mu_{\alpha}^{+}=-0.209897$ for bulk solids, but not for low-dimensional
materials. We realize that the increased band gap from the reduced screening
in monolayers and nanoribbons requires a different $\mu_{\alpha}^{+}$ and
therefore different $h_{X}^{0}$. Lacking more exact constraints for the ground
state, with any choice for ${F_{X}\leq h}_{X}^{0}\neq 1.174$, we have to
sacrifice the tight Lieb-Oxford bound. In this work, we choose $F_{X}\leq
h_{X}^{0}=1.29$. This choice is far from the Lieb-Oxford bound $F_{X}\leq
1.804$ known for GGA’s. This condition is explained in Ref. Perdew  _et al._
(2014) for meta-GGA functionals. To make the exchange enhancement factor a
smooth and monotonically decreasing function of $\alpha$ for any value of s,
we choose the limiting condition $f_{X}\left(\alpha\rightarrow\
\infty\right)=\ -3.5$. From Figure 1, it can be seen that $\frac{\partial
F^{mTASK}_{X}}{\partial\alpha}\leq\frac{\partial
F^{TASK}_{X}}{\partial\alpha}$ for any value of s. We expect the inclusion of
more nonlocality in the mTASK exchange-correlation potential compared to TASK
provides better band gaps and band structures for 2D materials and
nanoribbons. In addition, when choosing the value of $h_{X}^{0}=1.804$ known
for GGA’s, we have observed that the exchange enhancement factor is not a
smooth function of $\alpha$.

Figures 1 and 2 demonstrate the changes in mTASK compared to TASK. Figure 1
presents $F_{X}$ for s and $\alpha$ respectively. $F_{X}$ of mTASK starts out
at $h_{X}^{0}=$ 1.29, compared to $h_{X}^{0}=$ 1.174 of TASK. The larger
$\alpha$\- and s- dependence of mTASK is evident in Figure 1. The three-
dimensional contour plot of $F_{X}$ (Figure 2) allows showing all s and
$\alpha$-dependence.

Figure 1: The enhancement factors $F_{X}$ for the TASK (upper panel) and
modified TASK (lower panel) as functions of s and $\alpha$, respectively.
Figure 2: The three-dimensional contour plot of the exchange enhancement
factor $F_{X}$ as a function of $\alpha$ and s for the original TASK (purple,
$h_{X}^{0}=$ 1.174, $\mu_{\alpha}=-$ 0.209897,
$f_{X}\left(\alpha\rightarrow\infty\right)=-$ 3 ) and the modified TASK
(yellow, $h_{X}^{0}=$ 1.29, $\mu_{\alpha}=-$ 0.231993,
$f_{X}\left(\alpha\rightarrow\infty\right)=-$ 3.5).

##  III Computational Details

All calculations were performed in the Vienna ab initio simulations package
(VASP) Kresse and Furthmüller (1996); Kresse and Joubert (1999). The valence
electrons of all elements, are treated by the projector augmented wave (PAW)
pseudo-potential method Blöchl (1994), which is recommended in the VASP
manual. The pseudopotential for tungsten atom with valence electron
configuration 6S15D5 was utilized. The plane-wave energy cut-off is set as 450
eV for all calculations and it leads to converged results. The Brillouin zone
is sampled by a Gamma centered mesh of 8$\times$1$\times$1 for both the
hexagonal armchair TMD nanoribbons and trigonal TMD nanoribbons. The Gamma
centered k-point mesh of 18$\times$18$\times$1 and 20$\times$20$\times$20 were
used for TMD monolayers and bulk solids, respectively. A vacuum layer of more
than 15 Å was inserted along both the width and thickness directions to
prevent interactions between periodic images due to the long-range Coulomb
interactions. The bent nanoribbons were constructed by fixing the distance
between the two edge metal atoms (i.e M in MX2). Two hydrogen atoms are
attached to every edge metal atom and one hydrogen atom is attached to every
edge X atom in the MX2 nanoribbon. The edge atoms for the bent nanoribbons
were only allowed to relax in the periodic direction. All other atoms were
fully relaxed until the force on each atom became less than 0.01 eV/Å. For the
bulk solids we used experimental lattice constants in the calculations, in
order to directly compare results with references. The PBE functional was used
to optimize (or relax) the nanoribbon structures, and we utilized them to
calculate the band structures using various DFT approximations. The lattice
parameter along the periodic direction of the nanoribbon was also relaxed.

##  IV Results

###  IV.1 Bulk crystals

Bulks solids and low-dimensional materials often exhibit dramatically
different physical properties, especially for optical response. To obtain a
general picture about the applicability of the approximations considered in
this work, first we have assessed several density functional approximations
including TASK and mTASK for bulk solids. Figure 3 represents a correlation
between the experimental and calculated band gaps of the same set of bulk
solids from Ref Aschebrock and Kümmel (2019). The underestimation of band gaps
from PBE is not surprising. As a semilocal density functional approximation,
PBE has only the ingredients of the local density and the gradient of the
local density, without explicit inclusion of nonlocal exchange effects, the
latter being important for an accurate description of band gaps. The SCAN
meta-GGA yields a systematic and slight improvement. As a meta-GGA functional,
SCAN can include some nonlocal exchange effects through the orbital dependent
ingredient $\alpha$. Although SCAN is largely accurate for equilibrium
structures and energies for various bonds, the nonlocality derived from
$\alpha$ may not be pronounced due to cancellation within the exchange and
correlation parts. TASK proves to be very accurate, almost in all regions of
band gaps except for the Ar crystal that has the largest band gap in this set
and is underestimated by TASK. The enhanced screening in mTASK results in a
slight overestimation of band gaps for bulk solids.

The effect of the correlation is exemplified through the SXPW and TXSC
approximations. In the former method SX refers to SCAN exchange, while the
correlation component is from the PW92 local spin density approximation Perdew
and Wang (1992). In the latter, TX is the TASK exchange, and SC is the
correlation of SCAN. SXPW and TXSC can be directly compared with SCAN and
TASK, respectively, because they have the same exchange but different
correlation approximations. In general, both these approximations open the
band gaps slightly more than SCAN, but they underestimate them compared to
TASK. This can serve as an evidence for the cancellation effects in SCAN’s
exchange and correlation parts for the nonlocality needed for band gap
description. The improvement of the band gaps from TASK and mTASK is more
consistent for the large-gap crystals such as MgO, LiCl, Kr, LiF, and Ar, as
can be seen in Figure 3 a and b.

Figure 3: Comparison of calculated and experimental band gaps of bulk solids
with various density functional approximations (a) PBE, SCAN, TASK, and mTASK
(b) SCAN, SXPW, TXSC, and TASK. TXSC is a functional with the exchange
component of TASK and the correlation from SCAN. SXPW refers to the exchange
part of SCAN and correlation of PW92.

###  IV.2 Band gaps of single layers of 2D TMDs

Figure 4: (a) Top view and (b) side view of monolayer 3$\times$3$\times$1
1H-TMDs. (c) Top view and (d) side view of monolayer 3$\times$3$\times$1 1T
TMDs. The larger balls represent M, while the smaller ones represent X in MX2.

Next, we have extended the application of our methods to monolayers of
transition metal dichalcogenides (TMDs). Two forms of structure of TMD
monolayers, namely, monolayer-hexagonal (1H) and monolayer-trigonal (1T), are
considered, as shown in Figure 4. We have assessed a range of meta-GGA
approximations with increasing nonlocality for the band gaps of monolayers of
TMDs. Beside the available experimental band gap values, the available results
from the screened hybrid functional HSE06 and the many-body perturbation G0W0
are also included for comparison. Table 1 shows the band gap results from
different methods for MoS2, MoSe2, WS2, and WSe2 1H monolayers. As can be
seen, PBE generally underestimates the band gaps of TMD monolayers, due to the
same reason mentioned for the case of bulk solids. Unlike for bulk solids,
SCAN remains almost the same quality as PBE for the 1H monolayers for band
gaps, as shown in Table 1. A noticeable difference is however, observed for
the TASK meta-GGA. TASK improves the values of band gaps upon PBE and SCAN for
both 1H and 1T monolayers. mTASK predicts a better band gap results than TASK
compared to PBE and SCAN. In mTASK, the delicate control is better realized
with the increased slope of $\frac{\partial F_{X}}{\partial\alpha}$ for the
two-dimensional layers, leading to an improved description of band gaps. HSE06
yields a noticeable overestimation for the band gaps of MoS2 and MoSe2. HSE06
is based on PBE0 for the short-range. TASK and mTASK both have a more delicate
set of ingredients than the short-range PBE component of HSE06 so that these
meta-GGA’s can be fine tuned to find the band gap.

Table 1: Band gaps (eV) of 1H TMD monolayers from PBE, SCAN, TASK, and mTASK density functional approximations in comparison with those from HSE06 and experiments. TMDs | PBE | SCAN | TASK | mTASK | HSE | Expt.
---|---|---|---|---|---|---
MoS2 | 1.66 | 1.64 | 1.79 | 1.80 | 2.02 Kang  _et al._ (2013) | 1.88 Mak  _et al._ (2010)
MoSe2 | 1.44 | 1.56 | 1.62 | 1.61 | 1.72 Kang  _et al._ (2013) | 1.57 Lu  _et al._ (2014)
WS2 | 1.81 | 1.79 | 1.94 | 1.98 | 1.98 Kang  _et al._ (2013) | 2.01 Tongay  _et al._ (2014)
WSe2 | 1.53 | 1.53 | 1.66 | 1.69 | 1.63 Kang  _et al._ (2013) | 1.67 Tonndorf  _et al._ (2013)

1T monolayers have not been synthetized, therefore no experimental band gaps
are available for these materials. Therefore, the hybrid HSE06 and G0W0@PBE
results have been utilized as references for comparison. Table 2 displays the
band gaps of HfS2, HfSe2, ZrS2 and ZrSe2 monolayers. The conclusion for 1T
monolayers is different from the one of 1H monolayers. The SCAN meta-GGA opens
the band gaps significantly compared to PBE. The TASK meta-GGA adds an even
much larger amount of nonlocality to SCAN, and opens the band gaps even
further. mTASK provides more refinement in slightly increased band gaps from
TASK. mTASK is not only matching HSE06’s accuracy, but predicts even better
band gaps for 1T TMD monolayers.

Table 2: Band gaps (eV) of 1T TMD monolayers from PBE, SCAN, TASK, mTASK density functional approximations in comparison with the available results from HSE06 and G0W0@PBE. TMDs | PBE | SCAN | TASK | mTASK | HSE06Zhao  _et al._ (2017) | G0W0Zhuang and Hennig (2013)
---|---|---|---|---|---|---
HfS2 | 1.30 | 1.61 | 2.22 | 2.35 | 2.40 | 2.45
HfSe2 | 0.63 | 0.93 | 1.51 | 1.59 | 1.32 | 1.39
ZrS2 | 1.17 | 1.53 | 2.00 | 2.17 | 2.16 | 2.56
ZrSe2 | 0.48 | 0.82 | 1.26 | 1.38 | 1.07 | 1.54

###  IV.3 Bulk and edge-state band gaps of flat and bent nanoribbons

Figure 5: The cross-section views of the 1H nanoribbons with different bending
curvatures (1/R). (a) flat ribbon, R = $\infty$, (b) R = 15 Å, (c) R = 10 Å.
The periodical direction of the nanoribbon is perpendicular to the page of
paper. All nanoribbons are hydrogen passivated on the two ribbon edges and
their geometrical structures are fully relaxed as shown in (a) through (c).
The blue ball represents the M atom, the green ball represents X in MX2, and
the white ball represents the H atom.

Mechanical bending can effectively control the conductivity of 2D
semiconducting nanoribbons, since bending can induce non-uniform strain in
nanoribbons. The magnitude of the induced strain by bending can be much larger
than that of uniaxial strain Yu  _et al._ (2016); Nepal  _et al._ (2019).
Beside the delocalized states, the localized edge states gain a high relevance
in the electronic structures of nanoribbon systems, while in monolayer systems
of the same constituents, the edge states are usually absent. Some of the
authors found earlier that the donor-like in-gap edge-states of armchair MoS2
nanoribbon and their associated harmful Fermi-level pinning can be removed to
some extent by bending Yu  _et al._ (2016).

The nature of edge states differs significantly in 1T and 1H nanoribbons Nepal
_et al._ (2019). The 1T nanoribbons only have the edge states below the Fermi
level while both the edge states above and below the Fermi level are present
in the 1H nanoribbons Nepal  _et al._ (2019). The position of these edges
states in the band structure can be critical for water splitting reactions
Chakrapani  _et al._ (2007); Zhuang and Hennig (2013), and any accurate
simulation of these band gaps can be critical to guide experiments.

1H semiconductors have a different band structure. The band structures of bent
1H TMDs nanoribbons evaluated with mTASK are similar to those from the hybrid
HSE06. Figure 6 shows the band structures of armchair 1H MoS2 nanoribbon
obtained from mTASK and HSE06 at different radii of curvature. As can be seen,
by increasing the bending curvature from flat (R = $\infty$) to R = 15 Å, the
upper edge band shifts towards the conduction band continuum slightly, and the
lower edge band begins to touch (or merge) with the valence band continuum at
R = 15 Å. Further increase in the curvature from R =15 Å to R = 10 Å makes the
upper edge state shift downward slightly. Also note that the minimum of the
conduction band continuum is shifted downwards. This makes an overall effect
that the upper edge band is gradually closer to the conduction band continuum
with an increase in curvature. However, the upper edge bands never merge into
the conduction band continuum. The position of the upper edge band is
comparable to HSE06 level calculation for mTASK where PBE, SCAN, and even TASK
fail to predict the correct position of the upper edge band.
As can be seen in figure 6, the lower edge state begins to merge into the bulk
valence band continuum at R=15 Å, and it is eventually completely embedded at
R = 10 Å. When the lower edge states merge into the valence band continuum,
the TMD nanoribbons become n-type semiconductors, given only the upper edge
bands are in the gap region and above the Fermi level (or the chemical
potential level). It is important to accurately predict the positions of those
edge states when the harmful Fermi-level pinning needs to be evaluated for
contact engineering, in applications of the semiconducting nanoribbons.
Furthermore, in photocatalytic applications, such as water splitting, it is
crucial that the shift of the edge states under bending is in favor of
efficiency, and preserves the photocatalytic properties of the semiconductors.
This feature appears in MoS2 and WS2 where it is highly preferred to keep the
band edges in the positions where they straddle the water redox potentials
Zhuang and Hennig (2013).
Figure 7 shows the edge band gaps of four 1H TMD nanoribbons (MoS2, MoSe2,
WS2, and WSe2) at different bending curvatures, calculated from PBE, SCAN,
TASK, mTASK and hybrid HSE06. Here, we define the band gap involving only the
delocalized-states as the “non-edge” or “bulk” band gap, otherwise, it is
simply referred as “edge band gap” or simply “band gap”. As shown in figure 7
(c) and (d), the more localized edge band gaps are very accurately captured by
TASK and mTASK for WS2, and WSe2. TASK and mTASK produce much better band gaps
of MoS2 and MoSe2 than PBE and SCAN, when compared with a reference of HSE06.
The accuracy and the computational gain are apparent from mTASK, as long as
sizeable nonlocality is found even at the meta-GGA level, with a reasonable
agreement with HSE06. Table 3 summarizes the bulk band gaps for all 1H
nanoribbons with all methods. The more localized edge-state band gaps of
Figure 7 are very accurately captured by TASK and mTASK for WS2, and WSe2, but
the agreement becomes less accurate for the more delocalized band gaps.

Figure 6: Band structures of armchair 1H MoS2 nanoribbon obtained from the mTASK density functional at different radii of curvature (a) flat R = $\infty$ (b) R = 15 Å (c) R = 12.5 Å (d) R = 10 Å and the hybrid HSE06 functional at different radii of curvature (e) flat R = $\infty$ (f) R = 15 Å (g) R = 12.5 Å (h) R = 10 Å. The edge states are indicted by solid red lines. Figure 7: The band gap of the edge states of 1H TMDs nanoribbons (a) MoS2. (b) MoSe2, (c) WS2, and (d) WSe2 at different bending curvatures with different density functional approximations. Bending curvature $\kappa=\frac{1}{R}$. Table 3: Delocalized (non-edge) band gaps of armchair 1H TMD nanoribbons at different radii of curvature (R in Å) from different density functional approximations. MX2 | R | PBE | SCAN | TASK | mTASK | HSE06
---|---|---|---|---|---|---
MoS2 | $\infty$ | 1.76 | 1.86 | 1.91 | 2.00 | 2.32
40 | 1.79 | 1.90 | 1.95 | 2.03 | 2.34
20 | 1.69 | 1.79 | 1.84 | 1.92 | 2.24
15 | 1.50 | 1.59 | 1.66 | 1.73 | 2.03
12.5 | 1.36 | 1.45 | 1.50 | 1.58 | 1.88
10 | 1.14 | 1.19 | 1.25 | 1.31 | 1.77
MoSe2 | $\infty$ | 1.57 | 1.70 | 1.77 | 1.76 | 2.05
40 | 1.60 | 1.72 | 1.80 | 1.79 | 2.08
20 | 1.55 | 1.67 | 1.74 | 1.74 | 2.03
15 | 1.42 | 1.53 | 1.57 | 1.62 | 1.90
12.5 | 1.14 | 1.27 | 1.33 | 1.37 | 1.66
10 | 0.90 | 0.99 | 1.02 | 1.07 | 1.34
WS2 | $\infty$ | 2.00 | 2.03 | 2.16 | 2.23 | 2.54
40 | 2.03 | 2.05 | 2.20 | 2.25 | 2.56
20 | 1.94 | 2.01 | 2.09 | 2.18 | 2.49
15 | 1.75 | 1.82 | 1.90 | 1.96 | 2.39
12.5 | 1.48 | 1.59 | 1.74 | 1.78 | 2.06
10 | 1.34 | 1.40 | 1.46 | 1.54 | 1.83
WSe2 | $\infty$ | 1.73 | 1.76 | 1.89 | 1.92 | 2.22
40 | 1.74 | 1.77 | 1.90 | 1.94 | 2.24
20 | 1.70 | 1.73 | 1.86 | 1.89 | 2.19
15 | 1.60 | 1.63 | 1.75 | 1.78 | 2.08
12.5 | 1.21 | 1.33 | 1.44 | 1.47 | 1.74
10 | 1.07 | 1.18 | 1.22 | 1.28 | 1.54

Similarly, Figure 8 plots the band gaps calculated with four semilocal density
functional approximations (PBE, SCAN, TASK, and mTASK) and hybrid HSE06 for
four 1T flat nanoribbons, namely, HfS2, HfSe2, ZrS2 and ZrSe2. It is
noticeable that TASK and mTASK open the band gaps of both non-edge and edge-
states more than SCAN or HSE06. This is consistent with the results of the 1T
TMD monolayers shown in Table 2.

Figure 8: Band gaps of four 1T flat nanoribbons with hydrogen passivation,
calculated from different density functional approximations. (a) displays the
non-edge band gaps, which are measured as the gap between the maximum of the
valence band continuum (excluding edge bands) and the minimum of the
conduction band continuum (excluding edge bands), at the $\Gamma$ point in the
band structure. (b) displays the edge band gaps, which are measured as the gap
between the maximum of valence bands (including edge bands) and the minimum of
conduction bands (including edge bands), at the $\Gamma$ point in the band
structure.

##  V Conclusions

We have systematically assessed the fundamental band gaps of semiconductors
starting from bulk crystals to 1D nanoribbons. For low-dimensional materials,
we have modified the TASK meta-GGA within the generalized Kohn-Sham framework
Yang  _et al._ (2016) in order to increase the nonlocality in the exchange-
correlation potential. The mTASK greatly improves the band gaps and is in
close agreement to the experimental or the hybrid level HSE06 density
functional for 2D single-layer and nanoribbon systems. This improvement for
the band gap is related to the slope of the exchange enhancement factor with
respect to the dimensionless orbital indicator ingredient. The physics of
choosing the nonlocality in mTASK is analogous with choosing the exchange
mixing factor in hybrid functionals. The global or screened hybrids PBE0 and
HSE06 apply the exchange mixing globally while the TASK and mTASK meta-GGA’s
change their $\frac{\partial F_{X}}{\partial\alpha}$ semilocally. The limits
of the global exchange mixing factor are the well known values of zero and
one, but except its upper bound of zero, no in-principle lower limit can be
determined for $\frac{\partial F_{X}}{\partial\alpha}$. Based on its orbital-
dependence, $\alpha$ of meta-GGA’s is more than a local ingredient, and meta-
GGA’s are more than just semilocal density functionals. TASK was already found
to yield excellent band gaps of bulk solids. Low-dimensional materials such as
monolayers and nanoribbons are known to exhibit strongly reduced screening by
their reduced dimensionality. First we have investigated to what extent the
accuracy of the TASK meta-GGA is transferable to 2D materials. The overall
performance of mTASK for the band gaps of low-dimensional materials is better
than that of TASK because of the increased nonlocality within the generalized
Kohn-Sham scheme. The agreement of the band gaps and the band structures for
the 2D single layers and nanoribbons quantitatively and qualitatively is the
solid evidence for the accuracy of the mTASK functional for band gaps.
Furthermore, the band structure obtained from mTASK for flat and bent
nanoribbons is comparable to the accuracy of HSE06 calculations. The mTASK
meta-GGA accurately predicts the position of edge-state bands for flat and
bent nanoribbons compared to the hybrid HSE06 functional.

As pointed out, TASK and mTASK have a great potential for optical response
properties. The same prediction applies to electronic structure simulations
with structures that are different from their equilibrium. The optical
response of strained and bent low-dimensional systems within time dependent
density functional theory (TDDFT) is of great interest for device
applications. TDDFT is an affordable tool to evaluate optical properties.
Although this statement is true in principle, the application of TDDFT is not
so straightforward at this time, especially for strained or bent 2D materials.
When TDDFT is applied to a bulk crystal, the optical spectrum can be evaluated
on a scissor-operator-corrected DFT band gap Levine and Allan (1989); Del Sole
and Girlanda (1993). This approach is simpler than applying the more expensive
GW approximation, but the scissor correction is often taken from experiments
and the correction is very unlikely to be available for the bent structures at
different bending curvatures. In bulk crystals, GW yields a nearly constant
shift to the fundamental gap of semilocal DFT even for strained structures,
but this is not necessarily true any more for low-dimensional systems such as
nanoribbons. This difference is a consequence of enhanced many-body effects
present in low-dimensional materials Wagner  _et al._ (2019). HSE06 can be
feasible for the bent nanoribbons, but the HSE06 reference makes TDDFT more
expensive. TASK or mTASK meta-GGA’s, are promising references for TDDFT that
can make TDDFT feasible for bent or strained structures where the scissor
correction may not be available. These methods can then become practical for
the evaluation of exciton binding energies on strained nanoribbons.

##  VI Acknowledgment

The authors thank Prof. John P. Perdew for useful comments on the manuscript.
This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Basic Energy Sciences, under Award Number DE-
SC0021263. The calculations were carried out on HPC resources supported in
part by the National Science Foundation through major research instrumentation
grant number 1625061 and by the US Army Research Laboratory under contract
number W911NF-16-2-0189.

## References

  * Ruppert  _et al._ (2014) C. Ruppert, O. B. Aslan, and T. F. Heinz, Nano Letters 14, 6231 (2014).
  * Ugeda  _et al._ (2014) M. M. Ugeda, A. J. Bradley, S.-F. Shi, H. Felipe, Y. Zhang, D. Y. Qiu, W. Ruan, S.-K. Mo, Z. Hussain, Z.-X. Shen, _et al._ , Nature Materials 13, 1091 (2014).
  * Castellanos-Gomez  _et al._ (2013) A. Castellanos-Gomez, R. Roldán, E. Cappelluti, M. Buscema, F. Guinea, H. S. van der Zant, and G. A. Steele, Nano Letters 13, 5361 (2013).
  * Perdew  _et al._ (2017) J. P. Perdew, W. Yang, K. Burke, Z. Yang, E. K. Gross, M. Scheffler, G. E. Scuseria, T. M. Henderson, I. Y. Zhang, A. Ruzsinszky, _et al._ , Proceedings of the National Academy of Sciences 114, 2801 (2017).
  * Levine and Allan (1989) Z. H. Levine and D. C. Allan, Physical Review Letters 63, 1719 (1989).
  * Del Sole and Girlanda (1993) R. Del Sole and R. Girlanda, Physical Review B 48, 11789 (1993).
  * Spataru and Léonard (2013) C. D. Spataru and F. Léonard, Physical Review B 88, 045404 (2013).
  * Wagner  _et al._ (2019) C. Wagner, J. Schuster, and A. Schleife, Physical Review B 99, 075140 (2019).
  * Heyd  _et al._ (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 118, 8207 (2003).
  * Heyd  _et al._ (2006) J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 124, 219906 (2006).
  * Perdew  _et al._ (1982) J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz Jr, Physical Review Letters 49, 1691 (1982).
  * Perdew and Levy (1983) J. P. Perdew and M. Levy, Physical Review Letters 51, 1884 (1983).
  * Sham and Schlüter (1983) L. J. Sham and M. Schlüter, Physical Review Letters 51, 1888 (1983).
  * Perdew (1985) J. P. Perdew, International Journal of Quantum Chemistry 28, 497 (1985).
  * Perdew  _et al._ (2005) J. P. Perdew, A. Ruzsinszky, J. Tao, V. N. Staroverov, G. E. Scuseria, and G. I. Csonka, The Journal of Chemical Physics 123, 062201 (2005).
  * Sun  _et al._ (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Physical Review Letters 115, 036402 (2015).
  * Aschebrock and Kümmel (2019) T. Aschebrock and S. Kümmel, Physical Review Research 1, 033082 (2019).
  * Seidl  _et al._ (1996) A. Seidl, A. Görling, P. Vogl, J. A. Majewski, and M. Levy, Physical Review B 53, 3764 (1996).
  * Nazarov and Vignale (2011) V.U. Nazarov and G. Vignale, Physical Review Letters 107, 216402 (2011).
  * Tran and Blaha (2009) F. Tran and P. Blaha, Physical Review Letters 102, 226401 (2009).
  * Qiu  _et al._ (2013) D. Y. Qiu, F.H. daJornada, and S. G. Louie, Physical Review Letters 111, 216805 (2013).
  * Ramasubramaniam (2012) A. Ramasubramaniam, Physical Review B 86, 115409 (2012).
  * He  _et al._ (2014) K. He, N. Kumar, L. Zhao, Z. Wang, K. F. Mak, H. Zhao, and J. Shan, Physical Review Letters 113, 026803 (2014).
  * Kümmel and Perdew (2003) S. Kümmel and J. P. Perdew, Physical Review Letters 90, 043004 (2003).
  * Eich and Hellgren (2014) F. Eich and M. Hellgren, The Journal of Chemical Physics 141, 224107 (2014).
  * Arbuznikov and Kaupp (2003) A. V. Arbuznikov and M. Kaupp, Chemical Physics Letters 381, 495 (2003).
  * Perdew  _et al._ (2014) J. P. Perdew, A. Ruzsinszky, J. Sun, and K. Burke, The Journal of Chemical Physics 140, 18A533 (2014).
  * Levy (1991) M. Levy, Physical Review A 43, 4637 (1991).
  * Pollack and Perdew (2000) L. Pollack and J. Perdew, Journal of Physics: Condensed Matter 12, 1239 (2000).
  * Brack  _et al._ (1976) M. Brack, B. Jennings, and Y. Chu, Physics Letters B 65, 1 (1976).
  * Perdew  _et al._ (1986) J. P. Perdew, V. Sahni, M. K. Harbola, and R. K. Pathak, Physical Review B 34, 686 (1986).
  * Boyd (1987) J. P. Boyd, Journal of Computational Physics 70, 63 (1987).
  * Press  _et al._ (1992) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, The Art of Scientific Computing 2, 1002 (1992).
  * Gygi (2009) F. Gygi, Physical Review Letters 102, 166406 (2009).
  * Bickelhaupt and Baerends (2000) F. M. Bickelhaupt and E. J. Baerends, Reviews in Computational Chemistry 15, 1 (2000).
  * Yanai  _et al._ (2004) T. Yanai, D. P. Tew, and N. C. Handy, Chemical Physics Letters 393, 51 (2004).
  * Weintraub  _et al._ (2009) E. Weintraub, T. M. Henderson, and G. E. Scuseria, Journal of Chemical Theory and Computation 5, 754 (2009).
  * Adamo and Barone (1999) C. Adamo and V. Barone, The Journal of Chemical Physics 110, 6158 (1999).
  * Skone  _et al._ (2014) J. H. Skone, M. Govoni, and G. Galli, Physical Review B 89, 195112 (2014).
  * Refaely-Abramson  _et al._ (2015) S. Refaely-Abramson, M. Jain, S. Sharifzadeh, J. B. Neaton, and L. Kronik, Physical Review B 92, 081204(R) (2015).
  * Pozun and Henkelman (2011) Z. D. Pozun and G. Henkelman, The Journal of Chemical Physics 134, 224706 (2011).
  * Conesa (2012) J. C. Conesa, The Journal of Physical Chemistry C 116, 18884 (2012).
  * Alkauskas  _et al._ (2008) A. Alkauskas, P. Broqvist, F. Devynck, and A. Pasquarello, Physical Review Letters 101, 106802 (2008).
  * Alkauskas  _et al._ (2011) A. Alkauskas, P. Broqvist, and A. Pasquarello, Physica Status Solidi (b) 248, 775 (2011).
  * Broqvist  _et al._ (2010) P. Broqvist, A. Alkauskas, and A. Pasquarello, Physica Status Solidi (a) 207, 270 (2010).
  * Hedin (1965) L. Hedin, Physical Review 139, A796 (1965).
  * Becke and Edgecombe (1990) A. D. Becke and K. E. Edgecombe, The Journal of Chemical Physics 92, 5397 (1990).
  * Savin  _et al._ (1992) A. Savin, O. Jepsen, J. Flad, O. K. Andersen, H. Preuss, and H. G. von Schnering, Angewandte Chemie International Edition in English 31, 187 (1992).
  * Savin  _et al._ (1997) A. Savin, R. Nesper, S. Wengert, and T. F. Fässler, Angewandte Chemie International Edition in English 36, 1808 (1997).
  * Sun  _et al._ (2012) J. Sun, B. Xiao, and A. Ruzsinszky, The Journal of Chemical Physics 137, 051101 (2012).
  * Sun  _et al._ (2013) J. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao, A. Ruzsinszky, G. I. Csonka, G. E. Scuseria, and J. P. Perdew, Physical Review Letters 111, 106401 (2013).
  * Mejia-Rodriguez and Trickey (2018) D. Mejia-Rodriguez and S. B. Trickey, Physical Review B 98, 115161 (2018).
  * Mejía-Rodríguez and Trickey (2020) D. Mejía-Rodríguez and S. B. Trickey, Physical Review B 102, 121109(R) (2020).
  * Mejia-Rodriguez and Trickey (2017) D. Mejia-Rodriguez and S. B. Trickey, Physical Review A 96, 052512 (2017).
  * (55) A. D. Kaplan and J. P. Perdew, work in progress .
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999).
  * Blöchl (1994) P. E. Blöchl, Physical Review B 50, 17953 (1994).
  * Perdew and Wang (1992) J. P. Perdew and Y. Wang, Physical review B 45, 13244 (1992).
  * Kang  _et al._ (2013) J. Kang, S. Tongay, J. Zhou, J. Li, and J. Wu, Applied Physics Letters 102, 012111 (2013).
  * Mak  _et al._ (2010) K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Physical Review Letters 105, 136805 (2010).
  * Lu  _et al._ (2014) X. Lu, M. I. B. Utama, J. Lin, X. Gong, J. Zhang, Y. Zhao, S. T. Pantelides, J. Wang, Z. Dong, Z. Liu, _et al._ , Nano Letters 14, 2419 (2014).
  * Tongay  _et al._ (2014) S. Tongay, W. Fan, J. Kang, J. Park, U. Koldemir, J. Suh, D. S. Narang, K. Liu, J. Ji, J. Li, _et al._ , Nano Letters 14, 3185 (2014).
  * Tonndorf  _et al._ (2013) P. Tonndorf, R. Schmidt, P. Böttger, X. Zhang, J. Börner, A. Liebig, M. Albrecht, C. Kloc, O. Gordan, D. R. Zahn, _et al._ , Optics Express 21, 4908 (2013).
  * Zhao  _et al._ (2017) Q. Zhao, Y. Guo, K. Si, Z. Ren, J. Bai, and X. Xu, Physica Status Solidi (b) 254, 1700033 (2017).
  * Zhuang and Hennig (2013) H. L. Zhuang and R. G. Hennig, The Journal of Physical Chemistry C 117, 20440 (2013).
  * Yu  _et al._ (2016) L. Yu, A. Ruzsinszky, and J. P. Perdew, Nano Letters 16, 2444 (2016).
  * Nepal  _et al._ (2019) N. K. Nepal, L. Yu, Q. Yan, and A. Ruzsinszky, Physical Review Materials 3, 073601 (2019).
  * Chakrapani  _et al._ (2007) V. Chakrapani, J. C. Angus, A. B. Anderson, S. D. Wolter, B. R. Stoner, and G. U. Sumanasekera, Science 318, 1424 (2007).
  * Yang  _et al._ (2016) Z.-h. Yang, H. Peng, J. Sun, and J. P. Perdew, Physical Review B 93, 205205 (2016).

Supplementary Materials

For slowly varying densities with s $\approx$ 0, and $\alpha\approx$ 1, one
can get the condition for $\mu_{\alpha}$ from the fourth-order gradient
expansions (GE4) as:

$\mu_{\alpha}^{\pm}=-\frac{97+3h_{X}^{0}\pm\sqrt{9\left(h_{X}^{0}\right)^{2}+74166h_{X}^{0}-64175}}{1200}$
(10)

For $h_{X}^{0}=$ 1.29 , $\mu_{\alpha}^{+}\approx-0.231993$ and
$\mu_{\alpha}^{-}\approx 0.063877$. The more negative solution is chosen in
mTASK (similar to TASK) to get sizeable nonlocality,
$\mu_{\alpha}^{+}\approx-0.231993$. The interpolation functions
$h_{X}^{1}\left(s\right)$ and $f_{X}\left(\alpha\right)$ are written as
Chebyshev expansions.

$h_{X}^{1}\left(s\right)=\sum_{\nu=0}^{2}a_{\nu}R_{\nu}\left(s^{2}\right),\ \
f_{X}\left(\alpha\right)=\sum_{\nu=0}^{4}b_{\nu}R_{\nu}\left(\alpha\right)$
(11)

For mTASK, we choose, $f_{X}\left(\alpha\rightarrow\infty\right)=-3.5$. We
found that the eight coefficients are $a_{0}\approx 0.924374$,
$a_{1}\approx-0.09276847$, $a_{2}\approx-0.017143$, $b_{0}\approx\ -0.639572$,
$b_{1}\approx-2.087488$, $b_{2}=-0.625$, $b_{3}\approx-0.162512$, and
$b_{4}\approx 0.014572$.

Figure S1:  Exchange enhancement factor of mTASK for more positive value of
$\mu_{\alpha}^{-}=0.0638768$ as a function of $\alpha$ for fixed value of s.

The mTASK enhancement factor is not smooth for $\mu_{\alpha}^{-}=0.0638768$,
and is not a monotonically decreasing function of $\alpha$. This condition is
associated with a convergence issue. We had convergence problems while testing
for band gaps of bulk solids.

Table S1: The width of flat nanoribbons used in the calculations (Å). The nanoribbons considered here have 20 and 12 MX2 formula units in super-cell respectively, for the 1H and 1T phases.  Phase | MX2 | Width (Å)
---|---|---
1H | MoS2 | 30.43
MoSe2 | 31.74
WS2 | 30.44
WSe2 | 31.74
1T | HfS2 | 36.57
HfSe2 | 37.88
ZrS2 | 37.05
ZrSe2 | 38.19
Figure S2: The band gap of the non-edge states of 1H TMDs nanoribbons (a)
MoS2. (b) MoSe2, (c) WS2, and (d) WSe2 at different bending curvatures with
different density functional approximations.  Figure S3:  Band structures of
flat 1T HfS2 nanoribbon in eVs from (a) PBE, (b) SCAN, (c) TASK, (d) mTASK,
and (e) HSE06 density functional approximations. The pattern of the band
structures is similar for all 1T nanoribbons that we have reported in this
paper.