Introduction

Since the discovery of long-range magnetic ordering in atomically thin CrI3 and CrGeTe3 magnets1,2, there have been incredible efforts to unearth more compounds displaying 2D magnetism1,2,3,4. This discovery revitalized the discussion on the role of magnetic anisotropy in the stabilization of magnetism in 2D, which as recently demonstrated5, is not a requisite for magnetic ordering in low dimensions. As a result, interest in the van der Waals (vdW) class of materials has sky-rocketed3,4,6,7,8,9, as there has been a great success in uncovering ground-breaking findings in the fundamental physics of anisotropic low-dimensional magnetism.

2D materials are of particular interest because they create exciting opportunities for the realization of functional devices. It was found that bombarding graphene with protons induced ferromagnetism via single-atom defects, opening the door for opportunities in defect engineering of atomically thin vdW devices10,11,12. Since then, and after discovering a true 2D ferromagnet10,11, there have been achievements in tuning magnetism, spin-phonon coupling, exchange interactions, superconductivity, and charge and spin transport properties in vdW crystals and heterostructures stemming from proximity effects, high pressures, and mechanical defects5,13,14,15,16,17. While developments in ultrathin vdW magnets are exciting, the understanding of 2D limit characteristics stems from the bulk properties. Indeed, in Mn3Si2Te6 and sister compound Fe2.7GeTe2 (FGT), we have previously reported strong changes in the magnetization, spin-phonon coupling, and exchange interactions of vdW crystals after proton irradiation, photoexcitation, and pressure-proving that external perturbations create an exciting environment for advantageously tuning the properties of vdW materials as the monolayer limit is approached16,17,18,19,20.

In bulk, FGT is a soft, layered ferromagnet within a vdW crystal structure21,22,23. It is part of a family of ternary transition metal tellurides, which are well-known magnets for displaying 2D magnetism due to their strong structural stability after exfoliation. FGT displays a Curie temperature TC ~ 154 K, while being relatively high in comparison to the transition temperatures of similar compounds, is still far too low for any practical commercial application. Hence, efforts to increase the TC of FGT have included the intercalation of external molecules into its van der Waals gap. One report successfully observed a strong ferromagnetic ordering at high temperatures (>300 K) upon intercalation with sodium24. While exciting, that method created decomposition phases, effectively overpowering the intrinsic magnetism of FGT. Moreover, it did not provide results as transformative as that of tetrabutylammonium intercalation into Cr2Ge2Te625. In that report, there was an increase in TC of nearly 150 K, and a strengthening of the exchange mechanism25. Indeed, the literature proves that intercalation is a manageable and accessible method of altering the properties of layered materials25. Hence, we set out to intercalate FGT single crystals with tetrabutylammonium cations (TBA+) (Fig. 1a) and quantify its effects on the magnetic properties.

Fig. 1: Intercalation of FGT with TBA+ molecules and a comparison of the temperature-dependent magnetization.
figure 1

a Atomic scheme of FGT and TBA+ molecules intercalation. b Magnetization vs. Temperature curves for pristine (green) and intercalated (black) FGT. Zero-field-cooled (open symbols) and field-cooled (solid symbols) curves are also shown. The data were taken in the out-of-plane direction.

Full size image

In this work, we show that intercalated FGT exhibits room temperature ferromagnetism. By using temperature-dependent Raman measurements, in combination with vdW-corrected density functional theoretical calculations, we report TC values up to 350 K, which are robust against magnetic fields and temperature variations. The intercalation also modifies the intrinsic properties of the FGT crystals, including the spin-phonon coupling, the anisotropy field, and the uniaxial anisotropy coefficient. Our findings indicate that intercalated TBA+ molecules engineer the magnetic properties of FGT towards practical implementations.

Results and discussion

Sample preparation and characterization

As presented and discussed in the Supplementary Information (see Supplementary Fig. 1 and the discussion), from the discharge curve of (TBA)x FGT, the maximum intercalation amount of TBA+ in FGT is estimated as ~1.11. From the calculations, the amount of the intercalant was found to be 0.83, which is similar to the value estimated from the experimental method. As discussed in the Supplementary Information (see Supplementary Fig. 2 and the discussion), X-ray photoelectron spectroscopy (XPS) analysis of the blank substrate confirms -NH2 moieties as being present in the adhesive, along with a smaller fraction of more electronegative double/triple-bonded nitrogen, but ammonium is not detected. Conversely, spectra of the electrochemically intercalated FGT show positively charged nitrogen species, which we fit with a peak at ~401.8 eV, separate from the background amines, found at ~399.4 eV according to our spectral fitting. We ascribe this observation to our successful introduction of TBA+ molecules into the Fe3-xGeTe2 crystal matrix.

Temperature-dependent magnetization

Figure 1a shows the atomic scheme of TBA intercalated FGT. Figure 1b shows a comparison of the temperature-dependent magnetization of pristine and intercalated FGT. These measurements were conducted along the magnetic easy axis, where the field is applied perpendicular to the plane of the crystal, referred to as the out-of-plane direction (H//c). For the zero-field-cooled (ZFC) curves, the samples were cooled to 50 K in the absence of a magnetic field and allowed to thermally stabilize before applying a measuring field of 79.5 kA/m and sweeping the temperature back up to 300 K. For the field-cooled (FC) curves, a cooling field of 79.5 kA/m is applied as the sample is cooled and is maintained through the temperature sweep. In the case of ZFC measurements, the intercalated sample (denoted by black open circles) demonstrated a dramatic increase in magnetization when compared to the pristine crystal, possibly due to the intercalation-induced carrier doping (discussed below). While the precise origin for the difference in the ZFC and FC behavior of these compounds is unclear at this point, the effect of thermal-cycling history may be at play. Although the FC measurements between pristine and intercalated FGT (shown in solid green and black circles, respectively) did not show much difference, the change in behavior upon applying a cooling field is a clear indication that the magnetic response of intercalated FGT is altered upon intercalation.

The temperature-dependent magnetization data were used to estimate TC by taking the first derivative (see Supplementary Information Fig. 3). As shown in Fig. 1b, there is a negligible change in TC despite the changes observed in the magnetization. In fact, the values of the TC for both pristine and intercalated FGT agree well with the reported value of 153 K found in the literature21,22,23. In other words, the introduction of TBA+ did not affect the intrinsic TC of FGT, as is shown by the coinciding sharp drops in magnetization as the temperature is increased. Despite this, there is yet another striking difference in the temperature-dependent magnetization of intercalated FGT. In pristine FGT, after heating the sample through its TC, there is a weaker magnetization as thermal fluctuations overpower the energy provided by the applied magnetic field. However, this is not true of intercalated FGT, which shows a significant magnetization even up to 300 K. This is a striking indication of long-range magnetic ordering sustained up to 300 K, which is atypical of iron deficient FGT such as the pristine samples we present21,22,23. Noting an initial drop in magnetization at TC ~153 K, we attribute the high-temperature behavior of FGT to a secondary ferromagnetic phase induced by intercalation with TBA+ molecules.

Magnetic field-dependent magnetization

This assertion is reinforced by isothermal magnetization data taken at low temperatures below and through TC ~153 K up to 350 K. Shown in Fig. 2, these five-quadrant magnetization vs. magnetic field curves were taken at several temperature plots up to a maximum field of 2.38 MA/m (3 T) in both polarities. As expected of soft-ferromagnetic FGT, there is small coercivity even for the intercalated sample. In addition, the relatively low saturation field indicates that the easy axis is indeed along the stacking direction (the crystallographic c-axis). Furthermore, this shows that the easy axis is stable upon intercalation. The most striking comparison between pristine and intercalated FGT is present at room temperature, where a clear ferromagnetic hysteresis is observed at room temperature only for the intercalated sample. In addition, these data present a “wasp-waisted” hysteresis which may be an indication of competition between two phases with firmly distinct coercivities. This kind of hysteresis behavior is due to the negative (antiferromagnetic) exchange related to the interlayer antiferromagnetic coupling, as widely reported26. Indeed, this effect is only seen after intercalation, since pristine FGT is clearly paramagnetic at 300 K, as shown in Fig. 2a (shown in black). Further evidence of a secondary phase is shown in Fig. 2b, which presents a bifurcation of the first and fifth quadrant data from which a critical field of approximately 1.27 MA/m (1.6 T) is extracted. This feature is not present in the data for the pristine FGT sample. Most importantly, Fig. 2c presents the isothermal magnetization behavior of the intercalated FGT up to the temperature of 350 K. The strong hysteretic behavior confirms the ferromagnetic nature of this intercalated compound retains at least up to 350 K. The isothermal magnetization collected along the magnetic hard axis of intercalated FGT is shown in Supplementary Fig. S5.

Fig. 2: The isothermal magnetic field-dependent magnetization at various temperatures.
figure 2

a shows a comparison between fully paramagnetic pristine (black) and the room temperature ferromagnetic behavior of intercalated (red) FGT. b, c display the ferromagnetic ordering of intercalated FGT up to 350 K.

Full size image

Raman spectroscopy and theoretical calculations

To gain deeper insights, we recorded the temperature evolution of the Raman spectra for the two phonon modes \({E}_{2g}^{1}\) and \({A}_{1g}^{1}\), which are related to the magnetic ordering of FGT22,27. First, in agreement with the literature, the wavenumber of these two modes shifts to lower values after intercalation (see Supplementary Information Fig. 6, also see below)22,27. We propose that the effect is a result of not only weakened interlayer coupling stemming from increased separation distance of the FGT layers but also, most importantly, by electron doping caused by intercalation. Indeed, Fig. 3 includes the simulation results of the effect of TBA+ molecules on the charge rearrangement of FGT layers. A substantial number of electrons from TBA+ moved into the Fe and Te atoms, as calculated in the Bader charge analysis (Table 1). We also observed that during the relaxations, the TBA+ may stabilize in some geometrical positions that might not induce similar charge rearrangement in the FGT layers composing the supercell. That is, depending on how the atoms present in the TBA+ molecules (N, C, H) interact with the FGT atoms (Fe, Ge, and Te), the interactions locally may be different, which might induce slightly distinct charge density differences on the FGT layers (Fig. 3b). For instance, the far left FGT layer on Fig. 3a gained +0.608 electrons, whereas the far right FGT layer gained +0.622 electrons. Similar arguments apply to the TBA+ molecules. In addition, we hypothesize that the changes in magnetic behavior may correlate to alterations in the spin-phonon coupling of the intercalated sample. To investigate this possibility, Raman spectroscopy data in the parallel polarization configuration were collected at several temperatures and are shown in Supplemental Information Fig. 6a, b for pristine and intercalated FGT, respectively. The main phonon modes under investigation are the \({E}_{2g}^{1}\) and \({A}_{1g}^{1}\) modes with a Raman shift of approximately 125 and 142 cm−1, respectively, when exciting with 514.5 nm excitation. After modeling the spectra as a superposition of Lorentzian signals, the peak centers (in cm−1) for each mode are plotted versus temperature on a logarithmic scale in Fig. 4. There is a discontinuity in these values near TC ~153 K, which serves as an indication of spin-phonon coupling22,27. Hence, the points above the intrinsic transition temperature are fitted with a Boltzmann sigmoidal model to simulate the temperature-dependent behavior of the phonon modes as though they were unperturbed by the ordering of the spins. This process is outlined in the previous report on the sister compound CrSiTe328. The clear deviation of our experimental data with the simulation at low temperatures is additional evidence of spin-phonon coupling across both modes in pristine and intercalated FGT. This separation is quantified by taking the largest difference in the Raman shift (Δω) present at the lowest temperature (see Fig. 4). This parameter is used to estimate the spin-phonon coupling parameter in the system.

Fig. 3: Theoretical simulations of the charge density difference in intercalated FGT.
figure 3

a, b Charge density difference (\(\triangle {\rm{\rho }}={\rm{\rho }}\left({\rm{FGT}}+{\rm{TBA}}\right)-{\rm{\rho }}\left({\rm{FGT}}\right)-{\rm{\rho }}({\rm{TBA}})\), where \({\rm{\rho }}\left({\rm{FGT}}+{\rm{TBA}}\right)\) is the total charge density of FGT plus TBA, \({\rm{\rho }}\left({\rm{FGT}}\right)\) is the charge density of the isolated FGT system, and \({\rm{\rho }}({\rm{TBA}})\) is the charge density of the TBA molecules) through iso-surface and planar average along the z direction, respectively. In a, green and blue iso-surfaces correspond to positive and negative charge, respectively.

Full size image
Table 1 Values obtained from theoretical analysis of each atom in the system.
Full size table
Fig. 4: Peak centers for both phonon modes as a function of temperature.
figure 4

The temperature evolution of these modes for pristine (closed symbols, a, b) and intercalated (open symbols, c, d) FGT is modeled using a Boltzmann sigmoidal model (dashed line). The deviation between this simulated unperturbed mode and the experimental data is readily apparent.

Full size image

Equation (1) is routinely used to extract the spin-phonon coupling parameter (λ’)28. However, Eq. (2) is the shorthand notation of Eq. (1). Finally, λ’ is algebraically extracted using Eq. (3).

$${\omega }^{2}={\omega }_{0}^{2}+\lambda \left\langle {{\boldsymbol{S}}}_{i}\cdot {{\boldsymbol{S}}}_{j}\right\rangle$$
(1)
$$\omega \approx {\omega }_{0}+{\lambda }^{{\prime} }\left\langle {{\boldsymbol{S}}}_{i}\cdot {{\boldsymbol{S}}}_{j}\right\rangle$$
(2)
$${\lambda }^{{\prime} }=\frac{\triangle \omega }{\left\langle {{\boldsymbol{S}}}_{i}\cdot {{\boldsymbol{S}}}_{j}\right\rangle }$$
(3)

First, we take Δω as the value extracted from the deviation of the sigmoidal model and the experimental data as the numerator. Then, we assume that all the magnetic contribution is a result of the Fe3+ cation29, allowing us to set Si = Sj = 3/2 and make the denominator equal to 9/4. From these values, the spin-phonon coupling parameters for both the \({E}_{2g}^{1}\) and \({A}_{1g}^{1}\) phonon modes are calculated for both pristine and intercalated FGT. The values are shown in Table 2. Our analysis shows that the spin-phonon coupling for the \({E}_{2g}^{1}\) mode decreases after intercalation. On the other hand, the opposite is true for the \({A}_{1g}^{1}\) mode, as λ’ slightly increases after intercalation. Our analysis shows that there is a clear change in the spin-phonon coupling of FGT after intercalation. It is likely that this change is caused by alterations to the exchange interaction stemming from changes in the metal-ligand bond angles25. It is known that the phonon vibrations are strongly dependent on the spin ordering in magnetically active materials30,31. The breaking of magnetic symmetry changes the phonon spectrum upon intercalation. Interestingly, as recently reported by Weber et al. the simultaneous presence of two magnetically ordered phases can lead to the strengthening of spin-phonon coupling32. We expect a similar situation can occur here as well.

Table 2 Parameters obtained through Raman spectroscopy.
Full size table

To explain the presence of the room temperature ferromagnetic phase, we look at the role of the intercalant molecule in the overall structure of FGT. To reiterate, we do not expect any substitution of TBA+ into the iron vacancies present in the sample. Therefore, we expect all TBA+ to reside in the van der Waals gap, creating an increased interlayer separation as reflected in the Raman data. It was also found that electron doping has significant effects on the magnetic exchange within FGT24,33,34,35. In a previous report on TBA intercalation of Cr2Ge2Te6, the electrochemical reaction responsible for the intercalation of TBA into the crystal, and an increase in the saturation magnetization show evidence of electron doping in the sample25. Because we have intercalated our FGT with the same organic molecule, we assume that one full electron is added to FGT upon intercalation, as supported by our theoretical prediction (see Fig. 3). It was shown that high levels of electron doping cause a drastic increase in TC by suppressing the AFM arising from the geometric frustration of inequivalent Fe sites. Suppression of the AFM coupling between FGT layers results in the most dramatic changes to TC35. The temperature-dependent Raman data coupled with theoretical calculations point to the fact that there is an electron doping upon intercalation. This could potentially cause the appearance of high-temperature ferromagnetic phases associated with the enhancement in magnetization. Such carrier-tunable magnetization has been observed in many vdW magnets25. Another possible explanation lies in the electronic structure of intercalated FGT. Cr2Ge2Te6 saw a dramatic increase in TC because of the occupation of conducting dxz/dyz orbitals by the extra electrons introduced into the system upon intercalation25. This was reported due to a shift in exchange mechanism from a relatively weak superexchange across the ligand Te atom, to a robust direct exchange between Cr atoms25. The Stoner model readily describes the magnetic ordering of these systems, and the itinerant ferromagnetism of FGT is no exception21,29,36. Electron doping experiments on FGT have shown that the extra electrons begin to occupy the dxz/dyz orbital sub-bands, and alter the electronic density of FGT, similar to intercalated Cr2Ge2Te621,29,34. As the Stoner model implies, changes in the density of states near the Fermi level would cause considerable changes to the ferromagnetic ordering of FGT. Therefore, we expect that intercalation-induced electron doping would explain room temperature ferromagnetism by altering the electronic density of states of FGT in a way similar to intercalated Cr2Ge2Te6, where there was a shift in the exchange mechanism to a strong Fe-Fe direct exchange facilitated by additional exchange channels provided by intercalants34.

Magnetic anisotropy

Changes in the electronic density of states in Cr2Ge2Te6 resulted in the complete switching of the magnetic easy axis25. While we did not detect such drastic changes in FGT, there is evidence that the anisotropy of the intercalated sample is decreased, further cementing the idea that the density of states of FGT near the Fermi level is affected in a similar manner to that of Cr2Ge2Te6 after intercalation. An estimation of the anisotropy field and the uniaxial anisotropy constant of FGT reveal a decrease in both of these quantities after intercalation, which is detected as a decrease in the aforementioned anisotropy parameters. The results are shown in Fig. 5. The mechanism at play may be the contribution of a planar orbital moment component allowing the spins to align along the in-plane direction more easily25. Indeed, intercalated FGT shows a smaller saturation field along the magnetic hard axis (H//ab), meaning that less energy is needed to overcome the anisotropy barrier and align the moments along their preferred direction. The anisotropy constant Ku is related to the saturation magnetization MS and saturation field Hsat along the in-plane direction by the Stoner–Wolfarth model:

$$\frac{2{K}_{u}}{{M}_{s}}={\mu }_{0}{H}_{{sat}}$$
(4)

where µ0 is the vacuum permeability (unity in cgs units)21,29. The change in anisotropy is more readily correlated to a decrease in the saturation field, as the values of magnetization at high magnetic fields are comparable after intercalation. As a function of temperature, Ku is consistently lower in the intercalated sample throughout the ferromagnetic phase. Theoretical calculations on electron-doped FGT have also pointed to contributions from the dxz/dyz orbitals as mediators of the anisotropy energy while having only a minimal effect on the magnetic moment per Fe atom29. In our intercalated FGT, we observe an appreciable decrease in the anisotropy, but no drastic changes in the high-field magnetization. Therefore, we can expect that intercalated FGT sees a contribution from previously unoccupied states near the Fermi level and behaves similarly to Cr2Ge2Te6 after intercalation25. Thus, a shift in the exchange mechanism from a relatively weak superexchange to a stronger direct exchange may be what allows FGT to retain ferromagnetic ordering to such high temperatures. To note, we observed such intercalation-induced magnetic properties in similar compounds such as ferromagnetic metallic Cr2Ge2Te6 and ferrimagnetic insulating Mn3Si2Te6.

Fig. 5: Magnetic anisotropy coefficients for pristine and intercalated FGT.
figure 5

The anisotropy field (black, left y-axis) and the uniaxial anisotropy coefficient (blue, right y-axis) are shown for both pristine (solid symbols) and intercalated (open symbols) in the low-temperature region. It is shown that the anisotropy parameters of FGT decrease noticeably after intercalation.

Full size image

To summarize, tetrabutylammonium cations were electrochemically intercalated into the vdW gap of Fe2.7GeTe2 bulk crystals. The effects were investigated via low-temperature magnetometry and Raman spectroscopy measurements. Temperature-dependent magnetization data show a non-zero magnetization for intercalated FGT up to 350 K, confirmed by isothermal field-dependent magnetization data to be the result of an induced-ferromagnetic phase present far beyond the intrinsic TC ~154 K of pristine FGT. Despite this, a transition is still observed at the intrinsic Curie temperature, and there are no changes in the saturation magnetization in either crystallographic orientation at low temperatures. However, above room temperature, the intercalated FGT shows a clear ferromagnetic hysteresis above at least 350 K.

The temperature-dependent Raman measurements coupled with our density functional theoretical calculations suggest that charge transfer could be the cause for the room temperature ferromagnetism upon intercalation. These findings lend strength to the idea of electrochemical intercalation as a viable method of realizing room temperature ferromagnetism for the purpose of incorporation into practical devices.

Supplementary Information contains experimental methods, theoretical simulations, acknowledgements, and additional figures.

Methods

Experimental measurements

FGT crystals were grown through the self-flux technique, as outlined in previous reports, resulting in millimeter-sized single crystals22.

Electrochemical Intercalation of FGT crystals with TBA+

Initially, the FGT crystals were rinsed with ethanol many times. Then, single crystals were hard-pressed onto an indium plate and adjusted to work as positive electrodes, where a 0.1 mm silver piece was used as a negative electrode. A saturated solution of tetra butyl ammonium bromide powder (Sigma Aldrich, 99.0%) dissolved in DMF (Sigma Aldrich, 99.9%, was used for the electrochemical intercalation process. To perform homogenous intercalation, a potentiometric process was performed with a constant current of 2 μA. The whole intercalation process was carried out in an argon medium using an airtight cell. By applying the constant current, the negative electrode (Silver) will start to lose electrons and the FGT crystals will start to accept those electrons. The amount of intercalated material is then determined from the discharge curve and the specific criteria set by Wang et al.25.

X-ray diffraction data (see Supplementary Information Fig. 4) were taken for the pristine and intercalated samples. This XRD data confirm that the intercalant molecules indeed reside in the van der Waals gap of the crystal, and retains its crystal structure. This is of particular concern as this FGT sample is iron deficient, leaving a wide opportunity for molecules to readily substitute into any vacancies instead of intercalating into the van der Waals gap. In other words, we are confident that any changes to the magnetic behavior of our samples intercalated with TBA+ molecules are not due to the destruction of the intrinsic crystal structure of FGT. As reported, if the amount of intercalant exceeds an optimal value, then the structure of the crystal would be completely disrupted. Due to the sharp peaks and correspondence to previous literature reports, we can assert that intercalation did not affect the structure or crystallinity of FGT22,23,37.

Low-temperature bulk magnetometry measurements for intercalated FGT were conducted using a Quantum Design Physical Property Measurement System with the AC Measurement System option. DC moment data were extracted in a temperature range extending from 5 to 350 K and a maximum magnetic field of ±30 kOe. Orientation-dependent measurements were conducted for both the in-plane (H//ab) and out-of-plane (H//c) direction of the magnetic field with respect to the crystallographic c-axis of the crystal.

Raman spectroscopy measurements were performed in the parallel polarization configuration on the bulk crystals. A fresh surface was exfoliated for the measurements. To ensure a reduction of noise in the signals, the samples were exfoliated. Measurements were conducted in the temperature range of 3.2 to 295 K using 514.5 nm excitation.

X-ray photoelectron spectroscopy (XPS) was performed using a Kratos Axis Ultra system with a monochromatic Al Kα excitation source operating at 15 kV and 10 mA. An intercalated FGT crystal was mounted on tape and pumped to a base pressure of <2.0 × 10−9 Torr (2.7 × 10−7 Pa) before measurement. The detector was at a photoelectron take-off angle of 90° to the surface and a pass energy of 20 eV was used. Data analysis was conducted with Shirley’s background using CasaXPS.

vdW-corrected ab initio simulations

First-principles simulations were done using the Vienna Ab initio Simulation Package (VASP),38 with the PAW pseudo-potentials using the Perdew–Burke–Ernzerhof (PBE) exchange-correlation and an energy cutoff of 600 eV. We needed to use a 2x2x1 FGT supercell to intercalate the TBA. This led to a supercell containing 154 atoms. Such a large supercell to accommodate all atoms into the system became computationally impractical if more TBA+ molecules are considered (e.g., stacking of several TBA+ molecules between the FGT layers). That’s, the study of how the TBA+ molecules may organize themselves inside of the FGT layers is beyond the scope of this work. We used the DFT-D3 corrections39 to account for the van der Waals interactions between the FGT layers and TBA molecules. The relaxations were done using convergence criteria of 10−8 eV and a k-mesh of 1 × 1 × 1 at the Gamma point. The single-shot calculations were done using convergence criteria of 10−9 eV and a k-mesh of 1 × 1 × 1 at the Gamma point to compute the charge density. The single-shot calculations were done on the full system (TBA + FGT), and on two subsystems constructed by removing either the FGT or TBA from the full system so that we could compute the charge density difference via \(\triangle \rho =\rho \left({{\mathrm{FGT}}}+{{\mathrm{TBA}}}\right)-\rho \left({{\mathrm{FGT}}}\right)-\rho ({{\mathrm{TBA}}})\). We noticed during these simulations that both the interlayer distance and the charge transfer are strongly coupled to each other without any possibility to separate them. We also checked that in simulations where the FGT layers are kept apart beyond some threshold, full convergence is not achieved. This indicates that the interlayer distance and the charge transfer play together to stabilize the system.