## 2 Principle and Method  
The numerical simulation of the transient bubble radial response is carried out by using the Yang-Church model [18], which considers the influence of medium viscosity on the vibration of a single microbubble. The vibration of a spherical bubble in an infinite range of incompressible thermoviscous medium is solved by the Keller-Miksis dynamic equation [19; 20]. \[(1-\frac{\dot{R}}{c})R\ddot{R}+\frac{3}{2}(1-\frac{\dot{R}}{3c})\dot{R}^{2}=(1+ \frac{\dot{R}}{c})\frac{P}{\rho}+\frac{R}{\rho c}\frac{dP}{dt} \tag{1}\]  
where \(R\) is the bubble radius, \(\dot{R}\) is the bubble wall velocity, \(\ddot{R}\) is the bubble wall acceleration, \(c\) is the sound velocity in the surrounding medium, \(\rho\) is the density of the surrounding medium. The bubble pressure \(P\) can be expressed as follows [21]:  
\[P=p_{g}-p_{0}-\frac{2\sigma}{R}+3\int_{R}^{\infty}\frac{\tau_{rr}}{r}dr+p_{a} \tag{2}\]  
where \(p_{g}\) is the gas pressure in the bubble, \(p_{0}\) is the ambient pressure, \(p_{a}\) is the excitation signal sound pressure, \(\sigma\) is the surface tension, \(\tau_{rr}\) is the shear stress, \(r\) is the radial radius of the bubble. The gas pressure \(p_{g}\) in the bubble can be expressed as follows [22]:  
\[p_{g}=p_{g0}(R_{0}/R)^{3c} \tag{3}\]  
where \(R_{0}\) is the initial bubble radius, \(\kappa\) is the polytropic index, which is determined by the pressure and temperature in the bubble, reflecting the relationship between the air density and the pressure change in the bubble.\(p_{g0}\) is the initial gas pressure in the bubble, defined as:  
\[p_{g0}=p_{0}+2\sigma/R_{0} \tag{4}\]  
The shear stress can be expressed as [23]:  
\[\tau_{rr}=2(G\gamma_{rr}+\mu\dot{\gamma}_{rr}) \tag{5}\]  
where \(\gamma_{rr}\) is the radial strain, \(\dot{\gamma}_{rr}\) is the strain rate, \(\mu\) is the viscosity, \(u\) is the velocity, \(u=(R^{2}/r^{2})\dot{R}\) and \(\dot{\gamma}_{rr}=\partial u/\partial rG\) is the shear modulus. In the near field near the bubble surface, \(\gamma_{rr}=(2/3r^{3})(R_{0}^{3}-R^{3})\),\(\dot{\gamma}_{rr}=-(2R^{2}/r^{3})\dot{R}\), the shear stress integral is:  
\[3\int_{R}^{\infty}\frac{\tau_{rr}}{r}dr=\frac{4G}{3R^{3}}(R_{0}^{3}-R^{3})- \frac{4\mu\dot{R}}{R} \tag{6}\]  
In summary, \(P\) can be expanded as follows:  
\[P=p_{g}-p_{0}-\frac{2\sigma}{R}+p_{a}-[\frac{4G}{3R^{3}}(R^{3}-R_{0}^{3})+ \frac{4\mu\dot{R}}{R}] \tag{7}\]  
and  
\[ \begin{array}{l}\frac{dP}{dt}=\frac{d}{dt}(p_{g}-p_{0}-\frac{2\sigma}{R}+p_ {a})+\frac{d}{dt}(3\int_{R}^{\infty}\frac{\tau_{rr}}{r}dr)\\ =\frac{dp_{g}}{dt}+\frac{2\sigma\dot{R}}{R^{2}}+\frac{dp_{a}}{dt}-4G\frac{R^{3 }\dot{R}}{R^{4}}-4\mu(-\frac{\dot{R}^{2}}{R^{2}}+\frac{\ddot{R}}{R})\end{array}  \tag{8}\]  
Combining Equation (3) and Equation (7), we can get:  
\[\frac{dP}{dt}=\left(\frac{2\sigma}{R}-3\kappa p_{g}\right)\frac{\dot{R}}{R}+ \frac{dp_{a}}{dt}-4G\frac{R_{0}^{3}\dot{R}}{R^{4}}+4\mu(\frac{\dot{R}^{2}}{R^ {2}}-\frac{\ddot{R}}{R}) \tag{9}\]  
For multifrequency ultrasound irradiation, the excitation signal can be expressed as follows:  
\[p_{a}=p_{1}\sin(2\pi f_{1}t)+\phi p_{2}\sin(2\pi f_{2}t+\alpha)+\phi p_{3}\sin( 2\pi f_{3}t+\beta) \tag{10}\]where \(f_{1}\), \(f_{2}\), and \(f_{3}\) are the frequency of excitation signal, \(p_{1}\), \(p_{2}\), and \(p_{3}\) are the amplitude of sound pressure, \(0\leq\varphi\leq 1\), \(0\leq\phi\leq 1\). \(\alpha\) and \(\beta\) are the phase angle of the excitation signal, \(0\leq\alpha\leq 2\pi\), \(0\leq\beta\leq 2\pi\). Equations (7), (9) and (10) are substituted into Equation (1) to extract the equation with \(\ddot{R}\), and then the ordinary differential equation solver (ode45) [24] in MATLAB (R2018b) software is used to solve it to obtain the bubble radius \(R\). In this paper, the initial radius of bubble \(R_{0}\) varies from 0.1 \(\upmu\)m to 10 \(\upmu\)m, and the initial interface velocity between the bubble and the surrounding medium is zero (\(\dot{R}=0\)).