Estimation of Vibration Amplitude and Elastic Properties of Extra-capillary Tissue with Ultrasound Driven Vibration of Intra-capillary Gas Bubbles

Abstract

Estimation of vibration amplitude of intra-capillary micro-bubbles driven to vibrate with an incident ultrasound wave with amplitude and frequency to adjust the drive amplitude of the incident wave to obtain specified vibration amplitude of extra-capillary tissue. Estimation uses transmission of M groups of pulse complexes having low frequency pulse (LF) at bubble drive frequency, and high frequency (HF) pulse with angular frequency ω.sub.H> ~ 5 ω.sub.L, and pulse duration shorter than π/4ω.sub.L along HF beam. The phase between HF and LF pulses is ω.sub.Lt.sub.m for each group, where t.sub.m varies between the groups. Within each group, LF pulse varies between pulse complexes in amplitude and/or, where the LF pulse can be zero for a pulse complex, and LF pulse is different from zero for pulse complex within each group. HF receive signals are processed to obtain a parameter relating to bubble vibration amplitude when the HF pulse hits bubble.

Claims

1. A method for estimating parameters related to the surface or volume displacement vibration amplitude of at least one micro-bubble brought to vibration by an ultrasound wave of amplitude P.sub.L and angular frequency ω.sub.L transmitted along a LF beam with a given propagation direction along the unit vector e.sub.L, comprising a) transmitting towards said at least one micro-bubble at least two high frequency (HF) pulses with ω.sub.H > 5ω.sub.L along one HF beam along the unit vector e.sub.H at different times, b) directing at least one HF receive beam towards the region of the at least one micro-bubble with a direction e.sub.HR, the LF transmit and the HF transmit and receive beams are arranged so that they overlap at least in an observation region that includes said at least one micro-bubble, c) arranging the transmit timing of each of the transmitted HF pulses so that the pulses hits said at least one micro-bubble at different phases the LF oscillation, d) receiving at least two HF receive signals scattered from each said at least two transmitted HF pulses, and e) processing the at least two HF receive signals, to for each said at least one micro-bubble to produce an estimate of a bubble vibration parameter that gives one or both of the relative vibration of one or both of i) the bubble volume vibration relative to the resting bubble volume, and ii) a parameter related to relative bubble volume vibration.

2. The method according to claim 1, where the HF transmit beam is directed in close to the same direction as the LF transmit beam and the LF transmit and HF transmit and HF receive beams are positioned so that the three beams overlap in a region around said at least one micro-bubble, to produce estimates of said bubble vibration parameter for each said at least one micro-bubble.

3. The method according to claim 2, where a) a set S.sub.N of N≥1 groups G.sub.n of combined LF and HF pulse complexes are transmitted, each group G.sub.n comprising at least two pulse complexes, where for each group the LF pulse amplitude can be zero for one pulse complex, and the LF pulse amplitude is non-zero for at least one pulse complex, and at least one of the polarity, amplitude, or phase of the LF pulse varies between pulse complexes, and where b) for each group G.sub.n the HF pulse is transmitted with a timing relative t.sub.n to the LF pulse so that the HF pulses hit said at least one micro-bubble with a phase relation of ω.sub.Lt.sub.n relative to the transmitted LF oscillation, where the time t.sub.n is different for different groups G.sub.n, and c) the processing includes combining the HF receive signals from the scattered HF pulses from said at least two HF transmit pulses for each group G.sub.n to produce a detection parameter D.sub.n for group G.sub.n and each said at least one micro-bubble, and d) for each said at least one micro-bubble, the detection parameter D.sub.m for the group G.sub.m that gives maximal detection parameter for all G.sub.n ∈ S.sub.N is used as an estimate to obtain an estimate of the amplitude of said bubble vibration parameter.

4. The method according to claim 2, where a) a long LF pulse of duration T.sub.L is transmitted towards said at least one micro-bubble, where T.sub.L is longer than T.sub.R, preferably longer than N*TR, where T.sub.R is the propagation time from HF transmit transducer to said at least one micro-bubble, and b) a set S.sub.N of N ≥ 1 groups G.sub.n, comprising two HF pulses are transmitted at a reference time t.sub.n for each group, and internal timing within each group is selected so that at least two pulse complexes hits the bubble at opposite polarities of the LF oscillation, and c) the processing includes combining the scattered HF pulses from said at least two HF pulses for each group G.sub.n to produce a detection parameter D.sub.n for group G.sub.n and each said at least one micro-bubble, and d) for each said at least one micro-bubble, the detection parameter D.sub.m for the group G.sub.m, that gives maximal detection parameter for all G.sub.n ∈ S.sub.N is used as an estimate to obtain an estimate of the amplitude of said bubble vibration parameter.

5. The method according to claim 2, where the HF receive signals are compensated by one or both of i) a nonlinear propagation delay, and ii) a nonlinear pulse form distortion in the processing chain to produce estimates of said bubble vibration parameter for each said at least one micro-bubble.

6. The method according to claim 1, for an adequately dense group of micro bubbles, where a) the HF transmit beam defined by e.sub.H, is directed at an angle to the direction of the LF transmit beam defined by e.sub.L, and b) at least one HF receive beam are positioned together with the LF and HF transmit beams so that the three beams overlap at the region of said dense group of micro-bubbles, and c) transmitting a single group G.sub.n comprising at least two HF pulses at a reference time t.sub.n, and one of ci) the LF pulse polarity is switched between at least two of the HF pulses within the group, and cii) internal timing of the pulses within the group is selected, so that at least two HF pulse complexes hits at least one bubble within said group of micro-bubbles at opposite polarities of the LF oscillation, and d) processing that includes combining the received HF signal from said at least two HF pulses to produce a detection parameter D.sub.n for each bubble n in said group of micro-bubbles, and e) selecting estimate of a bubble vibration parameter the detection parameter D.sub.m for the bubble n = m that has the highest value.

7. The method according to claim 1, where said parameter related to relative bubble volume vibration is the average surface normal vibration amplitude relative to a radius of the bubble.

8. The method according to claim 7, where the value of the average surface normal is used to set the transmit amplitude of the LF wave.

9. The method according to claim 1, where the LF beam is focused onto in the therapeutic region.

10. The method according to claim 1, where the HF beam is broad with width less than the LF beam in the observation region, and several multi-focused parallel HF receive beams are generated to provide spatial resolution in the estimates of said bubble vibration parameter for a group of micro-bubbles.

11. The method according to claim 1, where scanning in 2 or 3 dimensions of one or more of i) the HF receive beam direction, and ii) the HF transmit beam direction, and iii) the LF beam direction to produce images of vibrating micro bubbles and the vibration amplitudes of the bubbles for a given transmitted amplitude of the LF wave.

12. The method according to claim 3, where the phase of the 2.sup.nd order transfer function between incident LF pressure oscillation and the amplitude bubble volume vibration is estimated from the timing t.sub.m that gives the maximal detection parameter D.sub.m.

13. The method according to claim 12, where the maximal detection parameter D.sub.m is obtained for at least two different frequencies of the LF beam, and estimates of the following parameters of the 2.sup.nd order approximation to the transfer function from incident pressure amplitude to bubble volume vibration i) the angular resonance frequency ω̂.sub.0 of the bubble ii) the absorption coefficient ζ̂̂ of the 2.sup.nd order approximation to the incident LF pressure to bubble volume transfer function, iii) the gain K.sub.v of the 2.sup.nd order approximation to the incident LF pressure to bubble volume transfer function.

14. The method according to claim 13, where for a bubble filling parts of a vessel estimates of at least one of the following vibration parameters for the surrounding tissue is obtained i) the ratio of the tissue shear wave velocity c.sub.s to the vessel radius a.sub.v is obtained from the estimate of the angular resonance frequency ω̂.sub.0 of the bubble, ii) for a known radius of the vessel the real part of tissue shear stiffness is obtained from the estimate of the angular resonance frequency ω̂.sub.0 of the bubble, iii) for a known radius of the vessel the imaginary part of tissue shear stiffness is obtained from the estimate of the angular resonance frequency ω̂.sub.0 of the bubble.

15. An instrument for obtaining estimates of parameters related to the surface or volume displacement vibration amplitude of at least one micro-bubble brought to vibration by an ultrasound wave of amplitude P.sub.L and angular frequency ω.sub.L, comprising a) means for transmitting towards said at least one micro-bubble ai) at least one LF pulse along a LF transmit beam with a transmit beam direction along a unit vector e.sub.L, and aii) at least two high frequency (HF) pulses with ω.sub.H> 5ω.sub.L along a HF beam along the unit vector e.sub.H at different times, and iii) means for arranging the transmit timing of each of the transmitted HF pulses so that the pulses hit said at least one micro-bubble at different phases the LF oscillation, and b) means for receiving that includes i) means for directing at least one HF receive beam with direction e.sub.HR to produce an overlap region (OR) with the LF and HF transmit beams that includes the at least one micro-bubble, and ii) means for receiving at least two HF receive signals scattered from each said at least one micro-bubble from said at least two transmitted HF pulses, and c) means for processing the at least two HF receive signals, to for each said at least one micro-bubble to produce an estimate of a bubble vibration parameter that gives one or both of the relative vibration of one or both of i) the bubble volume vibration relative to the resting bubble volume, and ii) a parameter related to relative bubble volume vibration.

16. The instrument according to claim 15, where a) said means for transmitting include means for transmitting said HF transmit beam in close to the same direction as the LF transmit beam, and b) said means for receiving include means for directing the HF receive beam so that the LF transmit and HF transmit and HF receive beams overlap at the position of said at least one micro-bubble.

17. The instrument according to claim 16, where a) said means for transmitting includes i) means for transmitting a set S.sub.N of N≥1 groups G.sub.n of combined LF and HF pulse complexes, each group G.sub.n comprising at least two pulse complexes, where for each group the LF pulse amplitude can be zero for one pulse complex, and the LF pulse amplitude is non-zero for at least one pulse complex, and ii) means for varying at least one of the polarity, amplitude, or phase of the LF pulse between pulse complexes, and iii) means that for group G.sub.n for sets the transmit time t.sub.n of the HF pulse relative to the LF pulse so that the HF pulses hit said at least one micro-bubble with a phase relation of ω.sub.Lt.sub.n relative to the LF pulse oscillation when the HF pulse hits the at least one micro-bubble, where the time t.sub.n is different for different groups G.sub.n. and b) said means for processing includes i) means for combining the scattered HF pulses from said at least two HF pulses for each group G.sub.n to produce a detection parameter D.sub.n for group G.sub.n and each said at least one micro-bubble, and ii) means for selecting for each said at least one micro-bubble, the detection parameter D.sub.m for the group G.sub.m, that gives maximal detection parameter for all G.sub.n ∈ S.sub.N, said D.sub.m is used as a basis to obtain an estimate of the amplitude of said bubble vibration parameter.

18. The instrument according to claim 16, where a) said means for transmitting include ai) means for transmitting a long LF pulse of duration T.sub.L towards said at least one micro-bubble, where T.sub.L is long, preferably longer than NT.sub.R where T.sub.R is the propagation time from HF transmit transducer to said at least one micro-bubble, and ii) a set S.sub.N of N ≥ 1 groups G.sub.n comprising two HF pulses transmitted at a reference time t.sub.n for each group, and selecting internal timing within each group so that at least two pulse complexes hits the bubble at opposite polarities of the LF oscillation, and b) said means for processing includes bi) means for combining the received scattered HF pulses from said at least two HF transmit pulses for each group G.sub.n to produce a detection parameter D.sub.n for group G.sub.n and each said at least one micro-bubble, and bii) means for selecting for each said at least one micro-bubble the detection parameter D.sub.m for group G.sub.m, that gives maximal detection parameter for all G.sub.n ∈ S.sub.N, Said G.sub.m is used to obtain an estimate of the amplitude of said bubble vibration parameter for each said at least one micro-bubble.

19. The instrument according to claim 16, where said processing means include means for compensating the HF receive signals by one or both of i) a nonlinear propagation delay, and ii) a nonlinear pulse form distortion in the processing chain to produce estimates of said bubble vibration parameter for each said at least one micro-bubble.

20. The instrument according to claim 15, for an adequately dense group of micro bubbles, where a) said means for transmitting includes ai) means for directing the HF transmit beam defined by e.sub.H, at an angle to the direction of the LF transmit beam defined by e.sub.L, and aii) means for transmitting a single group G.sub.n comprising at least two HF pulses at a reference time t.sub.n, and one of a) the LF pulse polarity is switched between at least two of the HF pulses within the group, and ii) internal timing of the pulses within the group is selected, so that at least two HF pulse complexes hits at least one bubble within said group of micro-bubbles at opposite polarities of the LF oscillation, and b) said means for receiving includes means for directing at least one HF receive beam together with the LF and HF transmit beams so that the three beams overlap at a region of said dense group of micro-bubbles, and c) said means for processing includes ci) means combining the received HF signal from said at least two HF pulses to produce a detection parameter D.sub.n for each bubble n in said group of micro-bubbles, and cii) means for selecting the detection parameter D.sub.m for the bubble n = m that has the highest value, and use D.sub.m in the processing chain to produce and estimate of a bubble vibration parameter.

21. The instrument according to claim 15, where the LF beam is focused onto in the therapeutic region.

22. The instrument according to claim 15, where the HF transmit beam is broad with width less than the LF beam in the observation region, and several multi-focused parallel HF receive beams are generated to provide spatial resolution in the estimates of said bubble vibration parameter for a group of micro-bubbles.

23. The instrument according to claim 15, where said parameter related to relative bubble volume vibration is the average surface normal vibration amplitude relative to a radius of the bubble.

24. The instrument according to claim 23, comprising means to use the value of the average surface normal to set the transmit amplitude of the LF wave.

25. The instrument according to claim 15, where said transmit and receive means comprises means for scanning in 2 or 3 dimensions of one or more of i) the HF receive beam direction, and ii) the HF transmit beam direction, and iii) the LF beam direction to produce images of vibrating micro bubbles and the vibration amplitudes of the bubbles for a given transmitted amplitude of the LF wave.

26. The instrument according to claim 17, where said processing means include means to estimate the phase of the 2.sup.nd order transfer function between incident LF pressure oscillation and the amplitude bubble volume vibration from the timing t.sub.m that gives the maximal detection parameter D.sub.m.

27. The instrument according to claim 26, where the processing means includes a) means for detection of the maximal detection parameter D.sub.m for at least two different frequencies of the LF beam, and b) means for estimating the following parameters of the 2.sup.nd order approximation to the transfer function from incident pressure amplitude to bubble volume vibration i) the angular resonance frequency ω̂.sub.0 for the bubble ii) the absorption coefficient ζ of the 2.sup.nd order approximation to the incident LF pressure to bubble volume transfer function, iii) the gain K.sub.v of the 2.sup.nd order approximation to the incident LF pressure to bubble volume transfer function.

28. The instrument according to claim 27, where the processing means comprises means to estimate at least one of the following vibration parameters for the tissue surrounding a bubble that is filling a part of a vessel i) the ratio of the tissue shear wave velocity c.sub.s to the vessel radius a.sub.v from estimates of the angular resonance frequency ω̂.sub.0 of the bubble, and ii) the real part of tissue shear stiffness for a known radius of the vessel from estimates of the angular resonance frequency ω̂.sub.0 of the bubble, and iii) the imaginary part of tissue shear stiffness for a known radius of the vessel from the estimate of the angular resonance frequency ω̂.sub.0 of the bubble .

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0015] FIG. 1a Illustrates a spherical micro-bubble with radius a in a larger vessel.

[0016] FIG. 1b Illustrates a micro-bubble filling a part of a capillary used to vibrate the capillary wall and extra capillary tissue

[0017] FIG. 2a Illustrates interface vibration patterns of a capillary surface vibrated by a micro-bubble as shown in FIG. 1b.

[0018] FIG. 2b Illustrates a numerical simulation of the displacement vibration amplitude as a function of position in the tissue surrounding the vibrating micro-bubble that fill s part of the capillary.

[0019] FIG. 3 Illustrates approximation of a 2nd order transfer function Eq.(5, A14) to the transfer function in Eq.(A12) from incident wave amplitude P.sub.i at a frequency f to the bubble volume complex vibration amplitude.

[0020] FIG. 4 illustrates two typical combined low frequency (LF) and high frequency (HF) ultrasound pulse complexes used to estimate the LF driven vibration amplitude and transfer function parameters of the micro-bubble.

[0021] FIGS. 5a - c Illustrates an ultrasound dual frequency array with LF transmit and HF transmit and receive beams typically used for identification of bubble vibration amplitude and tissue elastic properties.

[0022] FIG. 6a Illustrates the variation of H.sub.b(K) for a spherical bubble as a function of frequency when the bubble diameter 2a becomes comparable to the incident wave length.

[0023] FIGS. 6b - d Illustrates the variation of the H.sub.b(K) for a bubble that fills part of a capillary of radius a.sub.v= 5 .Math.m as a function of angular direction and incident frequency for a cylindrical length 2/.sub.c.

[0024] FIGS. 7 Shows the detection sensitivity function of the average relative radial vibration amplitude for a bubble that fills part of a capillary with relative cylindrical half length l.sub.c/a.sub.v.

[0025] FIG. 7b illustrates the logarithmic amplitude and phase of a 2.sup.nd order volume vibration transfer function in Eqs.(5, A14), with illustration of a linear approximation to parts of the phase function.

[0026] FIG. 8a Illustrates the variation of the relative resonance x.sub.d of the relative volume transfer function F.sub.d in Eq.(A11) as a function of relative cylindrical half length l.sub.c/a.sub.v and 3 values of the tissue absorption parameter.

[0027] FIG. 8b Illustrates estimated ratio of the tissue shear wave velocity to the capillary radius, c.sub.s/a.sub.v, Eq.(32), as a function of observed bubble resonance frequency f.sub.d.

[0028] FIG. 8c Illustrates a logarithmic plot of the real component of the tissue shear stiffness as a function of observed bubble resonance frequency f.sub.d and 4 capillary radii.

[0029] FIG. 9 Illustrates an overview block diagram of a system for carrying through the methods according to the invention.

DETAILED DESCRIPTION OF EMBODIMENTS

[0030] This section gives a more detailed description of example embodiments of the invention, where the examples by no means represents limitations of the invention, which in its broadest aspect is defined by the claims appended hereto.

A. Fundamental Bubble Vibration

[0031] The invention presents methods and instrumentation for extracting information about ultrasound driven vibrations of micro-bubbles in tissue, and also to use this information for estimation of tissue properties, such as elastic properties, and properties relating to elastic properties.

[0032] The invention is relevant for micro-bubbles with diameter less than the capillaries, and also bubbles that in large arteries have a spherical shape with radius a ~ 5 - 30 .Math.m that gives a bubble volume V.sub.b = 4πα.sup.3/3, while in capillaries with radii a.sub.v ~ 3 - 5 .Math.m the larger bubbles fill at least part of the capillary. Examples of the larger bubbles are shown in FIG. 1, where 101 shows a spherical bubble in a larger artery 100, while 111 shows the a bubble filling a smaller capillary 110 with radius a.sub.v for a cylindrical length 2/.sub.c with two semi-spherical endings 112 of radius a.sub.v.

[0033] There is spatially distributed elasticity and mass density of the tissue, the blood, and the gas, that implies wave propagation with reflection at boundaries. The tissue and the blood have both a high volume compression stiffness λ.sub.t = 1/κ.sub.t(κ.sub.t- volume compressibility) that upholds pressure waves with propagation velocity c.sub.p ≈1500 m/s for both tissue and blood. The tissue has in addition a low shear/deformation stiffness .Math..sub.t that produces tissue shear/deformation waves with no volume compression and low propagation velocity c.sub.s ~ 2 - 20 m/s. The blood and the gas have negligible shear stiffness, and hence no shear waves. The bubble gas has a low volume compression stiffness, but the mass density is also so low that the gas pressure can be approximated as constant across the bubble volume. The bubble gas pressure depends on the bubble volume and its variation as

[00001]$\begin{array}{l}{P}_g=P_0{\left(\frac{V_b}{V_b+\delta V}\right)}^{\gamma }\approx P_0-\gamma P_0\frac{\delta V}{V_b}\\ \gamma =\frac{C_p}{C_V}\approx 1\\ \gamma =C_p/C_v\approx 1\end{array}$

where P.sub.0 ~100 kPa is the environmental pressure, V.sub.b is the equilibrium bubble volume, δV is the change in bubble volume, C.sub.p and C.sub.vare the heat capacities at constant pressure and volume. Typical gases have large molecules that makes C.sub.p ≈ Cv and γ≈ 1. The volume-pressure relation is inherently nonlinear for larger changes in the volume, but we shall in our work use the approximate linear variation with δV/V.sub.b, to study fundamental aspects of the bubble vibration.

[0034] The low tissue shear stiffness allows large deformations of the cylindrical bubble surface 111, with large amplitude, short-range shear deformation of the surrounding tissue with comparatively low volume compression. The large vibration amplitude of the close surrounding tissue produces an outward acoustic radiation force (ARF) that has interesting applications for drug transport both across the capillary wall, through the interstitial space between the tissue cells, and across membranes of the cells and other tissue structures. In this context it is important to obtain adequate bubble/tissue vibration amplitude (~ 1 .Math.m) and one target of this invention is to use ultrasound scattering from vibrating bubbles to measure the bubble vibration amplitude produced by an incident ultrasound wave. The invention further uses the measured bubble vibration amplitude to set the amplitude of the incident vibration wave to achieve desired bubble vibration amplitude in the surrounding tissue to achieve the desired therapeutic effect.

[0035] A bubble in a large vessel, as 101 in 100, has a practically spherical shape with a spherical vibration pattern. A drive pressure P.sub.d that is the difference between the gas pressure P.sub.g and the pressure P.sub.∞ at a point so distant from the bubble that one can neglect the tissue shear vibration, produces a complex volume vibration amplitude □V.sub.c of the bubble as [1]

[00002]$\begin{array}{c}\begin{array}{ll}\frac{\delta V_c}{V_b}=H_d\left(\omega \right)P_d\hfill & H_d\left(\omega \right)\approx \frac{-K_d}{{\omega }^2-{\omega }_d^2-i2\zeta {\omega }_d\omega }\hfill \end{array}\\ \begin{array}{lll}{K}_d=\frac{3}{\rho a^2}\hfill & {\omega }_d^2=\frac{4E_{st}}{\rho a^2}\hfill & 2\zeta {\omega }_d=\frac{4{\eta }_b}{\rho a^2}\hfill \end{array}\end{array}$

where ρ.sub.b and η.sub.b are the mass density and coefficient of viscosity of the surrounding blood, and E.sub.st is the Youngs module of the blood/gas surface layer. We refer to H.sub.d as the direct drive transfer function, and ω.sub.d as the direct drive angular resonance frequency produced by the interaction between the co-oscillating mass of the surrounding blood and the elasticity of the blood/gas surface layer. Viscous friction in the blood and the surface layer is given by η.sub.b.

[0036] An incident acoustic wave with pressure amplitude P.sub.i and angular frequency ω produces a variation in the bubble volume and the gas pressure that produces drive pressure amplitude P.sub.d as

[00003]$P_d=-\gamma P_0\delta V_c/V_b-P_i=-\gamma P_0{H}_d\left(\omega \right)P_d-P_i$

[00004]$P_d=\frac{-P_i}{1+\gamma H_d\left(\omega \right)P_0}$

which gives the complete transfer function H.sub.V(ω) from the incident pressure amplitude P.sub.i to the complex relative volume amplitude as

[00005]$\frac{\delta V_c}{V_b}=-H_V\left(\omega \right)P_i{H}_V\left(\omega \right)=\frac{H_d\left(\omega \right)}{1+\gamma H_d\left(\omega \right)P_0}$

[0037] With the 2.sup.nd order form of the H.sub.d(ω) in Eq.(2), we get the following form of H.sub.V(ω)

[00006]$\begin{array}{l}{H}_V=\left(\omega \right)=\frac{-K_d}{{\omega }^2-{\omega }_0^2-i2\zeta {\omega }_0\omega }{\omega }_0^2={\omega }_d^2+\gamma K_d{P}_0\\ \zeta ={\zeta }_d{\omega }_d/{\omega }_0\end{array}$

[0038] For a bubble that partially fills a capillary with a large cylindrical region as 111 in FIG. 1b, the vibration pattern is more complex than the spherical vibrations behind the developments in Eqs.(2-5). The difference in material parameters across the vessel wall, i.e. between gas/tissue for the cylindrical bubble surface region (S.sub.c), and blood/tissue for the blood surface region (S.sub.v), produce interface waves along the vessel wall with propagation velocity c.sub.i ~ c.sub.s in the vascular region and c.sub.i ~ 2c.sub.s where the bubble contacts the vessel, as illustrated in FIG. 2a.

[0039] The FIG. 2a shows simulations of the spatial variation of the vibration amplitude of a bubble that fills part of a capillary (201 diam 2a.sub.v = 10 .Math.m) together with interface waves along the blood/tissue interface outside the bubble. The upper panel 200 shows the vibration amplitudes at a frequency close to the volume vibration resonance, and the lower panel 210 shows a frequency 6 times the volume resonance frequency. At equilibrium the cylindrical part of the bubble fills the capillary and ends with the semi-spheres 203. Close to the resonance frequency in 200 the bubble surface amplitude 202 is magnified 5 times [5X] while the blood tissue interface wave amplitude 204 is magnified 50 times [50X]. At 6 times the resonance frequency in 210 the bubble surface amplitude 212 is magnified 7.5 times [7.5X] while the blood tissue interface wave 204 is magnified 75 times [75X] interface wave.

[0040] FIG. 2b gives as 221 a grey-scale display of the simulated tissue displacement amplitude around the bubble (220) close to the bubble volume resonance frequency. The volume amplitude can be expressed by the average normal vibration displacement amplitude across the bubble surface S.sub.b with area A.sub.b as

[00007]$\delta V_c=A_b{\overline{\psi }}_n{\overline{\psi }}_n=\frac{1}{A_b}{\displaystyle \underset{S_b}{\int }{d}^2{r}_0{\psi }_n\left({\underset{\_}{r}}_0\right)}$

[0041] Eqs.(12,13,19) shows that for a bubble with dimensions adequately smaller than the incident wavelength the scattered signal is determined by the bubble volume V.sub.b and bubble volume compressibility κ.sub.b. We note from FIG. 2b that the normal displacement amplitude on the semi-spherical ends is lower than the cylinder radial displacement amplitude, and these two observations stimulates an approximation of the bubble with cylindrical length 2/.sub.c and semispherical ends of radius a.sub.v, with a cylinder of fixed length 2 L where the radius a.sub.v vibrates. The bubble volume V.sub.b for this cylindrical bubble is the same as for the original spherical bubble 101 in a large artery, and the bubble that fills part of a capillary in 111. We have the relation

[00008]$\begin{array}{l}{V}_b=\frac{4}{3}\pi a^3=\pi a_v^{2}2l_c+\frac{4}{3}\pi a_v^3=\pi a_v^{2}2LL=l_c+\frac{2}{3}{a}_v\\ \begin{array}{cc}\frac{l_c}{a_v}=\frac{2}{3}\left({\left(\frac{a}{a_v}\right)}^3-1\right)& \frac{L}{a_v}=\frac{l_c}{a_v}+\frac{2}{3}=\end{array}\frac{2}{3}{\left(\frac{a}{a_v}\right)}^3\end{array}$

where 2 L is the length of a cylindrical bubble in a capillary with the same volume V.sub.b. The vibrating area of the cylindrical bubble is

[00009]$A_c=4\pi a_{v}L=4\pi a_v{l}_c+\frac{8\pi }{3}{a}_v^2=A_b-\frac{4\pi }{3}{a}_v^2$

[0042] We notice that Ac is

[00010]$4\pi a_v^2/3$

which is ⅓ of the surface area of the end spheres. However, the effect of this error is reduced because the displacement of the semi-spherical end surfaces is less than the radial vibration of the cylinder surface as shown in FIG. 2b. We will therefore use the simplifying approximation

[00011]$\begin{array}{ccc}\begin{array}{cc}\delta V_c=2\pi a_v{\displaystyle \underset{-L}{\overset{L}{\int }}dz\cdot {\psi }_r\left({\underset{¯}{r}}_0\right)=A_c{\overline{\psi }}_r}& \end{array}& {\overline{\psi }}_r=\frac{1}{2L}{\displaystyle \underset{-L}{\overset{L}{\int }}dz\cdot {\psi }_r\left({\underset{¯}{r}}_0\right)}& A_c=4\pi a_{v}L\end{array}$

[0043] Analysis of vibration and scattering from the cylindrical bubble is carried through in Appendix A. We note from 210 in FIG. 2a that at the higher frequency range we get multiple interface wave lengths along the bubble cylindrical surface which complicates the approximation with the 2.sup.nd order transfer function in Eq.(5) over a larger frequency range. FIG. 3a shows as the solid lines a simulation of the magnitude |H.sub.v(ω)| (300) and the phase

[00012]${\theta }_V\left(\omega \right)=\angle {\text{H}}_V\left(\omega \right)\left(301\right)\text{of}{\text{H}}_V\left(\square \right)$

from Eq.(A10, A12) of Appendix A for 1.sub.c/a.sub.v=5. The dotted lines show the magnitude |Ĥ.sub.v(ω)| (302) and the phase θ̂v(ω) (303) of the 2.sup.nd order approximation |Ĥ.sub.v(ω)| of Eq.(A14). The approximation of the 2.sup.nd order function to Eq.(A12) is done in the following steps: [0044] 1. Re{H.sub.v(ω.sub.0)}=0 of Eq.(A14) gives the angular resonance frequency ω.sub.0. [0045] 2. ζ is determined from the differential of the ohase of H.sub.V as in Eq.(25). [0046] 3. ω.sub.m and K.sub.V is determined from Eq.(27)

[0047] In the example in FIG. 3a we see that f.sub.0 ≈ 340 kHz and we can say that the approximation is very acceptable in a range larger than Ω.sub.0 = 100 - 1000 kHz = (f.sub.0/3, 3f.sub.0).

[0048] To observe the LF vibration amplitude of a bubble driven by the incident LF wave with angular frequency ω.sub.L , the invention presents methods and instrumentation building on U.S. Pat. 8,550,998; 9,291,493; 9,939,413; 10,578,725 and 7,727,156; 10,879,867.

[0049] An LF pulse shown by example as 401 or 405 in FIG. 4 is typically transmitted along an LF transmit beam (LT - beam) shown by example as 501 of FIG. 5a and b, where 502 illustrates a therapy region. We have illustrated the LF transmit beam 501 with a large aperture and focused onto the therapy region 502 which allows higher LF pressure in the therapy region than in front of this region. For adequate therapy we want a certain surface displacement amplitude of the intra-capillary bubbles (FIG. 1b). For hard solid tumors, for example with large ingrowth of connective tissue and/or increased intra-tumoral pressure, adequate displacement amplitude requires increased LF pressure, e.g. ~ ~ 1 Mpa which gives an MI ~ 1.6 at LF ~ 400 kHz. This is seen as a risk MI for an intra-capillary contrast agent bubble with diam ~ 3 .Math.m which could get excessive displacement collapse. However for a bubble that fills the capillary, the vibration displacement is more controlled by the surrounding tissue, and for the corresponding spherical bubbles in the larger arteries and arterioles the artery wall is much stronger so that dangerous damage is less likely. In the normal tissue in front of the tumor where the extra-capillary tissue has much lower shear stiffness, one should have lower LF pressure amplitude. The focused LF beam provides such spatial variation of the LF amplitude.

[0050] In many situations, for example trans-abdominal and trans-costal applications, the access window to the therapeutic region is limited, and it is desirable to have a HF transmit and receive beam arrangement with the same direction as the LF beam as illustrated in FIG. 5a, illustrated as 503 for the HF transmit beam and 504 for the HF receive beam. One would then typically use a wide HF transmit beam (503) with well multi-depth focused parallel HF receive beams 504 illustrated by the dots 505. Such multi-depth focusing is today obtained by software receive beam-formers utilizing parallel processing of Graphic Processing Units (GPUs). Hardware receive beamformers with such parallel receive beams also exists to-day.

[0051] With adequate ultrasound access windows, it can be useful to use crossing multi/depth focused HF receive beams exemplified by 514 and the dots 515 in FIG. 5b, with the forward HF transmit beam as 503 in FIG. 5a. This reduces the length of the scattering wave vector K=k(e.sub.i-e.sub.s) as defined in Eq.(12). As discussed after Eq.(12) and in relation to FIGS. 6a - d, this length reduction opens for the us higher HF frequencies with controlled reduction in the magnitude of H.sub.b and ΔH.sub.b/H.sub.b value produced by LF variation of the scatterer dimension.

[0052] In a particular solution, one transmits a set S.sub.N of N≥ 1 groups G.sub.n pulse complexes where one for each group G.sub.n transmits at least two pulse complexes comprising a LF manipulation pulse and a HF observation pulse where the LF pulse varies in one of amplitude, and polarity, and phase between each pulse complexes, where the LF pulse might be zero for one pulse complex, and the LF pulse is non-zero for at least one pulse complex. The LF pulses 401 and 405 in FIG. 4 are example LF pulses for two pulse complexes with opposite polarity of the LF pulses for one group G.sub.n. The position of the HF pulse relative to the LF pulse is for the group G.sub.n defined by the timing t.sub.n of the HF pulse transmission. The timing varies between each group indicted by the time points ..., t.sub.n-1, t.sub.n, t.sub.n+1, ... shown as 404 in FIG. 4. The envelopes of the LF pulses varies slowly and shown as 404 in FIG. 4. The HF is typically > ~ 5*LF. The LF and the HF pulses can have different phase propagation with depth so that the positions t.sub.n of the HF pulses shown in FIG. 4 indicates the position of the HF pulses relative to the LF pulse at the depth of the bubble. This change in phase position can for example be obtained by numerical simulations of the combined LF and HF propagation to the bubble, or through measurements in a water or oil tank. For each group G.sub.n defined with position t.sub.n of the HF pulse, at least two pulse complexes are transmitted with different polarities and/or amplitudes of the LF pulse, where one LF pulse can have zero amplitude and at least one LF pulse has non-zero amplitude. In a preferred embodiment one transmits two pulse complexes with opposite polarities of the LF pulse shown as 401 and 405 in FIG. 4.

[0053] The time between transmit of the HF pulses is determined so that it is larger than the round trip propagation time T.sub.R from HF transmit transducer to the micro-bubble and to the HF receive transducer. In particular solutions where one wants long duration LF transmit vibration of the bubble we can use LF pulses with long duration time T.sub.L, for example longer than NT.sub.R where N is the the number of groups G.sub.n. Several groups G.sub.n of HF transmit pulses are then transmitted with time distance larger than T.sub.R and t.sub.n defining the detailed phase relation between the HF pulse and the LF oscillation of the bubble when the HF pulse hits the bubble.

[0054] It is also interesting the use a HF transmit beam 513 that crosses the LF beam 501 in the same direction as the set of HF receive beams 514/515. FIG. 5c illustrates regions of constant phase relation between the HF and LF pulses in the case of crossing HF and LF transmit beams as in FIG. 5b, where said regions are along lines normal to the sum of the unit direction vectors of the HF and LF transmit beams. The advantage with the crossing HF transmit beam is that for a fixed phase transmit timing t.sub.n between the HF pulse and the phase of the LF oscillation, there will be a spatial variation of said phase relation through the therapy region, and there will be line regions exemplified by 516 in FIG. 5c along which said phase relation is constant. Said line regions have direction normal to e.sub.L + e.sub.H, where e.sub.L is the unit vector along the LF transmit beam axis and e.sub.H is the unit vector along the HF transmit beam axis. With adequate density of intra-capillary bubbles in the therapy region 502 there will be intra-capillary bubbles with adequate closeness to these line regions to produce the full amplitude of the bubble volume oscillation and full average bubble surface displacement amplitude Ψ as defined in Eq.(9). Selecting the detection signal D.sub.cm according to Eqs.(18,20) from these bubbles with maximal D.sub.cn provides detection signals with the sought after phase relation between the HF pulses and the bubble vibrations when the HF pulses hits the bubble.

[0055] For therapy we generally drive the bubble with the low frequency (LF) pressure wave with angular frequency ω.sub.L, typically in the neighborhood of the bubble resonance frequency to produce large therapeutic cylindrical radius vibration amplitude (~ 1 .Math.m). The LF pulse has an amplitude envelope pP.sub.Le(t), where p is a polarity parameter of the LFwave. p = 1 implies positive LF amplitude (e.g. 401), p = -1 implies inverted LF amplitude (e.g. 405), while p = 0 implies zero LF amplitude. Other values/sequences of the polarity parameter can be used to enhance the signal processing, where the cited patents show examples of p = 0, ±1; p = ±1, ±2, for example to improve suppression of tissue noise and/or compensate for bubble movements between the HF pulses.

[0056] The LF pulse produces a vibration in bubble volume as

[00013]$\begin{array}{l}\delta {\text{V}}_b\left(t,{\omega }_L\right)=p\delta {\text{V}}_{Le}\left(t\right)\cos \left[{\omega }_{L}t+\theta \left({\omega }_L\right)\right]\\ \\ \delta {\text{V}}_{Le}\left(t\right)\approx \left|H\left({\omega }_L\right)\right|P_{Le}\left[t+{\tau }_e\left({\omega }_L\right)\right]\approx \left|H\left({\omega }_L\right)\right|P_{Le}\left(0\right)\end{array}$

[0057] Where θ(ω.sub.L)=∠H(ω.sub.L) is the phase of H(ω.sub.L), and τ.sub.e(ω.sub.L)≈θʹ(ω.sub.L) is an approximate delay of the envelope produced by the transfer function H(ω) for an adequately narrow band HF pulse. The last approximation of δV.sub.Le(t) assumes that he envelope varies slowly relative to the ω.sub.L oscillation. To measure the bubble vibration amplitude we transmit HF observation pulses as described in FIG. 4. The bubble volume vibration when the HF pulse hits the bubble is then

[00014]$p\delta V_n=\delta V_b\left(t_n,{\omega }_L\right)=p\delta V_{Le}\left(t\right)\cos \left[{\omega }_L{t}_n+\theta \left({\omega }_L\right)\right]$

[0058] The scattered pressure field of an incident ultrasound pressure wave with angular frequency ω and amplitude P.sub.i from a bubble with volume V.sub.b and adequately less extension than the incident wavelength λ, is by the Born approximation obtained as

[00015]$\begin{array}{l}{p}_s\left(\underset{\_}{r},\omega \right)\approx k^2\left(\frac{{\kappa }_b\left(\omega \right)}{{\kappa }_t\left(\omega \right)}-{\underset{\_}{e}}_s{\underset{\_}{e}}_i\right)V_b{H}_b\left(\underset{\_}{K}\right)\frac{e^{-ikr}}{4\pi r}{P}_i\\ {\kappa }_b\left(\omega \right)=-\frac{1}{P_i}\frac{\delta V_b}{V_b}{\kappa }_t\left(\omega \right)=-\frac{1}{P_i}\frac{\delta V_t}{V_t}\\ H_b\left(\underset{\_}{K}\right)=\frac{1}{V_b}{\displaystyle \underset{Vb}{\int }{d}^3{r}_0{e}^{-i\underset{\_}{K}{\underset{\_}{r}}_0}}\underset{\_}{K}={\underset{\_}{k}}_i-{\underset{\_}{k}}_s=k\left({\underset{\_}{e}}_i-{\underset{\_}{e}}_s\right)\\ k=\frac{2\pi }{\lambda }=\frac{\omega }{c_p}\end{array}$

where [0059] r is the radial distance from the bubble center of mass, [0060] k is the wave vector magnitude of both the incident and scattered wave, and λ = 2π/k is the wave length of both the incident and scattered waves. [0061] e.sub.i is the unit direction vector of the incident wave approximated as plane wave across the bubble with amplitude P.sub.i, [0062] e.sub.s is the unit direction vector of the observed scattered wave, [0063] c.sub.p is the pressure wave propagation velocity, [0064] δV.sub.b is the vibration amplitude of the bubble volume, V.sub.b, and δV.sub.t is the vibration amplitude of the surrounding tissue volume, V.sub.t, both produced by the incident wave, [0065] κ.sub.b (ω) is the bubble volume compressibility, and [0066] κ.sub.t (ω) is the volume compressibility of the surrounding material.

[0067] The term κ.sub.b/κ.sub.t represents scattering due to the deviation of the bubble volume compressibility from the tissue volume compressibility, and the term e.sub.se.sub.i represents scattering due the relative difference between the bubble mass density (≈ 0) and the tissue mass density ρ.sub.t. In a typical situation |κ.sub.b(ω)/κ.sub.t|>>|e.sub.se.sub.i|, and the term e.sub.se.sub.i can be neglected. Using

[00016]$c_p^2=1/{\rho }_t{\kappa }_t$

we modify Eq.(6) to

[00017]$\begin{array}{l}{p}_s\left(\underset{\_}{r},\omega \right)\approx {\rho }_t{\omega }^2{\kappa }_b\left(\omega \right)V_b{H}_b\left(\underset{\_}{K}\right)\frac{e^{-ikr}}{4\pi \text{}r}{P}_i\\ {\kappa }_b\left(\omega \right)=-\frac{1}{P_i}\frac{\delta V_{cb}}{V_b}=H_V\left(\omega \right)\end{array}$

[0068] Both for contrast agent bubbles and bubbles for therapy, we are generally interested in vibrating the bubbles in the low frequency (LF) range Ω.sub.0 around the bubble volume vibration resonance frequency, where the 2.sup.nd order transfer functions of Eqs.(5,A14) are useful. To measure the bubble vibration amplitude at ω.sub.L we transmit in addition high frequency (HF) observation pulses with angular center frequency ω.sub.H ~10ω.sub.L >> ω.sub.0 as illustrated in FIG. 4. We define this HF range as Ω.sub.H The 2.sup.nd order transfer function of Eq.(5) is still valid at ω.sub.H for the spherical bubble in FIG. 1a. For the bubbles that fills part of a capillary in FIG. 1b one gets a multipolar variation of the bubble/tissue interface vibration amplitude Ψ.sub.r(z,ω), shown as 212 in FIG. 2a. We therefore use the transfer function Hv given in Eqs.(A16,A17) to calculate the bubble volume compressibility κ.sub.b at ω.sub.H >> ω.sub.0 for these bubbles. This gives different formulas for the spherical and the cylindrical bubbles as

[0069] Spherical Bubble:

[00018]$\begin{array}{l}{\kappa }_b\left(\omega \right)\approx -\frac{K_d}{{\omega }^2}=-\frac{3}{\rho a^2{\omega }^2}\\ \\ p_s\left(\underset{¯}{r},\omega \right)\approx -a{H}_b\left(\underset{¯}{K}\right)\frac{e^{-ikr}}{r}{P}_i{Y}_s\left(a;\omega \right)~a{H}_b\left(\underset{¯}{K}\right)\end{array}$

[0070] Cylindrical Bubble:

[00019]$\begin{array}{l}{\kappa }_b\left(\omega \right)\approx -\frac{2K_d\left(L/a_v\right)}{\rho a_v^2{\omega }^2}{K}_{pH}\left(L/a_v\right)=\frac{L}{a_v}{K}_d\left(L/a_v\right)=\\ \frac{1}{\pi }{\displaystyle \int dq\frac{q{K}_1\left(q\right)}{K_0\left(q\right)}\frac{{\sin }^{2}qL/a_v}{q^2}}\\ p_s\left(\underset{¯}{r},\omega \right)\approx -a_v{K}_{pH}\left(L/a_v\right)H_b\left(\underset{¯}{K}\right)\frac{e^{-ikr}}{r}{P}_i\\ Y_c\left(a_v;\omega \right)~a_v{K}_{pH}\left(L/a_v\right)H_b\left(\underset{¯}{K}\right)\end{array}$

where it is shown the major dependencies of the receive signal Y(D;ω) on the bubble dimensions D. For the spherical bubble D is the bubble radius a, while for the bubble filling parts of the capillary D is (a.sub.v,L) or (a.sub.v, a) by Eq.(7). A change in the bubble dimension/volume at the HF/LF time reference t.sub.m when the HF pulse hits the bubble, produces a change in the amplitude of the HF receive signal. For polarity p of the LF pulse the major change in bubble dimensions are

[00020]$\text{Sphere:}a=a_0+p\psi \text{Cylindrical:}{a}_v=a_{v0}+p\overline{\psi }$

[0071] The spatial average Ψ of the normal surface displacement is found for the cylinder as described in Eq.(9). The received signal from the HF pulse has angular frequency components in a band Ω.sub.H centered around ω.sub.H. The propagation velocity c.sub.p of forward transmitted HF pulse observes a modification by the co-propagating LF wave as described in U.S. Pats. 9,291,493 ; 11,280,906 ; 11,275,006.

[00021]$c_p=c_0\left(1+p\beta P_L\right)$

where c.sub.0 is the low amplitude linear propagation velocity, P.sub.L is the pressure of the low frequency wave at the location of the HF pulse, β is a nonlinearity parameter of the tissue bulk compression, and p is a scaling parameter defined in relation to Eqs.(10,11). This nonlinear propagation effect gives a non-linear propagation delay (NPD) to the HF pulse, pτ.sub.n(r), and a weaker pulse form distortion (PFD) that increases with the propagation distance r is before scattering, as described in the cited patents. When the transmitted HF pulse is at the crest or through of the LF wave, the NPD has a larger effect than the PFD that generally can be neglected. When the transmitted HF pulse is at a zero crossing of the LF wave, the PFD has a larger effect than the NPD that generally can be neglected.

[0072] The back scattered HF signal will generally in addition to the scattering from the bubbles contain considerable scattering from tissue surrounding the bubble, and it is devised to suppress the tissue signal, for example by correcting the received signal for adequate NPD and/or PFD according to the methods in the cited patents, or using the 2.sup.nd harmonic component of the received HF signal, or a combination of both. We notice, however that the transmitted LF amplitude to vibrate bubbles is generally lower than for imaging, and both the NPD and the PFD have less effect, and corrections for NPD and/or PFD can in many cases be neglected. Before further processing, we assume that the received HF signals are adequately corrected with the NPD and PFD, as described in the cited patents, and in the notation for the received HF signals below, we assume that that adequate correction have been done.

[0073] For illustration of a simplest form, we linearize the receive signals and form a detection signal D.sub.m with p = ±1 as

[00022]$\begin{array}{l}{Y}_{np}\left(\omega \right)\approx Y_0\left(\omega \right)+p\delta Y_n\left(\omega \right)\hfill \\ D_n\left(\omega \right)=\frac{Y_{n-}-Y_{n+}}{Y_{n-}+Y_{n+}}=\frac{1}{Y_0}\frac{\partial Y_n}{\partial a_v}{\psi }_n\hfill \end{array}\delta Y_n=\frac{\partial Y_n}{\partial a_v}{\psi }_n$

where the subscript n indicates the time point t.sub.n in FIG. 4 the measurements are done for, and the subscript n± implies p = ±1. Other values for the polarity parameter p can be used as described in relation to FIG. 4 and Eqs.(10,11). We do in most of the equations operate in the angular frequency domain. As understood by any-one skilled in the art, Eq.(18) and other signal processing can through inverse Fourier transform be carried through in the time domain for a large range of depths.

[0074] Components for Y.sub.s and Y.sub.c in Eqs.(14,15) that have low variation with bubble dimension variations are shorted in the ratio of D.sub.m. In particular has H.sub.b low sensitivity to variations in bubble dimensions that are small compared to the wave length. H.sub.b(K) is the Fourier transform of the bubble volume with Fourier vector K=k(e.sub.i-e.sub.s). The largest magnitude K.sub.max=2ke.sub.i=4π/λ)e.sub.i is found for back-scattering where e.sub.s= –e.sub.i, which the situation in FIG. 5a where the HF transmit (503) and receive (504) beams have the same. With a HF receive cross beam as 514 in FIG. 5b and the HF transmit beam as 503 in FIG. 5a we get K=k|e.sub.i-e.sub.s|<K.sub.max=4π/λ. For a spherical bubble we get

[00023]$\begin{array}{l}\text{Spherical}\text{bubble:}{H}_b\left(\underset{\_}{K}\right)=\frac{3}{{\left(Ka\right)}^3}\left[\sin Ka-Ka\cos Ka\right]\left(a\right)\\ \\ \text{Cylindrical}\text{bubble:}{H}_b\left(\underset{\_}{K}\right)=\frac{2\sin \left(KL\cos \phi \right)}{KL\cos \phi }\frac{J_1\left(K{a}_v\sin \phi \right)}{K{a}_v\sin \phi }\left(b\right)\end{array}$

[0075] We notice that for the sphere we have for back scattering Ka=(4πa/c.sub.p).Math.f where f is the frequency and c.sub.p ≈ 1500 m/s is the propagation velocity of the pressure wave. FIG. 6a shows ΔH.sub.b/H.sub.b as a function of frequency for the spherical bubble where the curves 601,...,605 shows the curves for bubble diameters 2a = 10, 15, 20, 25, 30 .Math.m, and Ψ/a=0.01.ΔH.sub.b/H.sub.b0<0.02 is accepted as a tolerable error, which gives an upper frequency limit where H.sub.b can be accepted as a constant in the processing.

[0076] Scattering from the cylindrical bubble in FIG. 1b gives more complex dependency of H.sub.b with bubble dimensions when the bubble long axis is larger than the capillary diameter. H.sub.b then becomes a function of both the short axis and long axis bubble dimensions, and the direction of K in relation to the bubble, for example described by the angle φ to the bubble long axis in FIG. 1b. FIG. 6b shows H.sub.b as a function of φ. The calculations are done for a capillary radius of a.sub.v = 5 .Math.m, with a bubble cylindrical length 2/.sub.c = 20 .Math.m in 610, 2/.sub.c = 30 .Math.m in 610, 2/.sub.c = 40 .Math.m in 620, 2/.sub.c = 50 .Math.m in 630. Within each of the circular sub-figures we have a frequency from the outer to the inner curve of f = 2.5, 5, 7.5, 10, 15 MHz in each of the Figures. We notice that for K normal to bubble long axis there is negligible change in H.sub.b in the actual frequency range, while for K along the bubble long axis the variation is much larger. To approximate H.sub.b to constant, we first must both have H.sub.b adequately large to get adequate signal in Eq.(14,15,18), while at the same time we must have ΔH.sub.b/H.sub.b0<0.025 to assume H.sub.b constant.

[0077] ΔH.sub.b/H.sub.b0 is analyzed in FIGS. 6c and d for the frequencies f = 2.5, 5, 7.5, 10, 15 MHz that are displayed from the inner to the outer curve in each circular sub-figure. The capillary radius a.sub.v = 5 .Math.m in all Figures, and the cylinder half-length l.sub.c is 10, 15, 20, 25 .Math.m as listed under each circular sub-figure. FIG. 6c shows ΔH.sub.b/H.sub.b0 for for Δa.sub.v=1.Math.m, i.e. Δa.sub.v/a.sub.v=0.2.

[0078] For this situation ΔH.sub.b/H.sub.b0<0.02 is obtained for f = 2.5 and 5 MHz. Reducing □a.sub.v to 0.5 □m we have ΔH.sub.b/H.sub.b0<0.02 up to 7.5 MHz for all situations. FIG. 6d shows ΔH.sub.b/H.sub.b0 for Δl.sub.c=1 .Math.m for all values of l.sub.c. For this situation ΔH.sub.b/H.sub.b0<0.02 is obtained for f = 2.5 and 5 MHz. Reducing Δl.sub.c to 0.5 .Math.m we have ΔH.sub.b/H.sub.b0<0.02 up to 7.5 MHz for all situations.

[0079] We expect capillaries and hence also bubbles to have random position within the tumor, and if we select the bubbles with the strongest back-scatter, i.e. incident wave normal to the capilary length and φ close to ± π/2 eg , we can assume that the variation of H.sub.b for the variations of the bubble dimension produced by the LF pulse is negligible for all the given bubble lengths and frequencies HF < 15 MHz.

[0080] We note that when the incident HF beam hits close to the short axis of the bubble, both H.sub.b0 and ΔH.sub.β/H.sub.b0 are approximately independent of the dimension variations for HF < 15 MHz When we can neglect the changes in H.sub.b with the bubble dimension variations produced by the polarity change in the incident LF pulse, H.sub.b is essentially canceled in the fraction of D.sub.m in Eq.(18). Neglecting H.sub.b in the expression for Y.sub.s, Y.sub.c in Eq.(14,15) we see that D.sub.m of Eq.(18) gives

[00024]$\begin{array}{l}\text{Spherical:}{D}_{sm}=\frac{Y_{sm-}-Y_{sm+}}{Y_{sm-}+Y_{sm+}}=\frac{{\psi }_{sm}}{a}\left(a\right)\\ \text{Cylindrical:}{D}_{cm}=\frac{Y_{m-}-Y_{m+}}{Y_{m-}+Y_{m+}}=A_{cm}^{-1}\frac{{\psi }_{cm}}{a_v}\frac{{\psi }_{cm}}{a_v}=A_{cm}\frac{Y_{m-}-Y_{m+}}{Y_{m-}+Y_{m+}}\left(b\right)\\ A_{cm}^{-1}\approx \frac{1}{K_{pH}\left(L/a_v\right)}\left(K_{pH}\left(L/a_v\right)+a_v\frac{\partial K_{pH}\left(L/a_v\right)}{\partial a_v}\right)\end{array}$

where K.sub.pH(L/a.sub.v) is given in Eq.(15, A17). We note that ψ.sub.cm is the spatial average of the normal displacement ψ.sub.r for the cylinder, as shown in Eq.(A12), while ψ.sub.sm is the maximal displacement of the spherical bubble radius, as the bubble displacement is constant across the bubble surface. A.sub.cm as a function of l.sub.c/a.sub.v is shown as 700 in FIG. 7a. The relation between L and l.sub.c is given in Eq.(7). We notice that for l.sub.c = 0 we have a sphere, and FIG. 7a shows A.sub.cm = 1, which is the value given for the sphere in Eq.(21 a). This gives trust in the curve in the curve 700 for A.sub.cm. Newly evaporated ACT bubbles have a cylindrical length l.sub.c/a.sub.v ~ 3 - 20, depending on capillary radius and degree of gas absorption in the blood. The curve 700 shows that A.sub.cm ≈ 0.55 - 0.65 is a good assessment to obtain estimates of Ψ.sub.cm/a.sub.v with interesting accuracy.

[0081] The relative volume variation δV/V.sub.b can be approximated from the relative displacement for the spherical and cylindrical bubbles with the linear approximations as

[00025]$\text{Spherical:}\frac{\delta V_s}{V_b}\approx 3\frac{{\psi }_{sm}}{a}\text{Cylindrical:}\frac{\delta V_c}{V_b}\approx 2\frac{{\psi }_{cm}}{a}$

[0082] Considering measurements errors and relative displacement magnitudes, these approximate expressions have good practical value. More accurate nonlinear expressions between δV/V.sub.b and the displacements can be calculated by anyone skilled in the art.

[0083] For the 2.sup.nd order transfer functions in Eqs.(5,A14) it is possible to measure/estimate the three major parameters, w.sub.0, ζ, and K.sub.V, for example in the following way. We get a phase estimate θ(ω.sub.L)=∠H(ω.sub.L) of the transfer function at ω.sub.L for the t.sub.n = t.sub.m that gives maximum of δV.sub.b in Eq.(11), i.e.

[00026]${\omega }_L{t}_m+\widehat{\theta }\left({\omega }_L\right)=0.Math.\widehat{\theta }\left({\omega }_L\right)=-{\omega }_L{t}_m$

[0084] FIG. 7b shows an example amplitude (701) and phase (702) of the 2.sup.nd order transfer function H.sub.V(ω). Due to the - sign in front of H.sub.V(ω) in Eqs.(5,A14) the resonance frequency ω.sub.0 the phase of the volume vibration amplitude δV.sub.cb(ω) is – π/2 behind the driving pressure amplitude - P.sub.i, while for ω adequately less than ω.sub.0 the phase-lag is ≈ 0 behind - P.sub.i, i.e. the vibration amplitude has close to opposite polarity of P.sub.i. For ω adequately larger than ω.sub.0, the phase-lag is ≈ - π behind - P.sub.i, i.e. the δV.sub.cb(ω) has the same polarity as P.sub.i. The amplitude and phase of the bubble vibration hence depend on the details of the transfer function (Gain K, resonance frequency ω.sub.0, and absorption ζ).

[0085] Measuring θ̂(ω.sub.L1) and θ̂(ω.sub.L2) from Eq.(22) for two low frequencies ω.sub.L1 and ω.sub.L2 where we get the maximal detection signal for t.sub.m1 and t.sub.m2, we can for example though linear interpolation, indicated as 703 in FIG. 7b, write the phase variation between ω.sub.L1and ω.sub.L2 as

[00027]$\widehat{\theta }\left(\omega \right)=\widehat{\theta }\left({\omega }_{L1}\right)+\frac{\widehat{\theta }\left({\omega }_{L2}\right)-\widehat{\theta }\left({\omega }_{L1}\right)}{{\omega }_{L2}-{\omega }_{L1}}\left(\omega -{\omega }_{L1}\right)$

[0086] At the resonance frequency □.sub.0 we shall have θ(ω.sub.0) = -π/2 which gives an estimated resonance frequency from Eq.(23) as

[00028]${\widehat{\omega }}_0={\omega }_{L1}-\frac{\pi /2+\widehat{\theta }\left({\omega }_{L1}\right)}{\widehat{\theta }\left({\omega }_{L2}\right)-\widehat{\theta }\left({\omega }_{L1}\right)}\left({\omega }_{L2}-{\omega }_{L1}\right)$

[0087] Matching the differential of the phase in Eq.(5,A14) at ω= ω.sub.0 to Eq.(24) we get an estimate for the absorption losses in the tissue vibrations as

[00029]$\zeta =-\frac{1}{{\omega }_0{\theta }^{\prime }\left({\omega }_0\right)}.Math.\widehat{\zeta }=\frac{{\omega }_{L2}-{\omega }_{L1}}{{\widehat{\omega }}_0\left(\theta \left({\omega }_{L2}\right)-\theta \left({\omega }_{L1}\right)\right)}$

[0088] Another estimate of the absorption factor ζ is from Eq.(5) also found for a single low frequency ω.sub.L as

[00030]$\widehat{\zeta }=\frac{1}{2}\left(\frac{{\omega }_L}{{\widehat{\omega }}_0}-\frac{{\widehat{\omega }}_0}{{\omega }_L}\right)\tan \widehat{\theta }\left({\omega }_L\right)$

[0089] The parameter K.sub.V can for example be found by the angular frequency ω.sub.m that gives maximal amplitude of H.sub.V(ω). From the estimates of ω̂.sub.0 and ζ̂ we can obtain an estimate of

[00031]${\widehat{\omega }}_m^2={\widehat{\omega }}_0^2\left(1-2{\widehat{\zeta }}^2\right){\widehat{K}}_v=\left|H_v\left({\widehat{\omega }}_m\right)\right|\sqrt{{\widehat{\omega }}_0^4-{\widehat{\omega }}_m^4}$

where we measure |H.sub.V| at the estimated angular frequency ω.sub.m. If we do not have an estimate of ζ̂, we can search for the ω.sub.m that gives max|H.sub.V| and do the following calculations

[00032]${\widehat{K}}_V=\left|H_V\left({\widehat{\omega }}_m\right)\right|\sqrt{{\widehat{\omega }}_0^4-{\widehat{\omega }}_m^4}\widehat{\zeta }=\sqrt{\frac{{\omega }_0^2-{\omega }_m^2}{2{\omega }_0^2}}$

[0090] Estimating the phases θ̂.sub.n = θ̂(ω.sub.Ln) at more than two frequencies ω.sub.Ln, n=1, 2, ... , N, we can improve the estimate of ω̂.sub.0 with more complex interpolations, as known by any-one skilled in the art. By example, for n = 1, ... , N we the set of measured frequencies by the vector ωL = {ω.sub.L1,...ω.sub.Ln,...ω.sub.LN} and the set of estimated phases by ωL = {ω.sub.L1,...ω.sub.Ln,...ω.sub.LN}, we can then form an estimate of ω̂.sub.0 and ζ̂ω̂.sub.0 through minimization of a mean square functional as

[00033]$\begin{array}{c}j\left({\omega }_0,\zeta ;{\underset{\_}{\omega }}_L\right)={\left[\widehat{\underset{\_}{\theta }}\left({\underset{\_}{\omega }}_L\right)-\underset{\_}{\theta }\left({\omega }_0,\zeta ;{\underset{\_}{\omega }}_L\right)\right]}^T\left[\widehat{\underset{\_}{\theta }}\left({\underset{\_}{\omega }}_L\right)-\underset{\_}{\theta }\left({\omega }_0,\zeta ;{\underset{\_}{\omega }}_L\right)\right]\\ \left({\widehat{\omega }}_0,\widehat{\zeta }\right)=\underset{{\omega }_0,\zeta }{\text{min}}J\left({\omega }_0,\zeta ;{\underset{\_}{\omega }}_L\right)\end{array}$

[0091] Other forms of functionals for minimization, also including added information, can be formed by any-one skilled in the art.

[0092] From Eq.(5-7) we notice that for θ̂(ω.sub.L)+ω.sub.Lt.sub.m≈ 0 the amplitude of the received HF signal takes the form

[00034]$Y{}_{np}=Y_{n0}+p\delta Y_n\sim p\left|H_V\left({\omega }_{Ln}\right)\right|P_{Le}\left(t_n\right)$

[0093] and we can for example include the amplitude of the received signal ~|H.sub.v(ω.sub.Ln)| for the different {ω.sub.L1,...ω.sub.Ln} in addition to the phase in the functional for parameter estimation in Eq.(29).

[0094] The extra capillary tissue is determined by the mass density, ρ.sub.t, and the real and imaginary part, .Math..sub.tr and .Math..sub.ti of the complex shear stiffness of the tissue, as shown in Eq.(A6, A7). We can approximate ρ.sub.t≈ 1100 kg/m.sup.3 for most solid tumor tissues, and hence it is .Math..sub.tr and .Math..sub.ti which are the most sought after tissue parameters. From the measured/estimated ω̂.sub.0 and k̂.sub.v, we see from Eq.(A18) that we can obtain an estimate of ω̂.sub.d as

[00035]${\widehat{\omega }}_d=\sqrt{{\widehat{\omega }}_0^2-\gamma {\widehat{K}}_V{P}_0}{\zeta }_d=\frac{{\omega }_0}{{\omega }_d}\zeta {\widehat{x}}_d=\frac{{\widehat{\omega }}_d}{{\omega }_s}{\omega }_s=\frac{c_s}{a_v}=\frac{1}{a_v}\sqrt{\frac{{\mu }_{tr}}{\rho }}$

where γ ≈ 1 and P.sub.0 ≈ 101 kPa.

[0095] From F.sub.d(x) in Eq.(A11) we obtain estimated x.sub.d(l.sub.c/a.sub.v) from Re{F.sub.d(x.sub.d)=0. FIG. 8a shows x.sub.d(l.sub.c/a.sub.v) for a set of absorption parameters α= 0.5, 1.0, 1.5, defined from η= .Math..sub.ti/.Math..sub.tr=2αx.sub.d. We notice that the variation of the curves with α. Is barely visible, only for l.sub.c/a.sub.v<4. We also notice from Eq.(A10) that

[00036]${\omega }_s=c_s/a_v={\omega }_d/x_d\text{and}{\mu }_{tr}=\rho a_v^2{\omega }_s^2$

which allows us to calculate

[00037]$\frac{c_s}{a_v}\left(l_c/a_v\right)={\omega }_s=\frac{{\omega }_d}{{\widehat{x}}_d\left(l_c/a_v\right)}{\mu }_{tr}\left({\omega }_d\right)=\rho a_v^2\frac{{\omega }_d^2}{{\widehat{x}}_d^2\left(l_c/a_v\right)}$

[0096] FIG. 8b shows ω.sub.s = c.sub.s/a.sub.v as a function of ω̂.sub.d for of l.sub.c/a.sub.v=0,4, and 10 and FIG. 8c shows .Math..sub.tr(ω.sub.d) for a.sub.v=2,4,5, and 6 .Math.m.

[0097] The power dissipated by the surface vibration at a surface area element d.sup.2r is

[00038]$dW=-\frac{1}{2}\text{Re}\left\{d^{2}ri\omega {\psi }_n{\left(\underset{\_}{r}\right)}^{*}{P}_i\right\}=-\frac{\omega }{2}\text{Im}\left\{d^{2}r{\psi }_n\left(\underset{\_}{r}\right)P_i\right\}$

[0098] We assume P.sub.i approximately constant across the bubble surface, and integrated across the whole bubble surface gives the following power dissipation that provides a heating energy per unit time

[00039]$W=-\frac{\omega }{2}\text{Im}\left\{P_i{\displaystyle \int d^{2}r{\psi }_n\left(\underset{\_}{r}\right)}\right\}=-\frac{\omega }{2}\text{Im}\left\{P_i\delta V_{cb}\right\}$

Inserting the transfer function from Eq.(5), we get

[00040]$W\left(\omega \right)=\frac{\omega }{2}{V}_b{P}_i^2\mathrm{Im}\left\{H_V\left(\omega \right)\right\}=\frac{{\omega }_0\zeta K_V{V}_b{\omega }^2/{\omega }_0^2}{\left(1-{\omega }^2/{\omega }_0^2\right)+4{\zeta }^2{\omega }^2/{\omega }_0^2}{P}_i^2$

[0099] Hence, from estimation of ω.sub.0 and ζ and assessment of the bubble volume and the incident pressure amplitude, we can assess the dissipated power from the bubble on the surrounding tissue at ω= ω.sub.L. Because ω.sub.L is low, we can estimate the low frequency amplitude P.sub.i = P.sub.L using an assessment of the low absorption. From the ultrasound image (2D or 3D) we can automatically estimate the number of bubbles N per unit volume and assess the total heat generation as

[00041]$W_T\left({\omega }_L\right)=NW\left({\omega }_L\right)$

that can be used to estimate temperature increase in the tissue to avoid uncontrollable results.

[0100] FIG. 9 shows a block diagram of an example instrument for carrying out the methods according to the invention. 900 shows by example a 3D ultrasound probe comprising a dual frequency linear array 901 with a set of M LF elements and N HF elements in an azimuth direction indicated by the arrows 906. The dual frequency band linear array can be made according to known methods, for example as described in U.S. Pat. 7,727,156; 10,879,867. The LF and HF elements of the array are via a cable 902 connected to a transmit/receive unit 903 that connects each LF array element to LF transmit circuits, and each HF element to HF transmit/receive circuits comprising a HF transmit amplifier and a HF receive amplifier where the output of the HF receive amplifier is further connected to an analog to digital converter (A/D) presenting a digital representation of the HF received signals from all HF receive elements, according to known methods. The AD converter can in a modified embodiment present digital representations of the I-Q components of the HF receive signals from each HF element that represents the same information as the radio frequency (RF) HF signal, according to known methods. In a modified version of the instrument where one wants to receive the scattered LF pulse from the bubble, the LF elements are also connected to LF transmit/receive circuits.

[0101] For 3D scanning of the ultrasound beams, the linear array 901 can in this example embodiment be rotated around the long axis 904 that provides a mechanical scanning of the LF/HF beam in an elevation direction, indicated by the arrows 905. For each elevation position of the array, one does electronic scanning of LF/HF transmit beams in an azimuth direction indicated by the arrows 906, through electronic selection of transmitting LF and HF elements, and transmitting combined LF/HF pulse complexes for example as shown in FIG. 4, with selected beam directions and foci. An example HF transmit beam is illustrated schematically as 907 within a 2D azimuth plane 908 with given elevation position within a total 3D scan volume 909. Alternative elevation movements of the array, like side-ways movement can equivalently be done according to known methods, depending on the space available for such movement, and the shape of the object. Pure 2D scanning of the beam, and even measurements along a single beam direction is also within the scope of the invention. The array might also have a 2D matrix distribution of elements for 3D focusing and scanning of the beams, as known to anyone skilled in the art.

[0102] At least two pulse complexes with different LF pulses, for example as illustrated in FIG. 4, are transmitted for each transmit beam direction. The LF pulse might be zero in one pulse complex per HF transmit beam, but must be non/zero in at least one pulse complex for each HF transmit beam direction.

[0103] Two versions of the instrument are useful, where the first version 903 comprises beam former for HF receive cross-beams, illustrated as 914 in the 2D scan plane 908, and HF back scatter receive beams with the same axis as the HF transmit beam 907, shown in further detail in FIG. 5. In a preferred embodiment the HF back-scatter receive beam is equal to the HF transmit beam as this improves suppression of multiple scattering noise in the HF back-scatter receive signal, as discussed in U.S. Pat. 8,550,998; 9,291,493; 9,939,413; 10,578,725, and synthetic depth optimized focus is obtained through lateral filtering of the HF receive signals. During the scan, the HF cross-beam and back-scatter receive signals are via the high speed bus 910 transferred to the processor 911 for storage and further processing.

[0104] The processor 911 typically comprises at least one multicore central processing unit (CPU) and at least one graphics processor unit (GPU) that are SW programmable, according to known methods. The processor receives user inputs from a user/operator input unit 913 that operates according to known methods, and displays image data and other information necessary for communication with the user/operator through a combined display and audio unit 912, according to known methods.

[0105] In the second version of the instrument, the digital HF receive signals from each HF receive elements and each transmitted pulse complex are via the high-speed bus 910 transferred to the processor 911 for storage and further processing. For 2D imaging in the second version, a SW program in the processor unit 911 combines HF receive signals from multiple HF receive elements to produce both HF backscatter receive beams and HF receive cross-beams crossing each HF transmit beam in the 2D set, for example as shown as 914 and in more detail in FIG. 5b and c. A set of HF back-scatter receive signals from HF back-scatter receive beams with the same axis as the HF transmit beams, and preferably also equal to the HF transmit beams 907, also shown in FIG. 5, can be produced by a software program.

[0106] One use of the cross-beam (914) HF receive signals is to estimate the nonlinear propagation delay at a multitude of depths along the HF transmit beam with low influence of multiple scattered noise, as described in U.S. Pat. Appl. 2020/0405268. A set of cross beams is also useful to provide better suppression of signal from tissue, that is noise for detecting bubbles, vibration amplitudes, and tissue parameters. When we have imaging situations where the level of multiple scattering noise in the back-scatter HF receive signal is low compared to the 1st order scattered signals, the nonlinear propagation delay at multitude of depths along the HF transmit beam (907) can be obtained from the back-scatter HF receive signal, as described in the above cited US Patents, and the beam former for HF receive cross-beams, has less importance and can be omitted.

[0107] The processor 911 further comprises means, for example implemented as SW programs, to carry through the processing as described above to estimate [0108] i) the relative variation of the bubble volume or displacement amplitude when driven by the LF signal, as for example described in Eqs.(19,21) and following equations, and [0109] ii) the bubble vibration resonance frequency and absorption, as for example described in Eqs.(22-25, 30) and related equations, and [0110] iii) parameters of shear elasticity of the extravascular tissue, for example as shown in FIGS. 8b and c.

[0111] The method and instrument is useful for estimating vibration amplitudes of bubbles of any dimension at an angular frequency we define as the low frequency, ω.sub.L. For spherical bubbles we have a good analytic relationship of a 2.sup.nd order transfer function, Eq.(5) between the bubble radius and the scattered wave over a large frequency range from well below to well above the resonance frequency. For Bubbles that fill a capillary, like in FIG. 1b, we have in Appendix A develop an analytic expression for H.sub.v(w) in Eq.(A12) with an approximate 2.sup.nd order transfer function for frequencies close to the volume resonance frequency in Eq.(A14). Eqs.(22 - 32) gives methods to identify the parameters of said model and the tissue surrounding a bubble that fills part of a capillary.

[0112] Thus, while there have shown and described and pointed out fundamental novel features of the invention as applied to preferred embodiments thereof, it will be understood that various omissions and substitutions and changes in the form and details of the devices illustrated, and in their operation, may be made by those skilled in the art without departing from the spirit of the invention.

[0113] It is also expressly intended that all combinations of those elements and/or method steps which perform substantially the same function in substantially the same way to achieve the same results, are within the scope of the invention. Moreover, it should be recognized that structures and/or elements and/or method steps shown and/or described in connection with any disclosed form or embodiment of the invention may be incorporated in any other disclosed or described or suggested form or embodiment as a general matter of design choice. It is the intention, therefore, to be limited only as indicated by the scope of the claims appended hereto.

APPENDIX A

[0114] For the vibration field outside the cylindrical capillary from the vibrating bubble in FIG. 1b, we use the cylindrical (r,z) coordinates as shown in the Figure. The incompressible field in the vabrating co-oscillating mass outside the capillary can for vibration of a cylindrical bubble with length 2L at an angular frequency ω, be approximated by the Laplace equation for a potential A(r,z;ω)

[00042]${\nabla }^2\text{A}\left(r,z;\omega \right)+k^2\text{A}\left(r,z;\omega \right)\underset{\lambda >>2L}{.fwdarw.}{\nabla }^2\text{A}\left(r,z;\omega \right)=0$

where λ = 2πc.sub.p/ω is the wavelength of the pressure wave with propagation velocity c.sub.p. The particle vibration displacement ψ and the pressure P in the co-oscillating mass are obtained from A as

[00043]$\underset{\_}{\psi }\left(r,z;\omega \right)=\nabla \text{A}\left(r,z;\omega \right)P\left(r,z;\omega \right)=\rho {\omega }^{2}A\left(r,z;\omega \right)$

[0115] Fourier transform along the z-axis gives

[00044]$\begin{array}{l}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial A\left(r,z\right)}{\partial r}\right)+\frac{{\partial }^{2}A\left(r,z\right)}{\partial z^2}\underset{F{T}_z}{.fwdarw.}\\ \frac{{\partial }^{2}A\left(r,k_z\right)}{\partial r^2}+\frac{1}{r}\frac{\partial A\left(r,k_z\right)}{\partial r}-k_z^{2}A\left(r,k_z\right)=0\end{array}$

[0116] Setting x = k.sub.zr and A(r,k.sub.z) .sub.= y(k.sub.zr) we get the modified Bessel equation with the solution K.sub.0(x)

[00045]$x^2{y}^{″}\left(x\right)+xy\prime \left(x\right)-x^{2}y\left(x\right)=0x=k_{z}ry\left(x\right)=K_0\left(x\right)$

[0117] With A defined for r= a.sub.v and a tissue mass density ρ, we get the following field amplitudes in the tissue outside the capillary

[00046]$\begin{array}{l}\begin{array}{cc}A\left(r,k_z\right)=\frac{A\left(a_v,k_z\right)}{K_0\left(k_z{a}_v\right)}{K}_0\left(k_{z}r\right)& P\left(r,k_z\right)=p{\omega }^{2}A\left(r,k_z\right)\end{array}\\ {\psi }_r\left(r,k_z\right)=\frac{\partial A\left(r,k_z\right)}{\partial r}=-k_z\frac{A\left(a_v,k_z\right)}{K_0\left(k_z{a}_v\right)}{K}_1\left(k_{z}r\right)\\ {\psi }_z\left(r,k_z\right)=-i{k}_{z}A\left(r,k_z\right)=-i{k}_z\frac{A\left(a_v,k_z\right)}{K_0\left(k_z{a}_v\right)}{K}_0\left(k_{z}r\right)\end{array}$

[0118] Inside the co-oscillating mass the incompressible tissue deformation produces shear stresses, which balances to zero inside the tissue, but at the boarder to the bubble at the capillary wall we get a stress contribution to the tissue surface as

[00047]$\begin{array}{l}{\sigma }_{\mu r}\left(a_v,k_z\right)=2{\mu }_t\frac{\partial {\psi }_r}{\partial r}=2{\mu }_t\frac{{\partial }^{2}A}{\partial r^2}=-2{\mu }_t\left(\frac{1}{r}\frac{\partial A}{\partial r}-k_z^{2}A\right)\\ {\sigma }_{\mu r}\left(a_v,k_z\right)=2{\mu }_t\left(k_z^2+\frac{k_z{K}_1\left(k_z{a}_v\right)}{a_v{K}_0\left(k_z{a}_v\right)}\right)A\left(a_v,k_z\right)\end{array}$

where .Math..sub.t = .Math..sub.tr + i.Math..sub.ti = .Math..sub.tr(1 +iη) is the complex shear stiffness of the tissue where the imaginary component produces absorption of oscillation power. At the bubble-tissue surface interface this stress adds to the pressure in Eq.(A5) to the following pressure from the cylindrical bubble surface towards the tissue

[00048]$\begin{array}{l}{P}_d\left(a_v,k_z\right)=\rho {\omega }^{2}A\left(a_v,k_z\right)-{\sigma }_{ur}\left(a_v,k_z\right)=\rho H_{\omega }\left(\omega ;k_z{a}_v\right)A\left(a_v,k_z\right)\\ H_{\omega }\left(\omega ;k_z{a}_v\right)={\omega }^2-2{\omega }_c^2\left({\left(k_z{a}_v\right)}^2+\frac{k_z{a}_v{K}_1\left(k_z{a}_v\right)}{K_0\left(k_z{a}_v\right)}\right)\end{array}$

[00049]${\omega }_c^2=\frac{{\mu }_t}{\rho a_v^2}{\mu }_t={\mu }_{tr}\left(1+i\eta \right){\omega }_c^2={\omega }_s^2\left(1+i\eta \right){\omega }_s^2=\frac{{\mu }_{tr}}{\rho a_v^2}$

[0119] When the bubble cylindrical length 2/.sub.c is adequately larger than the capillary radius a.sub.v, we can approximate the bubble with semi-spherical end surfaces with radius a.sub.v with a cylinder with length 2L as shown in Eqs.(7, 8) above. P.sub.d (a.sub.v,k.sub.z is the z-Fourier transform of a spatially constant drive pressure amplitude P.sub.d relative to infinity in the interval (- L, L), and zero outside

[00050]$\begin{array}{l}{P}_d\left(a_v,k_z\right)=2\frac{\sin k_{z}L}{k_z}{P}_{d}A\left(a_v,k_z\right)=\frac{P_{d\left(a_v,k_z\right)}}{\rho H_{\omega }\left(\omega ;k_z{a}_v\right)}=\frac{2}{\rho H_{\omega }\left(\omega ;k_z{a}_v\right)}\frac{\sin k_{z}L}{k_z}{P}_d\\ {\psi }_r\left(a_v,k_z\right)=-\frac{k_z{K}_1\left(k_z{a}_v\right)}{K_0\left(k_z{a}_v\right)}A\left(a_v,k_z\right)=-\frac{k_z{K}_1\left(k_z{a}_v\right)}{K_0\left(k_z{a}_v\right)}\frac{2}{\rho H_{\omega }\left(\omega ;k_z{a}_v\right)}\frac{\sin k_{z}L}{k_z}{P}_d\\ {\psi }_r\left(a_v,z\right)=-\frac{1}{2\pi }{\displaystyle \int d{k}_z\frac{k_z{K}_1\left(k_z{a}_v\right)}{K_0\left(k_z{a}_v\right)}\frac{2}{\rho H_{\omega }\left(\omega ;k_z{a}_v\right)}\frac{\sin k_{z}L}{k_z}}{e}^{-i{k}_{z}z}{P}_d\end{array}$

[0120] The ψ.sub.r vibration amplitude at the cylindrical part of the bubble amplitude, gives a volume change of the cylinder

[00051]$\begin{array}{l}\frac{\delta V_c}{V_b}=\frac{4\pi \text{}{a}_{v}L}{2\pi \text{}{a}_v^{2}L}\frac{1}{2L}{\displaystyle \underset{-L}{\overset{L}{\int }}dz{\psi }_r}\left(a_v,z\right)=\\ -\frac{2}{\rho a_v^2}\frac{1}{\pi }{\displaystyle \int d{k}_z\frac{k_z{a}_v{K}_1\left(k_z{a}_v\right)}{K_0\left(k_z{a}_v\right)}}\frac{1}{H_{\omega }\left(\omega ;k_z{a}_v\right)}\frac{{\sin }^2{k}_{z}L}{k_z^{2}L}{P}_d\end{array}$

[0121] Changing integration coordinates q = k.sub.za.sub.v, k.sub.z = q/a.sub.v , dk.sub.z = dq/a.sub.v we get the direct drive transfer function H.sub.d as

[00052]$\frac{\delta V_c}{V_b}=H_d\left(\omega ;L/a_v\right)\frac{P_d}{{\mu }_{tr}}$

[00053]$\begin{array}{l}{H}_d\left(\omega ;L/a_v\right)=-\frac{1}{\pi }{\displaystyle \int dq}\frac{q{K}_1\left(q\right)}{K_0\left(q\right)}\frac{2{\omega }_s^2}{H_{\omega }\left(\omega ;q\right)}\frac{{\sin }^{2}qL/a_v}{q^{2}L/a_v}\\ {\omega }_s^2=\frac{{\mu }_{tr}}{{}^{\rho a_v^2}}=\frac{c_s^2}{a_v^2}{c}_s^2=\frac{{\mu }_{tr}}{\rho }\end{array}$

where c.sub.s is the shear wave velocity and ω.sub.s is the shear wave angular frequency. Introducing x = ω/ ω.sub.s we obtain the scaled direct drive transfer function F.sub.d as

[00054]$\begin{array}{l}\frac{\delta V_c}{V_b}=F_d\left(x;L/a_v\right)\frac{P_d}{{\mu }_{tr}}\\ F_d\left(x;L/a_v\right)=-\frac{1}{\pi }{\displaystyle \int dq}\frac{q{K}_1\left(q\right)}{K_0\left(q\right)}\frac{2}{H_x\left(x;q\right)}\frac{{\sin }^{2}qL/a_v}{q^{2}L/a_v}\\ H_x\left(x;q\right)=x^2-2\left[q^2+q{K}_1\left(q\right)/K_0\left(q\right)\right]\left(1+i\eta \right)\end{array}$

[0122] Note that for the cylindrical bubble in a capillary we have above put the surrounding tissue real shear stiffness .Math..sub.tr outside H.sub.d and F.sub.d. This gives F.sub.d the nice property that it does not depend on .Math..sub.tr and ω.sub.s, only on η. We assume a linear variation η with frequency, i.e. η = 2αx. The resonance frequency of F.sub.d is defined as the value x.sub.d that gives the real part of F.sub.d equal zero, i.e. Re{Fd(xd;L/av}. FIG. 3b shows xd(lc/av) as a function of lc/av. The angular resonance frequency of H.sub.d is ω.sub.d = ω.sub.s,x.sub.d . The variation of xd with lc/av is quite low and we can approximate xd ≈ 1 for lc/av in the range 1 to 8. This allows us to approximate the direct drive resonance function □This allows us to approxlc/av. For lc/av = 0 we get the direct drive values for the spherical bubble as xd ≈ 1.67 and □≈ 1.67 and

[0123] Including gas elasticity according to Eq.(3) we get the complete transfer function from incident pressure to relative volume amplitude

[00055]$\begin{array}{l}\frac{\delta V_c}{V_b}=H_d\left(\omega ;L/a_v\right)\frac{P_d}{{\mu }_{tr}}=-H_v\left(\omega ;L/a_v\right)P_i\\ H_V\left(\omega ;L/a_v\right)=\frac{H_d\left(\omega ;L/a_v\right)/{\mu }_{tr}}{1-\gamma H_d\left(\omega ;L/a_v\right)P{}_0/{\mu }_{tr}}\end{array}$

[0124] FIG. 3 shows that we for a practical interval around the resonance frequency can approximate the transfer function with a 2.sup.nd order transfer function with ω.sub.d = ω.sub.sx.sub.d

[00056]$\begin{array}{l}{\widehat{F}}_d\left(x;L/a_v\right)\approx \frac{-K_F}{x^2-x_d^2-i2{\zeta }_d{x}_{d}x}\\ {\widehat{H}}_d\left(\omega ;L/a_v\right)\approx \frac{-K_d}{{\omega }^2-{\omega }_d^2-i2{\zeta }_d{\omega }_d\omega }{\omega }_d^2={\omega }_s^2{x}_d^2{K}_d={\omega }_s^2{K}_F\end{array}$

[0125] Due to the frequency variation of the displacement along the bubble surface as shown in 212 of FIG. 2a, K.sub.F has to be adapted numerically, for example so that the maximal amplitude of the F.sub.d approximation in Eq.(A13) matches to the maximal amplitude in Eq.(A11). From Eq.(A12) we see that with a given incident amplitude P.sub.i we can include the variations in the gas pressure on the bubble volume to get the complete 2.sup.nd order transfer function as

[00057]$\begin{array}{l}{\widehat{H}}_V\left(\omega ;L/a_v\right)=\frac{{\widehat{H}}_d\left(\omega ;L/a_v\right)/{\mu }_{tr}}{1-\gamma {\widehat{H}}_d\left(\omega ;L/a_v\right)P_0/{\mu }_{tr}}\approx \frac{-K_V}{{\omega }^2-{\omega }_0^2-i2\zeta {\omega }_0\omega }\\ K_V=\frac{K_d}{{\mu }_{tr}}=\frac{K_F}{\rho a_v^2}{\omega }_0^2={\omega }_d^2+\gamma K_V{P}_0={\omega }_s^2\left(x_d^2+K_F\left(L/a_v\right)P_0/{\mu }_{tr}\right)\\ \zeta ={\zeta }_d{\omega }_d/{\omega }_0{K}_V={\omega }_s^2{K}_F/{\mu }_{tr}\end{array}$

[0126] With the 2.sup.nd order approximation of H.sub.v in Eq.(A14) we can from Eq.(13) approximate the scattered field in the low frequency range as

[00058]$\begin{array}{l}{\kappa }_b\left(\omega \right)=-\frac{1}{P_i}\frac{\delta V_{cb}}{V_b}=H_V\left(\omega \right)p_s\left(\underset{\_}{r},\omega \right)=-H_p\left(\omega \right)H_b\left(\underset{\_}{K}\right)\frac{e^{-ikr}}{r}{P}_i\\ H_p\left(\omega \right)=\frac{K_{pL}{\omega }^2}{{\omega }^2-{\omega }_0^2-i2\zeta {\omega }_0\omega }{K}_{pL}=\frac{1}{2}{\rho }_t{a}_v^{2}L{K}_V=\frac{1}{2}L{K}_F=\frac{1}{2{\omega }_s^2}L{K}_d\end{array}$

[0127] For detection of the low frequency LF vibration amplitude we transmit a high frequency (HF) pulse, HF ~ 10*LF or more, i.e. ω>> ω.sub.0, as described in relation to Eqs.(12-19) and FIG. 4. Eq.(A12) includes these HF variations of the bubble surface displacement amplitude as illustrated in 212. We can then approximate the bubble volume compressibility from Eq.(A12) for ω>> ω.sub.0

[00059]$\begin{array}{l}{\kappa }_b\left(\omega \right)=H_V\left(\omega \right)\frac{1}{{\mu }_{tr}}\approx H_d\left(\omega \right){\mu }_{tr}\approx -\frac{2K_d\left(L/a_v\right){\omega }_s^2}{{\omega }^2}{\mu }_{tr}=-\frac{2}{\rho a_v^2{\omega }^2}{K}_d\left(L/a_v\right)\\ K_d\left(L/a_v\right)=\frac{1}{\pi }{\displaystyle \int dq\frac{q{K}_1\left(q\right)}{K_0\left(q\right)}}\frac{{\sin }^{2}qL/a_v}{q^{2}L/a_v}\end{array}$

[0128] Inserted into Eq.(13) we then get the approximation to the scattered field for the Ω.sub.H frequency range as

[00060]$p_s\left(\underset{\_}{r},\omega \right)\approx -a_v{K}_{pH}\left(L/a_v\right)H_b\left(\underset{\_}{K}\right)\frac{e^{-ikr}}{r}{P}_i$

[00061]$K_{pH}\left(L/a_v\right)=\frac{L}{a_v}{K}_d\left(L/a_v\right)=\frac{1}{\pi }{\displaystyle \int dq}\frac{q{K}_1\left(q\right)}{K_0\left(q\right)}\frac{{\sin }^{2a}qL/a_v}{q^2}$

[0129] The received HF amplitude is then proportional to

[00062]$Y_r\left(a_v,L,\underset{\_}{K}\right)\sim a_v{K}_{pH}\left(L/a_v\right)H_b\left(\underset{\_}{K}\right)$

[0130] With the 2.sup.nd order approximation of H.sub.v in Eq.(A14) the angular frequency ω.sub.m for max|H.sub.v (ω)| is given by

[00063]$\begin{array}{ccc}\zeta =-\frac{1}{{\omega }_0\theta \prime \left({\omega }_0\right)}& {\omega }_m={\omega }_0\sqrt{1-2{\zeta }^2}& \text{}\\ \zeta =\sqrt{\frac{{\omega }_0^2-{\omega }_m^2}{2{\omega }_0^2}}& \left|H_v\left({\omega }_m\right)\right|=\frac{K_V}{\sqrt{{\omega }_0^4-{\omega }_m^4}}& K_V=\frac{K_d}{{\mu }_{tr}}=\frac{K_F}{\rho a_v^2}\end{array}$

REFERENCES

[0131] 1. Lars Hoff: “Acoustic Characterization of Contrast Agents for Medical Ultrasound Imaging” Springer Science + Business Media. 2001