Friction as Contrast Mechanism in Heterodyne Force Microscopy

The nondestructive imaging of subsurface structures on the nanometer scale has been a long-standing desire in both science and industry. A few impressive images were published so far that demonstrate the general feasibility by combining ultrasound with an Atomic Force Microscope. From different excitation schemes, Heterodyne Force Microscopy seems to be the most promising candidate delivering the highest contrast and resolution. However, the physical contrast mechanism is unknown, thereby preventing any quantitative analysis of samples. Here we show that friction at material boundaries within the sample is responsible for the contrast formation. This result is obtained by performing a full quantitative analysis, in which we compare our experimentally observed contrasts with simulations and calculations. Surprisingly, we can rule out all other generally believed responsible mechanisms, like Rayleigh scattering, sample (visco)elasticity, damping of the ultrasonic tip motion, and ultrasound attenuation. Our analytical description paves the way for quantitative SubSurface-AFM imaging.


INTRODUCTION
Many fields of research are in need of a nondestructive way of imaging nanometer-sized subsurface features. To this end, ultrasound was combined with Atomic Force Microscopy to invent (Waveguide-) Ultrasound Force Microscopy [1,2] and Heterodyne Force Microscopy (HFM) [3]. HFM makes use of two ultrasound waves at slightly different frequencies, one of which is sent through the sample and the other through the cantilever. The mixed, heterodyne signal (amplitude and phase) at their frequency difference contains possible subsurface information at an experimentally accessible frequency [4].
Using HFM, subsurface images with remarkable contrast and resolution have been reported [3,[5][6][7][8][9][10][11][12][13][14][15][16][17][18], like the detection of 17.5 nm large gold nanoparticles buried at a depth of 500 nm in a polymer [10]. Surprisingly, the generated contrast clearly exceeds the background variations in these images, although the size of the nanoparticles is only a fraction of the sample thickness, and the lateral fingerprint on the surface (resolution) is equal to the diameter of the nanoparticles. Both observations are hard to understand, if one considers the wavelengths of the ultrasonic excitations, which is in the order of mm's and therefore much larger than both the size of the nanoparticles (nm's) and their depth below the surface (up to µm's). Unfortunately, none of the published HFM experiments provides quantitative information on the measured amplitude and phase range, on the applied contact force during the measurement, and on the precise excitation scheme in combination with the resonance frequencies of the cantilever.
To pave the way for quantitative subsurface measurements, it is of crucial importance to understand the physical contrast formation mechanism [19]. This requires a detailed, quantitative understanding of the ultrasound propagation within the sample [20], the cantilever dynamics [21][22][23][24][25], nonlinear mixing [26][27][28], the explicit excitation scheme, the resonance frequency spectrum of the cantilever [29,30], resonance frequency shifting [29], and the response to variations in the tip-sample interaction [29,31,32] that are determined by the local elasticity and adhesion of the sample. All these factors can significantly change the heterodyne signal leading to a measurable contrast. Published HFM experiments that provide (at least some) quantitative information are scarce [11] and the actual depth of the subsurface features is confirmed independently only in Ref. [9].
In this paper, we present a full quantitative analysis that addresses all thinkable, physical contrast mechanisms to explain our experimental observations on a well characterized sample. We show that Rayleigh scattering [20] would produce a contrast that is orders of magnitude smaller than in the experiment. By calculating the cantilever dynamics for different tip-sample interactions, we show that variations in sample elasticity indeed can lead to contrasts that are, in magnitude, comparable to the experiments. However, we can also rule out this mechanism, as the contrast is inverted with respect to the experimentally observed one. The only remaining possibility is dissipation! As we can also exclude tip damping and ultrasound attenuation, we finally conclude that friction at shaking nanoparticles is the responsible physical contrast mechanism. Additional evidence for this comes from an estimate of the involved energy dissipation.
Our analysis shows that the contrast strongly depends on the applied contact force and the precise ultrasonic excitation scheme with respect to the resonance frequencies (and their shifts) of the cantilever.

RESULTS
To enable a quantitative analysis of our measurements, we carefully prepared a sample with 20 nm large gold nanoparticles embedded 82 nm below the surface, see Fig.  1. The preparation as well as the independently determined characterization of the sample with AFM, RBS, and SEM is described in Supplementary Notes 1 and 2.
As the explicit excitation scheme is of crucial importance for the measured HFM contrast, Fig. 2 shows our particular experimental choice, called experimental scheme, with an off-off resonance excitation scheme (see Methods for the definition of their excitation schemes). Figure 3 shows the actual HFM experiment with simultaneously measured height, amplitude A dif f and phase φ dif f of the difference frequency f dif f for various contact forces F c . Feedback was performed in contact mode operation. The contact force F c is decreased from top to bottom: 163 nN, 115 nN, 67 nN, and 2.4 nN. The gold nanoparticles are visible in all channels at F c = 163 nN. The observed density of 1.2 particles/µm 2 fits the independently determined density. Most of the gold nanoparticles are still visible at F c = 115 nN, although the contrasts are significantly reduced. At lower forces, we do not (or just barely) detect any nanoparticles in any of the channels, which supports the RBS measurements that the gold nanoparticles are indeed fully buried under a 82 nm thick PVA layer. Considering the tip indentation depths (note that this is different from the total height variation, see information at the left side in Fig. 3) and the thickness of the PVA top layer, we have to conclude that the gold nanoparticles are only visible by poking hard enough into the sample, although it is striking that we see them at all in the height images. At F c = 2.4 nN, we start probing the attractive part of the tip-sample interaction and recognize that we have damaged the surface, while measuring earlier at higher contact forces.
At F c = 2.4 nN, both subsurface channels (amplitude and phase) show a clear correlation with the height. As the cantilever mainly probes the attractive part of the tip-sample interaction during an oscillation, the effective contact area of the tip depends on the height variations of the sample: it is much smaller on a mountain than in a valley. Adhesion is directly proportional to the contact area and a variation of it indeed leads to a variation in both the amplitude and the phase of the subsurface signal [26]. We conclude that variations in the adhesion do generate a contrast in the subsurface channels.
To quantify the contrasts of the gold nanoparticles in Fig. 3, we extract from cross sectional lines, as shown in Fig. 4, the average values above the nanoparticles for the different channels with respect to their background, see Tab. 1.
Let us first compare the experimental values with the expected contrast based on Rayleigh scattering [20], for which we have to normalize the amplitudes A dif f with respect to their background amplitudes A b . At F c = 163 nN, we measure a normalized amplitude contrast, A c , of -0.44 and a phase φ dif f of 7.2 degrees. At F c = 115 nN, the normalized amplitude contrast is -0.11 and the phase contrast is 2.9 degrees. Based on Rayleigh scattering, the expected normalized amplitude contrast is 10 −6 and the phase contrast is 0.1 millidegree for a gold particle with a diameter of 20 nm buried 50 nm deep under a polymer (PMMA) [20]. As the experimentally observed normalized amplitude contrast is 5 orders of magnitude larger (and the phase contrast 4 orders of magnitude) than the theoretically predicted ones, we have to conclude that From left to right measured simultaneously: the height and both the amplitude A dif f and the phase φ dif f of the difference frequency. The contact force Fc as well as the resulting average indentation into the sample are indicated at the left in the height images. The gold nanoparticles are only visible at a contact force of 163 nN and 115 nN. At these forces, they are not only visible in the subsurface channels, but also in the height image. We 'loose' the nanoparticles in all three channels with decreasing force. At a Fc = 2.4 nN, we observe that we damaged the surface, while measuring at higher forces. All images within one channel do have the same (color) range such that the contrast for different contact forces can be compared directly. We provide typical cross sections with absolute values of the three channels at the positions of the nanoparticles in Fig. 4.
Rayleigh scattering does not form a major contribution to the physical contrast mechanism (at least not at MHz frequencies).
Recently, it was elucidated how the heterodyne signal is generated: its magnitude strongly depends on both the applied contact force and the specific characteristics of the tip-sample interaction [21,22,26]. In Supplementary Notes 3 and 4 we show, both experimentally and analytically, that the heterodyne signal depends on the elastic properties of the sample, which is characterized by   its Young's modulus E. For sufficiently soft samples, the amplitude A dif f is proportional to E. Let us, in the following, consider elasticity variations in the sample, due to the presence of the nanoparticles, as a possible contrast mechanism.
From an analytical 1D model, we estimate that the Young's modulus above a gold nanoparticle is ∼ 10% higher than the Young's modulus of PVA, which is 2.4 GPa, see Supplementary Note 5. To determine the contrast formation based on these elasticity variations, we numerically calculated the motion of the cantilever for different tip-sample interactions using the method outlined in [21]. The result is shown in Fig. 5, in which we, for reasons of clarity, only show the approach curves. To receive an upper bound on the contrast and to elucidate the contrast formation effect on the basis of small elasticity variations, we consider Young's moduli between 2 and 6 GPa. As the specific vibration spectrum of the cantilever has great influence on the results, we first matched the spectrum used in the calculations to that of our experiment, see Supplementary Note 6. We call the particular off-off resonance excitation scheme that we used in this experiment (see Fig. 2), experimental excitation. The graphical result, see Fig. 5, show the corresponding tip-sample interactions and, as a function of the applied contact force, the indentations as well as the amplitudes A dif f and phases φ dif f of the heterodyne signal at the difference frequency. The contrasts at a certain contact force can now be evaluated from the difference in the signals stemming from different elasticities (colors in the graphs). The indentation contrast decreases with decreasing contact force. The amplitude contrast stays almost constant over a large range (and even increases slightly), before it collapses to zero at very small contact forces. The phase contrast strongly depends on the specific excitation scheme, but always collapses to zero at very small contact forces. The extracted height, amplitude and phase values are listed in Tab. I. In addition, to elucidate the effect of different ultrasonic excitation schemes, we also considered an off-off resonance excitation, in which both ultrasound signals are midway between two resonance frequencies, as well as an off-on resonance excitation, see Supplementary Note 7. These results are, in addition, tabulated in Tab. I for comparison.
The experimental scheme with 2.4 GPa (PVA) to 2.6 GPa (effective elasticity above the nanoparticles, see Supplementary Note 5) perfectly reflects both the sample and the measurement conditions. To receive clear upper bounds, we determined further all excitations schemes from the differences between a sample with 2 GPa and 6 GPa. Starting with the height contrast, we find comparable values between the experiment and the calculated excitation schemes, except for the experimental scheme 2.4 → 2.6 GPa. The decrease in height contrast for smaller contact forces F c is reproduced for all cases. Considering the amplitude contrast ∆A dif f , the absolute values in the experiment are up to 100 times larger than  we calculated the tip-sample interaction and, as a function of the applied contact force, the corresponding sample indentation as well as the amplitude A dif f and phase φ dif f of the heterodyne signal for different sample elasticities: 2 GPa (black), 3 GPa (red), 4 GPa (magenta), 5 GPa (green), and 6 GPa (blue). The inset in the lower left panel shows A dif f for 6 GPa plotted as a function of the height of the cantilever's base, z b , such that a comparison becomes possible with other calculations [21,22,26].   the calculated ones. One notices three striking issues when performing a more detailed comparison. Firstly, the values of the off-on resonance case are significantly lower than most other values. This is due to its particular excitation scheme (see Supplementary Note 7), in which the ultrasonic tip amplitude significantly decreases when the cantilever gets into contact with the surface, due to the related frequency shift of the 4 th resonance. The size of this shift and, therefore, also of the amplitude reduction of the ultrasonic tip vibration, increases both with increasing contact force and with sample stiffness. This excitation scheme with its particular behavior is special, as the frequency shift acts as an amplifier/attenuator to the measured signal. Secondly, in contrast to the experiment, both the off-off resonance case and the experimental schemes show a larger contrast at lower force. Although this already indicates a problem, the most striking issue is the sign of the contrast, which is inverted in comparison with the experiment! In any case, the (visco)elasticity above the nanoparticle is for sure increased, which theoretically leads to a higher amplitude A dif f (see Fig. 5 and Supplementary Note 4) and, therefore, to a positive amplitude contrast ∆A dif f . We have to conclude that, although elasticity variations produce a contrast with similar magnitude than in the experiments, they cannot explain the inverted contrast. Consequently, a different physical mechanism must be present.
Please note that the amplitude contrast inversion of ∆A dif f in the off-on resonance case is due to its particular excitation scheme with the frequency shift of the 4 th mode. Above the nanoparticle, the amplitude reduction of the ultrasonic tip vibration A t is significantly larger than the reduction on the PVA without nanoparticles (see Supplementary Note 7). This indicates the importance of the precise excitation scheme and the spectrum of the cantilever for each published HFM measurement, in order to understand it quantitatively.
For the sake of completeness, we shortly turn our attention also to the phase behavior. The magnitude of the experimentally observed phase contrast ∆φ dif f is only comparable to the special case of the off-on resonance excitation scheme. The large phase shift in this scheme is due to the frequency shift of the 4 th resonance: the particular off-on resonance excitation scheme makes the tip vibration especially sensitive to phase changes based on frequency shifts. Although much smaller in magnitude, a similar argument holds also for the phase shifts in the offoff resonance and experimental excitation schemes. Since the ultrasonic tip excitation in the experimental scheme is closer to the 4 th resonance frequency of the cantilever, we observe a larger phase contrast than in the off-off resonance scheme where the excitation of the tip is midway between resonance frequencies.
Summarizing this part, we conclude that the contrast from (small) variations in the sample elasticity results in a much larger contrast than Rayleigh scattering: the order of magnitude is comparable to the experiments. However, variations in sample elasticity cannot be the physical contrast mechanism in our HFM experiment, as it would imply an opposite sign.

DISCUSSION
Ruling out both variations in the tip-sample interaction (elasticity and adhesion) and Rayleigh scattering, the remaining physical contrast mechanism must lead to a significant reduction of the tip amplitude A t or the sample amplitude A s above the nanoparticles, as . These reductions can be described as tip or sample damping. Tip damping can also be excluded, as it has been surprisingly shown that A t keeps 99.7% of its amplitude at a contact force of 25 nN even on a hard sample like Si [22]. Please note that the damping of the resonance frequencies of a cantilever that is in contact with a sample, is generally assumed to be directly proportional to the Young's modulus of the sample [37]. Without significant tip damping, the contrast must be due to a reduction in the sample amplitude. Since a reduction of A s is expected to occur also on the polymer without nanoparticles, and since A dif f is larger above the nanoparticle due to the increase in the effective Young's modulus, we need a mechanism that leads to a strong decrease of A s only above the nanoparticle to overcompensate the increase in A dif f such that it effectively leads to a contrast inversion (holes in A dif f , see Fig. 3).
Let us start with a possible vertical motion of the nanoparticles in the polymer matrix. At low ultrasonic sample frequencies, this motion is surely in phase with the excitation. However, if the ultrasonic excitation is above the resonance frequency of the system "nanoparticle in polymer", the motion will be out of phase leading to a significant reduction of A s only above the nanoparticles. The problem is, however, that the sample excitation is at 2.5 MHz and that we estimate the resonance frequency of the "nanoparticle in polymer" system to be ∼ 2.2 GHz (see Supplementary Note 8). The nanoparticles should, therefore, simply follow the ultrasonic displacements of the polymer.
Another mechanism worth considering is sample damping (reduction of A s ) by energy dissipation at the nanoparticles. Next to contrast formation based on attenuation or friction, a temperature effect might additionally enhance the contrast, especially if the elasticity of the polymer would have a strong temperature dependence. Therefore, we estimate the energy dissipation from the experiment. We determine the sample amplitude A s (far away from the nanoparticle) in analogy to the method described in [21] (see also Methods). With A s ∼ 0.22 nm at F c = 163 nN, we need a reduction of ∼ 41% to explain the observed contrast. A s is ∼ 0.29 nm at F c = 115 nm, which corresponds to a reduction of ∼ 13%. Both estimations deliver a similar value: 0.83 and 0.37 pW, respectively. This breaks down to an energy dissipation at the nanoparticles of less than 2.07 eV/oscillation of the ultrasonic sample excitation. This value is so small that we can rule out also any temperature effects. The only remaining physical mechanism that might cause this energy dissipation is ultrasound attenuation within the nanoparticles as well as friction at the interface between the nanoparticles and the polymer.
The ultrasound attenuation for gold is ∼ 150 times smaller than the attenuation for PVA. Therefore the total energy dissipation is less at the positions measured above the nanoparticles than at the positions far away from them. This effect results, in comparison to the experiment, again in a wrong sign of the contrast, as A s should be larger above the nanoparticles. We estimate this resulting energy 'gain' based on a smaller ultrasound attenuation at the nanoparticles to be 0.45 eV/oscillation. The dissipation that causes the contrast, must be increased with this value to overcompensate it and lead to contrast inversion.
This means that we are left with friction at the interface between the nanoparticles and the PVA. Due to a weak (chemical) bonding between the gold and the PVA, the nanoparticles might (slightly) slip instead of following all displacements of the PVA. One might even consider a small cavity around the nanoparticles such that they are shaken up and down. Both effects would lead to a significant amount of friction at the interface. Considering shaking nanoparticles, we are able to explain our observed contrast with a total energy dissipation of 2.52 eV/oscillation at the nanoparticles. This value is comparable (and definitively in the right order of magnitude) with the energy dissipation derived from atomic scale friction experiments of a sharp tip in contact with a surface [38]. Note that the tip radius in these experiments is comparable to the radius of the nanoparticles.
Pinpointing the physical mechanism to friction at shaking nanoparticles, we can consider the consequences for the lateral resolution. If one assumes that the propagation in amplitude reduction obeys a scattering-like behavior, the 'fingerprints' of the nanoparticles at the surface should show a significantly larger diameter than the diameter of the nanoparticles. Moreover, as we are measuring in near-field, the size of the 'fingerprints' should be in the order of the depth of the nanoparticles. The deeper the nanoparticle is, the larger should be its image at the surface. These considerations stand in clear contrast to experimental observations: nanoparticles with a diameter of ∼ 17.5 nm, buried 500 nm deep, are imaged with a diameter of only 20 nm [10], and the imaged fingerprint is even decreasing with increasing depth of the nanoparticles [9]. A solution to this might be found by considering a combination of a stress field that is introduced on the nanoparticle by the indenting tip [31], a resulting shaking that is no longer parallel to the initial ultrasonic displacements of the PVA, and a highly anisotropic propagation of the amplitude reduction such that there is a significant enhancement in the direction of the shaking movement. wrote the manuscript together, which was carefully read and improved by all authors.

METHODS
As a quantitative analysis of the contrast mechanism is impossible without a well-characterized sample, we carefully prepared a stack consisting of the following layers (from bottom to top): a Si wafer with native oxide, a ∼ 97 nm thick PMMA layer, a 30 nm thick PVA layer with embedded gold nanoparticles (diameter 20 nm), and a 82 nm thick PVA top layer. The density of the gold nanoparticles was determined via AFM and SEM to be 0.7 ± 0.6 particles/µm 2 . The precise sample preparation as well as its detailed characterization, in which we even determined the depth of the Au nanoparticles with an independent measurement based on Rutherford backscattering, is described in detail in Supplementary Note 1 and 2.
In our HFM experiment, we chose the ultrasonic excitation frequencies of both the tip and the sample as well as the difference frequency off resonance, i.e. not on (or within the width) of a resonance peak of the cantilever. We call this excitation scheme off-off resonance. The first on/off indication describes whether f dif f (heterodyne signal) is tuned to a resonance frequency of the cantilever, whereas the second on/off indication describes whether f t (ultrasonic tip excitation) is tuned to a resonance. This leads to four different excitation schemes, of which we evaluate also the off-on scheme in more detail in Supplementary Note 7. Figure 2 shows the excitation scheme and the vibration spectrum of the free hanging cantilever, of which we calibrated the spring constant to be 2.7 N/m using the thermal noise method [33]. Using the method described in [22,26], we determined the ultrasonic tip amplitude to be A t = 1.34 nm and the ultrasonic sample amplitude to be A s = 0.37 nm.

Sample Preparation
Inspired by the sample with buried gold nanoparticles of Shekhawat and Dravid [1], we set out to produce comparable ones. We decided to use gold nanoparticles with a diameter of 20 nm (±10%), which we got from BBI Solutions [9]. A schematic cross section of the final sample that we used in the current study, is shown in Fig. 1 in the main text. In the following, we describe important issues of the sample preparation and provide the recipe. As a substrate, we used a freshly with acetone cleaned Silicon (100) wafer that was covered with a native oxide. The polymer layers (including the suspension with the nanoparticles) were deposited by means of a spin coater, see also the recipe below. We decided to use two different polymers: polymethylmethacrylaat (PMMA) with a degree of polymerization of 970 and polyvinyl alcohol (PVA) with a degree of polymerization of 2700. The degree of polymerization is the number of monomers in the molecule and it characterizes the length of a single polymer molecule. This information is important, as the material properties of the polymer layers strongly depend on the molecule length. As we faced some problems with clustering of the nanoparticles as well as with their density, we describe these issues shortly in the following. Our first attempt to create a layer of gold nanoparticles on top of a spin coated PMMA layer was to let a suspension of pure (Milli-Q) water with gold nanoparticles evaporate at ambient conditions. This led to large "mountains" of clustered nanoparticles, which we measured with an Atomic Force Microscope (AFM). In our second attempt, we tried to embed the gold nanoparticles within a PVA layer by dissolving them in a PVA solution before spin coating on top of the PMMA. The gold nanoparticles did stick out with their "heads" just above the PVA layer, with which they were simultaneously spin coated, such that we easily could verify the density, again, with AFM. This approach resulted in an (for our research) unsuitable low density of nanoparticles of less than 0.1 nanoparticle/µm 2 . By increasing the concentration of the nanoparticles in the PVA solution, we were able to increase the density to 0.7 ± 0.6 nanoparticle/µm 2 . We derived this distribution from AFM measurements, see Fig. 6A. Finally, we buried the nanoparticles by spinning another PVA layer on top of this structure. As we considered that the solvent, which is present while spinning the additional PVA layer, might (partially) dissolve the thin nanoparticles/PVA layer that is to be buried, leading to a possible redistribution of the nanoparticles, we counterchecked the density with a Secondary Electron Microscope (SEM) on the final sample with the top PVA layer. Due to the different electron emissions between gold and PVA, the SEM is capable of imaging the nanoparticles, even if they are buried under a 82 nm thick PVA layer, see Fig. 6B.
The final recipe for the sample production is as follows: 1. solution: 30 mg PMMA / mL Toluene This results in a ∼ 97 nm thick PMMA layer, which was confirmed independently with an AFM measurement [10].
2. solution: 250 µL of 2 mg PVA / mL water + 750 µL suspension of pure water and gold nanoparticles This leads to a PVA layer with embedded gold nanoparticles with a diameter of 20 nm. The thickness of this PVA layer is less than 30 nm (∼ 10 nm), as we verified with AFM that the "heads" of the gold nanoparticles are sticking out. After burying this PVA layer with the top PVA layer, we find an effective thickness of 30 nm for this layer that contains the nanoparticles.
3. solution: 2mg PVA / mL water This step leads to a ∼ 82 nm thick PVA layer.
Each step in the recipe represents an individual spin coating procedure. In each spin coating step, a droplet of the corresponding solution was put onto the sample by means of a pipet before the spin coater started to rotate for 5 s at 2000 rpm immediately followed by 90 s at 4000 rpm.
Although we did not apply explicitly a curing (baking) step of the final sample after the preparation, the complete sample was baked for approximately 3 minutes at ∼ 140 0 C to glue it with crystalbond 509 onto the ultrasonic transducer of the sample. This was always done within 24 hours after the spin coating procedure. We assume that, during this baking procedure, most of the remaining solvents in the sample were evaporated.

Independent Verification of the Nanoparticle Depth
In order to quantify HFM experiments, it is of great importance to have a well defined sample, in which the depth of the subsurface particles (or features) is counterchecked with an independent technique. To this end, we performed a Rutherford Backscattering Spectrometry (RBS) measurement on the sample that we used for experiments, to quantify the exact depth of the gold nanoparticles as well as the thickness of the individual layers. To deduce quantitative data from an RBS measurement, it is necessary to perform a simulation [12]. Figure 7 shows both the RBS measurement (black) and the corresponding result of the simulation (red). The surface channels of the different elements in our sample (Carbon, Oxygen, Silicon, and Gold) are indicated in blue. Although almost at the detection limit of the RBS setup, the inset clearly shows a signal obtained from the buried gold nanoparticles: it is a sharp distribution, which indicates a well defined depth of the nanoparticles, with a clear shift away from the surface channel of Au, from which we can determine the thickness of the top PVA layer. The RBS spectrum in Fig. 7 shows that Si is present just below the sample's surface, see the rise (and the tiny  plateau) in the spectrum almost at the Si surface channel as well as layer 2 in Tab. II. We can explain this with the presence of air bubbles in our sample and/or holes in some of the spin coated polymer layers. As a consequence, the best simulation result contains 4 layers on top of the Si wafer (see Tab. II). The combined thickness of layers 1 and 2 is 82 nm. Therefore, the gold nanoparticles are buried approximatly 82 nm below the surface. The thickness of the underlying PMMA layer is approximately 97 nm. We verified the total thickness of 209 nm by scanning over scratches on the sample with an AFM. The minimum thickness that we found in all AFM heightlines is ∼ 250 nm, which confirms the RBS analysis.
It is striking that the sample contains more C than expected, but less O. From the simulation, we find the following composition for the PVA: C 2 H 1.6 O 0.2 , which has to be compared to C 2 H 4 O 1 . The lack of oxygen can be explained either by the formation of water during the baking procedure at ∼ 140 degrees 0 C after the spin-coating or by a decomposition of the polymer layers during the RBS measurements (a clear spot on the sample surface was visible after the experiment). For the PMMA layer we find a composition of

Experimental Dependence of the Difference Frequency Amplitude A dif f on the Sample Elasticity
To experimentally address the dependence of the amplitude A dif f of the heterodyne signal at the difference frequency on the elasticity of the sample, which is characterized by its Young's modulus E, we present results for the difference frequency generation on both a soft sample (∼ 97 nm thick PMMA, E ∼ 2.4 GPa) and a hard sample (Si(100) wafer, E ∼ 179 GPa). The HFM experiment was performed with a similar cantilever as described in the main text.
We obtained the Young's modulus on PMMA by fitting an experimentally obtained tip-sample interaction F ts with the Derjaguin-Muller-Toporov (DMT-) model [6]. A parameter called λ, which is related to the elasticities of the tip and the sample, is usually used to differentiate between the applicabilities of different models that describe the tip-sample interaction [7]. As λ = 0.63 in our case, one should use the Maugis-Dugdale model [8]. Nevertheless, our approach with the DMT-model is fully justified, as we have demonstrated in [3] that it does not matter at all for the numerical simulations which of the models describes the tip-sample interaction, as long as the fit perfectly matches the (experimentally obtained) tip-sample interaction. The only thing that matters is the particular shape (form) of F ts (z) and not the model that is used to describe this particular interaction.
The cantilever has a spring constant of 2.0 ± 0.4 N/m, which was calibrated using the thermal noise method [5]. We applied an off-off resonance excitation scheme with an ultrasonic tip frequency of 2.870 MHz and an ultrasonic sample frequency of 2.871 MHz leading to a heterodyne signal at a difference frequency f dif f of 1 kHz. The ultrasonic vibration amplitudes of both the tip A t and the sample A s were slightly different for the two experiments: A t = 0.94 nm and A s = 0.32 nm on Si, whereas A t = 1.23 nm and A s = 0.18 nm on PMMA. The tip amplitudes were determined using the procedure outlined in [3] and below we describe how we determined the sample amplitudes from the measurements. We measured the amplitude A dif f of the difference frequency as a function of the cantilever's base position z b on both the Si and the PMMA layer. z b is defined such that z b = 0, if the deflection δ = 0 during the approach cycle of the cantilever to the surface. This is exactly the point, at which the effective interaction on the tip changes sign from an attractive interaction to a repulsive interaction. Figure 8 shows the amplitude A dif f of the heterodyne signal as well as the corresponding deflection of the cantilever as a function of the cantilever's base position z b on both a Si sample (red) and a ∼ 97 nm thick PMMA layer (black). To determine the ultrasonic vibration amplitude A s of the sample, we can estimate A s , for the case of Si, from the height of the plateau, using the method described in [3], to be 0.32 nm. Without the existence of a clear plateau in the PMMA case (note the slight decrease of A dif f for negative z b ), we instead use the maximum amplitude of the As t , see below and [4], and as the vibration amplitudes are slightly different for the measurements on Si and PMMA, one has to multiply the amplitude A dif f for the PMMA case with a scale factor of 1.78 to accommodate for a valid comparison. Please note that, even with this correction factor, A dif f is significantly larger on the hard Si surface (179 GPa) than on the soft ∼ 97 nm thick PMMA layer (2.4 GPa). The lower panel shows the corresponding deflection δ of the cantilever.
difference frequency A dif f (0.18 nm) for the estimation and find A s = 0.18 nm.
To enable a valid comparison between the measurements on the two different samples, one has to multiply the amplitude A dif f for the PMMA case with a correction factor of 1.78, as A dif f ∝ A s · A t / A 2 s + A 2 t (see below and [4]) and as the vibration amplitudes are slightly different for the measurements on the different samples. Taking this correction factor into account, one still observes that, for the same contact force, A dif f is significantly larger on the hard Si surface than on the soft ∼ 97 nm thick PMMA layer. The peaks in the attractive regime are larger for the soft PMMA sample, because the adhesion is larger on the PMMA sample than on the Si sample (please note the difference in deflection in the attractive part of the tip-sample interaction). Thus we conclude that the amplitude A dif f significantly depends on the elasticity of the sample and increases with increasing Young's modulus E.

Analytical Dependence of the Difference Frequency Amplitude A dif f on the Sample Elasticity
Recently, an analytical theory has been developed that completely describes the generation of the heterodyne signal at the difference frequency for HFM experiments [4]. The signal is characterized by the following analytical expressions: (1) , in which A dif f and φ dif f are the amplitude and the phase, respectively, of the signal at the difference frequency, and A s and A t are the ultrasonic vibration amplitudes of the sample and the tip with corresponding phases φ s and φ t . |H −1 (ω dif f )| represents the absolute value of the inverse transfer function and its corresponding phase shift Λ, F ts is the tip-sample interaction as a function of the tip-sample distance z, z b is the position of the cantilever's base, and δ is the deflection of the cantilever. The integrals I 1 and I 2 completely determine the generation of the signal at the difference frequency and they both depend on the particular tip-sample interaction F ts . As the tip-sample interaction in the experiment can be best described by the DMT-model [3,6], F ts can be expressed by , in which R is the radius of the cantilever's tip, H the Hamaker constant, a 0 the distance at which the repulsive part of the tip-sample interaction is first felt by the cantilever (∼ at the minimum of F ts ), and E f is an effective Young's modulus describing the effective tip-sample stiffness. This effective Young's modulus E f is determined by the elasticities (E t and E) as well as the Poisson ratio's (µ t and µ) of the cantilever and the sample, respectively: Since we probe our final sample that consists of several polymers layers, of which E ∼ 2.4 GPa, with a hard Silicon cantilever with E t ∼ 179 GPa, we can neglect (1 − ν 2 t )/E t and receive that the effective elasticity E f is directly proportional to the elasticity E of the sample. The repulsive part of the tip-sample interaction, see Eq. 4, is, therefore, also directly proportional to the elasticity E of the sample. As a consequence, this is valid also for the integrals I 1 and I 2 described by Eqs. 2 and 3. Using these proportionality relations in Eq. 1, we find a simple expression for the elasticity dependence of the amplitude A dif f of the heterodyne signal: , in which γ is a complex constant. If the cantilever is completely in the Hertzian contact regime (z < a 0 ) during its oscillation, gamma can be written as , in which α is the normalized indentation given by: We can evaluate an lower estimate for γ by setting the normalized indentation α = 1 and noticing that for smaller α, the integral in the expression for γ would become smaller, and γ, therefore, larger. Using the ultrasonic amplitudes of both the tip and the sample, we appraise A 2 s + A 2 t = 1.39 nm. For the tip radius we assume R = 5 nm. The inverse transfer function |H −1 (ω dif f )| can be derived as described in [4], in which we take a spring constant of 2.5 N/m and set Λ to zero. This leads to the following estimates for γ and A dif f : , in which E has to be inserted in GPa. Equation 9 describes an analytical dependence of A dif f on the sample elasticity E. For soft samples, in which E is smaller than 0.5 GPa, A dif f is approximately proportional to E. Therefore, we also expect analytically that a harder sample results in a higher amplitude A dif f of the heterodyne signal, especially above the nanoparticles, where the effective elasticity is slightly increased with respect to the soft polymer. On very hard samples, with E 0.5 GPa, A dif f approaches a constant value and becomes independent of E. Please note that we have neglected the influence of the elasticity on both the deflection of the cantilever and the transfer function of the cantilever. However, this is a valid approximation, as we never saw a decrease in the amplitude A dif f of the difference frequency at a given contact force while the elasticity E of the sample was increased.

Effective Sample Elasticity above the Nanoparticles
In this section, we derive an upper bound for the effective sample elasticity, measured at the sample surface, that is increased by the presence of the buried nanoparticles in the polymer. Figure 9 shows a schematic cross section of the sample, which consists of a PVA layer (PVA), the gold nanoparticles (Au), and a PMMA layer (PMMA). The Au is buried at the depth d, and has a radius R. The total thickness of the sample is denoted with t. The sample is compressed by a stress σ, which is equal to the force F per unit area A. From linear elasticity theory, we know that an applied external stress is negatively proportional to the relative change in thickness, in which the proportionality factor is given by the Young's modulus E of the material. For our sample this reduces to the following equations: , in which E P V A , E Au , and E P M M A are the Young's moduli of PVA, gold, and PMMA, respectively, δd is the variation in depth of the nanoparticle, δR is the variation in radius of the nanoparticle, and δt is the variation in thickness of the sample.
It is straightforward to derive the solutions for δd, δR, and δt, from Eq. 10: If one introduces an effective Young's modulus E ef f , the complete sample with all three layers can be regarded also as a sample consisting of one layer with a thickness t of an isotropic material such that By substituting Eq. 13 in Eq. 14, we find an expression for the effective Young's modulus E ef f : Let us now discuss the two limits of this equation. Firstly, if the diameter 2R of the nanoparticle is equal to the thickness t of the sample (and thus d = 0), we find that E ef f = E Au . Secondly, if the radius R of the nanoparticle is equal to zero and the sample is infinitely thick (t d), we find that E ef f = E P M M A . Thirdly, if the radius R of the nanoparticle is equal to zero and d = t, we find that E ef f = E P V A . These results reflect correct expectations, as the sample consists only of Au in the first case, only of PMMA in the second case, and only of PVA in the third case.
Equation 15 provides an upper bound on the elasticity on the surface above a nanoparticle. In reality, the variation in elasticity due to a nanoparticle should be derived from a 3D calculation, as the stress is spread out also laterally through the sample [13]. As a consequence, the rise in elasticity caused by the presence of the nanoparticle decreases with increasing depth of the nanoparticle. This effect is comparable to a stone underneath a pillow: if one just touches the pillow, the stone is not felt, but if one pushes harder into the pillow, the presence of the stone is clearly noticed.
Let us now calculate the expected effective elasticity increase for our samples. PMMA and PVA, both have a similar Young's modulus: E P M M A ∼ E P V A = 2.4 GPa. Under this assumption, Eq. 15 reduces to: We assume that the Young's modulus of the gold nanoparticle is equal to that of bulk gold, which is 78 GPa [15] and consider the total thickness to be t = 209 nm. For the radius, we take R = 10 nm of the gold nanoparticles, as this is the average of their radii distribution [9,10]. Using this value for the radius R, we find the effective Young's modulus E ef f to be equal to 2.65 GPa. Therefore, the surface directly above the nanoparticle has (at maximum) a 10% higher Young's modulus than that of the bulk polymer of 2.4 GPa.

Setting up the Numerical Calculations
As the amplitude and phase contrast highly depend on both the exact excitation scheme and precise resonance frequency spectrum of the cantilever, which can even result in a contrast inversion, it is of uttermost importance to match the spectrum of the cantilever in the numerical calculation (numerical cantilever) to the spectrum of the cantilever used in the experiment (experimental cantilever). In this section, we describe the matching procedure.
On the basis of the resonance frequencies f num i of the numerical cantilever and the corresponding resonance frequencies f exp i of the experimental cantilever, we defined a normalized, relative error e i for each resonance frequency: We took into account the first 5 modes of the cantilever and used the average e i as a measure for the quality of our fit. We optimized the fit by varying the elasticity E t of the cantilever, the length L of the cantilever, the tip mass m e , the moment of inertia I e of the tip, and the density ρ s of the cantilever. As a best fit, with an average error of 1.397%, the cantilever is described by following parameters: E t = 222 GPa, L = 207 µm, m e = 5.76 · 10 −15 kg, I e = 3.51 · 10 −22 kg m 2 , and ρ s = 3207 kg m −3 . We did not fit the width and the thickness of the cantilever. Instead we have chosen them to be 20 µm and 2.7 µm, respectively, to set the spring constant of the numerical cantilever to 2.5 N/m such that it is comparable to the spring constant of 2.7 ± 0.4 N/m of the experimental cantilever.
The Q-factors that describe the widths of the resonance peaks, were chosen such that the widths of the resonance peaks match between the numerical and the experimental cantilever.  10: The vibration spectrum of both the experimental and the numerical cantilever. The bottom two panels show the amplitude and the phase of the experimental cantilever. The phase extremely decreases almost linear with the frequency (notice the phase change from 0 0 to −4200 0 ): this is due to the phase change in the fixed cables, which deliver the electronic drive signal to the cantilever. The red lines indicate the experimentally determined resonance frequencies, which are indicated at the top together with their corresponding Q-factors that describe the widths of the resonance peaks. The second panel from the top shows the spectrum of the numerical cantilever for comparison. The blue lines indicate the frequencies of this particular excitation scheme with ft = 2.50 MHz (cantilever) and fs = 2.52 MHz (sample), as well as the difference frequency f dif f = 20 kHz. Fig. 10 shows the vibration spectrum of both the experimental cantilever and the numerical cantilever. The bottom two panels show the amplitude and the phase of the experimental cantilever. The phase extremely decreases almost linear with the frequency (notice the phase change from 0 0 to −4200 0 ): this is due to the phase change in the fixed cables, which deliver the electronic drive signal to the cantilever. The red lines indicate the experimentally determined resonance frequencies, which are indicated at the top together with their corresponding Q-factors that describe the widths of the resonance peaks. The second panel from the top shows the spectrum of the numerical cantilever. The blue lines indicate the frequencies of this particular excitation scheme with f t = 2.50 MHz (cantilever) and f s = 2.52 MHz (sample), as well as the difference frequency f dif f = 20 kHz.

Complete Overview of the Results of the Numerical Calculations
In the main text, we describe the results of three different schemes for the ultrasonic excitations: off-off resonance, in which the ultrasonic excitations are chosen halfway between the 3 rd and 4 th resonance of the cantilever; off-on resonance, in which the ultrasonic excitation frequencies are on the 4 th resonance of the cantilever; experimental excitation, in which the ultrasonic excitation frequencies are equally far away from the nearest resonance frequency as in the experiment. In this section, we present the full numerical results of both the off-off resonance and the off-on resonance scheme. FIG. 11: The top panel shows the vibration spectrum in the off-off resonance scheme: blue lines indicate the excitation frequencies, whereas red lines indicate the resonance frequencies of the cantilever. The bottom three panels show the numerical results of the off-off resonance scheme for different sample elasticities: 2 GPa (black), 3 GPa (red), 4 GPa (magenta), 5 GPa (green), and 6 GPa (blue). The panels display from left to right the indentation (∼ inverted height), the amplitude A dif f , and the phase φ dif f . We observe an instability in the cantilever's motion while indenting into the sample (see jumps at ∼ 70 nN). Note that, after the instabilities, we receive smooth motions of the cantilever and the phase differences approach values that correspond to the numerical error of the lock-in. Figure 11 shows the results for the off-off resonance excitation scheme. The top panel shows the vibration spectrum: blue lines indicate the excitation frequencies and red lines the resonance frequencies. The bottom three panels show the numerical results for different sample elasticities: 2 GPa (black), 3 GPa (red), 4 GPa (magenta), 5 GPa (green), and 6 GPa (blue). The panels depict from left to right the indentation (∼ inverted height), the amplitude A dif f , and the phase φ dif f . The cantilever's motion is unstable for some contact forces while indenting in the sample (see jumps at ∼ 70 nN). Note that after the instability the motion is stable again, A dif f is smooth, we receive smooth motions of the cantilever, and the phase differences approach values that correspond to the numerical error of the lock-in. Figure 12 shows the results for the off-on resonance excitation scheme. The top panel shows the vibration spectrum: blue lines indicate the excitation frequencies and red lines the resonance frequencies. The bottom three panels show the numerical results for different sample elasticities: 2 GPa (black), 3 GPa (red), 4 GPa (magenta), 5 GPa (green), and 6 GPa (blue). The panels depict from left to right the indentation (∼ inverted height), the amplitude A dif f , and the phase φ dif f . This time, we do not observe any instabilities, as the value of the transfer function of the cantilever decreases with the shifting of the resonance frequencies towards higher frequencies. This is also the reason, why we do not see an instability in the results of the experimental excitation scheme, which is presented in the main text. Let us, in the following, have a closer look to the implications on the contrasts, if applying a the specific excitation scheme.
We start with the height contrast. In both the off-off resonance and the off-on resonance excitation scheme, we observe that a softer sample (2 GPa) leads to a deeper indentation at a given contact force. Since we consider measurements that are performed with the feedback operating in contact mode, the contact force is held constant and a variation in elasticity results in different indentations, which translates into a measurable height signal: a harder material appears to be higher. This consideration holds for all excitation schemes including also the experimental excitation scheme, as discussed in the main text.
Considering the contrast in the amplitude A dif f that results from parts of the sample with different elasticities, we observe opposite behavior between the off-off resonance scheme and the off-on resonance scheme. In the off-off FIG. 12: The top panel shows the vibration spectrum in the off-on resonance scheme: blue lines indicate the excitation frequencies, whereas red lines indicate the resonance frequencies of the cantilever. The bottom three panels show the numerical results of the off-on resonance scheme for different sample elasticities: 2 GPa (black), 3 GPa (red), 4 GPa (magenta), 5 GPa (green), and 6 GPa (blue). The panels display from left to right the indentation (∼ inverted height), the amplitude A dif f , and the phase φ dif f . resonance excitation scheme, we see that a hard surface leads to a higher amplitude A dif f than a soft surface. This additionally supports both the experimental results of Sect. 4 and the analytical result of Sect. 5. In contrast, in the off-on resonance excitation scheme, we observe that for large contact forces (> 110 nN, see Fig. 12), a soft surface generates a higher amplitude A dif f than a hard one. We trust this result of our simulation at large contact forces, as the cantilever is completely vibrating in the Hertzian contact regime of the tip-sample interaction at these contact forces: the cantilever does not feel any attractive forces during is motion. This is not the case at lower forces (< 110 nN), where the contrast is inverted. Further evidence for a contrast inversion as a function of the applied contact force comes from the fact that the retract curves (not shown here) show exactly the same characteristics. The contrast inversion between the off-on resonance case and the off-off resonance case is caused by the frequency shift of the 4 th resonance frequency of the cantilever, which is explicitly excited in the off-on resonance excitation scheme.

Resonance Frequency of the "Nanoparticle in Polymer" system
For the estimation of the resonance frequency of the system gold nanoparticle in polymer, we need the spring constant of the PVA as well as the PMMA layer that are above and below the nanoparticle, respectively. From the stress equations of the polymer layers (see Eq. 10), we find the following spring constants: for which we used the same physical values as in Supplementary Note 6: R = 10 nm, d = 97 nm, t = 209 nm, E P V A = E P M M A = 2.4 GPa. Next we need the total mass M of the spherical nanoparticle. As the mass density, ρ Au , of gold is 19300 kg/m 3 , we find: M = 4 3 πR 3 ρ Au = 8.1 · 10 −20 kg (21) By assuming a simple harmonic oscillator, in which two springs are attached to a mass, we find an estimation for the resonance frequency: ω 0 = k P V A + k P M M A M = 1.4 · 10 10 rad/s.
This results in a resonance frequency of 2.2 GHz.