Room temperature multiferroicity in a transition metal dichalcogenide

The coexistence of multiple ferroic orders, i.e. multiferroicity, is a scarce property to be found in materials. Historically, this state has been found mainly in 3-dimensional complex oxides, but so far this state has still been elusive for the most widely studied and characterized family of 2-dimensional compounds, the transition metal dichalcogenides. In this study we report the experimental realization of multiferroic states in this family of materials, at room temperature, in bulk single crystals of Te-doped WSe 2 . We observe the coexistence of ferromagnetism and ferroelectricity, evidenced in the presence of magnetization and piezoresponse force microscopy hysteresis loops. These findings open the possibility of widening the use and study of van der Waals-based multifunctional devices for new nanoelectronics and spintronics applications.


Introduction
Multiferroics are materials that exhibit different coexisting ferroic orders such as ferroelectricity, ferromagnetism, or ferroelasticity.Due to the coupling among the different degrees of freedom leading to these ordered states, the order parameters of one state can be controlled by tuning parameters different from their conjugate variable.For example, magnetoelectric coupling allows the control of magnetization via an electric field, or conversely, the control of electric polarization via magnetic field, for the case of ferromagnetic/ferroelectric heterostructures and multiferroics.Multiple applications in spintronics, data storage, actuators, among others (1-4) have been envisaged and already realized.The search for new and better multiferroic materials is arduous: ferroelectrics tend to be insulators as they need to preserve electric polarization (free charges in metals screen this effect) while ferromagnets are metals in their majority (5).In this sense, a guiding strategy that has been used to attain profitable multiferroicity is chemical doping of either a ferroelectric host or conversely, a ferromagnetic host.For example, replacement of the transition-metal ion in ferroelectric LiNbO3 or ZnO leads to magnetic order.On the other hand, chemical doping in complex oxides multiferroics has been effective in optimizing their properties (6,7).
Historically, these phenomena have been found mainly in 3-dimensional complex oxides and perovskites, such as Cr2O3, YMnO3, BiFeO3, among others, or in heterostructures of ferroelectric/ferromagnetic thin films (8)(9)(10).The 3D character of the crystal structure of these compounds poses challenges in their incorporation into nanoscale multiferroic devices.Their nanostructuration requires advanced fabrication techniques, as well as induces physical effects inherent to quantum confinement that typically weakens the multiferroic states (11).In this sense, the use of multiferroics with an intrinsic 2D-crystalline structure is highly desirable for obtaining affordable and robust multiferroic nanodevices.
Recent advancements in the realization and study of intrinsically 2D-multiferroics have been done (12,13), for example in the transition metal halide compounds such as NiI2 (14) and CrI3 (15), and in doped In2Se3 (16).However, multiferroicity has been elusive, until now, for the most widely studied and characterized family of 2D compounds, the transition metal dichalcogenides (TMD's), in spite of numerous theoretical works in this respect (17,18).In recent years, for this family of materials, individual ferroic orders have been experimentally observed at room temperature (19)(20)(21)(22)(23)(24).Ferromagnetism has been attained in MoS2 and WS2 due to defects (21,23), and in WSe2 via chemical doping in the transition metal site (24), among other examples.Ferroelectricity has been reported in bulk-WTe2 (20), while for an odd number of layers, WSe2 and MoS2 are reported to be piezoelectric (25,26).Despite the presence of these individual ferroic orders, there is no experimental evidence for simultaneous ferromagnetism and ferroelectricity in TMDs to date.
In this study we report the experimental realization of multiferroic states in TMDs, specifically in Te-doped WSe2.The multiferroic state is observed at room temperature and in bulk single crystals, for Te doping above 7%.The multiferroic states are revealed in the observation of ferromagnetic hysteresis loops in the magnetization, and ferroelectric hysteresis loops in the phase and amplitude of piezoresponse force microscopy (PFM) measurements.In addition, the highest Te-doping samples studied (15%) show very large effective piezoelectric coefficients, comparable to the commercially used 3Dpiezoelectric LiNbO3 (27), and one order of magnitude higher than in WSe2 and other pure TMD's nanolayers (25,26,28).In addition, we present density functional theory (DFT) calculations which provide insights about the possible mechanisms leading to the multiferroic state.These findings open the possibility of widening the use and study of van der Waals (vdW)-based multifunctional devices, and provide an affordable platform for the study of fundamental physics of ferroic orders on low-dimensional compounds.

Raman spectroscopy and X-rays techniques
Figure 1(a) shows a picture of representative single crystals of Te-doped WSe2, with typical lateral sizes around 1mm, and a hexagonal morphology.A level of chalcogen vacancies is naturally developed during the crystal growth process, and crystals with a general stoichiometry W(Se1-xTex)2(1-δ) are produced.x represents the fraction of Se atoms that were replaced by Te atoms, and δ is the fraction of chalcogen sites that are not occupied by neither a Se nor a Te.The values of x and δ for the different crystals studied were determined through X-ray fluorescence spectroscopy (see methods, and table S1 for values).Samples with four different Te-doping values, up to a maximum of x=15%, plus undoped samples (x=0) with three different levels of chalcogen (Se) vacancies were characterized in this study.Two different types of X-ray diffraction measurements (powder diffraction, presented in figs.1(d-f) and fig.S2; and single crystal diffraction, presented in table S2 and crystallographic information file -CIFavailable in supplementary files) reveal that the crystalline structure of all the pure and Tedoped compounds is the 2H-structure (space group P63/mmc, # 194) (29).The most intense diffraction peaks in fig.1(d) come from the {0 0 1} family of planes, given that the a-b plane of the crystals tends to align parallel to the sample holder (see methods).Close inspection to the {0 0 1} diffraction peaks (fig.1(e) for the (0 0 2) peak) reveals a continuous shift to the left, and therefore, an increase of interplanar clattice parameter with increasing Te-doping.Single crystal diffraction confirms the increase in c-lattice parameter and, in addition, reveals an increase of the in-plane a-lattice parameter with Te-doping.Table S2 shows a summary of a and c lattice constants found by X-ray diffraction.
The 2H-polytype is also confirmed by Raman spectroscopy.Figure 1(c) shows the Raman spectra for the different Te concentrations studied (x=0-15%), around the E2g and  1 modes for 2H-WSe2.A continuous decrease in the Raman shifts of both modes as tellurium concentration increases is observed, which also reveals a 2H-structure with continuously increasing a-and c-lattice parameters.Raman peaks associated with the 1Td-polytype, which for pure WTe2 appear around 150 cm -1 and 200 cm -1 , are not present, implying that secondary phase domains with this structure are not detected.Our results of the crystallographic structure align with a previous study in W(Se1-xTex)2 which reports a pure 2H-structure for Te compositions below 50% (30).

Magnetic response
Magnetization measurements as a function of magnetic field at room temperature are shown in figures 2(a) and S3.These reveal two contributions to the magnetic signal: a diamagnetic component and a paramagnetic or ferromagnetic component, depending on the chalcogen vacancy level.The paramagnetic behavior has been previously reported by other groups in WSe2 and other TMDs (21,22), and the paramagnetic or ferromagnetic contribution have been theoretically predicted (31) and experimentally realized in this family of materials (21)(22)(23)(24)(32)(33)(34).Measurements of samples with different times of exposure to ambient conditions suggests that the diamagnetic component is highly dependent on the level of oxidation of the samples.For instance, freshly synthesized samples showed no or minimum diamagnetic contribution, whereas the same samples after being stored in soft vacuum for a couple of days showed a marked diamagnetic contribution (fig.S4).The paramagnetic/ferromagnetic contribution appears to be more robust to sample oxidation (figs.S4 and S5, and table S3), and therefore the diamagnetic contribution of all curves in fig.2(a) was subtracted in order to focus on the paramagnetic/ferromagnetic component.Samples with no Te doping and low concentration of chalcogen vacancies (black curve in the inset to fig.2(a)) show a paramagnetic behavior, with a saturating magnetization above fields of the order of 150 Oe.Te-free samples with an increased number of Se vacancies showed weak hysteresis loops with small coercive fields, indicating weak ferromagnetism (see red curve in the inset to fig.2(a), and fig.S4(a-c) for other vacancy levels).For all the Te-doped samples studied, the ferromagnetic hysteresis loops became more notorious, reaching a maximum coercive field of 700 Oe for the sample with x=7.6% and δ=38%.Figures 2(b,c) plot the coercive field for the different samples studied, as a function of Te doping (x) and chalcogen vacancies (δ), respectively.The coercive field as a function of chalcogen vacancies grows monotonically, in contrast to its behavior as a function of Te doping.This suggests that magnetism is driven by the amount of chalcogen vacancies.Nevertheless, the apparent strengthening of the ferromagnetism in the Te-doped samples, in contrast to the undoped samples, suggests some role of the Te atoms.We argue that this is due to the fact that the presence of Te in the synthesis promotes the formation of chalcogen vacancies (Table S1).For instance, all the Tefree synthesis reached maximum vacancy levels of 12%, even though for some of them nominal vacancy values of 50% were aimed.Whereas for the Te-doped synthesis, measured vacancy levels above 12% could be easily reached, with a maximum obtained vacancy level of 38%.

Piezo-and ferro-electric response
Piezoelectricity in these compounds was investigated using an atomic force microscope (AFM) through the Dual AC Resonance Tracking (DART) piezo force microscopy (PFM) technique (see methods).Freshlycleaved samples were measured in an inert N2 atmosphere in order to prevent surface oxidation and minimize electrostatic effects.This technique measures the amplitude of elongation for a given applied AC bias, and therefore it is a direct measurement of the piezoelectric response.Figures 3(a-c) show DART piezoresponse maps for different values of the AC bias voltage, in single crystals with three different doping levels.The color in each plot represents the amplitude of out-of-plane deformation, which spatial average increases with the AC applied bias.The spatiallyaveraged piezoresponse as a function of AC bias, presented in figure 3(d), is linear for all samples.The slope in each curve is equivalent to the effective piezoelectric constant d33 (see methods, and section S4).The two lowest Teconcentration samples exhibit comparable values for their effective d33, of the order of 4 pm/V.This value is similar to the piezoelectric constants previously reported for TMD monolayers (26,35), and higher than the d11 coefficient reported for a WSe2 monolayer (25,28).For the highest Te-doped sample, with x=15% and δ=21%, d33 = 26±4 pm/V.This remarkably high value is comparable to the d33 of PPLN -a device based on LiNbO3, a widely used piezoelectric material (27) -and it is the highest reported or predicted d33 value among the TMDs in any configuration (36,37).It is important to mention that previously reported piezoelectric properties in the TMD's have been measured in nanostructures of these compounds, generally in single layers in which centrosymmetry is broken, and recently, in thin films of TMD alloys (38).Our measurements are performed in macroscopic bulk single crystals, and as such, constitute, to the best of our knowledge, the first report of piezoelectricity in the bulk of TMDs, with piezoelectric coefficients that are comparable to the materials used in commercial piezoelectric devices.
In order to further characterize the electromechanical properties in these compounds, switching spectroscopy (SS)-PFM measurements were performed (see methods and section S5 for experimental details).The phase (figs.3(e,g,i)) and amplitude (figs. 3(f,h,j)) of the SS-PFM piezoresponse were recorded as a function of the DC bias voltage at different random locations over the (insets to figures 3(e,g,i)).The (x=1.4%, δ=9%) Tedoped crystal (figs.3(e,f)) reveals a piezoelectric response without hysteresis.Interestingly, the amplitude tends to saturate for |VDC|>1V, suggesting a possible saturation of the electric dipole moment response.The (x=7.6%, δ=38%) and (x=15%, δ=21%) Tedoped samples show clear hysteresis loops both in the amplitude and phase of the SS-PFM piezoresponse (figs.3(g,h) and 3(i,j), respectively), with shapes equivalent to the ones shown in prototypical ferroelectric materials such as BaTiO3 and PZT (39).These loops indicate the presence of domains of electric polarization that can be aligned, with a coercive voltage of 0.2 V and 0.5V for the x=7.6% and 15% Te-doped bulk samples, respectively, and therefore reveal that these materials are ferroelectric.
In addition, electrical current vs voltage (IV) measurements taken in the Te-doped samples reveal asymmetric and hysteretic curves (fig.S8), with similar characteristics to those measured in the prototypical ferroelectric BaTiO3 and multiferroics such as BiFeO3 (40).This makes Te-doped WSe2 suitable for a variety of multifunctional applications in electronic devices (41).
As presented in a previous section, the crystalline structure of our samples is the centrosymmetric 2H polytype.Previous reports of piezoelectric behavior in monolayer TMD's, and ferroelectric behavior in WTe2 have attributed the non-centrosymmetric structures in those systems as the underlying mechanism for the existence of electric dipole moments.Therefore, these results invite us to search for alternative mechanisms for ferroelectricity in the TMD's, as will be presented in the computational modeling and discussion sections.

Computational modeling results
In order to gain further insight into the mechanisms that can lead to the multiferroic properties of W(Se1-xTex)2(1-δ), first principles calculations were performed using density functional theory (DFT) (see methods).For each concentration random configurations were selected to perform the calculations.Resulting total magnetizations for different levels of Te doping and chalcogen vacancies can be found in Table S4, for one particular realizations of each composition.Interestingly, the system with the highest magnetic moment corresponds to the one with the highest value of chalcogen vacancies, δ, independent of the value of Te doping, x.
Figure 4 shows different views of two layers of a simulated crystal structure with x≃10.26% and δ≃39%, very close to the values of one of our experimentally characterized samples that shows the strongest ferromagnetism and intermediate ferroelectricity.Arrows in figure 4(c) represent the magnetic moment contribution per atom for the atoms that contribute the most magnetic moment in this configuration.In this figure, in order to facilitate visualization, bonds between W atoms are drawn if the interatomic distance is equal or less than the bond distance in crystalline tungsten (2.74 Å); and bonds between W and chalcogen atoms are drawn if the interatomic distance is less than 2.72 Å and 2.747 Å for W-Se and W-Te bonds, respectively.The W atoms that are not bonded to a chalcogen element in the figure (labeled as W24 in fig.4(a), and W1 and W5 in fig.4(b)) are the ones that contribute the most magnetization of all atoms in the structure (as represented by the length of the gray arrows in fig.4(c)).The next-higher contributions come from W atoms that have a reduced number of bonds with chalcogen atoms (this is, W atoms labeled as 17, 21, 23, 25, 28 and 29 in fig.4(a), and 4, 6 and 16 in fig.4(b), displayed as red arrows in fig.4(c)).The magnitude of the contributions seems to be unaffected by whether the W is bonded to Te or Se atoms.Therefore, the relevant parameter that drives the magnitude of the magnetic moment is the number of bonds of each W with neighboring chalcogen atoms.Although the criteria for drawing a bond in figure 4 is rather arbitrary, a different number of bonds for each W ion is an indication of a different local environment for each.Therefore, our calculations indicate that the magnitude of the W magnetic moment in this system is driven by the local environment of the W atoms, and the strength of its hybridization to chalcogen atoms (or alternatively, its bonding to a chalcogen vacancy).Figure 4(d) shows the local spin density, ns(r)=nu(r)-nd(r), where nu (nd) is the density of spin-up (spin-down) electrons, plotted on two planes containing W1, W5 and W24 (atoms without bonds to chalcogen atoms).It is important to notice that throughout the unit cell the spin density keeps low values ~±5x10 -4 μB (pink) while on these 3 atoms the spin density peaks at ~8.8x10 -2 μB (violet).For the atoms not bonded to a chalcogen element in figure 4(a-c), atomic magnetization (the integral of the spin density up to the covalent radius) takes values as high as 0.56 μB.Furthermore, the overall ferromagnetic character of the system is corroborated by means of the total density of states (DOS) for spin-up and -down (fig.4(e)), which shows that at the Fermi energy there is a larger contribution from the spin-up channel (red) than from the spin-down channel (blue).The projection of the DOS on atomic orbitals allows for the identification of tungsten d-orbitals as the electronic states associated with the found magnetic order, as shown in figs.4

(f) and S9(a).
In order to corroborate if magnetic order is dependent on the positions of chalcogen vacancies and/or Te substitutions, calculations were performed for a different configuration keeping the same concentration of vacancies and Te doping.In this second configuration (fig.S9(b)) there is a more homogeneous distribution of vacancies which causes that all W-atoms are bonded to at least one chalcogenatom, either Se or Te.The resulting supercell magnetization is null in this case, with almost negligible atomic magnetizations (around 10 -4 μB) with opposite directions uniformly distributed throughout the cell.
In order to find possible sources for the ferroelectric state, electric polarization and ferroelectric behavior were analyzed in the framework of the modern theory of polarization (42), which allows to obtain the unit-cell polarization from the Berry phase of the system.The electric polarization vector was calculated in the same systems for which magnetization was studied.Results show a dependency on vacancy concentration and Te doping as can be seen in Table S5.Additionally, Hirshfeld charge analysis (43) was carried out in order to identify both the atoms responsible for polarization and the possible mechanism of ferroelectricity.
In the case of pristine WSe2 it is found that chalcogen atoms have positive Hirshfeld charges while W atoms gain a negative charge which is approximately equal to twice the charge of the chalcogens.This is expected from the coordination numbers of Se and W, three and six respectively (Table S5).This distribution of charge makes each layer of 2H-WSe2 non-polar due to its centrosymmetry, which is broken by the presence of vacancies.As can be seen in Figure 5, the absence of a Se atom develops electric polarization directed from the negative W plane towards the remaining Se atom below or above the vacancy.Since the sign of the polarization depends on the position of the chalcogen atom, this result suggests that polarization could be reverted by moving the chalcogen atom (or equivalently, the chalcogen vacancy) from one face of a layer to the other face, passing through the triangular space defined by the neighboring W atoms.This establishes a possible mechanism for the presence of switchable electric dipole moments in WSe2.Interestingly, even though for this configuration there is a high concentration of vacancies (δ≃39%), there is only one place with a non-mirrored vacancy located at the top layer.The remaining vacancy sites are located on both sides of both layers of the unit cell, therefore not breaking its centrosymmetry, as can be seen in the bottom layers of the unit cell in figure 5.
To study the effect of Te doping on the electric polarization, Se5 (as labeled in Fig. 5) was exchanged by a Te atom.The resulting Hirshfeld charge of this Te is ~2.61 times larger than that of Se5 and the total polarization of the supercell presents a ~20 times enhancement if compared to the previous system.The larger separation of charge seen with Te can be ascribed to the difference in its electronegativity (χ) with respect to both Tungsten and Selenium.In Pauling units χ Se =2.55, χ W =2.36 and χ Te =2.1, which accounts for the larger migration of charge from Te to W atoms.An additional configuration with a higher concentration of Te and less vacancies was studied (x=16% and δ=21.88%).In this case there are 13 non-mirrored vacancies, and chalcogen atoms are distributed on both faces of the two layers of the unit cell, 9 on the superior face and 4 in the inferior face.In general, Hirshfeld charges of Te range from 2 to 3 times the charges of Se, which gives larger local electric dipole moments.Since several local dipoles have opposite directions depending on the position of the vacancies, the computed unit cell polarization is about the same value found for the configuration with only one non-mirrored vacancy.In the case where all local dipoles point in the same direction, a large unit cell polarization is expected.

Discussion
Ferromagnetic states have been predicted theoretically (31) and experimentally realized in TMD through defects formation (21), vacancies formation (23), proton irradiation (33), nanostructuration (34), and doping with magnetic impurities (22,24).In the case of this study, the dependency of the coercive fields on Te doping and chalcogen vacancies suggests that the magnetic properties of bulk W(Se1-xTex)2(1-δ), are driven by the amount of chalcogen vacancies, rather than by the presence of Te.This result aligns with the mechanism for ferromagnetic order derived from our DFT calculations, which reveal that the spin density is strongly localized onto the dorbitals of the W-atoms less bonded with chalcogen atoms (or equivalently, more bonded to chalcogen vacancies), as presented in figures 4(f) and S9.This suggests that unpaired/localized d-electrons in W are responsible for the observed magnetic behavior.The key role of chalcogen vacancies in the magnetic order of TMDs has been previously predicted and experimentally observed in systems such as MoS2 (44), WS2 (23) and monolayer V-doped WSe2 (32).For instance, the presence of vacancies in 2H-MoS2 nanosheets has been shown to transform the crystalline structure around the vacancies into a 1T-structure (44).In this environment the Mo 4+ ions have a net magnetic moment, in contrast to stoichiometric 2H-MoS2 in which the Mo 4+ ions are non-magnetic.Data suggests that the origin of the measured ferromagnetism in this system is intimately connected to an enhanced Mo 4+ -Svacancy exchange interaction.Similar conclusions have been recently reached in monolayer V-doped WSe2, in which ferromagnetic properties are significantly enhanced by promoting Se-vacancies in the structure (32).In the case of W(Se1-xTex)2(1-δ), the role of chalcogen vacancies can be completely analogous to their TMD's counterparts.This is, a local magnetic moment is generated in the d-orbitals of transition metal atoms surrounded by vacancies due to their modified local environment.This can result in a strengthening of the exchange interaction through their coupling with chalcogen vacancies, leading to the ferromagnetic state.
On the other hand, ferroelectric states have been commonly associated with the noncentrosymmetry inherent to certain crystal structures.Among the TMD's, piezoelectric and ferroelectric properties have been previously reported in crystals with the noncentrosymmetric 1Td-phase, like the case of pure WTe2 (20), or in nanostructures of the 2Hpolytype with an odd number of layers, for which centrosymmetry is broken (25,26).However, alternative mechanisms, other than the underlying bulk crystalline symmetry, have been proposed for the appearance of electric polarizations in ferroelectric materials.For example, in type-II multiferroics such as TbMnO3 (45), Ca3CoMnO (46) and the 2D material NiI2 (14), local electric polarization is proposed to be induced by the strong spin-orbit coupling in a magnetic host, which can lead to inversion symmetry breaking, and to a ferroelectric state.In addition, and for the specific case of TMD's, inequivalent Te-W distances in adjacent 1T-WTe2 layers has been theoretically shown to induce a charge imbalance, leading to an out-of-plane electric polarization (47).In the case of this work, our DFT calculations identify a feasible mechanism for the existence of electric dipole moments in the 2H unit cell that can lead to a non-zero total electric dipole moment in the crystal.In our model, the electric dipoles are originated in the charge imbalance created by a chalcogen vacancy aligned with a chalcogen atom in the c-direction.Our calculations indicate that such imbalance is an order of magnitude stronger if the chalcogen is a Te-atom, which is consistent with our experimental results: no ferroelectric state is found in the pure WSe2(1-δ), in which, although the charge imbalance mechanism due to vacancies is also present, it is possibly not strong enough to result in a ferroelectric state.Whereas ferroelectricity grows stronger with Te doping, being well established for the largest Te-doped crystals.
For W(Se1-xTex)2(1-δ), the role of chalcogen vacancies seems to be crucial for both ferromagnetic and ferroelectric states.Interestingly, each of these ferroic orders in this material can exist independent from the other: WSe2(1-δ) is only ferromagnetic and WTe2 is only ferroelectric.This, together with the fact that multiferroicity in W(Se1-xTex)2(1-δ) is observed at room temperature, suggests that, although part of the origin of both ferroic orders can be common, this material is not a type-II multiferroic.In this class of multiferroics the mechanisms for both ferroic orders are highly intertwined through spin-orbit coupling, and generally show cryogenic critical temperatures (45).Nevertheless, the strong influence of chalcogen vacancies on both ferroic orders in W(Se1-xTex)2(1-δ) makes possible, not only the simultaneous presence of these states, but also their magnetoelectric coupling -a crucial ingredient in several envisaged applications of multiferroic materials.These ingredients, combined in an intrinsically 2D vdW layered material, and at room temperature, opens the door to a wide use of nanostructured multiferroic devices.

Synthesis of single crystals
Tellurium-doped tungsten diselenide single crystals, W(Se1-xTex)2(1-δ) (with x=0-15% and δ=1-38%), were grown by chemical vapor transport (CVT) with iodine as transporting agent.Powders with the desired stoichiometry were produced from elemental W, Se and Te through a sintering process.For each batch, stoichiometric amounts of tungsten (No. 357421, Sigma Aldrich), selenium (No. 36208, Alfa Aesar) and tellurium (No. 266418, Sigma Aldrich) were thoroughly ground, and then cold-pressed at 1000 psi to form a pellet.The pellet was sealed in a quartz tube, leaving a small pressure of argon inside the quartz ampoule.The ampoule was then heated from room temperature to 500°C in 24h and maintained at this temperature for 72h.The obtained sintered powder and the iodine were then encapsulated in vacuum in a long quartz tube, and then placed in a two-zone tube furnace with 1010°C in the hot-zone and 900°C in the cold-growth-zone, for 144 hours.The single crystals were collected from the growthzone of the tubes.

X-ray diffraction
The doping values x and δ of our W(Se1-xTex)2(1δ) single crystals were determined by triplication in X-ray fluorescence (XRF) using a ZSX Primus RIGAKU spectrometer.
X-ray powder diffraction patterns of W(Se1-xTex)2(1-δ) single crystals for all compositions were obtained using an Empyrean PANalytical series 2 diffractometer with a CuKα radiation source (λ = 1.5405980Å), operated at 40 mA, 45kV, and with a step size of 0.0262606° over a 2θ range from 10° to 70°, in an Eulerian-Cradle geometry.This characterization was performed in a collection of single crystals of each batch, which were cut in tiny fragments and spread in the sample holder to emulate a powder diffraction experiment.This was done with the purpose of obtaining a statistically meaningful characterization of many crystals of each batch.The most intense diffraction peaks for all batches come from the {0 0 1} family of planes, given that the a-b plane of the crystals tend to align parallel to the sample holder.Few peaks corresponding to different families of planes can still be recognized.
In order to resolve the exact crystal structures of the x=1.4% and 7.6% Te-doped samples, single crystal X-ray diffraction was also measured.The intensities were measured at room temperature, 298 (2) K, using CuKα radiation (λ = 1.54184Å), and ω scans in an Agilent SuperNova, Dual, Cu at Zero, Atlas four-circle diffractometer equipped with a CCD plate detector.The collected frames were integrated with the CrysAlis PRO software package (CrysAlisPro 1.171.39.46e,Rigaku Oxford Diffraction, 2018).Data were corrected for the absorption effect using the CrysAlis PRO software package by the empirical absorption correction using spherical harmonics, implemented in the SCALE3 ABSPACK scaling algorithm (CrysAlisPro 1.171.39.46e,Rigaku Oxford Diffraction, 2018).The structures were solved using an iterative algorithm (48), and then completed by difference Fourier map.The crystal structures were refined using the program SHELXL2018/3 (49).

Raman spectroscopy
Raman spectra for crystals of all the compositions studied were taken in a HORIBA Scientific XPLORA Xl041210 Raman spectrometer.The excitation wavelength used for all measurements was 532 nm with a grating of 2400 lines/mm, in the range of 70 cm -1 to 1000 cm -1 .

Magnetic measurements
Isothermal magnetic hysteresis loops at 300K and 80K were measured in a Lakeshore TM 7400 Series vibrating-sample magnetometer (VSM).Cryogenic temperatures were reached through a flow-cryostat operated with liquid nitrogen.

Piezoresponse Force Microscopy Measurements
The piezoelectric effect can be measured through piezoresponse force microscopy (PFM) technique, which is a variation of the atomic force microscopy (AFM) technique.The piezoresponse is acquired when an AC electrical voltage is applied using a conductive tip in contact with the sample surface.The resulting oscillatory deformation of the sample under the AC voltage is detected through the AFM cantilever deflection.We use two kinds of PFM measurement: Contact resonance (Dual AC Resonance Tracking-DART) and Switching Spectroscopy -SS-PFM -, using a Cypher ES Environmental AFM operated at room temperature and in an inert N2 atmosphere.
In the DART-PFM mode the cantilever is operated with a AC voltage frequency near the contact resonance frequency, resulting in an driven harmonic oscillator, which enhances the piezoresponse signal, even if it is very small.Similarly, SS-PFM uses a sinusoidal voltage at contact resonance frequency overlapped with square voltages of smaller periods.In ferroelectric materials, it allows the determination of the electromechanical hysteresis loops and their switching parameters such as coercive and nucleation voltages.
To avoid misunderstandings in piezoresponse amplitude data, arbitrary units (a.u.) are used when piezoresponse amplitude is reported without treatment, whereas picometers (pm) are used when amplitude data is normalized by a single harmonic oscillator (SHO) model (50).This model provides an accurate way to calculate the real piezoelectric amplitude in the small damping regimen.Moreover, to demonstrate the accuracy of this model and its use in a methodology to obtain the effective piezoelectric constant, d33 in our single crystals, an standard sample of periodically poled lithium niobate (PPLN) was characterized by DART-PFM and SS-PFM using silicon probes (AC240TM-R3) (sections S4 and S5 for more details).Afterwards, the same methodology was used in our doped single crystals (section S5 for data without data analysis).Experimental spring constants were between 1.71-2.95N/m,and contact-mode resonance frequencies were between 219-241kHz.

Transport properties
Electrical properties of single crystals were measured by a four-probe technique for resistance vs temperature/field measurements, and two-probe technique for current vs voltage measurements, using a Keithley 2400 Source Meter.Gold pads were evaporated into the single crystals in order to reduce contact resistance and capacitive effects.Cryogenic temperatures were achieved using an Oxford Instruments liquid helium variable temperature insert cryostat equipped in an IntegraAC Recondensing Helium Magnet System.

DFT calculations
First principles calculations were performed using Density Functional Theory (51) (DFT) as implemented in the Quantum ESPRESSO package (52).Core electrons were treated within the pseudopotential approximation by means of norm-conserving ultrasoft pseudopotentials.
Generalized gradient approximation (GGA) in the parametrization of Perdew-Burke-Ernzerhof (PBE) (53), was chosen for the exchange-correlation (XC) energy functional.Grimme's semiempirical DFT-D3 method (54) was used to include van der Waals interactions between the layers, ignoring three-body terms.The number of plane-waves in the basis set correspond to a kinetic energy cutoff of 50 Ha.Structural optimizations were performed until forces reached a tolerance of 1x10 -4 a.u.
Starting with the simulation of pure WSe2, ionic positions and unit cell lattice parameters were optimized.Resulting structural parameters are in good agreement with the measured ones and with previously reported DFT calculations (62), as can be seen in table S4 of section S7.Inclusion of both, chalcogen vacancies and Tedoping, was achieved in a 4x4x1 supercell, which allowed the simulation of systems with Te-doping and chalcogen vacancy levels similar to the measured single crystals.The effect of vacancies and Te doping on the magnetic behavior was analyzed by performing spin-polarized calculations with and without spin-orbit coupling.For each concentration a random configuration was selected to perform the calculations.In all cases a 4 x 4 x 3 Monkhorst-Pack grid (55) was used for sampling the first Brillouin zone.
Berry phase calculations of electric polarization were performed with the implementation found in the ABINIT(56) program using a finer reciprocal space grid of 5 x 5 x 3 points and ONCVPSP (57) pseudopotentials.Hirshfeld charge analysis was performed with the postprocessing tool cut3d from ABINIT.and secured funds.All authors participated in the discussion of the data and manuscript writing.
Competing Interests: The authors declare no competing financial interests.(29).(e, f) Close-up look to the (0 0 2) and (1 0 3) peaks, which reveal a shift to the left for the Te-doped samples, indicative of an increase of 2.5% in the c lattice parameter           Fig. S3(a) shows the total magnetization in emu/g at 300K for x=0-15% tellurium-doped compositions.
Here, the magnetic response reveals two components: a diamagnetic and a ferromagnetic contribution.The ferromagnetic component, obtained by subtracting the linear-diamagnetic contribution of each curve, has a slight temperature dependence, particularly in the saturation magnetization, as evidenced in fig.S3(b) which shows measurements at 80 K and 300 K for different Te-compositions.
On the other hand, the diamagnetic component has a strong dependence with the level of oxidation of the samples.Figs.S4(a,b) show the total magnetization curves for two types of samples, comparing the results for two different levels of exposure to ambient conditions: short after synthesis (short exposure-time to ambient conditions), and after several days of storage in soft vacuum.For the curves it is clear that the diamagnetic slope is developed only after a prolonged exposure to ambient conditions, whereas the ferromagnetic contribution of the magnetization is present with and without exposure to oxygen.Moreover, the ferromagnetic contribution seems to be more robust to slight oxidation of the samples, as shown in figs.S4(c,d).These curves reveal that the coercive field of the ferromagnetic hysteresis loops present only slight variations with the level of exposure to ambient conditions.Table S3 shows the variation in coercive field and saturation magnetization of the ferromagnetic contribution to magnetization, for different times of exposure to oxygen, for four different types of samples.Although the saturation magnetization presents variations with the time of exposure to oxygen, the coercive field is very robust to the oxidation level, and this is the reason why this quantity was the one chosen to study the evolution of the magnetism with Tedoping and chalcogen vacancy level.The slope of the linear regression to this curve corresponds to the piezoelectric coefficient, d33.For our sample of PPLN, the experimental result d33 is 19.00 ± 1.18 pm/V.This value is similar to the reported values for d33 in multiple piezoelectric studies for PPLN (27).Therefore, this demonstrate that the presented method is reliable for the determination of effective piezoelectric constants.

Section S5. Piezoresponse force microscopy hysteresis loops
Switching Spectroscopy (SS) PFM hysteresis loops for all tellurium-doped compositions, and an undoped sample (x=0%, δ=4%), are compared in fig.S7.Given that the undoped sample -which is not piezoelectric -was measured in the same conditions as the Te-doped samples, it can be used to identify the contribution of electrostatic effects during the measurement.It is worth noting that those contributions are negligible when compared to the experimental signal in our Te-doped compositions.Significantly, the amplitude of the SS-PFM signal increases enormously with increasing tellurium doping, which confirms the large enhancement of the piezoelectric coefficient for the highest Te-doping sample.Another indication of ferroelectric properties in our single crystals can be found in the current vs voltage (IV) curves.Fig. S8 shows I-V measurements taken in samples from x=0% to 15%, at 300K and 50K.Section S9.List of suporting files.
The crystallographic information filed -CIF -of the resolved crystal structure found by single crystal diffraction, for two of the studied compositions, can be found in the files:

Figure 1 .
Figure 1.Crystalline Structure.(a) Representative single crystals of W(Se1-xTex)2(1-δ) (this case, for x=15% and δ =21%) synthesized by chemical vapor transport, and characterized in this work.Scale bar is 1 mm long.(b) Crystal structure of a 2H-polytype of Te doped WSe2 with space group P63/mmc (#194).(c) Raman spectroscopy measurements for W(Se1-xTex)2(1-δ) single crystals of the compositions studied in this work.Data shows a close up look to the Raman shift around the characteristic E2g and A1g peaks of the 2H phase, which shift lower values with increasing doping.Dashed gray lines indicate the peak positions for the pure x=0 compound.(d) Powder X-ray diffraction data for a collection of undoped (x=0, δ=4%), and highest Te-doped (x=15%, δ=21%) samples.Diffraction peaks are consistent with the 2H polytype, as compared with the reported peak positions of pure 2H-WSe2, shown by the blue vertical lines(29).(e, f) Close-up look to the (0 0 2) and (1 0 3) peaks, which reveal a shift to the left for the Te-doped samples, indicative of an increase of 2.5% in the c lattice parameter

Figure 2 .
Figure 2. Magnetic properties.(a) Magnetization hysteresis loops at 300 K for samples of W(Se1-xTex)2(1-δ) with different x and δ values.For all curves, the diamagnetic components were subtracted, and their values normalized by the saturation magnetization, Ms. Inset in (a) shows the hysteresis loops for undoped (x=0) WSe2(1-δ) samples with different amounts of chalcogen vacancies, δ.(b) and (c) show coercive fields Hc for all samples presented in (a), as a function of (b) Te-doping, x, and (c) chalcogen vacancies, δ.

Figure 4 .
Figure 4. Structural and magnetic calculations of W(Se1-xTex)2(1-δ) for x=10.26% and δ=39.06%.W-W and W-Te bonds are drawn if the interatomic distances are less than 2.74 Å and 2.81 Å, respectively (see text).Top views of the (a) superior and (b) inferior unit cell monolayers and (c) perspective view of the optimized structure.Arrows indicate resulting atomic magnetic moments more than 0.02 μB in magnitude.(d) Spin density plotted on planes containing atoms W1, W5 and W24.To facilitate visualization, atom W24 from a neighboring unit cell is depicted.(e) Total spin-up (red) and spin-down (blue) density of states.(f) Projected density of states on p orbitals of Se and Te atoms, and on d orbitals of W.

Figure 5 .
Figure 5. Electric polarization calculations.(a) Different views of a (0 1 0) plane of a W(Se1-xTex)2(1δ) simulated structure with x=10.26% and δ=39.06%.The top view is a 3D-perspective view of the atoms around a chalcogen vacancy located in the upper side of the top-layer of the unit cell.The bottom view shows both layers of the unit cell (marked in the dashed square), and the color scale represents a slice of the total electronic density (color bar scale in a.u.).Arrows indicate the direction of the electric dipole moment around the vacancy, in the c-direction.(b) Same representations as in (a), but now with the chalcogen vacancy located in the bottom side of the top layer of the unit cell.The direction of the electric dipole moment is flipped in this case, as indicated by the arrow.

Section S2 .
Raman and XRD spectraIn fig.S1, Raman spectra are shown for all Te-doped compositions and an undoped one, for a wide range of Raman shifts.A decrease in the Raman shift of the E2g and A1g modes of the 2H-structure, with increasing Te-concentration is observed.

Figure S1 .
Figure S1.Raman spectra for all tellurium-doped single crystals involved in this study.

Figure S2 .
Figure S2.XRD patterns for all tellurium-doped single crystals involved in this study.

Figure S4 .
Figure S4.(a,b) Total magnetization hysteresis loops, for two different times of exposure to ambient conditions, for samples without Te-doping and vacancy level of (a) δ=12%, and (b) δ=8.7%.(c,d) Ferromagnetic contribution to the magnetization for two different times of exposure to ambient conditions, for (c) the undoped sample with δ=8.7% and (d) the Te doped sample with x=15%, δ=21%.

Section S4 .
Figure S5.(a) DART-PFM piezoresponse normalized by SHO model as a function of the AC applied bias for PPLN material, and (b) their resultant effective d33.Scale bar: 1 μm.

Fig. S6 shows
Fig.S6shows additional PFM measurements on the same PPLN sample.DART-PFM measurements (figs.S6(a-c)) allow the identification of poled domains in this material, which are not visible on the topography image.These domains were investigated by SS-PFM over the marked points in fig.S6(c).Firstly, both, +z poled and -z poled regions were measured using DC pulses bias from 0 to ±1V.The SS-PFM phase shows that these domains are oriented in opposite directions.Their SS-PFM amplitudes reveal an opposite behavior on the sample deformation: a contraction for -z poled and an extension for +z poled (fig.S6(d)).

Figure S7 .
Figure S7.Switching Spectroscopy (SS) PFM hysteresis loops at 300K under an inert atmosphere for all tellurium-doped compositions.An undoped composition (black) is measured to illustrate contributions to the amplitude caused by electrostatic effects.

Table S5 .
Computed values for Hirshfeld charge of different atoms in the unit cell, and total unit cell electric polarization, for different values of Te-doping and chalcogen vacancies.

Figure S10 .
Figure S10.[1 0 0] views of the total electronic densities as obtained for the vacancy on top (left) and below (right) a Se atom.

Table S1 .
Chemical composition of single crystals, determined through X-ray Fluorescence Spectroscopy (XRF).

Table S2
presents the obtained lattice constants, obtained through powder x-ray diffraction and single crystal diffraction, as explained in the methods section.For the powder x-ray diffraction, the values shown in the table are the average over different peaks of the same family.

Table S2 .
Lattice constants for all compositions, found by X-ray diffraction measurements.

Table S3 .
Variation in ferromagnetic properties after soft vacuum storage.
S4_checkcif report for Te-7percent.pdfComputed atomic coordinates, lattice constants as well as Hirshfeld charge analysis for several configurations can be found in the following files: