Energy spectrum of localized quasiparticles renormalized by multi-phonon processes at finite temperature

The theory of renormalized energy spectrum of localized quasi-particle interacting with polarization phonons at finite temperature is developed within the Feynman-Pines diagram technique. The created computer program effectively takes into account multi-phonon processes, exactly defining all diagrams of mass operator together with their analytical expressions in arbitrary order over the coupling constant. Now it is possible to separate the pole and non-pole mass operator terms and perform a partial summing of their main terms. The renormalized spectrum of the system is obtained within the solution of dispersion equation in the vicinity of the main state where the high- and low-energy complexes of bound states are observed. The properties of the spectrum are analyzed depending on the coupling constant and the temperature.


Introduction
Obviously, there are many reasons that despite of the long period of its successful development, the theory of the interaction of quasi-particles (electrons, excitons and so on) with quantized fields (phonons, photons) remains in the center of physicists attention. One of them is that the new phenomena permanently discovered in physics (superfluidity, superconductivity, high-temperature superconductivity and so on) in the classic 3d and 2d structures [1][2][3] and, further [4,5] in low-dimensional heterostructures (in particular, nano-structures) demanded a deep understanding of the specific interaction between quasi-particles proper or with different fields and urged the development of mathematical methods in theoretical physics.
A fundamental progress in the theory of interaction of quasi-particles-quantized fields in condensed materials has been achieved late in the twentieth century when the methods of quantum field theory were developed [6,7], in particular, the Feynman diagram technique in the method of thermodynamic Green's functions [2,8,9]. However, this universal and powerful mathematical method is not devoid of the difficulties at the stage of its direct application to the specific problems where the perturbation method or the variational one is not applicable. For instance, the difficulties associated with the peculiar problem of sign [10] appear when the excited hybrid states of the system [11] are studied and it is necessary to take into account the multi-phonon processes. In any representation, except the Matsubara one, the expanded Green's function in the representation of interaction contains complex terms with the sign that changes in general case. Hence, the mathematical structure of mass operator (MO) ranges becomes complicated. Therefore, in order to reliably define the spectrum of a quasi-particle interacting with phonons in a wide range of energies, it is not enough to account for the finite number of diagrams, but it is necessary to perform a partial summing of an infinite number of diagrams, if not all of them, then at least the main ones.
The long period of a lack of reliable and accurate data on the main and excited states in macroscopic structures has stimulated a search for new approaches to the mathematical methods of their research. This, in particular, contributed to the creation of diagrammatic Monte Carlo method, which together with the method of stochastic optimization made it possible to calculate the Green's functions of quasi-particles in different models [10,[12][13][14][15][16][17][18][19][20][21] almost exactly and to get reliable data where they were incomplete or very controversial [10], like in the concept of the relaxed excited states. This method does not require any explicit "circumcision" of the orders of the ranges of Green's functions [10] and satisfies the central limit theorem, since, using a sufficient reserve of memory of modern computers, one can always achieve an insignificantly small effect of the system fluctuations and obtain an almost exact result.
Modern computers significantly expanded the scope of applying the Feynman diagram technique in the method of Green's functions for the problems of quasi-particles-phonons interaction, especially in the cases where it is necessary to consider the multi-phonon processes [3,10,11]. Recently, using a modified Feynman-Pines diagram technique, computer calculations of the leading classes of MO diagrams were performed [22]. Herein, a physically correct picture was established, which was missing earlier, for the ground and excited states of Frohlich polaron in the regime of a weak electron-phonon coupling at A renormalized spectrum of a localized quasiparticle weakly interacting with phonons at T 0 K was obtained in [23] using the Feynman-Pines diagram technique in the general model with the Frohlich Hamiltonian. It was shown that the new complexes of bound states existed in the system besides the ground and bound states which were known earlier [1,5] from the simplified model with the additional condition for the operators of quasi-particles. The energy spectrum was calculated using the MO, which contained all diagrams till the third order. Of course, this was not enough to identify the leading diagrams for their partial summing. However, it was impossible to hope for the analytical calculation of high-order diagrams without a computer, since their number in the n-th order was proportional to n · 2 n (2n − 1)!!.
In this paper, we use a new computer program for analytical calculation of MO in such a high order over the coupling constant, which is limited by the computer resource only. Although we take into account all diagrams till the tens order, using a personal computer, this is enough for a partial summing of the main MO diagrams. As a result, we confirm the existence of previously known and new states of a quasiparticle interacting with phonons and analyze the temperature dependence of the spectrum.
2. Hamiltonian of the system. Multi-phonon processes and structure of complete MO at T 0 K Localized dispersionless quasiparticles (excitons, impurities and so on) interacting with dispersionless polarization phonons are described by Frohlich Hamiltonian, like in [23] where E and Ω are the energies of quasi-particles and phonons, respectively, which can be arbitrary in the typical ranges (E ≈ 100 ÷ 1000 meV, Ω ≈ 20 ÷ 100 meV) for the solid states and nano-structures. It turns out that the form of ϕ(q) binding function is irrelevant because the energy of quasi-particle-phonon interaction is uniquely characterized by a constant α independent of ì q , also arbitrary in the natural range, since the problem becomes zero-dimensional. The operators of second quantization of quasi-particles ( a ì k , a + ì k ) and phonons ( b ì q , b + ì q ) satisfy the Bose commutative relationships. The renormalized spectrum of quasi-particles interacting with phonons at an arbitrary temperature (T), like in [23], is obtained from the poles of Fourier image of quasi-particle Green's function G( ì k, ω ′ ), which through the Dyson equation is related with MO M( ì k, ω ′ ). It is calculated using the Feynman-Pines diagram technique when the 43706-2 concentration of localized quasi-particles is small. Introducing the dimensionless variables and values the Dyson equation (2) is rewritten as As far as the problem becomes "zero-dimensional", the complete MO m(ξ) is defined by the infinite sum of indexed diagrams (figure 1). The simple rules for these diagrams presented in [23], allow their definite programming and numerical calculation of all topologically unequal diagrams and their analytical expressions till such a high order over the powers of a coupling constant (α), which is limited by a computer resource only (appendix A). We should note that in this paper all diagrams and their analytical expressions for the MO terms till the tenth order are obtained. Although in figure 1, all first terms of MO, including the third order, are presented for illustration and for further understanding.

43706-3
The rules of conformity between indexed diagrams and analytical expressions are simple: where number n is the sum of phonon lines with the arrows directed to the right (with the sign "−") and with the arrows directed to the left (with the sign "+") over the quasi-particle line (figure 1), ν = (e Ω/kT − 1) −1 is an average phonon occupation number. Thus, the analytical expression for the arbitrary diagram with indices is written as a product of contributions of all its lines and tops. For example, the analytical expressions are as follows: .
A detailed analysis of the structure of complete MO, figure 1, till high orders over the powers of α, ν, (1+ν) shows that it can be written in an analytical form which contains three types of terms. The terms with the symbol ">": m > k, 2k+n (ξ) ∼ α 2k+n ν k (1 + ν) k+n (k = 0, 1, . . . , n = 1, 2, . . . ) describe the processes of quasi-particle scattering accompanied by the prevailing creation of phonons, the terms with the symbol "<": m < k, 2k+n (ξ) ∼ α 2k+n ν k+n (1 + ν) k describe the processes of quasi-particle scattering accompanied by the prevailing annihilation of phonons and m (e) 2k (ξ) ∼ α 2k [ν(1 + ν)] k . Hence, the infinite ranges of complete MO terms can be written as In this formula, the upper m > u (ξ) and the lower m < u (ξ) terms of the complete MO are defined by the respective infinite ranges of the terms m ≷ 0,n (ξ), to which the diagrams without crossing phonon lines with arrays directed only to the right (>) and only to the left (<) correspond. These diagrams describe the unmixed (u) processes of quasi-particles scattering accompanied either only by the creation (>) or by annihilation (<) of phonons, respectively. After the partial summing, an exact analytical representation for these terms in the form of infinite chain fractions was obtained [23].
The rest terms of MO (7)

43706-4
Energy spectrum of localized quasiparticles describe the mixed (m) processes of quasi-particle scattering accompanied by an equal (e) number of the created and annihilated phonons by a prevailing number of the created phonons rather than annihilated ones and by a prevailing number of annihilated phonons rather than the created ones The functions ϕ (e)2k (ξ) = −ϕ (e)2k (−ξ) and ϕ > k,2l+n (ξ) = −ϕ < k,2l+n (−ξ) are polynomials, diagrammatic and analytical forms of which are numerically calculated within the computer using the diagram technique. Till the sixth order over the coupling constant (α), the expressions for these functions are presented in appendix B, though in this paper we took into account all diagrams till the tenth order.
Finally, the complete MO (7) now has the form where m > Contrary to the unmixed terms (9), the partial summing in the mixed ones (11)-(13) is, in general, impossible. However, using the newly revealed analytical structure, this becomes possible for the arbitrary energies, highlighting in these MO terms the pole and non-pole factors and expanding them in Taylor series.

Partial summing in mixed MO using the non-pole factors expanded in Taylor series
For a system having a weak quasi-particle-phonon interaction (α < 1) at the temperature when ν ≪ 1, we observe the range of energies (−2 < ξ < 2) where, according to physical considerations, one should expect the renormalized energies of the main state and its nearest satellite bound-with-phonon states. In order to obtain them, both in high-and low-energy regions, we perform a partial summing in mixed MO (11)- (13).
First of all, we observe m > (ξ) in the range 0 < ξ < 2, selecting pole and non-pole factors in it. Formula (11) proves that MO terms MO m > k, 2k+n (ξ) in this energy range can be submitted as a product

43706-5
where one can see a multiplier independent of ξ, pole multiplier and non-pole fractional-rational functions which in the vicinity of ξ = 1 are expanded into the Taylor series where F >(s) k,2n are the known numerical coefficients.
Substituting (17), (18) in (14) we see that now MO m > k, 2k+n (0 < ξ < 2) can be written as a sum of Similar analytics gives where F >(s) Substituting formulae (20)- (22) in (11)-(13) we select pole and non-pole terms in the mixed MO m (e) (0 < ξ < 2) and m > (0 < ξ < 2), m < (0 < ξ < 2) in the range (0 < ξ < 2). Taking into account formula (15), they are obtained in the range (−2 < ξ < 0). Thus, according to (10) The expression for m m (ξ) obtained in the form of separated pole and non-pole terms gives an opportunity to perform a partial summing of their infinite ranges. It is clear that since the maximal order (k max ) of functions f k,2k+n is limited by the maximal order of the accounted diagrams over α, then in the complete m m (ξ) one should take into account not more than k max (in our case k max = 10) mixed The sums over k quickly converge because k-th terms of this MO contain small multipliers α k and ν k . Avoiding sophisticated analytics, considering the symmetric relationship (15) and a similar structure of m k at all k values, we demonstrate the method of partial summing by the example of MO m 1 (ξ) in the range (0 < ξ < 2).
In table 1 where The main idea of partial summing of the main terms of pole MO is as follows: using the s-th column, The renormalized m (p)(s) 1 (ξ) of a higher order over s are obtained in the same way. Still relatively not complicated MO at s = 1 and s = 2 have the form k, 2k+n (ξ) linearly decreases. Thus, beginning from certain k and s, it is impossible to obtain an exact analytical expression for an arbitrary term of the respective infinite range and, thus, to perform a partial summing. Usually, this problem is solved with the help of rather powerful computers, which can essentially increase k max . However, the presence of two small parameters (α and ν) in the theory, brings us to the fact that the contribution of pole MO terms rapidly decreases when k and n increase. Thus, as further numerical calculations prove, the above mentioned circumstances have a weak effect on the magnitudes of the renormalized energies without changing the number of complexes of bound states in the region of energies studied.

43706-7
The non-pole terms of m (np) k (ξ) can be calculated using two methods: either within the respective tables (like table 1), performing a partial summing over s at fixed k and n, or, which is easier, using the obvious formulae

43706-8
Energy spectrum of localized quasiparticles  For example, the first three terms of non-pole MO at k = 1 are the following ones m (np) The rest terms of m (np) k,2k+n (ξ) (at k k max , n max ) are easily obtained using the formula (30), but we do not present them due to their sophisticated analytical form.
The developed theory defines a complete MO m(ξ), which is a real function of ξ. Thus, the poles of g(ξ) function (4) are found from the equation which gives a renormalized spectrum of the system in the range (−2 < ξ < 2).

Analysis of quasi-particle spectrum renormalized due to the multiphonon processes in the vicinity of its main state at finite temperature
The theory developed in the previous section makes it possible to calculate the renormalized energy spectrum of localized quasiparticle interacting with polarization phonons at T 0 K in a wide range of energies taking into account both unmixed and mixed multi-phonon processes. The proposed mathematical approach of MO analytical calculation differs from the one presented in [23], where the MO terms describing the unmixed processes [m u (ξ)] were completely taken into account using an exact partial summing, while the mixed ones [m m (ξ)] -in the first three orders over the coupling constant only. The "new", with respect to the simplified model [1], satellite bound states of the system were revealed in [23]. However, further, in this paper we should show that some properties of their energy spectrum turned out to be an artifact of the approximated MO [m m (ξ)] describing the mixed processes.
The renormalized energies of the main state (e 0 ) and complexes of bound states (e < 1 −1 , e >1 −1 , e < 1 +1 , e +1 , e >1 +1 , e >2 +1 ) in the range (−2 < ξ < 2) were obtained from dispersion equation (34). In order to describe its properties and establish the mechanisms of formation, we calculated MO and its main terms as functions of ξ at a different coupling constant (α) and average phonons occupation number [ν(T)].
The dependences of MO m(ξ) and its main partial terms m u (ξ) and m m (ξ) are presented too, in order to reveal the role of unmixed and mixed processes in the formation of the respective energy levels of the system. A more detail contribution of different multi-phonon processes into the formation of the spectrum, for instance, is shown in table 3. Analysis of figure 2 and table 3 gives a clear understanding of the formation of the spectrum of a localized quasiparticle renormalized due to multi-phonon processes in the regime of a weak coupling at a finite temperature (ν 0) and in the range of energies (−2 < ξ < 2). Figure 2 (a) and table 3 prove that at α = 0.16, ν = 0.05, the two stationary states, i.e., the ground one with the energy e 0 and its first high-energy repeat (e +1 ) are almost completely formed by the unmixed [m u (ξ) ≫ m m (ξ)] processes of quasi-particle-phonon interaction with the prevailing creation of phonons [m > u (ξ) ≫ m < u (ξ)]. These energies are shifted into the low-energy region by the magnitude ∼ α(1+ ν), so that the difference between them (∆e +1 = e +1 − e 0 = 0.982742) is close to the energy of one phonon.
In the high-energy region, one can see three stationary bound states, which are satellite to the first main one: the lower one with the energy e < 1 +1 , smaller than e +1 , and two upper states with the energies e >1 +1 , e >2

+1
bigger than e +1 . These satellite states are mainly formed by the mixed [m m (ξ) m u (ξ)] processes of quasi-particle-phonon interaction with the prevailing creation of phonons [m > m (ξ) ≫ m < m (ξ)]. We should note that the lower satellite state with the energy e < 1 +1 is also a stationary one, taking into account the multi-phonon processes. Table 3. Energy levels of the system and contributions of MO terms into their formation at α = 0.16;  The state with the energy (e < 1 +1 ) does not degenerate with the main one with the energy (e +1 ) and, thus, it does not decay contrary to the results of [23], where the positive contributions only of two-and three-phonon mixed processes (m m ) were taken into account. It was not sufficient to compensate the negative contributions produced by unmixed processes (m u ) into the formation of both energy levels when ν increased. Now, from table 3 one can clearly see that, among all high-energy levels, only the one with e < 1   figure 3. It is clear that in the regime of weak coupling, the temperature dependence of the spectrum is qualitatively similar. When ν (temperature) increases, the upper satellite low-and high-energy levels (e >1 −1 , e >1 +1 , e >2 +1 ) smoothly (quasi-linearly) shift into a high-energy region while the main level (e 0 ), its first phonon repeat (e +1 ) and both low-energy levels (e < 1 −1 , e < 1 +1 ) smoothly shift into the low-energy region of the spectrum. The magnitudes of the shifts of all levels slightly increase at a bigger α.

Main results and conclusions
The main result of this paper is that the theory of renormalized energy spectrum of a localized quasiparticle weakly interacting with polarization phonons at a finite temperature is developed within the Feynman-Pines diagram technique using the Frohlich model and taking into account the multi-phonon processes. Now it is possible to obtain the renormalized energies of the main state and the complexes of bound states.
The created computer program exactly calculates the analytical expressions for the diagrams of MO in  an arbitrary order over the coupling constant, which is limited by a computer resource only. The analysis of the diagrams till the tenth order made it possible to separate its pole and non-pole terms and perform a partial summing of their main terms. The analytical expressions for the general terms in the ranges of main diagrams are obtained and partially summed. As a result, a complete MO is obtained in the range −2 < ξ < 2 taking into account the multi-phonon processes. It is shown that a sufficient consideration of multi-phonon processes in MO provides a stationary energy spectrum of the system. Its properties are similar, in certain areas of the spectrum, and somewhere differ from that obtained both in the model with an additional condition for the quasi-particles operators [1] and in the Frohlich model, recently observed in [23], where the unmixed processes were accounted exactly while the mixed ones -only in the second and third orders of MO diagrams. It is established that, like in [23], contrary to the model used in [1], the spectrum of the system is not equidistant and contains "new" complexes of bound states. The latter are stationary and do not decay in a high-energy region as in [23] due to the above mentioned approximation in the mixed MO.
The renormalized energy of the main state and its high-energy phonon repeat are mainly formed by interaction accompanied by the creation of phonons in unmixed processes, while all satellite lowand high-energy states are mainly formed by the interaction accompanied by annihilation of phonons in mixed processes.
It is shown that in the regime of a weak coupling, when temperature increases, the upper satellite low-and high-energy levels (e >1 −1 , e >1 +1 , e >2 +1 ) smoothly quasi-linearly shift into a high-energy region while the main level (e 0 ), its first phonon repeat (e +1 ) and low-energy levels (e < 1 −1 , e < 1 +1 ) smoothly shift into the low-energy region of the spectrum.
The final conclusion is that using a computer program for analytical calculation of the high orders of MO diagrams, we effectively took into account the multi-phonon processes of interaction between the quasi-particle and polarization phonons, thus essentially clarifying the properties of the complexes of bound states near the main one. The developed theoretical approach together with a proper computer support can be successfully used for the calculation of the renormalized spectrum of this system in a wider range of energy and coupling constant. It can be used for the solution of the other problems where the quasi-particles interacting with phonons are investigated. On each recursive call we select the first available point as a starting point and enumerate all points on its right as the end points. For every variation received from the first level we shall have 2(N − 1) − 1 variations at the second level and for each of them, as before, we call the function recursively again.
At the third level in our N = 3 example, we shall always have only one possible option. Finally, once we have fixed this single possible option, we have fully finished the diagram for a given N.
In total, the third level of the recursive call will be called 2N − 1 + 2(N − 1) − 2(N − 2) + . . . 1 times giving us all possible and unique shapes of the diagrams.

A.2. Filtering out non-fully covered diagrams
Non-fully covered diagrams are not of interest to us and we just neglect them: 1 2 4 3 6 5

A.3. All possible directions
For each fully covered diagram, we enumerate all possible directions of dashed lines. This is implemented as a single loop from 0 to 2 N − 1. Each number in binary representation represents one possible variation of the lines direction where 0 means the left-hand direction and 1 means right-hand direction. In our example of N = 3, we iterate from 0 to 7 where: a. Binary representation of 0 is 000 which means all lines are directed to the left.