text:
# Mean parameter model for the Pekar-Fröhlich polaron in a multilayered heterostructure

## I Introduction

Quasi-two-dimensional (2D) systems have attracted a lot of attention during the last decade because of their practical realization. If a heterostructure is made of polar materials such as $`GaAs`$, $`InAs`$ etc., the electron-phonon interaction modifies the properties of the electron confined to a 2D-structure resulting in a shift of the binding energy and the effective band mass.

The polaron effects in a 2D electron gas have extensively been studied. At earlier stages the attention was paid to the properties of a polaron confined to an infinite thin 2D layer. The binding energy and the effective mass were calculated for a $`GaAs/Al_xGa_{1x}As`$ infinitely deep quantum well of a finite length. In these papers only the interaction with the bulk LO-phonon mode has been taken into account. Actually, LO-phonon modes are modified in a 2D layer (the so-called confined slab LO-phonon modes). Besides, there exist interface optical-phonon modes as well as half-space LO-phonon modes in a barrier material. For the review of these modes (also in complicated multi-layer structures) see the book by Pokatilov, Fomin and Beril and also more recent publications  of this group. The influence of the mentioned modes on polarons were studied in Refs. .

While different phonon modes were studied in details, the quantum well potential was supposed to be infinitely deep in the cited papers. On the other hand, the properties of the system would be quite different if a confining potential had a finite depth. Indeed, for an infinitely deep confining potential the binding energy is the monotonous function of the potential width which varies between limiting values $`E_{3D}^{(in)}=\alpha _{in}\mathrm{}\omega _{in}`$ for the three-dimensional (3D) space and $`E_{2D}^{(in)}=(\pi /2)\alpha _{in}\mathrm{}\omega _{in}`$, where $`\alpha _{in}`$ is the standard Fröhlich electron-phonon coupling constant and $`\omega _{in}`$ is the LO-phonons frequency for the quantum well material. If a particle is confined to a finite potential well, the limiting value of the binding energy should be the same at large width of the well. But when the width becomes too small, the energy level approaches the edge of the well, so that effectively the particle is spread over the 3D space. Thus, the limiting value of the binding energy should coincide with the 3D limiting value rather than with the 2D one. This means, the binding energy takes the value $`E_{3D}^{(out)}=\alpha _{out}\mathrm{}\omega _{out}`$ at small widths where the parameters $`\omega _{out}`$ and $`\alpha _{out}`$ are now related to the barrier material. The binding energy evidently has a peak at some intermediate value of the width if $`E_{3D}^{(out)}E_{3D}^{(in)}`$. If this is not the case, the existence of the peak should be checked in more detail.

Different rectangular quantum wells of a finite height have been investigated by Hai, Peeters and Devreese and Shi, Zhu et al. in the scope of the second order perturbation theory in powers of the electron-phonon coupling constant $`\alpha `$ with all phonon modes being incorporated. Peaks of the phonon induced energy shift and the polaron effective mass were found for some values of the confining potential widths.

In principle, the same approach can be used while dealing with a quantum well constructed of layers of different materials. But the problem becomes then too complicated because one has to take into account interface phonon modes at each frontier of different materials as well as quantized phonon modes inside each of the layers. The main goal of the present paper is to formulate a simplified model to take these effects into account and to deal with the effective confining potential and only one bulk phonon mode. We calculate polaron characteristics for the same rectangular quantum well as in Refs.  to compare the results. Another example is given of a quantum well of the finite depth for which the second-order correction due to the electron-phonon interaction can be calculated explicitly. Namely, we take the Rosen-Morse potential to confine electrons to a 2D-multilayered heterostructure and calculate the shift of the ground-state energy and the effective mass perturbatively, that is, in the weak-coupling limit. In contrast with the rectangular potential we should not worry about the correct including of all virtual states because the Green function is known analytically for this system.

## II Formulation of the model

Let us consider a quantum well in the $`z`$ direction constructed of the $`xy`$ plane layers of $`GaAs/Al_xGa_{1x}As`$. That is, the $`AlAs`$ mole fraction $`x`$ depends on the coordinate $`z`$: $`x=x(z)`$. The energy gap between different materials forms the confining potential $`V(z)`$ which serves us as the main entity. Given the potential $`V(z)`$, one can find the corresponding mole fractions $`x(z)`$ and a dependence on $`z`$ of any of the medium parameters (such as the electron band mass $`m(z)`$, phonon frequencies $`\omega (z)`$, dielectric constants $`\epsilon _0(z),\epsilon _{\mathrm{}}(z)`$, Fröhlich coupling constants $`\alpha (z)`$, etc.).

To avoid difficulties with mass mismatch in different layers we suggest to use a mean band mass $`m`$ which is common for all layers. Then we start with the following expression for the electronic part of the Hamiltonian:

$`H_{el}`$ $`=`$ $`H_{el,}+H_{el,},`$ (1)
$`H_{el,}`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{p}_{}^{\mathrm{\hspace{0.17em}2}}}{2m}},H_{el,}={\displaystyle \frac{p_z^{\mathrm{\hspace{0.17em}2}}}{2m}}+V(z),`$ (2)
where the electron mean band mass is defined by the relation

$`{\displaystyle \frac{1}{m}}={\displaystyle 𝑑z\frac{|\psi _1(z)|^2}{m(z)}}`$ (3)
and the ground state wave function $`\psi _1(z)`$ for the electron motion in $`z`$ direction is a solution to the Schrödinger equation

$`H_{el,}\psi _1=E_1\psi _1`$ (4)
with $`E_1`$ being a ground state energy. As the wave function $`\psi _1`$ also depends on the mean band mass $`m`$, the latter can be found as a self-consistent solution of Eqs. (3), (4).

In a similar way we define the free LO-phonon Hamiltonian

$`H_{ph}=\mathrm{}\omega _{\text{LO}}{\displaystyle \underset{\stackrel{}{k}}{}}a_\stackrel{}{k}^{}a_\stackrel{}{k},`$ (5)
where $`a_\stackrel{}{k}^{}(a_\stackrel{}{k})`$ are creation (annihilation) operators of a phonon with a wave vector $`\stackrel{}{k}`$, and mean frequency $`\omega _{\text{LO}}`$ can be found from the expression

$`\omega _{\text{LO}}={\displaystyle 𝑑z\omega (z)|\psi _1(z)|^2}.`$ (6)
Evidently, we have to address why the free phonon Hamiltonian is averaged with respect to the electron wave function. Our motivation is based on the fact that we are going to apply our model to calculate polaron effects. That is, our effective phonons will be considered only in a cloud around the electron, and the properties of this cloud depend on the electron position. So, in our model the effective phonons replace numerous phonon modes whose frequencies depend on the coordinate $`z`$ of the electron.

Finally, we accept the conventional form of the Hamiltonian describing the interaction of the electron with effective phonons:

$`H_{int}={\displaystyle \underset{\stackrel{}{k}}{}}\left(a_\stackrel{}{k}V_\stackrel{}{k}e^{i\stackrel{}{k}\stackrel{}{r}}+a_\stackrel{}{k}^{}V_\stackrel{}{k}^{}e^{i\stackrel{}{k}\stackrel{}{r}}\right),`$ (7)
where the Fourier transforms of the electron-phonon interaction potential are specified as follows:

$`V_\stackrel{}{k}=i\mathrm{}\omega _{\text{LO}}\left({\displaystyle \frac{4\pi \alpha }{Vk^2}}\sqrt{{\displaystyle \frac{\mathrm{}}{2m\omega _{\text{LO}}}}}\right)^{1/2}.`$ (8)
Here the use is made of a mean Fröhlich coupling constant $`\alpha `$ which can be found from the relation

$`\sqrt{\alpha }={\displaystyle 𝑑z|\psi _1(z)|^2\frac{\omega (z)}{\omega _{\text{LO}}}\left(\alpha (z)\sqrt{\frac{m\omega _{\text{LO}}}{m(z)\omega (z)}}\right)^{1/2}}.`$ (9)
Note that we define the mean parameters in Eqs. (3), (6), (9) according to the way they enter the Hamiltonian.

Thus, we describe a complicated multilayered heterostructure by the Hamiltonian

$`H=H_{el}+H_{ph}+H_{int}`$ (10)
with the bulk phonon mode only which inhabits an effective medium with mean characteristics defined above. The details of the heterostructure are taken into account by the confining potential $`V(z)`$.

Performing a unitary transformation $`HH^{}=U^1HU`$ with the operator

$`U=\mathrm{exp}\left[i\stackrel{}{r}_{}{\displaystyle \underset{\stackrel{}{k}}{}}\stackrel{}{k}_{}a_\stackrel{}{k}^{}a_\stackrel{}{k}\right],`$ (11)
we arrive at the Hamiltonian

$`H^{}=H_{el,}^{}+H_{el,}+H_{ph}+H_{int}^{},`$ (13)
$`H_{el,}^{}`$ $`=`$ $`{\displaystyle \frac{1}{2m}}\left(\stackrel{}{p}_{}\mathrm{}{\displaystyle \underset{\stackrel{}{k}}{}}\stackrel{}{k}_{}a_\stackrel{}{k}^{}a_\stackrel{}{k}\right)^2,`$ (14)
$`H_{int}^{}={\displaystyle \underset{\stackrel{}{k}}{}}\left(a_\stackrel{}{k}V_\stackrel{}{k}e^{ik_zz}+a_\stackrel{}{k}^{}V_\stackrel{}{k}^{}e^{ik_zz}\right).`$ (15)
The quantity $`\stackrel{}{p}_{}`$ is a c-number corresponding to the conserved momentum in the $`xy`$ plane and the Hamiltonians $`H_{el,},H_{ph}`$ are defined by Eqs. (2), (5), respectively.

Keeping in mind the smallness of the electron-phonon coupling constant $`\alpha `$ for most of the materials, we calculate the second-order correction to the unperturbed Hamiltonian $`H_0^{}=H_{el,}^{}+H_{el,}+H_{ph}`$ (note that the quantum-mechanical first-order correction is equal to zero). The unperturbed energy levels are given by the expression

$`E(\stackrel{}{p}_{},n_\stackrel{}{k},N)`$ $`=`$ $`{\displaystyle \frac{1}{2m}}\left(\stackrel{}{p}_{}\mathrm{}{\displaystyle \underset{\stackrel{}{k}}{}}\stackrel{}{k}_{}n_\stackrel{}{k}\right)^2+`$ (17)
$`\mathrm{}\omega _{\text{LO}}{\displaystyle \underset{\stackrel{}{k}}{}}n_\stackrel{}{k}+E_N,`$
where $`n_\stackrel{}{k}`$ is the number of phonons with the wave vector $`\stackrel{}{k}`$. The energy $`E_N`$ is the $`N`$-th energy level of the one-dimensional system $`H_{el,}`$ of Eq. (2). Here $`N`$ is the corresponding quantum number not necessarily a discrete one: it stands for both the quantum number $`n`$ which varies from 1 to $`n_{max}`$ and the wave vector $`q`$ of the continuous spectrum states.

The wave functions of the unperturbed Hamiltonian $`H_0^{}`$ are given by the direct product

$`|\stackrel{}{p}_{};n_\stackrel{}{k},N=|n_\stackrel{}{k}|N`$ (18)
of the corresponding wave functions of different terms in $`H_0^{}`$.

Because of the structure of the interaction term $`H_{int}^{}`$ only intermediate states with one phonon contribute to the second order correction to the ground-state energy. The latter is then given by the expression

$`\mathrm{\Delta }_2E(\stackrel{}{p}_{})=`$ (19)
$`{\displaystyle \underset{N,\stackrel{}{k}}{}}{\displaystyle \frac{|V_\stackrel{}{k}|^2|G(N,k_z)|^2}{E_N+\mathrm{}\omega _{\text{LO}}+{\displaystyle \frac{(\stackrel{}{p}_{}\mathrm{}\stackrel{}{k}_{})^2\stackrel{}{p}_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}}}{2m}}E_1}},`$ (20)
where

$`G(N,k_z)={\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑z\psi _N(z)\psi _1(z)e^{ik_zz}`$ (21)
and $`\psi _N(z)`$ are the wave functions of the Hamiltonian $`H_{el,}`$ in Eq. (2). The concrete application of these formulae is given in the following section.

## III Rectangular Potential

### A Medium mean characteristics

As an example we now consider the rectangular confining potential

$`V(z)=\{\begin{array}{cc}0,\hfill & |z|L/2\hfill \\ V_0,\hfill & |z|>L/2\hfill \end{array}`$ (24)
and

$`m(z)=\{\begin{array}{cc}m_{in},\hfill & |z|L/2\hfill \\ m_{out},\hfill & |z|>L/2\hfill \end{array}`$ (27)
with $`m_{in}(m_{out})`$ being the electron band masses in the well (barrier) material, respectively. For concreteness we assume $`GaAs`$ to be the quantum well material and $`Al_xGa_{1x}As`$ to be the barrier material.

Symmetrical wave functions of the discrete spectrum in the rectangular quantum well with the mean band mass $`m`$ take the form

$`\psi _{s,n}=N_{s,n}\{\begin{array}{cc}\mathrm{cos}q_nz,\hfill & |z|L/2\hfill \\ \mathrm{cos}{\displaystyle \frac{q_nL}{2}}e^{p_n(|z|L/2)},\hfill & |z|>L/2,\hfill \end{array}`$ (30)
where

$`p_n=\sqrt{q_{max}^2q_n^2},q_{max}^2={\displaystyle \frac{2mV_0}{\mathrm{}^2}}`$ (31)
and the normalization constant is given by

$`N_{s,n}=\sqrt{{\displaystyle \frac{2p_n}{p_nL+2}}}.`$ (32)
Antisymmetrical wave functions of the discrete spectrum take the form

$`\psi _{a,n}=N_{s,n}\{\begin{array}{cc}\mathrm{sin}q_nz,\hfill & |z|L/2\hfill \\ \text{sgn}(z)\mathrm{sin}{\displaystyle \frac{q_nL}{2}}e^{p_n(|z|L/2)},\hfill & |z|>L/2\hfill \end{array}`$ (35)
(36)
with the same normalization constant given by Eq. (32).

The total number $`n_{max}`$ of the discrete energy levels is given by the expression

$`n_{max}=1+\left[{\displaystyle \frac{q_{max}L}{\pi }}\right],`$ (37)
where $`[A]`$ is an integer part of a number $`A`$. The expression for the discrete energy levels reads as follows:

$`{\displaystyle \frac{q_nL}{2}}=\text{atan}\sqrt{{\displaystyle \frac{q_{max}^2}{q_n^2}}1}+{\displaystyle \frac{\pi (n1)}{2}},n=1,2,\mathrm{}.`$ (38)
Energies with odd (even) $`n`$ correspond to the symmetrical (antisymmetrical) wave functions.

The energy $`E_q=\mathrm{}^2q^2/2m`$ of the continuous spectrum state depends on the wave vector $`q`$. The corresponding symmetrical wave functions are as follows:

$`\psi _{s,q}={\displaystyle \frac{N_{s,q}}{\sqrt{L_z}}}\{\begin{array}{cc}p\mathrm{cos}qz,\hfill & |z|L/2,\hfill \\ p\mathrm{cos}{\displaystyle \frac{qL}{2}}\mathrm{cos}p(|z|L/2)\hfill & \\ q\mathrm{sin}{\displaystyle \frac{qL}{2}}\mathrm{sin}p(|z|L/2),\hfill & |z|>L/2,\hfill \end{array}`$ (42)
(43)
where

$`p=\sqrt{q^2q_{max}^2}`$ (44)
and $`L_z`$ is the (infinite) size of the system in the $`z`$ direction. The normalization constant is given by the expression

$`N_{s,q}=\sqrt{{\displaystyle \frac{2}{p^2\mathrm{cos}^2{\displaystyle \frac{qL}{2}}+q^2\mathrm{sin}^2{\displaystyle \frac{qL}{2}}}}}.`$ (45)
The antisymmetrical wave functions are as follows:

$`\psi _{a,q}={\displaystyle \frac{N_{a,q}\text{sgn}(z)}{\sqrt{L_z}}}\{\begin{array}{cc}p\mathrm{sin}q|z|,|z|L/2,\hfill & \\ p\mathrm{sin}{\displaystyle \frac{qL}{2}}\mathrm{cos}p(|z|L/2)+\hfill & \\ q\mathrm{cos}{\displaystyle \frac{qL}{2}}\mathrm{sin}p(|z|L/2),\hfill & \\ |z|>L/2,\hfill & \end{array}`$ (50)
(51)
where the normalization constant is given by the expression

$`N_{a,q}=\sqrt{{\displaystyle \frac{2}{p^2\mathrm{sin}^2{\displaystyle \frac{qL}{2}}+q^2\mathrm{cos}^2{\displaystyle \frac{qL}{2}}}}}.`$ (52)
The electron mean band mass is defined as

$`{\displaystyle \frac{1}{m}}`$ $`=`$ $`{\displaystyle \frac{W_{in}}{m_{in}}}+{\displaystyle \frac{W_{out}}{m_{out}}}`$ (53)
$`m`$ $`=`$ $`{\displaystyle \frac{m_{in}m_{out}}{W_{in}m_{out}+(1W_{in})m_{in}}},`$ (54)
where $`W_{in}`$ and $`W_{out}=1W_{in}`$ are probabilities to find the electron inside (outside) the quantum well. The expression for $`W_{in}`$ follows from Eq. (30)

$`W_{in}=2N_{s,1}^2{\displaystyle \underset{0}{\overset{L/2}{}}}𝑑z\mathrm{cos}^2q_1z=1{\displaystyle \frac{(q_1/q_{max})^2}{1+p_1L/2}},`$ (55)
where $`q_1`$ is a solution to Eq. (38) for the ground-state ($`n=1`$).

To finish this subsection, we note that the exact energy levels in the rectangular potential with different masses $`m_{in}`$ and $`m_{out}`$ calculated for the $`GaAs/Al_xGa_{1x}As`$ heterostructure practically coincide with the levels obtained with the electron mean band mass $`m`$. To obtain an inner criterion of the validity of the anzatz concerning the mean band mass we notice that the particle being on lowest energy levels is located mostly inside the well which means that its band mass is almost coincide with $`m_{in}`$.

One can await the largest discrepancy for a level near the potential edge. The $`n`$-th discrete level appears at the width $`L=L_n^{(av)}`$, where

$`L_n^{(av)}=\pi (n1){\displaystyle \frac{\mathrm{}}{\sqrt{2mV_0}}}={\displaystyle \frac{\pi (n1)}{q_{\mathrm{max}}}},`$ (56)
and the analogous width for the exact solution reads as follows:

$`L_n^{(ex)}=\pi (n1){\displaystyle \frac{\mathrm{}}{\sqrt{2m_{in}V_0}}}.`$ (57)
Thus, the ratio

$`{\displaystyle \frac{L_n^{(av)}}{L_n^{(ex)}}}=\sqrt{{\displaystyle \frac{m_{in}}{m}}}=\sqrt{W_{in}+(1W_{in}){\displaystyle \frac{m_{in}}{m_{out}}}}.`$ (58)
can serve us as the numerical criterion of the validity of the anzatz. The largest discrepancy happens at $`n=2`$ and in this case Eqs. (31), (38) (55) lead to the following expression:

$`{\displaystyle \frac{L^{(av)}}{L^{(ex)}}}=\sqrt{0.844+0.156{\displaystyle \frac{m_{in}}{m_{out}}}}.`$ (59)
Note that numerical coefficients here do not depend on the material parameters. For the $`GaAs/Al_{0.3}Ga_{0.7}As`$ quantum well we have $`m_{in}/m_{out}0.7`$ and the discrepancy is about $`2\%`$; in the worst possible case when $`m_{in}/m_{out}1`$ the discrepancy is still not large: $`100\%\sqrt{0.844}8\%`$.

### B Electron-phonon correction to the polaron energy and the effective mass

Summation over the wave vector $`\stackrel{}{k}`$ in Eq. (20) can be reduced to integration in a conventional way

$`{\displaystyle \underset{\stackrel{}{k}}{}}|V_\stackrel{}{k}|^2(\mathrm{})={\displaystyle \frac{V}{(2\pi )^3}}{\displaystyle 𝑑\stackrel{}{k}|V_\stackrel{}{k}|^2(\mathrm{})}=`$ (60)
$`(\mathrm{}\omega _{\text{LO}})^2\sqrt{{\displaystyle \frac{\mathrm{}}{2m\omega _{\text{LO}}}}}{\displaystyle \frac{\alpha }{2\pi ^2}}{\displaystyle \frac{d\stackrel{}{k}_{}dk_z}{k_{}^2+k_z^2}(\mathrm{})}.`$ (61)
Then, the integration over $`\stackrel{}{k}_{}`$ in Eq. (20) can be performed explicitly. As we are interested in corrections to the ground-state energy and the effective mass $`m_{eff}m+\mathrm{\Delta }_2m`$ of the polaron motion in the $`xy`$ plane, we expand $`\mathrm{\Delta }_2E(\stackrel{}{p}_{})\mathrm{\Delta }_2E{\displaystyle \frac{\mathrm{\Delta }_2m}{2m^2}}\stackrel{}{p}_{}^2`$ in powers of the conserved momentum $`\stackrel{}{p}_{}`$. Doing this the use is made of the integral

$`{\displaystyle \frac{d^2\stackrel{}{k}_{}}{(\stackrel{}{k}_{}^2+k_z^2)[\stackrel{}{k}_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}}2\stackrel{}{k}_{}\stackrel{}{p}_{}/\mathrm{}+b^2]}}\pi {\displaystyle \frac{\mathrm{ln}(k_z^2/b^2)}{k_z^2b^2}}+`$ (62)
$`\left({\displaystyle \frac{\stackrel{}{p}_{}}{\mathrm{}}}\right)^2\pi {\displaystyle \frac{k_z^4b^42k_z^2b^2\mathrm{ln}(k_z^2/b^2)}{b^2(k_z^2b^2)^3}}.`$ (63)
As the next step we use dimensionless “polaronic” units performing the scaling $`k_zk_z\sqrt{2m\omega _{\text{LO}}/\mathrm{}},zz\sqrt{\mathrm{}/2m\omega _{\text{LO}}}`$ and using the notation

$`l=L\sqrt{{\displaystyle \frac{2m\omega _{\text{LO}}}{\mathrm{}}}},\epsilon _N={\displaystyle \frac{E_N}{\mathrm{}\omega _{\text{LO}}}}.`$ (64)
In these units the correction to the ground-state energy takes the form

$`{\displaystyle \frac{\mathrm{\Delta }_2E}{\mathrm{}\omega _{\text{LO}}}}=`$ (65)
$`{\displaystyle \frac{\alpha }{\pi }}{\displaystyle \underset{N}{}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑k_z{\displaystyle \frac{\mathrm{ln}(k_z^2/b_N^2)}{k_z^2b_N^2}}\left(|G_s(N,k_z)|^2+|G_a(N,k_z)|^2\right),`$ (66)
(67)
where

$`b_N=\sqrt{\epsilon _N+1\epsilon _1}.`$ (68)
The correction to the effective mass reads as follows:

$`{\displaystyle \frac{\mathrm{\Delta }_2m}{m}}=`$ (69)
$`{\displaystyle \frac{\alpha }{\pi }}{\displaystyle \underset{N}{}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_z{\displaystyle \frac{k_z^4b_N^42k_z^2b_N^2\mathrm{ln}(k_z^2/b_N^2)}{b_N^2(k_z^2b_N^2)^3}}\times `$ (70)
$`\left(|G_s(N,k_z)|^2+|G_a(N,k_z)|^2\right).`$ (71)
Quantities $`G_j(N,k_z)`$ in Eqs. (67), (71) are given in dimensionless units by the same Eq. (21); the indices $`(a)s`$ are related to (anti)symmetrical wave functions being used in Eq. (21):

$`G_s(N,k_z)=2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑z\psi _{s,N}(z)\psi _{s,1}(z)\mathrm{cos}k_zz,`$ (72)
$`G_a(N,k_z)=2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑z\psi _{a,N}(z)\psi _{s,1}(z)\mathrm{sin}k_zz,`$ (73)
Evidently, the replacement $`Ll`$ should be done in the definition of the wave functions and their normalization constants; in addition $`L_zl_z`$ in Eqs. (43), (51) as well as in Eq. (38) for the energy levels of the discrete spectrum. Eq. (31) now reads as follows:

$`p_n=\sqrt{v_0q_n^2},v_0={\displaystyle \frac{V_0}{\mathrm{}\omega _{\text{LO}}}},q_{max}^2=v_0.`$ (74)
Eq. (37) takes the form

$`n_{max}=1+\left[{\displaystyle \frac{\sqrt{v_0}l}{\pi }}\right].`$ (75)
The relation of dimensionless energies of the discrete and continuous spectra with subsequent wave vectors takes the form $`\epsilon _n=q_n^2,\epsilon _q=q^2`$. All the changes mentioned should also be done in Eq. (55).

The final note of this section concerns summation over $`N`$ in Eqs. (67), (71):

$`{\displaystyle \underset{N}{}}(\mathrm{})={\displaystyle \underset{n=1}{\overset{n_{max}}{}}}(\mathrm{})+\underset{l_z\mathrm{}}{lim}{\displaystyle \frac{l_z}{2\pi }}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑p(\mathrm{}).`$ (76)
The replacing of the sum over the continuous spectrum by the integration over the wave vector $`p`$ follows from Eqs. (43, 51) in the limit $`L_z\mathrm{}`$. The wave vectors $`q`$ and $`p`$ are related to each other because of Eq. (44) which now takes the form $`q=\sqrt{p^2+v_0}`$. Note also that only $`G_s(N,k_z)`$ ($`G_a(N,k_z)`$) has to be taken into account for odd (even) $`n`$ in the sum over the discrete quantum number $`n`$.

The numerical results obtained are plotted in Fig. 1 for $`\mathrm{\Delta }_2E`$ and in Fig. 2a for $`\mathrm{\Delta }_2m/m`$. Because the mean effective mass $`m`$ depends on the potential width we also plotted in Fig. 2b the ratio of the mass shift to those in the well material, that is, the ratio

$`\delta _2m={\displaystyle \frac{\mathrm{\Delta }_2m}{\mathrm{\Delta }_2m_{in}}},\mathrm{\Delta }_2m_{in}=m_{in}{\displaystyle \frac{\alpha _{in}}{6}}.`$ (77)
The discussion of the numerical results is given in the last section.

## IV Rosen-Morse potential

### A Energy-dependent Green function

In this section we present another example — a multilayered heterostructure described by a confining potential $`V(z)`$ which is chosen in the form of the Rosen-Morse potential

$`V(z)`$ $`=`$ $`V_0\mathrm{tanh}^2\left({\displaystyle \frac{z}{L_{RM}}}\right),`$ (78)
$`V_0`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2mL_{RM}^2}}\kappa (\kappa +1).`$ (79)
where $`L_{RM}`$ is the parameter close to the half-width of the Rosen-Morse quantum well and $`\kappa `$ is the dimensionless parameter to govern the strength of the potential.

The summation (20) over the quantum number $`N`$ can be represented through the Green function which is known analytically for the Rosen-Morse potential. Namely, the second-order correction to the ground-state energy can be written in the form

$`\mathrm{\Delta }_2E`$ $`=`$ $`\mathrm{}\omega _{\text{LO}}\alpha {\displaystyle \frac{l_{RM}}{\sqrt{2}}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑k_{}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑z_a{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑z_be^{k_{}|z_az_b|}`$ (81)
$`\psi _1^{}(z_a)\psi _1(z_b)G(z_a,z_b;E),`$
where we made a scaling $`zzL_{RM},\stackrel{}{k}\stackrel{}{k}/L_{RM}`$ to use dimensionless variables $`z,\stackrel{}{k}`$ and integrated over $`k_z`$ and angles of $`\stackrel{}{k}_{}`$. The dimensionless parameter

$`l_{RM}=L_{RM}\sqrt{{\displaystyle \frac{m\omega _{\text{LO}}}{\mathrm{}}}}`$ (82)
is the width of the confining potential in polaronic units while the potential strength can now be written as follows:

$`V_0=\mathrm{}\omega _{\text{LO}}{\displaystyle \frac{\kappa (\kappa +1)}{2l_{RM}^2}}.`$ (83)
The quantity $`G(z_a,z_b;E)`$ is the Green function of the dimensionless Hamiltonian (2) which takes the form

$`H_{el,}^{\prime \prime }={\displaystyle \frac{1}{2}}{\displaystyle \frac{d^2}{dz^2}}+{\displaystyle \frac{\kappa (\kappa +1)}{2}}\mathrm{tanh}^2z,`$ (84)
that is $`G(z_a,z_b;E)=z_a|(H_{el,}^{\prime \prime }E)^1|z_b`$, while $`\psi _1(z)`$ is the ground-state wave function of the potential (84)

$`\psi _1(z)=\left[{\displaystyle \frac{\mathrm{\Gamma }(\kappa +1/2)}{\sqrt{\pi }\mathrm{\Gamma }(\kappa )}}\right]^{1/2}{\displaystyle \frac{1}{\mathrm{cosh}^\kappa z}}.`$ (85)
The ground-state energy of the Hamiltonian (84) is given by

$`E_1={\displaystyle \frac{\kappa }{2}}.`$ (86)
The energy $`E`$ in Eq. (81) reads as follows

$`E={\displaystyle \frac{k_{}^2}{2}}l_{RM}^2+{\displaystyle \frac{\kappa }{2}}.`$ (87)
The energy-dependent Green function of the system can be represented in the form:

$`G(z_a,z_b;E)=`$ (88)
$`{\displaystyle \frac{\mathrm{\Gamma }(\nu +\kappa +1)\mathrm{\Gamma }(\nu \kappa )}{\mathrm{\Gamma }^2(\nu +1)}}{\displaystyle \frac{1}{(4\mathrm{cosh}z_a\mathrm{cosh}z_b)^\nu }}\times `$ (89)
$`{}_{2}{}^{}F_{1}^{}(\nu \kappa ,\nu +\kappa +1;\nu +1;{\displaystyle \frac{1\mathrm{tanh}z_>}{2}})\times `$ (90)
$`{}_{2}{}^{}F_{1}^{}(\nu \kappa ,\nu +\kappa +1;\nu +1;{\displaystyle \frac{1+\mathrm{tanh}z_<}{2}}),`$ (91)
where $`z_>(z_<)`$ denotes the maximum (minimum) of $`z_a`$ and $`z_b`$. The parameter $`\nu `$ is defined by the relation

$`\nu =\sqrt{2\left(E{\displaystyle \frac{\kappa (\kappa +1)}{2}}\right)}=\sqrt{k_{}^2+\kappa ^2+2l_{RM}^2}.`$ (92)
The polaron effective mass can be represented in a similar way

$`{\displaystyle \frac{\mathrm{\Delta }_2m}{m}}`$ $`=`$ $`\alpha {\displaystyle \frac{l^3}{2\sqrt{2}}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑k_{}k_{}^2{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑z_a{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑z_be^{k_{}|z_az_b|}`$ (94)
$`\psi _1^{}(z_a)\psi _1(z_b){\displaystyle \frac{^2}{E^2}}G(z_a,z_b;E).`$
To simplify numerical calculations we may replace the derivative with respect to $`E`$ by the derivative with respect to $`\nu `$

$`{\displaystyle \frac{^2}{E^2}}={\displaystyle \frac{1}{\nu ^2}}{\displaystyle \frac{^2}{\nu ^2}}{\displaystyle \frac{1}{\nu ^3}}{\displaystyle \frac{}{\nu }}`$ (95)
and perform once the integration by parts. As the result, we arrive at the following representation equivalent to Eq. (94):

$`{\displaystyle \frac{\mathrm{\Delta }_2m}{m}}=\alpha {\displaystyle \frac{l^3}{2\sqrt{2}}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑k_{}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑z_a{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑z_b(1k_{}|z_az_b|)`$ (96)
$`e^{k_{}|z_az_b|}\psi _1^{}(z_a)\psi _1(z_b)\left[{\displaystyle \frac{1}{\nu }}{\displaystyle \frac{}{\nu }}\right]G(z_a,z_b;E).`$ (97)
Note that $`m,\alpha ,\omega _{\text{LO}}`$ in all these formulae stand for the mean characteristics of the medium. The wave function in their definitions is given by Eq. (85). The numerical results are plotted in Fig. 3 and discussed in the last Section.

### B Effective width

If we decide to compare the results for the rectangular and Rosen-Morse potentials, we have to define a parameter which plays the role of the effective width of the Rosen-Morse potential. That is, this parameter (for which we use a notation $`L`$) should be close to $`2L_{RM}`$ of Eq. (79) being also related to the rectangular potential. We accept the following definition: let us call the effective width of the Rosen-Morse potential the width $`L`$ of the rectangular well of the same height with the same ground-state energy in the absence of the electron-phonon interaction (that is, at $`\alpha =0`$). The advantage of this definition is that while calculating the polaron binding energy for the Rosen-Morse and rectangular potentials we subtract the same quantity in both the cases and can compare only energy shifts due to the electron-phonon interaction.

The ground-state energy of a rectangular potential with the height $`V_0`$ and width $`L`$ is given by the relations

$`E_{RC}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2k^2}{2m}},`$ (98)
$`\text{tan}{\displaystyle \frac{kL}{2}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{V_0}{E_{RC}}}1},`$ (99)
while the RM ground-state energy looks like

$`E_{RM}={\displaystyle \frac{\mathrm{}^2}{mL_{RM}^2}}{\displaystyle \frac{\kappa }{2}}`$ (100)
and the height $`V_0`$ of the potential is given by Eq. (79). With the equality $`E_{RM}=E_{RC}`$ we arrive at the relation between the parameter $`L_{RM}`$ of the Rosen-Morse potential and its effective width defined as has been discussed:

$`{\displaystyle \frac{L}{L_0}}`$ $`=`$ $`2\sqrt{\lambda }\text{arctg}\sqrt{\lambda 1},`$ (101)
$`\lambda `$ $`=`$ $`\kappa +1={\displaystyle \frac{1}{2}}\left[1+\sqrt{1+\left(2L_{RM}/L_0\right)^2}\right].`$ (102)
Here we introduce the distance scale

$`L_0=\sqrt{{\displaystyle \frac{\mathrm{}^2}{2mV_0}}}.`$ (103)
The relation to the other dimensionless parameter $`l_{RM}`$ of Eq. (82) is given by

$`{\displaystyle \frac{L_{RM}}{L_0}}=l\sqrt{{\displaystyle \frac{2V_0}{\mathrm{}\omega _{\text{LO}}}}}.`$ (104)
At small $`L_{RM}L_0`$ we obtain $`L2L_{RM}`$ from Eq. (102), that is indeed the parameter $`L_{RM}`$ plays a role of the half-width of the Rosen-Morse potential in this case. When $`L_{RM}L_0`$, it follows from Eq. (102) that $`L\pi \sqrt{L_{RM}L_0}`$.

The effective width $`L`$ defined in this subsection allows us to apply the results for the rectangular potential to the Rosen-Morse quantum well. The example is given in Fig. 3 where we plotted also the energy and the mass shifts for the rectangular potential vs. the parameter $`L_{RM}`$ related to $`L`$ as is described.

## V Numerical results and discussion

To proceed to the numerical calculations we need now the dependence of medium parameters on the $`AlAs`$ mole fraction $`x`$. At first we present the parametrization from the review by Adachi:

$`\alpha (z)=0.068+0.058x,`$ (106)
$`m(z)=m_e(0.0665+0.0835x),`$ (107)
$`\mathrm{}\omega (z)=(36.25+1.83x+17.12x^25.11x^3)\mathrm{meV},`$ (108)
which was used in numerical calculations by Hai, Peeters and Devreese. Here $`m_e`$ is the electron mass in vacuum and $`m(z)`$ — its band mass in the subsequent layer; the values of the electron-phonon coupling constant $`\alpha (z)`$ and the LO-phonon frequency $`\omega (z)`$ are also related to this layer.

Some comments are to the point. The expression for the electron band mass is nothing else but the linear interpolation between the values $`m=0.0665m_e`$ for $`GaAs`$ and $`m=0.150m_e`$ for $`AlAs`$. As to the LO-phonon frequency there are two phonon modes with different frequencies $`\omega ^{(G)}(z)`$ and $`\omega ^{(A)}(z)`$ for the $`GaAs`$-like and $`AlAs`$-like modes in $`Al_xGa_{1x}As`$ crystal. Experimental results of Ref.  are interpolated by the following formulae:

$`\mathrm{}\omega ^{(G)}(z)=(36.256.55x+1.79x^2)\mathrm{meV},`$ (110)
$`\mathrm{}\omega ^{(A)}(z)`$ $`=`$ $`(44.63+8.78x3.32x^2)\mathrm{meV}.`$ (111)
Because the exact theory of the two-phonon interaction in alloys where there are two-mode phonons present has not been reported, Adachi suggested to use the effective phonon frequency $`\omega =(1x)\omega ^{(G)}+x\omega ^{(A)}`$, that is the linear interpolation between these two modes. Inserting here the expressions (V) one arrives at the result (108).

As to the interpolation formula (106) for the Fröhlich coupling constant $`\alpha `$, the situation seems to be a bit inconsistent. Indeed, $`\alpha `$ depends on the values of the static $`\epsilon _0`$ and the high-frequency $`\epsilon _{\mathrm{}}`$ dielectric constants:

$`\alpha `$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}\omega }}{\displaystyle \frac{\overline{e}^2}{\sqrt{2}}}\sqrt{{\displaystyle \frac{m\omega }{\mathrm{}}}}\left({\displaystyle \frac{1}{\epsilon _{\mathrm{}}}}{\displaystyle \frac{1}{\epsilon _0}}\right)`$ (112)
$`=`$ $`116.643\left({\displaystyle \frac{1}{\epsilon _{\mathrm{}}}}{\displaystyle \frac{1}{\epsilon _0}}\right)\sqrt{{\displaystyle \frac{m}{m_e}}}\sqrt{{\displaystyle \frac{1\mathrm{meV}}{\mathrm{}\omega }}}.`$ (113)
Earlier measurements of $`\epsilon _0`$ of $`GaAs`$ have yielded widely different values ranging from 9.8 to 13.3 (see Ref.  and references therein). For instance, Kartheuser reports the values $`\epsilon _{\mathrm{}}=10.9`$ and $`\epsilon _0=12.83`$ and $`\mathrm{}\omega =36.75\mathrm{meV}`$ for $`GaAs`$. This leads to the result $`\alpha =0.068`$, which is widely known and used by many people.

On the other hand, Adachi used the more recent results for $`GaAs`$: $`\epsilon _0=13.18\pm 0.40`$ and $`\epsilon _{\mathrm{}}=10.89`$, and for $`AlAs`$: $`\epsilon _0=10.06\pm 0.04`$ and $`\epsilon _{\mathrm{}}=8.16\pm 0.02`$. This gives birth to his interpolation formulae:

$`\epsilon _0=13.183.12x,`$ (115)
$`\epsilon _{\mathrm{}}=10.892.73x,`$ (116)
Inserting formulae (107), (108) and (V) into Eq. (113) Adachi declared the result $`\alpha =0.126`$ for $`AlAs`$. Together with the value $`\alpha =0.068`$ reported in Ref.  this leads to the interpolation formulae (106). The problem is that both these values for $`\alpha `$ do not follow from the parametrizations mentioned above.

Taking the same values for $`AlAs`$ as Adachi did take ($`m=0.150m_e,\mathrm{}\omega =50.09\mathrm{meV},\epsilon _0=8.16,\epsilon _{\mathrm{}}=10.06`$) we arrive at the result $`\alpha =0.1477`$. Moreover, if one takes the same interpolation formulae (V) at $`x=0`$ one obtains the value $`\alpha =0.0797`$ for GaAs. That is, Adachi had to obtain the formula

$`\alpha (z)=0.0797+0.0680x`$ (117)
as a linear interpolation between the values of $`\alpha `$ in $`GaAs`$ and $`AlAs`$. Note, that this formulae can be presented in the form $`\alpha (z)=1.172(0.068+0.058x)`$. The expression between the brackets coincide (probably occasionally) with the Adachi interpolation formulae for $`\alpha `$ (cf. Eq. (106)). That is, the discrepancy of (106) and of our interpolation (117) is about 17% and do not depend on $`x`$. To be consistent we have to accept the parametrization (117) in what follows.

For the confining potential we take the expression derived from the band-gap energy fit of Ref.  and used in Ref. :

$`V(z)=600(1.155x+0.37x^2)\mathrm{meV}.`$ (118)
Thus, we use the parametrization (107), (108), (115), (116), (117) and the potential (118) in our numerical calculations.

The results of our study for the rectangular potential (which is formed by a layer of $`GaAs/Al_xGa_{1x}As`$) are shown in Fig. 1 for the polaronic energy shift and in Fig. 2 for the polaron effective mass at the $`AlAs`$ mole fraction $`x=0.3`$. The contribution of the discrete and continuous spectra are plotted separately for this potential. In Fig. 2a the relative mass shift $`\mathrm{\Delta }_2m/m`$ is shown where the mean mass $`m`$ also depends on the potential width $`L`$. Thus, the ratio $`\delta _2m=\mathrm{\Delta }_2m/\mathrm{\Delta }_2m_{in}`$ of the mass shifts in the potential and in $`GaAs`$ is presented also in Fig. 2b for the same $`AlAs`$ mole fraction. Evidently, the asymptotics of this curve is equal to the unity at large $`L`$ and to the ratio $`m_{out}\alpha _{out}/m_{in}\alpha _{in}`$ at $`L0`$.

We may conclude that the continuous spectrum dominates at small potential widths. At large widths its contribution could also be significant although it is smaller than the contribution of the discrete spectrum (especially in deep potential wells). We also confirm the conclusion of the preceding papers that the leading term approximation is not adequate to describe this system and leads to wrong asymptotics at both small and large potential widths (see the dashed lines in Figs 1, 2).

An example of a multilayered heterostructure is presented. The results for the energy and the effective mass for the polaron in the Rosen-Morse potential well are shown in Fig. 3. For the numerical calculations we fix the value $`V_0=227.9`$ meV in Eqs. (79), (83) which corresponds to the limiting mole fraction at large distances $`x_{\mathrm{}}=lim_z\mathrm{}x(z)=0.3`$. Thus, we obtain the dependence of the mole fraction $`x`$ on the coordinate $`z`$:

$`600(1.155x+0.37x^2)=227.9\mathrm{tanh}^2z.`$ (119)
Now Eqs. (107), (108), (117) allow one to define the dependence of parameters on the coordinate $`z`$ and to calculate the mean characteristics of the heterostructure.

The calculations were completely different in comparison with the rectangular potential: instead of the direct summation over all intermediate states we used the analytical expression for the Green function of the Rosen-Morse potential. The results obtained demonstrate a similar behavior which is also close numerically to the results for the rectangular potential. The polaronic energy and mass shifts for the rectangular quantum well are also plotted here (dashed line) vs. the Rosen-Morse width $`L_{RM}`$ obtained from $`L`$ as is described above. We see that both the energies almost coincide, which gives the opportunity to approximate different quantum wells by the rectangular potential. The discrepancy in the effective mass is larger but not so crucial. This serves also as an additional internal criterion of the validity of our calculations.

Thus, we obtained a monotonous behavior of $`\mathrm{\Delta }_2E`$ between the correct 3D limiting values $`\alpha _{in}\mathrm{}\omega _{in}`$ and $`\alpha _{out}\mathrm{}\omega _{out}`$ both for the rectangular and the Rosen-Morse potentials (see Fig. 1 and Fig. 3a). Actually the peaks are “hidden” and they reveal themselves if we plot the dimensionless energy shift $`\mathrm{\Delta }_2E/(\mathrm{}\omega _{\text{LO}}\alpha )`$ which has the same 3D limit (the unity) at both small and large potential widths. But in the “real” units (meV) the peaks are smoothed.

To compare our results with the calculations performed for the one-layer heterostructure we refer to the papers where the authors took into account the contributions of different phonon modes as well as mass and dielectric constant mismatches in the materials of the barrier and the well. Note that the analytical formulas of Ref.  contain a mistake — the wrong expression for the density of states. Namely, in some parts of the continuous spectrum contribution the integration is performed not over the wave vector $`p`$ but over the wave vector $`q`$ (that is $`_{V_0}^{\mathrm{}}𝑑E_z/\sqrt{E_z}(\mathrm{})`$ in the notations of that paper instead of the correct integration $`_{V_0}^{\mathrm{}}𝑑E_z/\sqrt{E_zV_0}(\mathrm{})`$). It is clear that this mistake results in lowering of the resulting curve for the energy, and the discrepancy is larger when the energy is closer to the potential edge, that is at small widths. This is just what we see in Fig. 4a comparing the result of Ref.  (the curve $`\mathrm{\Delta }_{HPD}E`$) with the new calculations of the same authors (the curve $`\mathrm{\Delta }_{HPD}E`$) which came to our knowledge when the present paper was already submitted for the publication.

Thus, our model does not reproduce the more complicated structure with the peak and the dip which was obtained in Ref. . Some hints on the existence of peaks can also be seen in our plots but the maximal values are so close to the asymptotics that the peaks are almost invisible. Probably, the dip appears because of the presence of several phonon modes (bulk, interface, etc.). At widths $`L50\AA `$ our results for the energy practically coincide with those of Ref. . The discrepancy at smaller widths seems to be more crucial. But the difference between the values in the peak and in the dip for the curve $`\mathrm{\Delta }_HE`$ in Fig. 4a is about 0.1 eV (3%). This phenomena hardly can be seen experimentally and this discrepancy is in the limits of the accuracy of our model estimated above. This gives indeed a strong support to our model and we may conclude that the latter provides us with the rather accurate approximation and can be used for more complicated calculations in multilayered heterostructures.

As to the shift of the electron band mass we found clear peaks for both the rectangular and the Rosen-Morse potentials (see Figs. 2, 3). As is seen in Fig. 2 the effective mass shift for the polaron in the rectangular quantum well has a peak at $`L20\AA (x=0.3)`$. Calculations show also that the larger is $`x`$ the smaller is the potential width corresponding to the peak. For the Rosen-Morse potential at $`x_{max}=0.3`$ the peaks in the effective mass occur at $`2L_{RM}20\AA `$. Note, that again the authors of Ref.  obtained curves with peaks and dips in contrast with our results (see Fig. 4b). The maximal discrepancy for the mass is about 11% at $`L3\AA `$ which is beyond the region available for experiments. Our results are very close to those of Ref.  at $`L10\AA `$ and practically coincde with them at $`L20\AA `$.

To compare our results with those of Ref.  we need now another parametrization used by these authors (although they refer also to the paper by Adachi). Namely, they took a slightly different expression for the confining potential:

$`V(z)=600(1.266x+0.26x^2)\mathrm{meV},`$ (120)
which follows from the band gap of Ref. . Furthermore, instead of the effective LO-phonon frequency they used the expression (110) for the energy of the $`GaAs`$-like phonons. The Fröhlich coupling constant $`\alpha `$ was calculated then using also the parametrization (107) and (V). Note, that these numerical calculation, as we found, can be approximated by the interpolation formula

$`\alpha (z)=0.0797+0.0772x+0.0295x^2.`$ (121)
The results of the comparison are shown in Fig. 5 (we used in our calculations for this plot the same parametrization as was used in Ref. ).

The curve $`\mathrm{\Delta }_{Ch}E`$ in Fig. 5a for the energy shift taken from Ref.  has also a small dip (qualitatively similar to this of Ref. ). But the discrepancy between energy shifts is much more drastic in this case, and we have no explanation for this.

It is clear that at large potential width only a bulk phonon mode inside the quantum well contributes so these curves should have the same limiting value $`\alpha _{in}\omega _{in}`$. Numerically we found $`\alpha _{in}=0.0797`$ and $`\omega _{in}=36.25`$ meV, so $`\alpha _{in}\omega _{in}=2.89`$ meV. Moreover, the behavior of the curves at large $`L`$ should be qualitatively and quantitatively the same which was the case when we compared our model with Refs. . In contrast with our model and the cited results by Hai, Peeters and Devreese the curve $`\mathrm{\Delta }_{Ch}E`$ in Fig. 5b approaches the asymptotics from below and the subsequent mechanism remains unclear. On the other hand there are some reasons why the curve have to approach its asymptotics from above. Indeed, at large potential width the particle does not feel yet the finite height of the potential, and the energy shift takes the same value as in the infinitely high potential which is a bit larger than the free polaron energy.

As to the opposite limit of the small width of the confining potential, it is surprising that the asymptotic value is not reached even at $`L0.3\AA ,`$ as is found in Ref. . Numerically we obtained $`\alpha _{out}=0.1014`$ in this scheme and $`\omega _{out}=34.72`$ meV, so $`\alpha _{out}\omega _{out}=3.52`$ meV.

Both asymptotic values coincide with what was obtained by the authors of Ref. . Looking at the behavior of the mass shift, we see that both curves coincide at large widths as it should be. At widths smaller than 100 Å the discrepancy becomes evident. But we may conclude that something is wrong with the numerical job of Ref.  because their curve approaches the wrong limit at $`L0`$. Indeed, in this limit the asymptotical value of the plotted ratio should be equal to $`\alpha _{out}m_{out}/\alpha _{in}m_{in}`$. As it follows from our analysis of the energy shift, we obtained the same values for the Fröhlich coupling constants. The values for the band masses follow from Eq. (107): $`m_{in}=0.0665m_e`$ and $`m_{out}=0.0874m_e`$ at $`x=0.25`$. Then, the asymptotical value of the plotted ratio should be equal to $`1.67`$, instead of 1.83 what was found in Ref.  That is, the discrepancy is about 10% in this limit and we cannot explain its origin as well.

It would be highly desirable to include a comparison of our results and the results by other authors with corresponding experiments. To the best of our knowledge, no such experiments do exist at the moment.

## VI Conclusions

To conclude, we suggested an approximate model to describe a multilayered $`GaAs/Al_xGa_{1x}As`$ heterostructure as an effective medium with one (bulk) phonon mode. The fundamental entity is the confining potential generated by these layers which we take into account explicitly. Then we calculate the mean characteristics of the electron in the effective medium (such as its band mass, phonon frequencies etc.) which depend on the form of the confining potential. With these parameters we calculated the energy and the effective mass of a polaron confined to a quasi-2D quantum well $`GaAs/Al_xGa_{1x}As`$ for different $`AlAs`$ mole fractions. The calculations include the full energy spectrum as intermediate states. Peaks are found for the effective mass at some potential widths while the energy demonstrates rather monotonous behavior between the correct 3D-limits. Finally, some discrepancies in the interpolation formulae for the experimental results are discussed as well as discrepancies with the previously obtained theoretical results. We demonstrated that our model gives practically the same (or very close) results as the explicit calculations of Ref.  for potential widths $`L10\AA `$.

###### Acknowledgements.

We thank J. T. Devreese, V. N. Gladilin, G. Q. Hai, H. Leschke, V. M. Fomin, F. M. Peeters, E. P. Pokatilov, and J. Wüsthoff for useful discussions and valuable remarks and advices. Special thanks are to the authors of the paper  for making their results available to us prior to publication and to the authors of the paper  who kindly provided us with their data-files. M.O.D. and M.A.S. are grateful to Dortmund University for the kind hospitality during their visits to Germany. Financial support of the Heisenberg-Landau program (Germany-JINR collaboration in theoretical physics) and Deutsche Forschungsgemeinschaft (Graduiertenkolleg GKP 50/2) is gratefully acknowledged.


format: text

text:
# Schematic phase diagram and collective excitations in the collisional regime for trapped boson-fermion mixtures at zero temperature

## 1 Introduction

The achievement of Bose-Einstein condensation in trapped atomic gases \[1 - 3\] has given new impulse to the study of dilute quantal fluids. Experimental attention has extended to the properties of double Bose condensates  and more recently to the attainment of quantal degeneracy in gases of fermionic atoms \[6 - 10\] and in dilute mixtures of bosons and fermions . A specific motivation for interest in boson-fermion mixtures is that the collisions between fermions and bosons can foster thermalization of the fermionic component and induce so-called sympathetic cooling . Elastic collisions between fermions in a single hyperfine state are inhibited by the Pauli principle  and even in multi-component Fermi gases the collisions between atoms in different hyperfine states become ineffective at very low temperature because of Fermi factors .

A number of theoretical studies have already been addressed to properties of dilute boson-fermion mixtures confined in harmonic traps. They have concerned the equilibrium density profiles of the two species both at zero and at finite temperature \[14 - 17\], the kinetic energy of the fermionic component  and collective excitations from sum rules in the collisionless regime . Our main focus here is on small-amplitude collective excitations for a mixture under harmonic confinement in the collisional regime. For a detailed discussion on how this regime may be achieved in dilute Fermi gases we refer to the work of Amoruso et al. .

We first set up a general background for our discussion of collective excitations in boson-fermion mixtures by constructing a schematic phase diagram from the earlier studies of the equilibrium state \[14 - 18\]. This helps to identify various dynamical behaviours, which may later be examined in detail by explicit numerical solution of the equations of motion. We then report analytic results for the eigenfrequencies of the mixture in limiting situations, corresponding to mixing of Bogolubov sound and fermionic first sound in a nearly homogeneous mixture and to surface modes in a regime of weak fermion-boson coupling.

We assume throughout a positive boson-boson scattering length, as is applicable to Bose condensates of 87Rb and 23Na. Mixtures of 6Li and 7Li should have special interest, once condensation is realized for 7Li in a hyperfine state where the scattering length is positive. According to calculations on the collisional properties of ultracold potassium , the 40K-41K mixture should also be an interesting system for study from the present viewpoint, if condensation of the rare isotope 41K can be achieved.

## 2 Equilibrium density profiles and phase diagram

We consider a system of $`N_f`$ fermions of mass $`m_f`$ and $`N_b`$ bosons of mass $`m_b`$ at $`T=0`$, confined by spherically symmetric external potentials $`V_{f,b}(r)=m_{f,b}\omega _{f,b}^2r^2/2`$ with frequencies $`\omega _f`$ and $`\omega _b`$. We assume that a single spin state is trapped for each component of the mixture and that all bosons are in the condensate. We can then omit the fermion-fermion interaction, which is inhibited by the Pauli exclusion principle. The boson-boson and boson-fermion interactions are described by contact potentials with scattering lengths $`a_{bb}>0`$ and $`a_{bf}`$, yielding coupling-strength parameters $`g=4\pi \mathrm{}^2a_{bb}/m_b>0`$ and $`f=2\pi \mathrm{}^2a_{bf}/m_r`$ where $`m_r`$ is the reduced boson-fermion mass.

The equations of motion for the partial particle densities $`n_\sigma (𝐫,t)`$ and partial current densities $`𝐣_\sigma (𝐫,t)`$, where $`\sigma `$ is an index for the atomic species ($`\sigma =f`$ or $`b`$) follow from those for the one-body density matrix upon projection on the main diagonal (see for instance ). They take the form of continuity equations,

$$_tn_\sigma (𝐫,t)=𝐣_\sigma (𝐫,t)$$
(2.1)
and of generalized hydrodynamic equations which, in a mean-field approximation and neglecting coupling to energy fluctuations, can be written as

$$m_f_t𝐣_f(𝐫,t)=\mathrm{\Pi }^{(f)}(𝐫,t)n_f(𝐫,t)\left[V_f(𝐫)+fn_b(𝐫,t)\mu _f\right]$$
(2.2)
and

$`m_b_t𝐣_b(𝐫,t)`$ $`=`$ $`\mathrm{\Pi }^{(b)}(𝐫,t)`$ (2.3)
$`n_b(𝐫,t)\left[V_b(𝐫)+gn_b(𝐫,t)+fn_f(𝐫,t)\mu _b\right]`$
In these equations $`\mathrm{\Pi }^{(\sigma )}(𝐫,t)`$ are the kinetic stress tensors and $`\mu _\sigma `$ are the chemical potentials of the two atomic species. In the following we treat the kinetic stress tensors by a local density approximation for fermions,

$$\mathrm{\Pi }_{ij}^{(f)}(𝐫,t)=\frac{2}{5}A[n_f(𝐫,t)]^{5/3}\delta _{ij}$$
(2.4)
with $`A=\mathrm{}^2(6\pi ^2)^{2/3}/2m_f`$ and by a Thomas-Fermi approximation for bosons,

$$\mathrm{\Pi }^{(b)}(𝐫,t)=0.$$
(2.5)
In these equations we have dropped terms which are quadratic in the velocity field.

The treatment given above assumes that local equilibrium has been established for the ground-state density profiles and that it is maintained during dynamic fluctuations of the particle densities. The dynamical behaviour that we shall discuss in Sect.3 will therefore correspond to a collisional regime .

### 2.1  Equilibrium density profiles

The equilibrium density profiles ($`n_\sigma ^{eq}(r)`$, say) follow from Eqs. (2.2) - (2.5) by setting the current densities to zero. The results are well known :

$$n_b^{eq}(r)=g^1\left[\mu _bV_b(r)fn_f^{eq}(r)\right]\theta (\mu _bV_b(r)fn_f^{eq}(r))$$
(2.6)
and

$$n_f^{eq}(r)=A^{3/2}\left[\mu _fV_f(r)fn_b^{eq}(r)\right]^{3/2}\theta (\mu _fV_f(r)fn_b^{eq}(r)).$$
(2.7)
The chemical potentials are to be determined by self-consistency conditions on the particles numbers, $`N_\sigma =n_\sigma ^{eq}(r)𝑑𝐫`$. On comparing the two density profiles in the limit of vanishing boson-fermion coupling, it is easily seen that the quantity $`(2A/3)[n_f^{eq}(r)]^{1/3}`$ may be viewed as an effective fermion-fermion coupling arising from the kinetic pressure of the Fermi gas . It is also useful to introduce the radii $`R_b=(2\mu _b/m_b\omega _b^2)^{1/2}`$ and $`R_f=(2\mu _f/m_f\omega _f^2)^{1/2}`$ for the two clouds in the absence of interactions, but taken at the true values of the chemical potentials.

For values of the boson-boson coupling and of the confinement frequencies which are typical of current experiments, and assuming comparable magnitudes for the boson-boson and the boson-fermion coupling, the fermionic cloud is considerably more dilute than the bosonic cloud. In this situation the boson cloud is essentially unaffected by the interactions with the fermionic component and the latter can be treated as confined inside a double-parabola effective potential $`V_f^{eff}(r)`$ . This is given by

$$V_f^{eff}(r)=\{\begin{array}{cc}m_f\omega _f^2(1\gamma )r^2/2+f\mu _b/g\hfill & \hfill \mathrm{for}r<R_b\\ m_f\omega _f^2r^2/2\hfill & \hfill \mathrm{for}r>R_b\end{array}$$
(2.8)
where

$$\gamma =(fm_b\omega _b^2)/(gm_f\omega _f^2).$$
(2.9)
Evidently, if $`\gamma >1`$ the effective potential (2.8) has negative concavity inside the Bose radius and the fermions will tend to be expelled from the centre of the trap. The result at strong coupling is phase separation in the mixture, with the fermions forming a shell around the boson cloud .

A second type of phase separation for the boson-fermion mixture was demonstrated in the work of Mølmer . In this case one takes $`N_fN_b`$ and considers very large values of the boson-boson coupling $`g`$, such that the two clouds acquire similar densities. One ultimately finds that the bosons are expelled from the centre of the trap and form a shell around the fermion cloud.

Finally, Eq. (2.8) is valid also in the case $`\gamma <0`$, i.e. for attractive boson-fermion interactions, as long as the fermionic cloud is more dilute than the boson cloud. However, the fermions now tend to draw the bosons towards the centre of the trap, ultimately leading to collapse once $`\gamma `$ becomes sufficiently negative .

### 2.2  Schematic phase diagram at zero temperature

As shown by Vichi et al. , the thermodynamics of the boson-fermion mixture at $`T=0`$ is well characterized by means of two scaling parameters, which are $`\gamma `$ in Eq. (2.9) and

$$x=\left(\frac{R_b}{R_f}\right)^{1/2}=\left(\frac{m_f\omega _f}{2m_b\omega _b}\right)^{1/2}\left[15\frac{a_{bb}}{a_{ho}}\frac{N_b}{(6N_f)^{5/6}}\right]^{1/5},$$
(2.10)
with $`a_{ho}=(\mathrm{}/m_b\omega _b)^{1/2}`$. The description of the system with only two scaling parameters, instead of the eight original ones entering the expression of the equilibrium density profiles, is a major simplification. The ratio $`R_b/R_f`$ determines the deviation of the kinetic energy of the fermion cloud from its ideal-gas value .

For the purposes of the present discussion we make the simplifying assumptions $`m_f=m_b`$, $`\omega _f=\omega _b`$ and $`N_f=N_b`$. In this case the two scaling parameters take the simple expressions $`x=(15a_{bb}k_f^{(0)}/48)^{1/5}`$ and $`\gamma =a_{bf}/a_{bb}`$, with $`k_f^{(0)}=(48N_f)^{1/6}/a_{ho}`$ being the Fermi wave number at the trap centre for a non-interacting Fermi gas. As remarked under Eq. (2.7) above, the quantity $`a_{ff}3\pi /4k_f^{(0)}`$ may be viewed as an effective fermion-fermion scattering length arising from the Pauli kinetic pressure. We can now proceed to construct a schematic phase diagram for the boson-fermion mixture in the plane defined by the variables $`\gamma =a_{bf}/a_{bb}`$ and $`y=a_{bb}k_f^{(0)}`$. Figure 1 shows the result for the case $`\gamma >0`$ (i.e. repulsive boson-fermion interactions).

We search for the boundary of phase separation by means of a simple condition of linear stability on the two-by-two matrix of scattering lengths $`a_{bb}`$, $`a_{bf}`$ and $`a_{ff}`$  (see also ). A mixed state is stable if the inequality $`a_{bb}a_{ff}a_{bf}^2>0`$ holds, i.e. when

$$k_f^{(0)}a_{bb}\frac{3\pi }{4}\left(\frac{a_{bb}}{a_{bf}}\right)^2$$
(2.11)
Within this region we can distinguish a region of appreciable overlap between the two clouds (region I in Figure 1), bounded by two regions of diminishing overlap: a region at $`\gamma >1`$ where the fermions are being expelled from the centre of the trap (region II in Figure 1) and by a region at $`y>4/3\pi `$ where the bosons are being expelled from the centre of the trap (region III in Figure 1). If the spherical symmetry of the trapping potential is preserved in the demixed state, these two regions end into spherical-shell structures with the fermions at large $`\gamma `$ (or the bosons at large $`y`$) forming the outer shell.

## 3  Small amplitude oscillations around the equilibrium profiles

The equations of motion for small-amplitude density fluctuations $`\stackrel{~}{n}_f(𝐫,t)`$ and $`\stackrel{~}{n}_b(𝐫,t)`$ in the collisional regime are found by combining Eq. (2.1) with the linearized form of Eqs. (2.2) and (2.3). They are

$$_t^2\stackrel{~}{n}_f(𝐫,t)=\frac{1}{m_f}\left\{n_f^{eq}(r)\left[\frac{2}{3}A(n_f^{eq}(r))^{1/3}\stackrel{~}{n}_f(𝐫,t)+f\stackrel{~}{n}_b(𝐫,t)\right]\right\}$$
(3.1)
and

$$_t^2\stackrel{~}{n}_b(𝐫,t)=\frac{1}{m_b}\left\{n_b^{eq}(r)\left[g\stackrel{~}{n}_b(𝐫,t)+f\stackrel{~}{n}_f(𝐫,t)\right]\right\}.$$
(3.2)
Again the quantity $`(2A/3)[n_f^{eq}(r)]^{1/3}`$ enters Eq. (3.1) as an effective fermion-fermion coupling.

It is easily checked that the form of Eqs. (3.1) and (3.2) is such as to satisfy the generalized Kohn theorem  for the centre-of-mass coordinate $`𝐱(t)`$ of the whole fluid when the two confinement frequencies coincide. In this mode of motion we have $`n_{b,f}(𝐫,t)=n_{b,f}^{eq}(𝐫𝐱(t))`$ and hence

$$\stackrel{~}{n}_{b,f}(𝐫,t)=𝐱(t)n_{b,f}^{eq}(r)$$
(3.3)
for small amplitude oscillations. By substituting Eq. (3.3) in Eqs. (3.1) and (3.2) and using Eqs. (2.6) and (2.7), we find

$$m_{b,f}_t^2𝐱(t)n_{b,f}^{eq}(r)=\underset{i,j}{}x_j(t)_i\left[n_{b,f}^{eq}(r)_i_jV_{b,f}(r)\right].$$
(3.4)
Therefore, the centre of mass of the system oscillates at the bare trap frequency.

Evidently, a general solution of Eqs. (3.1) and (3.2) can only be obtained numerically, since the equilibrium density profiles are not known analytically from Eqs. (2.6) and (2.7). The dynamics of the fluid as a function of the scaling parameters will be especially rich in the regions of the phase diagram in which the nature of the partial density profiles is changing on the approach to phase separation. Of course, this will be signalled by a softening of the modes associated with concentration fluctuations. Below we restrict ourselves to analytic results which are applicable to special regimes.

### 3.1  Homogeneous limit

The dynamics of a spatially homogeneous mixture is relevant to the situation in which the density profiles of both components are slowly varying in space, on the length scale set by the mean interparticle distances. Taking the equilibrium densities as constant, it is easily seen from Eqs. (3.1) and (3.2) that the eigenmodes are two sound waves, with linear dispersion relation $`\omega _{1,2}(q)=c_{1,2}q`$. The sound velocities are given by

$$c_{1,2}^2=\frac{1}{2}\left[c_b^2+c_f^2\pm \sqrt{(c_b^2c_f^2)^2+4f^2n_b^{eq}n_f^{eq}/m_bm_f}\right].$$
(3.5)
In Eq. (3.5) $`c_b=(gn_b^{eq}/m_b)^{1/2}`$ is the speed of Bogolubov sound and $`c_f=[2A(n_f^{eq})^{2/3}/3m_f]^{1/2}`$ is that of first sound in a Fermi gas.

This result may be useful in regard to measurements of the speed of sound in elongated systems, in experiments such as already carried out by Matthews et al.  on a single Bose condensate.

### 3.2  Surface modes for weak boson-fermion coupling

Assuming equal confinement frequencies for the two components of the mixture ($`\omega _{ho}`$, say), it is known from earlier results of Stringari  for a Bose condensate and of Amoruso et al.  for a one-component Fermi gas that in the limit $`f=0`$ the frequencies of the surface modes ($`n`$= 0) with angular momentum number are equal and given by $`\omega _l=l^{1/2}\omega _{ho}`$. A perturbative treatment of the effect of fermion-boson coupling is therefore especially simple for these modes.

We define the unperturbed equilibrium density profiles $`n_{b,f}^{(0)}(r)`$ and the $`𝐫`$-dependent amplitudes $`\stackrel{~}{n}_{b,f}^{(l,m)}(𝐫)`$ of the unperturbed density fluctuations, for each value of the angular momentum $`l`$ and of its $`z`$-component $`m`$. After Fourier-transforming Eqs. (3.1) and (3.2) and linearizing them in the coupling parameter $`f`$, we multiply Eq. (3.1) by $`[n_f^{(0)}(r)]^{1/3}\stackrel{~}{n}_f^{(l,m)}(𝐫)`$ and Eq. (3.2) by $`\stackrel{~}{n}_b^{(l,m)}(𝐫)`$, and then integrate both equations over $`𝐫`$. We obtain in this way a determinantal equation for the shift of the eigenmode frequencies $`\mathrm{\Omega }_l`$ due to the fermion-boson coupling. Writing

$$\mathrm{\Omega }_l^2=\omega _{ho}^2(l+\mathrm{\Delta }_l),$$
(3.6)
we find

$$\mathrm{\Delta }_l=\frac{a_{bf}}{2a_{ho}}\left[(E+B)\pm \mathrm{sign}(a_{bf})\sqrt{(EB)^2+4CD}\right],$$
(3.7)
this result being valid for either sign of $`a_{bf}`$. The quantities entering Eq. (6) are defined in the Appendix, where it is also shown that they can all be expressed in terms of the Gauss hypergeometric function .

The result given in Eq. (3.7) should be better applicable in the case $`R_b<R_f`$, since in deriving it we have assumed that the equilibrium density profiles for both components are perturbed in a symmetric manner. This corresponds to working in region I of the phase diagram (and in the corresponding region at negative $`\gamma `$).

## 4  Summary and discussion

The main features of the zero-temperature phase diagram that we have sketched in Figure 1 for boson-fermion mixtures arise from the competition between the kinetic energy of the Fermi gas and the repulsive boson-boson and boson-fermion interactions. The former disfavours localization and phase separation, while the latter favour the spatial separation of the two components. On varying these system parameters while preserving the spherical symmetry of the trapped fluid mixture, a spontaneous symmetry breaking occurs towards a situation where the two components are spatially separated along the radial direction. Evidently, axial separation of the two components would instead be observed in an anisotropic trap. One may think of exploiting Feshbach resonances to tune the scattering lengths towards such phase-separation regime.

Given the complex behaviour of the density profiles on the approach to phase separation, a complete investigation of the oscillatory eigenmodes of a boson-fermion mixture can only be carried out by numerical means. This will be well worth doing once a mixture of specific experimental relevance and with reasonably known scattering lengths is identified. Of course, the dynamics in the phase-separated region is simply related to that of the two pure components, aside for the presence of interfacial modes. On the other hand, we have seen that suggestive analytic results can be obtained in the quasi-homogeneous limit and in the weak boson-coupling regime, for either sign of the boson-fermion coupling. With increasingly large attractions between boson and fermions, however, the mixture is expected to undergo collapse .

We thank Dr Ilaria Meccoli for many useful discussions and for her help in the early stages of this work.

## Appendix A Calculation of frequency shifts for surface modes

We introduce the notations $`|F=\stackrel{~}{n}_f^{(l,m)}(𝐫)`$ and $`|B=\stackrel{~}{n}_b^{(l,m)}(𝐫)`$, define the operators $`P_{b,f}=(f/m_{b,f})[n_{b,f}^{(0)}(r)]`$ and $`D_f=(f/m_f)\{(n_f^{(0)}(r))^{1/3}[(n_f^{(0)}(r))^{1/3}]\}`$, and introduce the scalar-product notation $`F|\widehat{O}|F=\stackrel{~}{n}_f^{(l,m)}(𝐫)\widehat{O}n_f^{(l,m)}(𝐫)𝑑𝐫`$ etcetera. With these notations we find the following expressions for the quantities entering Eq. (3.7):

$`B`$ $`=`$ $`F|(n_f^{(0)}(r))^{1/3}D_f|F/F|(n_f^{(0)}(r))^{1/3}|F`$ (A.1)
$`=`$ $`{\displaystyle \frac{2\pi ^{1/2}l}{\alpha _{bb}X_f^{2l+2}}}{\displaystyle \frac{\mathrm{\Gamma }(l+3)}{\mathrm{\Gamma }(l+3/2)}}{\displaystyle _0^{X_b}}𝑑xx^{2l+2}(X_f^2x^2)^{1/2}(2X_f^2+X_b^23x^2),`$
$`C`$ $`=`$ $`F|(n_f^{(0)}(r))^{1/3}P_f|B/F|(n_f^{(0)}(r))^{1/3}|F`$ (A.2)
$`=`$ $`{\displaystyle \frac{2l}{\pi ^{5/2}X_f^{2l+2}}}{\displaystyle \frac{\mathrm{\Gamma }(l+3)}{\mathrm{\Gamma }(l+3/2)}}{\displaystyle _0^{X_b}}𝑑xx^{2l+2}(X_f^2x^2)^{1/2},`$
$`D`$ $`=`$ $`B|P_b|F/B|B`$ (A.3)
$`=`$ $`{\displaystyle \frac{(2l+3)}{2\alpha _{bb}X_b^{2l+3}}}{\displaystyle _0^{X_b}}dxx^{2l+2}\{2l(X_f^2x^2)^{1/2}+X_f^2(X_f^2x^2)^{3/2}(X_b^2x^2)`$
$`+(X_f^2x^2)^{1/2}[x^2+2(l+1)(X_b^2x^2)]\},`$
and

$`E`$ $`=`$ $`B|P_f|B/B|B`$ (A.4)
$`=`$ $`{\displaystyle \frac{l(2l+3)}{2\pi ^2X_b^{2l+3}}}{\displaystyle _0^{X_b}}𝑑xx^{2l+2}(X_f^2x^2)^{1/2}.`$
In these equations $`\mathrm{\Gamma }(n)`$ is the Gamma function and we have set $`X_{b,f}=R_{b,f}/a_{ho}`$ and $`\alpha _{bb}=a_{bb}/a_{ho}`$.

All the integrals in Eqs. (A.1) -(A.4) can be expressed through the Gauss hypergeometric function $`{}_{2}{}^{}F_{1}^{}(a,b;c;z)`$, using the relation

$${}_{2}{}^{}F_{1}^{}(\frac{1+\alpha }{2},\frac{n}{2};\frac{3+\alpha }{2};\frac{X_b^2}{X_f^2})=\frac{1+\alpha }{X_f^nX_b^{\alpha +1}}_0^{X_b}𝑑xx^\alpha (X_f^2x^2)^{n/2}.$$
(A.5)


text:
# 1 The law of conservation of 𝜎-number.

## 1  The law of conservation of $`\sigma `$-number.

It is known that upper and lower quarks are positioned asymmetrically on the ”charge axis”. With the aim to bring in some symmetry to this asymmetric disposition we will introduce the new additive quantum number $`\sigma `$ determined so as to result, in combination with quark electric charge $`q`$, the charge of the respective lepton.

Then for $`u`$ and $`d`$ quarks we will have, respectively, $`\sigma _u=1/3`$ and $`\sigma _d=2/3`$. In general, $`\sigma `$-numbers for all quarks and anti-quarks are determined by the following formula:

$$\sigma =q1/3,\stackrel{~}{\sigma }=\stackrel{~}{q}+1/3.$$
(1)
Now the quarks (and leptons) are positioned symmetrically with respect to $`\sigma =q`$ axis on $`(q,\sigma )`$ plane (see Fig.1). This brings about an idea that together with the law of conservation of electric charge there is realized the second law of conservation, the law of conservation of $`\sigma `$-number. As we will see later, this idea results in extremely valuable findings. In particular, the conservation of $`\sigma `$ -number results in conservation of the difference between baryon and lepton numbers obtained earlier by Georgi and Glashow from $`SU(5)`$ symmetry.

Let us first show that $`\sigma `$-numbers for baryons, mesons, leptons and photons are determined by the following expressions:

$$\sigma _B=Q_B1,\sigma _M=Q_M,\sigma _L=Q_L+1,\sigma _\gamma =0,\stackrel{~}{\sigma }=\sigma $$
(2)
where $`Q`$ is an electric charge of a particle. Indeed, formula (2) is obvious for baryons and mesons due to their quark structure. Similarly, $`\sigma `$-numbers can be obtained for exotic adrons. Using (2), we will obtain, in particular, $`\sigma (p)=0,\sigma (n^0)=1,\sigma (\pi ^0)=0,\sigma (\stackrel{~}{\mathrm{\Lambda }}^0)=1,\sigma (^{})=2`$ and so on. Taking into account that $`\sigma (\pi ^0)=0`$, we will obtain $`\sigma _\gamma =0`$from $`\pi ^02\gamma `$decay.

In order to find $`\sigma _L`$, we will give the same quantum number (by analogy with the electric charge) to all charged leptons $`(\sigma _e=\sigma _\mu =\sigma _\tau )`$and another quantum number to neutral leptons $`(\sigma _{\nu _e}=\sigma _{\nu _\mu }=\sigma _{\nu _\tau })`$ . These numbers are related to each other by the formula $`\sigma _{\nu _e}=1+\sigma _e`$that can be derived from $`n^0p+e^{}+\stackrel{~}{\nu }`$decay. Then from $`uuxe^+\stackrel{~}{d}`$reaction where the same boson may decay into antilepton + antiquark or quark pair  we derive $`\sigma _{e^+}=\sigma _e^{}=0`$ , and $`\sigma _\nu =1`$. As we can see, these results are in accordance with expressions (2). After we have derived formula (2), we can check out that the conservation of $`\sigma `$-number takes place for all the reactions observed so far . It does not prohibit also possible proton decay through channels $`pe^+\pi ^0;e^+\pi ^+\pi ^{}`$ where conservation of baryon and lepton numbers is violated.

The general formula for conservation of $`\sigma `$-number can be written as follows:

$$Q_BN_B+\stackrel{~}{Q}_B+\stackrel{~}{N}_B+Q_M+\stackrel{~}{Q}_M+Q_L+N_L+\stackrel{~}{Q}_L\stackrel{~}{N}_L=const$$
(3)
where $`N`$ is the number of particles. Here, taking into account the conservation of electric charge, we will finally obtain the desired result:

$$(N_B\stackrel{~}{N}_B)(N_L\stackrel{~}{N}_L)=const$$
(4)
It should be pointed out that the baryon number in formula (4) could be replaced by total number of quarks$`N_k`$ comprised of baryon and meson quark numbers. Indeed, taking into account that for baryons and mesons $`(N_k\stackrel{~}{N}_k)=3(N_B\stackrel{~}{N}_B)`$and $`(N_k\stackrel{~}{N}_k)=0`$and using formula (4), we will obtain the following relation between quarks and leptons

$$(N_k\stackrel{~}{N}_k)3(N_L\stackrel{~}{N}_L)=const$$
(5)
As we will see later, the law of conservation of $`\sigma `$-number predicts that electron type neutrino mass is exactly zero. It should be mentioned in this relation that two contradictory theories trying to resolve the problem of neutrino mass have been developed so far by Dirak-Weyl in 1929 and Majoran in1936. According to the former  electron neutrino mass is exactly zero and the conservation of lepton number takes place. According to the latter  the electron neutrino mass is not zero and therefore the conservation of lepton number does not take place. The Majoran theory allows neutrinoless double beta decay where two $`d`$-quarks decay through channel $`du+e^{}`$with violation of lepton number.

The contemporary theories basing on $`SU(5)`$ and $`SO(10)`$ symmetry result in the same contradictory predictions  ”The $`SO(10)`$ symmetry allows occurrence of some phenomena prohibited by $`SU(5)`$ symmetry. In particular, $`SU(5)`$ theory predicts conservation of difference $`BL`$ while conservation of baryon number $`B`$ and lepton number $`L`$does not take place. In $`SO(10)`$ theory, the difference $`BL`$ may not hold constant if sufficient number of Higgs fields are involved”. Approximately the same results are mentioned in the Gell-Mann report.

Thus, the main question relating to neutrino mass remains open. There have been a number of attempts for the past 60 years to find theoretical or experimental solution but the problem still exists. The recent experiments in many scientific centers worldwide have been aimed at observing of neutrinoless double beta decay and neutrino oscillations. Currently several new projects are proposed and huge scientific potential and funds mobilized to get a breakthrough on this fundamental problem. Twenty-three scientific centers are involved in the project MINOS (USA). The underground experiment is aimed at neutrino beams detecting after passing some 730km from Fermi-lab (Viskonsin) to Soudan (Minnesota). It is expected that some of the detected neutrinos will have changed their flavor due to neutrino oscillations. In Japanese project KEK neutrino beam will pass some 230km from an accelerator to the detector Super-Kamiokande.  The new Heidelberg project GENIUS is aimed at increasing the sensitivity of Majorana neutrino mass from the present level of at best 0.1ev down to 0.01 or even 0.001ev.

Summarizing, we can say that the newly proposed law of conservation of $`\sigma `$-number prohibits neutrinoless double beta decay and, in compliance with Dirac-Weyl and $`SU(5)`$ theories, predicts that electron type neutrino mass is exactly zero.

## 2 Quark-lepton symmetry.

Now we will focus on the problem of quark-lepton symmetry. The essence of the problem is to find such a symmetry among the elementary particles, which would explain why each quark-lepton generation is necessarily comprised of the pair of leptons and pair of quarks. Let us shift to the new system of coordinates $`(q^{^{}},\sigma ^{^{}})`$ with the point of origin in the center of a rectangle $`e\nu _eud`$ and axes directed along the axes of symmetry (see Fig1.). The relation with an original system of coordinates is given by the following expressions:

$$q^{^{}}=\frac{1}{\sqrt{2}}q+\frac{1}{\sqrt{2}}\sigma ,\sigma ^{^{}}=\frac{1}{\sqrt{2}}q+\frac{1}{\sqrt{2}}\sigma \frac{1}{3\sqrt{2}}$$
(6)
The coordinates of the vertices are now as follows:

$$e(\frac{1}{\sqrt{2}},\frac{\sqrt{2}}{3}),\nu _e(\frac{1}{\sqrt{2}},\frac{\sqrt{2}}{3}),u(\frac{1}{\sqrt{2}},\frac{\sqrt{2}}{3}),d(\frac{1}{\sqrt{2}},\frac{\sqrt{2}}{3});$$
(7)
The group that represents the symmetry transformation for the above mentioned rectangle, consists of the following elements: $`E`$-identity transformation; turn around axis $`z`$ for angle $`2\pi `$ ; $`A`$-turn around axis $`\sigma ^{^{}}`$ for angle $`\pi `$; $`B`$-turn around axis $`q^{^{}}`$ for angle $`\pi `$; $`C`$-turn around axis $`z`$ for angle $`\pi `$.

For these elements we will obtain the group multiplication table (see Table 1).

Under multiplication we understand the subsequent execution of the corresponding operations. The elements of the group have order 2 (except for identity element ), since $`\chi ^2=E`$and $`\chi ^1=\chi `$, where $`\chi `$is an arbitrary element of the group. Thus, the totality of the elements $`E,A,B,C`$ makes up the Abelian group. The ”self-transformation” of a regular polygon is expressed by means of the following matrices

$$D_k=\left(\begin{array}{cc}\mathrm{cos}\frac{2\pi k}{n},\hfill & \mathrm{sin}\frac{2\pi k}{n}\hfill \\ \mathrm{sin}\frac{2\pi k}{n},\hfill & \mathrm{cos}\frac{2\pi k}{n}\hfill \end{array}\right),U_k=\left(\begin{array}{cc}\mathrm{cos}\frac{2\pi k}{n},\hfill & \mathrm{sin}\frac{2\pi k}{n}\hfill \\ \mathrm{sin}\frac{2\pi k}{n},\hfill & \mathrm{cos}\frac{2\pi k}{n}\hfill \end{array}\right).$$
(8)
where $`k=0,1,2,,n1`$. These $`2n`$-dimensional matrices make up the group of order $`2n`$ known as the group of diedr. In case of n=2 we have the simplest case of the group with elements

$$E=\left(\begin{array}{cc}1\hfill & 0\hfill \\ 0\hfill & 1\hfill \end{array}\right),A=\left(\begin{array}{cc}1\hfill & 0\hfill \\ 0\hfill & 1\hfill \end{array}\right),B=\left(\begin{array}{cc}1\hfill & 0\hfill \\ 0\hfill & 1\hfill \end{array}\right),C=\left(\begin{array}{cc}1\hfill & 0\hfill \\ 0\hfill & 1\hfill \end{array}\right).$$
(9)
Taking into account that the group (9) is isomorphic to the above mentioned group, we come to conclusion that (9) is the matrix representation of the group. Thus, if particles of the same generation are located in the vertices of rectangle (7) on the plane $`(q^{^{}},\sigma ^{^{}})`$, this distribution is subject to the symmetry transformation described by the group (9). It is also easy to show that if distribution of particles on the plane $`(q^{^{}},\sigma ^{^{}})`$ is subject to the symmetry transformation described by the group (9) and if coordinates of any arbitrary particle from (7) coincide with one of the vertices, there should be three more particles (and only three, without taking into account the color of the quarks) whose coordinates coincide with the remaining vertices.

Summarizing this paragraph, we can say that the nature of the quark-lepton symmetry could be described by $`q\sigma `$-symmetry. It is reflected by the fact that particles of the same generation are subject to the symmetry transformation represented by 4-group of diedr.

## 3  Symmetry among colorless elementary particles.

It is natural to expect that there is some symmetry among colorless particles as well. We will consider all the baryons and mesons which can be comprised of eighteen known quarks as well as the leptons. Each point on the plane $`(q,\sigma )`$ is assumed to ”contain” the whole family of the particles (see fig.2). For example, point $`(1;0)`$ contains, besides electron, other charged leptons $`\mu ^{}`$and $`\tau ^{}`$. It is clear from the picture that there is symmetry with respect to $`\sigma ^{^{}}`$axis. The number of particles in symmetrical points is identical. It is easy to show that the same symmetry exists for multi-quark baryons as well. For example, if baryons with quantum numbers $`q=3,\sigma =2,B=5`$ consist of eight upper and seven lower quarks, then symmetrical baryons with quantum numbers $`q=2,\sigma =3,B=5`$ consist of seven upper and eight lower quarks. It is obvious that the number of colorless combinations is identical.

For baryons the same symmetry exists also with respect to other axis $`q^{^{}}`$ which is determined by equation $`\sigma =q3`$. It can be proved by direct count of the particles in corresponding points (see fig.3). The dashed lines parallel to $`q^{^{}}`$ axis are determined by the following formula:

$$\sigma =q(BL)$$
(10)
It should be pointed out that formula (2) can be represented by formula (10). For baryons, instead of (10), we have $`\sigma =qB`$. The point $`q=3,\sigma =3,B=6`$contains only one baryon comprised of all the eighteen known quarks. The distribution of colorless particles shows that, in order to have completed the symmetry, it is necessary to assume some particles not observed so far in the points marked by crosses. For example, there should be three particles in the point with coordinates $`q=4,\sigma =3`$because there are three particles (leptons) in the symmetrical point. Now the distribution of colorless particles will be subject to the symmetry transformation described by 4-group of diedr (9). The predicted particles should be located in five points symmetrically positioned to those of leptons and mesons. These particles, like leptons, should not have quark structure because for them $`B=6;7`$. They are supposedly ten times heavier than baryons. We will call them transbaryons. The number of transbaryons is less by one of the total number of leptons and mesons. If, by any selection rule, it is prohibited for some particles to be located in some arbitrary points, then, in compliance with above mentioned symmetry, the same number of particles should be missing out in symmetrically positioned points. But this does not pertain to transbaryons because symmetrical to them leptons and mesons (most of them) have already been observed. The newly proposed symmetry for colorless particles does not depend on generation. Based on this symmetry alone, we cannot judge about forth generation.

Concluding this paragraph, we can say that within the frame of the proposed symmetry for colorless elementary particles the existence of transbaryons seems undisputable .

Fig 1. Distribution of quarks and leptons of the first generation on $`q,\sigma )`$ plane.

Fig 2. Distribution of ordinary particles on $`(q,\sigma )`$ plane. The distance between neighbouring points along axes $`q`$ and $`\sigma `$ is taken as a measurement unit.

Fig 3. Distribution of colorless particles on $`(q,\sigma )`$ plane. Transbaryons are marked by crosses. The figures indicate the number of particles in each point.