Anisotropic features in the electronic structure of the two-dimensional transition metal trichalcogenide TiS$_3$: electron doping and plasmons

TL;DR

This study investigates the electronic structure and plasmon behavior of TiS$_3$ monolayers, revealing anisotropic features, effects of electron doping, and potential electronic instabilities at specific doping levels using first-principles calculations.

Contribution

It provides a detailed analysis of anisotropic electronic properties and plasmon spectra in TiS$_3$, highlighting the impact of electron doping and spin-orbit coupling on its behavior.

Findings

Electron doping can induce incommensurate charge or spin modulations.

Highly anisotropic plasmons are supported with opposite anisotropy for electrons and holes.

Doped systems with smaller x remain 2D metals without electronic instabilities.

Abstract

Analysis of the band structure of TiS$_3$ single-layers suggests the possibility of changing their physical behaviour by injecting electron carriers. The anisotropy of the valence and conduction bands is explained in terms of their complex orbital composition. The nature of the Fermi surface and Lindhard response function for different doping concentrations is studied by means of first-principles DFT calculations. It is suggested that for electron doping levels $x$ (number of electrons per unit cell) $\sim$ 0.18-0.30$e^-$ the system could exhibit incommensurate charge or spin modulations which, however, would keep the metallic state whereas systems doped with smaller $x$ would be 2D metals without any electronic instability. The effect of spin-orbit coupling in the band dispersion is analysed. The DFT effective masses are used to study the plasmon spectrum from an effective low energy model. We find that this material supports highly anisotropic plasmons, with opposite anisotropy for the electron and hole bands.

Tables

Table 1. Table 1: Effective masses for the top of the valence band and bottom of the conduction band of the TiS 3 single-layer in units of the electron rest mass ( m e − subscript 𝑚 superscript 𝑒 m_{e^{-}} ).

		V.B.		C.B.
$\Gamma$-B		0.3317		1.8718
Y-$\Gamma$		0.8333		0.5253

Equations

$$\chi (\mathbf{q})=-\sum _{i,j}\sum _{\mathbf{k}}\frac{f_F({ϵ}_i(\mathbf{k}))-f_F({ϵ}_j(\mathbf{k}+\mathbf{q}))}{{ϵ}_i(\mathbf{k})-{ϵ}_j(\mathbf{k}+\mathbf{q})},$$

$$\epsilon (\mathbf{q},\omega )=1-V(\mathbf{q}){\Pi }^0(\mathbf{q},\omega )$$

$$V(\mathbf{q})=\frac{2\pi e^2}{{ϵ}_{B}q}$$

$${\Pi }^0(\mathbf{q},\omega )=\frac{g_s}{\mathcal{V}}\sum _{\mathbf{k}}\frac{f_F(\mathbf{k})-f_F(\mathbf{k}+\mathbf{q})}{\xi (\mathbf{k}+\mathbf{q})-\xi (\mathbf{k})-\hslash (\omega +i\eta )}$$

$${\xi }^{c,v}(\mathbf{k})=\frac{{\hslash }^2{k}_1^2}{2m_1^{c,v}}+\frac{{\hslash }^2{k}_2^2}{2m_2^{c,v}}-\mu$$

$${\Pi }^0(\mathbf{q},\omega )=\frac{g_{2D}}{2}\sum _{j=\pm 1}\left(\prod _{l=\pm 1}\sqrt{1+j\frac{\hslash (\omega +i\eta )}{\frac{{\hslash }^2{q}^{2}f(\theta )}{2m_1}}-2l\sqrt{\frac{\mu }{\frac{{\hslash }^2{q}^{2}f(\theta )}{2m_1}}}}-1\right)$$

$$f(\theta )={\cos }^2(\theta )+\frac{m_1}{m_2}{\sin }^2(\theta ).$$

$$\hslash {\omega }_{\mathrm{pl}}(\mathbf{q})=\sqrt{\frac{g_s{e}^2\mu }{{\epsilon }_b}\left[{\left(\frac{m_2}{m_1}\right)}^{1/2}{\cos }^2\theta +{\left(\frac{m_1}{m_2}\right)}^{1/2}{\sin }^2\theta \right]q}.$$

$$\hslash {\omega }_{\pm }=\frac{{\hslash }^2{q}^2}{2m_1}f(\theta )\pm 2\sqrt{\frac{{\hslash }^2{q}^2}{2m_1}\mu f(\theta )}.$$

Full text

I Introduction

Since the discovery of graphene in 2004,Castro Neto et al. (2009) there has been a huge improvement in the fabrication and manipulation of layered materials.Cao et al. (2015) Recently, the discovery that a monolayer of MoS2 changes its electronic properties with respect to the bulkMak et al. (2010) brought much expectation in the scientific community towards the transition metal chalcogenides (TMCs). The MX2 transition metal dichalcogenides (TMDCs) have been thoroughly studied both experimentallyCao et al. (2015); Novoselov et al. (2005); Wang et al. (2012); Mak et al. (2010) and theoretically.Roldán et al. (2014); Guo et al. (2014); Zhou et al. (2012); Rösner et al. (2014) Nowadays, TMCs with different chemical stoichiometries such as the transition metal trichalcogenides (TMTCs) MX3 (Fig. 1) are also being intensely studied. Dai et al. (2016); Island et al. (2017) Interestingly, TiS3 Island et al. (2015) has shown to have cleavage energies close to that of graphite, showing that similar methods can be used to fabricate TiS3 monolayers.Jin et al. (2015a)

The electronic structure of bulk group IVB TMTCs, MX3 (M= Ti, Zr, Hf; X= S, Se, Te), has been the subject of several theoretical studiesBullett (1979); Myron et al. (1981); Canadell et al. (1988, 1989); Stowe and Wagner (1998); Felser et al. (1998) but those concerning slabs of different thicknesses are more scarce. Jin et al. (2015a); Dai and Zeng (2015); Island et al. (2015); Li et al. (2015); Iyikanat et al. (2015); Kang and Wang (2016) These works have mostly dealt with the variation of the band gap in single-layers when changing the nature of both the transition metal and chalcogenide atoms or applying strain, the mechanical properties, and the role of vacancies. Some interesting tendencies have already been pointed out.Dai et al. (2016) For instance, among the semiconducting members of this family (i) the single-layer band gap diminishes when the size of the chalcogenide increases, (ii) the bulk indirect band gap may change to direct in the single-layers (TiS3) or keep its indirect nature (TiSe3), and (iii) the indirect band gap of ZrS3 and HfS3 can undergo an indirect-to-direct band gap transition with increasing tensile strain whereas the corresponding selenides keep their indirect band gap.

From a more fundamental viewpoint the potential of slabs with a few layers of group IVB TMTCs is still far from being fully explored. The conduction band of the semiconducting systems originates from the $d$ orbitals and, not surprisingly, the carrier mobility associated with electrons in the pristine TiS3 single-layer seems to be very anisotropic, with better mobility along the trigonal prismatic chains.Jin et al. (2015a) Therefore, it is expected that under electron doping the system may become an anisotropic conductor and thus be subject to the electronic instabilities of pseudo-one-dimensional (pseudo-1D) conductors. If this is the case, it would be interesting to see how much 1D would be these doped single-layers, i.e. would they exhibit metal to insulator or metal to metal transitions? It may be expected that the trigonal prismatic chains become more strongly coupled when the size of the chalcogen atom increases. The doped system could be then less 1D so that the physical behaviour of the doped single-layers could be modulated by varying the chalcogen atoms. At this point, let us remind that the Zr and Hf tellurium compounds are already metallic even without doping.Sambongi (1986) In fact, bulk ZrTe3 exhibits a charge density wave (CDW)Takahashi et al. (1984) at 63 K (whose origin is still debated) and bulk superconductivity under pressure ($T_{c}$ = 4 K at 5 GPa).Yomo et al. (2005)

The family of TMTCs with the TiS3 structure may thus provide a challenging series of compounds whose physical behavior under the form of single-layers or slabs with few layers can provide interesting and new phenomena besides the tunable optical properties. In the present work we explore such possibilities by means of first-principles density functional theory (DFT) calculations for TiS3 single-layers with different doping degrees. We start with a detailed analysis of the electronic band structure, in which the origin of the band anisotropy is discussed in terms of the orbital character of the valence and conduction bands. The results are used to study the screening properties and possible sources of instabilities of doped samples. We finally analyse the nature of the plasmon modes for $n$ and $p$ doped samples, finding a strong anisotropy in the spectrum, which is opposite for the valence and the conduction band modes.

II Computational details

First principles calculations were carried out using a numerical atomic orbitals approach to DFT, Hohenberg and Kohn (1964); Kohn and Sham (1965) which was developed for efficient calculations in large systems and implemented in the Siesta code.Soler et al. (2002); Artacho et al. (2008) We have used the generalized gradient approximation (GGA) and, in particular, the functional of Perdew, Burke and Ernzerhof. Perdew et al. (1996) Only the valence electrons are considered in the calculation, with the core being replaced by norm-conserving scalar relativistic pseudopotentials Troullier and Martins (1991) factorized in the Kleinman-Bylander form.Kleinman and Bylander (1982) The non-linear core-valence exchange-correlation scheme Louie et al. (1982) was used for all elements. We have used a split-valence double-$\zeta$ basis set including polarization functions.Artacho et al. (1999) The energy cutoff of the real space integration mesh was set to 500 Ry. To build the charge density (and, from this, obtain the DFT total energy and atomic forces), the Brillouin zone (BZ) for the bulk and monolayer was sampled with the Monkhorst-Pack schemeMonkhorst and Pack (1976) using grids of (30$\times$30$\times$30) and (30$\times$30$\times$1) k-points, respectively.

The Lindhard response function:

[TABLE]

where $f_{F}$ is the Fermi function, was obtained from the computed DFT values of the band eigenvalues $\epsilon_{i}({{\bm{k}}})$. The integral over k-points of the BZ was approximated by a direct summation over a regular grid of points. As the Lindhard function is more sensitive to the accuracy of the BZ integration than the total energy, especially in very anisotropic systems, and/or in the presence of hot spots in the band structure (e.g. saddle points with the corresponding van Hove singularity in the DOS), the k-points grid used for its calculation must be more dense. The calculations are done, nevertheless, using the eigenvalues obtained in the DFT calculation for the coarser grid, and interpolating their values in the denser grid, using a post-processing utility available within the Siesta package. In this work, for the calculation of the Lindhard function of the single layers, the BZ was sampled using a grid of (90$\times$90) k-points in the layer plane. The two lower conduction bands, which are those becoming partially filled for the electron doping levels considered in this work, were those taken into account in the calculations.

For the structural relaxations, we maintain the known symmetry of the crystal structure. For the bulk structures, as the weak interaction between layers is known to be severely underestimated by the GGA functionals, we have maintained the experimental value for the $c$ lattice parameter, while allowing the in-plane lattice parameters and the internal atomic coordinates to vary.

We also performed calculations using the HSE06Heyd et al. (2006) functional implemented in the Vienna ab initio simulation package (VASP)Kresse and Furthmüller (1996) in order to obtain the corrected values of the gap. The plane wave cutoff was set to 340 eV. The Brillouin zone for these calculations was sampled using a grid of ($7\times 10\times 4$) and ($7\times 10\times 1$) k-points for the bulk and monolayer structures, respectively.

III Results and discussion

III.1 Structural aspects

The MX3 layers of group IVB TMTCs (MX3) are built up from trigonal prismatic MX3 chains in such a way that two rectangular faces of an MX6 trigonal prism are capped by X atoms of the neighboring chains (see Fig. 1). Thus, every transition metal atom is coordinated to eight chalcogen atoms, X. One of the three X-X distances of the triangles is short and compatible with a single bond X-X so that for electron counting purposes this X-X pair must be considered as (X-X)2-. The system can thus be formulated as M4+(X-X)2-X2- and should be a semiconductor. Furuseth and coworkersFuruseth et al. (1975) found that these MX3 systems may have two structural variations, named type A and type B (see Fig. 1). In the type A structure the triangles and the two capping M-X distances are practically symmetric with respect to a plane containing the M and “isolated” X atoms of a chain. This plane is not really a symmetry plane of the structure because of the monoclinic crystal structure. In the type B structure any trace of this symmetrical situation is lost. Some compounds of this family were found to exhibit the type A structure (ZrSe3, ZrS3) whereas others the type B structure (TiS3, ZrTe3).Furuseth et al. (1975) However, later work concerning one type B structure (ZrTe3) showed that the structure was really of type A, casting doubts about the existence of two different structure types. Furuseth and Fjellvåg (1991) The inter and intrachain S-S contacts are considerably different in the two structures, so that the structural type can influence the location of the top of the valence band which is mostly chalcogen-based. Hence, we first decided to revisit the structure of bulk TiS3 by looking in detail to this point. Geometry optimizations with different starting geometries including those of the type B structure led always to a type A structure with geometrical parameters similar to those reported by Jin *et al. *Jin et al. (2015a) Thus, we confirm that TiS3 does not exhibit type B structure but the more symmetric type A and this is most likely the case for all TMTCs of this group. The calculations reported in the following sections use the bulk and single-layer DFT optimised structure.

III.2 Nature of valence and conduction bands

Shown in Fig. 2 are the calculated band structures and density of states (DOS) for the optimized structures of TiS3 bulk (a) and single-layer (b). In our calculations we find an indirect gap of 225 meV from $\Gamma$ to Z for the bulk, and a direct gap of 246 meV at $\Gamma$ for the single layer. The calculated effective masses of the top of the valence band and bottom of the conduction band for the TiS3 single-layer are shown in Table 1. We find that the mass anisotropy ratio, obtained from the effective masses along the $a$ and $b$ crystallographic directions $m_{a}/m_{b}$, is opposite for the valence and conduction bands, with respective values of $\sim 0.4$ and $\sim 3.56$. All these results are in good agreement with previous plane wave-based DFT calculations.Dai and Zeng (2015); Li et al. (2015); Iyikanat et al. (2015); Kang and Wang (2016) We have also checked that inclusion of spin-orbit coupling (SOC) does not lead to any substantial variation, as it can be seen in Fig. 2b (inset), where we compare the electronic band structure in the presence (red) and in the absence (blue) of SOC. Nevertheless it is important to notice that SOC leads to several band avoided crossings. This is specially important for the two branches of the valence band along the $\Gamma$-B direction, for which SOC leads to a sizeable splitting of $\sim 18$ meV.

Because of the well known underestimation of these band gaps when using the PBE functionalDai and Zeng (2015) we have recalculated the bulk and single-layer band structure using the HSE06 functional.Heyd et al. (2006) The calculated gaps are 1.09 and 1.12 eV for the bulk and monolayer structures, respectively. These values are in good agreement with the experimental bulk value, $\sim$ 0.9 eV.Grimmeiss et al. (1961) The band structures calculated with the two functionals are practically identical except for an almost rigid shift of the empty bands so that from now on our analysis is based on the PBE type calculations.

It is important for our discussion to clearly grasp the relationship between the nature of the valence and conduction bands and the structural parameters. Shown in Fig. 3 is a fatband analysis of the single-layer band structure. The orbitals are specified according to a local system of axes such that the x-axis coincides with the $a$-direction and the z-axis goes along the direction of the chains, i.e., the $b$-direction. The top of the valence band is mostly made of the S 3$p_{x}$ orbitals and a smaller contribution of the Ti 3$d_{xy}$ orbitals. In addition, the S 3$p_{x}$ contribution comes only from the inner sulfur atoms (i.e. those not forming S-S bonds; see Fig. 3). Thus, the valence band concentrates in the inner part of the layer and is mostly directed along the $a$-direction. Consequently, the top of the valence band should exhibit very weak interlayer interactions when the number of layers increases. This is indeed what we obtain in our DFT calculations for multilayer samples, shown in Fig. 4. Notice that the edge of the upper valence band appears basically degenerated for multilayer samples, whereas the second valence band suffers a noticeable splitting into well separated subbands, the number of which depends on the number of stacked layers. Since the main orbital contribution to this band comes from S 3$p_{z}$ of both inner and outer sulphur atoms, which play an essential role in inter-layer hopping, this leads to a splitting of this band in multilayer samples. In addition, the effective mass of the hole carriers should be considerably smaller along the $a$- than the $b$-direction of the layer, as shown in Table 1. The conduction band is strongly based on the Ti 3$d_{z^{2}}$ orbitals which point along the $b$-direction and this confers a strong anisotropy to the electron carriers. However, the conduction bands have also an important participation of the Ti 3$d_{x^{2}-y^{2}}$ and to a lesser extent from the S 3$p_{y}$ orbitals of the sulphur atoms of the outer part of the layer (i.e. those implicated in the S-S bonds; see Fig. 3). Such S 3$p_{y}$ contribution is the responsible for the downward shift of the valence band along the $\Gamma$-Z direction in the bulk and it further leads to splitting of the conduction band edges in multilayer samples, as it can be seen in Fig. 4. Of course, the effective mass of the electron carriers is anisotropic and smaller along the $b$-direction (see Table 1).

III.3 Electron doping

The fatband analysis of the TiS3 band structure is challenging in that it highlights the possibility of altering the physical behavior of TiS3 single-layers by injecting electron carriers in the conduction band through electric field gating. The dominance of Ti 3$d_{z^{2}}$ orbitals in the conduction band points out toward a strongly anisotropic metallic behavior (i.e. better conductivity along the $b$-axis) with associated inherent electronic instabilities under such conditions. However, one must bear in mind that the number of carriers which can be injected by gating has some limitations. Assuming, for instance, doping levels ($x$) similar to those attained in MoS2Ye et al. (2012) it would lead to $x$ values of up to $\sim$ 0.3$e^{-}$ per unit cell (i.e. per two TiS3 units). This would result with Fermi levels occurring between the bottom of the conduction band and the energy where the two lower conduction bands become degenerate at the B point (see Fig. 2b). In view of the non-negligible dispersion of the lower conduction band along the $a^{*}$-direction, at least for low and intermediate doping levels which will not lead to open Fermi surfaces, it is not clear what will be the detailed topology of the Fermi surface. In order to explore this issue we have calculated the Fermi surface and Lindhard response function for different doping levels. The atomic positions were re-optimized for every value of $x$ although it is not expected that structural effects can play an important role for these low doping levels.

The calculated band structure, Fermi surfaces and Lindhard response function for several electron doped TiS3 single-layer systems with doping levels $x$ = 0.1 to 0.3$e^{-}$ are shown in Figs. 5a-d. In these figures we superpose the calculated Fermi surface and the Lindhard response function. The portions of the Fermi surface originating from the lowest and second lowest conduction bands are shown in red and green, respectively. As can be seen from the almost identical calculated band structures, electron doping with these $x$ values practically does not alter the band structure. If the conduction band was a perfect 1D system the Fermi surface would consist of two parallel lines perpendicular to the $b$-direction at $\pm$ 0.025$b$ ($x$ = 0.1$e^{-}$) and 0.05$b$* ($x$ = 0.2$e^{-}$). However, this is very far from the computed results. From $x$ = 0 to $x$ = 0.1$e^{-}$ the Fermi surface is closed (Fig. 5a), meaning that there are non-negligible inter-chain interactions. The Fermi surface is an ellipse with the long axis along $a$* essentially because the slope of the band is larger along $b$. The Lindhard function has a somewhat broad region with larger values for wave vectors corresponding to the nesting of the locally flat borders of the ellipse at 0.22$a$. However, these values are relatively small and it is not expected that they will lead to any instability. Thus, for this doping interval the system should be a 2D metal (except for very small values of $x$ where the system may exhibit activated conductivity because of the low density of carriers and the potential due to random sulfur vacancies).

When $x$ increases, the Fermi surface is still closed but clear 1D features already appear. The shape of the ellipse changes in such a way that long flat portions occur (see Fig. 5b for $x$ = 0.2$e^{-}$). This is due to the fact that the density of carriers has increased and the Fermi level is now reaching a region where the slope along the chains direction ($\Gamma$-Y) strongly dominates over the slope along the inter-chain direction ($\Gamma$-B), i.e. a region where the electron gas already exhibits a pseudo-1D behaviour. However, the Fermi surface is not made of two cosine-like lines as it could be expected for such a case. This is due to the fact that the coupling between chains is large enough to outweigh the effect of the carrier increase and the Fermi surface is ultimately closed even if the pseudo-1D character inherent to the valence band already shows up. The flat sections are well nested and lead to maxima of the Lindhard function which for $x$ = 0.2$e^{-}$ occur at $\pm$ 0.5$a$$\pm$0.085$b$ (see Fig. 5b). Thus, for values of $x$ approaching 0.2$e^{-}$, electronic instabilities partially destroying the Fermi surface originating from these nesting vectors can occur.

Following the previous reasoning one could expect that for larger values of $x$ the pseudo-1D character of the system would clearly appear as an open and well nested Fermi surface made of two cosine-like lines. Slightly above 0.20$e^{-}$ the Fermi surface indeed becomes open but the shape is not the expected one. That corresponding to 0.25$e^{-}$ is shown in Fig. 5c. Note that for doping levels just above $x$ = 0.22$e^{-}$ the Fermi level cuts also the second conduction band around the B point. For very low band fillings this band exhibits similar curvature along the B-$\Gamma$ and the B-A directions so that a closed component around the B point emerges (the green lines of the Fermi surface). The shape of the portion originating from the first band strongly reminds that of Fig. 5b but around $\pm$ 0.38$a$* exhibits two minima. In fact, there are these regions between $\pm$ 0.38$a$* and B which are responsible for the maxima of the Lindhard function, which exhibits two arcs with maxima in the center, at $\pm$0.05$b$. The Fermi surface keeps this general shape until $x$ values of 0.3$e^{-}$ (Fig. 5d) and above. Except for the increase in the separation of the warped lines, the only difference is that when the area of the closed portion around B increases, the circle becomes an ellipse. The Lindhard function of systems with doping levels between 0.22$e^{-}$ and 0.3$e^{-}$ is the same although the separation of the arcs increases with the doping (they occur at $\pm$ 0.02$b$* for $x$ = 0.22$e^{-}$ but $\pm$0.085$b$* for $x$ = 0.3$e^{-}$). This nesting vector occurs in the direction of the trigonal prismatic chains, $b$*, but the modulus is smaller than it would be if the system was really a pseudo-1D system (even if taking into account the contribution of the small closed part around B; see Section III.4). These results suggest that for values of $x$ between 0.22 and 0.3$e^{-}$ there could occur an instability associated with this nesting vector which is very different from that discussed for lower carrier concentrations. This nesting vector would be responsible for the destruction of the flat portion of the Fermi surface around the border of the Brillouin zone. Note than when the closed portion of the Fermi surface around B appears, the Lindhard function has an additional contribution around the center of the arcs due to partial nesting of this closed part (for similar reasons as those discussed for Fig. 5a). To summarize, we believe that systems doped with $x$ values around 0.18-0.2 $e^{-}$ and between 0.22-0.3$e^{-}$ could exhibit two different types of incommensurate modulations of the charge density wave (CDW) or spin density wave (SDW) type which, however, would keep the metallic state whereas systems doped with less than around 0.18-0.2 $e^{-}$ would be 2D metals without any electronic instability.

III.4 Inter-chain coupling and the conduction band

As discussed above, the Fermi surfaces for the doping levels considered here are different from those intuitively expected meaning that the interchain interaction is not a trivial one. As shown in Fig. 2b, the two lowest empty bands exhibit the typical behavior along $b$* for a highly dispersive band based on the Ti 3$d_{z^{2}}$ orbitals. The separation of the two bands at $\Gamma$ is around 0.3 eV which is non negligible but not exceptionally strong. Looking at the $\Gamma$–B direction (see Fig. 3) one notices that the two lowest bands essentially result from a strong mixing of the 3$d_{z^{2}}$ and 3$d_{x^{2}-y^{2}}$ orbitals, leading to a moderately dispersive and a dispersionless pair of bands, respectively. However, in the close vicinity of B there is an avoided crossing with a strongly dispersive 3$d_{xy}$ based pair of bands. Thus, it appears that the nature of the lowest conduction bands along the interchain direction ($\Gamma$–B) changes significantly, acquiring 3$d_{xy}$ character in the last part of the line and this is at the origin of the unexpected behavior.

As it is well known,Hoffmann et al. (1980); Canadell et al. (1989) the fact that the sulfur triangles are far from equilateral, does not only change the electron counting but also induces an important rehybridization of two of the three low-lying levels of the transition metal, 3$d_{z^{2}}$ and 3$d_{x^{2}-y^{2}}$. The reason is that these levels try to minimize the antibonding interactions with the S 3$p$ levels. When the Ti atoms become eight coordinated because of the two additional Ti-S bonds between chains, one of these two levels becomes well separated from the other Ti-based levels and spreads out to lead to the two lowest bands of the undoped system. Depending on the degree of mixing between the two orbitals, the resulting one can be described as 3$d_{z^{2}-y^{2}}$ (for moderate mixing) or 3$d_{y^{2}}$ (for a strong mixing). Analysis of the fatbands and the detailed crystal orbitals clarify the shape of the band structure in that zone. The two lowest empty crystal orbitals of the undoped system at $\Gamma$ are schematically shown in Fig. 6a and b. The lowest level is based on an out-of-phase combination of 3$d_{z^{2}-y^{2}}$ orbitals. These orbitals make moderately stabilizing lateral interactions with the 3$p_{y}$ orbital of the non-bonded S atoms of the two adjacent cells, which are at the same level (the interactions marked in orange in Fig. 6). Half of these interactions occur within the cell and half occur between cells. When the phase changes from $\Gamma$ to B, the intercell interactions become antibonding, so that globally this band is non bonding between chains and consequently the band goes up in energy. In principle, the interaction between the 3$p_{y}$ orbital of the inner S atom of one chain and those of the bonded outer S atoms of the two neighboring chains could also contribute to the energy raising since the intracell overlap is positive both in $\Gamma$ and B but the intercell one changes to antibonding in B. However, the distance between the two sulfur atoms is large because the S atoms are located in planes differing by half-$b$ and the effect is small. In fact, there is a shorter S–S interchain contact between the bonded S atoms in second-neighboring chains along $a$. This $\pi$-type interaction could also contribute to the energy raising since they are in-phase at $\Gamma$ but out-of-phase at B. However, this contribution turns out to be extremely weak as shown by the fact that it also occurs for the next band (see Fig. 6b) which is, nevertheless, very flat along $\Gamma$-B. Thus, for around two thirds of the $\Gamma$ to B traject, the lowest conduction band behaves as resulting from a series of slightly coupled strongly 1D systems located in the trigonal prismatic chains and the closed nature of the Fermi surface is simply due to the joint effect of the low density of carriers and the non-negligible interaction between chains of 3$d_{z^{2}-y^{2}}$ orbitals.

The next empty crystal orbital is based on an in-phase combination of 3$d_{y^{2}}$ orbitals. An important difference with the previous band is that now the 3$p_{y}$ orbital of the non-bonded S atoms practically does not mix into this crystal orbital so that the different chains become uncoupled and hence, the completely flat nature of this band. Finally, near the B point a strongly descending pair of bands based on the Ti 3$d_{xy}$ orbitals (see Fig. 6c), which make strong interactions with the 3$p_{x}$ orbitals of the adjacent non bonded S atoms near $\Gamma$ but not around B, undergoes an avoided crossing with the lower pair of bands. As a result, Ti 3$d_{xy}$ character is built up in the two lower bands near the region of the B point stabilizing the two lower bands. More importantly, the antibonding interaction between the bonded S atoms of two successive chains along $a$ at $\Gamma$ (see Fig. 6c) becomes bonding around B. Thus, the avoided crossing brings about an effective decrease of the inter-chain coupling in the region around B. The abnormal shape of these bands in this region is thus the result of the competition between destabilizing interactions between the Ti and capping S atoms (see Fig. 6a) and stabilizing direct S–S interactions of two S–S bonded pairs in adjacent chains along the interlayer direction. It is important to note that the slope of the strongly descending pair of bands will increase when the chalcogen atom is bulkier, i.e. for TiSe3 vs. TiS3, or under some contraction along the $a$ axis. In these cases the bottom of the conduction band will most likely occur around B. This will have a very strong influence both on the Fermi surface of the electron doped system and the optical properties.

III.5 Plasmons in Single-Layer ${\rm{\bf TiS_{3}}}$

One of the most promising applications of the new families of 2D materials is their potential for optoelectronics and nanoplasmonics.Grigorenko et al. (2012); Low et al. (2017) Plasmons are collective density oscillations of an electron liquid that occur in many metals and semiconductors. Because of their single or few layer structure, these collective modes in 2D materials are highly confined, and can be tuned and manipulated by external gate or chemical doping. Plasmons have been studied in different 2D materials like graphene,Shung (1986); Wunsch et al. (2006); Hwang and Sarma (2007); Stauber (2014); Roldán and Brey (2013) silicene,Tabert and Nicol (2014) TMDs,Scholz et al. (2013); Groenewald et al. (2016) or black phosphorus.Low et al. (2014) In this section we study the plasmon spectrum of doped single-layer TiS3. The zeros of the dynamical dielectric function $\varepsilon({\bm{q}},\omega)$ yield the excitation spectrum of the plasmon modes of the electron liquid. The loss function, defined as ${\cal L}({\bm{q}},\omega)=-\Im[1/\varepsilon({\bm{q}},\omega)]$, quantifies the spectral weight of the plasmon mode. In the limit of zero damping, ${\cal L}({\bm{q}},\omega)$ consists in a delta peak. The dielectric function can be calculated, within the random phase approximation (RPA) asAndo et al. (1982); Giuliani and Vignale (2005)

[TABLE]

where $V({\bm{q}})$ is the Coulomb interaction

[TABLE]

where $\epsilon_{B}$ is the background dielectric constant and $\Pi^{0}({\bf q},\omega)$ is the polarization function

[TABLE]

where $g_{s}=2$ is the spin degeneracy, $\cal V$ is the system size, $f_{F}({{\bm{k}}})=\{\exp[\xi({\bm{k}})/k_{B}T]+1\}^{-1}$ is the Fermi-Dirac distribution function, and $\eta$ is a phenomenological broadening. Here we are interested on low energy intra-band plasmon modes. Therefore, we neglect inter-band transitions in the polarization (4) and consider the energy dispersion within the effective mass approximation:

[TABLE]

where $\mu$ is the chemical potential, the superscript $c(v)$ indicates conduction (valence) band, and the corresponding effective masses have been obtained from the DFT band structure and given in Table 1.

As we have seen in the previous sections, single layer TiS3 is a highly anisotropic material, with energy bands dispersing very differently in the two crystallographic directions. Therefore, we follow the scheme used by Low et al. to study plasmons in black phosphorus,Low et al. (2014) which is another anisotropic 2D crystal that consists in P atoms arranged in a puckered honeycomb lattice. This procedure has been later generalized to study plasmons in rotated double layer systems,Rodin and Castro Neto (2015) Coulomb drag in anisotropic van der Waals heterostructuresSaberi-Pouya et al. (2016) and 2D systems with merging Dirac points.Pyatkovskiy and Chakraborty (2016) Denoting ${\bm{q}}=q(\cos\theta,\sin\theta)$, we can write the polarization function as

[TABLE]

where $g_{\rm 2D}=m_{d}/\pi\hbar^{2}$ is the 2D density of states (DOS), in terms of the DOS mass $m_{d}=\sqrt{m_{1}m_{2}}$, and

[TABLE]

Our results are shown in Fig. 7, where we plot the loss function ${\cal L}({\bm{q}},\omega)$ for the valence and conduction bands along the two crystallographic directions. As in the case of black phosphorus,Low et al. (2014); Jin et al. (2015b) we observe a clear anisotropy in the plasmon mode. However, contrary to black phosphorus, the plasmon anisotropy in TiS3 is opposite for the valence and conduction bands. Indeed, one can clearly observe that for electron (hole) doping, the plasmon dispersion is faster (slower) for wave-vectors ${\bm{q}}\parallel{\hat{\bm{q}}}_{2}$ (${\bm{q}}\parallel{\hat{\bm{q}}}_{1}$), as seen by the slope of the corresponding modes in Fig. 7. The anisotropy is more clearly seen in Fig. 8 where we show a cross section of the loss function in the $q_{1}-q_{2}$ plane for $\hbar\omega/\mu=0.2$. The origin of the opposite anisotropy for the valence and conduction bands originate from the opposite rate between the effective masses $m_{1}/m_{2}$ for each band (see Table 1). This can be well understood from the different orbital character of the valence and conduction band edges, as explained in Sec. III.2. In brief, the top of the valence band is mostly made of the 3$p_{x}$ orbitals of S coming only from the inner sulfur atoms. This leads to a valence band dispersing mostly along the $q_{1}-$direction ($m^{v}_{1}<m^{v}_{2}$). On the other hand, the conduction band is mainly based on the $d_{z^{2}}$ orbitals of Ti directed along the $q_{2}-$direction. Therefore the effective masses of the electron carriers are anisotropic with $m^{c}_{1}>m^{c}_{2}$. This opposite rate between the effective masses for the valence and conduction bands leads at the end of the day to the different behaviour of the collective plasmon modes that appear after electron or hole doping.

The plasmon dispersion at long wavelengths presents the standard $\sim\sqrt{q}$ behaviour typical for 2D electron gases,Ando et al. (1982) with a correction that accounts for the anisotropyLow et al. (2014); Rodin and Castro Neto (2015); Saberi-Pouya et al. (2016); Pyatkovskiy and Chakraborty (2016) and it can be well approximated by the analytical expression

[TABLE]

The dispersion obtained from the approximation (8) is shown by dashed grey lines in Figs. 7 and 8. One should notice that the the plasmons are coherent outside the particle-hole continuum, which is defined as the region where $\Im[\Pi^{0}({\bm{q}},\omega)]\neq 0$. The boundaries of the continuum are defined by

[TABLE]

When the plasmon touches the continuum threshold, the mode starts to be damped by decaying into electron-hole pairs. One of the consequences of the anisotropy of the spectrum is that the damping start to occur at different $\omega-{\bm{q}}$ rates depending on the direction of propagation. We finally notice that the plasmon frequency presents the standard $\sqrt{\mu}$ dependence with the chemical potential (see Fig. 9), as expected for a plasmon in a 2D crystal, like in the well known case of graphene.Shung (1986); Wunsch et al. (2006); Hwang and Sarma (2007) However, contrary to graphene for which $\mu\propto n^{1/2}$ where $n$ is the carrier density, in TiS3 we have $\mu\propto n$, leading to the expected scaling relation $\omega_{pl}\propto n^{1/2}$, as in standard 2D electron gases,Giuliani and Vignale (2005) while $\omega_{pl}\propto n^{1/4}$ in graphene due its low energy linear dispersion relation.

IV Conclusions

In summary, we have used first-principles methods to study the electronic properties of single layer TiS3. Careful analysis of the DFT band structure points out the possibility of changing their physical behaviour by injecting electron carriers. A study of the Fermi surface and Lindhard response function for different electron doping levels shows that for electron doping levels $x$ $\sim 0.18-0.30e^{-}$ ($x$ being the number of electrons per unit cell) the system could exhibit incommensurate charge or spin modulations which, however, would keep the metallic state. In contrast, systems doped with smaller $x$ would be 2D metals without any electronic instability. The inter-chain coupling between the different trigonal prismatic TiS3 chains influences the bottom states of the conduction bands throughout a complex mixing of Ti 3$d_{z^{2}}$, 3$d_{xy}$ and 3$d_{x^{2}-y^{2}}$ orbitals as well as 3$p_{x}$ orbitals of the bonded S-S pairs. This leads to an unexpected shape of the Fermi surface and Lindhard response function for a large part of the carrier concentrations examined. We further find that, while spin-orbit interaction does not produce any significant modification of the valence and conduction band states, it leads to several avoided band crossings in the spectrum. Finally, we discuss the main features of the plasmon spectrum in doped single layer TiS3 samples. By using a low energy effective mass theory model, we find that plasmons in TiS3 are highly anisotropic. Interestingly, the anisotropy is opposite for the electron and hole branches, which can be understood from the different orbital nature of the valence and conduction bands. We find that the mass anisotropy ratio is $\sim 0.4$ for the valence band and $\sim 3.5$ for the conduction band. This anisotropy, which is opposite for the valence and conduction bands, might be compared to the same directional anisotropy of black phosphorus, with ratios of $\sim 6.2$ and $\sim 5.1$ for the valence and conduction bands respectively.Yuan et al. (2015) Another measure of the anisotropy can be obtained from the mobility ratios. In this case, it has been experimentally shown that few layers TiS3 present a highly temperature dependent ratio of $\sim 2.3$ to $7.6$ from room temperature to 25 K respectively,Island et al. (2016) while the anisotropy in the mobility of black phosphorus is of $\sim 1.5$ and only weakly dependent on temperature.Xia et al. (2014)

Our results show that TiS3 nanofilms are promising platforms for future optoelectronics and nanoelectronics applications, including field effect transistors or infrared photodetectors. In particular, the strong anisotropy of the electrical and optical properties of this material can lead to novel functionalities for devices based on TiS3 like high-performance transistors built along the light effective mass direction, or the directional focus of photocarriers, as demonstrated in black phosphorus,Quereda et al. (2016) but with the additional feasibility in TiS3 to reverse the main direction of propagation by switching from $n-$ to $p-$doping and viceversa.

V Acknowledgments

We thank A. Castellanos-Gomez for useful discussions and R. Robles for assistance with the calculations. This work has received funding from the European Union’s Seventh Framework Programme (FP7/2007-2013) through the ERC Advanced Grant NOVGRAPHENE (GA No. 290846), European Commission under the Graphene Flagship, contract CNECTICT-604391, the Spanish MINECO (through Grants No. FIS2015-64886-C5-3-P, FIS2015-64886-C5-4-P, FIS2014-58445-JIN and the Severo Ochoa Centers of Excellence Program under Grants SEV-2013-0295 and SEV-2015-0496), and the Generalitat de Catalunya (2014SGR301), as well as from the Comunidad Autónoma de Madrid (CAM) MAD2D-CM Program (S2013/MIT-3007).