Nikolai Botkin
email: botkin@ma.tum.de

Varvara Turova
email: turova@ma.tum.de

Simulation of acoustic wave propagation in anisotropic media using dynamic programming technique

Keywords: WKB-approximation, Hamilton-Jacobi equations, viscosity solutions, differential game, unification

It is known that the Hamiltonian of the eikonal equation for an anisotropic medium may be nonconvex, which excludes the application of Fermat's minimum-time principle related to minimum-time control problems. The idea proposed in this paper consists in finding a conflict control problem (differential game) whose Hamiltonian coincides with the Hamiltonian of the eikonal equation. It turns out that this is always possible due to Krasovskii's unification technique. Having such a differential game allows us to apply dynamic programming methods to computing the value function of the game, and therefore to describe the propagation of wave fronts. This method is very appropriate for the simulation of wave patterns in surface acoustic wave biosensors. Numerical computations given in this paper prove the feasibility of the approach proposed.

Introduction

The paper concerns the development of methods for modeling the propagation of acoustic waves in anisotropic media. This investigation is very important for many applications such as acoustic sensors whose operating principle is based on the excitation and detection of acoustic waves of very high frequency in piezoelectric crystals.

For anisotropic media, the WKB (Wentzel-Kramers-Brillouin) approximation yields eikonal equations whose Hamiltonians are neither convex nor concave in the impulse variable. Therefore, the well-known Fermat principle of wave propagation fails in this case. Moreover, the propagation occurs in such a way as if an antagonistic opponent aims to slow down the movement of the wave fronts. Thus, we arrive at the idea to use methods of differential games in the analysis of wave propagation. If the negative of the Hamiltonian of a differential game approximates the Hamiltonian of the eikonal equation, then the value function of Corresponding author.

the game approximates the phase function satisfying the eikonal equation. The authors have developed effective and precise algorithms for solving Hamilton-Jacobi equations arising from differential games, which yields an effective tool for the numerical investigation of eikonal equations.

If a differential game is chosen appropriately, the level sets of its value function represent the wave fronts, and optimal trajectories are associated with the propagation of rays. Thus, it makes possible to describe very complicated behavior of rays using game-theoretic classification of the so-called singular surfaces that can attract, repulse, and break the trajectories. For example, the caustic-like behavior of rays can be interpreted as the attraction of neighboring optimal trajectories to a singular surface. The monograph [START_REF] Melikyan | Generalized Solutions of First Order PDEs[END_REF] extends the classical method of characteristics by introducing the so-called generalized characteristics that are related to the above mentioned trajectories and singular surfaces.

The main objective of this paper is the numerical simulation of propagation of bulk and surface acoustic waves in anisotropic monocrystals and multi-layered structures used in surface acoustic wave sensors. It is demonstrated that the propagation fronts can be found very precisely even in the case of very complicated geometry of wave emitters. Numerical results are presented for the case of bulk and surface waves characterized by non-convex slowness surfaces. This investigation is inspired by the cooperation with professor A. A. Melikyan (deceased) from the Institute for Problems in Mechanics, Moscow, Russia.

Wave velocity in piezoelectric crystals and Eikonal equation

Assume that the indexes i, j, k, l run from 1 to 3 and use the summation convention over the repeated indexes. Let u 1 , u 2 , and u 3 be the displacements in x 1 , x 2 , and x 3 directions, respectively; ϕ is the electric potential such that the electric field E i is given by the relation E i = ∂ϕ/∂x i . Electro-elasticity equations for a piezoelectric anisotropic crystal read:

ρu itt -C ijkl ∂ 2 u l ∂x j ∂x k + e kij ∂ 2 φ ∂x k ∂x j = 0, ( 1
) ij ∂ 2 φ ∂x i ∂x j + e ikl ∂ 2 u l ∂x i ∂x k = 0 . (2)
where ρ, ij , e ikl , and C ijkl denote the density, the material dielectric tensor, the stress piezoelectric tensor, and the elastic stiffness tensor, respectively. The WKB (high frequency) approximation (see e.g. [START_REF] Barles | Remarks on a Flame Propagation Model[END_REF]) uses the ansatz

u j = u 0 j (t, x) • εe ıS(t,x)/ε , φ = φ 0 (t, x) • εe ıS(t,x)/ε , (3)
where ε = ω -1 is a small parameter (ω is the frequency), S(t, x) is the phase function, and u 0 j (t, x) and φ 0 (t, x) are functions defining the wave polarization. The symbol ı in the exponent denotes the imaginary unit. Substitution of (3) into (1) and ( 2) and collection of the terms of order 1/ε yield the equations

-ρS 2 t δ il + C ijkl ∂S ∂x j ∂S ∂x k u 0 l -e kij ∂S ∂x k ∂S ∂x j φ 0 = 0, (4)
-e ikl ∂S ∂x i ∂S ∂x k u 0 l -ij ∂S ∂x i ∂S ∂x j φ 0 = 0 . (5)
The condition of nontrivial solvability of the system (4) and ( 5) leads to the eikonal equation det

     1 ρ      C ijkl ∂S ∂x j ∂S ∂x k -e kij ∂S ∂x k ∂S ∂x j -e ikl ∂S ∂x i ∂S ∂x k -ij ∂S ∂x i ∂S ∂x j      - δ il 0 3×1 0 1×3 0 S 2 t      = 0,
which can be rewritten as

S t -|∇S|c α ∇S |∇S| = 0, (6)
where

c α (n), α = 1, 2, 3, are eigenvalues of the problem det   1 ρ   C ijkl n j n k -e kij n k n j -e ikl n i n k -ij n i n j   - δ il 0 3×1 0 1×3 0 c 2   = 0.
Here, n 1 , n 2 , n 3 , |n| = 1, are components of the normalized wave vector (the direction of propagation). Therefore, for each vector n, there are three types of waves propagating in this direction. Each of them has its own velocity c α and the corresponding nontrivial solutions, u 0 l and ϕ 0 , of (4) and (5) defining the wave polarization.

Figure 1a shows the phase velocity surface for LiT aO 3 piezoelectric crystals. This surface is obtained as the set of points of the form c α (n) • n, where n belongs to a grid on the surface of the unit sphere. The index α corresponds to a quasi shear wave where the displacements are near orthogonal to the propagation direction. Figure 1b presents the so-called slowness surface consisting of points of the form c -1 α (n) • n. It is easy to prove that the slowness surface can be described as {p ∈ R 3 : c α (p/|p|)|p| = 1}. The nonconvexity of the slowness surface shows that the Hamiltonian, c α (p/|p|)|p|, of equation ( 6) is non-convex in p.

Surface acoustic wave biosensors

Biosensors serve for the measurement of small amounts of biological substances in liquids. Usually a biosensor can be considered as a multi-layered structure (see Fig. 2) whose bottom layer is the ST-cut of piezoelectric α-quartz. Acoustic shear waves are excited here by means of a high-frequency voltage applied to electrodes placed on the ST-cut surface. The waves are transmitted into an isotropic guiding layer deposited on the top of the quartz substrate. The top gold layer is covered by DNA or RNA molecules, aptamers, that are able to specifically bind protein molecules from the contacting liquid. Binding protein molecules results in additional mass loading, which causes a phase shift in the electric signal measured by the output electrodes. A high sensitivity regarding to the added mass is achieved due to the usage of shear horizontally polarized guided waves (Love waves) because of their low interaction with the contacting fluid. The input and output electrodes are located between the substrate and the guiding layer. To obtain purely shear polarized modes, the direction of the wave propagation is chosen to be orthogonal to the crystalline X-axis (see Fig. 3). The next three subsections consider a mathematical model of the biosensor and two methods of numerical investigation of acoustic Love waves including their phase velocity, decay with depth, polarization, etc.

Mathematical model of biosensor

The governing equations for the displacements and the electric potential in the quartz substrate are given by formulae (1) and [START_REF] Barles | Remarks on a Flame Propagation Model[END_REF].

The gold layer is conductor so that there is no electric field inside it. The electric field inside the guiding layer is also neglected because of its low dielectric permeability. Therefore, the electric potential vanishes, and the gold and guiding layers are described by the equation of the form

¯ u i tt -Cijkl ∂ 2 u l ∂x j ∂x k = 0. (7)
In the fluid layer, the Stokes and mass conservation equations hold:

0 v i t -ν∆v i -(ζ + ν 3 ) ∂ ∂x i div v + ∂ ∂x i P = 0, ( 8
)
γP t + ∂ ∂x i v i = 0,
where v i are components of the velocity, P is the pressure, 0 is the fluid density at a reference pressure P 0 , ν and ζ are the dynamic and volume viscosities of the fluid, respectively, and γ = 1

0 ∂ ∂P P0
is the compressibility of the fluid.

A special homogenization technique developed in [START_REF] Hoffmann | Homogenization of Interfaces Between Rapidly Oscillating Fine Elastic Structures and Fluids[END_REF] is used to treat the aptamer-fluid structure. This bristle structure is replaced by an averaged material whose properties are derived as the number of bristles goes to infinity, their thickness tends to zero, and the height remains constant. The resulting new layer whose thickness is equal to the height of the aptamer emulates the aptamer-fluid structure. The governing equation for this layer is given by the relation (see [START_REF] Hoffmann | Homogenization of Interfaces Between Rapidly Oscillating Fine Elastic Structures and Fluids[END_REF] and [START_REF] Botkin | Dispersion Relations for Acoustic Waves in Heterogeneous Multi-Layered Structures Contacting with Fluids[END_REF])

ˆ u i tt -Ĉijkl ∂ 2 u l ∂x j ∂x k -Pijkl ∂ 2 u lt ∂x j ∂x k = 0, (9)
where the term containing the tensor P describes the viscous damping coming from the liquid part of the aptamer-fluid structure. The term containing Ĉ represents elastic stresses. The density ˆ is a weighted combination of the densities of the fluid and the aptamer. The tensors P and Ĉ are computed with FE-method using an analytical representation of solutions of the so-called cell equation arising in homogenization theory.

The conditions on the interfaces between the layers are carefully considered in [START_REF] Botkin | Mathematical Models of a Biosensor[END_REF] and [START_REF] Botkin | Two Complementary Approaches in Modelling a Biosensor[END_REF]. Briefly, the continuity of the displacements and the equilibrium of the normal pressures must hold on the interface between every two neighboring solid layers (the averaged aptamer-fluid layer is considered as solid). Moreover, the electric displacement and the tangent component of the electric field in the substrate must be zero on the interface between the quartz substrate and the guiding layer. The conditions on the interface between the aptamer layer and the fluid include the no-slip assumption and the equilibrium of the pressures.

Finite element modeling

The FE-model extends the above described basic model by accounting for two alternated groups of electrodes (see Figure 2a) and a damping area around the side and bottom faces to suppress the wave reflection thereon.

The electrodes are typically made of gold. Therefore, they can be accounted for by the linear elasticity equation of type [START_REF] Botkin | Numerical Computation of Dispersion Relations for Multi-Layered Anisotropic Structures[END_REF].

Accounting for the damping is done by adding the term -div(β(x)∇u i t ) to equations ( 2), [START_REF] Botkin | Numerical Computation of Dispersion Relations for Multi-Layered Anisotropic Structures[END_REF], and [START_REF] Royer | Elastic Waves in Solids I: Free and Guided Propagation[END_REF], where β(x) is a piecewise-linear function which is equal to zero outside of the damping region and grows up to some value β 0 > 0 towards the side and bottom faces.

The FE-approach provides accurate results because of accounting for the exact parameters of the sensor such as the shape of the electrodes, their position, mass, electro-conductivity properties. This allows us to estimate important characteristics of the biosensor and effects caused by scattering of waves (see [START_REF] Botkin | Mathematical Models of a Biosensor[END_REF] and [START_REF] Botkin | Two Complementary Approaches in Modelling a Biosensor[END_REF] for simulation results).

The main difficulty of this approach is very high resource-consuming because of a very small wavelength. A large number of finite elements in x 1 -direction is required to resolve the wave structure. The number of degrees of freedom lies in the range of 10 6 -10 7 , which makes impossible, e.g. to compute the phase wave surface with appropriate accuracy.

Harmonic analysis (Dispersion relations)

The approach related to the harmonic analysis is developed in [START_REF] Botkin | Dispersion Relations for Acoustic Waves in Heterogeneous Multi-Layered Structures Contacting with Fluids[END_REF] and provides a method for the construction of travelling wave solutions feasible in the biosensor structure under the assumption of its unboundedness in the lateral and downward directions. This assumption is very realistic because real biosensor chips are imbedded up to the surface in very viscous damping media that suppresses the reflection of waves on the the side and bottom faces.

The algorithm described here quite briefly (see [START_REF] Botkin | Numerical Computation of Dispersion Relations for Multi-Layered Anisotropic Structures[END_REF] and [START_REF] Botkin | Dispersion Relations for Acoustic Waves in Heterogeneous Multi-Layered Structures Contacting with Fluids[END_REF] for more details). It is assumed that all the layers are infinite in x 1 and x 2 directions, the (top) fluid layer and the (bottom) substrate layer are semi-infinite in x 3 direction. The electrodes are not taken into account.

We are looking for solutions describing plain waves propagating in x 1 direction. This means that the displacements in the solid layers, the velocities in the fluid, and the electric potential in the substrate are of the form:

u i (x 1 , x 3 ) = a i (x 3 ) cos(κx 1 -ωt) + b i (x 3 ) sin(κx 1 -ωt), (10)
v i (x 1 , x 3 ) = c i (x 3 ) cos(κx 1 -ωt) + d i (x 3 ) sin(κx 1 -ωt), (11)
ϕ(x 1 , x 3 ) = f (x 3 ) cos(κx 1 -ωt) + g(x 3 ) sin(κx 1 -ωt), (12)
where κ is the wave number, and ω is the circular frequency which is equal to the frequency of the voltage applied to the input electrodes in our case. Substitute [START_REF] Crandall | Viscosity Solutions of Hamilton-Jacobi Equations[END_REF] and ( 12) into (1) and (2) for the substrate; (10) into ( 7) and ( 9) for non-piezoelectric layers and for the aptamer layer; and (11) into [START_REF] Auld | Acoustic Fields and Waves in Solids[END_REF] for the fluid layer. Collecting all coefficients of cos and sin yields a system of ordinary linear differential equations for the coefficients a i , b i , c i , d i , f , g in each layer. Solving these systems for every layer, we obtain the representation of the functions a i , b i , c i , d i , f , g in the following form (only the expression for the function a = (a 1 , a 2 , a 3 ) is given here because the form of the other functions is similar):

a(x 3 ) = j D j h j e λ j κx3 , (13)
where D j are arbitrary coefficients, λ j and h j are eigenvalues and eigenvectors of the matrix of the corresponding system of differential equations. For the semiinfinite fluid and substrate layers, only terms decreasing towards x 3 for the fluid and towards -x 3 for the substrate, i.e. terms with negative Reλ j for the fluid and positive Reλ j for the substrate, are kept. Every layer has its own set of coefficients D j , eigenvalues λ j , and eigenvectors h j . To find any particular travelling wave solution in the whole structure we need to determine the coefficients D j for each layer, which is being done by substituting the expressions of the form (13) for the functions a i , b i , c i , d i , f , g into ( 10)-( 12) and then the resulting functions u i , v i , ϕ into the interface conditions outlined at the end of subsection 3.1 (see [START_REF] Botkin | Two Complementary Approaches in Modelling a Biosensor[END_REF] for exact description of the interface conditions). Since all of the interface conditions are linear relations, the computation yields a homogeneous system of linear equations for the unknown coefficients D j . Denote by G(ω, κ) the matrix of this system. Fix the circular frequency ω and denote the unknown phase velocity by V = ω/κ to consider G as a function of V . The phase velocity is feasible if and only if the system has a nontrivial solution, which is equivalent to the condition det G T (V )G(V ) = 0, where G T (V ) the conjugate transpose of G(V ). The last equation can be easily solved because the computation of the left-hand-side runs very quickly even on an ordinary computer. Usually, there are several roots corresponding to different types of waves propagating with different phase velocities. Concerning the biosensor, the root corresponding to a shear wave, i.e. only u 2 = 0, is to be chosen.

Figure 4a shows the computed phase velocity contour for surface acoustic waves exited in the biosensor structure using the excitation frequency of 96 MHz. Figure 4b presents the slowness contour scaled by 10 3 . It is seen that the slowness contour is not convex (see remark at the end of Section 2).

Description of wave propagation

This section addresses the question how to describe the propagation of waves if the velocity surface (contour in the case of surface waves) is known. Let us first recall the classification of surfaces related to the wave propagation. Then, the applicability of Fermat's minimum time principle and the direct usage of eikonal equations will be discussed.

Characteristic surfaces

In acoustics, three characteristic surfaces are used to characterize the wave propagation (see [START_REF] Auld | Acoustic Fields and Waves in Solids[END_REF] and [START_REF] Royer | Elastic Waves in Solids I: Free and Guided Propagation[END_REF]).

Wave surface. The wave surface (or the group velocity surface) describes the propagation of acoustic energy. This surface is involved in the formulation of minimum time principles, e.g. Fermat's law. The wave surface is the locus of points traced by the energy velocity vector V e , drawn from a fixed point O, as the propagation direction varies. The propagation direction n (the normalized wave vector) is orthogonal to the wave surface (see Fig. 5a). It should be noticed that V e is as a rule not collinear to the wave vector in the case of anisotropic media.

Phase velocity surface. The phase velocity surface (see Figures 1a and4a) describes the propagation of wave fronts. It defines the Hamiltonian of the eikonal equation. The phase velocity surface is obtained from the wave surface by projecting the vector V e onto the wave propagation direction n (see Fig. 5b) so that the phase velocity vector V is given by V = (V e • n) n. It should be noticed that the phase velocity surface can be constructed independently on V e , e.g. as shown in Section 2 and Subsection 3.3, and the wave surface can then be defined through the phase velocity surface.

Slowness surface. The slowness surface (see Figures 1b and4b) indicates the local convexity/concavity properties of the Hamiltonian of the eikonal equation.

The slowness surface is related to the phase velocity surface by the inversion through the origin (see Fig. 5c), i.e. m = n/|V |. The energy velocity, V e , is normal to the slowness surface at all points. Local concavities on the slowness surface can cause formation of cusps ("swallow tails") on the wave surface as it is shown in Fig. 5c: The arc (acb) is mapped into a "swallow tail" on the wave surface. This points out to the intersection of characteristics of the eikonal equation.

Fermat's principle

The Fermat principle describes how a ray, trajectory orthogonal to the wave front at all time instants, propagates from point A to point B. The principle says that the propagation time should be minimal. To express this, consider the minimization problem Thus, feasible rays are solutions of the Euler equation

T = B A dt = B A ds |V e | → min
L x - d dτ L ẋ = 0 .
This approach works well if the wave velocity V e (x, n) is well-defined for all directions n. It holds if the slowness surface is convex, which is as a rule violated in the case of anisotropic media. The next subsection discusses the method of direct solving eikonal equations.

Eikonal and Hamilton-Jacobi equations

Let c(x, n) be the phase velocity depending on the spatial position x and the propagation direction n, |n| = 1, Let S(t, x) be the phase function that shows the phase of the wave at the time instant t and at the point x. According to the results of Section 2, cf. equation ( 6), the eikonal equation reads

S t -|∇S|c x, ∇S |∇S| = 0. (14)
If the Hamiltonian c(x, p/|p|)|p| is convex in p, the method of characteristics can be used for solving equation ( 14) in the case of convex wave emitter. Under these conditions, the characteristics representing the rays do not intersect each other. If the convexity property is violated, equation ( 14) may not have classical solutions. Nevertheless, it is always uniquely solvable in the sense of viscous solutions (see [START_REF] Crandall | Viscosity Solutions of Hamilton-Jacobi Equations[END_REF] and [START_REF] Subbotin | Generalized Solutions of First Order PDEs: The Dynamical Optimization Perspective[END_REF]). Moreover, a unique viscous solution of ( 14) is the valid phase function. Thus, we arrive at the idea to use numerical methods of finding viscosity solutions of Hamilton-Jacobi equations. It should be noticed that common Lax-Friedrichs methods (see e.g. [START_REF] Mitchell | Application of Level Set Methods to Control and Reachability Problems in Continuous and Hybrid Systems[END_REF]) are not applicable in this case because they smooth solutions very strong. On the other hand, the authors have developed numerical methods that do not contain any smoothing (see e.g. [START_REF] Botkin | Stable Numerical Schemes for Solving Hamilton-Jacobi-Bellman-Isaacs Equations[END_REF] and [START_REF] Botkin | Approximation Schemes for Solving Disturbed Control Problems with Non-Terminal Time and State Constraints[END_REF]). These methods assume that the Hamilton-Jacobi equations arise from conflict control problems (differential games). Therefore, their application requires solving the following problem: Given an eikonal equation, it is required to construct a differential game whose Hamiltonian coincides (up to the sign) with that of the eikonal equation. The next section shows how to do that using an unification technique proposed in [START_REF] Krasovskii | On the Problem of the Unification of Differential Games[END_REF].

Usage of differential games

Assume that a set M ⊂ R d , d = 2 or 3, represents the shape of the acoustic wave emitter. For example, M is a ball in the case of bulk crystal, and M is the two dimensional area of the input electrodes in the case of biosensor. Let the game dynamics be described by the following system of ordinary differential equations:

ẋ = f (x, u, v), x ∈ R d , t ∈ (-∞, 0], u ∈ P ⊂ R a , v ∈ Q ⊂ R b , (15)
where u and v are control parameters of the first and second players, respectively. Introduce the signed distance, σ, to the set M as follows:

σ(x) = dist(x, M ) if x ∈ M , and σ(x) = -dist(x, R d \M ) if x ∈ M .
Consider the objective functional, γ, defined on the trajectories of ( 15) as follows:

γ(x(•)) = min τ ∈[t,0] σ(x(τ )). ( 16
)
The game is formalized using the concept of feedback strategies (see [START_REF] Krasovskii | Game-Theoretical Control Problems[END_REF]). The value function is defined by the relation

Ψ (t, x) = max V min x(•)∈X(t,x,V) γ(x(•)),
where V is a feedback strategy of the second player, and the set X(t, x, V) expresses the actions of the first player. This set consists of all limits of Euler trajectories of [START_REF] Krasovskii | On the Problem of the Unification of Differential Games[END_REF] which are obtained when the second player chooses v ≡ V(t i , x(t i )) on each interval [t i , t i+1 ) of partitions of [t, 0], and the first player uses admissible controls u(ξ), ξ ∈ [t, 0]. In doing that, all possible partitions whose diameter tends to zero and all admissible controls of the first player are exhausted. All Euler trajectories start at t from the initial state x.

The value function is locally bounded and Lipschitzian (see e.g. [START_REF] Subbotin | Generalized Solutions of First Order PDEs: The Dynamical Optimization Perspective[END_REF]).

Define the Hamiltonian

H(x, p) = max v∈Q min u∈P p, f (x, u, v) , p ∈ R d , (17)
and consider the Hamilton-Jacobi-Bellman-Isaacs equation

Ψ t + H(x, Ψ x ) = 0, Ψ (0, x) = σ(x). (18)
It is proven in [START_REF] Botkin | Stable Numerical Schemes for Solving Hamilton-Jacobi-Bellman-Isaacs Equations[END_REF] that the value function of the game [START_REF] Krasovskii | On the Problem of the Unification of Differential Games[END_REF] with the objective functional ( 16) is a viscosity solution of (18). Therefore, the following proposition holds: This proposition opens the way to use numerical methods of finding viscosity solutions of Hamilton-Jacobi equations. The only question consists in constructing an appropriate differential game whose Hamiltonian satisfies the condition of Proposition 1. The next subsection discusses this task.

Unification

Denote E(x, p) = -c x, p/|p| |p|, p ∈ R d . To find a differential game whose Hamiltonian coincides with E, the technique of unification (see [START_REF] Krasovskii | On the Problem of the Unification of Differential Games[END_REF]) can be used. Consider the following conflict control system

ẋ = E(x, p)p + q, x, p, q ∈ R d , |p| = 1, |q| = λ, p, q ≥ 0. (19)
Here, q is the control parameter of the first player who strives to minimize the objective functional ( 16), whereas p is the control parameter of the second player who maximizes the objective functional. The parameter λ is a constant which is greater than the Lipschitz constant of the function E in p. It should be noticed that the proof of this proposition essentially uses the positive homogeneity of the function p → E(•, p). Therefore, the unification procedure cannot be applied in the case of absence of positive homogeneity.

Proposition 2. If |E(x, p 1 ) -E(x, p 2 )| < λ|p 1 -p 2 |, p 1 , p 2 ∈ R d , |p 1 | = 1, |p 2 | = 1,
Assume now that |p| = 1, then E(x, p)p = -c(x, p)p. Taking into account that the points c(x, p)p exhaust the phase velocity surface V surf (x) when varying p, the game (19) can be rewritten as ẋ = -p + q, p ∈ V surf (x), |q| = λ, p, q ≥ 0.

Assuming that only the velocity magnitude depends on the spatial position yields the game ẋ = -a(x)p + q, p ∈ V surf , |q| = λ, p, q ≥ 0,

where the phase velocity surface V surf does not depend on x. Notice that the coefficient a(x) is necessary to take into account the damping area of the biosensor, the phase velocity is strongly reduced there. Therefore, a(x) ≡ 1 outside the damping area, and a(x) → 0 towards the outer boundary. Numerical methods developed by the authors (see [START_REF] Botkin | Stable Numerical Schemes for Solving Hamilton-Jacobi-Bellman-Isaacs Equations[END_REF] and [START_REF] Botkin | Approximation Schemes for Solving Disturbed Control Problems with Non-Terminal Time and State Constraints[END_REF]) provide an effective tool for finding the value function of the game (20) with the objective functional [START_REF] Krasovskii | Game-Theoretical Control Problems[END_REF]. The next section uses these methods and demonstrates a good reconstruction of wave patterns.

6 Simulation of wave propagation using the unified differential game

Now, the differential game (20) with the objective functional ( 16) is used for the computation of wave fronts. First, consider waves in LiT aO 3 piezoelectric crystals. Assume that the wave emitter, the set M , is a ball of radius 0.1, and a(x) ≡ 1, i.e. there is no damping area. The phase velocity surface is shown in Figure 1a. Thus, all data required for the formulation of the differential game ( 16) and (20) are available. Application of a finite difference upwind scheme described in [START_REF] Botkin | Stable Numerical Schemes for Solving Hamilton-Jacobi-Bellman-Isaacs Equations[END_REF] yields an approximation, Ψ , of the value function of the game.

According to Proposition 1, the wave front at any time instant t ≥ 0 is given by the relation {x ∈ R 3 : Ψ (-t, x) = 0}. Figure 6 shows the wave front at t = 0.25, 1, and 2 ms. Consider now wave fronts for surface shear waves in the biosensor structure (see Section 3). The phase velocity contour is shown in Figure 4a. Figure 7 shows the position of the wave front for a time sequence with a small sample time. The wave emitter and the damping area are easily recognizable in this figure.

Figures 8 and9 show the case of two wave emitters. Figure 10 shows the wave propagation in the presence of a hole. The hole is interpreted as an obstacle such that the trajectories of the game (20) can not penetrate therein. This case is numerically processed using a method for finding value functions in games with state constraints (see [START_REF] Botkin | Approximation Schemes for Solving Disturbed Control Problems with Non-Terminal Time and State Constraints[END_REF]).

Conclusion

The technique presented in this paper is also appropriate for the numerical treatment of arbitrary Hamilton-Jacobi equations with positive homogeneous Hamiltonians. As it was seen, the unification procedure described in Section 5.1 does not use any specific features of the Hamiltonian with the exception of the positive homogeneity which is necessary for the proof of Proposition 2. This opens new possibilities of investigation of physical processes related to optimality principles involving non-convex Lagrangians and Hamiltonians.

Fig. 1 .

 1 Fig. 1. Characteristic surfaces for lithium tantalate (LiT aO3): (a) Phase velocity surface; (b) Slowness surface. The coordinate frames show the crystalline axes X, Y , and Z.

Fig. 2 .

 2 Fig. 2. Structure of acoustic biosensor: (a) Spatial representation; (b) Vertical crosssection with thicknesses of the layers.

Fig. 3 .

 3 Fig. 3. Piezoelectric α-quartz crystal: (a) Direction of crystalline axes; (b) Orientation of ST-cut.

Fig. 4 .

 4 Fig. 4. Characteristic contours for surface acoustic waves: (a) Phase velocity contour; (b) Slowness contour.

Fig. 5 .

 5 Fig. 5. Schematic explanation to characteristic surfaces: (a) Wave surface; (b) Phase velocity surface; (c) Slowness surface.

Proposition 1 .

 1 Let c x, p/|p| |p| be the Hamiltonian of the eikonal equation. If H(x, p) = -c x, p/|p| |p|, p ∈ R d , then the wave front at any time instant t ≥ 0 is given by the relation {x : Ψ (-t, x) = 0}.

Fig. 6 .

 6 Fig. 6. Wave fronts in lithium tantalate LiT aO3 piezoelectric crystal (quasi shear wave): (a) t = 0.25 ms; (b) t = 1 ms; (c) t = 2 ms.

Fig. 7 .

 7 Fig. 7. Surface wave fronts in the biosensor structure. The wave emitter is an unclosed ring.

Fig. 8 .

 8 Fig. 8. Surface wave fronts in the biosensor structure with two emitters: a half-ring and a ring.

Fig. 9 .

 9 Fig. 9. Surface wave fronts in the biosensor structure with two half-rings as the emitters.

Fig. 10 .

 10 Fig. 10. Surface wave fronts in the biosensor structure in the presence of a round hole. The wave emitter is a half-ring.