text:
# Negative Magnetoresistance in Andreev Interferometers

## Abstract

We consider transport in a diffusive cross-shaped Andreev interferometer geometry, both theoretically and experimentally. A strong and unexpected modulation of the conductance with the superconducting phase is found. In particular, a reversed resistance vs. phase difference is predicted, where the resistance is decreased by a phase gradient. A comparison of our quasiclassical calculation with experimental data shows quantitative agreement if we account for diamagnetic screening.

Quantum mechanical interference effects can be ideally studied in heterostructures of normal metals and superconductors. Large conductance oscillations in structures containing loops threaded by a magnetic flux have been predicted and observed petrashov1 ; pothier ; petrashov2 ; vanwees ; courtois ; volkov ; stoof:96 ; GWZ . The origin of these oscillation is the phase coherence between Andreev coupled electron-hole pairs, which can be maintained over large distances of the order of $`\mu `$m, making their observation in microstructured devices accessible VZK ; raimondi ; suplattrev ; pannetier . Since the origin of these oscillation are superconducting correlations, they are usually suppressed by a magnetic field, i.e. these devices show a positive magnetoresistance. This trend can be reversed by means of tunnel junctions, in which case the conductance is strongly reduced due to a suppression of the low energy density of states leads. Obviously, the suppression of superconducting correlations then leads to a decreased resistance antonov . Thus, one would expect that a magnetic field drives the resistance towards its normal state value, no matter in what direction this actually is.

Here we report on a theoretical calculation and an experimental measurement of the magnetotransport of a purely diffusive heterostructure, which shows a negative magnetoresistance due to an subtle interference effect. This is surprising, since generally one would expect, that both finite energies and a phase gradient would lead to a destructive interference, diminishing the proximity effect. Thus, although the conductance shows the so-called reentrance effect as function of temperature, one would not expect a phase gradient to reverse the temperature effect. In the present Letter we study a geometry, in which a separation of the energy scales responsible for these two effects result in a negative magnetoresistance.

We will first describe the theoretical approach, which predicts a negative magnetoresistance for small magnetic flux. To obtain the experimentally observed fully inverted resistance oscillations, we will have to include diamagnetic screening in the superconducting loop imposing the phase difference. This allows us to obtain excellent agreement using parameters determined from the experiment independently. We note that this agreement is contrary to the claim in Ref. shaikaidarov1 , in which a similar experimental layout was studied. It was however claimed, that the experimental observations can not be explained within the ’standard theory’ of the mesoscopic diffusive proximity effect. Here, we show that these observations are in perfect agreement with the ’standard’ quasiclassical theory.

The system we study is shown in Fig. 1. The left part shows the experimental setup and the right shows the theoretical model structure. The resistance of the vertical piece of normal metal is to be measured. Each part (left and right) of this vertical wire has a resistance $`R_N`$ and a length $`L_N`$. The reservoirs inducing superconductivity (S1 and S2) are attached to the middle of the wire by normal resistors of length $`L_S`$ and resistance $`R_S`$. The actual measurement of the resistance is performed in a 4-terminal geometry by attaching current and voltage probes (N1-N4) through resistors of length $`L_I`$ and resistance $`R_I`$. Note, that these probes have to be included into the modeling in contrast to normal systems. In a proximity effect structure, the properties of these probes matter, since they modify the equilibrium properties (i.e. the density of states and therefore the energy- and space-dependent conductivity of the system).

The calculation of the resistance of the structure depicted in Fig.1 has been performed numerically, along the lines of Ref. nazarov:99 ; stoof:96  using the quasiclassical formalismquasicl . From the spectral Usadel equation the local energy-dependent conductivity is extracted via

$$\sigma (E,x)=\frac{\sigma _N}{4}\text{Tr}\left[\widehat{1}\widehat{\tau }_z\widehat{g}_R\widehat{\tau }_z\widehat{g}_A\right].$$
(1)
Here $`\sigma _N`$ is the normal state conductivity and $`\widehat{g}_{R,A}`$ are retarded and advanced Nambu Green’s functions obeying the Usadel equation (see Ref. nazarov:99  for a definition). The solution of the kinetic equation yields for the measured 4-terminal resistance

$$R^1=\frac{dE}{4T\mathrm{cosh}^2\left(\frac{E}{2T}\right)}\left(_0^{L_N}\frac{dx}{\sigma (E,x)}\right)^1.$$
(2)
Similarly we can find the supercurrent in equilibrium as

$$I_S=\frac{\sigma }{8}𝑑E\text{Im}\text{Tr}\left[\widehat{\tau }_z\widehat{g}_R\frac{}{y}\widehat{g}_R\right]\mathrm{tanh}\left(\frac{E}{2T}\right).$$
(3)
These expressions show, that the 4-Terminal configuration indeed measures the resistance of the vertical branch only. However, due to the proximity effect the spectral conductivity $`\sigma (E,x)`$ depends also on the properties of the current and voltage leads. As a consequence this will also lead to a strongly increased supercurrent for temperatures at which these branches are probed, as we demonstrate below.

We now discuss the concrete theoretical results. The parameters, that we have chosen correspond to the experimental setup shown in the left panel of Fig. 1. The dimensions are $`L_N=1000`$nm, $`L_S=225`$nm and $`L_I=800`$nm, where $`L_I`$ is taken a little bit larger than in the experiment to account for the proximity effect into the wider (”lead”) regions of the experiment. For the connectors to the superconductor we have taken the distance between the central point and the closest point of the superconductor. Taking the experimental value of the diffusion constant $`D80`$cm$`{}_{}{}^{2}/`$s, we find for the Thouless energy $`E_N=\mathrm{}D/L_N^2=57`$mK$`k_B`$. The branch connecting the two superconductors directly has correspondingly $`E_S=\mathrm{}D/4L_S^2=280`$mK$`k_B`$ (we took the full length, since this is mainly used in literature to characterize SNS junctions). In the calculation we have assumed a temperature independent superconducting gap corresponding to $`\mathrm{\Delta }/k_B=2`$K.

Let us now turn to the phase dependent resistance of the vertical branch, measured in 4-point configuration through the leads N1-N4. The resistance oscillations for a large range of temperatures are shown in Fig. 2. Following the temperature dependence of the resistance for zero phase difference, we observe the usual reentrance behaviour, with a maximal suppression of about 13% at a temperature around 200mK. At a phase difference of $`\pi `$ (or any odd multiple), the proximity effect is completely suppressed for all temperatures. At temperatures above 200mK the resistance vs. phase oscillation are similar to that predicted previously volkov ; stoof:96 ; GWZ  and found experimentallycourtois ; petrashov1 . However, below 200mK the oscillations change qualitatively. Whereas the resistance increases monotonically starting from zero phase difference for the larger temperatures, it decreases first for the lower temperatures. Only close to half-integer phases the resistance increases again. At half-integer flux the resistance is always equal to the normal state value, as it should be.

Obviously, such a strong change of the phase-dependence should also affect other transport properties, for example the temperature- and the phase-dependence of the supercurrent. The theoretical results are shown in the inset of the left panel of Fig. 3. Note that, for simplicity, we have not included the temperature dependence of the superconducting gap. Therefore, the supercurrent is finite above the experimental critical temperature of $`1.4K`$. The critical current strongly increases below a temperature close to $`E_S`$. At higher temperatures the critical current depends only weakly on the temperature (the exponential decrease $`\mathrm{exp}(\sqrt{k_BT/E_S})`$ sets in only at much larger temperatures). The strong increase is related to the influence of the vertical branch. As superconducting correlation become strong in this branch, the pair breaking influence of normal lead gets effectively suppressed. Accordingly the supercurrent increases drastically. The phase-difference, for which the critical current is reached, shifts from $`0.65\pi `$ up to $`0.8\pi `$ in the same temperature interval. In the left panel of Fig. 3 we show the current-phase relation in an temperature interval between 100mK and 2K. It is clearly non-sinusoidal and also changes qualitatively at low temperatures.

Summarizing these theoretical results, we have observed a rather surprising phase-dependence of transport properties in a multi-terminal diffusive heterostructure. However, the experimental results, shown in the left column of Fig. 4, are quite different. The reason, as we will quantify below, is the diamagnetic screening in the superconducting loop. The applied magnetic flux induces a supercurrent in the loop, which screens the magnetic flux. Therefore, the phase difference between the superconducting reservoirs is not directly given by the applied magnetic flux. This effect depends on the self-inductance of loop. We will show below, that accounting of this leads to a fully inverted phase dependence and a very good agreement of experimental results and theoretical calculations.

To address the screening effect we consider a superconducting loop interrupted by a weak link. For the system depicted in Fig. 1 the weak link is the normal metal structure. We therefore have to assume that the entire phase-drop occurs over the normal metal, which is equivalent to say that the critical current in the normal metal (induced by proximity) is smaller than the critical current in the superconducting ring connecting the two reservoirs S1 and S2. The energy of a superconducting loop, containing a weak link reads,

$$E(\phi )=\frac{1}{2}LI(\phi )^2+E_\text{J}(\phi ).$$
(4)
Here, $`L`$ is the selfinductance of the loop, which depends on the geometry, and $`E_\text{J}(\phi )`$ is the Josephson energy of the weak link. It has to be determined from the integration of the current-phase relation $`I(\phi )`$. Fluxoid quantization requires that $`\phi =2\pi \mathrm{\Phi }_{\text{tot}}/\mathrm{\Phi }_0`$ (we restrict ourselves here to $`\phi [0,\pi ]`$, for simplicity). The total flux $`\mathrm{\Phi }_{\text{tot}}`$ is the sum of the externally applied flux $`\mathrm{\Phi }_\text{x}`$ and the flux created by the supercurrent $`\mathrm{\Phi }_{\text{ind}}=LI(\phi )`$. Introducing $`\phi _x=2\pi \mathrm{\Phi }_x/\mathrm{\Phi }_0`$ and minimizing (4) with respect to $`\phi `$ we obtain an equation, which determines the phase difference as function of the external flux

$$\phi _x\phi =\frac{2e}{\mathrm{}}LI(\phi )\gamma I(\phi ).$$
(5)
In the limit of vanishing selfinductance ($`L0`$) the phase difference across the junction is equal to the externally applied flux. This is usually the desired result. On the other hand, if the left hand side of Eq. (5) is not negligible, the phase at the junction differs from the applied flux. If e.g. the current is $`I(\phi )=I_c\mathrm{sin}\phi `$ and $`(2e/\mathrm{})LI_c1`$, we find $`\phi \phi _x/(2eLI_c/\mathrm{})\phi _x`$. The phase is always much smaller than the applied flux. This holds as long as the applied flux is less than $`\pi `$. If the external flux exceeds this value, other solutions become important and the phase jumps. The main consequence of these jumps is, that certain intervals of phases around odd multiples of $`\pi `$ can not be reached by modulation of the external flux. Below, we will fully account for the screening, by calculating the full current phase-relation for each temperature. Then the solution of Eq. (5) determines the actual phase difference.

As discussed previously accounting for the screening effect renders a certain interval of phase differences unobservable. Thus, we expect that screening suppresses the resistance at half integer external flux, since the actual phase difference differs from an odd multiple of $`\pi `$. The screening effect is directly seen in the temperature dependence of the resistance for half-integer flux. This is depicted in the right panel Fig. 3 together with the zero-flux resistance. The value of the selfinductance parameter $`\gamma `$ is chosen to match the experimentally achieved minimal resistance around 100mK. Note, that without screening the resistance at half-integer flux would be always equal to the normal state resistance. Taking screening into account strongly suppresses the resistance below 400mK with a maximal suppression of $``$7%. At lower temperatures, the resistance reenters, similar as the zero-flux resistance, back to the normal states resistance at zero temperature. As can be seen the resistance for half-integer flux is also slightly below the zero-flux resistance, indicating a reversal of the resistance vs. flux oscillations. This is in agreement with the experimental data. The deviations at high temperatures can be attributed to the neglection of the temperature dependence of the superconducting gap.

To confirm, that the theoretical results are fully consistent with the experimental data, we compare the resistance oscillations for various temperatures in Fig. 4 using the screening parameter, determined from the fit in Fig. 3. The similarity to the experimental curves is quite striking, although not all details coincide. However, it is clearly seen, that the screening effect can lead to the inverted resistance oscillations (”$`\pi `$-shift”) at low temperatures, i. e. the system has a negative magnetoresistance.

Let us finally comment on the numerical value of the screening parameter, $`\gamma =0.005e/G_NE_N`$. Using the experimental values $`R_N3\mathrm{\Omega }`$ and $`E_N=57`$mK$`k_B`$, we calculate for the selfinductance of the loop $`L2`$pH. On the other hand, we find for a circular loop of radius $`R=2\mu `$m and cross section $`W^2=20\times 100`$nm2 a self-inductance $`L_{\text{loop}}=\mu _0R(\mathrm{ln}(8R/W)7/4)=9`$pH. Although the selfinductance of a loop with the experimental geometry is not easy to determine, this rough estimate agrees very well with the selfinductance determined from our fit.

In conclusion we have shown, that a diffusive Andreev interferometer can show a negative magnetoresistance and a non-sinusoidal supercurrent-phase relation. This effect has been calculated theoretically and found experimentally. We obtained excellent agreement with experimental data, if the diamagnetic screening in the superconducting loop imposing the phase difference is taken into account. This eventually leads to a fully inverted (’$`\pi `$-shifted’) resistance oscillations.

###### Acknowledgements.

W. B. was supported by the Swiss NSF and the NCCR Nanoscience.


format: text

text:
# hep-ph/0203177 Transport coefficients, spectral functions and the lattice

## 1 Introduction

In recent years there has been a growing interest in the calculation of transport coefficients in quantum field theories at finite temperature. With the advent of relativistic heavy-ion colliders, such as RHIC, a proper knowledge of transport coefficients has become relevant since hydrodynamical descriptions of heavy-ion collisions provide a useful tool to analyse the experimental data . In practice the extension of ideal relativistic hydrodynamics to include finite transport coefficients is far from straightforward and applications to heavy-ion physics have only just begun .

If the temperature is sufficiently high and the theory is weakly coupled, transport coefficients can be computed in a perturbative expansion, employing either kinetic theory or field theory using Kubo formulas. It turns out that the latter approach requires the summation of an infinite series of higher-order diagrams, known as ladder diagrams, which has been a serious drawback for its use. For a scalar theory, the higher-order contributions in the loop expansion have been identified and summed in Ref. , using an intricate diagrammatic analysis, and the leading-order results for the shear and bulk viscosities have been found. In Ref.  the equivalence of an effective Boltzmann equation and the field-theoretical calculation is shown. Using a more transparent analysis, the diagrammatic conclusions of Ref.  have been confirmed recently . So far, the only other transport coefficient for which the ladder series has been summed explicitly is the color conductivity in QCD . The viscosities in the scalar theory and color conductivity in QCD have the property that the one-loop contribution and the ladder contributions are of the same order in the coupling constant. However, for other transport coefficients in gauge theories, such as the shear viscosity or the electrical conductivity, the ladder contributions are in fact larger than the one-loop one . Only recently and using kinetic theory, a complete computation of the leading logarithmic order of these transport coefficients has appeared . Unfortunately, a full leading-order computation is still lacking. Ref.  provides a useful guide to the literature.

Euclidean lattice simulations offer in principle the possibility to compute transport coefficients completely nonperturbatively . However, transport coefficients are determined by the small frequency limit of zero-momentum spectral functions of appropriate composite operators (such as components of the energy-momentum tensor) and spectral functions or other real-time correlators are not readily available on a euclidean lattice, although recent progress has been made with the Maximal Entropy Method (MEM) .444Note that for classical field theories at finite temperature spectral functions can be computed nonperturbatively by numerical simulations directly in real time . For that reason it was proposed in Ref.  to introduce instead an ansatz for the spectral function and fit the result to the numerical data for the euclidean correlator, employing a standard dispersion relation between these two. This approach was pursued more recently in Refs. .

Motivated by these studies, our goal in this paper is to compute the spectral function relevant for the shear viscosity at high temperature in weakly-coupled scalar (Sec. 3) and nonabelian gauge theories (Sec. 4). In the Conclusions we compare our findings with the analysis carried out so far in Refs. . It is found that the ansatz used in these papers is inadequate and we suggest a better one. We also point out a potential problem in the calculation of spectral functions at very low frequencies ($`\omega 0`$) from euclidean lattice correlators using the MEM approach.

## 2 Correlation functions

We start with a summary of basic relations between transport coefficients, spectral functions and euclidean correlators, using the shear viscosity as an example. The relations presented in this section are quite general and valid for arbitrary transport coefficients.

The shear viscosity can be defined from a Kubo relation as

$$\eta =\frac{1}{20}\underset{\omega 0}{lim}\frac{1}{w}d^4xe^{i\omega t}[\pi _{kl}(t,𝐱),\pi _{kl}(0,\mathrm{𝟎})],$$
(1)
with $`\pi _{kl}`$ the traceless part of the spatial energy-momentum tensor. The brackets denote the equilibrium expectation value at temperature $`T`$. We define spectral functions of hermitian (composite) operators, such as $`\pi _{kl}`$, as the expectation value of the commutator,

$$\rho _{\pi \pi }(xy)=[\pi _{kl}(x),\pi _{kl}(y)],$$
(2)
and in momentum-space

$$\rho _{\pi \pi }(\omega ,𝐩)=d^4xe^{i\omega ti𝐩𝐱}\rho _{\pi \pi }(t,𝐱).$$
(3)
Spectral functions obey basic symmetry relations555 When most manipulations take place in real-space, it can be convenient to define spectral functions such that they are real instead of purely imaginary in real-space. Here we stick to the usual convention and spectral functions are real in momentum-space. $`\rho _{\pi \pi }^{}(x)=\rho _{\pi \pi }(x)=\rho _{\pi \pi }(x)`$ and $`\rho _{\pi \pi }^{}(\omega ,𝐩)=\rho _{\pi \pi }(\omega ,𝐩)=\rho _{\pi \pi }(\omega ,𝐩)`$ as well as the positivity condition $`\omega \rho _{\pi \pi }(\omega ,𝐩)0`$. The shear viscosity is determined by the slope at zero frequency:

$$\eta =\frac{1}{20}\frac{d}{d\omega }\rho _{\pi \pi }(\omega )|_{\omega =0},$$
(4)
where $`\rho _{\pi \pi }(\omega )=\rho _{\pi \pi }(\omega ,\mathrm{𝟎})`$. Since the spectral function is odd, we consider from now on positive $`\omega `$ only.

The euclidean correlator (at zero spatial momentum) is given by

$$G_{\pi \pi }^E(\tau )=d^3x\pi _{kl}(\tau ,𝐱)\pi _{kl}(0,\mathrm{𝟎})_E(\tau =it),$$
(5)
where the imaginary time $`\tau [0,1/T]`$ and $`G_{\pi \pi }^E(1/T\tau )=G_{\pi \pi }^E(\tau )`$. The euclidean correlator and the spectral function are related via an integral equation, originating from a dispersion relation,

$$G_{\pi \pi }^E(\tau )=_0^{\mathrm{}}\frac{d\omega }{2\pi }K(\tau ,\omega )\rho _{\pi \pi }(\omega ),$$
(6)
with the kernel

$$K(\tau ,\omega )=e^{\omega \tau }n(\omega )+e^{\omega \tau }\left[1+n(\omega )\right]=e^{\omega \tau }+2n(\omega )\mathrm{cosh}\omega \tau ,$$
(7)
obeying $`K(\tau ,\omega )=K(\tau ,\omega )=K(1/T\tau ,\omega )`$. The Bose distribution is

$$n(\omega )=\frac{1}{\mathrm{exp}(\omega /T)1}.$$
(8)
The low-frequency part of the spectral function contains all information on the transport coefficient and its effect on the euclidean correlator can be estimated quite easily. When $`\omega T`$, the kernel can be expanded as

$$K(\tau ,\omega )=\frac{2T}{\omega }+\frac{\omega }{T}\left[\frac{1}{6}\tau T(1\tau T)\right]+𝒪(\omega ^3/T^3),$$
(9)
and all except the first term are suppressed. As a consequence we find that the contribution to the euclidean correlator from low frequencies,

$$G_{\pi \pi }^{E,\mathrm{low}}(\tau )=2T_0^{\omega _\mathrm{\Lambda }}\frac{d\omega }{2\pi }\frac{\rho _{\pi \pi }(\omega )}{\omega },$$
(10)
is independent of $`\tau `$. The frequency cutoff $`\omega _\mathrm{\Lambda }T`$ is introduced here to justify the expansion of the kernel. We conclude that the dominant effect of the low-frequency region is a constant $`\tau `$-independent contribution to the euclidean correlator.

## 3 Scalar field

We consider a one-component massless scalar field with a quartic $`\lambda \varphi ^4/4!`$ interaction.666We assume that the temperature is sufficiently high such that a possible zero-temperature mass $`m_0^2\lambda T^2`$ plays no role. The one-particle spectral function is

$$\rho (xy)=[\varphi (x),\varphi (y)]=G^>(xy)G^<(xy).$$
(11)
The one-particle Wightman functions, $`G^>(xy)=\varphi (x)\varphi (y)`$ and $`G^<(xy)=\varphi (y)\varphi (x)`$, are related to the one-particle spectral function via the KMS condition

$$G^>(p)=\left[n(p^0)+1\right]\rho (p),G^<(p)=n(p^0)\rho (p).$$
(12)
The traceless part of the spatial energy-momentum tensor reads

$$\pi _{kl}=_k\varphi _l\varphi \frac{1}{3}\delta _{kl}_m\varphi _m\varphi .$$
(13)
The lowest-order skeleton diagram that contributes to the spectral function for the shear viscosity in Eq. (2) follows from simple Wick contraction,

$$[\varphi ^2(x),\varphi ^2(y)]=2\left[G^{>2}(xy)G^{<2}(xy)\right],$$
(14)
and we find

$$\rho _{\pi \pi }(\omega )=\frac{4}{3}\frac{d^4k}{(2\pi )^4}(𝐤𝐤)^2n(k^0)\rho (k^0,𝐤)\left[\rho (k^0+\omega ,𝐤)\rho (k^0\omega ,𝐤)\right],$$
(15)
where we used the KMS conditions (12) and $`n(\omega )=n(\omega )+1`$. The four $`𝐤`$’s in the integrand arise from the derivatives in $`\pi _{kl}`$. In the remainder of this section we compute the one-loop spectral function (15) as a function of the external frequency.

The spectral function $`\rho _{\pi \pi }`$ depends on the one-particle spectral functions $`\rho `$, which contain the quasiparticle structure of the theory at finite temperature. We describe this in some detail since similar, though more complicated, considerations appear in gauge theories. For hard momenta $`|𝐤|T`$, excitations are on-shell with energy $`|𝐤|`$. For softer momenta screening effects become important and hard thermal loop (HTL) resummation yields a temperature-dependent plasmon mass, $`m^2=m_{\mathrm{th}}^2(13m_{\mathrm{th}}/\pi T+\mathrm{})`$ with $`m_{\mathrm{th}}^2=\lambda T^2/24`$ . Finally, collisions in the plasma result in a finite (but narrow) momentum-dependent width $`2\gamma _𝐤mT`$, changing the one-particle spectral function from an on-shell delta function to a Breit-Wigner spectral function (see below) .

We may now discuss the one-loop expression (15). First we note that the integral is dominated by hard $`T`$ momenta. Therefore, for external frequencies that are not too small a simple on-shell delta function for the one-particle spectral functions suffices to find the dominant contribution. We will refer to this region as the high-frequency region. For smaller frequencies, however, the arguments of the delta functions come close, producing a so-called pinch singularity . The pinch singularity is screened by a finite external frequency or width, whichever one is the largest. For very small external frequencies, the inclusion of the width777And of ladder diagrams, see below. is essential . This situation is sketched in Fig. 1. The effect of these nearly pinching poles is to enhance the spectral function compared to naive estimates. We will refer to the region where pinching poles are important as the low-frequency domain.

We start with the region where the frequency is much larger than the thermal width and no pinch-singularity problems are encountered. The one-particle spectral function can be taken on-shell and is

$$\rho _0(k^0,𝐤)=2\pi ϵ(k^0)\delta (k_0^2\omega _𝐤^2),$$
(16)
where $`\omega _𝐤=\sqrt{𝐤^2+m^2}`$ and $`ϵ(x)`$ the sign-function. HTL effects are included in the mass parameter. Evaluating the integrals in Eq. (15) with the use of the delta-functions (the angular integrals are trivial) results in (recall that we take $`\omega >0`$)

$$\rho _{\pi \pi }(\omega )=\theta (\omega 2m)\frac{\left(\omega ^24m^2\right)^{5/2}}{48\pi \omega }\left[n(\omega /2)+\frac{1}{2}\right],$$
(17)
which is shown in Fig. 2. The physical processes are the decay of a zero-momentum excitation with energy $`\omega `$ into two on-shell particles with equal and opposite momentum, and the inverse process of creation. The decay process contributes also at zero temperature and makes the spectral function increase as $`\omega ^4`$ at large frequencies. The threshold at $`\omega =2m`$ arises from (simple) HTL resummation. This concludes the analysis of the spectral function in the high-frequency domain.

We continue with the low-frequency region where pinch singularities lead to a nontrivial enhancement of the spectral function. We replace the on-shell one-particle spectral functions with Breit-Wigner spectral functions:

$$\rho _{BW}(k^0,𝐤)=\frac{1}{2\omega _𝐤}\left[\frac{2\gamma _𝐤}{(k^0\omega _𝐤)^2+\gamma _𝐤^2}\frac{2\gamma _𝐤}{(k^0+\omega _𝐤)^2+\gamma _𝐤^2}\right].$$
(18)
The width $`2\gamma _𝐤`$ is determined by the imaginary part of the retarded self energy

$$\gamma _𝐤=\frac{\text{Im}\mathrm{\Sigma }_R(\omega _𝐤,𝐤)}{2\omega _𝐤},$$
(19)
and the dominant contribution at weak coupling arises from two-to-two scattering from the two-loop setting-sun diagram . A convenient way to write this damping rate is as

$$\gamma _𝐤=\gamma \frac{T}{\omega _𝐤}B(|𝐤|/T;m/T),$$
(20)
where

$$\gamma =\frac{\lambda ^2T}{1536\pi },$$
(21)
determines the parametrical behaviour. The function $`B`$ contains the nontrivial momentum dependence and is related to $`A`$ defined in Ref.  as $`B(|𝐤|/T;m/T)=(6/\pi ^2)A(|𝐤|/T;m/T)`$. In the limit of hard momenta and small mass

$$\underset{m0}{lim}\underset{|𝐤|\mathrm{}}{lim}B(|𝐤|/T;m/T)=1.$$
(22)
For analytical estimates we will neglect the momentum dependence and take $`B=1`$. In the results obtained by numerical integration the full momentum dependence is included.

We insert the Breit-Wigner functions into expression (15) for the spectral function and perform the $`k^0`$ integral by integrating around the poles in the complex plane. We preserve only the dominant contributions and discard all terms suppressed by (powers of) the coupling constant with respect to the leading order contribution. Breit-Wigner spectral functions have four poles at complex energy-arguments $`k^0=\pm (\omega _𝐤\pm i\gamma _𝐤)`$. From the residue of these poles we keep $`n(\omega _𝐤\pm i\gamma _𝐤)n(\omega _𝐤)`$. The Bose distribution $`n(k^0)`$ has poles along the imaginary axis at $`k^0=2\pi inT`$, $`n`$. However, the residues at these poles are subdominant compared to those from the poles of the Breit-Wigner functions. Hence we do not include these contributions. After performing the $`k^0`$ integral we find

$$\rho _{\pi \pi }(\omega )=\frac{4}{3}_𝐤\frac{|𝐤|^4}{2\omega _𝐤}\left\{[n(\omega _𝐤)n(\omega _𝐤\omega )]I(\omega ,𝐤)\frac{\omega _𝐤\omega }{\omega ^2+4\gamma _𝐤^2}[\omega \omega ]\right\},$$
(23)
with

$$_𝐤=\frac{d^3k}{(2\pi )^3},$$
(24)
and

$$I(\omega ,𝐤)=\frac{8\gamma _𝐤}{(\omega 2\omega _𝐤)^2+4\gamma _𝐤^2}.$$
(25)
For sufficiently large $`\omega `$ and in the limit of small coupling (width) we may take $`I(\omega ,𝐤)4\pi \delta (\omega 2\omega _𝐤)`$ and the result of the previous calculation in the high-frequency region is recovered, as it should.

Let’s now consider the low-frequency region, $`\omega m`$. Here we may approximate $`I(\omega ,𝐤)2\gamma _𝐤/\omega _𝐤^2`$ and expand the difference between the Bose distributions to find

$$\rho _{\pi \pi }(\omega )=\frac{8}{3}_𝐤\frac{|𝐤|^4}{\omega _𝐤^2}n^{}(\omega _𝐤)\frac{\omega \gamma _𝐤}{\omega ^2+4\gamma _𝐤^2}(0\omega m).$$
(26)
We note that the last factor controls the pinch singularity in precise agreement with Fig. 1. Although in principle $`\rho _{\pi \pi }(\omega )/T^4`$ may depend on the three dimensionless combinations $`\omega /T`$, $`\gamma /T`$, and $`m/T`$, we find that in practice it only depends on $`\omega /\gamma `$ and $`m/T`$. Since the integral is dominated by hard momenta the $`m/T`$ dependence is subdominant and the natural parameter on which the spectral function depends in this region is $`\omega /\gamma `$. It is easy to see that the spectral function has a local maximum at $`\omega \gamma `$ and $`\rho _{\pi \pi }(\omega \gamma )/T^41`$.

For very small frequencies $`\omega \gamma `$ we expand and obtain

$$\rho _{\pi \pi }(\omega )/T^4=a_1\left(\frac{\omega }{\gamma }\right)+\frac{a_3}{3!}\left(\frac{\omega }{\gamma }\right)^3+\mathrm{}(0\omega \gamma ),$$
(27)
with

$`a_1={\displaystyle \frac{2}{3T^4}}{\displaystyle _𝐤}{\displaystyle \frac{|𝐤|^4}{\omega _𝐤^2}}n^{}(\omega _𝐤){\displaystyle \frac{\gamma }{\gamma _𝐤}}{\displaystyle \frac{5!\zeta (5)}{3\pi ^2}}`$ (28)
$`a_3={\displaystyle \frac{1}{T^4}}{\displaystyle _𝐤}{\displaystyle \frac{|𝐤|^4}{\omega _𝐤^2}}n^{}(\omega _𝐤)\left({\displaystyle \frac{\gamma }{\gamma _𝐤}}\right)^3{\displaystyle \frac{7!\zeta (7)}{2\pi ^2}},`$ (29)
where the $``$ indicates that the final integrals are evaluated by neglecting the remaining momentum dependence of the damping rate (i.e. taking $`B=1`$) as well as the thermal mass (the first approximation has numerically the largest effect). From this the one-loop viscosity follows as

$$\eta _{1\mathrm{loop}}=\frac{1}{30}_𝐤\frac{|𝐤|^4}{\omega _𝐤^2}n^{}(\omega _𝐤)\frac{1}{\gamma _𝐤}\frac{2\zeta (5)}{\pi ^2}\frac{T^4}{\gamma }.$$
(30)
However, as is well-known  these one-loop results are not complete and a ladder summation is required to obtain the complete leading-order result (see Fig. 3). Due to nearly pinching poles, each additional rung in the ladder contributes with a factor $`\lambda ^2T/\gamma 1`$ and is therefore not suppressed. The effect of ladder summation is to change the coefficients $`a_1,a_3,\mathrm{}`$, but not the parametric dependence on the coupling constant.

In the region $`\gamma \omega m`$ we find

$$\rho _{\pi \pi }(\omega )/T^4=\frac{8}{3T^4}_𝐤\frac{|𝐤|^4}{\omega _𝐤^2}n^{}(\omega _𝐤)\frac{\gamma _𝐤}{\omega }\frac{8\zeta (3)}{\pi ^2}\frac{\gamma }{\omega }(\gamma \omega m).$$
(31)
In this frequency interval pinching-pole contributions from ladders are subdominant since each additional rung comes with a factor $`\lambda ^2T/\omega 1`$. The perturbative part of the three-loop ladder diagram (i.e. with a single rung) also contributes at this order and has the same $`\lambda ^2T/\omega `$ behaviour as we find above.888We thank Guy Moore for pointing this out. Fig. 4 shows the contribution to the spectral function from the nearly pinching poles in the low-frequency interval, obtained by numerical integration of Eq. (26) with the full momentum and mass dependence. Note that the natural dimensionless combinations in the low-frequency region differ from those in the high-frequency region (compare Figs. 2 and 4).

We may now combine the results obtained so far to construct the complete one-loop spectral function at high temperature in the weak-coupling limit. The spectral function can be written as the sum of the contributions discussed above:

$$\rho _{\pi \pi }(\omega )=\rho _{\pi \pi }^{\mathrm{low}}(\omega )+\rho _{\pi \pi }^{\mathrm{high}}(\omega ),$$
(32)
where $`\rho _{\pi \pi }^{\mathrm{low}}`$ represents the contribution from the nearly pinching poles in Eq. (26), dominating at low frequencies, and $`\rho _{\pi \pi }^{\mathrm{high}}`$ is the contribution from decay/creation processes in Eq. (17), dominating at higher frequencies.

We find that for very small frequencies the spectral function rises quickly as $`\omega /\gamma `$. It reaches a maximum of order 1 (in units of temperature) at $`\omega \gamma `$ and decays then slowly as $`\gamma /\omega `$. Around the thermal mass the contribution from decay/creation processes enter and the spectral function increases again. Note that the results obtained in both frequency domains smoothly match parametrically at $`\omega m`$, since $`\rho _{\pi \pi }(\omega 3m)/T^4\lambda ^{3/2}`$, both from the low- and the high-frequency calculation. For large $`\omega `$ the spectral function increases as $`\omega ^4`$, due to the zero-temperature decay process. Ladder diagrams do not affect this characteristic shape. In Fig. 5 we present the complete one-loop spectral function as a function of $`\omega /T`$.999 We used here that $`m/T=0.1`$ corresponds to $`\lambda 0.267`$ and $`\gamma /T1.4810^5`$. In order to combine the low- and the high-frequency contribution in one figure, we show $`\rho _{\pi \pi }(\omega )/\omega ^4`$, which enhances the contribution at lower frequencies.

We are now ready to calculate the euclidean correlator using Eq. (6). Because $`G_{\pi \pi }^E(\tau )`$ depends linearly on the spectral function, we write it as a sum of two contributions, $`G_{\pi \pi }^E=G_{\pi \pi }^{E,\mathrm{low}}+G_{\pi \pi }^{E,\mathrm{high}}`$, and discuss each term separately. We start with the contribution due to decay/creation processes which reads

$$G_{\pi \pi }^{E,\mathrm{high}}(\tau )=_{2m}^{\mathrm{}}\frac{d\omega }{2\pi }K(\tau ,\omega )\frac{\left(\omega ^24m^2\right)^{5/2}}{48\pi \omega }\left[n(\omega /2)+\frac{1}{2}\right].$$
(33)
It is easy to see that the mass plays only a subdominant role and finite-mass corrections are suppressed by $`m^2/T^2`$. Therefore we take $`m=0`$ which yields

$$G_{\pi \pi }^{E,\mathrm{high}}(\tau )=\frac{T^5}{96\pi ^2}_0^{\mathrm{}}𝑑xx^4\left[e^{sx}+e^{(1s)x}\right]n(x)\left[n(x/2)+\frac{1}{2}\right],$$
(34)
where $`x=\omega /T`$ and $`s=\tau T`$. The remaining integral can be performed and we find

$$G_{\pi \pi }^{E,\mathrm{high}}(\tau )=\frac{\pi ^2T^5}{3\mathrm{sin}^5u}\left\{(\pi u)\left[11\mathrm{cos}u+\mathrm{cos}3u\right]+6\mathrm{sin}u+2\mathrm{sin}3u\right\},$$
(35)
where $`u=2\pi \tau T`$. An approximate but illuminating expression for the euclidean correlator can be obtained by noticing that the integral in Eq. (34) is dominated by hard frequencies, $`x1/s`$ with $`0<s<1`$, such that the Bose distributions may be approximated with Maxwellian ones, $`n(x)e^x`$. In that case we find

$$G_{\pi \pi }^{E,\mathrm{high}}(\tau )\frac{1}{8\pi ^2}\left[\frac{1}{\tau ^5}+\frac{1}{(1/T\tau )^5}+\frac{2}{(3/2T\tau )^5}+\frac{2}{(1/2T+\tau )^5}\right].$$
(36)
This approximate expression differs less than 2% from the exact result (35). The dominant $`1/\tau ^5`$ behaviour of the correlator arises from decay at zero temperature. Finite temperature is manifested mainly through the reflection symmetry, $`\tau 1/T\tau `$. At the central point $`\tau T=1/2`$, we find the contribution to the euclidean correlator from the decay/creation processes to be

$$G_{\pi \pi }^{E,\mathrm{high}}(\tau =1/2T)=\frac{4\pi ^2}{45}T^5\left[1\frac{25}{8\pi ^2}\frac{m^2}{T^2}+\mathrm{}\right].$$
(37)
The effect of a finite mass is to lower the mimimal value of the correlator.

The contribution to the euclidean correlator from the nearly pinching poles in the low-frequency region can be found easily from Eqs. (6) and (26) by interchanging frequency and momentum integrals. Since $`\rho _{\pi \pi }^{\mathrm{low}}`$ gives the dominant contribution to the spectral function up to frequencies of order $`m`$, whereas $`\rho _{\pi \pi }^{\mathrm{high}}`$ dominates for higher frequencies, we can use the expansion (9) for the kernel and obtain

$`G_{\pi \pi }^{E,\mathrm{low}}(\tau )`$ $``$ $`{\displaystyle \frac{8}{3}}{\displaystyle _𝐤}{\displaystyle \frac{|𝐤|^4}{\omega _𝐤^2}}n^{}(\omega _𝐤){\displaystyle _0^{\omega _\mathrm{\Lambda }}}{\displaystyle \frac{d\omega }{2\pi }}{\displaystyle \frac{2T}{\omega }}{\displaystyle \frac{\omega \gamma _𝐤}{\omega ^2+4\gamma _𝐤^2}}`$ (38)
$`=`$ $`{\displaystyle \frac{4}{3}}{\displaystyle \frac{T}{\pi }}{\displaystyle _𝐤}{\displaystyle \frac{|𝐤|^4}{\omega _𝐤^2}}n^{}(\omega _𝐤)\mathrm{arctan}\left({\displaystyle \frac{\omega _\mathrm{\Lambda }}{2\gamma _𝐤}}\right)`$
$``$ $`{\displaystyle \frac{4\pi ^2}{45}}T^5\left[1{\displaystyle \frac{25}{8\pi ^2}}{\displaystyle \frac{m^2}{T^2}}+\mathrm{}\right],`$
with $`\omega _\mathrm{\Lambda }m`$. The error that is introduced by expanding the kernel is negligible (see Fig. 6). We note that in our one-loop calculation this result is at leading order independent of the coupling constant. Similarly, while ladder diagrams determine the precise shape of the spectral function, the effect on the euclidean correlator $`G_{\pi \pi }^{E,\mathrm{low}}(\tau )`$ appears only in subleading corrections in the coupling constant, as can be seen with the kinetic approach . Therefore the low-frequency contribution to the euclidean correlator is constant and of order one (in the appropriate units) and insensitive to details of the ladder summation.

Combining the low- and high-frequency contribution to the euclidean correlator at high temperature and weak coupling we find, to leading order in the coupling constant,

$$G_{\pi \pi }^E(\tau )=\frac{\pi ^2T^5}{3\mathrm{sin}^5u}\left\{(\pi u)\left[11\mathrm{cos}u+\mathrm{cos}3u\right]+6\mathrm{sin}u+2\mathrm{sin}3u\right\}+\frac{4\pi ^2}{45}T^5,$$
(39)
with $`u=2\pi \tau T`$. Corrections due to the finite thermal mass are suppressed by $`m^2/T^2`$. The euclidean correlator is minimal at $`\tau T=1/2`$, and $`G_{\pi \pi }^E(1/2T)`$ receives contributions of the same order from both the high- and the low-frequency region (at leading order in the coupling they are equal). A comparison between the analytical result at $`m=0`$ and the full result obtained by numerical integration of Eq. (6) with Eq. (32) in the presence of a finite mass is shown in Fig. 6.

## 4 Nonabelian gauge fields

We leave the scalar case and consider nonabelian $`SU(N_c)`$ gauge theory. The traceless spatial part of the energy-momentum tensor is

$$\pi _{ij}=F_i^{a\mu }F_{j\mu }^a\frac{1}{3}\delta _{ij}F^{ak\mu }F_{k\mu }^a.$$
(40)
The coupling vertex between the operator $`\pi _{ij}`$ and two gluons with incoming momenta $`P,K`$ and indices $`(\mu ,a),(\nu ,b)`$ respectively can be read from (40) and we find101010 In order to arrive at the basic one-loop expression (44) below we use here the imaginary-time formalism with $`P=(p_4,𝐩)`$, the Matsubara frequency $`\omega _n=p_4=2\pi nT`$ ($`n)`$ and $`PK=p_4k_4+𝐩𝐤`$.

$`\mathrm{\Gamma }_{ij,\mu \nu }^{ab}(P,K)`$ $`=`$ $`\delta ^{ab}[\delta _{\mu \nu }(p_ik_j+p_jk_i{\displaystyle \frac{2}{3}}\delta _{ij}𝐩𝐤)`$ (41)
$`+PK(\delta _{i\mu }\delta _{j\nu }+\delta _{i\nu }\delta _{j\mu }{\displaystyle \frac{2}{3}}\delta _{ij}\delta _{k\mu }\delta _{k\nu })(p_iK_\mu \delta _{j\nu }+p_jK_\mu \delta _{i\nu }`$
$`+P_\nu k_i\delta _{j\mu }+P_\nu k_j\delta _{i\mu })+{\displaystyle \frac{2}{3}}\delta _{ij}(P_\nu k_k\delta _{k\mu }+p_kK_\mu \delta _{k\nu })].`$
In the nonabelian theory $`\pi _{ij}`$ also couples to three and four gluons, which leads to diagrams as depicted in Fig. 7. However, these contributions are suppressed by powers of the coupling constant and will not be considered further.

We use the Coloumb gauge in which the gluon propagator reads

$$D_{\mu \nu }(P)=𝒫_{\mu \nu }^T(\widehat{𝐩})\mathrm{\Delta }_T(P)+\delta _{4\mu }\delta _{4\nu }\mathrm{\Delta }_L(P),$$
(42)
with $`𝒫_{ij}^T(\widehat{𝐩})=\delta _{ij}\widehat{p}_i\widehat{p}_j`$, $`𝒫_{4\mu }^T=𝒫_{\mu 4}^T=0`$. The transverse and longitudinal components have the following spectral representations

$$\mathrm{\Delta }_T(P)=_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{2\pi }\frac{\rho _T(\omega ,𝐩)}{i\omega _n\omega },\mathrm{\Delta }_L(P)=\frac{1}{𝐩^2}+_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{2\pi }\frac{\rho _L(\omega ,𝐩)}{i\omega _n\omega }.$$
(43)
The one-loop contribution to the spectral function reads, after evaluating the Matsubara sum and taking all frequencies real again,

$`\rho _{\pi \pi }(\omega )=`$ $`{\displaystyle \frac{2d_A}{3}}{\displaystyle }{\displaystyle \frac{d^4k}{(2\pi )^4}}[n(k^0)n(k^0+\omega )]\{V_1(k,\omega )\rho _T(k^0,𝐤)\rho _T(k^0+\omega ,𝐤)`$ (44)
$`+V_2(k)\rho _T(k^0,𝐤)\rho _L(k^0+\omega ,𝐤)+V_3(𝐤)\rho _L(k^0,𝐤)\rho _L(k^0+\omega ,𝐤)\},`$
where $`d_A=N_c^21`$ is the number of gluons and $`V_1(k,\omega )=7|𝐤|^410𝐤^2k^0(k^0+\omega )+7(k^0)^2(k^0+\omega )^2`$, $`V_2(k)=6𝐤^2(k^0)^2`$ and $`V_3(𝐤)=32|𝐤|^4`$. This expression is the equivalent of Eq. (15) in the scalar theory.

As in the scalar case, we start the analysis for external frequencies $`\omega `$ that are sufficiently large such that no pinch singularities are present. The collective (HTL) effects in a nonabelian plasma  can be incorporated in the one-particle spectral functions, which are however more complex than in the scalar case. The spectral function for transverse gluons can be written as

$$\rho _T(k^0,𝐤)=2\pi Z_T(|𝐤|)\left\{\delta [k^0\omega _T(|𝐤|)]\delta [k^0+\omega _T(|𝐤|)]\right\}+\beta _T(k^0,|𝐤|).$$
(45)
The delta functions describe propagating quasiparticles with a dispersion relation $`\omega _T`$ and a residue

$$Z_T(k)=\frac{\omega _T(k)[\omega _T^2(k)k^2]}{3\omega _{\mathrm{pl}}^2\omega _T^2(k)[\omega _T^2(k)k^2]^2},$$
(46)
where $`\omega _{\mathrm{pl}}^2=g^2T^2N_c/9`$ is the plasma frequency squared. For small and large spatial momentum one finds

$`\omega _T^2(k)\omega _{\mathrm{pl}}^2+{\displaystyle \frac{6}{5}}k^2,Z_T(k)1/(2\omega _{\mathrm{pl}})(k0),`$ (47)
$`\omega _T^2(k)k^2+m_{\mathrm{}}^2,Z_T(k)1/(2k)(k\mathrm{}),`$ (48)
where $`m_{\mathrm{}}^2=\frac{3}{2}\omega _{\mathrm{pl}}^2`$ is the asymptotic gluon mass squared. The $`\beta _T`$ function describes Landau damping and is nonzero below the light-cone only. The spectral function of longitudinal gluons $`\rho _L(k^0,𝐤)`$ has a similar form .

We start studying the contribution to Eq. (44) when both gluons are transverse. There are three parts, depending on whether we take from the one-particle spectral functions the delta functions of the quasiparticles (pole-contribution) or the $`\beta `$ function of the Landau damping (cut-contribution). When both gluons are quasiparticles the integrals can be done with the help of the delta functions and we find the pole-pole contribution for transverse gluons

$$\rho _{\pi \pi }^{\mathrm{pp}}(\omega )=\theta (\omega 2\omega _{\mathrm{pl}})\frac{2d_A}{3\pi }\frac{Z_T^2(f_0)f_0^2}{\omega _T^{}(f_0)}\left(7f_0^4+\frac{5}{2}f_0^2\omega ^2+\frac{7}{16}\omega ^4\right)\left[n(\omega /2)+\frac{1}{2}\right].$$
(49)
Here we use the notation $`f_0=f(\omega /2)`$, with $`f(u)`$ defined as the inverse of the transverse dispersion relation, $`f[\omega _T(k)]=k`$. The function $`f(u)`$ vanishes when $`u\omega _{\mathrm{pl}}`$, and

$`f(u)\sqrt{{\displaystyle \frac{5}{6}}(u^2\omega _{\mathrm{pl}}^2)}(u\omega _{\mathrm{pl}}),`$ (50)
$`f(u)u{\displaystyle \frac{3\omega _{\mathrm{pl}}^2}{4u}}(u\mathrm{}).`$ (51)
Again, as for the scalar case, the HTL resummation produces the threshold in the spectral function (49) for soft external frequency. We find that

$$\rho _{\pi \pi }^{\mathrm{pp}}(\omega gT)g^3T^4.$$
(52)
For large $`\omega `$ the spectral function behaves as

$$\rho _{\pi \pi }^{\mathrm{pp}}(\omega )=\frac{d_A}{4\pi }\omega ^4\left[n(\omega /2)+\frac{1}{2}\right](\omega \omega _{\mathrm{pl}}),$$
(53)
which is, up to the prefactor, what we found in the scalar case as well.

The contribution when one transverse gluon is a quasiparticle and the other undergoes Landau damping (pole-cut contribution) is

$`\rho _{\pi \pi }^{\mathrm{pc}}(\omega )={\displaystyle \frac{2d_A}{3\pi ^2}}{\displaystyle _{\omega _{\mathrm{pl}}}^{\mathrm{}}}𝑑u`$ $`{\displaystyle \frac{f^2(u)}{\omega _T^{}[f(u)]}}Z_T[f(u)][n(u\omega )n(u)]\beta _T[u\omega ,f(u)]`$ (54)
$`\times \left[7f^4(u)10u(u\omega )f^2(u)+7u^2(u\omega )^2\right].`$
For soft $`\omega gT`$ the dominant contribution arises when the energy $`u`$ is hard, and we may use

$$\beta _T[u\omega ,f(u)]\frac{3\pi }{4}\frac{\omega _{\mathrm{pl}}^2}{\omega u^3},$$
(55)
as well as other simplifications given above. We find

$$\rho _{\pi \pi }^{\mathrm{pc}}(\omega gT)\frac{d_A}{\pi }\omega _{\mathrm{pl}}^2_{\omega _{\mathrm{pl}}}^{\mathrm{}}𝑑uu^2n^{}(u)g^2T^4.$$
(56)
For soft external frequencies the pole-cut contribution dominates over the pole-pole contribution (52).

For hard frequencies $`\omega T`$, the $`\beta _T`$ function determines the lower integration limit to be $`u_0=\omega /2+3\omega _{\mathrm{pl}}^2/(4\omega )`$. As a result, $`u`$ is always hard and we can simplify the integrand to arrive at

$$\rho _{\pi \pi }^{\mathrm{pc}}(\omega T)=\frac{d_A}{3\pi ^2}_{u_0}^{\mathrm{}}𝑑uu^3\left(4u^24u\omega +7\omega ^2\right)[n(u\omega )n(u)]\beta _T(u\omega ,u).$$
(57)
It is convenient to substitute $`u=\omega (z+1)/2`$ such that

$$\rho _{\pi \pi }^{\mathrm{pc}}(\omega T)=\frac{d_A\omega ^6}{48\pi ^2}_a^{\mathrm{}}𝑑z(z+1)^3\left(z^2+6\right)[n(\omega z_{})n(\omega z_+)]\beta _T(\omega z_{},\omega z_+),$$
(58)
where $`z_\pm =(z\pm 1)/2`$ and $`a=3\omega _{\mathrm{pl}}^2/(2\omega ^2)`$. The dominant contribution comes from the lower integration limit ($`za`$) where the $`\beta _T`$ function can be approximated as

$$\beta _T(\omega z_{},\omega z_+)\frac{4\pi }{\omega ^2}\frac{za}{(z+a)^2}.$$
(59)
Using this expression we find that

$$\rho _{\pi \pi }^{\mathrm{pc}}(\omega T)g^2T^4\mathrm{ln}(1/g).$$
(60)
For hard frequencies the pole-cut contribution is therefore suppressed compared to the pole-pole contribution (53). Finally, we found that the remaining cut-cut contribution when both gluons are transverse is suppressed with respect to the pole-pole contribution when $`\omega `$ is hard and to the pole-cut contribution when $`\omega `$ is soft.

In a similar way we have analysed the remaining longitudinal-transverse and longitudinal-longitudinal contributions in Eq. (44) with the result that they do not modify the conclusions drawn from the transverse-transverse contribution analysed above. In particular, for hard frequencies $`\omega \omega _{\mathrm{pl}}`$ transverse gluons dominate and the spectral function is given by Eq. (53). We do not need to be more explicit about the region where $`\omega gT`$ because the dominant contribution in this region arises from the pinching poles, screened by a finite width and/or external frequency, as we will show now.

For small external frequencies $`0\omega gT`$ the loop integral is dominated by the region of hard momentum and we only need to consider transverse gluons, since the residue for longitudinal gluons vanishes exponentially. In order to avoid pinch singularities we follow the same steps as in the scalar case and substitute for the one-particle spectral function a Breit-Wigner function, see Eq. (18). For transverse on-shell gluons with hard momentum the dispersion relation is $`\omega _𝐤=\sqrt{𝐤^2+m_{\mathrm{}}^2}`$ and the leading (momentum-independent) contribution to the damping rate is

$$\gamma =\frac{g^2}{4\pi }N_cT\mathrm{ln}(1/g),$$
(61)
where the logarithm is sensitive to the magnetic mass, $`\mathrm{ln}(\omega _{\mathrm{pl}}/m_{\mathrm{mag}})\mathrm{ln}(1/g)`$ with $`m_{\mathrm{mag}}g^2T`$.111111In QED there is no magnetic mass which could regularize the logarithmic divergence of the leading contribution to the damping rate. It turns out that the electron damping rate is ill defined . The mass $`m_{\mathrm{}}`$ only plays a subdominant role and is neglected below. Evaluating the integral over $`k^0`$ exactly as in the scalar case, we find

$`\rho _{\pi \pi }(\omega )={\displaystyle \frac{d_A}{3}}{\displaystyle _𝐤}{\displaystyle \frac{1}{|𝐤|}}\{`$ $`[n(|𝐤|)n(|𝐤|\omega )]I(\omega ,𝐤)A(\omega ,𝐤){\displaystyle \frac{|𝐤|\omega }{\omega ^2+4\gamma ^2}}`$ (62)
$`[\omega \omega ]\},`$
where $`I(\omega ,𝐤)`$ was defined in Eq. (25) and $`A(\omega ,𝐤)=𝐤^2\left(4𝐤^24\omega |𝐤|+7\omega ^2\right)`$. For $`\omega \omega _{\mathrm{pl}}`$ this expression simplifies to

$$\rho _{\pi \pi }(\omega )=\frac{16d_A}{3}_𝐤𝐤^2n^{}(|𝐤|)\frac{\omega \gamma }{\omega ^2+4\gamma ^2}(0\omega \omega _{\mathrm{pl}}),$$
(63)
which is $`2d_A`$ times the scalar result. In the region $`\gamma \omega \omega _{\mathrm{pl}}`$ we find

$$\rho _{\pi \pi }(\omega )/T^4=\frac{16d_A}{3T^4}\frac{\gamma }{\omega }_𝐤𝐤^2n^{}(|𝐤|)=\frac{32\pi ^2}{45}d_A\frac{\gamma }{\omega }.$$
(64)
As in the scalar case, the three-loop ladder diagram with one rung contributes in this region at the same order. Note that the contribution from the pinching poles at $`\omega gT`$,

$$\rho _{\pi \pi }(\omega gT)gT^4\mathrm{ln}(1/g)$$
(65)
actually dominates over the contribution (56) from the collective (HTL) excitations in this region. The shear viscosity in the one-loop approximation follows from (63) as

$$\eta _{1\mathrm{loop}}=\frac{8\pi ^2d_A}{225}\frac{T^4}{\gamma }.$$
(66)
However, for the complete spectral function at small frequencies $`0\omega \gamma `$ the effects of ladder diagrams must be taken into account. As was mentioned in the Introduction, kinetic theory predicts that the shear viscosity is parametrically larger than the one-loop result (66) . Therefore, for vanishing frequency the spectral function due to ladder diagrams is expected to be larger than the one-loop contribution and to behave as $`\rho _{\pi \pi }^{\mathrm{ladder}}(\omega 0)=20\eta \omega `$, with

$$\eta \frac{N_c^21}{N_c^2}\frac{T^3}{g^4\mathrm{ln}(1/g)}.$$
(67)
Hence the slope of the spectral function is much steeper close to the origin, compared to the one-loop result. Up to frequencies $`\omega \gamma `$ the actual behaviour of the spectral function depends on the pinching-poles contribution of ladders diagrams. When $`\gamma \omega \omega _{\mathrm{pl}}`$ pinching-poles contributions of the ladders are subdominant by a factor $`g^2T/\omega `$ and the spectral function decreases as $`g^2T/\omega `$ until it meets with the rising contribution from Eq. (53). For large $`\omega `$ the spectral function increases as $`\omega ^4`$.

The euclidean correlator at high temperature can be computed as in the scalar case and we find, at leading order in $`g`$,

$$G_{\pi \pi }^E(\tau )/T^5=\frac{4\pi ^2d_A}{\mathrm{sin}^5u}\left\{(\pi u)\left[11\mathrm{cos}u+\mathrm{cos}3u\right]+6\mathrm{sin}u+2\mathrm{sin}3u\right\}+\frac{8\pi ^2d_A}{45},$$
(68)
with $`u=2\pi \tau T`$. The last constant term reflects the low-frequency region.

## 5 Conclusions

We have studied the spectral function relevant for the shear viscosity in scalar and nonabelian gauge theories at high temperature as a function of the external frequency. While for small frequencies ladder diagrams are important in the scalar case and essential in the nonabelian case, for higher frequencies a simple one-loop computation yields the dominant contribution.

We found that the spectral function has a characteristic shape: for small frequencies the spectral function rises, reaches a local maximum and decreases as $`1/\omega `$. This contribution is due to scattering processes in the plasma and is enhanced due to nearly pinching poles. We referred to this contribution as the low-frequency contribution. For higher frequencies, decay/creation processes dominate and the spectral function increases essentially as $`\omega ^4`$. This contribution is referred to as the high-frequency contribution.

In order to extract transport coefficients from euclidean lattice correlators, a simple ansatz, essentially a Breit-Wigner spectral function, was introduced in Ref.  to model spectral functions of components of the energy-momentum tensor. The resulting three-parameter ansatz was subsequently used in Refs.  to determine transport coefficients in hot gauge theories from lattice simulations. Unfortunately, as we have seen in this paper, at high temperature spectral functions of composite operators do not resemble simple Breit-Wigner functions at all. In order to improve this analysis, we propose therefore to use a different ansatz, which is written as the sum of two terms:

$$\rho _{\pi \pi }(\omega )=\rho _{\pi \pi }^{\mathrm{low}}(\omega )+\rho _{\pi \pi }^{\mathrm{high}}(\omega ).$$
(69)
The high-frequency part can be described by Eqs. (53) or (17) with $`m`$ as a possible free parameter. For the low-frequency part we note that the spectral function is odd, increases linearly with $`\omega `$ for small $`\omega `$ and decreases with $`1/\omega `$ for larger $`\omega `$. A simple ansatz reflecting this is

$$\rho _{\pi \pi }^{\mathrm{low}}(x)/T^4=x\frac{b_1+b_2x^2+b_3x^4+\mathrm{}}{1+c_1x^2+c_2x^4+c_3x^6+\mathrm{}},x=\omega /T,$$
(70)
with $`b_i=c_i=0`$, $`i>n`$ for given $`n`$. The viscosity is given by $`\eta /T^3=b_1/20`$. The ansatz (69), with free parameters $`m`$, $`b_i`$, and $`c_i`$, should be used in Eq. (6) to fit the corresponding euclidean correlator to the numerical results. When insisting on a three-parameter fit, one may take $`n=1`$ which leaves $`m`$, $`b_1`$ and $`c_1`$ to be determined.

Concerning the euclidean correlator, we found that the dominant $`\tau `$ dependence is determined by the high-frequency part ($`\omega T`$) of the spectral function. However, around $`\tau T=1/2`$ both the high- and the low-frequency regions of the spectral function contribute at the same order. The low-frequency contribution is of special importance since transport coefficients are determined by the slope of the spectral function at zero frequency and a precise calculation of the spectral function at low frequencies is therefore essential. We found that no matter how complicated the spectral function up to frequencies $`\omega gT`$ ($`g1`$) might be, its contribution to the euclidean correlator will be a $`\tau `$ independent constant. This latter feature poses a severe challenge for the MEM approach, since the spectral function at low frequencies should be reconstructed from the knowledge of a single constant alone. It turns out that this constant carries information on the transport coefficient (or more generally on ladder diagrams) only in subleading contributions. As a result euclidean correlators are remarkably insensitive to transport coefficients, which makes it extremely difficult to extract those at weak coupling.

We emphasize that these results are not specific for the correlator we considered here. For instance, in the first paper of Ref.  the current-current correlator relevant for thermal dilepton rates in QCD was studied on the lattice and an enhancement in the central value of the euclidean correlator compared to the free result was observed. This enhancement might be accounted for by the pinching-poles contribution in the low-frequency region of the spectral function.121212Note, however, that in the lattice study the spectral function reconstructed with the Maximal Entropy Method seems to vanish below $`\omega 3T`$.

Acknowledgements
It is a pleasure to thank Ulrich Heinz and Eric Braaten for discussions and comments. G. A. thanks Jan Smit and Nucu Stamatescu as well. We thank Guy Moore for a careful reading of and valuable comments on the manuscript. G. A. is supported by the Ohio State University through a Postdoctoral Fellowship and by the U. S. Department of Energy under Contract No. DE-FG02-01ER41190. J. M. M. R. is supported by a Postdoctoral Fellowship from the Basque Government.


text:
# Construction of particular solutions to nonlinear equations of Mathematical Physics by using matrix algebraic equation

## 1 Introduction

Analysis of nonlinear Partial Differential Equations (PDE) is severe problem in mathematical physics. Many different methods have been developed for analytical investigation of nonlinear PDE during last decades: Inverse Scattering Problem , Sato theory , Hirota bilinear method , Penlevé method , $`\overline{}`$-problem , with some generalizations , nonlinear $`\overline{}`$-problem , recent modification of the dressing method based on the algebraic matrix equation . A wide class of PDE (so-called completely integrable systems) has been studied better then others. Nevertheless, there are many methods which work in nonintegrable case as well: .

We represent the method for construction of families of particular solutions to multidimensional nonlinear PDE without relation to the complete integrability. This method is based on general properties of linear algebraic matrix equations. Essentially we develop some ideas represented in the ref. and recently in the ref. .

General algorithm is discussed in the Sec. 2. Sec. 3 deals with example of (2+2)-dimensional system of PDE together with family of particular solutions. By using an appropriate reduction, this system leads to Burgers equation. This example is simple demonstration of the method. Conclusions are given in the Sec.4.

## 2 General results

The algorithm represented in this section is based on the fundamental properties of linear matrix algebraic equation, which is written in the following form:

$`\mathrm{\Psi }U=\mathrm{\Phi },`$ (1)
where $`\mathrm{\Psi }=\{\psi _{ij}\}`$ is $`N\times N`$ matrix $`U=[u_1,\mathrm{},u_N]^T`$, $`\mathrm{\Phi }=[\varphi _1,\mathrm{},\varphi _N]^T`$. Let us recall these properties.

1. If $`\mathrm{\Psi }`$ is nondegenerate matrix, i.e.

$`det\mathrm{\Psi }0,`$ (2)
then equation (1) has unique solution which can be written in the following form:

$`U=\mathrm{\Psi }^1\mathrm{\Phi }.`$ (3)
Only nondegenerate matrices $`\mathrm{\Psi }`$ will be considered hereafter.
2. (consequence of the previous property) If $`\mathrm{\Phi }=0`$ and condition (2) is held, then the equation (1) has only the trivial solution

$`U0.`$ (4)
3. (superposition principle) Consider the set of $`K`$ matrix equations with the same matrix $`\mathrm{\Psi }`$:

$`\mathrm{\Psi }U_i=\mathrm{\Phi }_i,i=1,\mathrm{},K.`$ (5)
Then for any set of scalars $`b_k`$ ($`k=1,\mathrm{},K`$), function $`\stackrel{~}{U}=_{k=1}^Kb_kU_k`$ is solution of the following matrix equation

$`\mathrm{\Psi }\stackrel{~}{U}={\displaystyle \underset{k=1}{\overset{K}{}}}b_k\mathrm{\Phi }_k.`$ (6)
4. (consequence of properties 2 and 3) If columns $`\mathrm{\Phi }_i`$ are linearly dependent, i.e there are scalars $`a_k`$, $`k=1,\mathrm{},K`$, such that

$`{\displaystyle \underset{k=1}{\overset{K}{}}}a_k\mathrm{\Phi }_k=0,`$ (7)
then

$`{\displaystyle \underset{i=1}{\overset{K}{}}}a_iU_i=0.`$ (8)

Note that analogous properties of linear integral equation have been used in the classical dressing method based on the $`\overline{}`$-problem .

To relate the matrix equation (1) with the system of nonlinear PDE, one needs to introduce the set of additional variables $`𝐱=(x_1,\mathrm{},x_M)`$ in matrices $`\mathrm{\Psi }`$, $`\mathrm{\Phi }`$ and $`U`$, which are independent variables of succeeding system of nonlinear PDE. $`M`$ is dimension of $`𝐱`$-space. For this purpose we use the following system of linear differential equations

$`M_k\mathrm{\Phi }+N_k\mathrm{\Psi }=\mathrm{\Psi }A_k,k=1,\mathrm{},P,`$ (9)
where $`M_k`$ and $`N_k`$ are linear differential operators with matrix coefficients; $`A_k`$ are $`N\times N`$ matrices, which are constant for the sake of simplicity. The number of these equations, $`P`$, depends on situation. Equations (9) introduce variables $`𝐱`$ into the matrices $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$. Due to the eq.(3), elements $`u_i`$ of the column $`U`$ are functions of variables $`𝐱`$ as well.

Operators $`M_k`$ and $`N_k`$ may not be arbitrary. They have to satisfy two conditions:

1. Overdetermined system (9) is compatible.
2. Operators $`M_k`$ and $`N_k`$ have to provide existence of differential operator $``$ with non-constant matrix coefficients such that

$`(\mathrm{\Phi }\mathrm{\Psi }W)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{P}{}}}(M_k\mathrm{\Phi }+N_k\mathrm{\Psi })\stackrel{~}{U}_k+\mathrm{\Psi }\stackrel{~}{V}=\mathrm{\Psi }B(U)0,`$ (10)
$`B(U)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{P}{}}}A_k\stackrel{~}{U}_k+\stackrel{~}{V},`$
where $`\stackrel{~}{U}_k`$ and $`\stackrel{~}{V}`$ are some matrix functions of elements $`u_i`$ and their derivatives (see eqs.(23-28) as an example).

Since $`\text{det}\mathrm{\Psi }0`$, the last equation means:

$`B(U)=0,`$ (11)
which represents the system of nonlinear PDE on functions $`u_i`$. The structure of this system is defined by the equations (9).

One can see that the matrix equation (1) is generalization of the linear systems used in the refs.. In fact, let us introduce notations

$`\varphi _i\varphi ^{(i)},\psi _{ij}\psi _j^{(i)},`$ (12)
and rewrite the matrix equation (1) in the form of system of $`N`$ equations:

$`\varphi ^{(i)}={\displaystyle \underset{j=1}{\overset{N}{}}}u_j\psi _j^{(i)},i=1,\mathrm{},N.`$ (13)
Two evident reductions are possible (these reductions will not be used in the succeeding sections).

1. If

$`\psi _j^{(i)}=_{x_1}^j\psi ^{(i)},\varphi ^{(i)}=_{x_1}^{(N+1)}\psi ^{(i)},i,j=1,\mathrm{},N,`$ (14)
then system (13) is system of $`N`$ ordinary differential equations on the function $`\psi ^{(i)}`$ which is basic for the algorithm represented in the ref..

2. Introduce multi-index $`\beta =(\beta _1,\mathrm{},\beta _M)`$, $`\beta =_{k=1}^M\beta _k`$ instead of index $`j`$ ($`M`$ is dimension of $`𝐱`$-space), and use notations

$`\psi _\beta ^{(i)}`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{M}{}}}_{x_k}^{\beta _k}\psi ^{(i)}^\beta \psi ^{(i)},\beta =1,\mathrm{},N,`$ (15)
$`\varphi ^{(i)}`$ $`=`$ $`^\alpha \psi ^{(i)},\alpha =N+1.`$ (16)
Then the system (13) can be written in the form

$`^\alpha \psi ^{(i)}={\displaystyle \underset{\beta =0}{\overset{N}{}}}u_\beta ^\beta \psi ^{(i)},i=1,\mathrm{},Q,`$ (17)
where $`Q`$ is the number of terms in the right hand side of the eq.(17), $`\alpha `$ takes any value such that $`\alpha =N+1`$. Eq.(17) is basic in the ref..

Hereafter, we will use eq. (13) instead of matrix equation (1). Superscript $`i`$ runs the values from 1 to $`N`$, unless otherwise specified.

## 3 Example: (2+2)-dimensional generalization of Burgers equation

In this section we consider (2+2)-dimensional equations with following notations for independent variables:

$`xx_1,yx_2,tx_3,\tau x_4.`$ (18)
For convenience of construction, we separate the first term of the sum in the eq. (13):

$`\varphi ^{(i)}=u_1\psi _1^{(i)}+{\displaystyle \underset{k=2}{\overset{N}{}}}u_k\psi _k^{(i)}.`$ (19)
Introduce two differential operators $`L_1`$ and $`L_2`$:

$`L_1`$ $`=`$ $`l_1+\delta _{11}u_{1}^{}{}_{x}{}^{}+\delta _{12}u_{1}^{}{}_{y}{}^{},l_1=_t+_x^2+\gamma _1_x_y,`$ (20)
$`L_2`$ $`=`$ $`l_2+\delta _{21}u_{1}^{}{}_{x}{}^{}+\delta _{22}u_{1}^{}{}_{y}{}^{},l_2=_\tau +_y^2+\gamma _2_x_y.`$
Let us apply operator $`L_1`$ to the eq.(19) and investigate the result

$`l_1\varphi ^{(i)}\delta _{11}u_{1}^{}{}_{x}{}^{}\varphi ^{(i)}\delta _{12}u_{1}^{}{}_{y}{}^{}\varphi ^{(i)}+(L_1u_1)\psi _1^{(i)}+u_1(l_1\psi _1^{(i)})+`$ (21)
$`u_{1}^{}{}_{x}{}^{}(2\psi _{1}^{(i)}{}_{x}{}^{}+\gamma _1\psi _{1}^{(i)}{}_{y}{}^{})+\gamma _1u_{1}^{}{}_{y}{}^{}\psi _{1}^{(i)}{}_{x}{}^{}+`$
$`{\displaystyle \underset{k=1}{\overset{N}{}}}\left((L_1u_k)\psi _k^{(i)}+u_k(l_1\psi _k^{(i)})+u_{k}^{}{}_{x}{}^{}(2\psi _{k}^{(i)}{}_{x}{}^{}+\gamma _1\psi _{k}^{(i)}{}_{y}{}^{})+\gamma _1u_{k}^{}{}_{y}{}^{}\psi _{k}^{(i)}{}_{x}{}^{}\right)`$ $`=`$ $`0`$
Our purpose is to represent the equation (21) in the form of linear homogeneous equation

$`{\displaystyle \underset{j=1}{\overset{N}{}}}F_{1k}\psi _k^{(i)}=0,`$ (22)
where $`F_{1k}`$ are some expressions of $`u_k`$ and their derivatives. Owing to the property 4 (eqs.(7) and (8)), this would mean that $`F_{1k}0`$, $`k=1,\mathrm{},N`$. The simplest way to result in the eq. (22) is introduction of the system of differential equations (9) on the functions $`\psi _j^{(i)}`$ and $`\varphi ^{(i)}`$, which looks as follows:

$`\varphi ^{(i)}{\displaystyle \frac{1}{\delta _{11}}}(2\psi _{1}^{(i)}{}_{x}{}^{}+\gamma _1\psi _{1}^{(i)}{}_{y}{}^{})`$ $`=`$ $`{\displaystyle \underset{k=2}{\overset{N}{}}}\stackrel{~}{\alpha }_{1k}\psi _k^{(i)},`$ (23)
$`\varphi ^{(i)}{\displaystyle \frac{\gamma _1}{\delta _{12}}}\psi _{1}^{(i)}{}_{x}{}^{}`$ $`=`$ $`{\displaystyle \underset{k=2}{\overset{N}{}}}\alpha _{1k}\psi _k^{(i)},`$ (24)
$`\psi _{2}^{(i)}{}_{x}{}^{}`$ $`=`$ $`\beta _1\psi _1^{(i)}+{\displaystyle \underset{k=2}{\overset{N}{}}}\alpha _{3k}\psi _k^{(i)},\psi _{2}^{(i)}{}_{y}{}^{}=\beta _2\psi _1^{(i)}+{\displaystyle \underset{k=2}{\overset{N}{}}}\alpha _{4k}\psi _k^{(i)},`$ (25)
$`\psi _{j}^{(i)}{}_{x}{}^{}`$ $`=`$ $`{\displaystyle \underset{k=2}{\overset{N}{}}}\alpha _{(2j1)k}\psi _k^{(i)},\psi _{j}^{(i)}{}_{y}{}^{}={\displaystyle \underset{k=2}{\overset{N}{}}}\alpha _{(2j)k}\psi _k^{(i)},j>2,`$ (26)
$`l_1\psi _j^{(i)}`$ $`=`$ $`{\displaystyle \underset{k=2}{\overset{N}{}}}\beta _{1jk}\psi _k^{(i)},l_2\psi _j^{(i)}={\displaystyle \underset{k=2}{\overset{N}{}}}\beta _{2jk}\psi _k^{(i)},j1,`$ (27)
$`l_1\varphi ^{(i)}`$ $`=`$ $`{\displaystyle \underset{k=2}{\overset{N}{}}}\beta _{10k}\psi _k^{(i)},l_2\varphi ^{(i)}={\displaystyle \underset{k=2}{\overset{N}{}}}\beta _{20k}\psi _k^{(i)}.`$ (28)
Here parameters $`\stackrel{~}{\alpha }_{ij},\alpha _{ij},\beta _i,\beta _{ijk}`$ have to provide compatibility of the overdetermined system (23)-(28). Equation (23) can be simplified if one eliminates $`\varphi ^{(i)}`$ by using the eq.(24):

$`\psi _{1}^{(i)}{}_{y}{}^{}G\psi _{1}^{(i)}{}_{x}{}^{}={\displaystyle \underset{k=2}{\overset{N}{}}}\alpha _{2k}\psi _k^{(i)},G={\displaystyle \frac{\delta _{11}}{\delta _{12}}}{\displaystyle \frac{2}{\gamma _1}},\alpha _{2k}={\displaystyle \frac{\delta _{11}}{\gamma _1}}(\alpha _{1k}\stackrel{~}{\alpha }_{1k}).`$ (29)
The system (24)-(29) will be used below.

Now the equation (21) takes the form (22) with

$`F_{11}=u_{1}^{}{}_{t}{}^{}+u_{1}^{}{}_{xx}{}^{}+\gamma _1u_{1}^{}{}_{xy}{}^{}+\delta _{11}u_1u_{1}^{}{}_{x}{}^{}+\delta _{12}u_1u_{1}^{}{}_{y}{}^{}+`$ (30)
$`(2\beta _1+\beta _2\gamma _1)u_{2}^{}{}_{x}{}^{}+\beta _1\gamma _1u_{2}^{}{}_{y}{}^{}`$ $`=`$ $`0.`$
Expressions for $`F_{1k}`$ with $`k>1`$ will not be used.

Analogously, apply operator $`L_2`$ to the eq.(19):

$`l_2\varphi ^{(i)}\delta _{21}u_{1}^{}{}_{x}{}^{}\varphi ^{(i)}\delta _{22}u_{1}^{}{}_{y}{}^{}\varphi ^{(i)}+(L_2u_1)\psi _1^{(i)}+u_1(l_2\psi _1^{(i)})+`$ (31)
$`u_{1}^{}{}_{y}{}^{}(2\psi _{1}^{(i)}{}_{y}{}^{}+\gamma _2\psi _{1}^{(i)}{}_{x}{}^{})+\gamma _2u_{1}^{}{}_{x}{}^{}\psi _{1}^{(i)}{}_{y}{}^{}+`$
$`{\displaystyle \underset{k=1}{\overset{N}{}}}\left((L_2u_k)\psi _k^{(i)}+u_k(l_2\psi _k^{(i)})+u_{k}^{}{}_{y}{}^{}(2\psi _{k}^{(i)}{}_{y}{}^{}+\gamma _2\psi _{k}^{(i)}{}_{x}{}^{})+\gamma _2u_{k}^{}{}_{x}{}^{}\psi _{k}^{(i)}{}_{y}{}^{}\right)`$ $`=`$ $`0`$
Due to the relations (24)-(29) one gets analogous expression

$`{\displaystyle \underset{j=1}{\overset{N}{}}}F_{2k}\psi _k^{(i)}=0,`$ (32)
if only

$`\delta _{21}={\displaystyle \frac{\gamma _2}{\gamma _1^2}}(\delta _{11}\gamma _12\delta _{12}),\delta _{22}={\displaystyle \frac{1}{\gamma _1^2}}(\delta _{12}\gamma _1\gamma _2+2\delta _{11}\gamma _14\delta _{12}).`$ (33)
Coefficients $`F_{2k}`$ are expressed in terms of the functions $`u_k`$ and their derivatives with

$`F_{21}=u_{1}^{}{}_{\tau }{}^{}+u_{1}^{}{}_{yy}{}^{}+\gamma _2u_{1}^{}{}_{xy}{}^{}+{\displaystyle \frac{\gamma _2}{\gamma _1^2}}(\delta _{11}\gamma _12\delta _{12})u_1u_{1}^{}{}_{x}{}^{}+`$ (34)
$`{\displaystyle \frac{1}{\gamma _1^2}}(\delta _{12}\gamma _1\gamma _2+2\delta _{11}\gamma _14\delta _{12})u_1u_{1}^{}{}_{y}{}^{}+\beta _2\gamma _2u_{2}^{}{}_{x}{}^{}+(2\beta _2+\beta _1\gamma _2)u_{2}^{}{}_{y}{}^{}`$ $`=`$ $`0,`$
Eqs.(30) and (34) form the system of two nonlinear PDE on the functions $`u_1`$ and $`u_2`$.

One shell discuss the structure of the system (24-29). Operators $`L_1`$ and $`L_2`$, given by the formulas (20), determine the left hand sides of the equations in this system. The structure of the system of nonlinear PDE is defined by the operators $`L_1`$ and $`L_2`$ and by the right hand sides of the equations in the system (24-29). Only terms proportional to $`\psi _1^{(i)}`$ effect on the structure of the nonlinear system of PDE. At least one of the equations (24-29) must have this term (see eqs.(25)). Otherwise, only simplest solutions to the nonlinear system of PDE are available. Other terms enrich the family of available particular solutions.

Function $`\varphi ^{(i)}`$ may not be linear combination of functions $`\psi _j^{(i)}`$ (see eqs. (23,24)). Otherwise all coefficients $`u_j`$ in eq.(19) are zero. Equations (27) and (28) introduce dependence on variables $`t`$ and $`\tau `$. Right hand sides of these equations are zero in the examples below for the sake of simplicity.

### 3.1 Construction of particular solutions

This section is devoted to construction of some families of particular solutions to the system (30,34). We will use simplified form of the system (24-29):

$`\varphi ^{(i)}`$ $`=`$ $`{\displaystyle \frac{\gamma _1}{\delta _{12}}}\psi _{1}^{(i)}{}_{x}{}^{}+\alpha _{12}\psi _2^{(i)},`$ (35)
$`\psi _{1}^{(i)}{}_{y}{}^{}`$ $`=`$ $`G\psi _{1}^{(i)}{}_{x}{}^{},G={\displaystyle \frac{\delta _{11}}{\delta _{12}}}{\displaystyle \frac{2}{\gamma _1}},`$ (36)
$`\psi _{2}^{(i)}{}_{x}{}^{}`$ $`=`$ $`\beta _1\psi _1^{(i)}+\alpha _{32}\psi _2^{(i)},\psi _{2}^{(i)}{}_{y}{}^{}=\beta _2\psi _1^{(i)}+\alpha _{42}\psi _2^{(i)},`$ (37)
$`\psi _{j}^{(i)}{}_{x}{}^{}`$ $`=`$ $`{\displaystyle \underset{k=2}{\overset{N}{}}}\alpha _{(2j1)k}\psi _k^{(i)},\psi _{j}^{(i)}{}_{y}{}^{}={\displaystyle \underset{k=2}{\overset{N}{}}}\alpha _{(2j)k}\psi _k^{(i)},j>2,`$ (38)
$`l_k\psi _j^{(i)}`$ $`=`$ $`0,l_k\varphi ^{(i)}=\mathrm{\hspace{0.33em}0},k=1,2,j=1,\mathrm{},N`$ (39)
Below are two examples with $`N=2`$ and $`N=4`$. Example with $`N=3`$ does not have a principal difference with the case $`N=2`$. Solutions from both families considered below are parameterized by the arbitrary functions of single variable.

#### 3.1.1 One-dimensional kink

Let $`N=2`$, superscript $`i`$ takes the values 1 and 2. The system (19) has the following form:

$`\varphi ^{(i)}`$ $`=`$ $`u_1\psi _1^{(i)}+u_2\psi _2^{(i)},`$ (40)
with

$`u_1`$ $`=`$ $`{\displaystyle \frac{\mathrm{}_1}{\mathrm{}}},u_2={\displaystyle \frac{\mathrm{}_2}{\mathrm{}}},`$ (41)
$`\mathrm{}=\left|\begin{array}{cc}\psi _1^{(1)}& \psi _2^{(1)}\\ \psi _1^{(2)}& \psi _2^{(2)}\end{array}\right|,\mathrm{}_1=\left|\begin{array}{cc}\varphi ^{(1)}& \psi _2^{(1)}\\ \varphi ^{(2)}& \psi _2^{(2)}\end{array}\right|,\mathrm{}_2=\left|\begin{array}{cc}\psi _1^{(1)}& \varphi ^{(1)}\\ \psi _1^{(2)}& \varphi ^{(2)}\end{array}\right|,`$ (48)
We need only equations (35)-(37) and (39) (with $`j=1,2`$) to find functions $`\psi _j^{(i)}`$ ($`j=1,2`$) and $`\varphi ^{(i)}`$. They admit the following solutions

$`\psi _1^{(i)}`$ $`=`$ $`{\displaystyle c_1^{(i)}(k)e^{kx+qy+\omega t+\nu \tau }𝑑k},`$ (49)
$`\psi _2^{(i)}`$ $`=`$ $`{\displaystyle c_2^{(i)}(k)e^{kx+qy+\omega t+\nu \tau }𝑑k},`$ (50)
$`\varphi ^{(i)}`$ $`=`$ $`{\displaystyle c^{(i)}(k)e^{kx+qy+\omega t+\nu \tau }𝑑k},`$ (51)
$`c^{(i)}(k)`$ $`=`$ $`{\displaystyle \frac{c_2^{(i)}(k)}{\beta _1\delta _{12}}}(\gamma _1k^2\alpha _{32}\gamma _1k+\alpha _{12}\beta _1\delta _{12}),`$ (52)
$`c_1^{(i)}(k)`$ $`=`$ $`{\displaystyle \frac{c_2^{(i)}(k)}{\beta _1}}(k\alpha _{32}),`$ (53)
$`q`$ $`=`$ $`{\displaystyle \frac{\beta _2}{\beta _1}}k,`$ (54)
$`\omega `$ $`=`$ $`k^2kq\gamma _1={\displaystyle \frac{k^2}{\beta _1}}(\beta _1+\beta _2\gamma _1),`$ (55)
$`\nu `$ $`=`$ $`q^2kq\gamma _2={\displaystyle \frac{k^2\beta _2}{\beta _1^2}}(\beta _2+\beta _1\gamma _2),`$ (56)
$`\delta _{11}`$ $`=`$ $`{\displaystyle \frac{\delta _{12}}{\beta _1\gamma _1}}(\beta _2\gamma _1+2\beta _1),\alpha _{42}={\displaystyle \frac{\beta _2}{\beta _1}}\alpha _{32}`$ (57)
One can see from the eqs.(49-53) that functions $`\psi _2^{(i)}`$ are arbitrary functions of single variable (say, variable $`x`$). So that appropriate solutions (41) depend on two arbitrary functions of single variable.

As a simple example, let us take the following expressions for the functions $`\psi _j^{(i)}`$ and $`\varphi ^{(i)}`$:

$`\psi _1^{(1)}`$ $`=`$ $`{\displaystyle \frac{k_1\alpha _{32}}{\beta _1}}s_1E_1,`$ (58)
$`\psi _1^{(2)}`$ $`=`$ $`{\displaystyle \frac{k_2\alpha _{32}}{\beta _1}}s_2E_2+{\displaystyle \frac{k_3\alpha _{32}}{\beta _1}}s_3E_3,`$ (59)
$`\psi _2^{(1)}`$ $`=`$ $`s_1E_1,`$ (60)
$`\psi _2^{(2)}`$ $`=`$ $`s_2E_2+s_3E_3,`$ (61)
$`\varphi ^{(1)}`$ $`=`$ $`{\displaystyle \frac{s_1E_1}{\beta _1\delta _{12}}}(\alpha _{12}\beta _1\delta _{12}\alpha _{32}\gamma _1k_1+\gamma _1k_1^2),`$ (62)
$`\varphi ^{(2)}`$ $`=`$ $`{\displaystyle \frac{s_2E_2}{\beta _1\delta _{12}}}(\alpha _{12}\beta _1\delta _{12}\alpha _{32}\gamma _1k_2+\gamma _1k_2^2)+`$
$`{\displaystyle \frac{s_3E_3}{\beta _1\delta _{12}}}(\alpha _{12}\beta _1\delta _{12}\alpha _{32}\gamma _1k_3+\gamma _1k_3^2),`$
$`E_n=e^{k_nx+q_ny+\omega _nt+\nu _n\tau },n=1,2,3,`$
where

$`q_n={\displaystyle \frac{\beta _2}{\beta _1}}k_n,\omega _n={\displaystyle \frac{k_n^2}{\beta _1}}(\beta _1+\beta _2\gamma _1),\nu _n={\displaystyle \frac{k_n^2\beta _2}{\beta _1^2}}(\beta _2+\beta _1\gamma _2),`$ (64)
$`n=1,2,3,`$
which are in agreement with eqs.(49)-(57). Determinants (48) are expressed as follows:

$`\mathrm{}`$ $`=`$ $`{\displaystyle \frac{s_1E_1}{\beta _1}}((k_1k_2)s_2E_2+(k_1k_3)s_3E_3),`$ (65)
$`\mathrm{}_1`$ $`=`$ $`{\displaystyle \frac{\gamma _1s_1E_1}{\beta _1\delta _{12}}}((k_1k_2)(\alpha _{32}k_1k_2)s_2E_2+`$
$`(k_1k_3)(\alpha _{32}k_1k_3)s_3E_3),`$
$`\mathrm{}_2`$ $`=`$ $`{\displaystyle \frac{s_1E_1}{\beta _1^2\delta _{12}}}\times `$
$`((k_1k_2)(\alpha _{12}\beta _1\delta _{12}\alpha _{32}^2\gamma _1+\alpha _{32}\gamma _1(k_1+k_2)\gamma _1k_1k_2)s_2E_2+`$
$`(k_1k_3)(\alpha _{12}\beta _1\delta _{12}\alpha _{32}^2\gamma _1+\alpha _{32}\gamma _1(k_1+k_3)\gamma _1k_1k_3)s_3E_3).`$
Functions $`u_1`$ and $`u_2`$ have no singularities provided that $`\mathrm{}0`$ for all values of parameters $`x,y,t,\tau `$, which happens if only

$`\text{sign}((k_1k_2)s_2)=\text{sign}((k_1k_3)s_3).`$ (68)
Taking into account relations (57), system (30,34) gets the following form:

$`u_{1}^{}{}_{t}{}^{}+u_{1}^{}{}_{xx}{}^{}+\gamma _1u_{1}^{}{}_{xy}{}^{}+{\displaystyle \frac{\delta _{12}}{\beta _1\gamma _1}}(2\beta _1+\beta _2\gamma _1)u_1u_{1}^{}{}_{x}{}^{}+`$ (69)
$`\delta _{12}u_1u_{1}^{}{}_{y}{}^{}+(2\beta _1+\beta _2\gamma _1)u_{2}^{}{}_{x}{}^{}+\beta _1\gamma _1u_{2}^{}{}_{y}{}^{}`$ $`=`$ $`0,`$
$`u_{1}^{}{}_{\tau }{}^{}+u_{1}^{}{}_{yy}{}^{}+\gamma _2u_{1}^{}{}_{xy}{}^{}+{\displaystyle \frac{\gamma _2\delta _{12}\beta _2}{\gamma _1\beta _1}}u_1u_{1}^{}{}_{x}{}^{}+{\displaystyle \frac{\delta _{12}}{\gamma _1\beta _1}}(2\beta _2+\beta _1\gamma _2)u_1u_{1}^{}{}_{y}{}^{}+`$ (70)
$`\beta _2\gamma _2u_{2}^{}{}_{x}{}^{}+(2\beta _2+\beta _1\gamma _2)u_{2}^{}{}_{y}{}^{}`$ $`=`$ $`0,`$
Note that owing to the relations (54-56) functions $`u_1`$ and $`u_2`$ depend on two variables rather then four:

$`X=x+{\displaystyle \frac{\beta _2}{\beta _1}}y,T={\displaystyle \frac{\beta _1+\beta _2\gamma _1}{\beta _1}}t+{\displaystyle \frac{\beta _2(\beta _2+\beta _1\gamma _2)}{\beta _1^2}}\tau `$ (71)
Because of this fact, functions $`u_1`$ and $`u_2`$, given by the formulas (41,65-3.1.1), represent single kink, which is essentially one-dimensional. Moreover, in virtue of eqs.(71), the above system (69,70) can be written in the form of differentiated Burghers equation on the function $`u_1`$ with independent variables $`X`$ and $`T`$.

#### 3.1.2 Multidimensional kinks

In this section $`N=4`$, superscript $`i`$ takes values from 1 to 4, unless otherwise specified. The system (19) has the form

$`\varphi ^{(i)}=u_1\psi _1^{(i)}+u_2\psi _2^{(i)}+u_3\psi _3^{(i)}+u_4\psi _4^{(i)}.`$ (72)
Formulas (41) for the functions $`u_1`$ and $`u_2`$ are held with determinants in the following form:

$`\mathrm{}`$ $`=`$ $`\left|\begin{array}{cccc}\psi _1^{(1)}& \psi _2^{(1)}& \psi _3^{(1)}& \psi _4^{(1)}\\ \psi _1^{(2)}& \psi _2^{(2)}& \psi _3^{(2)}& \psi _4^{(2)}\\ \psi _1^{(3)}& \psi _2^{(3)}& \psi _3^{(3)}& \psi _4^{(3)}\\ \psi _1^{(4)}& \psi _2^{(4)}& \psi _3^{(4)}& \psi _4^{(4)}\end{array}\right|,\mathrm{}_1=\left|\begin{array}{cccc}\varphi ^{(1)}& \psi _2^{(1)}& \psi _3^{(1)}& \psi _4^{(1)}\\ \varphi ^{(2)}& \psi _2^{(2)}& \psi _3^{(2)}& \psi _4^{(2)}\\ \varphi ^{(3)}& \psi _2^{(3)}& \psi _3^{(3)}& \psi _4^{(3)}\\ \varphi ^{(4)}& \psi _2^{(4)}& \psi _3^{(4)}& \psi _4^{(4)}\end{array}\right|,`$ (81)
$`\mathrm{}_2`$ $`=`$ $`\left|\begin{array}{cccc}\psi _1^{(1)}& \varphi ^{(1)}& \psi _3^{(1)}& \psi _4^{(1)}\\ \psi _1^{(2)}& \varphi ^{(2)}& \psi _3^{(2)}& \psi _4^{(2)}\\ \psi _1^{(3)}& \varphi ^{(3)}& \psi _3^{(3)}& \psi _4^{(3)}\\ \psi _1^{(4)}& \varphi ^{(4)}& \psi _3^{(4)}& \psi _4^{(4)}\end{array}\right|`$ (86)
Functions $`\psi _1^{(i)}`$, $`\psi _2^{(i)}`$, $`\varphi ^{(i)}`$ are defined by the same formulas (35-37) and (39) ($`j=1,2`$) together with solutions (49)-(57). Because of this, solution (41) may depend on four arbitrary functions of single variable (compare with the paragraph below the eq.(57)).

Functions $`\psi _3^{(i)}`$ and $`\psi _4^{(i)}`$ are solutions of the equations (38,39) with $`j=3,4`$:

$`\psi _{3}^{(i)}{}_{x}{}^{}=\alpha _{53}\psi _3^{(i)}+\alpha _{54}\psi _4^{(i)}+\alpha _{52}\psi _2^{(i)},\psi _{3}^{(i)}{}_{y}{}^{}=\alpha _{63}\psi _3^{(i)}+\alpha _{64}\psi _4^{(i)}+\alpha _{62}\psi _2^{(i)},`$ (87)
$`\psi _{4}^{(i)}{}_{x}{}^{}=\alpha _{73}\psi _3^{(i)}+\alpha _{74}\psi _4^{(i)}+\alpha _{72}\psi _2^{(i)},\psi _{4}^{(i)}{}_{y}{}^{}=\alpha _{83}\psi _3^{(i)}+\alpha _{84}\psi _4^{(i)}+\alpha _{82}\psi _2^{(i)},`$ (88)
$`l_k\psi _j^{(i)}=0,j=3,4,k=1,2.`$ (89)
This is a system of linear nonhomogeneous equations with solution in the following general form:

$`\psi _3^{(i)}=\psi _{30}^{(i)}+\psi _{3p}^{(i)},\psi _4^{(i)}=\psi _{40}^{(i)}+\psi _{4p}^{(i)},`$ (90)
where $`\psi _{j0}^{(i)}`$ and $`\psi _{jp}^{(i)}`$ are general solution of homogeneous system, related with the system (87-89), and particular solution of nonhomogeneous system (87-89) respectively.

The general solution of homogeneous system reads:

$`\psi _{30}^{(i)}`$ $`=`$ $`p_1^{(i)}r_3^{(1)}e^{\stackrel{~}{k}_1x+\stackrel{~}{q}_1y+\stackrel{~}{\omega }_1t+\stackrel{~}{\nu }_1\tau }+p_2^{(i)}r_3^{(2)}e^{\stackrel{~}{k}_2x+\stackrel{~}{q}_2y+\stackrel{~}{\omega }_2t+\stackrel{~}{\nu }_2\tau },`$ (91)
$`\psi _{40}^{(i)}`$ $`=`$ $`p_1^{(i)}r_4^{(1)}e^{\stackrel{~}{k}_1x+\stackrel{~}{q}_1y+\stackrel{~}{\omega }_1t+\stackrel{~}{\nu }_1\tau }+p_2^{(i)}r_4^{(2)}e^{\stackrel{~}{k}_2x+\stackrel{~}{q}_2y+\stackrel{~}{\omega }_2t+\stackrel{~}{\nu }_2\tau },`$ (92)
$`r_4^{(i)}`$ $`=`$ $`{\displaystyle \frac{(\stackrel{~}{k}_i\alpha _{53})r_3^{(i)}}{\alpha _{54}}},`$ (93)
$`\stackrel{~}{q}_j`$ $`=`$ $`{\displaystyle \frac{1}{\alpha _{54}}}(\alpha _{74}\stackrel{~}{k}_j+\alpha _{54}\alpha _{73}\alpha _{53}\alpha _{74}),j=1,2,`$ (94)
$`\stackrel{~}{k}_{1,2}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\alpha _{53}+\alpha _{64}\pm `$
$`{\displaystyle \frac{1}{\sqrt{\alpha _{74}}}}\sqrt{\alpha _{53}^2\alpha _{74}2\alpha _{53}\alpha _{64}\alpha _{74}+\alpha _{64}^2\alpha _{74}+4\alpha _{54}^2\alpha _{83}}),`$
where $`\stackrel{~}{\omega }_j`$ and $`\stackrel{~}{\nu }_j`$ are related with $`\stackrel{~}{k}_j`$ and $`\stackrel{~}{q}_j`$ by the equations

$`\stackrel{~}{\omega }_j=\stackrel{~}{k}_j^2\stackrel{~}{k}_j\stackrel{~}{q}_j\gamma _1,\stackrel{~}{\nu }_j=\stackrel{~}{q}_j^2\stackrel{~}{k}_j\stackrel{~}{q}_j\gamma _2,`$ (96)
which follows from the eqs.(89).

The particular solution can be taken in the form (in virtue of the eq.(50))

$`\psi _{3p}^{(i)}`$ $`=`$ $`{\displaystyle \mathrm{\Gamma }_1^{(i)}(k)c_2^{(i)}(k)\mathrm{exp}((x+\beta _2/\beta _1y)k+\omega (k)t+\nu (k)\tau )𝑑k},`$ (97)
$`\psi _{4p}^{(i)}`$ $`=`$ $`{\displaystyle \mathrm{\Gamma }_2^{(i)}(k)c_2^{(i)}(k)\mathrm{exp}((x+\beta _2/\beta _1y)k+\omega (k)t+\nu (k)\tau )𝑑k},`$
with $`\omega (k)`$ and $`\nu (k)`$ given by the eq.(55,56). Substitution of the eqs.(97) into the eqs.(87,88) results in the following expressions for $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$ together with additional relations among parameters $`\alpha _{ij}`$:

$`\mathrm{\Gamma }_1(k)`$ $`=`$ $`{\displaystyle \frac{\alpha _{52}^2}{\alpha _{52}k\alpha _{52}\alpha _{53}\alpha _{54}\alpha _{62}}},\mathrm{\Gamma }_2(k)={\displaystyle \frac{\alpha _{62}}{\alpha _{52}}}\mathrm{\Gamma }_1(k),`$ (98)
$`\alpha _{52}`$ $`=`$ $`{\displaystyle \frac{\alpha _{54}\alpha _{62}(\alpha _{54}\beta _2\alpha _{74}\beta _1)}{(\alpha _{53}\alpha _{74}\beta _1\alpha _{64}\alpha _{74}\beta _1+\alpha _{54}\alpha _{84}\beta _1\alpha _{53}\alpha _{54}\beta _2)}},`$
$`\alpha _{72}`$ $`=`$ $`{\displaystyle \frac{\alpha _{52}\beta _2}{\beta _1}},\alpha _{82}={\displaystyle \frac{\alpha _{62}\beta _2}{\beta _1}},`$ (99)
$`\alpha _{83}`$ $`=`$ $`{\displaystyle \frac{\alpha _{74}(\alpha _{84}\beta _1\alpha _{74}\beta _2)}{\alpha _{54}(\alpha _{54}\beta _2\alpha _{74}\beta _1)^2}}\times `$
$`(\alpha _{53}\alpha _{74}\beta _1\alpha _{64}\alpha _{74}\beta _1+\alpha _{54}\alpha _{84}\beta _1\alpha _{53}\alpha _{54}\beta _2)`$
As an example, let us choose the following expressions for $`\psi _j^{(i)}`$ and $`\varphi ^{(i)}`$:

$`\psi _1^{(i)}`$ $`=`$ $`s_ic_1^{(i)}e^{k_ix+q_iy+\omega _it+\nu _i\tau },i=1,2,3,`$ (101)
$`\psi _1^{(4)}`$ $`=`$ $`s_4c_1^{(4)}e^{k_4x+q_4y+\omega _4t+\nu _4\tau }+s_5c_1^{(5)}e^{k_5x+q_5y+\omega _5t+\nu _5\tau },`$ (102)
$`\psi _2^{(i)}`$ $`=`$ $`s_ie^{k_ix+q_iy+\omega _it+\nu _i\tau },i=1,2,3,`$ (103)
$`\psi _2^{(4)}`$ $`=`$ $`s_4e^{k_4x+q_4y+\omega _4t+\nu _4\tau }+s_5e^{k_5x+q_5y+\omega _5t+\nu _5\tau },`$ (104)
$`\varphi ^{(i)}`$ $`=`$ $`s_ic^{(i)}e^{k_ix+q_iy+\omega _it+\nu _i\tau },i=1,2,3,`$ (105)
$`\varphi ^{(4)}`$ $`=`$ $`s_4c^{(4)}e^{k_4x+q_4y+\omega _4t+\nu _4\tau }+s_5c^{(5)}e^{k_5x+q_5y+\omega _5t+\nu _5\tau },`$ (106)
$`\psi _3^{(1)}`$ $`=`$ $`p_1e^{\stackrel{~}{k}_1x+\stackrel{~}{q}_1y+\stackrel{~}{\omega }_1t+\stackrel{~}{\nu }_1\tau }+\mathrm{\Gamma }_1(k_1)s_1e^{k_1x+q_1y+\omega _1t+\nu _1\tau },`$ (107)
$`\psi _3^{(2)}`$ $`=`$ $`p_2e^{\stackrel{~}{k}_2x+\stackrel{~}{q}_2y+\stackrel{~}{\omega }_2t+\stackrel{~}{\nu }_2\tau }+\mathrm{\Gamma }_1(k_2)s_2e^{k_2x+q_2y+\omega _2t+\nu _2\tau },`$ (108)
$`\psi _3^{(3)}`$ $`=`$ $`\mathrm{\Gamma }_1(k_3)s_3e^{k_3x+q_3y+\omega _3t+\nu _3\tau },`$ (109)
$`\psi _3^{(4)}`$ $`=`$ $`\mathrm{\Gamma }_1(k_4)s_4e^{k_4x+q_4y+\omega _4t+\nu _4\tau }+\mathrm{\Gamma }_1(k_5)s_5e^{k_5x+q_5y+\omega _5t+\nu _5\tau },`$ (110)
$`\psi _4^{(1)}`$ $`=`$ $`p_1r_4^{(1)}e^{\stackrel{~}{k}_1x+\stackrel{~}{q}_1y+\stackrel{~}{\omega }_1t+\stackrel{~}{\nu }_1\tau }+\mathrm{\Gamma }_2(k_1)s_1e^{k_1x+q_1y+\omega _1t+\nu _1\tau },`$ (111)
$`\psi _4^{(2)}`$ $`=`$ $`p_2r_4^{(2)}e^{\stackrel{~}{k}_2x+\stackrel{~}{q}_2y+\stackrel{~}{\omega }_2t+\stackrel{~}{\nu }_2\tau }+\mathrm{\Gamma }_2(k_2)s_2e^{k_2x+q_2y+\omega _2t+\nu _2\tau },`$ (112)
$`\psi _4^{(3)}`$ $`=`$ $`\mathrm{\Gamma }_2(k_3)s_3e^{k_3x+q_3y+\omega _3t+\nu _3\tau },`$ (113)
$`\psi _4^{(4)}`$ $`=`$ $`\mathrm{\Gamma }_2(k_4)s_4e^{k_4x+q_4y+\omega _4t+\nu _4\tau }+\mathrm{\Gamma }_2(k_5)s_5e^{k_5x+q_5y+\omega _5t+\nu _5\tau }.`$ (114)
where $`c^{(n)}c^{(n)}(k_n)`$, $`c_1^{(n)}c_1^{(n)}(k_n)`$ are given by the eqs.(52,53) with $`c_2^{(n)}=1`$; $`r_4^{(n)}`$ are given by the eq.(93) with $`r_3^{(n)}=1`$; parameters $`k_n,q_n,\omega _n,\nu _n`$ are mutually related by the eqs.(54-56); relationship among parameters $`\stackrel{~}{k}_n,\stackrel{~}{q}_n,\stackrel{~}{\omega }_n,\stackrel{~}{\nu }_n`$ is given by the eqs.(94-96); $`\mathrm{\Gamma }_n`$ are given by the eqs.(98); $`s_n`$ and $`p_n`$ are arbitrary parameters.

General expressions for the functions $`u_1`$ and $`u_2`$ are very complicated. As an example, let us fix the following list of parameters:

$`k_1`$ $`=`$ $`1,k_2=1,k_3=\mathrm{\hspace{0.33em}2},k_4=2,k_5=\mathrm{\hspace{0.33em}0},`$ (115)
$`\beta _1`$ $`=`$ $`1,\beta _2=3,\gamma _1=\gamma _2=\delta _{12}=\mathrm{\hspace{0.33em}1},`$ (116)
$`\alpha _{12}`$ $`=`$ $`\alpha _{32}=\alpha _{53}=\alpha _{62}=\alpha _{74}=\mathrm{\hspace{0.33em}1},`$ (117)
$`\alpha _{54}`$ $`=`$ $`\alpha _{84}=\mathrm{\hspace{0.33em}2},\alpha _{64}=1.`$ (118)
Then, formulas (41) read:

$`u_1`$ $`=`$ $`{\displaystyle \frac{45}{\stackrel{~}{\mathrm{}}}}(s_5e^{2(x+12\tau )}2s_4e^{2(3y+4t)})\times `$
$`\left(343s_2e^{\frac{6}{7}y+\frac{48}{49}t}+76p_2e^{\frac{2}{7}x+\frac{144}{49}\tau }\right),`$
$`u_2`$ $`=`$ $`{\displaystyle \frac{6}{\stackrel{~}{\mathrm{}}}}(15680s_2s_4e^{\frac{48}{7}y+\frac{440}{49}t}+8967s_2s_5e^{2x+\frac{6}{7}y+\frac{48}{49}t+24\tau }+`$
$`4560s_4p_2e^{\frac{2}{7}x+6y+8t+\frac{144}{49}\tau }+1140s_5p_2e^{\frac{16}{7}x+\frac{1320}{49}\tau }),`$
$`\stackrel{~}{\mathrm{}}=24010s_2s_4e^{\frac{48}{7}y+\frac{440}{49}t}+21609s_2s_5e^{2x+\frac{6}{7}y+\frac{48}{49}t+24\tau }+`$
$`6840s_4p_2e^{\frac{2}{7}x+6y+8t+\frac{144}{49}\tau }+3420s_5p_2e^{\frac{16}{7}x+\frac{1320}{49}\tau }.`$
This solution has no singularities if

$`\text{sign}(s_2s_4)=\text{sign}(s_2s_5)=\text{sign}(s_4p_2)=\text{sign}(s_5p_2)`$ (121)
The nonlinear system keeps the same form (69,70). The change of variables (71) is not effective in this case, so that the nonlinear system is essentially (2+2)-dimensional. If one of parameters $`x`$, $`y`$, $`t`$ or $`\tau `$ is fixed, then functions $`u_1`$ and $`u_2`$ represent (2+1)-dimensional inelastic 3-kinks interaction.

## 4 Conclusions

The represented method is another way of using the algebraic system of equations for analysis of nonlinear PDE. In comparison with method considered in , it supplies wide class of particular solutions for nonlinear PDE which may be solved by this method. In general, nonlinear PDE considered in this paper, are not completely integrable in classical sence. For instance, attempt to introduce commuting flows by the standard way leads, generally speaking, to new constraints on the manifold of available solutions to nonlinear PDE.

There is no formal restriction on the linear differential operators $`L_1`$ and $`L_2`$ (see eqs.(20)). This makes the algorithm fluxible, but its application field is not defined yet. An interesting problem is construction of commuting flows on submanifold of particular solutions to the given nonlinear PDE.