A new equation and exact solutions describing focal fields in media with modular nonlinearity

Brand-new equations which can be regarded as modifications of Khokhlov–Zabolotskaya–Kuznetsov or Ostrovsky–Vakhnenko equations are suggested. These equations are quite general in that they describe the nonlinear wave dynamics in media with modular nonlinearity. Such media exist among composites, meta-materials, inhomogeneous and multiphase systems. These new models are interesting because of two reasons: (1) the equations admit exact analytic solutions and (2) the solutions describe real physical phenomena. The equations model nonlinear focusing of wave beams. It is shown that inside the focal zone a stationary waveform exists. Steady-state profiles are constructed by the matching of functions describing the positive and negative branches of exact solutions of an equation of Klein–Gordon type. Such profiles have been observed many times during experiments and numerical studies. The non-stationary waves can contain singularities of two types: discontinuity of the wave and of its derivative. These singularities are eliminated by introducing dissipative terms into the equations—thereby increasing their order.


Introduction
Wave beams in quadratic nonlinear media are governed by the Khokhlov-Zabolotskaya (KZ) equation, published first in Ref. [1]: Here p is the acoustic pressure, x is the axial coordinate along which the wave propagates, r is the radial polar coordinate of the beam, τ = t − x/c is the time in the reference system moving with sound velocity c. The nonlinear parameter of the medium is ε. The history of construction of this equation and main physical results following from it are reviewed in Ref. [2]. The technique of derivation is based on a combination of the slowly varying profile method and the quasi-optical approximation and is described in detail in books [3,4]. The model (1) is widely used for underwater [3] and medical [4] applications. Modern ultrasonic technologies and medical devices are based mainly on high-intensity focused ultrasound (HIFU). Focusing is a natural way to concentrate energy, but also leads to strong nonlinear phenomena. When converging at the focal point, shock fronts form in the wave and a large nonlinear absorption appears. At the same time, the diffraction effects in the focal region are important. These diversified phenomena accompanying the nonlinear focusing in HIFU are described in the review [5] (Sect. 5) where a detailed list of references is given. Only those which contain results needed for our further analysis will be discussed here.
In Ref. [6], the coefficient of concentration of energy is calculated. It usually decreases with increasing intensity because of nonlinear decay of the wave (see the review by Naugol'nykh [7]). However, nonlinearity is sometimes able to increase the concentration due to the sharper focusing of higher harmonics [8,9]. Ostrovskii and Sutin analyzed this phenomenon [9] on the base of a stage-by-stage approach, which goes as follows.
Two stages can be distinguished at sharp focusing (see Fig. 1), at two ranges of the x-coordinate. In the first stage, which extends from the surface of concave source to the beginning of focal area, the wave is spherically converging. It is subjected to only nonlinear distortion-diffraction is negligible. At the second stage around the focal point within the short waist of length l * , the wave front becomes plane, and the wave profile distorts almost exclusively due to the phase shifts between harmonics, which appear due to the lowfrequency dispersion caused by diffraction.
As is shown in Fig. 1, the interesting waist region has a form which is close to cylindrical. For beams round in cross section, the length of this cylinder is l * , and its base radius is a * [3,4]: In this formula l d = ωa 2 /2c is the diffraction length, R is the radius of curvature of the source (or focal distance) and a is the initial value of the radius of the beam.
Hereafter the focusing will be considered strong, and The cylindrical volume is the focal ray tube in which the nonlinear beam with plane fronts propagates. Based on the ideas expressed above, the following mathematical model was derived in Ref. [10] for a wave in a tube: Equation (2) is valid within the waist. One can derive it straight from the KZ (1) by formally putting the parabolic radial dependence of the field near the axis (r → 0) to: Substituting (3) into (1) and limiting the analysis to the paraxial region (r → 0), one can derive the 1D model (2). This model dramatically simplifies the qualitative analysis of nonlinear low-frequency diffraction. Evidently, the 1D Eq. (2) is much simpler to solve than the 2D Eq. (1) both analytically and numerically. Equation (2) obtained in Ref. [10] for the focal wave field has a quite general meaning and is a popular model of nonlinear dynamics. It is often called the Ostrovsky-Vakhnenko (OV) Eq. [11] and was initially constructed for other physical systems, like for internal waves in rotating ocean [12], and for finding solutions to mathematical physics models similar to solitons [13]. An analysis of an even more general OV equation with varying coefficients was done in Ref. [14].
This current article is associated with Ref. [10], but here appears fundamentally different mathematics. The new equation that is solved here describes a medium with a modular nonlinearity, as distinct from the more common quadratic nonlinear medium.
The form of the wave in the focal area was calculated in Ref. [10] for a common medium. At high intensities of the focused wave, the profile is depicted by a periodic sequence of arcuate sections having singularities of the derivative at the points of matching-typical for waves in nonlinear systems with low-frequency dispersion. Shockwaves containing steep fronts vary their shape while passing the focal area and they cannot be stationary. The maximum possible focal fields were estimated to be 50 kW/cm 2 [10], which is in agreement with both experiments and computer simulations.
In this paper, strongly distorted nonlinear forms of waves are described. These may in some aspects remind of solitons which in conservative systems result from competition between nonlinearity and dispersion [15,16]. Furthermore, solitons are single pulses, while this article treats periodic waves.

The modified KZ-OV equation for a bimodular material
Studies of plane waves in media with modular nonlinearity described by modified Burgers-type equations started recently [17]. These media have different elastic constants for the tensile and compressional deformations. They are bimodular materials, and examples are armed polymers and concretes (see [18], Ch.1, Sec.10), as well as cracked solids, and gas-liquid and multiphase media. The difference in elasticity at change of sign of deformation can be 10-15% [18]. The corresponding equations are formally linear when the wave function conserves its sign, i.e., for purely p > 0 or for purely p < 0, and the nonlinear phenomena are displayed only for waves containing both compression and rarefaction areas. The analogue to Eq. (2) for the bimodular medium can be obtained by the replacement of the quadratic nonlinearity p 2 /2 in (2) by the module nonlinearity | p|: The scheme of derivation of (4) is standard [3,4]. For example, for a fluid the equations of motion (Euler) and continuity for weak disturbances of pressure and density can be reduced to the following differential equation: The equation of state is here: The density variation ρ is gotten rid of by combining (5) and (6). Thereafter the model is transformed to a coordinate system [x, τ = t − x/c], moving along x with the sound velocity c. The wave distortion is slow along the x axis (∼μ x) and in the transverse direction r (∼ √ μ r ). This slow variation is caused by weak influence of nonlinear and diffraction effects and these terms are small on the order of μ 1. Neglecting terms of second and higher orders of the small parameter μ, and restricting attention to the focal region, one obtains the evolution Eq. (4). By analogy, this equation can be derived for solid media from the equations of the theory of elasticity.
It is convenient to reduce (4) to dimensionless form. Use the notations: The meanings of the constants ω and p 0 are the typical frequency and amplitude of initial signal, and l NL is the nonlinear length. The number γ 2 is the similarity criterium. At γ 2 1, diffraction dominates, and at γ 2 1, nonlinearity is stronger. The notations (7) reduce (4) to the following equation:

Stationary solution to the modified KZ-OV
Equation (8) can be represented by a pair equivalent to Klein-Gordon's equations: The upper signs are for V > 0, and the lower signs are for V < 0. The linear equations (9) can be easily solved separately, but minor complications appear when linking the positive and negative branches of the solution to Eq. (9). The positive branch is described by one of the two formulas: For the negative branch, by analogy: where δ, C and D are constants.
In (10) and (11), one can see that the waves V + and V − can propagate with equal velocities. For example, equating the phases from the second formula in (10) and the first formula in (11) one obtains: It follows from (12) that: Consequently, for any given number γ 2 the constants δ 2 + and δ 2 − can be determined to fullfil relation (13). This means that the branches V + and V − being linked at an initial z = z 0 will be continuously sewn together during the propagation of the wave, i.e., at all subsequent z.
To calculate the form of stationary wave, it is convenient to seek for the steady-state solutions to Eq. (8) in the form: For beams having their maximum amplitude on the axis, diffraction increases the velocity of propagation, so it must be that β > 0 [19]. The substitution (14) reduces (8) to the ordinary differential equation: Its two solutions, which can be linked through the use of the formulas (10)- (13), can be written as: These solutions are valid for 0 < β < 1.
For the convenience to perform the matching procedure, we will combine the solutions (16) to construct an even function. One can do that by shifting along T because we are solving the autonomous equation whose coefficients have no T -dependence. The equation is nonlinear, but is homogeneous in the variable V , which means that all integration constants in (15) can be multiplied by the same number. Therefore, the solution (16) can be rewritten in a more simple form: Both branches of the solution (17) are continuously matched in the point T 1 . In addition, continuity of derivative in the same point is required and also that the mean-over-period value of the wave field is zero: The formula (18) follows directly from Eq. (8). The physical meaning of (18) is the conservation of linear momentum of a volume of the medium in which the wave propagates [4]. These additional conditions make it possible to calculate the constant D in formula (17), as well as the connection between parameters γ and β. The final form of this solution is: The solution (19) for one period of wave is shown in Fig. 2. Curves 1, 2, 3 correspond to different values of the parameter β = 0, 0.8, 0.99. One can see that the duration of the positive part is shorter than the negative. Moreover, the branch V + has a sharp maximum, and the derivative has a discontinuity at T = π . This singularity would evidently be smoothed by the presence of dissipation. The negative branch V − is less steep but has longer duration. The maximum of V + exceeds the maximum of module of the negative branch |V − |: This difference increases with increase in velocity of propagation, i.e., when β → 1. Wave profiles such as those shown in Fig. 2 have been observed many times during experiments and numerical studies of the KZ equation near the focus and have been described by approximate analytical methods as well. However, the exact analytical solution for modular nonlinearity is found only in this present work. The mean intensity of the wave (19) can be defined through the integral over the period of the square of the function: Another periodic solution can be constructed by matching the trigonometric functions for both positive V + and negative V − branches. The joining of such functions is possible at higher values of velocities β > 1. In this case the analogue of (17) is: This solution is continuous in the point T 1 . By means of the additional condition (18), one of the constants is calculated: The connection between constants β and γ follows from the requirement for the period to be 2π : With account for the connections (23), (24), the final form of the solution (22) is: By analogy with the main characteristics of the solution (19), namely its maximum value (20) and its mean intensity (21), the corresponding values for the solution (25) are: In this case too, the maximum of V + exceeds the maximum of the absolute value of the negative branch |V − |. This difference increases with decrease in velocity of propagation, i.e., when β → 1 (see Fig. 3).

Non-stationary solution
The case to consider now is when the wave arriving at the input of the waist does not satisfy the conditions of linking, periodicity and zero mean level (18). This wave can then not be stationary, i.e., it cannot propagate with an unchanged shape. As an example, consider an harmonic initial wave and start by using the formulas (10) and (11): The only requirements for this wave is to be continuous and monochromatic. Put the constants C + 2 = C − 2 = C and δ + = δ − = 1. Then inside one period, the initial wave (at z = 0) is: During propagation, at increasing z, the wave is described by the solution: The initial wave (28) is shown in Fig. 4 by the dashed curve, and the solution (29) by the solid curve. The diffraction is not too strong-as γ 2 < 1. As shown in Fig. 4, nonlinearity creates three singularities during each period: the break of the function V (θ ) itself (at θ = −γ 2 z) and the breaks of its derivative (at θ = π − (1 + γ 2 )z and at θ = π + (1 − γ 2 )z). These singularities are marked by solid circles.
A discontinuity singularity occurs in the case when V + follows V − . A shock front appears inside the intersection area and cuts two equally shaded areas (Fig. 4).
Another type of singularity, a discontinuity of the derivative occurs if V − follows V + . In this case the two parts move away from each other. One can see in Fig. 4 how a nonlinear decay of the wave takes place, and at the finite distance z = π the wave completely disappears.
The Fourier series expansion of solution (29) is: The harmonic amplitudes B n are for 0 < z < π given by π 2 B n = n n 2 − 1 sin z + (−1) n n 2 − 1 sin(nz), With increase in z, the fundamental harmonic (n = 1) fades monotonically, giving its energy to the higher harmonics, which are in turn transferred to the medium-a nonlinear dissipation process. At first, the amplitudes of the higher harmonics (n > 1) grow while receiving energy from the fundamental wave, but thereafter they decay as well-to zero at z = π .

Smoothing of singularities by dissipation
The singularities shown in Figs. 2, 3 and 4 would be smoothed by dissipation, which of course exists to some extent in any medium. To derive the smooth solutions, it is necessary to generalize the differential equations (4) and (8) by introducing a dissipative term of higher order. Starting with Navier-Stokes' equation of motion, one can, instead of (5), derive the following differential relation: Here b is the dissipative coefficient, and ρ 0 is the equilibrium density of the medium. Using the equation of state (6) and repeating the procedure of derivation described in Sect. 2, we obtain the generalized Eq. (4): It is reasonable to call this equation a modified Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation [3] because it differs from original KZK in the replacement of the quadratic nonlinearity p 2 /2 by a modular one | p|.
Using the dimensionless variables (7), Eq. (33) can be reduced to a normalized equation [generalizing (8)]: The new parameter Γ is connected with the typical dissipative length l DISS for a given medium through: To smooth the discontinuity of a function (see shock front in Fig. 4), the singularly perturbed problem must be considered with Γ being a small parameter in the highest (third) derivative in Eq. (34). As follows from Fig. 4, the symmetric shock is located at θ = −γ 2 z. If a new temporal variable T = θ + γ 2 z is introduced, the shock will be immovable at T = 0. With the new variable ξ = T /Γ , Eq. (34) takes a steady-state form: One can see that at weak dissipation Γ 1 the shock front structure is governed by the simplified equation: The solution to (35) satisfying the two boundary conditions at the infinity-V → V 0 , ξ → +∞ and V → −V 0 , ξ → −∞-has the form: It describes the steep shock front which transforms to an ideal discontinuity at Γ → 0. In accordance with the matched asymptotic expansions method described by Nayfeh [20], the result (38) is the main term of the inner expansion, and corresponding term of the outer expansion is determined by (29): A similar problem can be solved to smooth the derivative discontinuity. The final result, omitting the calculating details, is given as: Here Φ(x) = 2 √ π x 0 exp(−ν 2 ) dν is the error integral, and θ 1 (z) = π − (1 + γ 2 )z and θ 2 (z) = π + (1 − γ 2 )z are the positions of the singular points (see Figs. 4,5). With increase in the number Γ and the distance z, the smoothing is enhanced (see curves for Γ z = 0.25·10 −2 and 10 −2 in Fig. 5).
Note that the analytical smoothing described above is consistent with the results of numerical modeling performed by Nazarov and co-workers for other equations with modular nonlinearity [17,21]. The construction of smoothed solutions can be done using the standard method of matched asymptotic expansions. A detailed summary of these methods is given in the book by Nayfeh [20].

Conclusion
This article contains two types of results. On the one hand, new nonlinear equations are obtained which generalize the previously studied models of KZ, OV and KZK types. Solutions to these equations are found, describing the nonlinear dynamics of waves in focal areas. In media without dissipation the singularities that form are of two types: discontinuities of the wave itself and discontinuities of its derivative. Smoothing procedures which eliminate these singularities are thereafter indicated. This group of results is attributed to the field of mathematical physics and nonlinear dynamics.
On the other hand, the studies of the physics of nonlinear wave focusing, which was initiated in Ref. [10], are here continued for high-intensity waves in bimodular media. These materials exist in nature-one example is cracked or granular materials [22]. Modular and other related unusual nonlinearities exist also in composites and meta-materials, increasing the applicative interest of these studies [23], and the mathematical and physical aspects of the nonlinear dynamics in such materials will be expanded in the nearest future.
The analytical methods used in this work are standard ones, and the main difficulty encountered was that singularities exist as long as the models do not have damping. These singularities disappeared when dissipation was introduced. As any real system of this type always has a nonzero amount of dissipation, this also made the results being closer to the physical reality. The results here are new as they investigate the modular material model which is based on real materials' structures and the obtained results are qualitatively close to their behavior.
The biggest advantage with the modular equation investigated here is that exact analytical solutions were obtained. There is therefore no need for numerical results to make comparisons to and all solutions in this paper are purely analytical.