## Abstract

A new theoretical model of ultrasound propagation in soft biological media is presented based on an extended thermodynamics formalism. The long-standing experimental conjecture claiming that a continuous distribution of internal degrees of freedom can be used to model ultrasound in biological media is given theoretical justification. A strategy to derive a well-defined set of equations coupling the balance equations of mass, momentum, energy and entropy with relaxation kinetics of a medium characterized by a continuous distribution of internal states is presented. We demonstrate that new phenomenological coefficients of the proposed governing equations can be extracted directly from experimental data. Our theory successfully explains the anomalous attenuation law found in experiments with biological media that is inconsistent with the conventional models using a finite number of internal degrees of freedom. The results presented offer new possibilities for medical applications of high-intensity ultrasound and ultrasound emission methods to study matter with complex internal structure. These techniques include using pressure relaxation methods for accurate investigation of fast protein folding and a variety of other applications for media where irreversible thermodynamic simulations are essential.

## 1. Introduction

Interaction of ultrasound with soft biological tissues plays a fundamental role in biophysics and its biomedical applications [1]. The thermodynamics of this interaction has recently attracted much interest in connection with the mechanisms of thermal ablation of oncological tumours [2]. Currently, there is no reliable theoretical model describing the response of tissue to ultrasound. Although a variety of approaches have been proposed [3–5], none of them provides thermodynamic analysis of the rheological relations and the wave damping mechanisms used. Consequently, to the best of the author's knowledge, it is still an open question whether the earlier-mentioned approaches agree with the fundamental physical principle that the rate of entropy production must be positive. An explicit proof of this principle is essential for the self-consistency of a physical theory aiming to describe media with complex laws of sound attenuation, for example, the power law, because for these media, its fulfilment does not follow from the formalisms of conventional thermodynamics and continuum mechanics. As is known [6], violation of this principle can result in physically unrealistic thermodynamic instabilities. In addition to this, the knowledge of the entropy is often required in applications. Thus, it is of importance to explore how to construct a theory that results in a well-defined set of self-consistent nonlinear governing equations. The aim of this study is to examine rational formulations of the problem of ultrasound propagation in soft biological media that follow from the thermodynamic method.

A number of phenomenological approaches towards the construction of non-equilibrium thermodynamics have been documented in the literature (see Lebon *et al*. [6] for a detailed discussion). However, it still remains a challenge to understand the nature of the irreversible transformations accompanying the propagation of nonlinear ultrasound in tissue.

Experiments on ultrasound attenuation in tissues reveal a broadband statistical distribution of internal degrees of freedom over the whole range of relaxation times available for observations, with the major part of the attenuation arising at a molecular, rather than cellular, or higher, level of organization [7–9]. It was then suggested that a continuous distribution of internal parameters should be used to model the dynamics of the medium because a finite number of internal states is not sufficient to derive the correct law of attenuation. Although such an approach has found applications in the thermodynamics of proteins [10,11], an irreversible thermodynamic theory of ultrasound propagation in media with continuous distributions of internal states has not been formulated.

We address this question in the present study from the viewpoint of Prigogine's [10] extended thermodynamic formalism, whose core feature is to define the specific entropy *s* with the help of the equation
1.1where *T*, *e*, *V* =1/*ρ* and *p* are the absolute temperature, the specific internal energy, the specific volume and the pressure, *ρ* is the density; *t* is time and parameter *τ* characterizes possible different internal states. For determinacy, here *τ*≥0. *ξ*(*τ*) is the density of molecules in the state *τ*. Function *ξ*(*τ*) satisfies the kinetic equation , where is the number of molecules that change their state from *τ* to *τ*+d*τ* per unit time, *x*(*τ*) is the chemical potential.

In a more general thermodynamic context, one observes that variables *x* and d*ξ*/d*t* in formula (1.1) can be considered as a continuously distributed thermodynamic flux and its conjugate generalized force [6]. Prigogine's construction then becomes particularly suitable for the problem of sound propagation in tissue because it is in line with the formulated experimental conjecture and requires only one partial differential equation to capture the relaxation kinetics. It turns out, remarkably, that the strategy outlined in Prigogine's analysis results in consistent sets of governing equations for sound propagation in tissue.

In theories of sound attenuation with a finite number of internal degrees of freedom, the internal parameter *x* is traditionally defined either as the irreversible pressure *p*^{r}, or the irreversible portion of the internal energy *e*^{r} associated with the related internal relaxation process. The first situation is typical for conventional acoustics [12,13], the second one for physico-chemical hydrodynamics [14,15]. We consider the continuous versions of these two scenarios in what follows.

## 2. Pressure relaxation scenario

The problem of ultrasound attenuation in soft biological matter is very complex, and the details of physico-chemical events initiated in tissue by ultrasound are largely unknown. An in-depth discussion of this topic can be found in the book by Hill *et al.* [1]. Here, we present only some of the important findings that motivate the analytical approach proposed.

According to the available biochemical data [16,17], application of an external pressure to a protein medium results in rapid relaxation kinetics whose typical time scales are of the order of microseconds and are in the same range as the periods of ultrasound waves used in biomedical applications. Typical pressure wave amplitudes used in ultrasonics are also comparable with the forcing pressures required to initiate kinetic reactions of protein denaturation.

The existing evidence [18,19] suggests that the denatured protein is not a single uniform state, but is a large statistical ensemble of molecular conformation states with similar free energies. These conformation states can rapidly establish equilibrium with one another, and are highly sensitive to thermodynamic perturbations of the medium. Consequently, it has been proposed in the thermodynamic literature that a statistically large number of microscopic transitions in protein structures and their influence on the overall behaviour of the whole medium as a continuum can be described using a continuous distribution of internal states [10,11].

On the hypothesis that an ultrasound wave initiates the myriad of individual kinetic reactions, a continuous distribution of internal thermodynamic parameters can be used to explain the mechanism of ultrasound attenuation in soft tissue. As pointed out in §1, there is also abundant evidence in the literature that a finite number of internal degrees of freedom is not sufficient to explain the law of ultrasound attenuation in soft tissue. For this reason, a continuous set of internal parameters is introduced in the following, in order to construct the thermodynamic theory of soft tissue ultrasound.

In the approach developed in this section, these internal parameters are chosen to be the irreversible relaxation pressures that result in the medium as a consequence of the related internal relaxation kinetics. The set of internal parameters is assumed to be continuously distributed over the internal states *τ*>0. Here, is the position vector.

The remaining thermodynamic state variables used in this section are the density and the specific entropy . As opposed to the energy relaxation scenario considered in §4, here the specific internal energy is taken in the quasi-equilibrium form. It is assumed to be a real-valued function *e*=*e*(*ρ*,*s*) of the density and the entropy, but not the internal parameters , whereas the total pressure *p*_{Σ} is the sum of the pressure *p*(*ρ*,*s*) given by the equilibrium thermodynamic equation of state and the total relaxation pressure *I*_{p},
2.1Here, function *φ*_{A}(*τ*)>0 describes the distribution of relaxation pressures over the internal states *τ*. Only small relaxation corrections to the related equilibrium variables are considered in what follows, with *p*^{r}=0 at thermodynamic equilibrium.

For better visualization of the ideas that form the basis of the proposed extended thermodynamic theory and for clarity of presentation, only the case of scalar, rather than tensorial, internal parameters is considered in this study.

The following balance equations of mass, momentum and energy are assumed to hold in the medium:
2.2and
2.3where is the velocity vector, *J*_{NS} and *σ*_{NS} stand for the entropy flux and entropy production owing to heat conduction and viscosity, exactly as in the conventional Navier–Stokes–Fourier equations. These well-studied effects present no interest here and owing to shortage of space will be discarded.

In order to formulate a closed set of governing equations, the entropy balance equation and a transport equation for the internal parameter *p*^{r} which guarantees that the entropy production is positive should be added to equations (2.1)–(2.3).

In this scenario, the equation of entropy balance is defined by means of Prigogine's formula (1.1), in which we take *x*=*p*^{r}(*τ*). Formula (1.1) then yields the following balance equation for the specific entropy *s*:
2.4where *T* is the quasi-equilibrium temperature, that is, the derivative of the specific internal energy *e*(*ρ*,*s*) with respect to the entropy at constant density.

The term *ρσ*^{r} in equation (2.4) is the rate of entropy production owing to the relaxation motion,
2.5whereas the first term in formula (2.4) is zero by virtue of the energy balance equation (2.3).

The term in parenthesis in formula (2.5) is the *total* thermodynamic force related to the thermodynamic flux *p*^{r}. The last term in this parenthesis is its part associated with temporal variations of the density, d*ρ*/d*t*. The remaining part of the total thermodynamic force d*ξ*/d*t* does not depend on d*ρ*/d*t* because the variation of the state variable *ρ* with time has already been accounted for by the last term in parenthesis in formula (2.5). Consequently, d*ξ*/d*t* depends only on temporal variations of the internal variable *p*^{r}, i.e. d*p*^{r}/d*t*, and of the state variable d*s*/d*t*. We present the remaining part of thermodynamic force d*ξ*/d*t* by the following linear form:
2.6where *α* and *β* are phenomenological coefficients characterizing the medium.

The total thermodynamic force corresponding to the thermodynamic flux *p*^{r} is given by the expression
In line with the methodology of irreversible thermodynamics, the thermodynamic flux *p*^{r} and the related force *X*^{r} are identically zero in the state of thermodynamic equilibrium.

Thermodynamics conjectures the existence of a linear force–flux relation that we take in the following form:
2.7where *γ* is an unknown coefficient. Additional terms can, in principle, appear on the right-hand side of equation (2.7). However, at present, there is no experimental data to justify their presence.

The entropy production term may be transformed by a partial integration into the following form:
2.8provided *γ*⋅(*p*^{r})^{2}=0 at *τ*=0 and . In order for the entropy production to be positive, we require that
2.9

Function ∂*γ*/∂*τ* has the physical meaning of the rate of change of the adiabatic modulus of the medium in the configurational space *τ*. This change results from the excitation of the internal degrees of freedom as a consequence of the change in volume in response to the application of the external pressure.

Formula (2.7) yields the following kinetic equation for the function *p*^{r}:
2.10with the initial conditions *p*^{r}(*t*≤0)=0 and *p*^{r}(*τ*≤0)=0. We assume here that the coefficient *γ*/*α* is positive.

The material (phenomenological) function *m*_{ρ} determines the intensity of the pressure relaxation kinetics induced by the perturbation of the initial thermodynamic state via application of an external pressure field. These kinetic transformations are responsible for the dispersion and attenuation of sound in the medium. Consequently, the material function *m*_{ρ} is the measure of sound–speed dispersion in the reacting medium. It describes the overall effect of individual microscopic kinetic reactions taking place in the medium on the variation of the speed of sound and on the attenuation with the wave frequency. These effects are discussed in detail in Vilensky *et al*. [20].

The phenomenological coefficient *α* has the physical meaning of the change in the adiabatic compressibility with time as a consequence of the pressure-initiated relaxation kinetics in the material of the medium. The physical meaning of the function *β* is less straightforward. It can be shown to describe the simultaneous change in the coefficient of thermal expansion and the heat capacity ratio.

Equation (2.10) shows that the temporal change of relaxation pressure *p*^{r} results from the action of the external force that perturbs the medium from the state of equilibrium (the term on the right-hand side), and the pressure variation that results from the medium's transformations in the internal configurational space *τ*. The coefficient *γ*/*α* determines the rate of this variation. In principle, it can depend on *τ*, if the intensity of relaxation phenomena changes with *τ*. If this effect is negligible, than *γ*/*α*≈*constant* and, without loss of generality, we can then take *γ*/*α*=1 by normalizing *τ* to *γ*/*α*. In this case, parameter *τ* has the dimension of time and serves as the replacement for discrete relaxation times *τ*_{k} of the conventional relaxation theory with a finite number of internal states [12,13].

A practical implementation of the theoretical approach developed in this section requires the knowledge of function *φ*_{A}. Bearing in mind that true distribution of relaxation times in media with complex internal structure, such as tissue, is unknown, it makes sense to choose function *φ*_{A} in the simplest possible way that is sufficient for explanation and quantitative modelling of the available experimental results. We take *φ*_{A}=1/*τ*, in order for the integral term in formula (2.1) to be of the same dimension as the remaining terms.

In the particular problem of ultrasound propagation in soft tissue, equation (2.10) can be further simplified due to the fact that the motion is close to isentropic. Indeed, it follows from formulae (2.4) and (2.8) that d*s*/d*t* is *O*((*p*^{r})^{2}), where the relaxation pressure *p*^{r} itself is a small parameter. Thus, the last term in equation (2.10) is of higher order than the remaining ones and can be neglected to a good approximation. Then, the only phenomenological function that remains undetermined is *m*_{ρ}. It satisfies the restriction ∂(*τm*_{ρ})/∂*τ*>0 that follows from formula (2.9) and the fact that *γ*=*α*=*φ*_{A}/(*ρ*^{2}*m*_{ρ}).

In principle, it should be possible to determine the material function *m*_{ρ} from the pressure profile measurements obtained in the experimental biochemical studies of pressure-induced relaxation kinetics, such as those reported in Jacob et al. [16] and Pearson et al. [17]. This could be particularly interesting and important for biochemical applications in view of the fact that other kinetic measurements follow the pressure profiles. However, this approach requires the development of a rather complicated mathematical method of analysis of such experimental results, and it is not pursued in this study.

A more straightforward way to determine function *m*_{ρ} is to use the conventional approach based on the sound attenuation measurements that are reported in the literature for a wide range of soft tissues. An analytical procedure of its reconstruction from sound attenuation data is described in detail in Vilensky *et al.* [20]. It constitutes one of the central results of the cited paper and involves complicated mathematical techniques that are not given here owing to a shortage of space. For the reader's convenience, and in order to clarify how thermodynamic theory based on an infinite number of degrees of freedom can be used to describe power-law attenuation, we briefly illustrate its derivation from the proposed thermodynamic theory in §2*a*.

### (a) Power-law attenuation

Consider the following one-dimensional adiabatic version of the governing system (2.1), (2.2), (2.10) linearized about the steady state *ρ*=*ρ*_{o}, *p*_{Σ}=*p*_{o}, *u*=0, *p*^{r}=0:
2.11
2.12Here, the subscripts *t* and *x* stand for the partial derivatives with respect to time *t* and spatial coordinate *x*; *ρ*_{o}, *p*_{o} and *c*_{so} are the density, the pressure and the speed of sound in the medium at rest. Functions *U*, *R*=*ρ*−*ρ*_{o}, *P*_{Σ}=*p*_{Σ}−*p*_{o} and *P*^{r} describe the perturbations of the velocity, the density, the total pressure and the relaxation pressure *p*^{r}, respectively. As discussed earlier, we have taken *φ*_{A}=1/*τ* and *γ*=*α*=*φ*_{A}/(*ρ*^{2}*m*_{ρ}). The phenomenological function *m*_{ρ} is taken in the form
2.13where 0<*b*<1 and *C*_{1}>0, to ensure that the entropy production rate is not negative.

We assume that the medium remains at rest for all *t*≤0 and is set in motion by means of the time-periodic excitation of the velocity field on the left boundary *x*=0,
2.14where *ω* is the excitation frequency.

In order to simplify the formulated problem, it is advantageous to eliminate function *P*^{r} from the system of equations (2.11) and (2.12) to express the total relaxation pressure *I*_{p} in terms of the velocity *U*.

The Laplace transform of the second equation in (2.12) with respect to the configurational variable *τ* results in the following well-known first-order relaxation equation:
2.15where the Laplace transform variable *ϖ* plays the role of the relaxation frequency. However, unlike the classical relaxation theory, now *ϖ* is not restricted to a discrete set of parameters, but is distributed continuously over the positive half axis *ϖ*≥0. Consequently, a continuum of relaxation processes corresponding to relaxation times equal to 1/*ϖ* is considered here. Functions
2.16are the Laplace transforms of functions *P*^{r} and *m*_{ρ}, respectively; *Γ*(*z*) is the Euler gamma function. Variable *P*_{r} has the physical meaning of the relaxation pressure per unit relaxation frequency *ϖ*.

Because
it follows from the first formula in (2.12) that the total relaxation pressure *I*_{p} is simply the sum of individual relaxation pressures per unit relaxation frequency, *P*_{r}, over all permissible *ϖ*, that is,
2.17

According to equations (2.15), we can write
2.18where, in agreement with the first equation of system (2.11), the function *R*_{t} on the right-hand side of equation (2.15) has been expressed in terms of *U*_{x}. The total relaxation pressure *I*_{p} corresponding to the infinite number of degrees of freedom characterized by the distribution is then given by the formula
2.19which can be obtained by substitution of expression (2.18) in equation (2.17).

Substitution of formula (2.19) into the system of equations (2.11) and subsequent elimination of the functions *R* and *P*_{Σ} results in the following equation for the velocity *U*:
2.20is a non-dimensional parameter.

Mathematical analysis of equation (2.20) can be greatly simplified if parameter *ϵ* is assumed to be small. In this case, parameter *ϵ* can be used to construct the related formal asymptotic expansion of the solution *U*.

Because here we are not interested in the transient behaviour of the sound field associated with the process of ‘switching on’ of the excitation (2.14), the analysis of the problem can be further restricted to the study of the solution in the limit of sufficiently large time *t*. In this case, we can seek for the leading-order term of the asymptotic solution of equation (2.20) for *ϵ*≪1 and in the form
2.21where the symbol ‘Re’ denotes the real part of a complex expression, i is the imaginary unit and the dots stand for higher-order terms; *α*_{r} and *α*_{i} are unknown constants.

These constants can be determined by substituting expression (2.21) in equation (2.20), discarding all the terms of order greater than *O*(*ϵ*) and taking the limit . This procedure results in the following expressions for *α*_{r} and *α*_{i}:

Because, in agreement with formula (2.21), the attenuation coefficient *α*_{at} is equal to the product *α*_{r}*ϵ*, we see that the thermodynamic model with an infinite number of degrees of freedom described in this section exhibits attenuation with the following power-law frequency dependence:
2.22

## 3. Nonlinear ultrasound waves in soft tissues

To illustrate how the proposed theory can be used to describe ultrasound in soft biological media, we take the following simplified one-dimensional adiabatic version of the governing system of equations derived in §2:
3.1
3.2Here, *w*=*ρu* is the mass flux, *u* is the velocity along the *x* axis, *p*=*p*(*ρ*) is the adiabatic equation of state for the medium in thermodynamic equilibrium. We take Tait's adiabatic equation of state,
3.3where *n*−1 is the so-called nonlinearity parameter. Because, in soft tissue, *n* is usually much larger than unity, the equation of state *p*=*p*(*ρ*) is the main source of the nonlinearity in the governing system of equations. Further analysis of various equations of state currently available in the literature and their influence on the properties of a signal's wave form can be found in Vilensky *et al.* [21]. The derivation of expression (3.2) for the total relaxation pressure *I*_{p} can be found in Vilensky *et al.* [20].

Here, for the purpose of numerical solution, the system of governing equations has been written in the conservation form (3.1). In this form, it can be integrated numerically with the help of the well-documented fully discrete central difference total variation diminishing scheme described in Kurganov & Tadmor [22]. This numerical scheme was developed by its authors to construct solutions of nonlinear systems of conservation laws and of a wide range of nonlinear convection–diffusion equations. The detailed analysis of its stability and approximation properties can be found in the cited paper. These issues are not discussed here any further as our numerical tests of accuracy and stability have shown to be in agreement with the predictions of the above authors. In accordance with the methodology of Kurganov & Tadmor [22], the diffusive term represented by the partial derivative of the convolution integral ∂*I*_{p}/∂*x* in (3.1) is treated explicitly as the diffusion flux. At each discrete moment in time *t*_{n} function, *I*_{p} is approximated by the following quadrature sum:
where

As an illustration, we consider the propagation of a one-dimensional ultrasound wave in a biological system consisting of two layers of tissue with substantially different physical properties. An initially sinusoidal wave is excited at the left boundary of the computational domain and propagates through a strongly nonlinear and strongly absorbing layer of fatty tissue followed by a less nonlinear and less absorbing, though more dense layer of liver tissue. According to Duck [23], in fatty tissue, *n*−1=9.63, *ρ*_{o}=1.02 g cm^{−3}, *c*_{so}=1430 m s^{−1}, and *n*−1=7.75, *ρ*_{o}=1.093 g cm^{−3}, *c*_{so}=1537 m s^{−1} in liver tissue. In the simulation, we take the following power laws of attenuation: *α*_{at}=0.6*f*^{1.4} and *α*_{at}=0.4*f*^{1.14} dB cm^{−1} for fatty tissue and liver, respectively. Here, *f* is the wave frequency in megahertz. There is a sharp interface between the two layers, across which the unperturbed density, the speed of sound, the equilibrium equations of state and the attenuation change abruptly. The interface is 0.7 cm away from the left boundary. The ultrasound wave is generated by the periodic forcing applied at the left boundary of the computational domain. Here, *A*=14.3 m s^{−1} is the amplitude of the wall velocity *u*. This magnitude of *A* roughly corresponds to about 20 MPa of the wall pressure amplitude. The excitation frequency is *f*=1 MHz. This value is in line with those used in many typical applications of medical ultrasound.

Figure 1 depicts typical pressure and density profiles at an instant of time *t* equal to several wave periods after the system was set in motion. The solid line in figure 1 is the pressure profile obtained by solving numerically the full nonlinear system (3.1)–(3.3). The dotted line represents the corresponding density profile. For reference, the pressure profile of the related linear system of equations is also plotted in the same figure with a dash-dotted line.

It can be inferred from figure 1 that a system of nonlinear waves rapidly develops in both media. Comparison of the linear (dash-dotted line) and the nonlinear solution (solid line) shows that the nonlinearity is manifested in the steepening of the nonlinear pressure wave in the course of its propagation and a somewhat more rapid decay of its amplitude, owing to the damping effect of the nonlinearity [24]. Contrary to this, the rate of decay of the linear pressure wave is solely due to the tissue attenuation properties.

As the pressure wave passes through the tissue interface at *x*=0.7 cm, its amplitude grows. This fact has the following explanation. The characteristic acoustic impedance *Z*_{o}=*ρ*_{o}*c*_{so} of liver tissue is higher than the characteristic acoustic impedance of fatty tissue, with their ratio being *ς*=1.15. In agreement with the Fresnel equations for the reflection and transmission coefficients *R*=(*ς*−1)/(*ς*+1) and *T*=1+*R*, respectively, in a linear lossless medium, the amplitude of the transmitted wave must be larger than the amplitude of the incident wave for the chosen configuration of tissue types. Although attenuation and nonlinearity bias this effect to some extent, qualitatively it still manifests itself in figure 1 in the fact that the amplitude of the transmitted wave is slightly higher than that of the approaching wave.

Although the one-dimensional character of the computation precludes direct comparison of these results with measurement data, qualitatively they are in agreement with experimental observations.

## 4. Energy relaxation scenario

Finally, we conduct the analysis of the second scenario of continuous distribution of internal degrees of freedom, in order to extend the discrete hydrodynamic relaxation theory of physico-chemical type.

Following the existing tradition [14,25] the state variables are now taken to be the density *ρ*, the temperature *T* and a set of internal parameters . Here, variables are assumed to be continuously distributed over the internal states *τ*>0.

Function is the excess of the internal energy of the non-equilibrium relaxation process that corresponds to the internal state associated with the variable *τ* over its equilibrium value. The total internal energy *e*_{Σ} is the sum of the equilibrium internal energy defined by the equilibrium caloric equation of state *e*(*ρ*,*T*) and the total relaxation energy *I*_{e},
4.1where function *φ*_{A}(*τ*) now describes the distribution of energies *e*^{r} with respect to the variable *τ*. As above, only small relaxation corrections to the related equilibrium variables are considered here, and *e*^{r}=0 at thermodynamic equilibrium.

Following Zel'dovich & Raizer [14], the temperature *T* is defined as the equilibrium temperature.

In contrast to the pressure relaxation scenario, the total pressure is given by the equilibrium thermodynamic equation of state 4.2As in §2, the dynamics of the medium is described by the equations of continuity, momentum, energy and entropy, 4.3and 4.4and the two equations of state (4.1) and (4.2). In this scenario, we take .

To close this system, a kinetic equation for the internal parameters that describes the time evolution of the non-equilibrium internal processes in the medium must be added. Because these processes result in the production of entropy, the related kinetic equation must satisfy the principle of positive entropy production.

The required equation can be obtained using the same line of reasoning as in the pressure relaxation scenario. However, now the generalized force d*ξ*/d*t* in the entropy equation (4.4) is taken to be proportional to the rates of change of the state variables d*e*^{r}/d*t*, d*ρ*/d*t* and d*T*/d*t*,
4.5

Here, function depends on *τ* and the state variables *ρ* and *T*. It characterizes the intensity of the medium's dynamic response to perturbation of its mean thermodynamic state by the external forcing field; *β*^{r} is a phenomenological coefficient. Functions and *β*^{r}(*τ*,*ρ*,*T*) are determined in the following.

As in the previous scenario, functions d*ξ*/d*t* and *e*^{r} must be linked by a linear force–flux relation such that the resulting rate of entropy production is a positive quadratic functional in *e*^{r}. Here, we take
4.6

The physical meaning of the proportionality coefficients *L*, *β*^{r} and the quantity are discussed later.

Substitution of the third formula from (4.3) and formula (4.6) into formula (4.4) and subsequent partial integration in *τ* results in the equation
4.7provided that *L*⋅(*e*^{r})^{2}=0 at *τ*=0 and . We require that
4.8in order that the rate of entropy production is positive.

Formulae (4.5) and (4.6) give the sought kinetic equation for the unknown function *e*^{r},
4.9Together with formulae (4.1)–(4.3) and (4.7), it forms a closed set of equations, provided the phenomenological functions *L*, *β*^{r} and are known.

Here, the ratio *L*/*β*^{r} plays the same role as its counterpart *γ*/*α* in the previous scenario. Because the configurational space has been defined for non-negative values of parameter *τ*, we require that *L*/*β*^{r}>0. In this case, the slopes of the characteristic lines of equation (4.9) are always positive, and the system does not leave its configurational space in finite time. As in §2, below we consider the most simple situation when *L*=*β*^{r} and *φ*_{A}=1/*τ*. The exception is §4*a*, where the case of only one relaxation mode is analysed, and the distribution function *φ*_{A} is known exactly.

### (a) Physical interpretation of material functions

To gain a better insight into the physical meaning of the phenomenological parameters *L*(*τ*,*ρ*,*T*), *β*^{r}(*τ*,*ρ*,*T*) and , it is useful to compare them with the related parameters of the hydrodynamic relaxation theory of a non-equilibrium gas [14].

For this purpose, let us define function *E*^{r}(*τ*,*ρ*,*T*) as the sum of the quantity and the increment of the non-equilibrium internal energy *e*^{r}(*τ*,*ρ*,*T*),
4.10It follows from equation (4.9) that function *E*^{r}(*τ*,*ρ*,*T*) satisfies the kinetic equation
4.11whereas the expression for the total internal energy (4.1) and the entropy equation (4.4) take the following form:
4.12and
4.13

The structure of equations (4.12) and (4.13) is analogous to the structure of the related expressions of the relaxation theory of non-equilibrium gas [14]. In this theory function, *E*^{r}(*τ*,*ρ*,*T*) plays the role of the non-equilibrium-specific internal energy of the internal degree of freedom corresponding to the parameter *τ*, and the quantity is its equilibrium value. Function *e*_{o}=*e*_{o}(*ρ*,*T*) can be interpreted as the equilibrium internal energy less the equilibrium internal energy of the degrees of freedom of the relaxation processes.

For example, when only one non-equilibrium mode is present, the equations for gas can be obtained from formulae (4.11)–(4.13) by taking
4.14where *τ*_{o} is a constant proportional to some typical value of the relaxation time *t*_{rel}(*ρ*,*T*)>0, *δ*(*τ*) is the Dirac delta-function, and *T*^{r} are the specific heat capacity and the non-equilibrium temperature associated with the non-equilibrium mode.

The temperature *T*^{r} and the entropy are governed by the following well-known equations:
4.15whose detailed analysis can be found in Zel'dovich & Raizer [14]. Because only small deviations from the state of equilibrium are considered here, we assume that the temperature *T*^{r} of the non-equilibrium mode is close to the equilibrium temperature *T*, and is taken to be a constant, for simplicity.

This example shows that the phenomenological function plays the role of the equilibrium internal energy of the internal degree of freedom corresponding to the parameter *τ*, whereas the ratio *L*/*β*^{r} characterizes the deviation of the rate of the related microscopic kinetic process from its typical value owing to the changes in the temperature, the density and, possibly, the configurational variable *τ*. This follows from the fact that, according to the third formula in (4.14), the quantity *L*/*β*^{r} is equal to the ratio of the typical value of the relaxation time *τ*_{o} to the actual value of the relaxation time *t*_{rel}(*ρ*,*T*).

The quantity (−*β*^{r}*e*^{r}) is analogous to the affinity, which in the above example is equal to the term before the derivative d*T*^{r}/d*t* on the right-hand side of the entropy equation (4.15) [25]. It is worth noting, however, that although in chemical thermodynamics and kinetics, the affinities are expressed in terms of the concentrations of the chemical components or the microscopic reaction rates and, at least in principle, can be determined experimentally [26,27], already in the simplest case of ideal systems, the related analytical formulae are essentially empirical.

The earlier-mentioned example demonstrates that even in the most simple physical situation the determination of the physical meaning of the phenomenological parameters, their microscopic interpretation and evaluation requires the availability of an appropriate statistical mechanical theory and the information on the microscopic kinetic rates of the involved relaxation dynamics. In the presented example, these determine the values of and of the relaxation time, respectively.

In ultrasonics, the microscopic interpretation of the material (phenomenological) parameters *L*(*τ*,*ρ*,*T*), *β*^{r}(*τ*,*ρ*,*T*) and depends on the knowledge of the functional role of tissue biomolecules in the process of sound propagation, their spatial structure, types of chemical bonds between the atoms, interaction of protein structures with the surrounding molecules, etc. Unfortunately, advances in understanding of the detailed events underlying these very complex processes are still limited. Despite the awareness that the empirically observed features of tissue result from changes in its molecular structure under the action of the imposed ultrasound field, there is currently no satisfactory approach to analysing these changes and to constructing the quantitative theory. Even the qualitative understanding of the nature of the involved molecular transformations is fragmentary. Consequently, the methods that can be used to give the microscopic interpretation or to measure the material parameters of soft tissue are limited in number. This is particularly true in the case of soft tissue ultrasonics, where the typical physical time scale is defined by the sound wave period and is usually very small (of the order of 1 μs and less).

For these reasons, comparison of the present theory with the conventional theory of relaxation processes in gases was used in order to better appreciate the physical meaning of the involved phenomenological parameters.

Further insights into the structure of parameters *L*(*τ*,*ρ*,*T*) and *β*^{r}(*τ*,*ρ*,*T*) can be gained by noting that their values must depend on the way the relaxation degrees of freedom are distributed in the material and that both these functions have the dimension reciprocal of specific internal energy. Consequently, it can be expected that *L*(*τ*,*ρ*,*T*) and *β*^{r}(*τ*,*ρ*,*T*) are inversely related to either *E*^{r}(*τ*,*ρ*,*T*) or . Because it has been assumed that the relaxation correction to the equilibrium internal energy is small, either of the above functions can be taken so that, for example,
4.16

In principle, some experimental information about the phenomenological parameters *L*(*τ*,*ρ*,*T*), *β*^{r}(*τ*,*ρ*,*T*) and can be obtained from calorimetric measurements. However, the pressure-initiated kinetics are very rapid compared with the temperature-induced kinetics. Consequently, such measurements are unlikely to be of practical use for tissue ultrasonics, because the degrees of freedom characterized by fast relaxation kinetics either would have come to equilibrium or remained unexcited on the thermal time scale.

Because this work is concerned with the problems of ultrasound in tissue, and in the absence of other well-established measurement methods, the unknown functions *L*, *β*^{r} and will be determined following the traditional approach of ultrasonics based on the use of the experimental law of attenuation.

### (b) Determination of phenomenological coefficients

In order to obtain the expressions for the functions *L*, *β*^{r} and and also to compare the two scenarios, we express phenomenological parameters *L* and in terms of function *m*_{ρ} whose link with the law of attenuation of linear ultrasound waves has already been discussed. To do this, we exclude the equilibrium internal energies *e*(*ρ*,*s*) and *e*(*ρ*,*T*) from equation (2.3) and the third equation in (4.3) by making use of the thermodynamic identity
where *Γ* is the Gruneisen coefficient, *c*_{s} is the adiabatic speed of sound. After some simple algebra, this results in the equations
for the first and the second scenarios, respectively. It is straightforward to verify now that in the case of linear motion, the two approaches become identical, if we set
4.17where *ρ*_{o} and *Γ*_{o} stand for *ρ* and *Γ* in the medium at rest. Thus, phenomenological functions *L*, *β*^{r} and can be determined from the experimental law of attenuation.

Formulae (4.17) also show that the non-dimensional proportionality coefficient in formula (4.16), which was obtained earlier from dimensional considerations, is equal to the Gruneisen parameter *Γ*, that is,
4.18

## 5. Conclusions

We have developed a thermodynamic theory of ultrasound propagation in media where the simultaneous excitation of a continuous range of internal degrees of freedom can take place. We have shown that the long-lasting conjecture claiming that a continuous distribution of internal states can be used to model the dynamics of ultrasound waves in biological media has rigorous theoretical justification in terms of the extended thermodynamic formalism. More importantly, we have established that the two typical theoretical scenarios traditionally used to describe media with finite number of internal degrees of freedom have well-defined extensions for media with a continuous configurational space. We have shown that for linear motion, these new models are equivalent to each other and agree well with the experiments on ultrasound attenuation in soft tissue. In particular, it turns out, remarkably, that the well-known experimental power law of attenuation in soft tissue automatically guarantees that the rate of entropy production is positive within our theoretical approach. To the best of the author's knowledge, previous work on ultrasound in tissue has not systematically addressed this side of the problem, although the issue of attenuation is fundamental for biological media, and there exists a whole set of experimental attenuation laws in the literature [20]. Our theory offers a new possibility for accurate analysis of the consistency of experimentally derived material relationships for media with complicated internal structure as these must be in agreement with the theoretical restrictions imposed on the related phenomenological coefficients.

These observations suggest that the presented theory should be fundamental for ultrasound emission techniques used to study biological matter. These include new tools for biophysical studies of rapid kinetics of protein reactions such as pressure relaxation methods for fast protein folding [11,28] and thermodynamic simulations of media with complex internal structure [29].

## Acknowledgements

The author thanks Prof. V. Shrira for discussions and acknowledges the financial support of the EPSRC under grant no. EP/F0257501/1.

- Received July 6, 2012.
- Accepted December 14, 2012.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.