Spin–orbit and exchange proximity couplings in graphene/1T-TaS₂ heterostructure triggered by a charge density wave

The ability to stack two-dimensional materials offers to design ultrathin van der Waals heterostructures [1–3] that can foster spintronics [4–9], optoelectronics [9–11], and twistronics [12, 13] applications and pave the way to novel carbon-based and, hence, environmentally-friendly devices. Modifying the properties of van der Waals materials through proximity effects is an added degree of functionality so far very intensively explored, particularly, proximitizing graphene leads to the appearance of induced spin–orbit coupling (SOC) [10], or inducing magnetic exchange proximity effects [14, 15].

In order to add a novel tuning knob to the van der Waals heterostructures, we investigate the effect of charge density wave (CDW) that develops in 1T-TaS₂ and explore its effect on the spin degree of freedom in proximitized graphene. In addition to various CDW orderings that are emerging in layered 3D bulk structures of transition metal dichalcogenides (TMDCs) [16], also their 2D counterparts [17–19] offer a potential platform for studying the effects stemming from the dimensionality, and CDW-triggered commensurability and shape-topography [20, 21]. In particular, some of TMDCs in the 1T polymorph phase have the propensity to exhibit a variety of periodic lattice distortions [22]. Frequently, the Peierls instability gives rise to a stable low temperature periodic CDW phase with a triangular $\sqrt{13}\times\sqrt{13}$ superlattice. Particular examples count 1T-TMDCs such as TaS₂ [23–25], TaSe₂ [25–27] and TaTe₂ [25, 27] or NbS₂ [28] and NbSe₂ [29–31]. Structural CDW distortions have the form of David stars [32]—clusters of 13 transition metal atoms in which an excess electron is effectively trapped at the star center [33–38]. What is particularly interesting, the commensurate $\sqrt{13}\times\sqrt{13}$ CDW phase of TaS₂ forming the superlattice of David stars [39] possesses below 180 K [40] enough electronic correlations [41] to develop the Mott phase [42–46] and even superconductivity under pressure [32, 47]. Formation of CDW in 1T-TaS₂ results from the subtle interplay of lattice, orbital and spin degrees of freedom [48], along with the electron–phonon coupling [49, 50] which results in strongly localized electronic states near the Fermi level [51]. The anharmonic phonon–phonon interactions in 1T-TaS₂ are predicted to stabilize CDW at elevated temperatures [52]. The correlations at low-temperatures lead to spin polarization and transform a monolayer of 1T-TaS₂ into a magnetically polarized semiconductor [53] with in-plane magnetization [54]. Rich electronic properties of 1T-TaS₂ can be engineered by the temperature [55–58], pressure [59], electrical pulses [60, 61], electrical bias [62, 63], creation of domain electronic walls [64], insertion of adsorbates [65, 66], charge doping [67], substrate selection [68], strain [69–71], or laser light [45, 72–74]. All this gives the 1T-TaS₂ monolayer a potential for 2D spintronics applications counting ionic field-effect transistor [75], logic gates and circuits [76], non-volatile memory devices [58, 77, 78], and others.

In addition to many promising studies of 1T-TaS₂ in the bulk or monolayer form, also its various heterostructures have been considered in recent literature, including the experimental works on the heterostructures with 2H-TaS₂ [79], 2H-MoSe₂ and 2H-WSe₂ [80], 2H-MoS₂ [81], permalloy [82], h-BN [83] or black phosphorus [84], and also graphene [41, 85–87]. The latter is providing an ideal platform for investigating fine spectral features stemming from the proximity of graphene Dirac electrons to the correlated CDW phase.

The CDW degree of freedom present in 1T-TaS₂ and its spatial topology provides a channel to control the proximity exchange coupling and SOC in graphene without any physical change in a twist of the heterostructure components, as the CDW phase can be manipulated locally, for example, by optical means [88–90], electrical pulses [77, 78], or just temperature and in a fully reversible manner. Such a feature has been only achieved in optical lattice analog to bilayer graphene [91]. Twisted graphene/TMDC heterostructures have been recently intensively studied [92–96] and control of the twist angle is of paramount importance for possible applications based on the design of spin-charge conversion utilizing the Rashba–Edelstein effect [8, 90, 97–99]. Here we show that the same level of functionality can be expected also for graphene/CDW states where the 'role of the twist' is played by the spatial topology of CDW.

In the paper, we study the van der Waals heterostructure of graphene and 1T-TaS₂ monolayer by means of first-principles calculations and propose an effective tight-binding model for graphene bands near the Dirac point extracting relevant model parameters. We found that the model parameters reflect induced SOC and exchange interaction proximity effects for transitions of the 1T-TaS₂ monolayer from the normal to CDW state, and then to the ferromagnetic ground state achieved by gradually decreasing temperature.

Based on a given structural symmetry one can derive a specific effective model Hamiltonian for graphene π-bands [100, 101]. It has been noted [92, 93] that the most general form of the induced Rashba SOC term in twisted graphene/TMDC heterostructures that obeys time reversal and threefold rotation C₃ symmetries can be written as $U^\dagger H_{\mathrm R} U$, where $H_{\mathrm R}$ is the effective Rashba Hamiltonian, $H_{\mathrm R} = \lambda_{\mathrm R}(\kappa \sigma_x s_y - \sigma_y s_x)$ [102, 103] and $U = e^{-i s_z \phi_{\mathrm R}/2}$, where $\sigma_{x,y}$ are the Pauli matrices acting in sublattice space, while $s_{x,y,z}$ are the Pauli spin matrices for the spin degree of freedom. The Rashba angle $\phi_{\mathrm R}$ originates from SOC matrix elements between d-orbitals of metal atoms in TMDC [92] and in general, represents a sum of the geometric angle and the quantum phase [104]. The results provided below are given with respect to the ordered Bloch basis:

$$\begin{equation} |A(\mathbf{q}),\uparrow\rangle, |A(\mathbf{q}),\downarrow\rangle, |B(\mathbf{q}),\uparrow\rangle, |B(\mathbf{q}),\downarrow\rangle, \end{equation} \tag{ 1 }$$

where the states are sublattice ($A, B$) and spin ($\uparrow$, $\downarrow$) resolved, and the crystal momentum q is measured from the Γ point.

The corresponding effective Hamiltonian in the studied cases has C₃ symmetry, and its low energy expansion close to the K $(\kappa = 1)$ and K^ʹ $(\kappa = -1)$ point is given as follows:

$$\begin{align} H_{\kappa}=\left(\begin{matrix} E_0+\Delta+\kappa \lambda_{\mathrm I}^{ A} & 0 & \frac{\sqrt{3}}{2}\mathrm {a t} \left(\kappa k_x-ik_y\right) & i(1-\kappa)\lambda_{\mathrm R} e^{-i\phi_{\mathrm R}} \\[2pt] 0 & E_0+\Delta-\kappa \lambda_{\mathrm I}^{ A} & i(1+\kappa)\lambda_{\mathrm R} e^{i\phi_{\mathrm R}} & \frac{\sqrt{3}}{2}\mathrm {a t} \left(\kappa k_x-ik_y\right) \\[2pt] \frac{\sqrt{3}}{2}\mathrm {a t} \left(\kappa k_x+ik_y\right) & -i(1+\kappa)\lambda_{\mathrm R} e^{-i\phi_{\mathrm R}} & E_0-\Delta-\kappa \lambda_{\mathrm I}^{ B} & 0\\[2pt] -i(1-\kappa)\lambda_\mathrm {R} e^{i\phi_{\mathrm R}} & \frac{\sqrt{3}}{2}\mathrm {a t} \left(\kappa k_x+ik_y\right) & 0 & E_0-\Delta+\kappa \lambda_{\mathrm I}^{ B} \end{matrix}\right). \end{align} \tag{ 2 }$$

The orbital and spin proximity effects result in the following parameters: E₀ corresponds to the Dirac point offset with respect to the Fermi energy, Δ is the staggered potential discriminating electronic states located at sublattice A and B, $\lambda_{\mathrm I}^{ A}$ and $\lambda_{\mathrm I}^{ B}$ are the sublattice-resolved spin-conserving next-nearest neighbor SOCs, $\lambda_{\mathrm R}$ is the Rashba SOC parameter with phase angle $\phi_{\mathrm R}$. Moreover, orbital nearest neighbor hopping t and the lattice constant a complete the set of parameters for the Hamiltonian as those can be subjected to CDW distortion. We should mention that the sublattice-resolved SOC parameters can be expressed as $\lambda_{\mathrm I}^{ A} = \lambda_{\mathrm I}+\delta\lambda_{\mathrm I}$ and $\lambda_{\mathrm I}^{ B} = \lambda_{\mathrm I}-\delta\lambda_{\mathrm I}$, where $\lambda_{\mathrm I}$ is conventionally called intrinsic, or Kane–Mele SOC.

The above-described Hamiltonian (2) was used to model graphene band structure in the vicinity of the $K/K^{^{\prime}}$ points for normal and CDW phase of 1T-TaS₂. The model parameters were determined to fit the density functional theory (DFT) calculations. A partially filled band in the vicinity of the Fermi level tends to spin split [53] at low temperatures and strong correlations of the Ta d-electrons play an important role in the positioning of the lower and upper Hubbard bands [47, 56, 86, 105, 106]. This has an effect on the emergence of the in-plane magnetic polarization of the 1T-TaS₂ [54]. In order to describe the SOC and also the exchange proximity coupling of the graphene bands around the Dirac points—we model the magnetic interaction between graphene and 1T-TaS₂ Hubbard band in terms of effective Zeeman-like interactions captured by the following Hamiltonian:

$$\begin{align} H_{\mathrm {e}x}=\begin{pmatrix} 0&\Delta_{ A}{ e}^{-{ i}\phi_{ A}}&0&0\\ \Delta_{ A}{ e}^{{ i}\phi_{ A}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \Delta_{ B}{ e}^{-{ i}\phi_{ B}} \\ 0 & 0 & \Delta_{ B}{ e}^{{ i}\phi_{B}} & 0 \end{pmatrix}. \end{align} \tag{ 3 }$$

The effective Hamiltonian $H_{\mathrm {e}x}$, independent of the valley, represents the coupling of a spin with the sublattice-resolved (subscripts $A/B$) in-plane magnetization that is parameterized by the amplitude $\Delta_{ A/B}$ and the angle $\phi_{ A/B}$ between the x-axis and the magnetization direction.

We consider a heterostructure with a structural model formed by the supercell of $\sqrt{13}\times\sqrt{13}$ for 1T-TaS₂ and $5\times 5$ for graphene, see figure 1. We study a slab geometry with a vacuum of about 15 Å in the direction perpendicular to the planes in order to separate periodic images. The graphene lattice parameter was compressed by approximately 1.4% to appropriately fit the supercell, capturing the commensurate CDW in 1T-TaS₂, using the experimental lattice constant of 12.1323 Å [107]. The commensurate CDW can be described by a wavevector $\mathbf{q}_{\mathrm {CDW}}$ whose orientation stems from the emerging of $\sqrt{13}\times\sqrt{13}$ reconstruction and rotation angle of 13.9^∘ [107, 108], see figure 1(b).

Zoom In Zoom Out Reset image size

All the electronic structure calculations were performed using the Quantum Espresso suite [109, 110], implementing a plane wave basis for DFT calculations [111]. The Perdew–Burke–Ernzerhof exchange-correlation functional was utilized [112] for the projector augmented wave method [113]. For non-collinear DFT calculations including SOC, fully relativistic pseudopotentials were applied [114]. Small Methfessel–Paxton energy level smearing [115] of 1 mRy was used along with the kinetic energy cut-off of 53 Ry. The semiempirical van der Waals corrections were included [116, 117] and the dipole correction [118] was applied to properly determine Dirac point energy offset due to dipole electric field effects. The relaxation of atomic positions within the spin-unpolarized DFT approach was allowed with a threshold force of 10⁻⁴ Ry per Bohr radius. For spin-polarized DFT+U calculations [119], we used on-site Hubbard U = 3 eV for Ta d-orbitals, inspired by the calculations reported in [21, 87], with $9\times 9$ mesh of k-points used for sampling the first Brillouin zone.

We study two stacking configurations of graphene on 1T-TaS₂, particularly, the top position where the carbon atom is placed over the Ta atom of the David star center and the hollow position where the center of graphene ring and the David star are placed on top of each other, see insets in figure 3. The CDW formation lowers the heterostructure energy per supercell by 152 meV for the top and by 136 meV for the hollow stacking, respectively. We note that those energies reasonably correspond to the transition of 1T-TaS₂ to the commensurate CDW phase [32].

In figure 2, we show calculated band structures of the supercell heterostructure for the top stacking along the high symmetry lines in the first Brillouin zone for the case where 1T-TaS₂ is in the normal and in the CDW state (including the magnetic state, which we discuss in section 5). The fat bands shown by the dark gray circles correspond to the states originating from the carbon atoms. The graphene gets hole-doped and its Dirac cone is shifted about 0.25 eV and 0.30 eV above the Fermi level for the normal and CDW phase, respectively. The π-bands of graphene hybridize with the 1T-TaS₂ states in the vicinity of the Fermi level. We note that in the case of the CDW phase, the bands of 1T-TaS₂ significantly change and the Dirac states hybridize near the Fermi level with the partially filled Mott band originating predominantly from the 5d orbitals of Ta atoms. The change in hybridization has a direct effect on proximity-induced SOC in graphene. The band structure for the hollow stacking (not shown) possesses on the energy scale as shown in figure 2 almost the same bands' dispersion as the top stacked configuration. The bands differ noticeably in the vicinity of the Dirac point for two stackings in the CDW phase as will be seen from discussion of figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The figure 3 shows details of the band structure topology of graphene in the vicinity of the Dirac point and corresponding s_z , s_y , and s_x spin expectation values owing to proximity effects for the normal and CDW phase of 1T-TaS₂. The calculated DFT results near the K point (k = 0 corresponds to the K point) are shown by circles while the solid lines are tight-binding model fits. Four different colors of these lines mark distinct bands in all the panels of the figure. In each case, both the electron and hole bands are spin split. For the normal phase at the top and hollow stacking, the band dispersions are electron–hole-like symmetric with significant band inversion and avoided crossings due to proximity-induced Rashba SOC. A similar band topology has been observed in the case of graphene on 1H-WSe₂ leading to quantum spin Hall effect and pseudohelical states predicted in graphene nanoribbons [120, 121]. We note that in the normal phase, the s_x and s_y spin expectation values are identical for pairs of particle-hole symmetric bands, so that they overlap in the plots. Also, s_x and s_y spin expectation values vanish exactly at the K point for all cases considered in figure 3 due to the inherent C₃ symmetry. The emergence of CDW in the underlying 1T-TaS₂ monolayer leads to the strong band topology modification due to a change in the SOC proximity-induced effects.

The complex features of the proximity-induced effects on graphene electronic structure can be studied by analyzing an effective tight-binding model. We used the Hamiltonian given by equation (2) and fit its parameters to the DFT data. The parameter values are summarized in table 1 and provide a very good reproduction of the electronic band structure close to the Dirac point as well as the corresponding spin expectation values; compare DFT data (circles) and the tight-binding model (solid lines) in figure 3. The effect of the CDW phase in the case of top stacking leads to the band topologies mostly controlled by the sign change of the staggered potential Δ and significant spin expectation values variation. In the case of hollow stacking with CDW, the bands' dispersions are similar to the case of bare graphene in a strong transverse electric field which induces $|\lambda_{\mathrm R}|\gg|\lambda_{\mathrm I}|$ [122] and enforces the degeneracy of a pair of bands at the K point. However, in this case, the bands are nearly degenerate and the physical mechanism is more complex due to the presence of significant staggered potential and the sublattice-asymmetry of the intrinsic SOC (see the parameter values in table 1). The most striking effect of the CDW phase is the enhancement of the Rashba angle $\phi_{\mathrm R}$ in both stackings. This notably cants the in-plane spin components towards the radial direction, see figure 4 showing the behavior of the in-plane spin components along the circular contour in wavevector space at the distance of 10⁻³ Å⁻¹ from the K point. It is noteworthy that each arrow color corresponds to an individual energy band as in figure 3 and each band has different energy along the selected contour in wavevector space. A somehow similar effect regarding the presence of a radial component of the Rashba field was predicted in twisted graphene/TMDC structures [93, 96, 104]. In the context of the enhancement of the Rashba angle by the emergence of CDW, let us emphasize the fundamental change in the band structure in the vicinity of Dirac point focusing on the bands originating from TaS₂ (compare figures 2(a) and (b)). In the spirit of the model presented in [92] (discussing the case of a twisted heterostructure of graphene with TMDC), the change of the Rashba angle is a consequence of a pronounced difference between the TMDC energy band structure in the normal phase and the CDW phase.

Zoom In Zoom Out Reset image size

Table 1. Tight-binding model parameters for graphene on 1T-TaS₂ determined by fitting the model given by Hamiltonian in equation (2) to DFT data.

	Top stacking		Hollow stacking
	Normal	CDW	Normal	CDW
E0 (meV)	250.08	297.13	249.57	302.24
t (eV)	2.547	2.516	2.542	2.537
Δ (meV)	1.90	−1.30	1.91	2.13
$\lambda_{\mathrm{I}}$ (meV)	−0.02	0.62	−0.02	−1.02
$\delta\lambda_{\mathrm{I}}$ (meV)	−2.52	−1.85	−2.52	−2.06
$\lambda_{\mathrm R}$ (meV)	1.05	1.12	1.05	0.90
$\phi_{\mathrm R}$ (deg)	17.36	43.07	17.57	28.15

The partially filled band at the Fermi level in the CDW phase of the 1T-TaS₂ monolayer is suspected to indicate non-vanishing magnetic polarization according to recent first-principles studies [27, 54]. In order to study the presence of spontaneous magnetic polarization in graphene/1T-TaS₂ monolayer, we have performed spin-polarized DFT+U calculations considering the CDW phase for 1T-TaS₂. The stability of the magnetically ordered state can be evaluated by computing the energy difference between the non-magnetic and the magnetic ground state. Our DFT+U calculations yielded values (per supercell) equal to 34.7 meV and 33.1 meV for top and hollow stacking, respectively. On the other hand, analogous calculations performed for 1T-TaS₂ monolayers with the graphene layer removed resulted in values of 29.8 meV and 28.2 meV, respectively. The results support the view that the creation of a heterostructure increases slightly the stability of the magnetic state in 1T-TaS₂. For the calculations we have considered a preferential in-plane magnetization along the x-direction, as we found that the magnetic anisotropy energy, i.e. the out-of-plane minus in-plane configuration energy (per supercell) amounts to 0.56 meV and 0.76 meV for top and hollow stacking, respectively, supporting the easy-plane anisotropy. The same type of anisotropy is predicted to emerge for a bare 1T-TaS₂ monolayer with magnetic moments centered on David stars [54].

The band structure resulting from DFT+U calculations with CDW phase in 1T-TaS₂ is shown in figure 2(c), where the fat bands originating from Ta atoms are marked with circles with color scale representing the degree of spin-polarization, linear color scale from red color for spin-up polarization, through gray to blue color for spin-down polarization. The main feature distinguishing this case from unpolarized CDW phase shown in figure 2(b) is pronounced spin splitting of the Hubbard flat Ta band located close to the Fermi level to a pair of strongly polarized bands. The remaining bands for the plotted energy range exhibit significantly lower degree of polarization.

The effect of the in-plane magnetization on the graphene bands in the vicinity of the $K/K^{^{\prime}}$ point is well described by the effective Zeeman-like Hamiltonian $H_{\mathrm ex}$, given by equation (3), representing the coupling of a spin with the sublattice-dependent in-plane magnetization, having the interaction strength $\Delta_{ A/B}$ and an angle $\phi_{ A/B}$ with respect to the considered x-direction. Using the Hamiltonian $H_\kappa+H_{\mathrm ex}$ we were able to fit relativistic non-collinear DFT+U data. The obtained parameters are given in table 2, while the comparison between the DFT+U and fitted data in figure 5 shows that the addition of the $H_{\mathrm ex}$ is able to capture the dominant behavior of the graphene bands in the vicinity of the $K/K^{^{\prime}}$ point when the magnetic polarization emerges. In both the top and the hollow stacking, the sublattice magnetization direction $\phi_{ A/B}$ mostly points towards the x-direction. It resembles the situation found in monolayer 1T-TaS₂, in which the Mott-like band that lies on the Fermi level strongly contributes to the x-component of magnetization [54]. The in-plane spin pattern, shown in figure 6, differs significantly from the one depicted in figure 4, for the non-magnetic case with CDW. The tilt angle shows pronounced angular dependence and, also, differs significantly for all four bands involved.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Tight-binding model parameters for graphene on 1T-TaS₂ in CDW phase and in-plane magnetization determined by fitting the model given by equations (2) and (3) to DFT+U data.

	Top stacking	Hollow stacking
E0 (meV)	336.50	336.16
t (eV)	2.516	2.537
Δ (meV)	1.81	−2.27
$\lambda_{\mathrm{I}}$ (meV)	1.22	−1.83
$\delta\lambda_{\mathrm{I}}$ (meV)	1.98	2.47
$\lambda_{\mathrm R}$ (meV)	1.65	1.07
$\phi_{\mathrm R}$ (deg)	56.90	30.93
$\Delta_{ A}$ (meV)	−2.33	−0.85
$\Delta_{ B}$ (meV)	−1.66	−0.74
φA (deg)	−9.31	−5.90
φB (deg)	−3.86	−20.37

Our fitting shows that the presence of magnetism influences the graphene bands not only qualitatively through the appearance of the effective Zeeman-like Hamiltonian, but also quantitatively through the change of the nonmagnetic parameters of graphene Hamiltonian (equation (2)). This is not a surprise, since the exchange parameters $\Delta_{ A/B}$ of the order of 1 meV correspond to the giant built-in magnetic field of 17 T (see a comparison between the nonmagnetic data in table 1 and magnetic data in table 2), thus justifying the observed interference between the proximity-induced SOC and proximity-induced-magnetism. The effective model Hamiltonian we developed can be used in future studies probing interplays of several proximity-induced interactions, for instance, investigating spin-to-charge inter-conversion phenomena [123] or novel anisotropic Rashba–Edelstein effects that can be tuned by a twist [98] or by the engineering of CDW. The temperature in such cases can be used as an easy knob for switching between different phases and for turning on/off spin–orbit and exchange interactions.

In the paper, we report a study of the electronic properties of van der Waals graphene/1T-TaS₂ heterostructure focusing on the proximity effects in the vicinity of the graphene Dirac point. The underlying band structures were computed by DFT calculations and then modeled analytically employing an effective tight-binding Hamiltonian derived from the symmetry considerations. The model parameters were determined by fitting the tight-binding to the DFT data. A crucial finding is the fundamental influence of CDW formation on the proximity effect, resulting in a profound change in the parameters, spin textures, and fine-band topologies. As a consequence, controlling the CDW presence in a heterostructure consisting of graphene on 1T-TaS₂ provides a promising tool to manipulate the electronic properties of graphene near the Dirac point, in particular, a possibility to tailor the radial Rashba SOC, which cants the in-plane spin projections. In addition to SOC proximity, the in-plane magnetic polarization of the CDW phase superimposes the SOC, and exchange proximity effects offering to control the proximity effects by temperature.

Our findings demonstrate that the emergence of CDW in 1T-TaS₂ profoundly modifies the proximity effects in the studied heterostructure with graphene. The effect bears resemblance to the physical situation in twisted heterostructures, where modification of the twist angle allows control over the proximity effects. It should be strongly emphasized that in our system this effect takes place without physical modification of the relative twist angle of 13.9^∘ between the graphene layer and 1T-TaS₂ layer, as enforced by the lattice commensurability for $\sqrt{13}\times\sqrt{13}$ reconstruction. Therefore, a specific sort of spintronics/twistronics without a twist may be realized based on CDW control in our heterostructure.

The proximity-induced SOC controlled with the twist angle in van der Waals graphene/TMDC heterostructure might find application in spin-to-charge inter-conversion [98, 124]. In particular, the modification of the in-plane spin texture with the twist angle enables tunability of the Rashba–Edelstein effect [98]. In this context, we might mention that the controllability of the proximity effects in our heterostructure offers a route to master also the spin-to-charge inter-conversion but without the need to tune physically the twist angle and with full reversibility.

K Sz acknowledges financial support provided by the National Science Centre (Poland) under Grant No. 2015/19/B/ST3/03142 and by University of Łódź under Grant No. 1/IDUB/DOS/2021. M M acknowledges financial support provided by the Ministry of Education, Science, and Technological Development of the Republic of Serbia and DAAD Research Grant 57552336. Also, this project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under the Programme SASPRO 2 COFUND Marie Sklodowska-Curie grant agreement No. 945478. D K acknowledges partial support from the project SUPERSPIN funded by Slovak Academy of Sciences within the program IMPULZ IM-2021-26, as well as from Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within Project-ID 314695032-SFB 1277, and the VEGA Grant No. 2/0156/22—QuaSiModo awarded by the Ministry of Education, Science, Research and Sport of the Slovak Republic. M G acknowledges financial support from Slovak Research and Development Agency provided under Contract No. APVV SK-CZ-RD-21-0114 and APVV 20-0425, and by the Ministry of Education, Science, Research and Sport of the Slovak Republic provided under Grant No. VEGA 1/0105/20 and Slovak Academy of Sciences project IMPULZ IM-2021-42 and project FLAG ERA JTC 2021 2DSOTECH. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre (www.lrz.de).

The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.