## 2 Material and Methods
### Nonlinear Resonant Ultrasonic Spectroscopy Acquisition  
Nonlinearity was determined by analyzing the wave propagation velocity dependence on strain. To this purpose, Nonlinear Resonant Ultrasonic Spectroscopy (NRUS) was used. The method consists of measuring the resonance spectra and varying the amplitude of excitation and measuring the resonance frequency and resonance amplitude for each excitation.  
To determine the resonance spectra, a tapered chirp excitation was used to sweep over one longitudinal resonance peak:  
\[u(t)=A_{\text{drive}}\cos\left(2\pi\left(f_{\text{min}}+\frac{f_{\text{max}}- f_{\text{min}}}{\Delta t}t\right)t\right). \tag{1}\]  
The sweep was performed round a resonance mode (\(f_{\text{min}}\leq f_{\text{res}}\leq f_{\text{max}}\), where \(f_{\text{res}}\) is the resonance frequency). The bandwidth of the chirp spectra (\(f_{\text{max}}-f_{\text{min}}\)) was kept constant at 2.5 kHz, as was the duration \(\Delta t=50\) ms. This choice allowed for accurate resolution of the resonance peak. A discrete Fourier transform with 2 Hz spacing was computed for each response signal, followed by a deconvolution of the chirp spectrum, resulting in the elimination of spectral ripples and normalization of spectral amplitudes. This normalization ensured alignment between the spectral amplitude and the amplitude of the observed signal.  
The NRUS sequence consisted of selecting the drive amplitude \(A_{\text{drive}}\) (following the protocol in Figure 1c, i.e., alternating baseline and conditioning amplitudes). For each amplitude, the experimentally obtained spectrum (see Figure 2a) was fit with the MoDaNE solution [41]:  
\[B(f)=\frac{U_{0}}{\sqrt{\cosh^{2}[\alpha L]-\cos^{2}\left[2\pi\left(f+\psi \left(f-\frac{\eta v}{2L}\right)^{2}\right)L/v\right]}}, \tag{2}\]  
where \(f\) is frequency, \(n\) is the mode number, \(B\) is the amplitude and \(L\) is the sample length. The fit parameters \(U_{0}\), \(\alpha\), and \(v\) correspond to the uncalibrated source amplitude, damping coefficient, and wave velocity, respectively. The term proportional to \(\psi\) is a tilting term, which allows one to account for the sensitivity curve of transducers; \(\psi\) is thus not a material property. The resonance frequency of the sample is function of one of the fit parameters: \(f_{\text{res}}=nv/2L\).  
Repeating the procedure for different drive amplitudes, velocities, and damping was performed as a function of output amplitude, which is defined as the maximum of \(B\), i.e., the amplitude at the resonance frequency. Thus, it is derived as follows:  
\[A_{\text{out}}=\text{max}(B)=\frac{U_{0}}{\sqrt{\cosh^{2}(\alpha L)-1}}. \tag{3}\]  
The output amplitudes, from now on simply called amplitudes, are reported in volts, but, at least in the linear approximation, they are proportional to the strain amplitude in the center of the sample. Since our experiment is not calibrated, it was not possible to convert the voltage to strain, so whenever the term strain is used, it is intended as voltage.  
The amplitude \(A_{\text{out}}\) was derived for all drive amplitudes, both conditioning and baseline. Thus, it provided the definition of the conditioning strain and the baseline strain in the material. In the following, velocity variations will be plotted as a function of the output conditioning strain.