Misfit Layer Compounds as Ultratunable Field Effect Transistors: From Charge Transfer Control to Emergent Superconductivity

Misfit layer compounds are heterostructures composed of rocksalt units stacked with few-layer transition metal dichalcogenides. They host Ising superconductivity, charge density waves, and good thermoelectricity. The design of misfits’ emergent properties is, however, hindered by the lack of a global understanding of the electronic transfer among the constituents. Here, by performing first-principles calculations, we unveil the mechanism controlling the charge transfer and demonstrate that rocksalt units are always donor and dichalcogenides acceptors. We show that misfits behave as a periodic arrangement of ultratunable field effect transistors where a charging as large as ≈6 × 1014 e– cm–2 can be reached and controlled efficiently by the La–Pb alloying in the rocksalt. Finally, we identify a strategy to design emergent superconductivity and demonstrate its applicability in (LaSe)1.27(SnSe2)2. Our work paves the way to the design synthesis of misfit compounds with tailored physical properties.

: Exact mismatch ratios a 2 /a 1 . Lattice parameters of the considered compounds are experimental values.  mismatch direction are slightly strained in order to build a commensurate misfit structure.
We obtain a TiSe 2 =3.6Å (tensile strain of ≈ 2%, exp. value a TiSe 2 =3.54Å) and a SnSe 2 =3.8Å (compressive strain of ≈ 0.3%, exp. value a SnSe 2 =3.81Å), respectively. In the other in-plane direction, we find b TiSe 2 =6.0191Å and b SnSe 2 =6.5818Å, respectively. Starting from this cell, a supercell is built according to the misfit proportions ( 5 × 1 and a 3 × 1 supercell for the cases of TiSe 2 /LaSe and SnSe 2 /LaSe misfits, respectively). Regarding the Q-layer rocksalt, we optimized LaSe with centered orthorhombic cell and slightly strained in-plane lattice parameter in order to obtain some commensurability with the considered TMD. A 3 × 1 LaSe supercell with in-plane lattice parameters a LaSe ≈b LaSe =6Å (tensile strain of ≈ 0.5%, exp. value a LaSe ≈b LaSe =5.97Å) is considered in order to match with TiSe 2 , and a 2 × 1 × 1 with lattice parameter a LaSe =5.7Å (compressive strain of ≈ 4%, exp. value a LaSe =5.97Å) and b LaSe =6.5818 Å (tensile strain of ≈ 10%, exp. value a LaSe =5.97Å) is considered in order to match with SnSe 2 . The two cells (one for the TMD and one for the rocksalt) are then appropriately assembled to build the slab system, which possesses a P1 symmetry in the most general case. This is composed of a Q-layer LaSe sandwiched between two TiSe 2 (or SnSe 2 ) single layers. The final (LaSe) 1.18 (TiSe 2 ) 2 cell is composed of 84 atoms, the mismatch ratio is a 2 /a 1 =(6Å)/(3.6Å) = 1.66 ≃ 5/3, leading to lattice parameter of the misfit supercell equal to a = 18Å. Instead, (LaSe) 1.27 (SnSe 2 ) 2 is composed of 52 atoms, the mismatch ratio is a 2 /a 1 =(5.7Å)/(3.8Å) = 1.57 ≃ 3/2, leading to lattice parameter of the misfit supercell equal to a = 11.4Å. In the tables S1 and S2 we report the optimized atomic positions expressed in crystalline coordinates. The crystal cell of (LaSe) 1.18 (TiSe 2 ) 2 is orthorhombic with a=18Å, b=6.0191Å, c=30Å. The crystal cell of (LaSe) 1.27 (SnSe 2 ) 2 is orthorhombic with a=11.4Å, b=6.581789Å, c=30Å.

II. Technical details
For what concerns the convergence parameters of DFT calculations of surface of the misfits, we employ a 2 × 8 × 1 a Monkhorst-Pack k-points grid and a Gaussian smearing of 0.025 Ry for Brillouin Zone (BZ) sampling. We use the generalized gradient approximation in the Perdew-Burke-Ernzerhof 1 parametrization for the exchange-correlation functional.
In the calculation of surface properties, as the interaction among transition metal dichalcogenides layers is missing and only covalent bonds among the rocksalt and the transition metal dichalcogenides are present, we did not consider any Van der Waals correction. For the case of the bulk superconducting properties of (LaSe) 1.27 (SnSe 2 ) 2 , where the interaction among adjacent SnSe 2 layers is present, we included the Van der Waals corrections Grimme-D3 2 in the cell relaxation in order to carefully reproduce the interlayer distance.
Spin orbit coupling (SOC) is included in all the electronic structure calculations in the main text. We verified that the inclusion of SOC marginally modifies the band structure of misfit layer compounds containing Bi and Pb rocksalts. Relativistic effects are negligible for the electronic properties of (La x Pb 1−x Se) 1.18 (TiSe 2 ) 2 family of compounds. For this reason, SOC is neglected in the calculations made for these compounds in the SI.
In case of TiSe 2 componunds, we employ PBE+U method described in Refs. 3,4 in order to take into account the strong correlation effects due to the localized d orbitals of Ti. The Hubbard correction is set to U= 3.25 eV, consistently with previous work on bulk TiSe 2 where a good agreement with ARPES spectra was demonstrated. 5 We consider the following pseudopotential configurations taken from the Vanderbilt 6 and  Figure S4: The calculation of V(z), which is the planar averaged electrostatic potential along the stacking axis z, is depicted in this figure. The negative peaks represents the spacial region in which the material is located, then we can determine the value of V VAC in the vacuum region of the plot (red line). Then, the work function is

III. Band Alignment Calculation
We use the band alignment method in order to make a computational screening of MLCs heterostructures from their elemental constituents. The work function is the amount of energy needed to extract an electron from a solid so that it exits with zero kinetic energy. Therefore, the difference of work functions of two compounds A and B is defined as:  IV. Band unfolding method applied to (La x Pb 1−x Se) 1.18 (TiSe 2 ) 2 We use unfolding method based of effective band structure (EBS) 12 as implemented in the BandsUP software. 13 We report the full calculation of (LaSe) 1.18 (TiSe 2 ) 2 with partial substi-

V. Doping-induced Superconductivity
As detailed in the text, in order to recover the (LaSe) 1.27 (SnSe 2 ) 2 surface and bulk behaviour, we used the 2D material FET setup 14 to precisely dope SnSe 2 as it is inside the misfit. A Coulomb long range interaction cutoff is placed at z cut = c/2 with c being the unit-cell size in the direction perpendicular to the 2D plane: c is set opportunely for single and double layer SnSe 2 at 16.14Å and 25.83Å, respectively. For SnSe 2 monolayer we use a single gate setup, placing charged plate modelling a single gate electrode at z bot = −0.25c with a charge of +0.7, equal and opposite to the one of the single layer SnSe 2 to ensure charge neutrality.
For SnSe 2 bilayer we use a double gate configuration, with the bilayer sandwiched between two charged plates at z bot = −0.266c and z top = +0.266c each with a charge of ρ=+(-)0.7, such that ρ tot = ρ 2L + ρ bot + ρ top = 0. For all systems, potential barriers V with a height of V H = 2.5 Ry are placed before the gates at z V = z bot + 0.1 (z V = z top − 0.1) in order to prevent the ions from moving too close to the gate electrodes.
The harmonic phonon frequencies are evaluated within density-functional perturbation theory 15 on a 7×7×1 phonon momentum grid (q-grid) and 21×21×1 electron-momentum grid (k-grid). The slightly imaginary acoustic phonons (order of 8 cm −1 ) observable in the phonon dispersion are an artifact linked to the acoustic sum rule breaking in linear response.
To avoid spurious effects, we did not impose the acoustic sum rule manually since it has been shown that in field-effect gated 2D system one of the three acoustic modes remains finite at Γ (see 14 ). The correct evaluation of the electron-phonon coupling properties of doped SnSe 2 requires a precise knowledge of electron-phonon matrix elements for very dense electron and phonon momentum grids. Since the direct calculation of electron-phonon matrix elements over a ultradense q-and k-point grids is very time consuming in linear response, we perform a Wannier interpolation of the electron-phonon coupling as described in Ref. 16 We used the Wannier90 17 code to obtain the Bloch to Wannier transformation. We use as starting guess of the Maximally Localized Wannier Function procedure three p-like orbitals at every chalcogen site and 5 d-like orbitals at any molybdenum site. Within this approach, the electron-phonon matrix elements are first calculated on a coarse 7×7×1 q grid and 21×21×1 k-grid, and then Wannier interpolated to 96×96×1 q-and k-grids in order to evaluate the electron-phonon coupling parameter λ and the isotropic Eliashberg function α 2 F (ω). We employ a Gaussian smearing of 0.001 Ry for k-and q-summations in λ.
The superconducting gap is evaluated by solving the Migdal-Eliashberg equations 18 in the Wannier basis over the imaginary frequency axis and then by performing analytic continuation to the real axis using N-point Padé approximants. 19 In order to perform the k-