\section{Introduction}
Two-dimensional (2D) transition metal dichalcogenides (TMD) have recently attracted significant attention due to their promising opto-electronic and mechanical properties and ultrathin nature \cite{Manzeli2017}. As in graphene, the weak van-der Waals (vdW) coupling between adjacent layers enables scalability from thick TMDs to mono- or few-layer crystals \cite{Manzeli2017}. This reduction in dimensionality, from bulk to monolayer (ML), has been shown to affect the material's band structure, driving bandgap transitions, and modifying the bandgap width \cite{Mak2010, Jin2013, Komsa2012, Trainer2017}. TMDs at the atomic scale are of interest for device scaling \cite{Radisavljevic2011, Baugher2013, Kim2012, Desai2016, Jun2014, Subhamoy2014}. In particular, ML to few-layer molybdenum disulfide (MoS$_2$)-based transistors have demonstrated impressive device performances, characterized by high mobility, low sub-threshold slope, and excellent on/off current ratios \cite{Radisavljevic2011, Ghatak14}. Moreover, as the band structure of TMD MLs is strongly affected by mechanical deformations\cite{Naumis2017, Roldan2015, Guinea2009, Blundo2021}, the performance of TMD-based devices can be further tuned by strain \cite{Datye2022}. Given their high in-plane stiffness and low flexural rigidity \cite{Akiwande2017}, individual atomic sheets of TMDs can sustain much larger mechanical strains compared to conventional semiconductors \cite{Changgu2088}, thus allowing wider and finer tuning of electronic and photonic properties \cite{DiGiorgioREV}. However, their atomic scale thickness makes their intrinsic properties, and related device performance, very sensitive to structural defects, such as atomic vacancies, dislocations, and grain boundaries \cite{Hu2018, Yin2019}, which are unavoidably present and sometimes further induced by strain production techniques. Some of these defects can behave as charge trapping states, introducing mid-gap or band edge states \cite{Qiu2013, Hong2015, Liu2013}. A Shockley-Read-Hall (SRH)-like exchange of electrons or holes, can occur between the defect state and the majority carrier band, affecting the capacitive response of the material to an external electric field \cite{Brammertz2008, Zhu2014}. This a priori uncontrollable \textit{defect-state capacitance} adds to the intrinsic \textit{quantum capacitance} ($c_Q$), which accounts for the partial screening of electric fields in 2D materials, due to their low density of states (DOS) \cite{Luryi1988, Bera2019, Ma2015, Fang2018}.

In this article, we investigate both the defect states and the quantum capacitance of biaxially strained MoS$_2$ MLs in the form of `bubbles': low-energy hydrogen (H)-ion irradiation induces (sub-)micrometric protrusions where a ML becomes detached from the bulk \cite{Tedeschi2019, He2019, DiGiorgio2020, DiGiorgio2021, Blundo2020}. These bubbles are filled with highly pressurized H$_2$ gas \cite{DiGiorgio2020, DiGiorgio2021, Blundo2020, Blundo2021, Blundo2021PRL, Blundo2022}. Recently, low-energy H-plasma treatments have been shown to induce sulfur vacancies at the topmost interface of the MoS$_2$ ML, with the Mo and S atoms in both the underneath planes left untouched \cite{Ghorbani2017, Lee2021}. In addition, the partial passivation of sulfur vacancies via H-intercalation has been demonstrated in MoS$_2$ MLs exposed to atomic H-gas \cite{Pierucci2017}. Indeed, Li et al. \cite{Li2016} have shown that the H adsorption free energy of a MoS$_2$ surface is affected by the density of sulfur vacancies, as well as by strain, decreasing as these increase. The presence and/or passivation of atomic-scale defects is thus expected to contribute to the electrostatic properties of the MoS$_2$ bubbles.

Here, we study the electrostatic properties of MoS$_2$ bubbles as well as bulk MoS$_2$ via radio frequency-assisted electrostatic force microscopy (RF-EFM). EFM is a non-contact atomic force microscopy (AFM)-based technique. It is nondestructive and has high spatial resolution ($\sim$few tens of nm). By studying the long-range electrostatic interaction between MoS$_2$ and conducting Pt probes, biased with (DC+RF) voltage, we distinguish the capacitive contribution of the ML defect states from that of the finite DOS (quantum capacitance). Indeed, we demonstrate that at low frequency (sub-MHz range), the capacitive behavior of the MoS$_2$ ML is governed by the response of defect states. At higher frequencies, the role of the intrinsic quantum capacitance dominates. Compared to existing macroscopic capacitance-voltage measurements \cite{Brammertz2008, Fang2018}, which provide information on defect states and quantum capacitance averaged over the whole ML-based device, the method employed here allows for their local, sub-micrometric investigation. Moreover, while diverse scanning probe spectroscopies have been used in the past for the investigation of the quantum capacitance in 2D materials \cite{Giannazzo2009, Tino2013, Li2022}, we present here direct imaging of $c_Q$, over the bubble surface. Finally, from the quantum capacitance we derive the Fermi energy $E_F$. The knowledge of $E_F$, together with its spatial dependence, gives the local charge carrier density. Through our analysis, we  elucidate the correlation between defects, band structure and carrier density, over the ML-thick membrane.

The use of RF-EFM for finite frequency capacitance mapping opens up possibilities for investigating other phenomena. The nanoscale electron compressibility in quantum materials, the evolution of conductance domains at the metal-insulator transition \cite{Liu2020}, the defect states of perovskites-based harvesting systems \cite{Futscher2020, Yun2022}, the frequency variations of the relative dielectric constant in polymeric insulators (commonly implemented in metal-insulator-semiconductor structures \cite{Estrada2013}), are just a few examples where this technique can bring new insights. Moreover, as it is nondestructive both for the tip and for the sample surface, EFM is highly reproducible and reliable.

\section{Results and discussion}

Figure \ref{fig:asimmetry}(a) shows a representative ($\mathrm{13\ \mu m \times 9.2\ \mu m}$) AFM topography of the surface of an irradiated MoS$_2$ bulk crystal. Here, numerous H$_2$-filled bubbles can be seen. As discussed in Refs. \cite{Tedeschi2019, DiGiorgio2020, DiGiorgio2021, Blundo2020, Blundo2021, Blundo2021PRL} and in Supporting Information 1, the bubbles have an universal (\textit{i.e.}, not dependent on their size) ratio of maximum height $h_{max}$ to radius $R$, which determines their strain profile \cite{Khestanova2016}. Furthermore, the bubbles are separated by flat regions where the bulk crystal has remained intact after irradiation. Thus, both bulk MoS$_2$ and bubbles can be studied in the same sample. As shown in the inset of figure \ref{fig:asimmetry}(b), EFM is a non-contact imaging technique based on the electrostatic interaction between the sample and a conducting AFM probe, biased with a (DC+RF) voltage, $V=(V_{DC}-V_0)+V_{RF}cos(2 \pi f_{RF}t)$, with $f_{RF}$ the frequency of the applied RF voltage. The effective voltage drop $V$ is shifted from the applied $V_{DC}$ by $V_0$, the \textit{contact potential difference}. The measurement of $V_0$ gives the mismatch between sample and tip work functions, $W_S$ and $W_T$ respectively, being $-{eV}_0=W_T-W_S$ \cite{Melitz2011}. To investigate the MoS$_2$ electric field response, we determine the tip-sample electrostatic interaction $F_{el}$ by measuring the frequency of an oscillating AFM cantilever, $df=f-f_0$, with $f_0$ and $f$ being the free resonance frequency and its interaction-driven shifted value. $df$ is related to $F_{el}$ through $df=-\frac{f_0}{2k}\frac{\partial F_{el}}{\partial z}=-\frac{f_0}{4k}\frac{{\partial }^2}{{\partial z}^2}\left(c_T{V_T}^2\right)$. Here, $c_T$ and $V_T$ are the tip-surface capacitance and voltage drop, respectively, and $k$ is the cantilever spring constant. If the sample is a metal or a highly conducting semiconductor (very short penetration depth, or Debye length, $L_D\simeq 0$), all of the applied voltage drops between the tip and the surface and $V_T=V$ is $z$-independent. In contrast, when the sample conducts poorly, $V_T=V-V_S=\frac{c_S}{c_S+c_T}V$, with $V_S$ and $c_S$ being the sample to ground voltage drop and capacitance, respectively. In both cases $df$ is quadratically dependent on $V$, as can be seen in Supporting Information 2, Eqs. (\ref{asymDF}).

\begin{figure}
    \centering
    \includegraphics[width=16cm]{Figure1}
    \caption{(a) AFM topography, $\mathrm{13\ \mu m \times 9.2\ \mu m}$ in lateral size, of an irradiated MoS$_2$ bulk crystal. Scale bar is $ \mathrm{2\ \mu m}$. (b) Main: $df(V_{DC})$ spectroscopy of MoS$_2$ bulk (black scatters) and strained ML (blue scatters). The red solid line is the fit to Eq.(\ref{asymDF}a). Inset: scheme of RF-EFM. (c) Main: cantilever resonance curves in presence of tip-sample electrostatic interaction at $V_{DC}$=+3,-6 V. Black dashed lines indicate the cantilever free resonance frequency $f_0$=7.8 kHz. Inset: rescaled resonance curves; red solid lines are fit to Eqs. (\ref{amplitude}). (d) Main: cantilever oscillation amplitude vs supplied bias $V_{DC}$. Inset: circuital configuration equivalent to the EFM experiment on MoS$_2$ bubbles. }
    \label{fig:asimmetry}
\end{figure}

Figure \ref{fig:asimmetry}(b) shows $df(V_{DC})$ spectra of an irradiated MoS$_2$ sample in two different locations: the bulk (black dots) and a bubble (blue dots). As shown by the parabolic fit (red line), the bulk exhibits the expected quadratic voltage dependence of $df$ on $V_{DC}$. In contrast, the $df(V_{DC})$ of the bubble is asymmetric: at positive $V_{DC}$, the $df$ of bulk and bubble fully overlap, while, at negative $V_{DC}$, the $df$ of the strained ML is always lower than that of the bulk. This deviation from $V^2$ is explained below.

In addition to $df$, the amplitude of cantilever oscillations $A$ is also asymmetric in $V_{DC}$. Figure \ref{fig:asimmetry}(c) shows cantilever resonance curves acquired at $V_{DC}=+3$ V and $-6$ V, in a bulk region and on a bubble, in black and blue, respectively, while sweeping the excitation frequency $f$. The vertical black dashed lines indicate the cantilever free resonance frequency, $f_0$= 7.8 kHz. At $+3$ V: i) the bulk and bubble traces fully overlap, having the same $A\simeq 10$ nm and $df=-10$ Hz; at $-6$ V: i) $df_{bubble}$ is smaller (in modulus) than $df_{bulk}$, $-28$ and $-35$ Hz, respectively; ii) $A_{bubble}$ is smaller than $A_{bulk}$, $\sim 7$ and $\sim 10$ nm, respectively. In addition, figure \ref{fig:asimmetry}(d) demonstrates that, for any value of $V_{DC}<-1$ V, $A_{bubble}$ is affected by a significant reduction, compared to $A_{bulk}$. In contrast, for all positive $V_{DC}$, they are practically identical.

We modeled the observed phenomena by assuming that, while the MoS$_2$ bulk fully screens the electric field in the entire bias voltage range, the bubble is capable of full screening only for positive $V$ ($L_D\simeq 0$) and has a non-negligible capacitive contribution at negative $V$ (finite $L_D$). Thus, the bubble's capacitance, $c_B$, is voltage-sign dependent, having a finite value only for $V<0$. As $c_B$ is in series with the tip-surface capacitance $c_T$ (as shown in the inset of figure \ref{fig:asimmetry}(d)), this leads to voltage-sign dependence of $V_T$, $F_{el}$, and $df$. Thus, the electrostatic force on the negative bias branch, $F_{el,neg}$, is weaker than that on the positive branch, $F_{el,pos}$ (Supporting Information 3):
\begin{footnotesize}
\begin{equation} \label{elecForce}
F_{el, neg}\simeq \frac{1}{2}(V_{DC}-V_0)^2\frac{\partial c_T}{\partial z}\left[1-\frac{c_T\left(z_0\right)}{c_B}-\frac{4}{c_B}\left.\frac{\partial c_T}{\partial z}\right|_{z=z_0} \delta z\left(f,t\right)\right]=F_{el,pos}+F_{stat}\left(c_B\right)+F_{dyn}\left(c_B,\ f\right)
\end{equation}
\end{footnotesize}

Here $\delta z(f,t)=A(f)cos(2 \pi ft)$ is the displacement of the cantilever from its equilibrium position, while $f$ and $A(f)$ are the excitation frequency and the Lorentzian-like oscillation amplitude function. Thus, $F_{el,neg}$ differs from $F_{el, pos}$ by two additional, static and dynamic, negative, contributions -- $F_{stat}(c_B)$ and $F_{dyn}(c_B,f)$ respectively. The reduction of $F_{el,neg}$ compared to $F_{el,pos}$ explains the observed reduction in $df$ when $V<0$: given $df=-\frac{f_0}{2k}\frac{\partial F_{el}}{\partial z}$, $F_{el,neg}<F_{el,pos}$ leads to $df_{neg}<df_{pos}$ (figure \ref{fig:asimmetry}(b) and (c)).

Further, $F_{dyn}(c_B,f)$ explains the observed reduction in amplitude $A$ (figure \ref{fig:asimmetry}(c)). Indeed, modelling the AFM cantilever as a harmonic oscillator one can solve the equations of motion in two cases: i) a single driving force (external excitation due to the piezoelectric actuator, $F_{ex}$), as for bulk MoS$_2$ at all $V$ and bubbles at $V>0$; ii) two dynamic forces, $F_{ex}$ and $F_{dyn}(c_B,f)<0$, as for bubbles at $V<0$. One finds that (see Supporting Information 4), with the same damping coefficient and driving force, the oscillation amplitude at resonance is reduced in case (ii). The inset of figure 1(c) shows measured resonance curves (dots), offset along the frequency axis to match their resonance frequencies. The solid red lines are fits to Eq. (\ref{amplitude}a) -- case (i) above, for the MoS$_2$ bulk, at $+3$ V and $-6$ V, and the bubble at $+3$ V -- and to Eq. (\ref{amplitude}b) -- case (ii) above, for the bubble at $-6$ V. The fit to the bubble gives $c_B=4.5$ fF (Supporting Information 4).

The proposed capacitive contribution of the strained MoS$_2$ ML should not be a surprise. As mentioned above, the capacitance of low-dimensional semiconductors (as MoS$_2$ ML) has two main components: the \textit{quantum capacitance} ($c_Q$), and the \textit{defect-state capacitance} ($c_{SRH}$). Both are DOS-related: the former accounts for the intrinsically low DOS, which causes imperfect screening of the electric field in the material; the latter is due to the presence of charge trapping defects, possibly coming from sulfur vacancies, dislocations, and grain boundaries \cite{Qiu2013, Hong2015, Liu2013}, which give rise to states in the bandgap. The exchange of charge carriers, throughout SRH-like transitions, between the majority carrier band (conduction band in the case of MoS$_2$) and these additional states, affect the electrostatic response of TMDs. Furthermore, this exchange mechanism is frequency dependent, with SRH transitions occurring on a characteristic time scale $\tau \propto e^{{\Delta E}/{k_BT}}$, related to their position in the bandgap. Here $\Delta E$ is the energy difference between the majority carrier band edge and the defect state, $k_B$ the Boltzmann constant, and $T$ the temperature \cite{Brammertz2008, Zhu2014}. Multiple SRH defect states can occur at different $\Delta E$, each contributing to the carrier exchange on a different time scale \cite{Zhu2014}.

The H+ bombardment used in our sample preparation might create defect states, in particular sulfur vacancies. We can thus describe the bubble as a circuit with three elements in parallel (inset of figure \ref{fig:parabola}(a)-left): $c_Q$, and the frequency ($f_{RF}$)-dependent SRH capacitive and resistive branches ($c_{SRH}(f_{RF})=c^*_{SRH}\frac{arctg(2\pi f_{RF}\tau)}{2\pi f_{RF}\tau}$ and $R_{SRH}(f_{RF})=\frac{1}{c^*_{SRH}}\frac{\tau}{ln[1+(2\pi f_{RF})^2 \tau^2]}$, respectively, where $c^*_{SRH}=\frac{e^2 D_{SRH}}{2}$, with $D_{SRH}$ being the defect density) \cite{Zhu2014, Fang2018, Lehovec1966}. Multiple defect states can be represented as additional SRH capacitive and resistive branches in parallel to $c_Q$ \cite{Zhu2014, Lehovec1966}. The $f_{RF}$-dependent voltage drop at the tip-bubble surface is then given by:

\begin{equation} \label{VRF}
V_T(f_{RF})=\frac{1+i 2 \pi f_{RF}\left[c_Q+c_{SRH}\left(f_{RF}\right)\right]R_{SRH}\left(f_{RF}\right)}{1+i 2\pi f_{RF}\left[c_Q+c_{SRH}\left(f_{RF}\right)+c_T\right]R_{SRH}\left(f_{RF}\right)}V
\end{equation}

In the DC limit, $V_T(f_{RF}\rightarrow 0)=\frac{c_Q+c^*_{SRH}}{c_Q+c^*_{SRH}+c_T}V=\frac{c_B}{c_B+c_T}V$, with $c_B=c_Q+c^*_{SRH}$, which is precisely what we assumed initially (see figure \ref{fig:asimmetry}(d)). Here, the pre-factor $\frac{c_B}{c_B+c_T}$ now includes the asymmetric behavior of $c_Q(V)$ \cite{XU2020}. As MoS$_2$ MLs are intrinsically n-type, a small positive applied $V$ suffices to push the Fermi level $E_F$ of bubbles into the conduction band (CB), with $c_Q$ eventually saturating to its degenerate value ($c_{Q,deg}=e^2g_{2D}\approx 84 \mathrm{\frac{\mu F}{{cm}^2}}$, with $e$ and $g_{2D}=\frac{g_sg_vm^*}{2 \pi \hbar^2}$ the electron charge and the band-edge DOS, respectively \cite{Luryi1988, Bera2019, Ma2015, Fang2018}, $g_s$ and $g_v$ the spin and valley degeneracy factors ($\sim 2$ each \cite{Bera2019, Ma2015, Fang2018}), and $m^*\sim 0.4 m_0$ \cite{Phuc2018} the band-edge effective mass). In contrast, $c_Q$ remains negligible at all negative $V$ in the MoS$_2$ bandgap, reaching a finite value only when $E_F$ intercepts the valence band (VB) edge -- this is not achievable in our experiment and thus not addressed here. Thus, in the $V_{DC}$ range we use, we have:

\begin{table}
  \caption{Voltage sign dependence of $c_Q$, $c_B$, and $V_T$.}
  \label{tbl:signDep}
  \begin{tabular}{lllll}
    \hline
   $f_{RF}$ & $V_{DC}$ & $c_Q$ & $c_B$ & $V_T$  \\ \hline
    0 & Positive & $c_{Q,deg}>>c_T$\textsuperscript{\emph{a}} & $c_Q+c^*_{SRH}$ & $\simeq V$   \\
    0 & Negative & $\simeq 0$ & $\simeq c^*_{SRH}$, $c_T<<c^*_{SRH}$\textsuperscript{\emph{a}} & $\simeq V\frac{c^*_{SRH}}{c_T+c^*_{SRH}}\simeq V\left(1-\frac{c_T}{c^*_{SRH}}\right)$ \\
     \hline
  \end{tabular}

  \textsuperscript{\emph{a}} Supporting Information 5

\end{table}

Combining the results of Table 1 with the equation governing the dependence of $df$ on the electrostatic force (Eqs. (\ref{asymDF})), we get:

\begin{subequations} \label{dfasimmetry}
    \begin{empheq}[left={df=-\frac{f_0}{4k}(V_{DC}-V_0)^2\frac{\partial^2c_{T,eq}}{\partial z^2}=\empheqlbrace\,}]{align}
      & c_{T,eq}=c_T, \ \ \ V>0 \ \ \ (L_D\simeq 0) \\
      & c_{T,eq}=c_T{\left(1-\frac{c_T}{c^*_{SRH}}\right)}^2, \ \ \ V<0 \ \ \ (finite\ L_D)
    \end{empheq}
\end{subequations}

We thus expect the MoS$_2$ ML to behave as a metal for positive bias (Eq. (\ref{dfasimmetry}a) reproduces Eq. (\ref{asymDF}a)), and to have a finite capacitance for negative $V$, hence the observed asymmetry of $df(V_{DC})$.

\begin{figure}
    \centering
    \includegraphics[width=16cm]{Figure2}
    \caption{(a) Main: $df(V_{DC})$ spectroscopy on the top-most location of a MoS$_2$ bubble (left), MoS$_2$ bulk (center), and Au (right), as a function of the tip-surface separation, from 100 to 330 nm. The red solid lines are fits to Eqs. \ref{dfasimmetry}. The orange dashed line in left-hand panel is a representative parabolic fit. Inset: equivalent circuital configurations of the EFM experiments. Evolution of the contact potential difference $V_0$ (b) and defect states capacitance $c^*_{SRH}$ (c) as a function of the tip-surface separation. The red dashed lines work as a guide for the eye to follow the trend. Error bars are smaller than dot size.}
    \label{fig:parabola}
\end{figure}

The leftmost panel of figure \ref{fig:parabola}(a) shows $df(V_{DC})$ spectra acquired at the center of an MoS$_2$ bubble ($\mathrm{\sim 3\ \mu m}$ in diameter and $\sim$ 125 nm in height) as a function of the tip-surface separation, varying from 100 to 330 nm. The solid red lines are fits to Eqs. (\ref{dfasimmetry}a) and (\ref{dfasimmetry}b), for positive and negative $V_{DC}$, respectively. A parabolic fit to the positive voltage branch of the spectrum acquired at 100 nm (orange dashed line) underlines the asymmetry of the data. For comparison, we also acquired $df(V_{DC})$ spectra on a bulk region of the same MoS$_2$ sample and on an Au crystal (middle and right-hand panels of figure \ref{fig:parabola}(a), respectively). The corresponding circuits appear as insets. In the MoS$_2$ bulk, the electric field is screened over the Debye length $L_D$, giving rise to a capacitive contribution ($c_{dep}=\frac{\varepsilon_0 \varepsilon_R S}{L_D}$, with $S$ being the involved surface area, assumed as $\pi R^2_{tip}$). Here, $c_{dep}\gg c_T$, leading to $V_T\simeq V$ (Supporting Information 5). On the other hand, Au screens the electric field perfectly, and $V_T=V$. In conclusion, the tip-sample voltage drop, $V_T = V - V_S$, is equal to the applied bias, $V$, for both bulk MoS$_2$ and Au, so that Eq. (\ref{dfasimmetry}a) can be used to fit their $df(V_{DC})$ spectra in the entire $V$ range (red solid lines). Simulations of the electrostatic force acting between the AFM tip and both MoS$_2$ bulk and Au, based on their charge carrier densities, are reported in the Experimental section (see figure \ref{fig:simulation}).

Notice that the $df(V_{DC})$ spectra of figure \ref{fig:parabola}(a) have maxima which are are shifted by few-mV from $V_{DC}=0$. This can be explained by a non-zero contact potential difference $V_0$ between tip and sample, with $-eV_0=W_T-W_S$ \cite{Melitz2011}. The values of $V_0$ extracted from the fits as a function of distance are shown in figure \ref{fig:parabola}(b) as blue, black, and orange scatters for a bubble, bulk MoS$_2$, and Au, respectively. While the $V_0$ of Au and bulk are scattered around constant values ($-250\pm 3$ mV and $40\pm 9$ mV respectively), the bubble's $V_0$ decreases with increasing tip-sample separation. This can be explained by a sampling area which increases with the tip-sample distance. Indeed, as the distance increases, the electric field lines exiting the tip and entering the sample (or vice versa, depending on $V$ polarity) intercept a larger area, over which the electrostatic response is averaged. As expected, this does not affect the $V_0$ of Au and MoS$_2$ bulk. However, it does influence the measurement of $V_0$ on the bubble: averaging over regions with decreasing strain (moving from the top toward the edge of a bubble) leads to a decrease of $V_0$, which gradually reaches its bulk value. We see that the $V_0$ measured on the MoS$_2$ bubble is always higher than that of the bulk, given the same tip material (and thus the same $W_T$), in agreement with previous reports \cite{Tosun2015}. Moreover, the $V_0$ values measured on Au and MoS$_2$ bulk give the difference in their work functions: $W_{Au}-W_{MoS_2,bulk}=290 \pm 12$ meV. This is in agreement with $W_{Au}=5.4$ eV \cite{Uda1988} and $W_{{MoS}_2,bulk}\simeq 5.1$ eV \cite{Choi2014} from previous works. Note that $W_{MoS_2}$ can vary from $\sim 4$ (n-type) to $\sim 6$ eV (p-type), depending on doping \cite{Schauble2020}. Additionally, we mapped $V_0$ across different bubbles with Kelvin probe force microscopy, and found consistent behavior across them (Supporting Information 6). DFT calculations of the MoS$_2$ DOS suggests that the minimum of the CB, $E_c$, increases with reduced material dimensionality \cite{Trainer2017}, and decreases with tensile strain \cite{Hsu2017}. Increasing tensile strain also moves the vacuum level towards lower energies \cite{Postorino2020}. Taking all this into account, our results point to a location-dependent $E_F$, whose position inside the bandgap reflects a delicate equilibrium between strain-induced band structure changes and local doping.

Similarly to $V_0$, the asymmetry in $df(V_{DC})$ spectra is also separation-dependent, and gradually decreases as the tip moves farther from the surface of the bubble. As we explain the asymmetric behavior by the capacitance of ML MoS$_2$, its reduction means an increase in $c^*_{SRH}$ (see Eq. (\ref{dfasimmetry}b)). As shown in Figure \ref{fig:parabola}(c), when the separation increases from 100 to 330 nm, $c^*_{SRH}$ evolves from $3.9$ to $13.7$ fF. Given $c^*_{SRH}=\frac{e^2 D_{SRH}}{2}$, this could be due to an inhomogenous defect density $D_{SRH}$: minimal at the center of the bubble and increasing towards the edge.

To visualize the defect distribution across the bubble surface, we measured the location-dependent $df$ at a fixed separation of 153 nm and fixed negative $V_{DC} = -7$ V (figure \ref{fig:defectmap}(a)).
Figure \ref{fig:defectmap}(b) shows a representative $df(V_{DC})$ spectrum acquired at the center of the bubble shown in the inset. It is fitted to Eqs. (\ref{dfasimmetry}) (red line). From $df$, we can obtain the defect density by calculating $c^*_{SRH}$ for each pixel from Eq. (\ref{dfasimmetry}b), then converting to $D_{SRH}$. The defect distribution is strongly inhomogeneous and increases monotonically from the center of the bubble to the edge, consistent with our previous observation (color scale of figure \ref{fig:defectmap}(a)). A second bubble, partially visible in the bottom-right corner of the topography, can also be seen in the $D_{SRH}$ map.

\begin{figure}
    \centering
    \includegraphics[width=16cm]{Figure3}
    \caption{(a) Map of $df/D_{SRH}$, acquired at $V_{DC}$=-7 V, and $z$=153 nm. (b) Main: $df(V_{DC})$ spectroscopy measured in the center of a MoS$_2$ bubble. The red solid line is the fit to Eqs. \ref{dfasimmetry}. Inset: AFM topography of the selected bubble. Scale bars are 1 $\mathrm{\mu m}$.  }
    \label{fig:defectmap}
\end{figure}

Finally, we studied the frequency-dependent electrostatic response of bubbles to an RF bias, $V_{RF}cos(2\pi f_{RF}t)$, with $f_{RF} = 0.1 \div 300$ MHz. By combining Eqs. (\ref{asymDF}) and (\ref{VRF}), given $V=(V_{DC}-V_0)+V_{RF}cos(2\pi f_{RF}t)$, in the limit of  $V_{DC}=V_0$, and $\left\langle {\left[V_{RF}cos(2\pi f_{RF}t)\right]}^2\right\rangle =\frac{V_{RF}^2}{2}$, we get:

\begin{scriptsize}
\begin{subequations} \label{dfRF}
    \begin{empheq}[left={df_{RF}=-\frac{f_0}{8k}V_{RF}^2\frac{\partial^2c_{T,eq}}{\partial z^2}=\empheqlbrace\,}]{align}
      & c_{T,eq}=c_T \ \ \ \ \  L_D \simeq 0 \\
      & c_{T,eq}=c_T{\left(\frac{1+i2\pi f_{RF}\left[c_Q+c_{SRH}(f_{RF})\right]R_{SRH}(f_{RF})}{1+i2\pi f_{RF}\left[c_Q+c_{SRH}(f_{RF})+c_T\right]R_{SRH}(f_{RF})}\right)}^2 \ \ L_D \ge L^*_{D}
    \end{empheq}
\end{subequations}
\end{scriptsize}

Here $L^*_D$ is the distance between the plates of an equivalent plane capacitor, $c^*_T=\frac{\varepsilon_0 \varepsilon_r S}{L^*_D}$, such that $c^*_T$ is comparable to $c_T$. We see that, when the electric field is imperfectly screened, $df_{RF}$ becomes $f_{RF}$-dependent -- this is the case for bubbles at negative $V$.

\begin{figure}
    \centering
    \includegraphics[width=16cm]{Figure4}
    \caption{(a) Main: frequency shift ratio $R=\frac{df_{RF,bulk}}{df_{RF,bubble}}$ fitted, in red, to the ratio between Eq. (\ref{dfRF}a) and (\ref{dfRF}b). Different inner-scatter colors refer to two sets of data taken with smaller steps on a smaller frequency range (yellow), and with larger step in a larger frequency range (black). Inset: map of $R(x,y)$ at $f_{RF}$=300 MHz. Scale bar is $\mathrm{1 \mu m}$. (b) Black: theoretical plot of $c_Q(E,0)$ for an intrinsic semiconductor having $E_g$=1.65 eV, at $T$=300 K. Red: theoretical plot of $c_Q(E,E^*)$ for an n-doped semiconductor having $E_g$=1.65 eV, $E_C-E_F \simeq$ 225 meV, $E^* \simeq$ 0.6 eV. Blue: experimentally measured value of $c_Q(0\ V)$. CB and VB appear at positive and negative energy, respectively. The yellow region indicates $E_g$. The $E^*$ arrow points at the shift between the two curves. Black and red $E_C$ indicate the minimum of CB for intrinsic and n-doped semiconductor.}
    \label{fig:dfratio}
\end{figure}

Figure \ref{fig:dfratio}(a) shows the ratio between the measured $df_{RF,bulk} (f_{RF})$ and $df_{RF,bubbles} (f_{RF})$,  $R=\frac{df_{RF,bulk}}{df_{RF,bubble}}$. See Supporting Information 7 for individual $df_{RF,bulk} (f_{RF})$ and $df_{RF,bubbles}$ traces. The fit, in red, is to the analytical expression obtained by dividing Eq. \ref{dfRF}(a) by Eq. \ref{dfRF}(b). The resulting fit parameters, $c_Q$, $c_{SRH}(D_{SRH},\tau)$ and $R_{SRH}(D_{SRH},\tau)$, are shown in the figure. A good fit is achieved when assuming two types of defect states (A and B) at different energies in the bandgap, and with different emission times $\tau$, being 1 and 0.1 $\mathrm{\mu s}$, respectively (the fit to a single defect state is shown in Supporting Information 8). Thus, both lie in the MHz and sub-MHz range, where they dominate the electrostatic response of the material.

At higher frequencies, $R$ is instead mostly governed by $c_Q(E,E^*)=e^2g_{2D}\left[1+\frac{exp(E_g/{k_BT})}{2cosh((E+E^*)/{k_BT})}\right]^{-1}$ \cite{Ma2015, Fang2018}. The $c_Q(E,E^*)$ of intrinsic and extrinsic semiconductors differ by $E^*$, which is zero for the former ($E_F$ in the middle of the bandgap, $E_g$), and finite for the latter. The black curve in figure \ref{fig:dfratio}(b) is the theoretical $c_Q(E,0)$ for a 0.8\% strained MoS$_2$ ML, which has $E_g=1.65$  eV \cite{Lopez2016}. Here $E_F \equiv 0$, $E_g$ is the width of the shaded yellow region and $E_c$, the bottom of the conduction band (black arrow). Thermal broadening at the band edges is also observable, due to the input temperature of $T=300$ K. Overlaid on the simulated quantum capacitance, the $c_Q$ value measured in our experiment at $V=0$ V, and extracted from the fit of $R$, is plotted as the horizontal blue dashed line ($c_Q(0\ V)=7.7\times 10^{-3}\ \frac{\mu F}{cm^2}$). As no states are expected at the Fermi level in intrinsic semiconductors, the blue and black curves do not intercept each other at zero energy. The red curve is, instead, the plot of $c_Q(E,E^*)$, with $E^*=\frac{E_g}{2}-(E_C - E_F)$ measuring the shift of $c_Q$ toward the VB. It thus models the quantum capacitance expected in an n-type semiconductor, having $E_C -E_F = 225$ meV. This latter intercepts the measured value of $c_Q$ at zero energy. This match allows the experimental determination of the $E_F$ position inside the bandgap of the strained MoS$_2$ ML.

We determined $c_Q$ over the surface of a bubble by measuring $R(x,y)$ at  $f_{RF}=300$ MHz (inset of figure \ref{fig:dfratio}(a)). At this frequency, the effect of defects is negligible, and $c_Q$ can be extracted. $R$ evolves from a maximum of 2.4, at the center of the bubble, to 1, in the bulk region surrounding it. By using Eq. (\ref{dfRF}b) (at fixed $f_{RF}$, $c_{SRH}$ and $R_{SRH}$), we convert $R(x,y)$ to $c_Q(x,y)$ (figure \ref{fig:quantumcap}(a)). The quantum capacitance can be seen to be minimal in the center of the bubble and increase monotonically towards the edge, from $1$ to $\mathrm{3.2\times 10^{-2}\ \frac{\mu F}{cm^2}}$.
The evolution of $c_Q(E,E^*)$ from the center of the bubble (blue) to the edge (red) obtained from $E_g$ (extracted from Ref. \cite{Lopez2016} , see Supporting Information 1) is shown in Supporting Information 9.
We then convert $c_Q(0\ V)$ to $E^*$ (color bar of figure \ref{fig:quantumcap}(a)). The increase of $E^*$ is due to an increase in n-type doping as we move from the center of the bubble towards the edge (Supporting Information 9).

At the very edge of the bubble, a lateral heterojunction occurs. The equalization of bulk and bubble Fermi levels ($E_{F,bulk}$ and $E_{F,bubble}$ respectively) to a pinned level, assumed to be $E_{F,bulk}$, gives rise to a non-zero contact potential difference between the two. We measured this via DC-EFM (figure \ref{fig:parabola}(b)) and obtained $V_{0,ML}-V_{0,Bulk} \simeq 200$ mV at the center of the bubble, with $V_{0,ML}>V_{0,Bulk}$, consistent with previous reports \cite{Tosun2015}. The variation of $E_c-E_F$ (top panel of figure \ref{fig:quantumcap}(b)) across the bubble reflects that of the CB. A representative bubble profile is shown in the background of figure \ref{fig:quantumcap}(b)-bottom. $E_C$ increases by $\sim$30 meV, when going from the bubble's edge towards the center. This is compatible with an intrinsic $E_{F,bubble}< E_{F,bulk}$. Pierucci et al. \cite{Pierucci2017} have demonstrated that the passivation of sulfur vacancies via hydrogen incorporation in MoS$_2$ MLs alters the band structure, by lowering $E_F$ toward VB. Not only the incorporation of hydrogen atoms at the bubble site would justify $E_{F,bubble}<E_{F,}$, but would also establish a correct correlation between the inhomogeneous distribution of defects, shown in figure \ref{fig:defectmap}(a), and the measured CB profile. The lower concentration of defects (possibly sulfur vacancies) at the top, compared to the edge, would indicate a locally higher H passivation activity at the topmost location. Consequently, the CB minimum at the top is expected to be further, with respect to the edge, in qualitative agreement with our findings. Note that the preferential H-absorption at the bubble center (where strain is highest) agrees with the work of Li et al. \cite{Li2016}, which predicts a decreasing H-adsorption free energy with increasing tensile strain.

\begin{figure}
    \centering
    \includegraphics[width=16cm]{Figure5.png}
    \caption{(a) Map of $c_Q(0\ V)$ obtained from the inset of figure \ref{fig:dfratio}(a) by using the ratio between Eq. (\ref{dfRF}a) and (\ref{dfRF}b). The evolution of $c_Q(0\ V)$ can be converted into the evolution of $E^*$. Note that in the red region (MoS$_2$ bulk) $c_Q$ and thus $E^*$ are not defined. Scale bar is 1 $\mathrm{\mu m}$. (b) Evolution of $E_C - E_F$ (top), 2D and 3D carrier density, along the bubble profile (bottom). The gray region in the background shows the bubble profile ($h_{max}$ =127 nm). The vertical dashed line indicates the transition from ML to bulk (in yellow). The horizontal dashed line is the MoS$_2$ bulk carrier density.}
    \label{fig:quantumcap}
\end{figure}

Once $E_F$ is determined, one can calculate the 2D carrier density, $n_{2D}=g_{2D}k_BT ln\left[1+e^{{(E_F-E_C)}/{k_BT}}\right]$ \cite{Bera2019, Ma2015}. This can then be converted into the 3D carrier density, $n_{3D}$, by dividing by the thickness of the ML (0.65 nm), so as to compare it to the density of the MoS$_2$ bulk. The bottom panel of figure \ref{fig:quantumcap}(b) shows the evolution of 2D and 3D carrier densities across the bubble, whose height profile is shown in gray in the background. They increase from 1.6 (2.5) to 5.1 $\mathrm{\times 10^9\ cm^{-2}}$ (7.9 $\mathrm{\times10^{16}\ cm^{-3}}$), when moving from the bubble center to its edge. Here, the $n_{3D}$ we find tends to the value found for MoS$_2$ bulk (horizontal dashed line) via four-probe resistance measurements, and is confirmed by the EFM spectroscopy (see Experimental, figure \ref{fig:simulation}). Consistent with the lower value of $E_F$, the carrier density measured on the bubble ($n_{3D}\sim 10^{16}\ cm^{-3}$) is lower than that in the bulk, in agreement with local passivation \cite{Pierucci2017}.

\section{Conclusions}

We used RF-EFM to investigate the electric field response of biaxially strained MoS$_2$ MLs in the form of bubbles. These bubbles were produced via H-ion irradiation of bulk crystals. The irradiation method induces strain in thick TMD crystals without the need to first isolate MLs through exfoliation or CVD growth. These latter are known to induce defect states in an uncontrollable manner, causing diffused dangling bonds and atomic scale defects (vacancies or grain boundaries). In contrast,  low energy treatments performed in H-plasma have been demonstrated to induce only a selective desulphurization of the topmost interface of the ML \cite{Lee2021} and possible partial passivation of sulfur vacancies due to H intercalation \cite{Pierucci2017}. Importantly, defect states affect the electrostatic response of MLs as much as their intrinsic quantum capacitance. Both contributions are DOS-related, and their investigation gives insight into the material's band structure. By using RF-EFM we separated the contribution of defect states from that of quantum capacitance: the former dominates the electrostatic response in the sub-MHz range. Our data are compatible with the existence of defect states at two different energies in the bandgap, with characteristic emission times of 0.1 and 1 $\mathrm{\mu s}$, respectively, and corresponding defect densities of the order of $\mathrm{10^{-10}\ eV^{-1}cm^{-2}}$. While the emission times are consistent with previous reports, the density values are one to two orders of magnitude smaller than that measured in CVD-grown or mechanically-exfoliated MLs, testifying the low invasiveness of our ML production procedure and/or H-induced passivation. At higher frequency, the quantum capacitance dominates instead. We imaged $c_Q$ across a bubble surface, at $f_{RF}=300$ MHz, demonstrating the correlation between quantum capacitance (electrostatic properties, band structure) and bubble topography (strain). From the local value of $c_Q$ we obtained $E_C-E_F$, and showed that Fermi level and CB minimum move progressively closer together, when going from the bubble center towards the edge. Finally, from $E_C-E_F$ we obtained the charge carrier density ($n_{3D}$ or $n_{2D}$). This increases monotonically from the center to the edge of the bubble, where it reaches the $n_{3D}$ value extracted from four probe resistance measurements on a non-irradiated twin MoS$_2$ bulk crystal.
Consistent with the lower value of $E_F$, the carrier density measured on the bubble ($n_{3D}\sim \mathrm{10^{16}\ cm^{-3}}$) is lower than that in the bulk, suggesting possible local S-vacancy passivation. Moreover, the evolution of charge carrier density correlates with the measured distribution of defect density (possibly sulfur vacancies), both being minimal at the top and maximal at the edge of the bubble. This suggests greater H passivation at the bubble center, where the ML is subjected to the maximum strain, in agreement with theoretical calculations \cite{Li2016}: the higher the tensile strain, the lower the H adsorption free energy of a MoS$_2$ strained surface. Note that, as shown in Ref. \cite{Bennet2023} , the variation of carrier density can severely limit the on-state capacitance of a MoS$_2$ ML-based device, to values well below its degenerate capacitance, and thus deserves attention for practical applications. Finally, we underline that a favored formation of single photon emitters has been recently found at the dome edges in similar systems \cite{Cianci2023} . As they can come from optically active defects, our result may correlate with the one discussed in \cite{Cianci2023} from a complementary point of view.

In conclusion, imaging the quantum capacitance at finite frequency, through RF-EFM, allowed us to gain insight into the effect of H-irradiation on the local MoS$_2$ ML doping level and defect states, and its variation on a sub-micrometric scale. This new technique can be used to investigate diverse electrostatic phenomena on the mesoscopic scale over and up to microwave frequencies. Additionally, being nondestructive as a non-contact AFM technique, RF-EFM is robust and can be used for a wide range of applications.

\section{Experimental}

\subsection{Sample Preparation and Proton Irradiation}

Commercial MoS$_2$ bulk crystals, from 2D Semiconductors, were first mechanically exfoliated with tape. Once a fresh and flat surface was obtained, they were placed in a vacuum chamber, equipped with a Kaufman ion source. The proton irradiation was performed at a base pressure of $<1 \times 10^{-6}$ mbar, and at a temperature of 120-150 $^\circ$C, with the crystal being grounded along the entire procedure. Hydrogen ions were obtained in an ionization chamber and accelerated by a system of biased grids, thus irradiating the sample with an ion beam having energy in the range of 10--20 eV. The samples were irradiated with a total dose in the range of 6-7 $\times 10^{16} ions/cm^{2}$.

\subsection{RF assisted-Electrostatic Force Microscopy}

The EFM measurements were performed by using an attoAFM I, from Attocube, equipped with an interferometric sensor for the detection of cantilever deflection and/or oscillation. The data acquisition system is controlled by Nanonis electronics from Specs, interfaced with ACC100 and ANC350 modules, from Attocube, for powering the interferometric cavity, and driving step motors and piezo movements.

We used conductive solid platinum AFM probes from Bruker (RMN-25PT400B), having resonance frequency in the range $f_0=5 \div 11$ kHz, and spring constant in the range $k=6 \div 14$  N/m. With the tip holder solidly connected to an excitation piezo, we preliminarily characterized the cantilever resonance curve by monitoring the cantilever oscillation amplitude, throughout the interferometric sensor, while sweeping the frequency of the AC bias supplied to the piezoelectric actuator. Upon identifying the cantilever resonance frequency, we worked in a closed loop mode: the cantilever oscillation amplitude was kept constant ($\sim$10 nm, much smaller than the tip-surface distance) by adjusting, if needed, the amplitude of the AC bias supplied to the piezo-actuator (figure 1d).

We measured the tip-sample electrostatic interaction, and its spatial and bias dependence, throughout the variation of the frequency of the oscillating cantilever, $df=f-f_0$, with $f_0$ and $f$ being the free resonance frequency and its interaction-driven shifted value. The frequency shift $df$ is indeed related to the interaction force, electrostatic in our case, $F_{el}$, as $df=-\frac{f_0}{2k}\frac{\partial F_{el}}{\partial z}$. The topographic measurements were acquired in closed-loop mode, i.e., by keeping, pixel by pixel, a constant $df$ (constant tip-surface separation). To do so, we first set up a tip-surface interaction satisfying the condition of negative $df$ (attractive interaction), at fixed supplied bias $V_{DC}$ ($F_{el}\propto V_{DC}^2$). We chose $df=-10\ Hz$ and $V_{DC}=2\ V$, which guarantees a tip-surface distance of $\sim$100 nm (well outside the vdW force range). Subsequently, a feedback loop was used to counterbalance for the variation of tip-surface distance, occurring while scanning due to the topographic features. The feedback circuit biased the scanning piezo so as to induce displacements along the z-direction counteracting for the morphological variations, and thus preserving the strength of the electrostatic interaction and the initial value of the frequency shift $df$. Conversely, in spectroscopic mode, the feedback loop was disabled and the variations of $df$ (vs $V_{DC}$ or $f_{RF}$) were recorded at fixed tip-surface distance. We used the measured $df(V_{DC})$ spectra of MoS$_2$ bulk and Au, at tip-sample separations of 192 and 203 nm, respectively, and their fit to quantify $\frac{\partial^2 c_T}{\partial z^2}$. This latter was used to calculate the electrostatic force, $F_{el}$, acting between the AFM probe and both samples, at fixed separation, resulting in the solid lines of figure \ref{fig:simulation}. By describing the tip shape as an axially symmetric paraboloid, Xu et al. \cite{Xu2019} have numerically simulated the electrostatic force acting between a metallic AFM probe and a sample with known carrier density, $n$. Following their formalism, we have simulated $F_{el}$ by assuming a charge carrier density of $n_{Au}=5.9\times 10^{22}\ cm^{-3}$ \cite{AuDatabase} and $n_{MoS_2}=1.6\times 10^{17}\ cm^{-3}$ (the latter extracted from four-probe resistance measurement on a twin MoS$_2$ crystal). The numerical simulations (scatters in figure \ref{fig:simulation}) reproduce the experimental measurements (solid lines), thus demonstrating that EFM experiments can also lead to the quantitative evaluation of the charge carrier density.

\begin{figure}
    \centering
    \includegraphics[width=10cm]{Figure6.png}
    \caption{$F_{el}(V)$ spectroscopy (solid lines) and numerical simulations (scatters) of tip-Au and tip-MoS$_2$ electrostatic interaction, at 203 and 192 nm, respectively, obtained by using the values of charge carrier densities $n$ reported in the legend.}
    \label{fig:simulation}
\end{figure}

Finally, the maps of $df$(-7 V) - figure \ref{fig:defectmap}(a) - and $R$(300 MHz) - figure \ref{fig:dfratio}(c)-left -- were obtained by acquiring the value of $df$ (at $-7$ V or $300$ MHz), per each pixel, at the same tip-sample separation. To do so, the initial feedback conditions ($df=-10$ Hz and $V_{DC}=3$ V) were always restored before moving from one pixel to the next, so as to ensure that the same tip-surface distance ($\sim$153 nm) is set up per each point. Once moved, the feedback was disabled, the bias was swept to the desired conditions, and the local value of $df$ was measured (at fixed height).

 We used arbitrary wave function generators (AFG3252 from Tektronix or Hameg HM8134 from Rohde \& Schwarz) to supply the tip with a sinusoidal voltage signal, by keeping the sample grounded. The output signal was first amplified and then summed up to the DC bias supplied by Nanonis, throughout the use of a bias tee (ZFBT-6GW from Mini-Circuits - inset of figure 1b). Finally, it travelled to the tip within a microwave coaxial cable. As discussed in Supporting Information 7, we found a RF-frequency dependent attenuation of the RF power, which decreases as $f_{RF}$ increases (as also shown in the top panel of figure 4(a)). We also found that, differently from the case of MoS$_2$ bubble, the $df(f_{RF})$ spectra measured on Au and MoS$_2$ bulk crystals are almost undistinguishable, and thus only affected by the attenuation of the RF line.

Finally, as discussed in Supporting Information 2, we have adopted a sphere-plane (sphere-sphere) model to describe the electrostatic interaction between the tip and Au or MoS$_2$ bulk crystal (MoS$_2$ bubble). By doing so, we have analytically modelled the quantities involved in the problem, such as $F_{el}$ and $c_T$ (and its derivatives).

\begin{acknowledgement}

C.D.G. acknowledges financial support from MUR (Italian Ministery for University and Research) under the project PON-AIM (Programma Operativo Nazionale - Attraction and Internazional mobility).

\end{acknowledgement}

\setcounter{equation}{0}
\setcounter{figure}{0}

\begin{suppinfo}
\renewcommand{\theequation}{S.\arabic{equation}}
\renewcommand{\thefigure}{S.\arabic{figure}}

\section{1. Bubble characterization}

Figure \ref{fig:asimmetry}(a) of the main text shows a representative ($\mathrm{13\ \mu m \times 9.2\ \mu m}$) AFM topography of an irradiated MoS$_2$ bulk crystal. Numerous localized bubbles are observed on the crystal surface, having maximum height $h_{max}$ and radius $R$ evolving in a dependent manner (bigger protrusions being taller, and smaller being shorter), so as to preserve the universal law regulating the $\frac{h_{max}}{R}$ ratio \cite{Khestanova2016}. This ratio dictates the magnitude of the radial, $\varepsilon _r,$ and tangential, $\varepsilon_t$, strain components affecting the membrane, as \cite{Blundo2021PRL}:

\begin{subequations} \label{strain}
    \begin{empheq}[left=]{align}
      & \varepsilon_t(r)=f(\nu,q){\left(\frac{h_{max}}{R}\right)}^2\left[1-{\left(\frac{r}{R}\right)}^{2q-2}\right] \\
      & \varepsilon_r(r)=f(\nu ,q){\left(\frac{h_{max}}{R}\right)}^2\left[1-g(\nu ,q){\left(\frac{r}{R}\right)}^{2q-2}\right]
    \end{empheq}
\end{subequations}

Both components increase from the bottom toward the top, where the total strain becomes equally biaxial. Here, $0<r<R$ is the lateral coordinate of the bubble profile, while $f(\nu,q)$ and $g(\nu,q)$ are tabulated functions of the Poisson ratio $\nu$, and of the exponent $q$. This latter governs the analytical shape of the bubble as \cite{Blundo2021PRL}:

\begin{equation} \label{profile}
h(r)=h_{max}\left[1-{\left(\frac{r}{R}\right)}^q\right]
\end{equation}

The left-hand plot of figure \ref{fig:prof} shows a typical profile taken along one axial direction of the bubble surface (a white dashed line is drawn in the inset as a guide for the eye). The experimental profile (black scatters) is fitted, in red, by Eq. \ref{profile}, resulting in $h_{max}$=126 nm, $\mathrm{R=1.4\ \mu m}$, and $q=2$. These latter parameters have been used to evaluate the strain components $\varepsilon_t(r)$ and $\varepsilon_r(r)$, from Eqs. (\ref{strain}a) and (\ref{strain}ab), respectively. The total strain $\varepsilon(r)=\varepsilon_t(r)+\varepsilon_r(r)$ is plotted in the top panel of figure \ref{fig:prof}-right. Here, $\varepsilon(r)$ is the average strain calculated over the entire bubble surface, which may slightly differ ($<0.3\%$) from the local one (measured over a single profile), given a not always perfect bubble axial symmetry. The bottom panel, instead, shows the evolution of the strain-dependent bandgap across the bubble surface, being minimum at the maximum strain location (top of the bubble) and increasing as the strain decreases (extracted from \cite{Lopez2016}).

\begin{figure}
    \centering
    \includegraphics[width=16cm]{S1.png}
    \caption{Left: experimental profile taken along one axial direction of the bubble surface (black scatters) and its fit (red solid line) by Eq. (\ref{profile}). Right: evolution of the total strain (top) and $E_g$ (bottom), along the bubble profile. Inset: 3D bubble topography. Dashed white line is indicative of an axial profile direction. Yellow arrows indicate tangential and circumferential strain directions.}
    \label{fig:prof}
\end{figure}

\section{2. Analytical description of tip-sample electrostatic interaction}

By modelling tip and sample surface as the plates of an arbitrary shape capacitor, the electrostatic force acting between the two can be written as:

\begin{equation} \label{Fel}
F_{el}=\frac{1}{2}\frac{\partial}{\partial z}\left(c_TV_T^2\right)
\end{equation}

Here, $c_T$ is the tip-surface capacitance, strongly dependent on their separation $z$ and shape, and $V_T$ is the drop of potential between the two.

If the sample is either a metal, or a highly conducting semiconductor (very small penetration depth, or Debye length $L_D$), all the voltage $V$ supplied drops between the tip and the sample surface, thus $V_T=V$ is $z$-independent. Conversely, when the sample conducts poorly, the $V$ supplied is redistributed as $V=V_T(z)+V_S(z)$, with $V_S$ being the drop of potential inside the investigated sample. This latter can be described as an effective capacitance $c_S$, in series with $c_T,$ so that $V_T(z)=\frac{c_s}{c_S+c_T(z)}V.$ Given $V=V_{DC}-V0$, with $V_0$ being the contact potential difference, this translates into:

\begin{subequations} \label{Felasimmetry}
    \begin{empheq}[left={F_{el}=\frac{1}{2}(V_{DC}-V_0)^2\frac{\partial c_{T,eq}}{\partial z}\empheqlbrace\,}]{align}
      & c_{T,eq}=c_T, \ \ \  \ \ \ (L_D\simeq 0) \\
      & c_{T,eq}=\frac{c_Tc_S}{c_S+c_T} \ \ \ \ \ L_D\ge L^*_D
    \end{empheq}
\end{subequations}

Here $L^*_D$ is the distance between the plates of an equivalent plane capacitor, ${c_T}^*=\frac{\varepsilon S}{L^*_D}$, such that $c_T^*$ is comparable to $c_T$. If $c_S$ is comparable to $c_T^*$, or $c_S\le c_T^*$, its contribution into the series cannot be neglected. Notice that, as discussed in detail in Ref. \cite{Sadeghi2012} , no minus sign should appear explicitly in Eq. (\ref{Felasimmetry}). However, since  $\frac{\partial c_T}{\partial z}<0$ the force stays negative, reflecting its attractive nature.

The cantilever's frequency shift $df$ can be calculated from Eqs. (\ref{Felasimmetry}4), by using $df=-\frac{f_0}{2k}\frac{\partial F_{el}}{\partial z}$, with $f_0$ being the cantilever's free resonance frequency and $k$ its spring constant:

\begin{subequations} \label{asymDF}
    \begin{empheq}[left={df=-\frac{f_0}{4k}(V_{DC}-V_0)^2\frac{\partial^2 c_{T,eq}}{\partial z^2}\empheqlbrace\,}]{align}
      & c_{T,eq}=c_T, \ \ \  \ \ \ L_D\simeq 0 \\
      & c_{T,eq}=\frac{c_Tc_S}{c_S+c_T} \ \ \ \ \ L_D\ge L^*_D
    \end{empheq}
\end{subequations}

The analytical expression of $c_T$ depends on the shape of the capacitor's plates. By following a generalized Derjaguin approximation \cite{Schnitzer2015}, the electrostatic force acting between two metallic objects of arbitrary shape can be written as:

\begin{equation} \label{Deraj}
F_{el}\simeq 2\pi R^*\int{\Phi_{el}}dz
\end{equation}

Here, $R^*$ is the so-called effective radius - in the case of two spheres, having radius of curvature $R_1$ and $R_2$, respectively, $\frac{1}{R^*}=\frac{1}{R_1}+\frac{1}{R_2}$, and $\Phi_{el}=\frac{\varepsilon }{2z^2}(V_{DC}-V_0)^2$ is the electrostatic force per unit area acting between the parallel plates of a plane capacitor, with $\varepsilon $ being the absolute dielectric constant.

We used a plane-sphere interaction model for the analytical description of the electrostatic force between the tip and both Au and MoS$_2$ bulk surfaces, and a sphere-sphere model for describing the tip-bubble interaction. By doing so, in the case of total electric field screening:

\begin{equation} \label{effectFel}
F_{el}\simeq -\frac{\pi \varepsilon R^*}{z}(V_{DC}-V_0)^2
\end{equation}

with $R^*=R_{tip}$ in the first case (plane-sphere interaction) and $\frac{1}{R^*}=\frac{1}{R_{tip}}+\frac{1}{R_{bubble}}$ in the second one (sphere-sphere interaction). By comparing Eqs. (\ref{Felasimmetry}a) and (\ref{effectFel}), $\frac{\partial c_T}{\partial z}=-\frac{2\pi \varepsilon R^*}{z}$. Under this assumption, at sufficiently high tip-surface separation (negligible van der Waals forces), the cantilever's frequency shift $df$, only due to the electrostatic interaction, is:

\begin{equation} \label{effectdf}
df=-\frac{f_0}{2k}\frac{\pi \varepsilon R^*}{z^2}(V_{DC}-V_0)^2
\end{equation}

Figure \ref{fig:decay} shows the plots of $\frac{df}{V^2}$ vs $\Delta z$ for Au (left), MoS$2$ bulk (middle) and MoS$2$ bubble, at positive bias, (right). Here $\Delta z=z-z_0$ is the experimentally measured tip-surface separation rescaled with respect to their reciprocal initial distance, at the approach. Note that, to guarantee the lowest invasiveness possible, the approach was always performed in non-contact mode, with a non-zero bias applied, ensuring for high ($\ge$ 100 nm) tip-sample approach distance. The fits, in red, have been used to extract the effective radius $R^*$. By using a conductive solid platinum AFM probe from Bruker (RMN-25PT400B), having resonance frequency of $f_0=7.8$ kHz, and spring constant of $k=8$ N/m, we found $\mathrm{R^*\simeq 3\ \mu m}$ in the case of Au and MoS$2$ bulk, and $\mathrm{R^*\simeq 2.2\ \mu m}$ in the case of MoS$2$ bubble (positive voltage). While the former value of $R^*$is a direct measure of $R_{tip}$ (radius of the equivalent sphere used to model the tip/Au and tip/MoS$2$ bulk interaction as a sphere-plane interaction), the latter value is compatible with the convolution of $\mathrm{R_{tip}\simeq 3\ \mu m}$ and the bubble's radius, by using a sphere-sphere modelling.

\begin{figure}
    \centering
    \includegraphics[width=12cm]{S2.png}
    \caption{Plot of $\frac{df}{(V_{DC}-V_0)^2}$ vs $\Delta z$ for Au (left), MoS$_2$ bulk (middle) and MoS$_2$ bubble, at positive bias, (right).}
    \label{fig:decay}
\end{figure}

We underline that, given the geometry of the probe used, the electrostatic contribution from the cantilever is negligible. Indeed, by modelling this latter as the second plate of a plane capacitor:

\begin{equation} \label{dfcant}
df=-\frac{f_0\varepsilon}{2k}\left[\frac{\pi R^*}{z^2}+\frac{S}{{(z+h)}^{3}}\right](V_{DC}-V_0)^2
\end{equation}

 Here, $S$ and $h$ are the cantilever surface area and the height of tip's cone. Given $\mathrm{h\sim 100 \mu m}$ and $\mathrm{S\sim 24\times 10^3 \mu m^2}$ \cite{RMN}, the electrostatic contribution of the cantilever is several orders of magnitude smaller than the one of the tip, as also discussed in \cite{ChuMM2020}.

 Finally, more complex models of the tip shape have been considered to describe the electrostatic tip-sample interaction. However, while adding complexity to the calculations, they would only contribute with negligible corrections to the sphere-plane model.

 \section{3. Derivation of bias-sign dependent tip-MoS$_2$ bubble electrostatic interaction}

 As discussed in the main text, the electrostatic behavior of the MoS$_2$ bubbles is bias-sign dependent: for $V>0$, the electric field screening is complete ($c_B=\infty $); for $V<0$, the electric field penetration is non-zero ($c_B$ has a finite value). Here $c_B$ is the total MoS$_2$ bubble capacitance, which is a combination of defect states capacitance, $c^*_{SRH}$, and quantum capacitance, $c_Q$. Thus, in presence of the effective capacitance $c_B$, in series with $c_T$, $V_T=\frac{c_B}{c_B+c_T}V$. In the limit of $c_B\gg c_T$, $V_T\simeq \left(1-\frac{c_T}{c_B}\right)V$. It follows:

\begin{subequations} \label{Vasimmetry}
    \begin{empheq}[left={F_{el}=\frac{1}{2}(V_{DC}-V_0)^2\frac{\partial c_{T,eq}}{\partial z}\empheqlbrace\,}]{align}
      & c_{T,eq}=c_T, \ \ \  \ \ \ V>0 \\
      & c_{T,eq}=c_T\left(1-\frac{c_T}{c_B}\right)^2 \ \ \ \ \ V<0
    \end{empheq}
\end{subequations}

Equations (\ref{dfasimmetry}) of main text are obtained by deriving Eqs. (\ref{Vasimmetry}).

Thus, depending on the bias polarity:

\begin{subequations} \label{explicitFel}
    \begin{empheq}[left=]{align}
      & V>0 \ \ \ \ \ F_{el,pos}=\frac{1}{2}(V_{DC}-V_0)^2\frac{\partial c_T}{\partial z} \\
      & V<0 \ \ \ \ \ F_{el,neg}\simeq \frac{1}{2}(V_{DC}-V_0)^2\frac{\partial c_T}{\partial z}\left(1-\frac{4c_T}{c_B}\right)
    \end{empheq}
\end{subequations}

In the limit of small cantilever oscillations $\delta z$, around the equilibrium position $z_0$, the tip-sample separation $z$ and the capacitance $c_T\ $can be written as:

\begin{subequations} \label{oscillations}
    \begin{empheq}[left=]{align}
      & z=z_0+\delta z(f,t)=z_0+A(f)cos(2\pi ft) \\
      & c_T(z)=c_T(z_0+\delta z)\simeq c_T(z_0)+\left.\frac{\partial c_T}{\partial z}\right|_{z=z_0}\delta z(f)
    \end{empheq}
\end{subequations}

Here $f$ and $A(f)$ are the excitation frequency and the Lorentzian-like cantilever oscillation amplitude, respectively. By substituting Eq. (\ref{oscillations}b) in (\ref{explicitFel}b):

\begin{footnotesize}
\begin{equation} \label{finalFel}
F_{el, neg}\simeq \frac{1}{2}(V_{DC}-V_0)^2\frac{\partial c_T}{\partial z}\left[1-\frac{c_T\left(z_0\right)}{c_B}-\frac{4}{c_B}\left.\frac{\partial c_T}{\partial z}\right|_{z=z_0} \delta z\left(f,t\right)\right]=F_{el,pos}+F_{stat}\left(c_B\right)+F_{dyn}\left(c_B,\ f\right)
\end{equation}
\end{footnotesize}

 which reproduces Eq. \ref{elecForce} of main text.

 \section{4. Cantilever resonance curves modelling}

 We modelled: i) the cantilever's resonance curves acquired on MoS$_2$ bulk (at any $V$) and strained ML, at positive $V$, by using the solution of a classic harmonic oscillator in presence of a single driving force (external excitation due to the piezoelectric actuator, $F_{ex}=F_0cos(2\pi ft)$); and ii) the resonance curve acquired on the strained ML at negative $V$ as an harmonic oscillator in presence of two dynamic forces, $F_{ex}$ and $F_{dyn}(c_B,f)=F_{dyn,0}(c_B,f)cos(2\pi ft)$. At $z=z_0$, $F_{dyn,0}(c_B,f)=-\frac{2}{c_B}(V_{DC}-V_0)^2\left(\left.\frac{\partial c_T}{\partial z}\right|_{z=z_0}\right)^2A(f)$. Expected oscillation amplitudes are:

 \begin{subequations} \label{amplitude}
    \begin{empheq}[left=]{align}
      & Bulk,\ any\ V;\ Strained\ ML,\ V>0 \rightarrow A(f)=\frac{\alpha}{\sqrt{\left[(2\pi f_0)^2-(2\pi f)^2\right]^2+\frac{(2\pi f)^2}{\beta^2}}}\\
      & Strained\ ML,\ V<0 \rightarrow A^*(c_B, f)=\frac{\alpha^*(c_B, f)}{\sqrt{\left[(2\pi f_0)^2-(2\pi f)^2\right]^2+\frac{(2\pi f)^2}{\beta^2}}}
    \end{empheq}
\end{subequations}

Here $\alpha =\frac{F_0}{m}$ and $\alpha^*(c_B,f)=\frac{F_0+F_{dyn,0}(c_B,f)}{m}$, with $\alpha^*$ being frequency-dependent and smaller than $\alpha$ ($F_{dyn,0}(c_B,f)<0$). Finally, $m$ is the mass of the cantilever, $m=\frac{k}{f^2_0}$, and $\beta$ is the oscillator damping coefficient. By fitting the cantilever's resonance curves, acquired on MoS$_2$ bulk with equation (\ref{amplitude}a), we extracted a damping coefficient of $\beta =0.012$ s (quality factor of 588), and a driving force of $5.4$ nN in amplitude. By using the so found $\beta$, $F_0$ and $A(f)$, and by using the analytical expression  $\left.\frac{\partial c_T}{\partial z}\right|_{z=z_0}=-\frac{2\pi \varepsilon R^*}{z_0}$, we fitted the resonance curve measured on the strained ML at $-6V$, by using equation (\ref{amplitude}b) and keeping $c_B$ as the only free parameter.

 \section{5. Capacitance comparison}\label{cap}

\begin{figure}
    \centering
    \includegraphics[width=10cm]{S3.png}
    \caption{Capacitances involved in the EFM experiments discussed in the main text.}
    \label{fig:capacitances}
\end{figure}

The electrostatic tip-surface capacitance can be approximated as $c_T(z)=2\pi \varepsilon R^*ln\left(\frac{z+R^*}{z}\right)$ \cite{Li2022}, which, in the limit of $R^*\gg z$, reproduces $\frac{\partial c_T}{\partial z}=-\frac{2\pi \varepsilon R^*}{z}$, as obtained by comparing Eqs. (\ref{Felasimmetry}a) and (\ref{effectFel}). Figure \ref{fig:capacitances} shows the z-dependence of $c_T$, both in the case of sphere-plane and sphere-sphere interaction (given the $R^*$ values characterized in the previous section), and compares it to the other capacitances involved in the problem, such as the depletion capacitance of MoS$_2$ bulk, $c_{dep}$, the defect states capacitance, $c^*_{SRH}$, and the quantum capacitance $c_Q$ (both measured, and degenerate, $c_{Q,deg}$). While $c^*_{SRH}$ and $c_Q$ are experimentally measured, $c_{Q,deg}=e^2g_{2D}=84\frac{\mu F}{cm^2}$ \cite{Fang2018}, and $c_{dep}$ is defined as $c_{dep}=\frac{\varepsilon S}{L_D}$. Here, $S$ is the involved surface area, chosen as $\pi R^2_{tip}$, and $L_D\simeq$ 7.9 nm is the Debye length derived from the free carrier density estimated by four probe measurements on a twin MoS$_2$ crystal ($n_{MoS2}=\mathrm{1.6\times 10^{17}cm^{-3}}$).

 \section{6. Mapping the contact potential difference: Kelvin probe force microscopy}

The effective drop of voltage $V$ between the tip and the sample surface is shifted when compared to the bias voltage supplied, $V_{DC}$, by an amount equal to the \textit{contact potential difference} $V_0$, as $V=V_{DC}-V_0$. This latter leads the information on the mismatch between sample and tip work function, $W_S$ and $W_T$ respectively, given $-eV_0=W_T-W_S$ \cite{Melitz2011}. To map the local variation of $V_0$ across the bubble surface, we performed Kelvin probe force microscopy (KPFM), by using a Bruker MultiMode AFM, under ambient conditions. A Pt-Ir coated tip (SCM-PIT, from Bruker) was used to scan the surface in tapping mode, first. Subsequently, the local value of $V_0$ was acquired, in lift mode, at a constant tip height of 50 nm, with the tip being supplied an AC voltage of 1.5 V at its resonance frequency ($\sim$75 kHz).

\begin{figure}
    \centering
    \includegraphics[width=16cm]{S4.png}
    \caption{topography (a) and contact potential (b) maps acquired on the surface of a hydrogenated MoS$_2$ crystal.}
    \label{fig:KPFM}
\end{figure}

Figures \ref{fig:KPFM}(a)-(b) compare topography and contact potential maps acquired in a region of the hydrogenated MoS$_2$ bulk crystal surface. A clear correspondence arises between the two: in qualitative accordance with EFM, the bubbles support higher $V_0$ (white contrast), compared to the bulk (blue). Moreover, $V_0$ smoothly decreases from the bubble's top toward the edge, in agreement with the trend measured by EFM, and shown in figure \ref{fig:parabola}(b) of the main text, by increasing the averaging surface area (increasing tip-surface separation). Importantly, all bubbles manifest a similar behavior, with higher contact potential difference, compared to the bulk. Finally, while the absolute values of bulk and bubble's $V_0$ are tip-dependent, their tip-independent difference is slightly underestimated compared to the EFM measurements ($\sim$60 mV vs $\sim$205 mV). The reduced contrast can be attributed to a partial lack of metallic coating at the very tip apex, which can be often occur when scanning the topography in tapping-mode.

 \section{7. Calibration of RF power}

 Figure \ref{fig:DCRF} shows the $df$ spectroscopy as a function of $V_{DC}$ and $V_{RF}$, at fixed frequency $f_{RF}$ =20 MHZ. Under the assumption $\left\langle {\left[V_{RF}cos (2\pi f_{RF}t)\right]}^2\right\rangle =\frac{V_{RF}^2}{2}$, which is a good approximation for a RF frequency much higher than cantilever's resonance frequency:

\begin{equation} \label{RFdf}
df\simeq -\frac{f_0}{2k}\left[(V_{DC}-V_0)^2+\frac{V^2_{RF}}{2}\right]\frac{\partial^2c_T}{\partial z^2}
\end{equation}

At $V_{DC}=V_0$, $df$ should exhibit a quadratic dependence on $V_{RF}$. Figure \ref{fig:DCRF}(a) shows the $df(V_{RF})$ spectroscopy, measured on Au, and fitted in red by a parabolic function, confirming the expected behavior. Figures S\ref{fig:DCRF}(b) and (c) show the $df(V_{DC})$ spectroscopy, measured on Au and MoS$_2$ bulk, respectively, at $V_{RF}=0$ (DC), and at fixed $V_{RF}=3V$, $f_{RF}$=20 MHZ. The vertical shift between the two is due to the additional $\frac{V^2_{RF}}{2}$ term. From the quadratic fits, in red, we extracted the effective $V_{RF}$ bias at the tip-sample junction, resulting in $2.8\div 2.9$ V. Thus, the value obtained from the fit confirms that microwave attenuation from the cables is negligible at the working frequency of 20 MHZ.

\begin{figure}
    \centering
    \includegraphics[width=16cm]{S5.png}
    \caption{(a) $df(V_{RF})$ spectroscopy, acquired on Au at $V_{DC}=V_0$ and $f_{RF}$=20 MHZ. $df(V_{DC})$ spectroscopy, acquired on Au (b) and MoS$_2$ bulk (c), at $V_{RF}=0$ and $V_{RF}$=3 V, at $f_{RF}$=20 MHZ. Quadratic fits are reported in red.}
    \label{fig:DCRF}
\end{figure}

Finally, Figure \ref{fig:RF} shows the $df$ dependence on the RF frequency, as measured on Au, in orange, on MoS$_2$ bulk, in black, and on top of a Mos$_2$ bubble, in blue. The overall increasing trend of $df_{RF}$ on $f_{RF}$ is due to the frequency-dependent attenuation of the $RF$ power along the transfer line. The effective RF bias to be considered in Eqs.\ref{dfRF} is thus $V^*_{RF}=V_{RF}\chi (f_{RF})$, with $\chi (f_{RF})$ being the transfer function of the RF line. Since this latter is generally unknown and hard to calibrate, one can get rid of it by considering the ratio between Eqs. \ref{dfRF} (figure \ref{fig:dfratio}(a)). Importantly, as long as a good electric field screening is provided, as in the case of Au and MoS$_2$ bulk, the $df_{RF}(f_{RF})$ spectroscopy is sample-independent. Differently, when the electric field penetrates inside the material, as in the case of MoS$_2$ bubble, the effect of the sample-to-ground capacitance plays a non negligible role in the $df_{RF}(f_{RF})$ spectrum.

\begin{figure}
    \centering
    \includegraphics[width=10cm]{S6.png}
    \caption{$df(f_{RF})$ spectroscopy, at fixed $V_{RF}$, as measured on Au, in orange, on MoS$_2$ bulk, in black, and on MoS$_2$ in blue.}
    \label{fig:RF}
\end{figure}

\section{8. Fit of $df$ ratio vs $f_{RF}$ spectroscopy - single vs double defect states}

Figure \ref{fig:fit} compares the $df$ ratio vs $f_{RF}$ spectroscopy fitted by assuming that the DOS of the MoS$_2$ ML disclose a single (a) or double (b) defect state. As shown by the circuital representation in the inset, this latter configuration turns into adding a parallel of capacitance and resistance branch pair \cite{Lehovec1966}. The fitting function is obtained by dividing Eq. (\ref{dfRF}a) by Eq. (\ref{dfRF}b) of the main text:

\begin{equation} \label{dfratio}
\frac{df_{RF,bulk}}{df_{RF, bubble}}=Abs\left[\frac{\frac{\partial^2c_T}{\partial z^2}}{\frac{\partial^2}{\partial z^2}\left[c_T{\left(\frac{1+i2\pi f_{RF}(c_Q+c_{SRH})R_{SRH}}{1+i2\pi f_{RF}(c_T+c_Q+c_{SRH})R_{SRH}}\right)}^2\right]}\right]
\end{equation}

where, in the case of two defect states, $c_{SRH}(f_{RF})$ and $R_{SRH}(f_{RF})$ are the parallel of $c_{A,B}(f_{RF})=c^*_{A,B}\frac{arctg(2\pi f_{RF}\tau _{A,B})}{2\pi f_{RF}\tau_{A,B}}$, and $R_{A,B}(f_{RF})=\frac{1}{c^*_{A,B}}\frac{\tau_{A,B}}{ln\left(1+(2\pi f_{RF})^2\tau _{A,B}^2\right)}$ ($c^*_{A,B}=\frac{e^2D_{A,B}}{2}$). As shown by figure \ref{fig:fit}, given the same quantum capacitance $c_Q$, which sets the uppermost level of the $df$ ratio at high frequency, the fit including two defect reproduces better the data.

\begin{figure}
    \centering
    \includegraphics[width=16cm]{S7.png}
    \caption{$df$ ratio vs $f_{RF}$ spectroscopy, at fixed $V_{RF}$, fitted by single (a) and double (b) defect state. Different inner scatter colors (yellow and black) refer to to two sets of data taken with smaller steps on a smaller frequency range (yellow), and with larger step in a larger frequency range (black).}
    \label{fig:fit}
\end{figure}

\section{9. Quantum Capacitance}

Figure \ref{fig:plotcq} shows the evolution of $c_Q(E, E^*)$ when moving from the top of the bubble, in blue, toward its edge, in red, obtained once including the position-dependent values of the bandgap $E_g$ (extracted from \cite{Lopez2016} as a function of position, see Supporting Information 1). The location-dependent $E^*$ is point by point univocally defined by the energy difference $E_c-E_F$ such that the simulated
value of $c_Q(E_F,E^*)$ matches, at zero energy, the locally measured $c_Q(0V)$. As shown by the inset, such a match is achieved by shifting $E_c$ progressively closer to $E_F$ (enhancing the n-type doping of the strained ML), while gradually moving from the top of the bubble toward the edge.

\begin{figure}
    \centering
    \includegraphics[width=10cm]{S8.png}
    \caption{Main: evolution of $c_Q(E,E^*)$ while moving from the bubble top (blue) to its edge (red). Inset: zoom of $c_Q$ around $E_F$.}
    \label{fig:plotcq}
\end{figure}

\end{suppinfo}


\section{Introduction}\label{sec:intro}

In studies with a large number of observations, statistical testing procedures are prone to detect even minor departures from a null hypothesis yielding very small p-values. However, since significance does not automatically imply relevance, measures for the size of the effect associated with a possible rejection of the null hypothesis are often recommended as a useful tool, see e.g. \citet{wilkinson1999statistical}.

A well known measure in a regression context when a quantitative variable $Y$ depends on independent regressors (quantitative and/or qualitative) is Cohen's $f^2$, see \citet[Chapt. 9]{cohen1988statistical}. Under multivariate normality this measure is strongly related to an $F$ test of a linear hypothesis that a subset $B$ of independent variables does not substantially contribute to the explanation of $Y$, given a set $A$ of independent regressors already included in the model. One may argue that even if the regression coefficients associated with variables $B$ significantly differ from zero, their contribution may be assessed as relevant when some meaningful size of their effect is measured.

In this light studies may additionally want to report $f^2$ values for each variable in a regression model given the others, see for example Tables 2 to 5 in \citet{taylor2020reactions}. With regard to the procedure of calculating $f^{2}$, the authors refer to \citet{selya2012practical} who introduced a practical method in relation with SAS\textsuperscript{\textregistered} software. This approach even carries over to linear mixed models, i.e. the case that some independent variables are associated with random effects rather than with fixed effects. Hence, the original approach by \citeauthor{cohen1988statistical} is generalized to a certain extent, involving the additional estimation of a variance-covariance matrix.

In the following we reconsider and discuss the generalization of $f^2$ to linear mixed models. We also point to an alternative computational procedure which turns out to be especially useful in context with the statistical software environment {\sf R} \citep{Rsoftware} and well known package {\tt lme4}, see \citet{R:lme4}, for fitting linear mixed models. Our explanations are illustrated on the basis of an artificially generated data set.

\section{Linear Mixed Model}\label{sec:lmm}

Consider a linear mixed model (LMM) described by
\begin{equation}\label{E1}
\bm{y} = \bm{X} \bm{\beta} + \bm{Z} \bm{u} + \bm{e}\; ,
\end{equation}
where $\bm{y}$ is an $n\times 1$ observable random vector. It is assumed that the $n\times p$ model matrix $\bm{X}$
of full column rank $p$ can be partitioned as
\begin{equation}\label{E2}
\bm{X} = (\bm{1}_{n} : \bm X_{1}: \bm{X}_{2})\; ,
\end{equation}
where $\bm{1}_{n}$ denotes the $n\times 1$ vector of ones, while the $n\times p_{1}$ and $n\times p_{2}$ matrices $\bm{X}_{1}$ and $\bm{X}_{2}$ contain the values of $p_{1}+ p_{2} = p-1$ regressors. The $p\times 1$ vector $\bm{\beta}$ comprises $p$ unknown parameters $\beta_{0}, \beta_{1}, \ldots, \beta_{p-1}$ addressed as fixed effects. The $n\times q$ matrix $\bm{Z}$ contains the values of independent variables associated with an $q\times 1$ vector $\bm{u}$ of unobservable random effects. It is assumed that $\bm{u}$ has expectation $\bm{0}_{q}$ (the $q\times 1$ vector of zeroes) and variance covariance matrix $\mathop{\operator@font Cov}\nolimits(\bm{u}) = \sigma^2 \bm{D}$ with unknown parameter $\sigma^{2} >0$ and $q\times q$ matrix $\bm{D}$, which may depend on further unknown parameters.
For the $n\times 1$ vector $\bm{e}$ of unobservable random errors it is assumed that $\mathalpha{\operator@font E}(\bm{e})=\bm{0}_n$ and
$\mathop{\operator@font Cov}\nolimits(\bm{e}) = \sigma^2 \bm{T}$, where the $n\times n$ positive definite matrix $\bm{T}$ may also depend on unknown parameters. Moreover, the assumption $\mathop{\operator@font Cov}\nolimits(\bm{u}, \bm{e}) = \bm{0}_{q,n}$ (the $q\times n$ matrix of zeroes) implies
\begin{equation}\label{E3}
\mathop{\operator@font Cov}\nolimits(\bm{y}) =  \sigma^2 \bm V, \quad \bm{V} = \bm{Z} \bm{D} \bm{Z}^{T} + \bm{T}\; .
\end{equation}
Then the above model may also be represented by the triplet $\{\bm{y}, \bm{X}\bm{\beta}, \sigma^2 \bm{V}\}$ and can be considered as a special case of the general Gauss-Markov model, see e.g. \citet{gross2004general}. As a matter of fact, formulas useful under model \ref{E1}) carry over from  a classical regression model $\bm{Q}^{-1} \bm{y} = \bm{Q}^{-1} \bm{X} \bm{\beta} + \bm{\varepsilon}$ with $\mathalpha{\operator@font E}(\bm{\varepsilon}) = \bm{0}_{n}$ and $\mathop{\operator@font Cov}\nolimits(\bm{\varepsilon}) = \sigma^2 \bm{I}_{n}$, see \citet[Sect 2.7]{christensen2020plane}. Here $\bm{Q}$ denotes some nonsingular matrix satisfying $\bm{Q}\bm{Q}^{T} = \bm{V}$. Let us assume for the moment that  $\bm{V}$ is completely known. Then
\begin{equation}\label{E4}
\widehat{\bm{\beta}} = (\bm{X}^{T} \bm{V}^{-1}\bm{X})^{-1} \bm{X}^{T} \bm{V}^{-1} \bm{y}
\end{equation}
is the best linear unbiased estimator for $\bm{\beta}$, and
\begin{equation}\label{E5}
\widehat{\sigma}^{2} =  \frac{1}{\nu} (\bm{y} - \bm{X}\widehat{\bm{\beta}})^{T} \bm{V}^{-1}(\bm{y} - \bm{X}\widehat{\bm{\beta}}),\quad \nu = n - p\; ,
\end{equation}
is the usual unbiased estimator for $\sigma^2$. Consider the linear hypothesis $H_{0}: \bm{R}\bm{\beta} =\bm{r}$ versus
$H_{1}: \bm{R}\bm{\beta} \not=\bm{r}$ for a given $r\times p$ matrix $\bm{R}$ of full row rank and a given $r\times 1$ vector $\bf{r}$. The corresponding $F$ statistic in model (\ref{E1}) is
\begin{equation}\label{E6}
F = \frac{(\bm{R}\widehat{\bm{\beta}} -\bm{r})^{T} (\bm{R} \bm{B} \bm{R}^{T})^{-1}(\bm{R}\widehat{\bm{\beta}} -\bm{r})}{r \widehat{\sigma}^2}
= \frac{(\bm{R}\widehat{\bm{\beta}} -\bm{r})^{T}(\bm{R} \bm{B} \bm{R}^{T})^{-1}(\bm{R}\widehat{\bm{\beta}} -\bm{r})}{(\bm{y} - \bm{X}\widehat{\bm{\beta}})^{T} \bm{V}^{-1}(\bm{y} - \bm{X}\widehat{\bm{\beta}})} \cdot \frac{\nu}{r}\; ,
\end{equation}
where
\begin{equation}\label{E7}
\mathop{\operator@font Cov}\nolimits(\widehat{\bm{\beta}}) = \sigma^2 \bm{B}, \quad \bm{B} = (\bm{X}^{T} \bm{V}^{-1} \bm{X})^{-1}\; .
\end{equation}
Under multivariate normality the statistic $F$ follows a
$F_{r,\nu}$ distribution provided $H_{0}$ holds true.

\subsection{Effect Size}

Suppose that we are interested in the effect of the independent variables represented by the model matrix $\bm{X}_{1}$, given the variables represented by $\bm{X}_{2}$. The corresponding linear hypothesis reads $\bm{R}_{1} \bm{\beta} = \bm{0}_{p_1}$ with $\bm{R}_{1} = (\bm{0}_{p_{1}}: \bm{I}_{p_{1}}:\bm{0}_{p_{1},p_{2}})$. Hence, by reasoning similar to Cohen (1988), an appropriate measure for the size of the effect based on the above $F$ statistic is provided by
\begin{equation}\label{E8}
f^2  = \frac{(\bm{R}_{1}\widehat{\bm{\beta}})^{T}(\bm{R}_{1} \bm{B}\bm{R}_{1}^{T})^{-1}(\bm{R}_{1}\widehat{\bm{\beta}})}{(\bm{y} - \bm{X}\widehat{\bm{\beta}})^{T} \bm{V}^{-1}(\bm{y} - \bm{X}\widehat{\bm{\beta}})}\;,
\end{equation}
thereby removing the factor $\nu/r$ from the $F$ statistic. This generalizes Cohen's $f^{2}$ in the sense that if $\bm{D} = \bm{0}_{q,q}$, then the unobservable random vector $\bm{u}$ equals $\bm{0}_{q}$ with probability one, the LMM reduces to the usual linear model of fixed effects only, and $f^2$ becomes identical to the measure provided by formula (9.2.1) in \citet{cohen1988statistical}.

We note that it is also possible to compute $f^2$ as
\begin{equation}\label{E9}
f^2 = \frac{R_{A,B}^2 - R_{A}^2}{1 - R_{A,B}^2}
\end{equation}
for appropriately defined coefficients of determination $R_{A,B}^{2}$ derived under the full model and $R_{A}^{2}$ derived under a reduced LMM assuming that variables represented by $\bm{X}_{1}$ are not present at all. Such a formula is the basis for the widely applied computational procedure suggested by \citet{selya2012practical}.

\subsection{Operational Effect Size}

Formula (\ref{E8}) for $f^2$ is only operational when $\bm{V}$ is completely known, a condition not met in practical applications. According to \citet{harville1977maximum} a simple LMM is the ordinary mixed and random effects ANOVA model, also referred to as traditional variance component model, see \citet[Chapt. 5]{christensen2019advanced}.
Under this model, the variance-covariance matrix of $\bm{y}$ depends on a total of $m$ variance components. The $n\times q$ matrix $\bm{Z}$ is partitioned as $\bm{Z} = (\bm{Z}_{1}: \cdots: \bm{Z}_{m-1})$, where each $n\times q_{i}$ matrix $\bm{Z}_{i}$ is the design matrix corresponding to a qualitative variable with a certain number of levels. In such a case one may assume
\begin{equation}
\bm{D} = \text{diag}\left[(\sigma_{1}^{2}/\sigma^2) \bm I_{q_{1}}, \ldots, (\sigma_{m-1}^{2}/\sigma^2) \bm I_{q_{m-1}}\right]
\end{equation}
where $\sigma^2 >0$ and $\sigma_{1}^{2}, \ldots \sigma_{m-1}^{2}\geq 0$ are $m$ unknown variance components. Then
\begin{equation}
\bm{V} = \bm{I}_{n} + \sum_{i=1}^{m-1} (\sigma_{i}^{2}/\sigma^2) \bm{Z}_{i} \bm{Z}_{i}^{T}\; .
\end{equation}
If $\sigma_{i}^2 = 0$ for some $i$, then the corresponding $q_{i}\times 1$ random effect vector $\bm{u}_{i}$ equals $\bm{0}_{q_{i}}$ with probability one. Such a model, and also more general ones, may be fitted in {\sf R} with the package
{\tt lme4}, see \cite{R:lme4}.  From the fitting procedure it is possible to obtain an estimate for the  variance-covariance matrix
\begin{equation}
\widehat{\mathop{\operator@font Cov}\nolimits}({\widehat{\bm{\beta}}}) = \widehat{\sigma}^2 \widehat{\bm{B}}
\end{equation}
of the estimated fixed effects parameter vector $\bm{\beta}$. Then a corresponding estimate for $f^2$ is given as
\begin{equation}\label{E10}
f^2  = \frac{1}{\nu}(\bm{R}_{1}\widehat{\bm{\beta}})^{T}(\bm{R}_{1}\widehat{\mathop{\operator@font Cov}\nolimits}({\widehat{\bm{\beta}}}) \bm{R}_{1}^{T})^{-1}(\bm{R}_{1}\widehat{\bm{\beta}})\; .
\end{equation}
The application of this formula is illustrated in the following section.

\begin{figure}[hbt]
\centering
\includegraphics[width=25pc,angle=0]{response_distribution.pdf}\\
  \caption{Distribution of response variable $Y$ over all $n=1000$ observations (histogram) and within two groups of sizes $n_{1} =687$ and $n_{2} = 313$ indicated by $0$ and $1$ (boxplots)}\label{F1}
\end{figure}

\section{Example}

In the following we discuss the performance of measure $f^2$ for an artificially generated data set of $n=1000$ observations intended to further illustrate some computational aspects. For the sake of simplicity our model consists of two independent variables $X_{1}$ (categorical/binary) and $X_{2}$ (quantitative) associated with fixed effects and one variable $Z$ (categorical) associated with random effects. However, the same principles apply when $X_{1}$, $X_{2}$, and $Z$ are extended to possible sets of variables containing more than one element. Our setting corresponds to the above mentioned variance component model with $p_{1}=p_{2}=1$ and $m=2$.

\subsection{Variable of Interest}

Figure \ref{F1} shows a discernible location difference with respect to the distribution of the response variable $Y$ in the two groups indicated by the binary variable $X_{1}$. The Welch two-sample $t$ statistic for the null hypothesis of no difference in group means reads  $|t|= 6.0751$ implying a highly significant result. A corresponding effect size measure is Cohen's $d$ which may be computed from {\sf R} package {\tt effectsize}, see \citet{R:effectsize}, as
$|d| = 0.4122$. From \citeauthor{cohen1988statistical}, values $|d|=0.2$, $|d|=0.5$ and $|d|=0.8$ indicate a small, medium and large effect, respectively.

\begin{figure}[hbt]
\centering
\includegraphics[width=25pc,angle=0]{scatter_with_groups.pdf}\\
  \caption{Scatter plot of $n=1000$ pairs $(X_{2}, Y)$ with group markings according to $X_{1}$}\label{F2}
\end{figure}

\subsection{Additional Fixed Effects}

From Figure \ref{F2} one may conclude, however, that the difference in the two groups may to some extent be explained by the variable $X_{2}$, since there is a tendency for larger values of $X_{2}$ to come along with larger values of $Y$ and observations from group 1 of variable $X_{1}$. Therefore one might be interested in the size of the effect of $X_{1}$ when $X_{2}$ is held constant. This can be achieved by considering a regression model with $X_{1}$ and $X_{2}$ as independent variables and deriving the measure $f^2$ as explained in \cite[Sect. 9]{cohen1988statistical}.  From package {\tt effectsize} one gets
$f^2 = 0.0017767$. Here, values $f^2=0.02$, $f^2=0.15$ and $f^2=0.35$ are supposed to indicate  a small, medium and large effect, respectively.

Recently, \citet{gross2023note} considered a generalized version $d_{\ast}$ of $d$ as an effect size measure for a binary variable $X_{1}$ given further variables. It may be computed from $f^2$ as
\begin{equation}
d_{\ast} = \sqrt{f^2 (n-2-w) \gamma}\; ,
\end{equation}
where in our analysis $w=1$ is the number of additional independent variables incorporated in the model  and $\sigma^2 \gamma$ is the variance of the regression coefficient for $X_{1}$. For our data $\gamma = 0.0065821$, yielding $d_{\ast} = 0.108$ and thus confirming a less than small effect.

\begin{figure}[hbt]
\centering
\includegraphics[width=25pc,angle=0]{group_means.pdf}\\
  \caption{Scatterplot of $(\#\{X_{1} =1\}, \overline{Y})$ for the 15 groups of $Z$, where
   $\#\{X_{1} =1\}$ denotes the number of observations with $X_{1}=1$}\label{F3}
\end{figure}

\subsection{Additional Fixed and Random Effects}
As a next step one may take the categorical variable $Z$, admitting 15 groups in our data set, into account. From Figure \ref{F3} it is seen that there is a tendency for $Z$ groups with larger means of the response variable $Y$ to contain less observations marked as 1 (referring to variable $X_{1}$) than $Z$ groups with smaller means of $Y$. However, since observations from group 1 were earlier seen to reveal on average larger $Y$ values than observations from group 0, holding $Z$ constant is expected to contribute in favor of an effect of $X_{1}$ again.

As noted before, the variable $Z$ is meant to be associated with a random effects vector $\boldsymbol{u}$. The corresponding
design matrix $\bm{Z}$ has 15 columns where each row contains 0's and a single 1, indicating the membership of the observations to the respective $Z$ variable group. There are two variance components, the overall $\sigma^2 > 0$ and the random effects variance denoted by $\sigma_{u}^2\geq 0$ here. The model may also be reparameterized via an unknown $k\geq 0$ by setting $\sigma_{u}^2 = k \sigma^2$. This gives variance-covariance matrix
\begin{equation}\label{E11}
\mathop{\operator@font Cov}\nolimits(y) = \sigma^2 \bm{V}, \quad \bm{V} = k \bm{Z}\bm{Z}^{T} + \bm{I}_{m}\; ,
\end{equation}
and the LMM reduces to the usual fixed effects regression model in case $k=0$. For $\bm{V}$ from (\ref{E11}) one may
compute the effect size $f^2$  for $X_{1}$ as a function of $k$ from formula (\ref{E8}) when all variables $X_{1}$, $X_{2}$ and $Z$ are included in the LMM formulation. Figure \ref{F4} shows $f^2$ for different choices of $k$. As expected, the effect size is larger when both,  $X_{2}$ and $Z$ are incorporated into the model compared to the case when only $X_{2}$ is employed, corresponding to the choice $k=0$.

\begin{figure}[hbt]
\centering
\includegraphics[width=25pc,angle=0]{effectsizes.pdf}\\
  \caption{Values of $f^2$ from a LMM fit depending on $k$}\label{F4}
\end{figure}

The operational version of $f^2$ from (\ref{E10}) can easily be computed from fitting the model by the function {\tt lmer} from package {\tt lme4} as
{\tt fit <- lmer(Y \~{} 1 + X1 + X2 + (1|Z))}. The estimated variance components are $\widehat{\sigma}^2 =  393.4455$ and $\sigma_{u}^2 = 180.4234$ giving
$\widehat{k} = \sigma_{u}^2/\widehat{\sigma}^2 = 0.4586$ as an estimation for $k$, see the dashed vertical line in Figure \ref{F4}. Then the vector $\widehat{\beta}$ is obtained from
{\tt fixef(fit)} and $\widehat{\mathop{\operator@font Cov}\nolimits}({\widehat{\bm{\beta}}})$ is obtained from {\tt vcov(fit)}. By using
$\bm{R}_{1} = (0,\, 1,\, 0)$ and $\nu= n-p = 997$, formula \ref{E10} yields
$f^2 = 0.0946626$ indicating a small but not medium effect size of $X_{1}$ when $X_{2}$ and $Z$ are held constant.

\subsection{Coefficient of Determination} Finally, we confirm that $f^2$ may also be obtained from formula (\ref{E9}). For this, we define
\begin{equation}\label{E12}
R_{AB}^2 = \frac{(r/\nu)  F}{1 + (r/\nu)  F}, \quad r=p-1,\, \nu = n-p\; ,
\end{equation}
which is the proposed $R^2$ for linear mixed  models by \citet[Eq. (19)]{edwards2008r2}. Here, $F$ is the $F$ statistic from (\ref{E6}) for testing the hypothesis $H_{0}:\bm{R}\bm\beta = \bm{0}_{p-1}$ with $\bm{R} = (\bm{0}_{p-1}: \bm{I}_{p-1})$. The relationship  between $F$ and $R_{AB}^2$ from (\ref{E12} may be established similar to  Example 4.8 from \citet{seber2003linear}. For our data $\bm{R} = (\bm{0}_{2}: \bm{I}_{2})$, $r= 2$, and $\nu = n-p= 997$. The measure $R_{A}^{2}$ is defined in the same way, but for a reduced model with all variables present except for $X_{1}$. For this reduced model we have $\bm{R} =(0, \, 1)$, $\bm{r} =(0)$, $r=1$, $\nu = n-p_{2}-1=998$. This results in
\begin{equation}
f^2 = \frac{R_{A,B}^2 - R_{A}^2}{1 - R_{A,B}^2}
= \frac{0.1539263 - 0.0738348}{1 - 0.1539263} = 0.0946626\; ,
\end{equation}
which is identical to the value computed above.