The Effect of Janus Asymmetry on Thermal Transport in SnSSe

Several ternary “Janus” metal dichalcogenides such as {Mo,Zr,Pt}-SSe have emerged as candidates with significant potential for optoelectronic, piezoelectric, and thermoelectric applications. SnSSe, a natural option to explore as a thermoelectric given that its “parent” structures are SnS2 and SnSe2 has, however, only recently been shown to be mechanically stable. Here, we calculate the lattice thermal conductivities of the Janus SnSSe monolayer along with those of its parent dicalchogenides. The phonon frequencies of SnSSe are intermediate between those of SnSe2 and SnS2; however, its thermal conductivity is the lowest of the three and even lower than that of a random Sn[S0.5Se0.5]2 alloy. This can be attributed to the breakdown of inversion symmetry and manifests as a subtle effect beyond the reach of the relaxation-time approximation. Together with its low favorable power factor, its thermal conductivity confirms SnSSe as a good candidate for thermoelectric applications.


■ INTRODUCTION
Thermoelectric (TE) materials are reckoned to play a significant role in catering to the present global thrust for green and sustainable energy by enabling the conversion of waste heat into electricity. 1,2 Continued research in the field of thermoelectrics over the last decade has led to the discovery of several highly efficient bulk TE materials. 3−16 Theory shows that the key to optimizing the performance of TE materials lies in enhancing their TE figure of merit (zT), which is mainly achieved through alloying, doping, or nanostructuring. 17 zT is defined as S 2 σT/κ, where S is the Seebeck coefficient, σ is the electrical conductivity, κ is the the total thermal conductivity [a sum of electronic (κ e ) and lattice (κ l ) contributions], and T is the absolute temperature. Here, we use κ as a shorthand for κ l since we work with materials and regimes where thermal transport by charge carriers can be neglected as a first approximation. Due to the role of nanostructuring in impeding lattice thermal transport in bulk materials, the potential of lowdimensional materials has drawn a lot of attention in this field. Specifically, several two-dimensional (2D) materials like stanene, silicene, phosphorene, bismuthene, MoS 2 , MoSe 2 , WSe 2 , and SnSe have been reported to show low lattice thermal conductivities at room temperature. 18 −23 Among the 2D materials mentioned above, transition metal dichalcogenides (2D-TMDs) have attracted considerable interest due to their unique structure and semiconducting characteristics, which can lead to large zT. Shafique et al. 24 recently reported on the ultralow lattice thermal conductivity and high carrier mobility in monolayer SnS 2 and SnSe 2 using first-principles calculations. They showed that high thermoelectric performance could be realized in these simpler compounds that combine chemical stability, low toxicity, and earth abundance. A more complex class of 2D-TMD structures that are being intensely explored in this area are the so-called "Janus" monolayers where the two layers of chalcogen atoms belong to two different elements. The key to their usefulness lies in preserving the desirable characteristics of their "parent" binary structures but using the additional degrees of freedom to tune the material properties further. Applications of Janus monolayers have been proposed in electronics, optoelectronics, and piezoelectronics. Regarding thermal transport specifically, the lattice thermal conductivities of Janus MoSSe/ZrSSe/ PtSSe 2D-TMDs have been predicted to be much lower than those of the MoS 2 /ZrS 2 /PtS 2 monolayers. 25−27 This observation is in line with the traditional wisdom that a more complex unit cell leads to increased phonon scattering and lower thermal conductivity. 28−30 Guo et al. 31 recently reported a new mechanically stable Janus monolayer compound, SnSSe, which can be conceived as derived from SnS 2 and SnSe 2 and was shown to possess a very high power factor, making it a promising candidate for thermoelectric applications. However, the study focused on the (opto)electronic and piezoelectric properties of the SnSSe monolayer and there is, to our knowledge, no corresponding investigation of its thermal conductivity. Here, we undertake such a study, both as a necessary element to determine the suitability of the Janus SnSSe structure as a thermoelectric material and as a tool to assess the effect of symmetry in the family of monolayers formed by SnS 2 , SnSe 2 , and SnSSe. To be able to perform this comparison, we also calculate the phonon contribution to thermal transport in these compounds. We start with ab initio calculations but also employ a simpler predictive model developed by some of us 4,32 to enable a more direct comparison of the key ingredients for thermal transport between the three structures.
We find that along with its structure, the Janus SnSSe monolayer also inherits the main features of its phonon dispersion from its parent compounds. Nevertheless, the specific features of the phonon dispersion due to the broken symmetry, together with anharmonicity, make its lattice thermal conductivity lower than those of SnS 2 and SnSe 2 . The effect is not as drastic as in the Janus monolayers investigated in previous reports, and in fact, the order of the conductivities is incorrectly predicted by the commonly employed relaxation-time approximation. To understand it better, we compare the three chalcogenides in detail and also include the case of a Sn[S 0.5 Se 0.5 ] 2 alloy based on a SnS 2 structure where S sites are populated by S or Se with the same probability.

■ METHODS
The lattice thermal conductivity tensor is calculated by solving the linearized Boltzmann transport equation (BTE) for phonons as implemented in the almaBTE 33 software package. The final result is given by the following: 34 where α, β ∈ {x, y, z} are Cartesian coordinates and λ is the phonon mode comprising both the phonon branch index j and the wave vector q. v λ and C λ are the group velocity and the mode contribution to the volumetric specific heat, respectively. F λ is the solution of the linearized BTE written in the form In the above equation, τ λ is the lifetime of mode λ, and D λ accounts for the deviation of the population of a specific phonon mode from the prediction of the relaxation-time approximation. The total scattering rate τ λ is calculated considering the Taylor series of the potential energy up to the third-order term. It is written as a sum of the contributions from inelastic three-phonon (τ 3ph −1 ) and elastic two-phonon (τ 2ph −1 ) scattering processes as 34 In this study, two-phonon scattering is enabled by mass disorder at each atomic site, coming from either different isotopes of each element or from alloying. Both are treated in Tamura's approximation 35 where the contribution of a site to the scattering rate of a mode is proportional to the mass variance at that site. In real samples, boundaries 33,36 or defects 37−39 may have to be taken into consideration as well.
We also approximate the thermal conductivity using a model by Madsen et al. 40 that operates within the relaxation-time approximation. The relaxation time is approximated by the contribution from umklapp anharmonic scattering processes parameterized as 4,28 p T e 1 where p is an adjustable parameter given by

5)
V is the unit cell volume, s is the speed of sound, M is the average atomic mass, and θ D̃i s the Debye temperature defined as n k g g 5 3 ( )d In the above equation, n is the number of atoms per unit cell, g(ω) is the vibrational density of states, and γ̃2 is modeaveraged squared Gruneisen parameter given by, 29 The volume derivatives needed to obtain the Gruneisen parameter are calculated for a specific volume range using a finite-difference scheme.
We analyze thermal expansion in SnSSe in the framework of the quasi-harmonic approximation. For each given unit cell volume V, we approximate the free energy of the system by the sum of its potential energy when the atoms are frozen at their equilibrium positions and the free energy contribution from phonons obtained under the usual harmonic approximation for that volume, i.e., Here, the sums again run over all phonon modes, ω λ is the angular frequency of each mode, n λ (T) is its Bose−Einstein occupancy factor, and k B is the Boltzmann constant. The dependences on volume of E frozen , E phon , S phon , ω λ , and n λ have been omitted for brevity. Based on F(V, T), the equilibrium volume at each temperature, V(T), is defined as the value of V that minimizes F(V, T) at fixed T. In this context, fixing V is equivalent to fixing the unit cell area A, since the vacuum layer in the simulation box plays no role in contraction and expansion.
For SnSSe, we also obtain finite-temperature phonon spectra. Those are defined and calculated according to the prescription of the Quantum Self-Consistent Ab Initio Lattice Dynamics (QSCAILD) method. 41,42 A set of second-order interatomic force constants determines the real-space probability density function for atomic displacements in the canonical ensemble and also the effective harmonic force for each possible displacement. The QSCAILD method looks for the set of force constants that bring the harmonic estimate of the forces as close as possible to the DFT forces (in a leastsquares sense) when the atomic displacements are sampled from the aforementioned canonical distribution. In practice, this self-consistency condition is achieved iteratively: from an initial estimate of the force constants, the covariance matrix for the atomic displacements of all atoms in a supercell is obtained as The Journal of Physical Chemistry C pubs.acs.org/JPCC Article where i and j are atom indices, M i is the mass of atom i, and ϕ λ, i s is the element of the polarization vector (eigenvector) for mode λ corresponding to atom i and Cartesian direction α. Then, a sample of atomic displacements u is drawn from a Gaussian distribution with a probability density function proportional to , from the quantum canonical distribution corresponding to the effective harmonic potential built from those constants. Forces on the sampled configurations are calculated from first principles, and then a new set of interatomic force constants is obtained that provides the best fit to those forces in a least-squares sense from which a new covariance matrix is computed, and so on and so forth; the process is iterated until satisfactory convergence is achieved. The QSCAILD approach is able to explain finite-temperature effects like the stability of cubic oxide and fluoride perovskites 43 whose regular harmonic spectrum contains imaginary frequencies and can be shown to lead to the same finitetemperature spectrum as other approaches to the same problem that start from different principles, 42 like the stochastic self-consistent harmonic approximation (SSCHA). 44

■ COMPUTATIONAL DETAILS
All density functional theory (DFT) calculations are performed using the Vienna ab initio simulation package (VASP 45,46 ), the projector augmented wave (PAW) method, 47,48 and the Perdew−Burke−Ernzerhof (PBE 49 ) approximation to the exchange and correlation terms of the functional. The planewave cutoff is set at 336 eV, corresponding to the highest recommended cutoff among the PAW datasets employed in the calculation plus 30%.
The structures of monolayer SnS 2 , SnSe 2 , and SnSSe, oriented perpendicular to the OZ axis, are first relaxed using a 7 × 7 × 1 Monkhorst−Pack grid for integrations over the Brillouin zone (BZ) and a tolerance of 1 × 10 −7 eV for the total energy and the band energies, until the absolute value of all components of the force is lower than 1 × 10 −4 eV.
We then obtain the phonon spectrum of each structure using the Phonopy package 50 by calculating the harmonic force constants in a 7 × 7 × 1 supercell of each of the structures through a finite-difference scheme applying atomic displacements with an amplitude of 0.01 Å. Those force constants are then post-processed to enforce the continuous rotational symmetries of free space, 51 which has been shown not only to be necessary to obtain the right behavior of the out-of-plane phonon branches but also to have a large effect on the calculated thermal conductivity. The procedure consists in projecting the interatomic forces onto a space of internal coordinates, all of which are explicitly scalar and, in particular, rotationally invariant, and then transforming them back to the original Cartesian basis. The process separates interactions between pairs of atoms, which were originally aliased on top of one another because of the periodic boundary conditions; therefore, the post-processed constants do not fit in the original supercell, but require a 10 × 10 × 1 one instead. All DFT calculations for supercells are Γ-only and otherwise use all the same parameters as the initial relaxation.
Another set of supercell calculations is required in order to compute the anharmonic force constants leading to the amplitudes of the allowed three-phonon processes and ultimately to the ab initio phonon lifetimes. These are performed using the thirdorder.py script. 34 The supercell size is again 7 × 7 × 1, and after careful convergence tests, the cutoff for interactions is set at the 12th nearest neighbors, leading to more than 800 DFT runs for each structure after all symmetries have been harnessed. This is an unusually large cut-off, but related laminar 2D structures like SnSe have been shown 52 to require similar attention to this parameter.
In order to solve the phonon BTE and compute the thermal conductivity, phonon properties are sampled with almaBTE 33 on a dense 80 × 80 × 1 grid. That choice is again the result of a careful convergence study for all three structures.
We use a version of almaBTE modified so that each threephonon process allowed by the conservation of energy and momentum is tagged as normal, if all three phonons involved lie in the first Brillouin zone, or umklapp otherwise. This small change enables us to split the anharmonic scattering rates and all other quantities that can be expressed as sums over allowed processes in normal and umklapp contributions. The classification is implemented as part of the existing loop over pairs of phonon modes that the program uses to detect allowed three-phonon processes. An array holds a minimum-norm image of each wave vector in the regular grid used to sample the phonon modes; when a valid process is detected, we check whether the sum (for absorption processes) or difference (for emission processes) of the minimum-norm images of the wave vectors of the two phonons under consideration has a norm compatible with that of its own minimum-norm image.
To obtain the temperature-dependent equilibrium unit cell size of SnSSe in the quasi-harmonic approximation, we again use Phonopy. 53 We explore 10 different values of the volume, specifically 96, 97, ..., 105%, of the minimum-energy volume obtained from DFT. Expansions and contractions are isotropic in the in-plane directions, as required by the symmetry of the structure, and do not affect the vacuum layer along the direction perpendicular to the monolayer. For each volume, we fully relax the atomic positions and then obtain the phonon spectrum through a series of calculations completely analogous to those described above, including the supercell size, the DFT parameters, and considerations about rotational symmetry. To obtain a continuous F(V, T) from the discrete samples in V, Phonopy fits the results to the Vinet equation of state. 54 As part of the process of finding the self-consistent finitetemperature phonon dispersions of SnSSe in the QSCAILD framework, the hiPhive package 55 is used to generate displaced configurations according to the current best estimate of the effective harmonic potential and to fit a new set of secondorder force constants to the forces the system experiences in those configurations, as computed by DFT. These calculations also use 7 × 7 × 1 supercells. The QSCAILD iterative process is run for five iterations, in each of which 15 displaced configurations are generated. We choose recursive feature elimination 56 as the regularization strategy during the leastsquares fits to promote sparsity and avoid overfitting. The cutoff for the second-order force constants is set to 12 Å. Since hiPhive contains its own implementation of the sum rules that preserve the continuous rotational symmetry of the force constants, they do not require any post-processing.
To estimate the properties of a disordered Sn[S 0.5 Se 0.5 ] 2 alloy, we employ the virtual crystal approximation (VCA) functionality implemented in almaBTE, starting from the structures and force constants of the parent compounds. The VCA is very appealing because of its simplicity, but it has been shown 57 to not always deliver quantitative predictions because of structural distortions it cannot account for. We consider it suitable for this study because we only use the simplified alloy ■ RESULTS AND DISCUSSION Structural Properties. The dynamically stable 1T phases of monolayer tin dichalcogenides (SnX 2 , where X ∈ { S, Se }), belong to the P3̅ m1 space group and are shown in Figure 1a,b.
The optimized lattice parameters that we obtain are a = b = 3.69 Å for monolayer SnS 2 and a = b = 3.86 Å for monolayer SnSe 2 , in good agreement with the previously reported values. 24 This structure provides an octahedrally coordinated environment for Sn(IV), which is surrounded by chalcogen atoms, forming a hexagonal pattern. The vertical distances from the chalcogen atoms to the central plane of the monolayers, where the Sn atoms lie, are Δ S = 1.48 Å and Δ Se = 1.60 Å.
The Janus structure of the SnSSe monolayer belongs to the P3m1 space group [ Figure 1c] and is less symmetric than the SnS 2 and SnSe 2 monolayers due to the lack of inversion symmetry with respect to the Sn atoms in the central plane. The optimized lattice constant of the SnSSe monolayer is a = b = 3.78Å, larger than that of SnS 2 but smaller than that of SnSe 2 . In order to be able to provide the thermal conductivities in the usual units, we also need a layer thickness, which we obtain from the interlayer distance (c) in the corresponding layered bulk systems. We use the values 5.92 Å for SnS 2 , 6.38 Å for SnSe 2 , and 6.15 Å for SnSSe, which follow the same trend as the lattice constant. In fact, both a = b and c for SnSSe fall remarkably close to the arithmetic mean of their values for the parent structures. The distances from the S and Se atoms to the central Sn plane in the Janus structure are Δ S = 1.46 Å and Δ Se = 1.62 Å. Compared to the symmetric cases, Δ S is slightly shorter and Δ Se is slightly longer, partially compensating for the expanded (with respect to SnS 2 ) or contracted (with respect to SnSe 2 ) unit cell and keeping the Sn−S and Sn−Se bond lengths closer to their values in the symmetric monolayers.
Vibrational and Thermal Properties. Figure 2 shows a comparison of the harmonic properties of the three monolayer compounds, namely SnS 2 , SnSe 2 , and Janus SnSSe. No trace of imaginary frequencies is found in any of the dispersions or the density of states, confirming that all the structures are mechanically stable. In contrast, when we run the same calculations for other two possible tin-based Janus monolayers, SnSeTe and SnSTe, using the same DFT parameters and supercell sizes, both turn out to contain clear imaginary frequencies in the phonon spectra, pointing to their being mechanically unstable and making SnSSe singular in this regard. 58 A comparison between the phonon dispersions of the SnS 2 and SnSe 2 monolayers reveals a general trend toward lower frequencies in the latter. This is to be expected, since the atomic mass of Se (78.96) is more than double that of S (32.065). The presence of a common element in both structures (Sn, the heaviest of the three, with an atomic mass of 118.71) leads to different mass ratios in each of these two structures that manifest in another clear difference: a small gap between the acoustic and optical modes in SnS 2 that is absent from the spectrum of SnSe 2 . As another example of the role that the masses play, the highest frequencies computed for SnS 2 and SnSe 2 are 10.24 and 7.39 THz, respectively. Their ratio is 1.38, close to the value of M M / 1 . 5 6 S Se = . A secondorder contribution to the strong softening of the acoustic branches and the large downward shift in the optical branches comes from the unit cell sizes of each monolayer. The larger unit cell of SnSe 2 is typically associated with weaker, less stiff bonds.
On the most superficial level, the phonon branches of the Janus SnSSe, as shown in Figure 2, are simply intermediate between those of the parent symmetric monolayers. However, a closer look at points like K quickly shows this not to be the case: the lower symmetry of the SnSSe structure leads to broken degeneracies and large splits between branches that have the same frequency in the parent spectra. This is specifically the case with the two transverse acoustic bands and also with pairs of optical bands. Naturally, this lack of degeneracy at high-symmetry points like K affects the shape of the smooth bands throughout the BZ. The clearest example is perhaps the highest-lying optical band, which is contained in the same range of frequencies as the second highest in the symmetric SnS 2 structure but floats well above the rest of the spectrum in the Janus monolayer.
To look deeper into the contribution to the vibrational spectrum from each of the elements of the Janus structure, in Figure 3, we show the projected densities of states (PDOS) for each element in each structure. Note that in the symmetric structures the PDOS corresponds to a single atom of S or Se, to make it more easily comparable with the case of SnSSe.
Starting from the low frequencies, the density of states is dominated by the quadratic, non-degenerate ZA branch, as the only one whose contribution does not vanish as f → 0. The protagonists in this region are the more massive Sn atoms, and  The Journal of Physical Chemistry C pubs.acs.org/JPCC Article the differences between all three structures amount to a slight change in the frequencies of each mode, confirming what can be derived directly from the spectra in Figure 2. In fact, the PDOS from each chalcogen atom (S and Se) is very similar in the respective dichalcogenides and the Janus monolayer with a small change in frequency being explainable by the change in unit cell size. In other words, as far as the ZA branch is concerned, each of the sides behaves almost independently. This can be understood based on the fact that the polarization of the ZA branch is largely determined by the symmetry of the underlying structure. Moving on to higher frequencies, the behavior is far less universal because of the specificities outlined above, especially the overlap between acoustic and optical bands or the lack thereof. A feature can be identified in SnSSe that it inherits from SnS 2 , namely the very low participation of S in the lowfrequency region. A sort of reciprocal phenomenon is observed at the opposite end of the spectrum, i.e., the highest frequencies: there the role of Se is almost negligible, although as discussed above, the shape of the bands in SnSSe and SnS 2 is very different in this range.
To advance the discussion beyond the harmonic properties, the panels on the right in Figure 4 show the atom-projected mode Gruneisen parameters for each of the three structures as a descriptor of the contribution of each element to the anharmonicity of the crystal. It can be seen that at higher frequencies, S atoms have larger contributions than Sn atoms in SnS 2 . In Janus SnSSe, although the atomic contribution of S to the Gruneisen parameter dominates at higher energies, the contribution of Se atoms is larger for lower-energy phonon modes. The ZA branch stands out because of its negative Gruneisen parameter, a sign that its frequencies increase upon expansion of the crystal. The other branches show a similar behavior in all three compounds with respect to their frequency dependency but with the Sn contribution (the dominant one) at low frequencies being typically higher in SnS 2 and SnSSe.
Therefore, at this point, most factors point to SnSe 2 having a lower thermal conductivity than SnS 2 , specifically the reduced phonon group velocities caused by the phonon-softening and the increase in phase space for three-phonon processes introduced by the lower-lying optical branches. However, the relatively significant Gruneisen parameters of SnS 2 would hint in the other direction. As for SnSSe, with harmonic features sitting between those of its two parent structures, lower band degeneracy and high Gruneisen parameters, it could be suspected to have a lower thermal conductivity than either of the symmetric dichalcogenides.
Surprisingly, when calculated in the relaxation-time approximation, the thermal conductivity of SnSSe is found to be closer to but higher than that of SnSe 2 for all temperatures studied. This is shown in the bottom panel of Figure 5, comprising both ab initio results and those from the model explained in the previous section. They highlight how difficult it is to intuitively predict how much each of the factors mentioned above will influence the thermal conductivity. The model and the first-principles calculations agree well, and specifically at T = 300 K, they yield an SnSSe conductivity ∼1.5 times lower than that of SnS 2 and only ∼10% higher than that of SnSe 2 .
However, the full solution of the BTE (top panel of Figure  5) tells a different story. For the two symmetric 2D-TMDs, the difference between the relaxation-time results and those from the solution of the BTE is roughly a factor of two, but for the Janus monolayer, it is far less drastic. This difference is typically interpreted in terms of the relaxation-time approximation considering all processes as resistive, which is appropriate here in light of the good agreement with the semiempirical method. Hence, the results suggest that the breakdown of inversion symmetry opens up the possibility of more resistive processes and therefore decreases the relevance of normal ones for the final result.
To explore this hypothesis, in Figure 6 we plot the normal and umklapp contributions to the anharmonic scattering rates for all three systems under study with those corresponding to the lowest-lying branch (ZA) represented in a different color.  The Journal of Physical Chemistry C pubs.acs.org/JPCC Article As can be extrapolated from the case of graphene, studied in detail by Lindsay et al., 59 it is the ZA branch that is affected by the lifting of the selection rule arising from inversion symmetry. Comparing the cases of Janus SnSSe (center) and its parent structures (top and bottom) it is clear that both normal and umklapp processes are dramatically affected. However, while the normal contribution to scattering rates stays within the region where it can be considered comparatively modest, in general, the umklapp contribution to the scattering rates of the ZA branch in SnSSe is orders of magnitude higher than most total scattering rates for any of the three compounds. Scattering of phonons from the out-of-plane branch by umklapp processes is so intense in SnSSe that normal processes are comparatively irrelevant. With that branch almost out of the picture, the system behaves much more like a typical three-dimensional low-conductivity material where indeed the RTA and full BTE solutions tend to be closer to one another than in systems like graphene where the ZA branch makes an important contribution. The presence of strong anharmonicity calls into question the assumption, inherent in our treatment, that the phonon dispersions are not affected by temperature. We therefore examine that hypothesis by looking into two specific classes of possible effects in SnSSe: thermal expansion and renormalization of the phonon spectrum for a given unit cell size.
As indicated in the Methods section, we tackle the first category in the framework of the quasi-harmonic approximation, i.e., by calculating the phonon spectrum for a range of possible unit cell sizes, building a free energy F(V, T) that includes the DFT ground state and the vibrational contribution and taking its minimum at each given temperature as the equilibrium V(T). The results are shown in Figure 7, in terms of the volume itself and of the coefficient of thermal expansion where A representes the in-plane unit     60 the values of α are rather significant in the order of magnitude from 10 −5 to 10 −4 K −1 , larger than those of typical 3D materials.
Once the equilibrium unit cell size is known for each temperature of interest, we iteratively calculate the selfconsistent finite-temperature phonon spectrum in the QSCAILD framework. Figure 8 shows those spectra for two representative temperatures of 300 and 600 K, compared with the regular harmonic phonon spectrum of SnSSe that we use for all the calculations in this article. There are significant differences between the optical parts of the harmonic and finite-temperature spectra, especially for the higher temperature; in particular, the most energetic branch is softened (i.e., pushed to lower energies) and flattened, as one could expect from an expanded unit cell. However, the acoustic branches, which play a much bigger role in thermal transport, do not change so much, especially when taking into account that in systems where temperature renormalization of the phonon branches is really strong, it can change which phase is stable by pushing the vibrational frequencies from imaginary to real. 43 The only noticeable effect detected on closer inspection of Figure 8 is a softening of the ZA branch along the Γ → M line and even then only at 600 K. To put this in more quantitative terms, we calculate the small-grain thermal conductivity tensor κ SG αβ , defined as 61 Working under the relaxation-time approximation, if all the phonons in a material had the same mean free path Λ, the lattice thermal conductivity computed with eq 1 would be κ αβ = Λκ SG αβ . The small-grain thermal conductivity, therefore, summarizes the influence of all factors derived from the phonon dispersions (frequencies and group velocities) on the lattice thermal conductivity while leaving scattering aside. Using the regular harmonic phonon spectrum of SnSSe, the inplane component of the small-grain thermal conductivity is 0.781 W m −1 K −1 nm −1 at T = 300 K and 0.809 W m −1 K −1 nm −1 at T = 600 K; using the temperature-renormalized spectra computed at the volume determined by thermal expansion, the corresponding values are 0.760 W m −1 K −1 nm −1 at T = 300 K and 0.638 W m −1 K −1 nm −1 at T = 600 K. Those values represent decreases of 2.8 and 21% because of finite-temperature effects, respectively. We must therefore conclude that while phonon softening can be neglected at room temperature, the high-temperature κ values reported here are an upper bound, and finite-temperature effects can lead to even higher zT values for SnSSe in that regime.
In light of these results, SnSSe can be considered an ultralow thermal conductivity material in any context in which SnSe2 can merit that consideration. In fact, the predicted conductivity for Janus SnSSe likely makes it a suitable material for thermoelectric applications. In particular, with the values for the power factor reported earlier, 31 namely 65 μW cm −1 K −1 at 800 K, and with κ = 2.9 W m −1 K −1 as calculated in this report at the same temperature, the figure of merit for this monolayer is predicted to be around 1.8 at 800 K.
From the phenomenological model, we can extract that the low value of the thermal conductivity of SnSSe monolayer is due to its strong anharmonicity as evidenced by the high atomistic Gruneisen parameters. However, SnS 2 shares this feature and yet displays a significantly higher κ. Part of the explanation is that the Gruneisen parameter is descriptive only of the typical amplitudes of three-phonon processes. It conveys very little information about how easy it is for these processes to be allowed (or, more precisely, how big is the phase space for allowed three-phonon processes). Since SnS 2 has higher symmetry than SnSSe, such a phase space is relatively smaller than in the Janus structure, resulting in SnS 2 having a higher thermal conductivity. To show that phonon scattering is indeed more intense in the Janus monolayer as compared to SnS 2 , in the left-hand panels of Figure 4, we plot the total ab initio three-phonon scattering rates for each mode. Note that as discussed above, these scattering rates alone are not enough to explain why the thermal conductivity of SnSSe is also lower than that of SnSe 2 .
To explore whether a structure based on SnSSe can be built with an even lower thermal conductivity by exploiting disorder (and lead to a higher zT), in Figure 5, we also include an estimate of κ for a random Sn[S 0.5 Se 0.5 ] 2 alloy under the VCA. Given the way this approximation is constructed, such an alloy has the full symmetry of SnS 2 or SnSe 2 (each S/Se site being equivalent to any other) and the effect of disorder is felt through elastic phonon scattering instead. This is clearly reflected in the temperature dependence of the corresponding curve in the figure. Quite remarkably, the thermal conductivity of the alloy is higher than that of the Janus monolayer throughout the whole temperature range and more drastically at high temperatures. This is because the elastic scattering rates are temperature independent; therefore, the thermal conductivity of the alloy decays more slowly than the T −1 trend of the anharmonic rates. This shows how efficient the effect of removing inversion symmetry is and also that alloying does not appear to be a particularly attractive option for further reduction of the thermal conductivity. That said, other avenues for nanostructuring or doping remain open. In particular, the introduction of carefully chosen point defects can also cause a breakdown of symmetry leading to efficient mixing of phonon modes and to drastic decreases of κ, as has been shown in bulk systems. 38

■ CONCLUSIONS
In summary, our detailed study of the thermal conductivity of the Janus SnSSe monolayer points to low values that considered together with its favorable electric transport properties, can make it a promising thermoelectric material.
The key to it having a lower thermal conductivity than either of its parent structures is to be found in the lack of inversion symmetry with respect to the central plane of the structure and is captured by the Boltzmann transport equation but only partially by the relaxation-time approximation. The effect is subtler than in Mo-or W-based chalcogenides; however, since the conductivities of Sn-based chalcogenides are lower, SnSSe is a very interesting option. Our results should therefore encourage further experimental work toward the synthesis and thermal transport characterization of Janus SnSSe monolayers.