{'fraction_non_alphanumeric': 0.04310433973320465, 'fraction_numerical': 0.025455861569749404, 'mean_word_length': 4.394477500929714, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Group classification of (1+3)-dimensional Schr\"odinger equations with position dependent mass

Kinematical invariance groups of the 3d Schr\"odinger equations with position dependent masses (PDM) and arbitrary potentials are classified. It is shown that there exist 94 classes of such equations defined up to the generic equivalence group, and 70 classes defined up to the equivalence groupoid. The maximally extended kinematical invariance algebras of such equations appears to be eight dimensional. The specific symmetries connected with the presence of the ambiguity parameters are discussed and an extended class of systems which keep their forms for arbitrary or particular changes of these parameters is specified. The exact solution of the selected PDM Schr\"odinger equation is presented. This equation describes a deformed 3d isotropic harmonic oscillator and possesses extended continuous symmetries and hidden supersymmetries with two different superpotentials as well.


Introduction
Group classification of Schrödinger equations with arbitrary potentials had been carried out in the seventieth of the previous century in papers [1], [2], [3] and [4]. This classification is a corner stone in the group theoretical grounds of quantum mechanics. It presents a priori information about all admissible Lie symmetries of one particle quantum mechanical systems, and is an important (and necessary) step in investigation of various generalized symmetries of Schrödinger equations.
In contrary, the group classification of Schrödinger equations with position dependent mass (PDM) is still far from the completeness. It can be considered as a challenge for experts in group analysis of differential equations, taking into account the fundamental role played by such equations in modern theoretical and mathematical physics. There are many papers devoted to PDM Schrödinger equations with various particular symmetries, see, e.g., [5], [6], [7], [8]. But the basic symmetries of these equations with respect to continuous transformation groups are still unknown.
The first steps in the systematic study of symmetries of the PDM Schrödinger equations was made apparently in papers [9] and [10]. Namely, in paper [9] the stationary 3d PDM Schrödinger equations with different symmetry groups have been classified. In other words, all such equations admitting first order integrals of motion were presented there.
In paper [10] we discuss symmetry properties of non-stationary PDM Schrödinger equations with arbitrary number of spatial variables and present the completed group classification of such equations with two spatial variables. The number of the equations with non-equivalent symmetry groups appears to be rather restricted. More exactly, there exist only seven classes of such equations, three of which are defined up to arbitrary functions and four -up to arbitrary parameters.
The 3d PDM systems admitting second order integrals of motion were classified in paper [11]. The classification presented there is completed but restricted to the PDM Schrödinger equations invariant w.r.t. the rotation group.
The natural next step is to extend the analysis presented in [9] and [10] to the time dependent PDM Schrödinger equations with three spatial variables, which is the most interested from the physical point of view. Just such equations are discussed in the present paper. We make the completed group classification of them and specify 94 classes of equations with various symmetry groups. In other words, we classify all continuous symmetries which can be admitted by PDM systems and constructively describe all systems with non-trivial symmetries up to equivalence transformations.
The specificity of PDM Schrödinger equations is that they include the ambiguity parameters. We consider also the specific symmetries connected with the presence of the ambiguity parameters and specify an extended class of systems which keep their forms for arbitrary or special changes of these parameters.
In contrast with the case of constant mass, PDM Schrödinger equations have less extended symmetries. However, there is a lot of inequivalent equations admitting six, seven, or eight dimensional symmetry algebras. Such symmetry is sufficient to make equation separable and even exactly solvable.
Quantum mechanical systems with extended Lie symmetries in many cases are also supersymmetric. We do not study supersymmetric aspects of PDM Schrödinger equations here, but consider an example of a PDM system which possesses both the extended Lie symmetry and supersymmetry. This system which we call deformed isotropic oscillator is shape invariant, exactly solvable and has a more general discrete energy spectrum (including an additional spectral parameter) than the standard 3d isotropic oscillator. This system together with many other supersymmetric and exactly solvable ones are included as a particular cases in our classification.
Since the Schrödinger equation with a constant mass appears as a particular case in our analysis, we revised also the classical results presented in the Boyer paper [4] and find few systems which were not presented in this classical work.

Time dependent PDM Schrödinger equations
We will discuss PDM Schrödinger equations of the following generic form where H is the PDM Hamiltonian H = 1 4 (m α p a m β p a m γ + m γ p a m β p a m α ) +V .
Here p a = −i ∂ ∂xa , m = m(x) andV =V (x) are the mass and potential depending on spatial variables x = (x 1 , x 2 , x 3 ), and summation w.r.t. the repeating indices a is imposed over the values a = 1, 2, 3.
The first term in the r.h.s. of equation (2) is interpreted as a kinetic energy term. We use the standard representation with α, β and γ being the ambiguity parameters satisfying the condition α + β + γ = −1 [12].
The choice of values of the ambiguity parameters can be motivated by physical reasons, see a short discussion of this point in [10]. Mathematically, we can treat two Hamiltonians (2) with different values of these parameters as similar provided their potential terms are not fixed. Indeed, Hamiltonian (2) can be rewritten in the form including another values of the ambiguity parameters marked by tildes: (mαp a mβp a mγ + mγp a mβp a mα) + V providedV and V are connected by the following relation: were f = 1 m , f a = ∂f ∂xa and f aa = ∆f = ∂fa ∂xa . In particular, we can chooseα =γ = 0 and reduce Hamiltonian (2) to the following form: We will restrict ourselves to classification of equations (1) with Hamiltonians in the form (5). The classification results can be easily reformulated for systems with generic Hamiltonian (2) using relation (4) withα =γ = 0,β = −1. In other words, potentials of equations in generic form (1) can be expressed via potentials V by the following formula: An interesting special cases of relation (4) appears ifV = V orV = V + Const. If so, the potentials in Hamiltonians (2) and (3) coincide, exactly in the first case and up a constant term in the second. Since Hamiltonians (2) and (3) are equal one to another, their kinematical parts are equal as well. Then the related equation (4) represents equivalence relations for the ambiguity parameters, see Section 8 for discussion of this point.

Determining equations for symmetries
Let us search for symmetries of equations (1), (5) with respect to Lie groups of transformations of dependent and independent variables. Since these equations are linear, they will not be changed if we add to ψ arbitrary particular solutions ψ p . The changes keep the equations invariant and form the infinite group of their symmetries. Our problem is to find all symmetries additional to (7). It can be reduced to searching for the first order differential operators: associated with generators of the continuous symmetry groups of equation (1), (5). In (8) η = 1 2 ξ a a + iη, ξ 0 , ξ a and η are functions of independent variables, whose explicit form can be found from the following operator equation where a is one more unknown function of t and x. The condition (9) generates a system of differential equations for functions ξ 0 , ξ a , η, f,V and a:ξ Here δ ab is the Kronecker symbol, the dot denotes the derivation w.r.t. the time variable whilst derivations w.r.t. the spatial variables are denoted by subindices:ξ 0 = ∂ξ 0 ∂t , f a = ∂f ∂xa , etc.
We will not deduce the system of the determining equations (11)- (14) here since it is nothing but a special case of the system corresponding to the PDM Schrödinger equation with arbitrary number of spatial variables presented in [10].
Equation (11) defines the conformal Killing vector whose generic form is given by the following formula: where r 2 = x 2 1 + x 2 2 + x 2 3 , ε abc is the totally antisimmetric unit tensor while λ a , ω, θ c and ν a are parameters which do not depend on x but in general can depend on t. Notice that, in accordance with (10), function ξ 0 does not depend on x = (x 1 , x 2 , x 3 ).
Considering the remaining determining equations (12)- (14) and using (15), we reduce them to the following form: where the summation is imposed over the repeating indices a over the values a = 1, 2, 3.
Thus the group classification of equations (1) is reduced to finding inequivalent solutions of equations (16)- (18) where ξ a are given by formula (15).

Algebraic content of symmetry operators and equivalence relations
Let us discuss the algebraic structure of symmetry operators (8), (15). It happens these operators can be represented as: where The evident solutions of relations (10)- (14) valid for arbitrary f and V are ξ a = η = 0, ξ 0 = 1. They correspond to the following symmetry operator Q = P 0 : which generates shifts of the time variable. Operators (20) form a basis of the Lie algebra c(3) of conformal group defined in the 3d Euclidean space. This algebra is isomorphic to so(1,4), i.e., to the Lie algebra of the Lorentz group in 1+4 dimensional space. This isomorphism can be fixed by choosing the new basis of algebra c(3): where M µν = −M νµ with µ, ν = 0, 1, 2, 3, 4 are basis elements of algebra so (1,4), satisfying the following commutation relations: were g µν = diag(1, −1, −1, −1, −1). Let us also note that the conformal group C(3) is the equivalence group of Hamiltonian (5) [9]. Indeed, just conditions (11) are necessary and sufficient for keeping the general structure of Hamiltonian (5) while functions f and V can be changed. In other words, transformations of spatial variables belonging to group C(3) do not change the generic form of the Hamiltonian and can be used to simplify it and symmetries (8) as well.
To present a way of such simplification we expand symmetry operators via basis elements of algebra so (1,4): where τ µν are functions of time and the summation over the repeating indices is imposed over the values 0, 1, 2, 3, 4.
We will see that all functions τ µν should be proportional to some fixed function, and so it is possible to reduce the number of nonzero components of tensor M µν using its transformation properties with respect to the group SO (1,4). For example, let τ 0a = 0 for all a = 1, 2, 3, 4, then we can restrict ourselves to symmetries (24) with the only nonzero coefficient τ 12 and then obtain the generic case with nontrivial τ 13 , τ 14 , τ 23 , τ 24 and τ 34 , using a transformation belonging to SO(4)⊂ SO (1,4).
The completed description of all nonequivalent linear combinations of tensors M µν satisfying relations (23) was presented in paper [13]. We use these results to decouple the classification problem to the set of relatively simple subproblems, like it was done in [9].
We will search for potentials V up to constant terms C, i.e., potentials V and V ′ = V + C will be treated as equivalent. To reduce V ′ to V it is sufficient to make the transformation ψ → e −iCt ψ in equation (5). In addition, rescaling the time variable we can reduce to the unity any constant multiplier for f .
Transformations discussed in the above belong to the generic equivalence group of equations (5) and can be applied for any arbitrary elements f and V . However, for some particular masses and potentials there are additional equivalence transformations which are listed in the following formulae: x 1 →x 1 = x 1 cos (ln(ω)) − x 2 sin (ln(ω)) , x →x = e µ λ 2 x, 2 , Transformations (25)- (39) do not belong to the equivalence group C(3) which is valid for arbitrary mass and potential. However, they keep the generic form of equations (1), (5) for some particular f and V . To find these transformations we test the pairwise non-equivalence of all obtained equations, especially of those ones which have equivalent invariance algebras. As a result we discover the additional equivalence transformations (25)-(39) which extend the equivalence group of equations (1) to the equivalence gruppoid.
Let us stress that there are equivalence transformations additional to (25)-(33) which can change arbitrary functions in mass and potential terms to other arbitrary functions. We do not enumerate them here since they are not essential.

Dependence on time
We already know the generic dependence of functions ξ 0 and ξ a on x which is described by equation (19). More sophisticated speculations are requested to define their dependence on t. We will specify the following qualitatively different versions of symmetries.
1. Functions ξ a and ξ 0 are time dependent. In this case we have the generic determining equations (16)-(18) with a being a function of t, whilst η in general depends on t and x. Such symmetries will be called symmetries of class 1.
2. Vector ξ a and function η are time independent but ξ 0 depends on time. In this case function a in equations (17) and (18) is non-trivial. We will refer to the corresponding symmetries as symmetries of class 2.
3. For symmetries of class 3 vectors ξ a are time independent but ξ 0 and η can depend on time. In accordance with (10) in this case function a is nontrivial if ξ 0 t = 0. 4. Symmetries of class 4 -vectors ξ a and their counterparts ξ 0 and η are time independent. In this case the generic expression (15) for ξ a includes arbitrary constant parameters λ a , ν a , θ a and µ. In addition, functions η and ξ 0 in (19) are reduced to constants, and function a in determining equations (17) and (18) is zero.
The fourth class in fact was considered in detail in paper [9], were the symmetries of the stationary Schrödinger equations with position dependent mass have been classified. The only new feature in comparison with [9] is the additional symmetry with respect to shifts of time variables whose generator is given by equation (21).
The second and the third cases are a bit more complicated thanks to the presence of new arbitrary functions of t, i.e., a and η, in the determining equations. However, it can be handled in analogy with the fourth case considered in [9]. To do it we can again exploit the results of Patera and Winternitz [13] concerning the subgroup structure of group P(1,4), which is isomorphic to the equivalence group of equation (1). Any such subgroup can be confronted by a reduced version of the determining equations with a small number of arbitrary parameters λ a , ν a , θ a and ω. Such equations with a = 0 anḋ η = 0 have been integrated in [9] and can be solved also for a nonzero.
The first case looks as much more complicated. However, equation (16) with non-trivial ξ a appears to be a rather strong condition which essentially reduces the number of the corresponding symmetry operators.
Let us start with the analysis of the first class of symmetries. To evaluate the additional restrictions for f generated by equation (16) with non-trivialξ a we differentiate (17) with respect to t and obtain: wheref = ln(f ),f a = fa f . One more rather evident differential consequence of equations (16), (17) and (15) is: Thus if symmetries are time dependent we have the additional conditions for f given above. Formulae (40) and (41) present a system of four algebraic equations for three unknownsf a . Its compatibility condition reads: Equation (42) generates rather strong restrictions on coefficient functions λ a , ν a , θ a and ω. Up to equivalence it admits only three nontrivial solutions for functions ξ a (see Appendix): where Φ = Φ(t) is a function of time. The corresponding functions ξ 0 , and a are: where σ is a constant.
The generic solutions of the first of equations (47) forΦ can be represented in the following form: where a, b and λ are arbitrary constants and the condition ab = 0 reflects the fact that for a nonzero we can reduce b to zero by an appropriate shift of the time variable. We specify the possible time dependence of functions ξ 0 , ξ a , a and η for symmetries of the first type. For symmetries of the second type we should deal with functions ξ a of generic form (15) with constant parameters λ a , θ a , ν a and ω, while ξ 0 is a function of t. The corresponding relations (10) and (17) are consistent iff this function is linear, i.e., where σ is a constant. In accordance with (16) function η is a constant also. For symmetries of the third type we have time independent ξ a discussed in the previous paragraph, and time dependent η which, however, does not depend on x. To satisfy equation (18) function η should be linear in t, i.e., η = µt where µ is a constant.
Finally, for symmetries of the fourth type all functions satisfying the determining equations (10) and (16) Thus we specify the possible dependence of coefficient functions ξ 0 , ξ a and a on time and spatial variables and are in a position to start the procedure of direct solving the determining equations (16)- (18). The results of this solving are presented in the following section while calculational detail can be found in the Appendix.

Classification results
Let us present the results of the group classification of equations (1). The number of such equations with essentially different symmetries is rather extended. The related position dependent masses are either fixed or arbitrary functions of reduced number of variables including parameters.
The admissible fixed mass functions are presented in the following formulae: Just such functions correspond to the most extended symmetries of equation (5). Mass functions defined up to arbitrary parameters also correspond to highly symmetric PDM Schrödinger equations. They can have the following forms: where Θ = arctan( x 2 x 1 ), σ and ν are arbitrary nonzero parameters. Masses with m = r −2+σ and m = r −2−σ are equivalent.
In addition, in Table 6 we specify 24 mass terms defined up to arbitrary functions. Table 1 Systems with power inverse PDM function f =r 2 = x 2 1 + x 2 2 . Parameter λ takes nonzero values, µ can be reduced to zero applying transformation inverse to (31) for Items 1-6 and transformation (25) for Items 7-18. The invariance algebras include the presented symmetries together with generator (21) and the unit operator.
No Potential V Symmetries Invariance algebras There is a big variety of potentials corresponding to PDMs given by equations (52) and (53). The corresponding classification results are arranged in Tables 1-5. Equations (1) with any m and V admit symmetry operator (21). Additional symmetries are indicated in the third columns of the tables where M ab , P a , L a , and D are operators defined in (20) and (22), G(·) are arbitrary functions of the arguments fixed in brackets, and small Greek letters denote arbitrary parameters. Table 2 Systems with power inverse PDM function f =r 2 e σΘ , σ = 0, ω = 0. The invariance algebras include generator (21) and can be extended by adding the unit operator. Potentials presented in Items 6-10 and 11-15 can be obtained from the ones presented in Items 1-5 using transformations (36) and (37) correspondingly. Alternatively, parameter ω can be reduced to unity applying transformation (26). Table 3 Systems with power inverse PDM functions f = r 2 = x 2 1 + x 2 2 + x 2 3 (Items 1-9) and f = r 2+σ , σ = 0, ±1 (Items 10-18). Parameter µ can be reduced to zero using transformation inverse to (32) for Items 1-3 and transformation (27) for Items 4-9 . Potentials and symmetries presented in Items 13,14,15 and 16,17,18 can be obtained from the ones presented in Items 10, 11, 12 using transformations (38) and (39). In addition, parameter ω can be reduced to unity applying transformation (29). (3) Symmetries collected in Tables 1 and 2 are of type 1. The corresponding vectors ξ k have the type presented in (43), (48)-(50). The related symmetry algebras are fixed in the fourth columns, where n a,b and s a,b are nilpotent and solvable Lie algebras of dimension a. We use the notations presented in [14] for low dimension Lie algebras. A discussion of symmetry algebras which are admitted by equation (1) can be found in Section 7.
In Table 3 we present symmetries of type (8), (44). The maximally extended symmetry algebras are seven dimensional, see Items 3, 6 and 9 there. However, these algebras and all other algebras presented in Items 1, 2, 4, 5, 7, 8 ultimately include the unit operator which is accepted by any of the considered equations. The invariance algebras presented in Items 10-18 can be extended by adding the unit operator. Algebra n 1,1 includes the only basis element L 3 .
Notice that the low dimension algebras with dimension d ≤ 5 and a certain class of the algebras with dimension 6 had been classified by Mubarakzianov [15], see also a more contemporary and accessible papers [16], [17] and [18].
In Table 4 some exotic PDMs are specified. In particular, there are two systems which admit eight dimensional invariance algebras, see Items 6 and 7. However, the related symmetries are time independent and commute with Hamiltonians. In other words, they belong to integrals of motion for the stationary Schrödinger equation which have been classified in paper [9]. 2 Table 4 Systems with some particular mass functions. The invariance algebras include the presented symmetries together with generator (21) and the unit operator. Parameter µ can be reduced to zero using transformation inverse to (33).
No Inverse mass f Potential V Symmetries Invariance algebras All the invariance algebras presented in Items 1-5 of Table 4 are solvable. Like in the previous tables we use notations s a,b for solvable algebra of dimension a and classification number b, see [14].
One dimension algebra n 1,1 in Items 1, 2 is spanned on L 3 and the unit operator correspondingly. In Item 4 the basis elements of s 1,1 ⊕ s 1,1 are L 3 and the unit operator. In Items 6 and 7 the symbol s 1,1 ⊕ s 1,1 denotes the direct sum of one dimensional algebras spanned on P 0 and the unit operator.
Symmetries of type (8), (45) are collected in Table 5. The corresponding masses are exponentials in the only spatial variable x 3 . The invariance algebras again have dimensions 3, 4, 5 or 6. These algebras are either simple of form a direct sum of a simple and solvable algebra, and can be extended by the unit operator. By e(2)∼ s 3,3 we denote the Euclidean algebra in two dimensional space, whose generators are P 1 , P 2 and L 3 , while n 1,1 denotes the one-dimensional algebras whose generators are presented in the last positions of Column 3. Table 5 Systems with exponential inverse PDM function f = e x 3 . Potentials and symmetries presented in Items 5, 6, 7, 8 and 9, 10, 11, 12 can be obtained from the ones presented in Items 1, 2, 3, 4 using transformations (34) and (35). In addition, parameter ω can be reduced to unity applying transformation (29), and equivalence transformation (30) is available which generates multiplier σ for x 3 .
No Potential V Symmetries Invariance algebras In Table 6 rather generic systems are presented. They include potentials (and masses) defined up to arbitrary functions, but admit reduced numbers of symmetries belonging to class 3 specified in Section 5. In Table 6 the symbols F (·) and G(·) denote arbitrary functions of arguments fixed in brackets, D 1 F and D 2 F are the derivations of function F = F r 2 +1 r , νΦ + Θ with respect to the first and second argument correspondingly, F ′ is the derivation of function F = F r 2 ±1 r with respect to its combined argument r 2 ±1 r . In addition, to save a room we denote Φ = arctan r 2 −1 2x 3 and Ψ = 3νx 3r 2 (r 2 +1) 2 −4r 2 . Table 6. Systems with masses and potentials defined up to arbitrary functions No Inverse mass f Potential Symmetries We indicate symmetries additional to P 0 and I in the third column. For the systems presented in Items 1, 3, 5, 7, 9 and 2, 4, 6, 8, 10 the corresponding symmetry algebras are s 2,1 ⊕ n 1,1 and s 3,1 correspondingly. In Items 11, 13, 15, 17, 19 and 21 we have symmetry algebras isomorphic to s 2,1 ⊕ n 1,1 ⊕ n 1,1 if at least one of parameters λ or σ is nontrivial. If all these parameters are equal to zero, these symmetry algebras (and algebras admitted by systems presented in Items 12,14,16,18 and 20) are degenerated to the direct sums of four one dimensional algebras.
In Items 12, 14, 16, 18 and 20 we have symmetry algebras isomorphic to s 3,1 ⊕ n 1,1 if at least one of parameters λ or σ is nontrivial. Finally, the systems presented in Items 22, 23 and 24 admit the direct sums of symmetry algebras e(2) ⊕ n 1,1 , so(3)⊕n 1,1 and so(1,2)⊕n 1,1 correspondingly. A short discussion of the invariance algebras is presented in the following section.
Thus we find all 3d Hamiltonians (5) with PDM which correspond to equations (1) with nonequivalent symmetries. The presented list of equations (1) is completed up to equivalence transformations belonging to group C(3). The additional equivalence transformations are indicated in the tables headings. In all tables parameters µ and ω can be reduced to zero using additional equivalence transformations (25)-(39), like the isotropic harmonic oscillator can be reduced to the free particle Schrödinger equation. We keep these parameters in the classification tables for the readers convenience.
The number of such equations is rather extended. Namely, we specify 22 classes of equations defined up to arbitrary parameters and 70 classes of equations defined up to arbitrary functions, see Tables 1-6. In addition, there are two systems with the fixed mass and potential terms presented in Items 6 and 7 of Table 4.
We did not consider equations with constant mass terms, since they had been classified long time ago. However, some of such equations appears as particular cases of our analysis, and the comparison of these cases with well known results leads to a rather non-excepted conclusion: the classification results presented in [4] are incomplete. A discussion of this point is presented in Section 9.

Algebras of symmetry operators
The complete sets of Lie symmetries of a given partial differential equation should have a structure of a Lie algebra. In particular, it is the case for symmetries presented in Tables 1-6. Let us describe these structures explicitly.
First we present commutation relations between essentially time dependent symmetry operators and generator P 0 . For symmetries presented in Table 1 and the first half of  Table 3 we obtain: where I is the unit operator, a = 1, 2 and no sum with respect to the repeating indices. Thus for systems indicated in Items 1, 7, 13 of Table 1 and Items 1, 4, 7 of Table  3 we have four dimensional symmetry algebras, which are solvable and include the unit operator as a central element. In accordance with [14], we use the notations n 4,1 for algebras (54) and n 4,7 for algebras (55) and (56).
The additional symmetries which are presented in the remaining items of the mentioned tables amend the mentioned algebras to five dimensional algebras n 5,4 or s 5, 16 . In addition, in Items 6, 12 and 18 of Table 1 the direct sums of solvable five dimensional algebras with the simple algebra sl(2,R) are indicated. The latter is a linear span of basis elements K 3 , P 3 and D.
The additional symmetries presented in Tables 2, 3 and 5 commute with Q k a , N k a and S k a and satisfy the following commutation relations between themselves: Relations (60), (61) and (62) specify algebras sl(2,R), so(3) and e(2) correspondingly.

Equivalence relations for ambiguity parameters
The subjects of our classification are equations (1) including PDM Hamiltonians of special form (5). In other words, we restrict ourselves to the a priori fixed values of the ambiguity parameters Nevertheless the obtained results can be easily reformulated for arbitrary values of these parameters. To this effect it is sufficient to change the potentials presented in the tables to the corresponding effective potentialsV given by relation (6). Let us note that the ordering ambiguity for the kinetic energy term is an interesting problem which attracted attention of numerous researchers. The ordering (5) was proposed in [19] to ensure the current conservation. The ordering with α = γ = − 1 2 was applied in [20] for reformulating the connection rule problem on the two sides of an heterojunction. The same ordering was derived in [21] as the nonrelativistic limit of the Dirac Hamiltonian with PDM, and in [22] via path integral evaluation. There are other ordering versions suitable for particular physical models: α = γ [23], β = α = 0, γ = −1 [24], γ = 0, α = β = − 1 2 [25], α = γ = − 1 4 , β = − 1 2 [26]. All the corresponding Hamiltonians (2) are equivalent to (5) with effective potential (6) which in general differs from the initial potential V .
Let us consider an abstract Hamiltonian (2) and transform it to the form (3), i.e., make the change α →α, β →β, γ →γ,V → V (64) were potentialsV and V are connected in the way given by formula (6). We say Hamiltonian (2) is ambiguity invariant if the transformed potential V is equal to the initial potentialV up to a constant term:V = V + C. In accordance with (6) it happens if inverse mass f satisfies the following condition: Equation (65) gives the necessary and sufficient conditions of the complete equivalence of two generic Hamiltonian (2) and (3) with the same inverse mass functions f = 1 m and potentials but different ambiguity parameters. In particular, forβ = −1 andα =γ = 0 this condition is reduced to the following form while the second Hamiltonian is reduced to the form given by equation (5). For all the inverse masses f represented in Tables 1-5 there are multiparameter families of the ambiguity parameters satisfying condition (66) with C = 0. For the cases included into Tables 1, 2, 5 and position 2 of Table 4 the parameters equivalent to (63) have to satisfy the following condition 2αγ + α + γ = 0 (67) which is a consequence of (66) with C = 0. Formula (67) defines a family of the ambiguity parameters, one of which, say, α, can take arbitrary values different from − 1 2 , whereas the others are can be expressed via α in the following manner: Thus we have a one parametric family of the ambiguity parameters for which Hamiltonians (2) and (3) are ambiguity equivalent provided their mass and potential terms are presented in Tables 1 and 3. The same is true for the systems presented in Table 3, however the condition (67) should be changed by the following one: 2αγ(σ + 2) + (α + γ)(σ + 3) = 0 were σ can take both nonzero and zero values, the latter case corresponds to the systems presented in Items 1-9.
For the systems presented in Table 1 and Items 10-18 of Table 3 there exist a more general family of equivalent ambiguity parameters with β = −1 − α − γ and arbitrary α and γ. In this case condition (67) includes a non-trivial C and the corresponding Hamiltonians (2) and (3) differ by a constant term, i.e., V =V + C which is not essential and can be removed by the change ψ → e −iC ψ in equation (5).
It is interesting to note that some systems presented in the classification tables can be treated as purely kinematical. This means that the corresponding Hamiltonians can be represented in generic form (2) with suitable ambiguity parameters and trivial or constant potential V .
First, it is the case when potentials do not include parameters, which can play the role of coupling constants, see Items 6 and 7 of Table 4. The corresponding Hamiltonians can be rewritten in the form (3) withσ =γ = − 1 2 ,β = 0 andṼ = 0. The same property have Hamiltonians which include masses and potentials specified in Items 4 and 5 of this table, Item 5 of Table 2, Item 12 of Table 3 and Item 4 of Table 5, in spite of that they include coupling constant κ.
Moreover, all potential terms in Tables 2-5 including parameter κ can be transferred into the kinetic part of the Hamiltonians with some special values of the ambiguity parameters. These values can be easily found solving the corresponding equation (66) where C is changed by the potential term with the multiplier −κ. In particular, for the Hamiltonian whose potential term is specified in Item 5 of Table 2 we obtain the following condition for the ambiguity parameters: (σ + 2) (2αγ(σ + 2) + (α + γ)(σ + 3)) = −4κ. (68) All PDM Hamiltonians (2) with the ambiguity parameters satisfying (68) and trivial potentialV = 0 can be rewritten in form (5) with V = κr σ .

New symmetries of 3d Schrödinger equation with constant mass term
It is generally accepted to think that symmetry classification of Schrödinger equations with constant masses has been completed long time ago in Boyer paper [4]. Let us show that, in order to be complete, the Boyer classification should be extended to include some systems missing in [4]. We considered equations with essentially nonconstant masses and ignored the systems with f = Const. However, some of systems with constant masses are presented implicitly in the classification tables. Thus, the systems presented in Table 2 include PDMs defined up to arbitrary parameter σ = 0. However, they are well defined for σ = −1, but the corresponding mass functions used in Items 10-18 are reduced to constants.
All PDMs presented in Table 6 are defined up to arbitrary functions, which in particular can be chosen as constants. In this way, taking if necessary σ = −2, many classes of PDM Schrödinger equations specified in Table 6 can be reduced to equations with constant mass terms.
Thus Tables 3 and 6 specify implicitly a number of constant mass Schrödinger equations which admit non-equivalent symmetries, and it is possible to compare these results with the results of Boyer classification [4]. Doing this, it is possible to see that there are at least three cases missing in [4]. Namely, setting in Items 4, 8 and 12 of Table 6 F = 1 and substituting the corresponding expressions for f and V into (5) we come to the following Hamiltonians: Equation (1) with Hamiltonian (69) admits the symmetry: while for Hamiltonian (70) we have: Finally, equation (1) with Hamiltonian (71) admits tree symmetry operators, namely, (72) and (74): Symmetries Q,Q andQ are fixed in Table 6. Symmetry G is not present there since it appears only in the case of a constant mass.
Equations (69)-(71) are missing in Boyer classification results, the same is true for symmetries of type (72). However, such symmetries do exist, and this fact can be verified directly, calculating commutators of operators (72) and (73) with L defined in (1). Such commutators are identically equal to zero.

A little on supersymmetry
Some of the presented equations admit rather extended symmetries and so have good chances to be exactly solvable. In many cases the extended symmetry guaranties the existence of solutions in separated variables, i.e., the equation is separable or even multiseparable.
In this section we discuss one more nice property of some of the presented equations. This property is called shape invariance, which can be considered as a reason for the exact solvability. Just shape invariant systems can be solved algebraically using tools of SUSY quantum mechanics.  Table 3:

Sharp invariance with oscillator potential
where we denote σ = 2ν. Equation (75) admits rather extended symmetries being invariant w.r.t. six parametrical Lie group. We will see that, in addition, it admits a hidden supersymmetry.
In view of the rotational invariance of equation (75) and its symmetry with respect to shifts of time variable, it is possible to search for its solution in separated variables.
In particular, we can use the spherical variables and search for solutions in the following form: where θ and ϕ are angular variables and Y lm (θ, ϕ) are spherical functions, i.e., eigenvectors of L 2 = L 2 1 + L 2 2 + L 2 3 and L 3 . As a result we come to the following radial equation for R lm : This equation can be simplified using the Liouville transform: which reduce (77) to the following form: where we denote Equation (79) describes a deformed 3d harmonic oscillator which includes two deformation parameters ν and κ. To justify such interpretation consider first the special case when these parameters are constrained by the following condition: and equation (79) is reduced to the following form: Equation (81) is shape invariant. Indeed, the Hamiltonian H r can be factorized: On the other hand, the superpartnerĤ l of Hamiltonian (82) has the following property which fixes the shape invariance: Shape invariant equations can be solved algebraically using tools of SUSY quantum mechanics, see, e.g., [27]. In particular, for equation (81) the ground state energy level E = E 0 is equal to 1 2 C l , while for the n-th exited state we obtain E n = 1 2 C l+νn , or where δ = ν − 1. Equation (85) represents the spectrum of 3d isotropic harmonic oscillator deformed by the term proportional to δ. For integer δ the spectrum is degenerated, but this degeneracy is much less than in the case of the standard 3d oscillator which corresponds to δ = 0.
Equation (79) is shape invariant also in the general case when condition (80) is not imposed. The corresponding energy spectrum is: whereκ = 8(κ + 1) + ν(ν + 3). The related eigenvectors are expressed via the confluent hypergeometric functions F : where n is integer and E n is the corresponding eigenvalue given by equation (86). Thus we effectively solve equation (75) using tools of SUSY quantum mechanics. However, there is one more way to do it, as it will be shown in the next section.

Shape invariance with Morse potential
Let us return to equation (77) and solve it using more sophisticated approach. First we divide all terms by r 2ν and obtain: In equation (88) we change the roles of the coupling constant and energy and consider it as an eigenvalue problem for ε. Making the Liouville transform we simplify (88) to the following form: Equation (91) includes the familiar Morse potential and so is shape invariant. Indeed, denoting µ = ω(ν + ν 2 ) we can factorize hamiltonian H ν like it was done in (82) where index l should be changed to ν and W = ν − ωe −aρ , C ν = ν 2 (93) and the shape invariance is easy recognized. To find the admissible eigenvalues ε and the corresponding eigenvectors we can directly use the results presented in paper [27], see Item 4 of Table 4.1 there: where y = ω ν r −ν . Thus we find the admissible values ofε n . Using definitions (89) and (92) we can find the corresponding values of E which are in perfect accordance with (86).

Discussion
It is an element of common knowledge that PDM Schrödinger equations are not invariant with respect to Galilei transformations, see discussion of this point in [28]. A natural question is: what do we have instead, what kinds of symmetry do exist in PDM systems?
The results presented above in Sections 6 and 7 give a constructive answer to this question. They include the completed list of continuous symmetries which can be admitted by PDM Schrödinger equations. In particular, the systems presented in Items 1-4 of Table 5 admit Galilei transformation of variable x 3 , provided the time variable is changed also. Indeed, integrating the Lie equations corresponding to generators Q 3 1 , we obtain the following symmetry transformations for the independent variables: Solving the Lie equations, it is not difficult to find also the corresponding transformation low for the wave function ψ, but we omit the related cumbersome formula.
Thus we extend the well known results [4] to the case of Schrödinger equations with position dependent mass. The number of PDM systems with non-equivalent symmetries appears to be much more large than in the case of constant masses. However, the maximal admissible symmetries are less extended. In particular, we cannot find PDM systems invariant w.r.t. ten parametrical Galilei group or twelve parametrical extended Galilei group which are admissible by the constant mass Schrödinger equation with harmonic oscillator, linear and constant potentials. Nevertheless, there are PDM systems admitting eight, seven or six parameter symmetry groups, see Section 7 for discussion of this point. In addition, thanks to the linearity of the considered equations, there is the infinite symmetry group of transformations (7).
As other extensions of results of paper [4] we can mention the group classification of the nonlinear Schrödinger equations [29] and the analysis of its conditional symmetries [30].
To find the determining equations for symmetries of equations (2), (5) we use the traditional technics applied in papers [1]- [4], i.e., evaluate commutators of the searched symmetry operators with i∂ t − H, see Section 3. The more general approach which is applicable also to nonlinear equations was developed long time ago by Sophus Lie, see, e.g., the fundamental Olver monograph [31]. But we choose a more simple way which does not request a knowledge of the group analysis grounds from potential readers of the present paper. Let us note that this way is applicable also for group analysis of some non-linear equations [32], [33].
However, the deducing of the determining equations is a necessary but rather simple step in the group classification of a given class of equations. Much more efforts are requested to obtain the complete set of inequivalent solutions of the determining equations. To achieve this goal we apply the algebraic approach which includes the description of subgroups of the equivalence group of equations (1) and a priory analysis of the admissible invariance groups. The grounds of the algebraic approach were created in papers [34], but in fact its elements were used for the first time in paper [35].
To obtain the completed group classification of equations (1) it was necessary to describe also the equivalence transformations which keep the generic form of these equations but can change the explicit form of arbitrary functions f and V . It is well known that in general the equivalence transformations do not form a Lie group, but have a structure of groupoid [36]. This is the case also for equations (1), and for some subclasses of these equations there exist additional equivalence relations which do not belong to their equivalence group. These relations are specified in the tables heads and can be used to transform some of the presented systems to more simple ones, like the isotropic harmonic oscillator can be reduced to the Schrödinger equation with the trivial potential.
A special attention was paid to the equivalence relations of the ambiguity parameters. In Section 8 we present the notion of the ambiguity invariance of PDM Hamiltonians which is defined as the invariance w.r.t. changes of the ambiguity parameters up to a constant potential term. In this section extended classes of Hamiltonians having this property are specified.
Thanks to their extended symmetries many of the presented systems are exactly solvable. In Section 10 an example of a solvable system is discussed. In addition to the symmetry under the six parameter Lie group whose generators are given in Item 18 of Table 3 and formula (21), equation (75) possesses a hidden dynamical symmetry w.r.t. group SO (1,2). The effective radial hamiltonian is shape invariant, and its eigenvalues can be easily found with tools of SUSY quantum mechanics. We call this system deformed 3d isotropic harmonic oscillator. In spite on the qualitative difference of its spectra (85) and (86) of the standard 3d oscillator, it keeps the main supersymmetric properties of the latter.
We discuss SUSY aspects only of one selected system. But in fact the classification tables present a number of supersymmetric systems with hidden dynamical symmetries. In particular, there are other systems with supersymmetric radial equations, systems with supersymmetric equations in angular variables and systems with supersymmetries of both the mentioned types. To keep the reasonable paper size we do not discuss SUSY aspects of systems other than equation (85).
The very possibility to solve shape invariant Schrödinger equation with using more than one superpotential is well known, see, e.g., [37]. We show that it is possible to do it with the PDM system (85) by changing the roles of the energy eigenvalue and coupling constant.
In the present paper we did not consider Schrödinger equations with constant mass. However, such equations can appear as particular case of systems presented in the classification tables. Examining these cases we found some systems missing in Boyer classification [4], refer to Section 9 in the above.
It is necessary to stress that symmetries of the constant mass Schrödinger equation cannot be completely described in frames of the more general problem of group classification of PDM Schrödinger equations, since the latter includes much more extended equivalence group which reduces the number of non-equivalent symmetries. Thus the revision of classical results presented in [4] is requested, but it deserves a separate publication which is in preparation.
In the present paper the continuous symmetries of PDM Schrödinger equations are discussed. We do not consider discrete symmetries which present additional powerful tools in construction of supersymmetric and exactly solvable models with constant masses [38], [39]. Moreover, the discrete symmetry operators can serve as constructive elements of physically consistent potentials [40], [41]. A systematic study of discrete symmetries of PDM Schrödinger equations is one more interesting field which still waits for researchers attention.

A.1 Time dependence of symmetry operators
Here we analyze the compatibility condition (42) and prove that the time dependence of symmetry operators with non-trivialξ a is described by formulae (43)-(45) and (48)-(50).
Let us differentiate (15) w.r.t. t and substitute the obtained expression forξ a into (42). Then, equating the coefficients for the same powers of x a we obtain the following conditions for functionsλ a ,ν a ,θ a , a andω: λ a =ν a = 0,ω = 0, a, b = 1, 2, 3.
(A1) or, alternatively, Considering the case (A1) we come to the following form of the doted symmetry operator:Q Let conditions (A2) are valid then functionsλ a andν a can be represented in the following forṁ λ a = gμ a ,ν a = qμ a (A4) where g, q andμ a are some functions of t. In accordance with (15), (A2) and (A4) the corresponding vectorξ a is reduced to: and can be correlated to the following symmetry operator (8) differentiated with respect to time: where relations (22) were used. Making hyperbolic rotations on the plane 0-4, it is possible to reduce (A6) to one of the following forms: where µ a = 2μ a (g 2 + q 2 ) 1 2 . Thus our analysis of time dependent symmetries is reduced to considering the four independent cases presented by formulae (A3), (A7), (A8) and (A9).
Let us start with symmetries defined by relation (A3). The corresponding vectorξ a is: which leads to the following equations (40) and (41): Solving (A11) forf a we obtain: Two evident consequences of equation (A12) are: In order to systems of equations (A12) and (A13) be compatible, all functionsθ 1 ,θ 2 ,θ 3 andȧ have to be linearly dependent, i.e. proportional to the same function of t. Thus, up to rotation transformations, we can restrict ourselves to the case whenθ 1 =θ 2 = 0,ȧ = aθ 3 and reduce equations (A11) to the following form: 2 is a function of t. Thus the version (A3) can be effectively reduced to the case when vectorξ a is given by equations (43) where we denote θ =Φ.
Let us consider the next version of doted symmetry operators given by equation (A3). The corresponding vectorξ a looks as followṡ and so equations (40) and (41) are reduced to the following system: Multiplying (A17) byξ b and summing up w.r.t. index b we obtain: f aξ bξb =ξ a (ȧ + 2ω).
Differentiating equation (A18) w.r.t. t and taking into account time independence of f a , we obtain the following differential consequence: ξ bξb (ξ a (ȧ + 2ω) +ξ a (ä + 2ω)) = 2ξ a (ȧ + 2ω)ξ bξb . (A19) Substituting (A15) into (A19) and equating the coefficients for the same powers of independent variables x a we come to the following system of equations for functions µ a , ω and τ =ȧ + 2ω: τω =ωτ ;μ a τ = µ aτ . (A20) In accordance with (A20) functions µ a ,ω and τ should either be zero or proportional one to another. Since zero value of τ corresponds to constant masses (see equation (A18)), we suppose τ = 0. Then equations (A20) present the following alternatives: where C 0 and C a , a = 1, 2, 3 are constants. For both versions (A22) and (A23) the corresponding functions (A15) can be reduced to the form given in (44). To this effect it is sufficient to make the shifts x a → x a − Cȧ ω . If version (A21) be realized, then up to rotation transformations we can restrict ourselves to functionsξ a given by formula (45).
Consider now the remaining versions (A8) and (A9). In both cases the corresponding vectorsξ a (A5) can be represented in the following unified form: were δ = ±1. The corresponding equations (40) and (41) are reduced to the following system: Solving these equations forf a we obtain: Differentiating (A26) with respect to t and and taking into account time independence off a , we obtain the following condition: (μ a R + µ aṘ + x aṖ )ξ bξb = 2(µ a R + x a P )ξ bξb . (A28) Substituting expressions (A27) into (A28) and equating the coefficients for the same powers of x a we conclude that all functions µ a andω should be proportional to the same function of t which we denote asΦ, i.e., The corresponding operators (A8) and (A9) are: Using the rotation transformations we can reduce nontrivial linear combinations C a M 0a of vectors M 0a to CM 03 with C 2 = C 2 1 + C 2 2 + C 2 3 . Then making a planar rotation on the plane 3-4, it is possible to reduce the linear combination C 0 M 04 + CM 03 toCM 04 withC 2 = C 2 + C 2 0 . As a result we reduce (A30) to the case with trivial C a , and the corresponding vectorξ a again takes the form given by equation (44).
Analogous speculations can be applied to operator (A31) which can be reduced to the case of trivial C a provided C 2 0 >C 2 . If C 2 0 <C 2 we can reduce this operator to the form (A3), while for C 2 0 =C 2 it is reducible to the particular case (A7) withω = 0. Thus we prove that the classifications of symmetries with time dependent Killing vectors ξ a can be reduced to solving the determining equations (16)-(18) withξ a enumerated in formulae (43), (44) and (45).

A.2 Symmetries with time dependent derivative terms
Let us start with symmetries of the first class when functions ξ a are time dependent and have one of the forms presented in (43)-(45).
For functions (43) the determining equations (16)- (18) are reduced to the following form: The generic solution of equation (A33) looks as follows: f = e νΘ F (r), a = νΦ wherer = x 2 1 + x 2 2 , Θ = arctan x 2 x 1 , ν is a constant, and F is an arbitrary function ofr. This solution is compatible with equation (A32) iff F = µr 2 with some constant µ, which can be reduced to the unity by rescaling the time variable. In this way we obtain f =r 2 e νΘ , η = 1 νΦ e νΘ + φ(t), ν = 0.
Substituting (A35) into (A34) we obtain: This equation specifies admissible functions Φ and φ. Indeed, to make it consistent we have to seẗ where κ and ν are constants. Then substituting (A37) into (A36) we obtain the following equations for potential V : Since ν is nonzero, it is possible it is possible to reduce ν to zero by a constant shift of the potential. In addition, the generic solutions of the first of equations (A37) is given by formulae (48)-(50).
Formulae (A35), (43), (48)-(50) present explicitly all components of symmetry operators (8) together with the corresponding inverse mass function f . Solving equation (A38) consequently for zero, negative and positive κ we obtain the corresponding potentials: In the special case σ = 0 relations (A35) are transformed to the following form: Functions Φ and φ should satisfy relations (A37) , and so equation (A42) is solved by the following function: Thus we have specified all Hamiltonians which admits symmetries of special form (8), (43). Any of such Hamiltonians admits the following representation: where arbitrary parameter σ can take both the zero and nonzero values and possible potentials V are given by formulae (A39), (A40), (A41) and (A44). The corresponding functions ξ 0 present in symmetry operator (8) are easily calculated solving equation (10). Such specified Hamiltonians include arbitrary functions G(r, x 3 ). It is important to note that for some fixed G(r, x 3 ) the corresponding equation (1) Tables 1 and 2. In complete analogy with the above we can solve determining equations (16)-(18) for special cases of functions ξ 2 presented in (44) and (45). The corresponding classification results are presented in Tables 3 and 4

A.3 Symmetries of classes 2 and 3
For symmetries (8) with time independent coefficients ξ a we can directly use the results of paper [13] concerning the subgroup structure of group SO(1,4) and decouple the system of determining equations (16)- (18) to inequivalent subsystems corresponding to selected values of arbitrary parameters in vector (15), like it was done in [9]. The only new feature in comparison with [9] is the necessity to take into account additional time dependent functions ξ 0 and a which appear now in the determining equations.
Let us start with one dimensional subalgebras of algebra c(3)∼so (1,4). In accordance with [13] it is possible to specify five inequivalent subalgebras of the mentioned type spanned on the following basis elements: Subalgebra spanned on P 3 is associated with the following symmetry operator Q = ξ 0 + ∂ 3 + iη where ξ 0 and η are functions of t and t, x respectively. The corresponding determining equations (10) and (16)- (18) are reduced to the following system: f 3 = af, a = −ξ 0 , η a = 0, V 3 = aV +η. (A47) Taking into account time independence of f and V , we conclude that ξ 0 and η are linear functions of t, and a is a constant: ξ 0 = νt + µ, η = κt + ρ, a = −ν We can restrict ourselves to the case ρ = µ = 0 since these parameters are coefficients for symmetries accepted by any equation (1). In addition, if σ is nontrivial, we can set κ = 0 since this parameter can be removed transforming Q → e −i κ σ t Qe i κ σ t . Thus effectively there are only two versions: ξ 0 = νt, η = 0 and ξ 0 = 0, η = κt. (A48) Solutions of the corresponding equations (A47) are: Just these solutions are represented in Items 1 and 2 of Table 6. In complete analogy with the above we can solve the determining equations corresponding to the other one dimension algebras presented in (A46).
All algebras whose basis elements are presented in (A51) include one dimensional subalgebras fixed in (A46). Functions f and V which correspond to the one dimensional subalgebras are known and presented in Items 1-10 of Table 6. Let f (1) and V (1) correspond to the first elements of pairs given in (A51). Our task is to substitute these functions into the determining equations generated by the second elements and integrate the obtained system.
For algebras presented in the first line of equation (A51) we have functions f = f (1) and V = V (1) given by equations (A49) and (A50). The determining equations corresponding to the second basis element Q = P 2 can be obtained from (A47) by changing the subindex 3 to 1.
We have again two possibilities fixed in (A48), so the second symmetry operator can be reduced to one of the following forms: Q = P 1 + iµt∂ t , orQ = P 1 + νt.
There are four different pairs Q,Q with Q andQ given in equations (A49), (A50) and (A52). Moreover, passing to a new basis in algebras Q,Q and making the suitable linear transformations of variables x 1 and x 2 it is possible to reduce them to the following two representatives: P 3 + κt, P 1 + iσt∂ t and P 3 + iσt∂ t , P 1 + κt . functions f and V have to solve the following systems of equations: for functions (A50), and for functions (A49). As a result we obtain solutions represented in Items 15 and 16 of Table 6. Consider now the second algebra from the list (85). Its basis element L 3 is associated with the following symmetry operator Q = ξ 0 + x 1 ∂ 2 − x 2 ∂ 1 + iη and the following determining equations (17) and (18): where we use the angular variable Θ = arctan( x 2 x 1 ). Like in (A47) a andη should be constants, and, by definition, f and V are given by formulae (A49) or (A50). We rewrite these formulae using radial and angular variables: and so functions f and V are reduced to the form presented in Item 11 of Table 6. Then, substituting (A57) into (A55) we recover functions f and V presented in Item 12 of Table  6.
In analogous way we solve the determining equations corresponding to the remaining pairs of symmetries represented in (A51). As a result we obtain functions f and V enumerated in Items 13, 14 and 17-21 of Table 6.
The specification of systems admitting more extended symmetry algebras presented in Table 5 and Items 22-24 of Table 6 can be made in analogy with the above with using three dimensional and more extended subalgebras of algebra so(1,4) which can be found in paper [13], see also [9]. We will not present here the routine calculations requested to realize this programm.