text:
# Renormalization of the elementary excitations in hole- and electron-doped cuprates due to spin fluctuations

## I Introduction

Elementary excitations in the cuprates are of central interest in order to learn more about the correlations and the pairing mechanism for superconductivity. For example, it is well-known that the understanding of the elementary excitations in conventional superconductors like lead as measured by tunneling spectroscopy played the crucial role in accepting the picture of electron-phonon-mediated Cooper-pairing . In the high-$`T_c`$ cuprates one expects that due to the presence of antiferromagnetic spin fluctuations a strong renormalization of the spectral density and the corresponding energy dispersion may occur. It was shown by several groups  that the so-called ’dip-hump’ structure seen in tunneling experiments  can be explained in terms of the feedback effect of superconductivity arising from the structure in the gap function $`\mathrm{\Delta }(𝐤,\omega )`$. Moreover, it has been argued that this structure reflects the effective pairing interaction and points towards a spin-fluctuation-mediated Cooper-pairing. Recent developments in angle-resolved photoemission spectroscopy (ARPES) allow to study the elementary excitations directly and in much more detail. In particular, the anisotropy of the elementary excitations close to the Fermi energy in the different parts of the Brillouin Zone (BZ) are studied. This is important, since the coupling of the quasiparticles to spin fluctuations varies at different parts of the BZ. Furthermore, the understanding of the structures seen by ARPES and their doping, momentum, and temperature dependence will help to understand the role played by spin fluctuations in contrast to phonons regarding the formation of superconductivity in the high-$`T_c`$ cuprates.

1. Theoretical background: Similar to the electron-phonon case, the coupling between the quasiparticles and the spin excitations should influence characteristically the energies $`\omega _𝐤`$ of the hole- or electron carriers

$$\omega _𝐤=ϵ_𝐤+\text{Re }\mathrm{\Sigma }(𝐤,\omega ).$$
(1)
Here, $`ϵ_𝐤`$ refers to the bare energy dispersion of the quasiparticles assuming no interaction with the spin excitations or phonons, i.e. a tight-binding energy dispersion for the CuO2-plane. In general, the self-energy $`\mathrm{\Sigma }(𝐤,\omega )`$ results from the coupling of the particles to spin excitations, see Fig. 1 for an illustration.

Obviously, the self-energy is a functional of the spin susceptibility, $`\mathrm{\Sigma }=\mathrm{\Sigma }\{\chi \}`$. In its simplest form, the latter is of Ornstein-Zernicke form  that allows for a sharp enhancement of fluctuations near the antiferromagnetic wave vector Q$`=(\pi ,\pi )`$:

$$\chi (𝐪,\omega )=\frac{\chi _𝐐}{1+\xi ^2(𝐪𝐐)^2i\frac{\omega }{\omega _{sf}}}.$$
(2)
Here, $`\chi _𝐐`$ is the value of the static spin susceptibility at the wave vector Q, $`\xi `$ is the magnetic correlation length, and $`\omega _{sf}`$ is the characteristic frequency of spin fluctuations (roughly the peak position in the imaginary part of Eq. (2)). Due to the fact that $`\chi (𝐪,\omega )`$ has only structures around $`\omega _{sf}`$ and wave vector Q, their influence on the elementary excitations is expected to be very anisotropic at different parts of the Brillouin Zone.

In order to illustrate the anisotropy of the elementary excitations we show in Fig. 2 the calculated Fermi surface for hole-doped cuprates . For the calculation we use the tight-binding energy dispersion for the CuO2-plane

$$ϵ_𝐤=2t\left(\mathrm{cos}k_x+\mathrm{cos}k_y\right)+4t^{}\mathrm{cos}k_x\mathrm{cos}k_y\mu .$$
(3)
Here, $`t`$ and $`t^{}`$ refer to the hopping of a hole (electron) between nearest, next-nearest sites on the square lattice, and $`\mu `$ is the chemical potential that defines the doping. In Eq. (3) and in the following we set the lattice constant to unity. One clearly sees that the scattering of quasiparticles by spin fluctuations in cuprates is anisotropic. First, the antiferromagnetic wave vector $`𝐐=(\pi ,\pi )`$ connects exactly quasiparticles at the Fermi level close to the ($`\pi ,0`$) points of the first Brillouin Zone. These quasiparticles

experience the strongest coupling to antiferromagnetic spin fluctuations. Quasiparticles at the diagonals of the BZ are not connected by the $`𝐐`$ and thus have smaller scattering by spin fluctuations at the Fermi level. The corresponding point on the Fermi surface is called cold spot or nodal point. Furthermore, the quasiparticle states at the Fermi level close to the $`M`$ point in the first BZ are usually called antinodal points. Note that the quasiparticle states at cold and hot spots which are connected by $`𝐐`$ lie close to the Fermi level. This will be important later for discussing the kink feature. Therefore, we may define three characteristic regions at the Fermi surface regarding their sensitivity to coupling to antiferromagnetic spin fluctuations.

2. Experimental findings: Let us now discuss the experimental situation close to the Fermi level as measured by ARPES. The combined study of energy distribution curves (EDC) and momentum distribution curves (MDC) allows the study of the quasiparticle excitations close and well below the Fermi level up the energies of $`200`$ meV. Most importantly, recent ARPES studies by various groups  reveal a kink structure in the energy dispersion at energies about $`60\pm 15`$meV below the Fermi energy and along the nodal direction $`(0,0)(\pi ,\pi )`$ direction of the first BZ. Furthermore, this kink structure appears in the normal state and almost does not change if one enters the superconducting state.

We illustrate the formation of the kink in Fig. 3 where we show the unrenormalized tight-binding energy dispersion $`ϵ_𝐤`$ in the normal state along the route $`(0,0)(\pi ,0)(\pi ,\pi )(0,0)`$ of the first BZ.

Below Ef, the measured energy dispersion $`\omega _𝐤`$ along $`(0,0)(\pi ,\pi )`$ has changed compared to the bare tight-binding case due to the self-energy corrections as it follows from Eq. (1). We define the kink at the inflection point of the dashed curve where the renormalized dispersion tends to approach again the bare dispersion $`ϵ_𝐤`$. Recent studies by Dessau and co-workers  reveal no kink structure along $`(\pi ,0)(\pi ,\pi )`$ direction in the normal state, but only in the superconducting state.

3. Doping and $`𝐤`$-dependence of the kink: Previously we have shown that from the momentum and energy conservation one expects changes of the quasiparticle velocities $`v_𝐤=ϵ/𝐤`$ in the nodal direction for $`\omega \omega _{sf}`$ and $`𝐤𝐤^{}=𝐪𝐐`$. This leads to the kink in the $`(0,0)(\pi ,\pi )`$ direction that has been observed by angle-resolved photoemission spectroscopy (ARPES) . Here, we extend our previous analysis and discuss in more detail the doping dependence and the anisotropy of the kink feature as a function of ($`𝐤𝐤_F`$) and energy $`\omega _𝐤`$ in the normal and superconducting state, for different routes $`𝐤𝐤_F`$ in the Brillouin zone: $`(0,0)(\pi ,\pi )`$, $`(0,0)(0,\pi )`$, and $`(\pi ,0)(\pi ,\pi )`$. Due to the fact that the superconducting order parameter has nodes along $`(0,0)(\pi ,\pi )`$ and maxima for $`(0,0)(0,\pi )`$ one expects additional changes in the kink structure as observed recently . This should help to contrast renormalization due to spin fluctuations to what we get due to electron-phonon interaction. Furthermore, we investigate the interdependence of the elementary excitations with the resonance peak that is observed in hole-doped cuprates by inelastic neutron scattering (INS) below $`T_c`$. It follows also from Eq. (1) and Im$`\mathrm{\Sigma }(𝐤,\omega )\omega ^2\omega `$ that the temperature dependence of the kink feature should reflect characteristically the coupling to spin fluctuations and be different from the electron-phonon coupling case.

Note that in the case of electron-doped cuprates neither a resonance peak nor a kink feature is present. We will also show that the kink feature is not restricted to the cuprates, but is also expected for other novel superconductors where quasiparticles couple strongly to spin excitations. A particularly interesting case might be $`\text{Sr}_2\text{RuO}_4`$ with large anisotropic behavior of $`\chi _{zz}(𝐪,\omega )`$ and in-plane $`\chi _+(𝐪,\omega )`$ . Generally, a strong nesting behavior of the Fermi surface might yield pronounced kink features.

This paper is organized as follow: In Section II we present the theory, and in Section III we discuss our results for the kink structure, its doping- and temperature dependence for hole- and electron-doped cuprates. Finally, in Section IV we summarize our analysis and contrast renormalization due to spin fluctuations and phonons.

## II Theory

A. One-band Hubbard model: The elementary excitations and the spin excitations are key quantities determining the properties of the cuprates and other superconductors with strong correlations and magnetic activity. The quasiparticles, holes or electrons, interact strongly with spin fluctuations and also with phonons. However, phonons and spin fluctuations differ with respect to their doping dependence and anisotropy. This is clearly demonstrated by neutron scattering experiments, for example. Also the different behavior of hole- and electron-doped cuprates and the feedback of superconductivity on $`\chi (𝐪,\omega )`$ is important.

In this paper we employ an effective one-band Hamiltonian. This is justified because upon hole doping antiferromagnetism disappears due to Zhang-Rice singlet formation and quenching of Cu-spins. In this one-band picture the Coulomb interaction between the quasiparticles refers to an effective interaction (i.e. the Hubbard $`U`$) within the conduction band. Then, further doping increases the carrier mobility and a system of strongly correlated quasiparticles occurs. In the overdoped case less magnetic activity is present yielding usual Fermi liquid. We assume $`UW/2`$ ($`W=`$ bandwidth) independent of the doping concentration.

The main physics of a single CuO2-plane is the two-dimensional one-band Hubbard model given by

$$H=\underset{ij\sigma }{}t_{ij}\left(c_{i\sigma }^+c_{j\sigma }+c_{j\sigma }^+c_{i\sigma }\right)+U\underset{i}{}n_in_i,$$
(4)
where $`c_{i\sigma }^+`$ creates an electron with spin $`\sigma `$ on site $`i`$, $`U`$ denotes the on-site Coulomb interaction, and $`t_{ij}`$ is the hopping integral. After diagonalization of the first term, one arrives at the bare tight-binding energy dispersion given by Eq. (3). The description of the electron- and hole-doped cuprates within a one-band approximation is possible if one takes into account different parameters and quasiparticle dispersion. Note that in the case of electron doping the electrons occupy the copper $`d`$-band, while in the hole-doped case holes refer mainly to the oxygen $`p`$-states yielding different dispersion parameters. Furthermore, the energy dispersions for optimally hole-doped La2-xSrxCuO4 (LSCO) and electron-doped Nd2-xCexCuO4 (NCCO) behave differently around $`(\pi ,0)`$ point. While in the case of LSCO the flat band (leading to the van Hove singularity in the density of states) lies close to the Fermi level, in NCCO the flat band is approximately $`300`$ meV below the Fermi level. Then, using $`t=250`$meV and $`t^{}=0.1t`$, one describes the hole-doped LSCO dispersion, whereas $`t=138`$meV and $`t^{}=0.3t`$ are needed for the description of the electron-doped NCCO compound fitting earlier photoemission data. These parameters including an intermediate coupling $`U=4t`$, will be used for our calculations of various physical quantities in the normal and superconducting state of electron- and hole-doped cuprates .

B. Generalized Eliashberg equations: In this one-band model, we assume that the same electrons (holes) are participating in the formation of antiferromagnetic fluctuations and in Cooper-pairing due to the exchange of these fluctuations. Thus, both the magnetic susceptibility and the quasiparticle self-energy must be calculated self-consistently. This is possible in the FLEX approximation . In this approach the dressed one-electron Green’s function are used to calculate the charge and spin susceptibilities. These susceptibilities are then used to construct an effective Berk-Schrieffer-like  pairing interaction $`V_{\mathrm{eff}}`$ describing the exchange of charge and spin fluctuations. The generalized Eliashberg equations are derived in Appendix A. In order to demonstrate the significant role of $`V_{\mathrm{eff}}`$ in our work, we also show in Fig. 13 its corresponding diagramatic representation. Note, in general, if the Cooper-pairing and the effective pairing potential $`V_{\mathrm{eff}}`$ are generated by the same quasiparticles (solid lines in Fig. 13), strong self-energy and feedback effects on $`G(𝐤,\omega )`$ and $`\chi (𝐪,\mathrm{\Omega })`$ are expected .

To be more precise, we write down the quasiparticle self-energy components $`\mathrm{\Sigma }_\nu `$ ($`\nu =0,3,1`$) with respect to the Pauli matrices $`\tau _\nu `$ in the Nambu representation, i.e. $`\mathrm{\Sigma }_0=\omega (1Z)`$ (mass renormalization), $`\mathrm{\Sigma }_3=\xi `$ (energy shift), and $`\mathrm{\Sigma }_1=\varphi `$ (gap parameter). They are given by

$`\mathrm{\Sigma }_\nu (𝐤,\omega )`$ $`=`$ $`N^1{\displaystyle \underset{𝐤^{}}{}}{\displaystyle _0^{\mathrm{}}}𝑑\mathrm{\Omega }V_{\mathrm{eff}}(𝐤𝐤^{},\mathrm{\Omega })`$ (6)
$`\times {\displaystyle _{\mathrm{}}^+\mathrm{}}d\omega ^{}I(\omega ,\mathrm{\Omega },\omega ^{})A_\nu (𝐤^{},\omega ^{})`$
with

$$V_{\mathrm{eff}}=\left[P_s(𝐤𝐤^{},\mathrm{\Omega })(\delta _{\nu 1}\delta _{\nu 0}\delta _{\nu 3})P_c(𝐤𝐤^{},\mathrm{\Omega })\right].$$
(7)
This is a generalization of figure 1. $`P_s`$ and $`P_c`$ denote the spectral density of the spin and charge excitations, respectively, and are defined in Eqs. (A6) and (A7). The second part of Eq. (6) is given in Eqs. (A8) and (A9). It is interesting to remark that the above formulae remain valid even if the elementary excitations and the magnetic activity that controls $`V_{\mathrm{eff}}`$ would result from different quasiparticles.

The generalized Eliashberg equations allow us to calculate all the properties of the system self-consistently like superconducting phase diagram, elementary excitations, the superconducting order parameter, energy dispersion and dynamical spin susceptibility, for example . Note, in Eq. (6) self-energy effects due to phonons are neglected. Their contribution will be discussed in section IV and in Appendix B.

## III Results and Discussion

### A  Anisotropic Renormalization

1. Nodal direction: We start the discussion analyzing the spectral density of hole-doped superconductors in the normal state. The spectral density reveals the elementary excitations and in particular the renormalized energy dispersion. First, we present our results for the spectral density along the nodal $`(0,0)(\pi ,\pi )`$ direction in the first BZ.

In Fig. 4 we show the calculated spectral density $`N(𝐤,\omega )`$, i.e. the local density of states, as a function of frequency and momentum $`𝐤𝐤_F`$. The peak positions correspond to the renormalized energy dispersion. Due to coupling of holes to antiferromagnetic spin fluctuations the quasiparticle dispersion changes its slope and shows a pronounced kink feature at the energy $`\omega _{kink}75\pm 15`$meV.

How can one understand the kink feature in a simple way? At the first glance the occurrence of a kink in the nodal direction seems to be surprising, since the main interaction of the carriers with spin fluctuations occurs at the hot spots while the kink feature is present along the diagonal of the BZ close to the cold spots. This argument, however, considers only the quasiparticles exactly at the Fermi level. Away from the Fermi level but close to it (along $`(0,0)(\pi ,\pi )`$) the quasiparticles couple strongly to spin fluctuations, as can be seen from Fig. 2. Most importantly, as follows from Fig. 2, the largest scattering will occur at values of $`𝐤𝐤_F=𝐐`$ and $`\omega =\omega _{sf}`$. To be more precise, let us rewrite Eq. (6) as

$`\mathrm{\Sigma }(𝐤,i\omega _n)`$ $`=`$ $`T^2{\displaystyle \underset{\omega _m,\nu _m}{}}{\displaystyle \underset{𝐤^{},𝐪}{}}\stackrel{~}{\tau }_0G(𝐤𝐤^{},i\omega _ni\nu _m)\stackrel{~}{\tau }_0U^2`$ (8)
$`\times `$ $`{\displaystyle \frac{1}{2}}\text{Tr}\left[\stackrel{~}{\tau }_0G(𝐤+𝐪,i\omega _m+i\nu _m)\stackrel{~}{\tau }_0G(𝐪,i\omega _m)\right]`$ (9)
and approximate the Green’s function by its non-interacting part

$$G(𝐤,i\omega _n)G_0(𝐤,i\omega _m)=\frac{i\omega _n\stackrel{~}{\tau }_0+ϵ_𝐤\stackrel{~}{\tau }_3\varphi _𝐤\stackrel{~}{\tau }_1}{(i\omega _n)^2E_𝐤^2}.$$
(10)
with $`E_𝐤^2=ϵ_𝐤^2+\varphi _𝐤^2`$. Thus, after little algebra one obtains on the real axis

$`\mathrm{\Sigma }(𝐤,\omega ){\displaystyle \frac{U^2}{4}}`$ $`{\displaystyle \underset{𝐤^{}}{}}{\displaystyle _0^{\mathrm{}}}𝑑\omega ^{}{\displaystyle \frac{\text{Im}\chi _{\text{RPA}}(𝐤𝐤^{},\omega ^{})}{\omega \omega ^{}E_𝐤^{}}}`$ (11)
$`\times `$ $`\left[\mathrm{coth}\left({\displaystyle \frac{\omega ^{}}{2T}}\right)\mathrm{tanh}\left({\displaystyle \frac{\omega ^{}\omega }{2T}}\right)\right].`$ (12)
The imaginary part of the spin susceptibility may be obtained within the RPA and is approximately given by Ornstein-Zernicke expression (see Eq.(2)). This self-energy $`\mathrm{\Sigma }`$ now enters Eq. (1). It is important that the self-energy is mainly frequency-dependent, while the bare dispersion of the carriers is not. Already in the normal state, $`\mathrm{\Sigma }(𝐤,\omega )`$ has a maximum reflecting a corresponding maximum of $`\text{Im }\chi `$ at $`𝐪𝐐`$ and $`\omega ^{}\omega _{sf}`$. Note that $`\omega _{sf}`$ can be determined according to Moriya and Ueda and paramagnon theory from $`\chi _{RPA}=\chi _0/(1U\chi _0)`$ or equivalentlt from $`\chi _0^1(𝐪𝐐,\omega \omega _{sf})U=0`$ . Obviously, $`\omega _{sf}`$ is strongly doping dependent. This will be discussed later. Then, the kink position follows from the pole of the denominator of Eq. (12). This leads to the ’kink condition’

$$\omega _{kink}E_{𝐤𝐐}+\omega _{sf}(x).$$
(13)
This gives an estimate of the position of the kink. Furthermore, since the superconducting gap is zero for $`\omega =0`$, but not for $`\omega =\omega _{sf}`$, the kink feature along the nodal direction $`(0,0)(\pi ,\pi )`$ will change only slightly below $`T_c`$. This we have demonstrated previously .

2. $`(\mathrm{𝟎},\mathrm{𝟎})(\mathrm{𝟎},\pi )`$ direction: In order to see whether the kink feature is present in other directions of the Brillouin Zone, we show in Fig. 5 the evolution of the spectral density along the $`(0,0)(0,\pi )`$ direction. Despite of the fact that along this direction we do not cross the Fermi level, the kink feature is still present and is found at an energy similar to the one for the nodal $`(0,0)(\pi ,\pi )`$ direction. This indicates that the occurrence of the kink feature is not connected to some specific conditions which might be present only along the $`(0,0)(\pi ,\pi )`$ direction. Instead, the kink is characteristic for all direction where $`𝐤𝐤_F𝐐`$ and $`\omega \omega _{sf}`$. Also below $`T_c`$ we find that the kink feature is present in the $`(0,0)(0,\pi )`$ direction (not shown). Note, that our results are in fair agreement with experimental data .

3. Antinodal direction: In Fig. 6(a) we show our results for $`N(𝐤,\omega )`$ along the $`(\pi ,0)(\pi ,\pi )`$ route, i.e. the antinodal direction, of the first BZ in the normal state. Note, that the spectral density at the $`(0,\pi )`$ point is broader than at the antinodal point due to stronger coupling to spin excitations peaked at $`𝐪=𝐐=(\pi ,\pi )`$ as discussed in Fig. 2. Clearly, no kink is present. The absence of a kink structure can be explained with the flat structure of the $`\text{CuO}_2`$-plane around the $`M`$ point (see Fig. 3). Simply speaking, for a flat band the frequency dependence of $`\mathrm{\Sigma }`$ in Eq. (1) does not play a significant role and therefore no change of the velocity and no kink structure is present.

What does happen in the superconducting state? Below $`T_c`$ the superconducting gap $`\varphi (𝐤,\omega )`$ opens rapidly for decreasing temperature $`T`$ and becomes maximal in momentum space around the $`M`$ point reflecting the momentum dependence of the effective pairing interaction (see Eq. (6)). In addition, due to the frequency dependence of the gap the flat band around $`M`$ disappears.

In Fig. 6(b) we show results for $`N(𝐤,\omega )`$ at a temperature $`T=0.5T_c`$ where the superconducting gap has opened. A kink structure around $`\omega _{kink}50\pm 10`$ meV is present reflecting the magnitude of $`\varphi `$. Hence, in the $`(\pi ,0)(\pi ,\pi )`$ direction this kink feature is only present below $`T_c`$ and connected to the feedback effect of $`\varphi `$ on the elementary excitations. We will show later that this feedback is also important for the resonance peak seen in INS.

Note, that the superconducting gap $`\varphi (𝐤,\omega )`$ is calculated self-consistently in our theory and reflecting the underlying spin fluctuations which dominate the pairing potential $`V_{\mathrm{eff}}`$. Therefore, the occurrence of a kink structure only below $`T_c`$ in the antinodal direction is a direct fingerprint of the spin excitation spectrum . Furthermore, as we will discuss below, $`\text{Im }\chi (𝐐,\omega )`$ entering in Eq. (12) is peaked at the resonance frequency $`\omega _{res}`$ (roughly at $`\omega _{sf}+\mathrm{\Delta }`$). Therefore, the kink condition is given by

$$\omega _{kink}E_{𝐤𝐐}+\omega _{res}(x).$$
(14)
In Fig. 7(a)

the frequency dependence of $`\text{Re }\mathrm{\Sigma }(𝐤_a,\omega )`$ in the normal and superconducting state at the antinodal point $`𝐤=𝐤_a`$ is shown. Due to the occurrence of the resonance feature in $`\text{Im }\chi (𝐐,\omega )`$ and the related feedback of the superconducting gap $`\varphi (\omega )`$, $`\text{Re }\mathrm{\Sigma }`$ shows a pronounced structure below $`T_c`$ at energies of about $`\omega _{res}+\mathrm{\Delta }_0`$. Also the corresponding imaginary part, $`\text{Im }\mathrm{\Sigma }(𝐤=𝐤_a,\omega )`$, shows a peak below $`T_c`$ (see Fig. 7(b)). This pronounced behavior is responsible for the kink formation along $`(\pi ,0)(\pi ,\pi )`$ direction in the BZ. Therefore, while the kink features are present along $`(0,0)(\pi ,\pi )`$ and $`(\pi ,0)(\pi ,\pi )`$ directions in the superconducting state of hole-doped cuprates, their nature is qualitatively different. Along the nodal direction the superconducting gap is zero (for $`\omega =0`$) and thus the feedback effect of superconductivity on the elementary and spin excitations is small. Therefore, $`\omega _{sf}`$ determines mainly the formation of the kink feature. On the other hand, along the antinodal direction the gap is maximal and yields a strong feedback of superconductivity on $`\chi `$. Thus, in the superconducting state $`\omega _{res}`$ and $`\mathrm{\Delta }_0`$ yield the kink structure along $`(\pi ,0)(\pi ,\pi )`$ direction that is not present in the normal state.

### B Doping dependence of renormalization

The different reasons for the kink structures in hole-doped cuprates along different directions in the first BZ will be also reflected in their doping dependence. So far, the results we have shown were for optimal doping concentration $`x=0.15`$ that refers to a band filling of $`n=0.85`$. Note, the superconducting transition temperature $`T_c`$ behaves differently in the overdoped (OD) and underdoped (UD) regime:

$`T_c`$ $``$ $`\mathrm{\Delta }(T0),\text{OD}`$
$`T_c`$ $``$ $`n_s(T0),\text{UD}`$
where $`n_s`$ is the superfluid density that is calculated self-consistently from the generalized Eliashberg equations .

In the antinodal $`(0,\pi )(\pi ,\pi )`$ direction the kink is only present below $`T_c`$ due to the feedback of $`\varphi (\omega )`$. In the OD case, $`\varphi (\omega )`$ decreases reflecting a mean-field-like behavior. Thus, the energy where the kink occurs must decrease with overdoping:

$$\omega _{kink}(x)\mathrm{\Delta }_0(x).$$
(15)
This behavior is indeed observed by Dessau and co-workers . Note, the above argument remains true also in the strongly OD case where no resonance peak in $`\text{Im }\chi (𝐐,\omega )`$ occurs because the feedback effect of $`\varphi (\omega )`$ should always be present.

Regarding the kink along the nodal $`(0,0)(\pi ,\pi )`$ direction we note the following. In Fig. 8 we show the calculated doping dependence of Im $`\chi `$ at the antiferromagnetic wave vector $`𝐐`$ versus frequency in the normal state. One clearly sees the characteristic Ornstein-Zernicke behavior (see Eq. (2)) of $`\text{Im }\chi `$,

$$\text{Im }\chi (𝐪=𝐐,\omega )\frac{\omega \omega _{sf}}{\omega ^2+\omega _{sf}^2},$$
(16)
and that $`\omega _{sf}`$ increases with increasing doping from underdoped to overdoped cuprates. Since $`\omega _{sf}`$ determines the kink position along $`(0,0)(\pi ,\pi )`$ direction we expect

$$\omega _{kink}(x)\omega _{sf}(x)..$$
(17)
This is in qualitative agreement with experimental data  (for underdoped regime and optimally doped superconductors). On the other hand, the spectral weight of $`\text{Im }\chi (𝐐,\omega )`$ decreases drastically with overdoping. Therefore, the coupling of the quasiparticles to spin fluctuations is getting much weaker in the OD case. These two competing effects seem responsible for the non-monotonic and weak doping dependence of the kink position in the nodal direction .

In Fig. 9 we illustrate the kink feature resulting from the renormalization ($`\frac{d\mathrm{\Sigma }^{}}{d\omega }v^{}(1+\lambda )`$, $`vv_F`$) of the bare dispersion. We estimate $`2\lambda 3`$. This renormalization is doping dependent and stronger for underdoped hole-doped cuprates. Of course, we expect that the position of the kink as well as the change of the quasiparticle velocity ($`vv^{}`$) are important fingerprints of the coupling to spin fluctuations. Note, $`v^{}v`$ for frequencies $`\omega >\omega _{kink}`$ reflects mainly the width of the peak in $`\text{Im }\chi `$. Important is the slope ratio $`v^{}/v`$ for $`\omega <\omega _{kink}`$.

Another important behavior concerns the asymmetry between hole and electron-doped cuprates. Note that no kink feature has been reported in the electron-doped cuprates . It is believed that the electron-phonon coupling is much more pronounced in electron-doped cuprates than in hole-doped ones. This is indicated, for example, by the behavior of the resistivity $`\rho T^2`$ in the normal state at optimum doping and by the transition between $`d_{x^2y^2}`$-wave symmetry of the superconducting gap towards anisotropic $`s`$-wave as it has been observed in several experiments . Simply speaking, the spin fluctuations in electron-doped cuprates are weaker than in the hole-doped ones yielding a smaller $`T_c`$ and a smaller superconducting gap . Thus, no kink is present in the nodal direction and also no kink occurs in the $`(0,\pi )(\pi ,\pi )`$ direction below $`T_c`$. This is related to the fact that the flat band around $`(0,\pi )`$ lies in electron-doped cuprates well below the Fermi level and, therefore, it cannot be softened due to $`\varphi (\omega )`$.

### C Relation of kink and resonance peak

In Fig. 10 we show results for the spin susceptibility $`\text{Im }\chi _{RPA}(𝐐,\omega )`$ in the optimally electron(a)- and hole(b)-doped cuprates in the normal and superconducting state taking into account different tight-binding energy dispersions .

While in the normal state of hole-doped cuprates $`\omega _{sf}`$ is of order of $`25`$ meV, in the electron-doped ones its value is much larger ($`\omega _{sf}70`$ meV) and $`\text{Im }\chi `$ is much less pronounced. Therefore, antiferromagnetic spin fluctuations are much weaker in the electron-doped cuprates due to weaker nesting of the Fermi surface and less density of states due to the flat band well below the Fermi level. Regarding the superconducting state note, that in the hole-doped cuprates a strong renormalization of the spin fluctuation spectra occurs due to the feedback effect of superconductivity and that $`\mathrm{\Delta }_0\omega _{sf}`$ leading to a resonance peak at $`\omega =\omega _{res}`$ (see Fig. 10(b)). To be more precise, a resonance condition

$$\frac{1}{U_{cr}}=\text{Re }\chi _0(𝐪=𝐐,\omega =\omega _{res}),$$
(18)
which signals the occurrence of a spin-density-wave collective mode, must be fulfilled in order to observe a resonance peak . In electron-doped cuprates, the spin excitations do not obey Eq. (18) and thus only a rearrangement of spectral weight occurs below $`T_c`$, but no resonance peak. Therefore, the kink feature is intimately connected with the resonance peak. As we see from Fig. 10(a) there is only a small feedback of superconductivity below $`T_c`$ on $`\text{Im }\chi `$ in the electron-doped cuprates due to $`\omega _{sf}>>\mathrm{\Delta }_0`$. Thus, we find also no kink feature in the superconducting state of electron-doped cuprates in the antinodal direction.

### D Anisotropic scattering rates

Finally we discuss the anisotropy of the scattering rate $`\tau ^1(\omega )`$ of hole-doped cuprates at different points on the Fermi surface. In Fig. 11 we show our results for $`\tau ^1(\omega )`$ at the antinodal point and the nodal point, respectively, for optimal doping (a) and for the overdoped case (b) for various temperatures. In Fig. 11(a) one clearly sees that the scattering rate is very anisotropic on the Fermi surface reflecting the anisotropy of the coupling of elementary excitations to spin fluctuations. In particular, $`\tau ^1(\omega )`$ in the normal state is almost three times larger at the antinodal point than at the nodal point. This agrees with recent ARPES experiments .

Furthermore, we find that $`\text{Im }\mathrm{\Sigma }\omega `$ demonstrating a non-Fermi liquid behavior in the hole-doped cuprates. Below $`T_c`$ at the antinodal point $`\tau ^1(\omega )`$ reveals a strong feedback of superconductivity at energies $`\omega _{res}+\mathrm{\Delta }_0`$. At the nodal point the effect of superconductivity is rather weak. In the overdoped cuprates the anisotropy between nodal and antinodal points is strongly reduced and for $`\omega 0`$ almost disappeared. Most importantly the system then behaves more Fermi-liquid-like. The latter is seen from Fig. 11(b) where one observes a crossover from the $`\text{Im }\mathrm{\Sigma }\omega `$ to the $`\text{Im }\mathrm{\Sigma }\omega ^2`$ behavior. This is also in agreement with experimental observation .

## IV Summary

In summary, we have analyzed the elementary excitations in hole- and electron-doped cuprates and the fingerprints of spin fluctuations on them. The quasiparticles around the antinodal points of the BZ experience the strongest scattering on spin fluctuations yielding a non-Fermi liquid behavior. In agreement with experimental data, we find that coupling of holes to spin fluctuations yields a kink feature in the renormalized energy dispersion.

Possible phonon contribution to the kink feature One of the interpretation of the kink structure in hole-doped cuprates has been the electron-phonon interaction suggested by Lanzara et al.  Indeed, it is clear that phonons would also cause a kink structure in the energy dispersion if one assumes that Eliashberg function $`\alpha ^2F(𝐪,\omega )`$ has the same features as $`\chi (𝐪,\omega )`$, namely peaked at the wave vector Q and at the Debye frequency $`\omega _D`$, i.e. $`\omega =\omega _D\omega _{sf}`$. By analyzing Fig. 2 it is clear that both spin fluctuations and electron-phonon coupling can cause a kink structure. However, in general, one would expect that its position and doping dependence might be different in both cases. For example, only in the case of dominant spin fluctuation coupling the kink structure can be related to INS experiments, i.e. $`\text{Im }\chi (𝐐,\omega )`$, and, furthermore, the kink position is given by $`\omega _{kink}E_{𝐤𝐐}+\omega _{sf}(x)`$. As discussed earlier, the kink feature along the antinodal $`(0,\pi )(\pi ,\pi )`$ direction results from the structure in $`\varphi (\omega )`$. Thus, additional structure in $`\varphi (\omega )`$ due to the electron-phonon interaction (EPI) may also contribute. Therefore, the question remains: How to distinguish between spin fluctuations and phonons as a reason for the kink formation? To answer this question one has to understand how consistent are both scenarios with available experimental data. For example, as was shown Zeyher and Greco , the value of electron-phonon coupling extracted from the kink analysis yields the value of the electron-phonon coupling $`\lambda `$ and Tc that are too low to account for high-$`T_c`$ superconductivity in the cuprates. Furthermore, assuming that the kink structure arises only from the EPI, it is difficult to understand the $`d_{x^2y^2}`$-symmetry of the superconducting order parameter and related observed anisotropy of the kink structure (see Appendix B). Note, only the spin fluctuation scenario yields $`T_c70`$K , a $`d_{x^2y^2}`$-wave order parameter, and a kink feature in qualitative agreement with experiment. Also the doping dependence of the kink is difficult to explain within the phonon scenario. In contrast to Eq. (13) one would expect $`\omega _{kink}E_{𝐤𝐐}+\omega _D(x)`$ in the case of electron-phonon coupling.

Kink structure in the triplet superconductor $`\text{Sr}_2\text{RuO}_4`$? Finally we want to emphasize that the formation of the kink feature due to spin fluctuations should not be restricted to cuprates. For example, the quasi-two-dimensional triplet superconductor $`\text{Sr}_2\text{RuO}_4`$ (isostructural to $`\text{La}_2\text{CuO}_4`$)  reveals pronounced incommensurate antiferromagnetic spin fluctuations at the wave vector $`𝐐_i=(2\pi /3,2\pi /3)`$ and frequency $`\omega _{sf}6`$meV that originates from the nesting properties of the quasi-one-dimensional $`\alpha `$ and $`\beta `$-bands  (see Fig. 12 for an illustration). On general grounds one would expect a kink structure in the renormalized energy dispersion of the quasiparticles. Although the correlation effects are weaker in $`\text{Sr}_2\text{RuO}_4`$ ($`U`$ is smaller), and $`𝐐_i`$ is an incommensurate wave vector, similar conditions as in cuprates are present. Note, the kink feature should occurs at smaller energies than in cuprates due to a lower value of $`\omega _{sf}`$ in the ruthenates. Further experimental studies should test our suggestion.

###### Acknowledgements.

We thank J. Fink, D. Dessau, and D. Fay for helpful discussions. The work of I. E. is supported by the ’Alexander von Humboldt’ foundation and CRDF Rec.007. We further thank the DFG via the SFB 290 and INTAS (Grant No. 01-0654) for financial support.

## A Generalized Eliashberg equations

The interdependence of elementary excitations with spin excitations leads to strong self-energy effects. The corresponding Dyson equation yields the dressed 2x2 matrix Green’s function $`G`$ in terms of the bare Green’s function $`G_0`$ and the self-energy $`\mathrm{\Sigma }`$:

$`G^1(k)`$ $`=`$ $`G_0^1(k)\mathrm{\Sigma }(k)`$ (A1)
$`=`$ $`i\omega _nZ(k)\tau _0\left[ϵ(k)+\xi (k)\right]\tau _3\varphi (k)\tau _1,`$ (A2)
where $`k=(𝐤,i\omega _n)`$. In the FLEX approximation for the Hubbard Hamiltonian the self-energy $`\mathrm{\Sigma }`$ is determined by the following generalized Eliashberg equations:

$`\mathrm{\Sigma }(k)=`$ (A4)
$`{\displaystyle \underset{k^{}}{}}\left[P_s(kk^{})\tau _0G(k^{})\tau _0+P_c(kk^{})\tau _3G(k^{})\tau _3\right]`$
$`=`$ $`{\displaystyle \underset{k^{}}{}}V_{\mathrm{eff}}(kk^{})G(k^{}).`$ (A5)
In order to provide a better understanding of our numerical procedure we show the corresponding Feynman diagrams for $`V_{\mathrm{eff}}`$ in Fig. 13. Within RPA the spin and charge fluctuation interaction are given by

$$P_s=(2\pi )^1U^2\text{Im }(3\chi _s\chi _{s0}),$$
(A6)
with $`\chi _s=\chi _{so}(1U\chi _{s0})^1`$ and

$$P_c=(2\pi )^1U^2\text{Im }(3\chi _c\chi _{c0}),$$
(A7)
with $`\chi _c=\chi _{c0}(1+U\chi _{c0})^1`$. Therefore, the kernel $`I`$ and the spectral functions of the one-particle Green’s function in Eq. (6), $`A_\nu `$, read

$`I(\omega ,\mathrm{\Omega },\omega ^{})={\displaystyle \frac{f(\omega ^{})+b(\mathrm{\Omega })}{\omega +i\delta \mathrm{\Omega }\omega ^{}}}+{\displaystyle \frac{f(\omega ^{})+b(\mathrm{\Omega })}{\omega +i\delta \mathrm{\Omega }\omega ^{}}},`$ (A8)
$`A_\nu (𝐤,\omega )=\pi ^1\text{Im }\left[a_\nu (𝐤,\omega )/D(𝐤,\omega )\right],`$ (A9)
and

$`D=[\omega Z]^2[ϵ_𝐤^0+\xi ]^2\varphi ^2,`$ (A10)
$`a_0=\omega Z,a_3=ϵ_𝐤^0+\xi ,a_1=\varphi .`$ (A11)
In Eq. (A8) $`f`$ and $`b`$ are the Fermi and Bose distribution function, respectively. Finally, the bare susceptibility is calculated from

$`\text{Im }\chi _{s0,c0}`$ $`=`$ $`{\displaystyle \frac{\pi }{N}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega ^{}\left[f(\omega ^{})f(\omega ^{}+\omega )\right]`$ (A14)
$`\times {\displaystyle \underset{𝐤}{}}[N(𝐤+𝐪,\omega ^{}+\omega )N(𝐤,\omega ^{})`$
$`\pm A_1(𝐤+𝐪,\omega ^{}+\omega )A_1(𝐤,\omega ^{})],`$
where we assume that the same itinerant carriers are responsible for the elemenatry excitations and, at the same time, generate the spin excitations. In Eq. (A14) we use $`N(𝐤,\omega )=A_0(𝐤,\omega )+A_3(𝐤,\omega )`$, and the real parts are calculated with the help of the Kramers-Kronig relation. The subtracted terms in $`P_s`$ and $`P_c`$ remove a double counting that occurs in second order. Note that $`V_{\mathrm{eff}}`$ in Eq. (6) is dominated by the exchange of spin fluctuations due to the fact that the system is in the vicinity of an antiferromagnetic phase transition, but the above equations still remain valid in the case where $`\chi _c`$ becomes more important.

Our numerical calculations are performed on a square lattice with $`256\times 256`$ points in the first Brillouin Zone and with $`200`$ points on the real $`\omega `$-axis up to 16$`t`$ on a logarithmic scale. Within our self-consistent procedure the full momentum and frequency

dependence of quantities is kept.

## B Phonons and $`𝐝_{𝐱^\mathrm{𝟐}𝐲^\mathrm{𝟐}}`$-wave superconductivity

In this Appendix we analyse how the magnetic mode which is mainly peaked at $`𝐪=𝐐=(\pi ,\pi )`$ leads to a $`d_{x^2y^2}`$-wave order parameter that is maximal around $`(\pi ,0)`$ and, in particular, to which extend phonons contribute to this result. In general, the generalized Eliashberg equations read after the inclusion of attractive phonons (branch $`i`$) via their spectral function $`\alpha ^2F_i(𝐪,\mathrm{\Omega })`$:

$`\mathrm{\Sigma }_\nu ^{(i)}(𝐤,\omega )=`$ (B3)
$`N^1{\displaystyle \underset{𝐤^{}}{}}{\displaystyle _0^{\mathrm{}}}𝑑\mathrm{\Omega }V_{\mathrm{eff}}(𝐤𝐤^{},\mathrm{\Omega })\alpha ^2F_i(𝐤𝐤^{},\mathrm{\Omega })`$
$`\times {\displaystyle _{\mathrm{}}^+\mathrm{}}d\omega ^{}I(\omega ,\mathrm{\Omega },\omega ^{})A_\nu (𝐤^{},\omega ^{}).`$
For $`\alpha ^2F_i(𝐪,\mathrm{\Omega })`$ we employ a Lorentzian in frequency $`\mathrm{\Omega }`$ around $`\mathrm{\Omega }_0\omega _D`$ (Debye frequency), and a normalized form factor $`F_i(𝐪)`$ peaked at $`𝐪=𝐪_{\mathrm{pair}}`$ as indicated in Fig. 14. The spin fluctuations that are dominating $`V_{\mathrm{eff}}(𝐪,\omega _{sf})`$ are peaked at $`𝐪=𝐐_{\mathrm{pair}}`$.

It is instructive to write down the weak-copuling limit of the $`\widehat{\tau _1}`$-component of Eq. (B3) that reads ($`T=0`$)

$$\mathrm{\Delta }(𝐤)=\underset{𝐤^{}}{}\frac{\left[V_{\mathrm{eff}}(𝐪)\alpha ^2F_i(𝐪)\right]}{2E_k}\mathrm{\Delta }(𝐤),$$
(B4)
where again $`E_k=\sqrt{\mathrm{\Delta }^2(𝐤)+ϵ_k^2}`$ is the dispersion of the quasiparticles in the superconducting state. Note that the contribution to the pairing potential is repulsive for spin

fluctuations and attractive for phonons, respectively. In the case where no phonons would contribute to the Cooper-pairing ($`\alpha ^2F_i(𝐪)=0`$), $`V_{\mathrm{eff}}(𝐪)`$ bridges parts of the Fermi surface where the superconducting order parameter has opposite signs. This momentum dependence of the pairing interaction is indeed required for solving Eq. (B3) and is typical for unconventional superconductivity. Note that for a repulsive and momentum-independent pairing potential, $`V_{\mathrm{eff}}(𝐪)=const`$, no solution of Eq. (B3) can be obtained.

How is the kink related to the pairing mechanism? Physically speaking, the interdependence of elementary excitations that dominate $`V_{\mathrm{eff}}(𝐪)`$, leads to $`d_{x^2y^2}`$-wave Cooper-pairing as well as to the kink structure as observed by ARPES experiments. In other words, the quasiparticles around the hot spots couple strongly to spin fluctuations that leads (a) to a $`d_{x^2y^2}`$-wave order parameter, and (b) the same coupling leads to the kink in the nodal direction that occurs close to the Fermi level where $`𝐐_{\mathrm{pair}}=(\pi ,\pi )`$ as indicated in Fig. 2.

It follows also from Eq. (14) that attractive phonons with a corresponding spectral function $`\alpha ^2F(𝐪)`$ peaked at $`𝐪=𝐪_{\mathrm{pair}}`$ contribute constructively to $`d_{x^2y^2}`$-wave pairing as long as the main pairing interaction is provided by spin fluctuations. However, the kink close to the antinodal points occurs only below $`T_c`$ and is a result of $`\varphi (\omega )`$ that is maximal around $`(0,\pi )`$. Therefore, the kink structure in the antinodal direction is mainly connected to spin excitations peaked at $`𝐐_{\mathrm{pair}}=(\pi ,\pi )`$ and not to the phonon branch peaked at $`𝐪_{\mathrm{pair}}`$.

Note, in the case where no spin fluctuations would be present, i.e. $`V_{\mathrm{eff}}(𝐪)=0`$, the attractive phonon contribution will cancel the minus sign on the RHS of Eq. (B3) yielding an order parameter with $`s`$-wave symmetry. Thus, we safely conclude that both, $`d_{x^2y^2}`$-wave Cooper-pairing and the anisotropy of the kink feature in the elementary excitations are hardly to recoincile within the same physical picture.


format: text

text:
# STREAMING CAUSED BY IMPULSE ACOUSTIC WAVE.

## 1 INTRODUCTION

Traditionally, streaming is referred as a localized mean flow caused by an acoustic source. This flow is created by the Reynolds forces, non-zero time-averaged values of the quadratic acoustic terms arising in the equation of conservation of momentum. There exist extensive reviews on this subject \[1-3\].

Steady acoustic streaming is described by the mean Eulerian velocity $`\stackrel{}{v}`$, or the mean mass flow velocity $`\stackrel{}{U}`$, where $`\stackrel{}{U}`$ is defined by $`\stackrel{}{U}=\frac{\rho \stackrel{}{v}}{\rho _0}=\stackrel{}{V}+\frac{\rho _a\stackrel{}{v}_a}{\rho _0}.`$Here, $`\stackrel{}{V}`$ is a slow component of velocity, namely the Eulerian velocity of the streaming, index $`a`$ labels acoustic wave, $`\rho _0`$ is undisturbed density. Square brackets designate averaging over time whose interval is much shorter than the transient time of streaming and much longer than a period of a sound wave. It is pointed, that the difference between $`\stackrel{}{U}`$ and $`\stackrel{}{V}`$ is roughly proportional to the acoustic intensity vector and is small compared with any of mean velocities . From the continuity equation it follows $`\stackrel{}{U}=0`$ if we assume $`\frac{\rho }{t}=0`$.

The recent discussion of acoustic streaming  suggest that there is important unresolved issue concerning acoustic streaming, namely the effect of fluid incompressibility. The inconsistency of this point is that fluid is supposed to be incompressible while the acoustic wave casing streaming may propagate only over compressible medium. The next, the conservation of energy usually is not considered. It is reasonable since excluding of one variable (density) reduces an initial system. It has been proved that the effect of the heat conductivity could not be discarded in a study of temperature variation associated with the streaming . Actually, this approximation is well understood in a typical liquid like water but should be revised for other liquids . That is a reason to start from the full system of conservation laws including the energy balance equation and the continuity equation.

The procedure of the temporal averaging needs a continuous acoustic wave of many periods as a source, thus the well-known results on acoustic streaming apply only for this type of acoustic source though quite realistic. There are clear inconsistency concerning the very procedure of the temporal averaging. As it has been mentioned before, an equality $`\stackrel{}{U}=0`$ is correct only if $`\frac{\rho }{t}=0`$ . Actually, the overall field includes not only acoustic and vortical modes, but also of the entropy mode. The entropy mode caused by acoustic wave is known as a slow decrease of density and isobaric increase of temperature of the background of acoustic wave propagation. Recently, the importance of change in density of surrounding (not only increase in temperature) was underlined in the review . If nonlinear interactions are accounted in a proper way, for density of the entropy mode one gets an equation :

$`\frac{\rho _e}{t}+\frac{\eta }{\rho _0}\mathrm{\Delta }\rho _e=S_a`$,

where $`\eta `$ is shear viscosity, $`S_a`$ is an acoustic source proportional to a gradient of acoustic intensity, index $`e`$ concerns the entropy mode. An attenuation of acoustic wave leads to slow but continuous decrease of density so $`\frac{\rho _e}{t}=0`$ becomes incorrect as well as $`\frac{\rho }{t}=0,`$ since $`\rho =\rho _0+\rho _a+\rho _e.`$ Therefore, if the acoustic heating caused by periodic acoustic wave is accounted, $`\stackrel{}{U}=0`$ becomes incorrect.

These inconsistencies call new ways to theoretical investigation of acoustic streaming, including streaming caused by non-periodic acoustic wave. We plan to start with all conservation equations for compressible thermoviscous flow, then to define modes by establishing of the relations between perturbations (section 2), to demonstrate how the well-known linear and nonlinear evolution equations may be obtained by projecting (sections 3,4), and, finally, to apply projectors to calculating of the streaming caused by pulse acoustic wave (section 5).

The idea to decompose the linear flow into specific modes is not novel and has been exploited for a long time , see  and referred papers, where homogeneous background with sources of heat and mass are considered. The concrete ideas of doing it automatically by using of projectors in the wide variety of problems (also flow over inhomogeneous media like bubbly liquid , flows affected by external forces including the gravitational one which changes the background density and pressure , both one- and multi-dimensional problems, see  ). The principle advance is an expansion of the ideas into area of nonlinear flow: to get nonlinear coupled evolution equations for interacting modes and to solve the system approximately. The deriving of the final nonlinear system is algorithmic, one just acts by projectors at the initial system of nonlinear equations. Some ideas concerning separation and interaction of acoustic-gravity and Rossby waves may be found in .

## 2  Modes and projectors for the flow with equations of state in the general form.

The mass, momentum and energy conservation equations for the thermoviscous flow read:

$$\frac{\rho }{t}+\stackrel{}{}(\rho \stackrel{}{v})=0$$

$$\rho \left[\frac{\stackrel{}{v}}{t}+(\stackrel{}{v}\stackrel{}{})\stackrel{}{v}\right]=\stackrel{}{}p+\eta \mathrm{\Delta }\stackrel{}{v}+\left(\varsigma +\frac{\eta }{3}\right)\stackrel{}{}(\stackrel{}{}\stackrel{}{v})$$

$$\rho \left[\frac{e}{t}+(\stackrel{}{v}\stackrel{}{})e\right]+p\stackrel{}{}\stackrel{}{v}\chi \mathrm{\Delta }T=\varsigma \left(\stackrel{}{}\stackrel{}{v}\right)^2+\frac{\eta }{2}\left(\frac{v_i}{x_k}+\frac{v_k}{x_i}\frac{2}{3}\delta _{ik}\frac{v_l}{x_l}\right)^2$$
$`(1)`$
Here, among already mentioned variables, $`p`$ is pressure, $`e,T`$\- internal energy per unit mass and temperature; $`\varsigma ,`$ $`\chi `$are bulk viscosity and thermal conductivity coefficient (both supposed to be constants as well as $`\eta `$ ), $`x_i`$ \- space coordinates. Except of the dynamic equations (1), the two thermodynamic relations are necessary: $`e(p,\rho ),T(p,\rho )`$ . To treat a wide variety of substances, let us use the most general form of these relations as expansion in the Fourier series:

$$\rho _0e^{}=E_1p^{}+\frac{E_2p_0}{\rho _0}\rho ^{}+\frac{E_3}{p_0}p^2+\frac{E_4p_0}{\rho _0^2}\rho ^2+\frac{E_5}{\rho _0}p^{}\rho ^{}+\mathrm{}$$

$$T^{}=\frac{\mathrm{\Theta }_1}{\rho _0C_v}p^{}+\frac{\mathrm{\Theta }_2p_0}{\rho _0^2C_v}\rho ^{}+\mathrm{}$$
$`(2)`$
The background values are marked by zero, perturbations are primed, $`C_V`$ means specific heat per unit mass at constant volume, $`E_1,\mathrm{}\mathrm{\Theta }_1,\mathrm{}`$are dimensionless coefficients,

$`\mathrm{\Theta }_1=\frac{\rho _0C_vk}{\beta },`$ $`\mathrm{\Theta }_2=\frac{\rho _0C_v\beta }{p_0},`$ (3)

where $`k=\frac{1}{\rho _0}\left(\frac{p}{\rho }\right)_{T=T(p_0,\rho _0)},\beta =\frac{1}{\rho _0}\left(\frac{\rho }{T}\right)_{\rho =\rho _0}.`$

For the tree-dimensional flow $`\stackrel{}{v}=(v_x,v_y,v_z)`$ in the Cartesian coordinates $`\stackrel{}{x}=(x,y,z)`$ over homogeneous stationary background without mean flow $`\stackrel{}{v}_0`$ $`=0`$ , the system (1), (2) looks:

$$\rho _0\frac{v_x}{t}+\frac{p}{x}\eta \mathrm{\Delta }v_x\left(\varsigma +\frac{\eta }{3}\right)\frac{}{x}(\stackrel{}{}v)=\rho _0(\stackrel{}{v}\stackrel{}{})v_x+\frac{\rho ^{}}{\rho _0}\frac{p}{x}\eta \rho ^{}\mathrm{\Delta }v_x\left(\varsigma +\frac{\eta }{3}\right)\rho ^{}\frac{}{x}(\stackrel{}{}v)$$

$$\rho _0\frac{v_y}{t}+\frac{p}{y}\eta \mathrm{\Delta }v_y\left(\varsigma +\frac{\eta }{3}\right)\frac{}{y}(\stackrel{}{}v)=\rho _0(\stackrel{}{v}\stackrel{}{})v_y+\frac{\rho ^{}}{\rho _0}\frac{p}{y}\eta \rho ^{}\mathrm{\Delta }v_y\left(\varsigma +\frac{\eta }{3}\right)\rho ^{}\frac{}{y}(\stackrel{}{}v)$$

$$\rho _0\frac{v_x}{t}+\frac{p}{x}\eta \mathrm{\Delta }v_z\left(\varsigma +\frac{\eta }{3}\right)\frac{}{z}(\stackrel{}{}v)=\rho _0(\stackrel{}{v}\stackrel{}{})v_z+\frac{\rho ^{}}{\rho _0}\frac{p}{z}\eta \rho ^{}\mathrm{\Delta }v_z\left(\varsigma +\frac{\eta }{3}\right)\rho ^{}\frac{}{z}(\stackrel{}{}v)$$

$$\frac{p^{}}{t}+c^2\rho _0\stackrel{}{}\stackrel{}{v}\frac{\chi }{E_1}\left(\frac{\mathrm{\Theta }_1}{\rho _0C_v}\mathrm{\Delta }p^{}+\frac{\mathrm{\Theta }_2p_0}{\rho _0^2C_v}\mathrm{\Delta }\rho ^{}\right)=\frac{1}{E_1}\left(\varsigma \left(\stackrel{}{}\stackrel{}{v}\right)^2+\frac{\eta }{2}\left(\frac{v_i}{x_k}+\frac{v_k}{x_i}\frac{2}{3}\delta _{ik}\frac{v_l}{x_l}\right)^2\right)+\left(Qp^{}+S\rho ^{}\right)\stackrel{}{}\stackrel{}{v}\stackrel{}{v}\stackrel{}{}p$$

$$\frac{\rho ^{}}{t}+\rho _0\stackrel{}{}\stackrel{}{v}=\rho ^{}\stackrel{}{}\stackrel{}{v}\stackrel{}{v}\stackrel{}{}\rho ^{}$$
$`(4),`$
where Q,S are constants depending on the equations of state (2) :

$$Q=\frac{1}{E_1}(1+2\frac{1E_2}{E_1}E_3+E_5),$$

$$S=\frac{1}{1E_2}(1+E_2+2E_4+\frac{1E_2}{E_1}E_5).$$
$`(5)`$
For an ideal gas, $`Q=\gamma =C_P/C_v,`$ $`S=0.`$

In order to simplify calculations, the geometry of beams will be considered. The equivalent system in the dimensionless variables :

$$\stackrel{}{v}_{},\stackrel{}{x}_{},\rho _{},p_{},t_{}:\stackrel{}{v}=\epsilon c\stackrel{}{v}_{},p^{}=\epsilon c^2\rho _0p_{},\rho ^{}=\epsilon \rho _0\rho _{},\stackrel{}{x}=(\lambda x_{}/\sqrt{\mu },\lambda y_{},\lambda z_{}/\sqrt{\mu }),t=\lambda t_{}/c,$$
$`(6)`$
($`c=\sqrt{\frac{p_0(1E_2)}{\rho _0E_1}}`$is adiabatic sound velocity, $`\lambda `$ means characteristic scale of the longitudinal perturbations, $`\mu `$ is small parameter expressing the relation the longitudinal and transverse scales of perturbation, $`\epsilon `$ is small amplitude parameter) looks as follows (asterisks for dimensionless variables will be later omitted):

$$\frac{}{t}\psi +L\psi =\epsilon \phi +\epsilon \phi _{tv},$$
$`(7)`$
where $`\psi `$ is column of field perturbations

$$\psi =\left(\begin{array}{ccccc}v_x& v_y& v_z& p& \rho \end{array}\right)^T,$$
$`(8)`$
$`L`$ is a linear matrix operator :

$$L=\left(\begin{array}{ccccc}\delta _1^1\mu \frac{^2}{x^2}\delta _1^2\mathrm{\Delta }& \delta _1^1\sqrt{\mu }\frac{^2}{xy}& \delta _1^1\mu \frac{^2}{xz}& \sqrt{\mu }/x& 0\\ \delta _1^1\sqrt{\mu }\frac{^2}{xy}& \delta _1^1\frac{^2}{y^2}\delta _1^2\mathrm{\Delta }& \delta _1^1\sqrt{\mu }\frac{^2}{yz}& /y& 0\\ \delta _1^1\mu \frac{^2}{xz}& \delta _1^1\sqrt{\mu }\frac{^2}{zy}& \delta _1^1\mu \frac{^2}{z^2}\delta _1^2\mathrm{\Delta }& \sqrt{\mu }/z& 0\\ \sqrt{\mu }/x& /y& \sqrt{\mu }/z& \delta _2^1\mathrm{\Delta }& \delta _2^2\mathrm{\Delta }\\ \sqrt{\mu }/x& /y& \sqrt{\mu }/z& 0& 0\end{array}\right)$$
$`(9)`$
with dimensionless parameters

$`\delta _1^1=\frac{\left(\zeta +\eta /3\right)}{\rho _0c\lambda },\delta _1^2=\frac{\eta }{\rho _0c\lambda },\delta _2^1=\frac{\chi \theta _1}{\rho _0c\lambda _0C_vE_1},\delta _2^2=\frac{\chi \theta _2}{\rho _0c\lambda _0C_v(1E_2)}.`$

Dimensionless operators $`\stackrel{}{},\mathrm{\Delta }`$ look : $`\stackrel{}{}=\left(\begin{array}{ccc}\mu /x& /y& \mu /z\end{array}\right),`$ $`\mathrm{\Delta }=\mu ^2/x^2+^2/y^2+\mu ^2/z^2,`$

$`\phi `$ is the quadratic nonlinear column:

$$\phi =\left(\begin{array}{c}(\stackrel{}{v}\stackrel{}{})v_x+\sqrt{\mu }\rho p/x\\ (\stackrel{}{v}\stackrel{}{})v_y+\rho p/y\\ (\stackrel{}{v}\stackrel{}{})v_z+\sqrt{\mu }\rho p/z\\ [Qp+S\rho ](\stackrel{}{}\stackrel{}{v})(\stackrel{}{v}\stackrel{}{})p\\ \rho (\stackrel{}{}\stackrel{}{v})(\stackrel{}{v}\stackrel{}{})\rho \end{array}\right)$$
$`(10),`$
and $`\phi _{tv}`$ is the quadratic nonlinear column $`O(\beta )`$ appearing in the viscous flow:

$$\phi _{tv}=\left(\begin{array}{c}\delta _1^2\rho \mathrm{\Delta }v_x\delta _1^1\rho \frac{}{x}(\stackrel{}{}v)\\ \delta _1^2\rho \mathrm{\Delta }v_y\delta _1^1\rho \frac{}{y}(\stackrel{}{}v)\\ \delta _1^2\rho \mathrm{\Delta }v_z\delta _1^1\rho \frac{}{z}(\stackrel{}{}v)\\ \frac{1}{E_1}\left(\left(\delta _1^1\delta _1^2/3\right)\left(\stackrel{}{}\stackrel{}{v}\right)^2+\frac{\delta _1^2}{2}\left(\frac{v_i}{x_k}+\frac{v_k}{x_i}\frac{2}{3}\delta _{ik}\frac{v_l}{x_l}\right)^2\right)\\ 0\end{array}\right)$$
$`(11)`$
For a linear flow defined by the linearized version of system (7)

$$\frac{}{t}\psi +L\psi =0,$$
$`(12)`$
a solution may be found as a sum of planar waves: $`v_x=\stackrel{~}{v_x}(\stackrel{}{k})\mathrm{exp}(i\omega ti\stackrel{}{k}\stackrel{}{x}),\mathrm{}.`$where

$$\stackrel{}{k}=(k_x,k_y,k_z)$$

is the wave vector. It is convenient to go to the Fourier transforms marked by tilde for the flows over the homogeneous background. In the Fourier space,$`ik_x`$ means $`/x,`$ $`i\omega `$ means $`/t,`$and (6) yields in the five roots of dispersion relation $`(\beta =\delta _1^1+\delta _1^2+\delta _2^1+\delta _2^2)`$:

$`\omega _1=\mathrm{\Omega }+i\beta \mathrm{\Omega }^2/2,`$ $`\omega _2=\mathrm{\Omega }+i\beta \mathrm{\Omega }^2/2,`$ $`\omega _3=i\delta _2^2\mathrm{\Omega }^2,`$ $`\omega _4=i\delta _1^2\mathrm{\Omega }^2,`$ $`\omega _5=i\delta _1^2\mathrm{\Omega }^2,`$where

$$\mathrm{\Omega }=k_y+\frac{\mu (k_x^2+k_z^2)}{2k_y}.$$
$`(13)`$
Also, an operator $`\left(ik_y\right)^1`$ represents in $`\stackrel{}{k}`$-space operator $`𝑑y`$ . Constant of integration should be chosen accordingly to the concrete physical problem. The first two roots relate to progressive (acoustic) modes of different directions of propagation, the third one relates to the entropy mode, and the fourth and fifth ones - to the vortical ones. For the real substances, $`\beta >0,`$and, $`\delta _2^2<0,`$that corresponds physically correct sign of imagine parts of all frequencies. The modes of linear flow are determined by relations of amplitudes of plane waves $`\stackrel{~}{v_x}(k_x,k_y,k_z)`$ ,$`\mathrm{}`$ Modes as eigenvectors of a linear problem in the Fourier space look:

$`\stackrel{~}{\psi }_1=\left(\begin{array}{c}\stackrel{~}{v_x}_1(k_x,k_y,k_z)\\ \stackrel{~}{v_y}_1(k_x,k_y,k_z)\\ \stackrel{~}{v_z}_1(k_x,k_y,k_z)\\ \stackrel{~}{p}_1(k_x,k_y,k_z)\\ \stackrel{~}{\rho }_1(k_x,k_y,k_z)\end{array}\right)=\left(\begin{array}{c}\sqrt{\mu }k_x/k_y\\ 1\mu (k_x^2+k_z^2)/(2k_y^2)+i\beta k_y/2\\ \sqrt{\mu }k_z/k_y\\ 1+i(\delta _2^1+\delta _2^2)k_y\\ 1\end{array}\right)\stackrel{~}{\rho }_1=\stackrel{~}{M}_1\stackrel{~}{\rho }_1,`$

$$\stackrel{~}{\psi }_2=\left(\begin{array}{c}\sqrt{\mu }k_x/k_y\\ 1+\mu (k_x^2+k_z^2)/(2k_y^2)+i\beta k_y/2\\ \sqrt{\mu }k_z/k_y\\ 1i(\delta _2^1+\delta _2^2)k_y\\ 1\end{array}\right)\stackrel{~}{\rho }_2=\stackrel{~}{M}_2\stackrel{~}{\rho }_2,$$
$`(14)`$
$`\stackrel{~}{\psi }_3=\left(\begin{array}{c}0\\ i\delta _2^2k_y\\ 0\\ 0\\ 1\end{array}\right)\stackrel{~}{\rho }_3=\stackrel{~}{M}_3\stackrel{~}{\rho }_3,`$ $`\stackrel{~}{\psi }_4=\left(\begin{array}{c}ik_y\\ i\sqrt{\mu }k_x\\ 0\\ 0\\ 0\end{array}\right)\stackrel{~}{\phi }_4=\stackrel{~}{M}_4\stackrel{~}{\phi }_4,\stackrel{~}{\psi }_5=\left(\begin{array}{c}0\\ i\sqrt{\mu }k_z\\ ik_y\\ 0\\ 0\end{array}\right)\stackrel{~}{\phi }_5=\stackrel{~}{M}_5\stackrel{~}{\phi }_5.`$

All calculations of modes and projectors have accuracy up to the terms of order $`\mu `$, $`\beta `$. As basic variable for the first three modes the perturbation of density is chosen, and for the last two- the stream function, since for the both vorticity modes density keeps unperturbed. So, there are five independent values indeed, one for every mode, that define other wave perturbations of the mode uniquely. Any field of the linear flow as a solution of a linearized equation (12) may be presented as a sum of independent modes.

The next step is to get projectors that decompose a concrete mode from the overall field $`\stackrel{~}{\psi }.`$ Let us define a matrix M in the following way:

$$\stackrel{~}{M}\left(\begin{array}{c}\stackrel{~}{\rho }_1\\ \stackrel{~}{\rho }_2\\ \stackrel{~}{\rho }_3\\ \stackrel{~}{\phi }_4\\ \stackrel{~}{\phi }_4\end{array}\right)=\stackrel{~}{\psi }=\underset{i=1}{\overset{5}{}}\stackrel{~}{\psi }_i=\stackrel{~}{M}_1\stackrel{~}{\rho }_1+\stackrel{~}{M}_2\stackrel{~}{\rho }_2+\stackrel{~}{M}_3\stackrel{~}{\rho }_3+\stackrel{~}{M}_4\stackrel{~}{\phi }_4+\stackrel{~}{M}_5\stackrel{~}{\phi }_5.$$
$`(15)`$
so that matrix $`\stackrel{~}{M}`$ consists of five columns,

$$\stackrel{~}{M}=(\begin{array}{ccccc}\stackrel{~}{M}_1& \stackrel{~}{M}_2& \stackrel{~}{M}_3& \stackrel{~}{M}_4& \stackrel{~}{M}_5\end{array}).$$
$`(16)`$
and the inverse matrix $`\stackrel{~}{M}^1`$ consists of five lines

$`\stackrel{~}{M}^1=\left(\begin{array}{ccccc}\stackrel{~}{M}_1^1& \stackrel{~}{M}_2^1& \stackrel{~}{M}_3^1& \stackrel{~}{M}_4^1& \stackrel{~}{M}_5^1\end{array}\right)^T,`$

where lines $`\stackrel{~}{M}_1^1,..\stackrel{~}{M}_5^1`$ are as follows:

$`\stackrel{~}{M}_1^1=(\begin{array}{ccccc}\sqrt{\mu }\frac{k_x}{k_y}& \frac{1}{2}\left(1i(\delta _2^1+\delta _2^2)k_y\mu (k_x^2+k_z^2)/(2k_y)\right)& \sqrt{\mu }\frac{k_z}{k_y}& \frac{1}{2}\left(1i\beta k_y/2i\delta _2^2k_y\right)& i\delta _2^2k_y/2\end{array})`$,

$`\stackrel{~}{M}_2^1=(\begin{array}{ccccc}\sqrt{\mu }\frac{k_x}{k_y}& \frac{1}{2}\left(1i(\delta _2^1+\delta _2^2)k_y+\mu (k_x^2+k_z^2)/(2k_y)\right)& \sqrt{\mu }\frac{k_z}{k_y}& \frac{1}{2}\left(1+i\beta k_y/2+i\delta _2^2k_y\right)& i\delta _2^2k_y/2\end{array}),`$

$`\stackrel{~}{M}_3^1=(\begin{array}{ccccc}0& \frac{1}{2}i(\delta _2^1+\delta _2^2)k_y& 0& 1& 1\end{array}),`$

$`\stackrel{~}{M}_4^1=(\begin{array}{ccccc}i\frac{1\mu k_x^2/k_y^2}{k_y}& i\sqrt{\mu }\frac{k_x}{k_y^2}& i\mu \frac{k_xk_z}{k_y^3}& 0& 0\end{array}),`$

$`\stackrel{~}{M}_5^1=(\begin{array}{ccccc}i\mu \frac{k_xk_z}{k_y^3}& i\sqrt{\mu }\frac{k_z}{k_y^2}& i\frac{1\mu k_z^2/k_y^2}{k_y}& 0& 0\end{array}).`$

$`(17)`$
Accordingly to definition of matrix $`\stackrel{~}{M}^1`$ projectors may be determined, so as

$`\stackrel{~}{P}_1\stackrel{~}{\psi }=\stackrel{~}{\psi }_1,\mathrm{},\stackrel{~}{P}_5\stackrel{~}{\psi }=\stackrel{~}{\psi }_5.`$ (18)

Relations (14), (15) written for the first mode as an example $`\stackrel{~}{M}_1^1\stackrel{~}{\psi }=\rho _1,\stackrel{~}{\psi }_1=\stackrel{~}{M}_1\rho _1`$ lead to $`\stackrel{~}{P}_1=\stackrel{~}{M}_1\stackrel{~}{M}_1^1`$ accordingly to (18) and so on for the all other modes:

$$\stackrel{~}{P}_i=\stackrel{~}{M}_i\stackrel{~}{M}_i^1,i=1,\mathrm{..5}.$$
$`(19)`$
At any moment of evolution a concrete mode is distinguished from the overall field by the correspondent projector. Projectors calculated with accuracy of order $`\mu ,\beta `$ look:

$`\stackrel{~}{P}_{1,2}=\left(\begin{array}{ccccc}\mu \frac{k_x^2}{2k_y^2}& \sqrt{\mu }\frac{k_x}{2k_y}& \mu \frac{k_xk_z}{2k_y^2}& \pm \sqrt{\mu }\frac{k_x}{2k_y}& 0\\ \sqrt{\mu }\frac{k_x}{2k_y}& \frac{1}{2}\left(1\pm \frac{i\beta }{2}k_yi(\delta _2^1+\delta _2^2)k_y\mu \frac{k_x^2+k_z^2}{2k_y^2}\right)& \sqrt{\mu }\frac{k_z}{2k_y}& \frac{1}{2}\left(\pm 1i\delta _2^2k_y\mu \frac{k_x^2+k_z^2}{2k_y^2}\right)& \frac{i\delta _2^2k_y}{2}\\ \mu \frac{k_xk_z}{2k_y^2}& \sqrt{\mu }\frac{k_z}{2k_y}& \mu \frac{k_z^2}{2k_y^2}& \pm \sqrt{\mu }\frac{k_z}{2k_y}& 0\\ \pm \sqrt{\mu }\frac{k_x}{2k_y}& \pm \frac{1}{2}\left(1\mu \frac{k_x^2+k_z^2}{2k_y^2}\right)& \pm \sqrt{\mu }\frac{k_z}{2k_y}& \frac{1}{2}(1\frac{i\beta }{2}k_y\pm i\delta _2^1k_y)& \pm \frac{i\delta _2^2k_y}{2}\\ \pm \sqrt{\mu }\frac{k_x}{2k_y}& \frac{1}{2}\left(\pm 1i(\delta _2^1+\delta _2^2)k_y\mu \frac{k_x^2+k_z^2}{2k_y^2}\right)& \pm \sqrt{\mu }\frac{k_z}{2k_y}& \frac{1}{2}(1\frac{i\beta }{2}k_yi\delta _2^2k_y)& \pm \frac{i\delta _2^2k_y}{2}\end{array}\right),`$

$$\stackrel{~}{P}_3=\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& i\delta _2^2k_y& i\delta _2^2k_y\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& i(\delta _2^1+\delta _2^2)k_y& 0& 1& 1\end{array}\right),$$
$`(20)`$
$`\stackrel{~}{P}_4=\left(\begin{array}{ccccc}1\mu \frac{k_x^2}{k_y^2}& \sqrt{\mu }\frac{k_x}{k_y}& \mu \frac{k_xk_z}{k_y^2}& 0& 0\\ \sqrt{\mu }\frac{k_x}{k_y}& \mu \frac{k_x^2}{k_y^2}& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right),\stackrel{~}{P}_5=\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& \mu \frac{k_z^2}{k_y^2}& \sqrt{\mu }\frac{k_z}{k_y}& 0& 0\\ \mu \frac{k_xk_z}{k_y^2}& \sqrt{\mu }\frac{k_z}{k_y}& 1\mu \frac{k_z^2}{k_y^2}& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right),.`$

Matrix projectors satisfy common properties of orthogonal projectors:

$`\underset{i=1}{\overset{5}{}}\stackrel{~}{P}_i=\stackrel{~}{I},`$ $`\stackrel{~}{P}_i\stackrel{~}{P}_n=\stackrel{~}{0}`$ if $`in,\stackrel{~}{P}_i\stackrel{~}{P}_i=`$ $`\stackrel{~}{P}_i`$ $`if`$ $`i=n,`$ (21)

where $`\stackrel{~}{I},\stackrel{~}{0}`$ are unit and zero matrices. The inverse transformation of formulae (20) to the $`(\stackrel{}{x},t)`$ space may be easily undertaken.

## 3  Linear flow and evolution equations.

A linear flow may be decomposed to modes uniquely accordingly to (18). Linear evolution equations for every mode may be originated from the linear version of (7) by acting of the corresponding projector: $`P\left(\frac{}{t}\psi +L\psi \right)=0.`$ Note that all projectors do commute with both $`/t`$ $`I`$ and $`L`$ ($`I`$ is unit matrix). When one acts by projector at the linear system (12), five equations for all the components of every mode appear. For example, an evolution equation for the dimensionless perturbations of density for both acoustic modes are as follows:

$`\frac{\rho _{1,2}}{t}\pm \frac{\rho _{1,2}}{y}\pm \frac{\mu }{2}\mathrm{\Delta }_{}\rho _{1,2}𝑑y\frac{\beta }{2}\frac{^2\rho _{1,2}}{y^2}=0,`$ (22)

where $`\mathrm{\Delta }_{}=\frac{^2}{x^2}+\frac{^2}{z^2}.`$

An evolution equation for the axial velocity of the entropy mode looks:

$`\frac{v_{y3}}{t}\delta _2^2\frac{^2v_{y3}}{y^2}=0,`$ (23)

and for the axial velocity of two vortical modes is:

$`\frac{v_{y4,5}}{t}\delta _1^2\frac{^2v_{y4,5}}{y^2}=0.`$ (24)

These linear evolution equations are well-known .Equations (21)-(23) are explicit up to terms of order $`O(\beta \sqrt{\mu },\mu ^2)`$ due to accuracy of evaluated projectors. In the paper , equations (22)-(24) are presented as three equations for the three basic modes, since there was no subdivision of acoustic and vortical modes into two independent branches everyone. Also, the evolution equations from  involves sources of mass, momentum and heat that easily may be added in to the right-hand side of the initial system (7) as well as to the equations (22)-(24).

## 4 Nonlinear flow and coupled evolution equations.

At this point, we declare that relations inside every mode become fixed in weakly nonlinear flow. Acting by projectors at the original system (7) with non-zero nonlinear part leads to coupled equations for interacting modes: $`P\left(\frac{}{t}\psi +L\psi \right)=\epsilon P\left(\phi _1+\phi _{1tv}\right)`$. Modes become separated in the linear side, and one should account that the overall field is a sum of all modes inputs in the nonlinear side. The proper returning to the $`\stackrel{}{x}`$ space should be proceeded as well. To illustrate the possibilities of the projecting method, some illustrations concerning the famous evolution equations deriving are presented below.

### 4.1 One-dimensional nonlinear flow.

Let $`k_x=0,`$ $`k_z=0,`$ that corresponds to one-dimensional flow along y-axis. For the planar geometry, only three modes exist: two acoustic and the entropy one. Acting by $`P_1`$ at the system (5), one get an evolution equation for the rightward acoustic mode consisting of three equations for , $`p_1,\rho _1,`$ and $`v_1`$. An evolution equation for the non-dimensional density perturbations of the rightward acoustic mode is as follows:

$$\frac{\rho _1}{t}+\frac{\rho _1}{y}\frac{\beta }{2}\frac{^2v_1}{y^2}=\frac{\epsilon }{2}\left(v\frac{}{y}v+\rho \frac{}{y}pv\frac{}{y}p+\left[Qp+S\rho \right]\frac{}{y}v\right)$$
$`(25),`$
Here, $`p,\rho ,`$ and $`v`$ in the right-hand part represent a sum of all modes:

$`v=v_1+v_2+v_3,`$ $`p=p_1+p_2+p_3,`$

$$\rho =\rho _1+\rho _2+\rho _3.$$
$`(26)`$
All cross nonlinear-viscous terms are neglected, as well as higher order nonlinear ones. If the nonlinear terms of the first mode are kept only, we get the known Burgers evolution equation, going to the famous Earnshow one in the limit $`\beta 0`$ :

$$\frac{\rho _1}{t}+\frac{\rho _1}{y}+\epsilon \frac{QS+1}{2}\rho _1\frac{}{y}\rho _1\frac{\beta }{2}\frac{^2\rho _1}{y^2}=0.$$
$`(27)`$
Physically, this case relates to a dominant rightward acoustic mode. For an ideal gas, again, $`Q=\gamma =C_P/C_v,`$

$`S=0.`$

### 4.2 Dynamics of beams.

In the same way, the famous Khokhlov-Zabolotskaya equation  (KZ, in the limit $`\beta 0`$) and Khokhlov-Zabolotskaya-Kuznetsov equation (KZK) for rightward (leftward) propagating beams follows when acting by $`P_{1,}P_2`$ at the system (7):

$$\frac{\rho _{1,2}}{t}\pm \frac{\rho _{1,2}}{y}\pm \frac{\mu }{2}\mathrm{\Delta }_{}\rho _{1,2}𝑑y\frac{\beta }{2}\frac{^2\rho _{1,2}}{y^2}\pm \epsilon \frac{QS+1}{2}\rho _{1,2}\frac{}{y}\rho _{1,2}=0,$$
$`(28)`$
This evolution equation may be understood as a limit for the progressive mode self-action, meaning that only this mode nonlinear terms are kept. It seems be suitable definition though the ’self-action’ is usually reserved for the phenomenon of acoustic heating when acoustic mode causes the entropy mode during its propagation thus changing the surrounding .

## 5 ACOUSTIC STREAMING CAUSED BY IMPULSES.

There are some reasons to find another approach to acoustic streaming. First, the procedure of temporal averaging is a basis of the modern theory of acoustic streaming and all corresponding results relate to averaged fields. All non-periodic acoustic sources thus are left of account. Second, the problem of modes interaction looks more extended. Actually there exist three types of independent modes. The initial amplitudes of these modes are determined by initial conditions of a concrete problems and may be evaluated by projecting of the overall initial field into every modal field. The acoustic streaming imposes an ability of predominate acoustic mode. The acoustic field causes slow varying with time flow which may grow with time since nonlinear effects can store, that is namely the acoustic streaming. So, the acoustic streaming relates to one type (though important) of possible nonlinear interactions in flow, namely to the beginning of acoustic mode evolution. Later, an amplitude of the vortical mode (and the entropy one also) may become so large that one can not treat these modes as secondary and other approach for nonlinear interaction of modes should be developed. Let us mark, that even traditional subdivision of overall flow into ’slow’ and ’quick’ components eliminates the entropy mode unclearly. It is also ’slow’ but essentially possesses non-zero density perturbation which slowly varies with time ( $`\rho _3/t=0`$ in linear approach), so the procedure of temporal averaging of continuity equation fails with storing of nonlinear effects. Also, there are many more complicated flows over inhomogeneous media and/or with background flows that may be algorithmically solved by the pointed method, see conclusion.

We will show, at first, that projecting results in the well-known equations for an acoustic streaming in terms of $`\stackrel{}{V},`$ and, the second and more important, that there appear the new equations for acoustic streaming caused by impulse (non-periodic) acoustic signals.

To consider interactions of vortical and acoustic modes, one acts by projector $`P_4`$ $`+`$ $`P_5`$ at the system (7) to get an evolution equation for vorticity $`\phi _4`$ $`+`$ $`\phi _5`$ . The right-hand nonlinear vector involves, again, all modes inputs. Then, to calculate the nonlinear interactions of only vortical and rightward propagating acoustic mode, the right-hand vector should be thought as a sum of vortical and the first acoustic mode:

$$\psi =\left(\begin{array}{c}\sqrt{\mu }\frac{}{x}𝑑y\\ 1\frac{\mu }{2}\mathrm{\Delta }_{}𝑑y𝑑y\frac{\beta }{2}\frac{}{y}\\ \sqrt{\mu }\frac{}{z}𝑑y\\ 1(\delta _2^1+\delta _2^2)\frac{}{y}\\ 1\end{array}\right)\rho _1+\left(\begin{array}{c}v_{x4}\\ v_{y4}\\ 0\\ 0\\ 0\end{array}\right)+\left(\begin{array}{c}0\\ v_{y5}\\ v_{z5}\\ 0\\ 0\end{array}\right).$$
$`(29)`$
An operator $`\mathrm{\Delta }_{}𝑑y𝑑y`$ represents in $`\stackrel{}{x}`$-space operator $`\frac{k_x^2+k_z^2}{k_y^2}.`$Acting by $`P_4+P_5`$ at the system (5) leads to an evolution equation (27) for vorticity $`\phi (\stackrel{}{x},t),`$ which we rewrite going to velocity components $`V_x=v_{x4}=\phi _4/y,`$ $`V_y=v_{y4}+v_{y5}=\sqrt{\mu }\phi _4/x+\sqrt{\mu }\phi _5/z,`$ $`V_z=v_{z5}=\phi _5/y:`$

$$\frac{V_x}{t}\delta _1^2\frac{^2V_x}{y^2}=\epsilon (\stackrel{}{V}\stackrel{}{})V_x+\epsilon F_{1x}$$
$`(30)`$
,

where in right-hand side of (30) only two term are left: the first relating to the vortical mode, and the second $`F_1`$ consisting of quadratic inputs of the first progressive acoustic mode:

$$F_{1x}=\left(\begin{array}{ccccc}1\mu \frac{^2}{x^2}𝑑y𝑑y& \sqrt{\mu }\frac{}{x}𝑑y& \mu \frac{^2}{xz}𝑑y𝑑y& 0& 0\end{array}\right)(\phi _1+\phi _{1tv}),$$
$`(31)`$
where $`\phi _1,\phi _{1tv}`$ are supposed to consist of specific perturbations for the rightward dominant acoustic mode, see (10). For example, $`\phi _1`$ looks

$`\phi _1=\left(\begin{array}{c}(\stackrel{}{v}_1\stackrel{}{})v_{1x}+\sqrt{\mu }\rho _1p_1/x\\ (\stackrel{}{v}_1\stackrel{}{})v_{1y}+\rho _1p_1/y\\ (\stackrel{}{v}_1\stackrel{}{})v_{1z}+\sqrt{\mu }\rho _1p_1/z\\ [Qp_1+S\rho _1](\stackrel{}{}\stackrel{}{v}_1)(\stackrel{}{v}_1\stackrel{}{})p_1\\ \rho _1(\stackrel{}{}\stackrel{}{v}_1)(\stackrel{}{v}_1\stackrel{}{})\rho _1\end{array}\right).`$ There are all other possible cross nonlinear terms in the right-hand side of (30) which are out of interest in the present investigation since the rightward acoustic mode supposed to be dominant at least at the initial stage of the evolution. Writing on the terms of order not higher then $`\beta \mu `$ , one goes to the expected result:

$$\left(\begin{array}{ccccc}1\mu \frac{^2}{x^2}𝑑y𝑑y& \sqrt{\mu }\frac{}{x}𝑑y& \mu \frac{^2}{xz}𝑑y𝑑y& 0& 0\end{array}\right)\phi _1=0,$$
$`(32)`$
that means that acoustical streaming may exist only in the thermoviscous flow. Finally, $`F_{1x}`$ is as follows

$`F_{1x}=\left(\begin{array}{ccccc}1\mu \frac{^2}{x^2}𝑑y𝑑y& \sqrt{\mu }\frac{}{x}𝑑y& \mu \frac{^2}{xz}𝑑y𝑑y& 0& 0\end{array}\right)\phi _{1tv}=`$

$`\sqrt{\mu }\left[\frac{\beta }{2}\left(\frac{\rho _1}{x}\frac{\rho _1}{y}\frac{}{x}𝑑y(\rho _1\frac{^2\rho _1}{y^2})\frac{}{x}𝑑y(\frac{\rho _1}{y})^2\right)+(\delta _2^1+\delta _2^2)\left(\frac{}{x}𝑑y(\rho _1\frac{^2\rho _1}{y^2})\rho _1\frac{^2\rho _1}{xy}\right)\right]+`$

$$\sqrt{\mu }\left[\beta \left(\frac{}{x}𝑑y\left(\rho _1\frac{^2\rho _1}{y^2}\right)\rho _1\frac{^2\rho _1}{xy}\right)+(\delta _1^1+\delta _1^2)\left(\frac{}{x}𝑑y(\rho _1\frac{^2\rho _1}{y^2})\rho _1\frac{^2\rho _1}{xy}\right)\right].$$
$`(33)`$
To show that both effects play role: (1) presence of the thermoviscous terms in $`\psi _1,`$ and (2) presence of thermoviscous terms in $`\phi _{1tv},`$ corresponding terms are included in square brackets, first and second ones though they partially compensate each other. Finally, (33) goes to the following:

$$F_{1x}=\sqrt{\mu }\beta \left(\frac{\rho _1}{x}\frac{\rho _1}{y}0.5\rho _1\frac{^2\rho _1}{xy}\frac{}{x}𝑑y(\frac{\rho _1}{y})^2\right),$$
$`(34)`$
and the constant of integration of the last term should be chosen accordingly to the real initial perturbations and /or boundary regime. It is remarkable that only three first elements of $`\phi _{1tv}`$ participate in the final radiation force, namely that three ones are responsible for the losses of momentum in the flow. The forth element expresses the thermoviscous loss of energy. This term (and also the fifth one) plays role in the forming of the heating, see the last line of $`\stackrel{~}{P}_3`$ given by (20). Relations of the specific perturbations of eigenvector $`\psi _1`$ (14) are used to express $`F_{1x}`$ through $`\rho _1`$.

The transverse component of the force caused by sound $`F_{1x}`$ given by formula (34) does not possesses terms depending on the other transverse coordinate z because these terms are of order $`\beta \mu ^{3/2}`$ . In all calculations, the properties of specific modes and related projectors applied only. The temporal averaging was not proceeded. It has been proved that the acoustic force should be rotational . In the referred papers, the rotational part of the averaged is just selected from the general expression. Projectors $`P_{4,}P_5`$ due to their properties yield in the vortical forces automatically.

### 5.1 Limit of the quasi-periodic acoustic source

To compare formula (30) with this caused by periodic at the tranducer source, let us take the transverse acoustic force (radiation force) given by Gusev, Rudenko

$$\mathrm{\Phi }_{1x}=\sqrt{\mu }\rho _0^2\frac{}{x}(\theta ^2/2)\mathrm{exp}(\beta y)$$
$`(35)`$
rewritten in our variables; note also that longitudinal coordinate there is $`x`$ instead of $`y`$. The radiation force (35) relates to the streaming evolving in the sound field of a highly attenuated beam with a quasi-periodic source as follows:

$$\rho _1(x,y,t)=\rho _0\theta (x)\mathrm{exp}(\frac{\beta }{2}y)\mathrm{sin}(ty)$$
$`(36)`$
see  for more details. Calculating based on formula (34) yields in the following:

$$F_{1x}=\sqrt{\mu }\rho _0^2\frac{}{x}(\theta ^2/2)\mathrm{exp}(\beta y)\left[1+\frac{\beta }{4}\mathrm{sin}(2(ty))+O(\beta ^2)\right]$$
$`(37)`$
Averaging over integral number of period of the sound wave in the leading order gives $`F_{1x}=\mathrm{\Phi }_{1x}.`$ In contrast to the known results, (34) fits to calculate force caused by any acoustic source including non-periodic ones. Moreover, evolution of streaming given by temporal averaging is not suitable in principle for detail tracing with small temporal step because of the very initial procedure of averaging over time whose interval is much longer than a period of a sound wave.

### 5.2 Radiation force cased by mono-polar acoustic source

As an appropriate acoustic source satisfying the limitations of , let us take a mono-polar two-dimensional wave as follows :

$$\rho _1(x,y,t)=\sqrt{\frac{2\beta }{\pi }}\mathrm{exp}(x^2)\frac{\mathrm{exp}\left(\tau ^2/2\xi \right)}{ϵ\sqrt{\xi /\beta }(CErf(\tau /\sqrt{2\xi })},$$
$`(38)`$
where $`\xi =\beta y`$, $`\tau =ty,`$ $`ϵ=\frac{QS+1}{2}`$ is a parameter of nonlinearity, see (5). For ideal gases, $`ϵ=\frac{\gamma +1}{2}.`$ Constant C is responsible for the shape of the curve: large C provides a curve close to the Gauss one. The function $`F_{1x}`$ achieves negative minimum at the transversal distance $`x_m=\sqrt{2}/2`$ and is an odd function as follows from (34), (38). Results of the calculating of $`\frac{F_{1x}}{\sqrt{\mu }\beta }`$ in accordance to formula (34) are shown at the figures 1,2 at the transversal point $`x=x_m`$. Figure 1 shows the dependence of $`\frac{F_{1x}}{\sqrt{\mu }\beta }`$ on $`y`$ at $`t=1`$, and figure 2 shows the same at $`t=3.`$ Constants of mono-polar source (35) used in calculations are: $`C=2,ϵ=1.2,\beta =0.1.`$ Also, the constant of numerical integration of the last term of (31) is chosen to the radiation force be zero at $`y>\mathrm{}.`$

All pictures show that streaming develop with some delay after the source passing. The reason is non-local relation between the source and induced streaming.The velocity field may be calculated accordingly to the formula (30) completed with equations for the other components of $`\stackrel{}{V}`$. At the first stage of evolution, the nonlinear term $`(\stackrel{}{V}\stackrel{}{})\stackrel{}{V}`$ may be neglected and the linear equation like this of thermal conductivity with viscous coefficient cased by shear viscosity $`\delta _1^2`$ should be solved. The acoustic source in the right-hand side of (30) is namely the radiation force. The temporal behavior of the radiation force at the two longitudinal coordinates $`y=1`$ and $`y=3`$ is presented at the figures 3,4, $`x=\sqrt{2}/2.`$ The mono-polar source is naturally attenuated during its propagation along y-axis. Note that the radiation force is multiplied by the large value $`\frac{1}{\sqrt{\mu }\beta }`$ at all pictures, so the real scale is much less.

Fig.1 Longitudinal distribution of the radiation force $`\frac{F_{1x}}{\sqrt{\mu }\beta }(y)`$ (thin line) and the acoustic source $`\rho _1(y)`$ (bold line) at $`t=1,x=\sqrt{2}/2.`$

Fig.2 Longitudinal distribution of the radiation force $`\frac{F_{1x}}{\sqrt{\mu }\beta }(y)`$ (thin line) and the acoustic source $`\rho _1(y)`$ (bold line) at $`t=3,x=\sqrt{2}/2.`$

Fig.3 The temporal distribution of the radiation force $`\frac{F_{1x}}{\sqrt{\mu }\beta }(t)`$ (thin line) and the acoustic source $`\rho _1(t)`$ (bold line) at $`y=1,x=\sqrt{2}/2.`$

Fig.4 The temporal distribution of the radiation force $`\frac{F_{1x}}{\sqrt{\mu }\beta }(t)`$ (thin line) and the acoustic source $`\rho _1(t)`$ (bold line) at $`y=3,x=\sqrt{2}/2.`$

## 6 Conclusions

The basic idea of the paper is to separate modes accordingly to their properties in the weakly nonlinear flow. At first, modes as eigenvectors of the linear flow should be defined. In the other words, the relations of the specific perturbations inside every mode should be established. Both homogeneous and inhomogeneous backgrounds ( see paper on interacting modes in bubbly liquid ), media affected by external forces including the gravitational one which changes the background density and pressure  are treated in the algorithmic way. The definition of modes is unique, determined by the linearized differential conservation equations only. Proper initial or/and boundary conditions are given by the proper superposition of modes.

The second step is to get coupled nonlinear equations for the interacting modes and to solve it approximately. If one of the acoustic modes is dominant the generation of the entropy mode and the other acoustic mode in the plane geometry (namely heating and the reflected wave) is governed by the simple analytical formulae. The evolution equations may be corrected up to the higher order nonlinear terms due to increasing influence of the other generated modes .

The principal advance in comparison to the widely used approaches consists in possibility to treat any sources: initial mixture of all possible types of motion, non-periodic acoustic source and so on. In the present paper, radiation force caused streaming is calculated for any acoustic source, an example of the monopolar pulse source is considered, the limit of periodic source is traced.

REFERENCES

 Zarembo, L.K.(1971). Acoustic streaming. Rozenberg, L.D., editor, High-intensity ultrasonic fields.Part 3, Plenum Press, NewYork,137-191.

 Rudenko, O.V., and Soluyan, S.I. (1977). Theoretical foundations of nonlinear acoustics. Consultants Bureau, New York, English edition.

 Nyborg, W.L. (1997). Acoustic streaming. Hamilton M.F., D.T.Blackstock, D.T., editors, Nonlinear Acoustics, Academic press, New York, 207-231.

 Tjotta,S, and Tjotta, J.N. (1993). Acoustic streaming in ultrasound beams.Proc.13 Int. Symp. Nonlinear Acoustics, 601-606.

 Qi,Q.(1993).The effect of compressibility on acoustic streaming near a rigid boundary for a plane travellimg wave. J. Acoust. Soc. Am. 94, 1090-1098.

 Menguy,L.,Gilbert, J.(2000). Non-linear acoustic streaming accompanying a plane stationary wave in a guide, Acta Acustica, 86, 249-259.

 Kamakura,T., Kazuhisa M., Kumamoto. Y., and Breazeale, M.A. (1995). Acoustic streaming induced in focused Gaussian beams, J.Acoust. Soc. Am. 97(5), Pt. 1, 2740-2746.

 Makarov,S., and Ochmann, M.(1996). Nonlinear and thermoviscous phenomena in acoustics, part I., Acustica, 82, 579-606.

 Chu, B.-T., and Kovasznay, L.S.G. (1958) Nonlinear interactions in a viscous heat-conducting compresible gas, Journ. Fluid. Mech.,3, 494-514.

 Perelomova, A.A.(2000) Projectors in nonlinear evolution problem: acoustic solitons of bubbly liquid. Applied Mathematical Letters 13, 93-98.

 Perelomova, A.A.(1998) Nonlinear dynamics of vertically propagating acoustic waves in stratified atmosphere, Acta Acustica, 84 , 1002-1006.

 Perelomova, A.A. (2001) Directed acoustic beams interacting with heat mode: coupled nonlinear equations and modified KZK equation, Acta Acustica, 87, 176-183.

 Leble, S.B. (1991) Nonlinear waves in waveguides with stratification, Springer-Verlag, Berlin.

 Zabolotskaya,E.A, Khokhlov, R.V , Quasiplane waves in the nonlinear acoustic field of confined beams, Sov.Phys.Acoust., 15(1969),35-40

 Bakhvalov, N.S.,Zhileikin,Ya.M and Zabolotskaya E.A, Nonlinear theory of sound beams (American Institute of Physics, New York, 1987).

 Gusev, V.E., and Rudenko, O.V. (1979) Non-steady quasi-one-dimensional acoustic streaming in unbounded volumes with hydrodynamic nonlinearity, Sov. Phys. Acoust. 25,493-497.

 Gusev, VE., and Rudenko, O.V. (1980) Evolution of nonlinear two-dimensional acoustic streaming in the field of a highly attenuated sound beam, Sov. Phys. Acoust. 27(6), 481-484.

 Rudenko, O.V., Soluyan S.I. (1977) Theoretical foundations of nonlinear acoustics, Plenum, New York.