Alena A Dobretsova

Alexei D Chepelianskii

Nikolay N Mikhailov

Don Ze

Kvon

Z D Kvon

Spin mixing between subbands and extraordinary Landau levels shift in wide HgTe quantum wells

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The recent discovery of topological insulators, a new state of mater, in materials with strong spin-orbit interaction [1][2][3][4][5][6] has opened an exciting research direction. In most topological insulators the strong relativistic effects and particularly spin orbit coupling lead to an inverted spectrum which is at the origin of their topological properties. Since relativistic effects strength grows with increasing atomic mass, investigations have focused on a broad range of materials with Bi and Sb [7][8][9][10][11], InAs/(In)GaSb [12][13][14] and HgTe [3,15,16] quantum wells, all containing elements with a large atomic numbers in Mendeleev periodic table. The quantum wells based on HgTe and InAs / GaSb can be grown epitaxially and the carrier density can be tuned with a gate in a broad range, these systems are thus particularly suited for fundamental investigations on topological insulators. One of the principal feature of HgTe quantum wells is that their properties drastically depend on their thickness. By changing the width d of the quantum well [17] one can get: an ordinary insulator like GaAs quantum wells at d < 6.3 nm [3,15,18]; a linear spectrum of so called Dirac fermions at the critical thickness d ≈ 6.3 nm [19,20]; inverted spectrum of 2D topological insulator at d > 6.3 nm [3,15,16,18]; 2D semimetallic spectrum with a small overlap between the conduction and valence bands at d ≈ 20 nm [21,22]; the spectrum of 3D topological insulator at width d > 70 nm [23][24][25].

In addition to its impact on the band structure, spin orbit interaction can appear directly in magneto-transport experiments through the shift of Landau levels (LLs). Previous works [26][27][28] focused on Landau levels shift due to the Rashba splitting arising as a consequence of the asymmetric deformation of the quantum well at high gate voltages. In this work we present magneto-transport experiments and a thorough theoretical analysis showing that the shift of Landau levels can also occur due to mixing of quantum well subbands with different spin. Indeed the average spin of the quantum-well wavefunction is directly connected to Berry phase since both describe the transformational properties of the wavefunction under plane rotation. A change of the wavefunction spin thus directly leads to an anomalous shift of Landau levels due to the contribution of Berry phase in the semiclassical quantization rule for Landau levels [29][30][START_REF] Berry | The quantum phase, five years after in geometric phases in physics[END_REF][START_REF] Littlejohn | [END_REF][33][34][35][36][37]. Following these theoretical works we develop a quantitative semi-classical model which reproduces the Landau levels obtained numerically in a realistic 6-band k • p theory calculation, these theoretical results are then compared to the experiment.

EXPERIMENT

Our experiments were carried out on 20 ÷ 22 nm undoped Cd 0.6 Hg 0.4 Te/HgTe/Cd 0.6 Hg 0.4 Te heterostructures with (013) surface orientation, grown by molecular beam epitaxy. The cross section of the structure is illustrated in Fig. 1 (a). The detailed description of samples preparation can be found in [22,38]. These quantum wells have a separation of 65meV between p-like well conduction subband [H1] with total spin (angular momentum) 3/2 and the first excited s-like subband [E2] with spin 1/2. This splitting is sufficiently large to resolve a large number of Landau levels in the conduction band at few Tesla magnetic fields but still within the range covered by changes in the gate potential, 20 ÷ 22nm wide quantum wells were thus optimal for our experiments. The evolution of the conduction band spin as a function of momentum k is shown in the inset in Fig. 1 (b), the average spin in z direction indeed changes from 3/2 at low energies to around 1 near the bottom of the subband [E2] as a consequence of their mixing. The previous works on similar HgTe samples [21,[39][40][41][42][43][44] were mostly concentrated on the energy region when the conduction and valence bands overlap (see Fig. 1 (b) leading to coexistence of electrons and holes. It has been shown in InAs/GaSb quantum wells [45] that mixing between plike electrons and s-like holes states can already lead to a non-trivial LLs behavior. In our case both the first conduction and valence subbands are p-like, in the work we thus concentrated on the energy region highlighted in red in Fig. 1 (b), which correspond to a metal regime where the two dimensional electron gas is formed only by electrons of the conduction band [H1]. To perform magnetotransport measurements, the structures were patterned into Hall bars with width 50 µm and the distance between the contacts 100 and 250 µm. The structures had also a metallic top gate allowing to reach the metallic regime and to change electron densities from 10 10 cm -2 to 8 • 10 11 cm -2 . Electron mobility in this electron concentration range varied from 5 to 30 m 2 /Vs and measurements were performed by standard lock-in technique at temperature T = 0.2 ÷ 2 K and magnetic field up to 9 T.

The typical experimental results are shown in Fig. 2, where as example results for 22 nm HgTe quantum well are presented. In Fig. 2 (a) 2D color map of longitudinal resistivity ρ xx versus magnetic field B and gate voltage V g is shown. Here the energy axis calculated from V g by using 6-band Kane model is also present. In the figure one can see that some non-degenerate neighbor Landau levels (resistivity maxima) coincide with each other with magnetic field increase. One of the example of such LLs behavior is highlighted by dots: "red" level (highlighted by red dots) at small finite magnetic field B ∼ 1 T is degenerate with "blue" and separated from "black" ones however with B increase it comes closer to the "black" one till their degeneration at B = 2.4 while the 'blue" level become separate.

The same pattern can be also seen in Fig. 2 (b), where the fixed magnetic field cross-sections of ρ xx from Fig. 2 (a) and Hall resistance ρ xy are shown for B from 1.0 to 2.4 T. Here Landau levels highlighted by dots in Fig. 2 (a) are also highlighted by dots of the same color. In the dependence of longitudinal resistivity on gate voltage the described effect is seen as the transformation of two peaks at B = 1 T: one corresponding to the "blue" and "red" levels and one separated to "black" one, into two different peaks at B = 2.4 T: one corresponding to "red" and "black" levels and now "blue" one is separated. To simplify the visualization we have joined the dots of the same color. In Hall resistance dependence on gate voltage one can also see that at B = 1 T there is Hall plateau corresponding to filling factor 14 while at 2.4 T it disappears. Similar results were obtained in [46], but the microscopic explanation in term of Berry phase change due to spin mixing between subbands was not presented.

THEORY

To describe the unusual Landau levels behavior observed in the experiment we have calculated the energy spectrum of (001) 20 nm HgTe quantum well with infinite wall boundaries in magnetic field by means of 6band Kane model as in [17] with Kane parameters taken from [47] (see Fig. 1 (b)). In order to calculate the position of Landau levels as function of the perpendicular magnetic field, we used an isotropic model for the quantum well. This assumption enabled us to find eigenstates using only a small number of nearby oscillator wavefunctions and thus to solve the problem without diagonalizing very large matrices, we checked by full diagonalization that the anisotropic terms deform the spectrum only weakly. Our numerical results for the Landau levels of the first conduction band versus magnetic field are shown in Fig. 3. The calculation gives the energy levels (highlighted by yellow, green and red lines) with the same behavior as was observed in the experiment: yellow and green levels non-degenerate at small magnetic field come closer to each other at large B. Note that in this calculation we didn't take into account the well asymmetry arising as a consequence of the applied gate voltage, the calculated Landau levels shift are thus an intrinsic property of the quantum well spectrum. The inclusion of an asymmetric contribution is discussed in the Appendix D of the manuscript and doesn't change the qualitative picture.

The physical origin of this unusual shift of the Landau levels can be understood from the semiclassical quantization rules. In semiclassical approximation Landau level energy can be found from the well-known Lifshitz-Onsager equation [START_REF] Kittel | Quantum theory of solids[END_REF] A(E n ) = (n + γ)2πeB/c, where A is an area of a cyclotron orbit cross-section in k-space (constant-energy and k z = const surfaces intersection, z -the direction of magnetic field B); n -orbital Landau level number; γ -some phase from 0 to 1. In the case of a free electron γ is simply equal to 1/2, while in general V g (V) case in solids:

γ = 1 2 + γ s = 1 2 - 1 2π Γ(C) + 1 2π Γ N G .
We denote by γ s the deviation of γ from the free electron value, we will show that it is connected to the average spin of electrons. This quantity has several contributions, the first term is the Berry phase Γ(C) = i C u k | ∂u k ∂k dk -written as a function of the envelope of the electron wave functions across the well |u k as function of the in plane momentum k (the semiclassical analysis uses wavefunctions without magnetic field). The second term Γ N G is a phase related to the electron orbital magnetic moment [29-34, 36, 37]. We decided to designate Γ N G as "non-geometric" as in [START_REF] Littlejohn | [END_REF], since it depends not only on the shape of the semiclassical trajectory in phase space but also on the rate at which the orbit is traversed. For completeness we provide an elementary derivation of the general expression for Γ N G in the Appendix C, while here only its expression for the case of an isotropic spectrum will be discussed.

In the isotropic approximation, the envelope wavefunction |u(k, θ) as function of the momentum amplitude k and of its polar angle θ can be expressed as function of |u(k, θ = 0) by application of the total angular momentum operator Ĵz which describes rotations along z direction perpendicular to the two dimensional electron gas: |u(k, θ) = e -3iθ/2 e -iθ Ĵz |u(k, 0) . The phase factor e -3iθ/2 is introduced to make the wavefunctions single valued upon a full rotation θ → θ + 2π, it is necessary since J z is a direct sum of spin 1/2 and 3/2 operators. This relation can be formally established from the identity e -iθ Ĵz Ĥ(k)e iθ Ĵz = Ĥ(R z (θ)(k)) where R z (θ) is the rotation operator around z-axis by an angle θ. Using this expression for |u(k, θ) , we can connect the Berry phase accumulated along a cyclotron trajectory in magnetic field as a function of the average spin J z in z direction:

Γ(C) = i dθ u(k, θ)|∂ θ |u(k, θ) = 2π J z + 3π. (1)
We see that for an isotropic quantum well the Berry phase is directly connected to the spin projection along z direction.

In isotropic approximation the non-geometric phase can be simplified to

Γ N G = 2π d /dk Im ∂ θ u| Ĥ(k, θ) -(k)|∂ k u , (2)
where Ĥ(k, θ) is the Hamiltonian of the quantum well as a function of the polar coordinates of the momentum k

and (k) the dispersion relation of the conduction band.

Due to the isotropy of the quantum well this expression is actually independent on θ and was numerically evaluated at θ = 0. For a symmetric quantum well (without Rashba) spin up and spin down carriers have the same energy, making the absolute value of the spin projection in the conduction band a well defined function of energy with |J z | changing from 1.5 to 0.9 in the explored energy range (see inset of Fig. 1 (b)). Both Berry and the non-geometric phases change sign for opposite spins. To make γ s a well defined function of energy we take its positive branch for the spin-up orientation keeping in mind that the sign is reversed for the spin-down branch, we also take the values of all phases in the (0, 2π) interval. For the following discussion it is convenient to introduce the dimensionless ratios ∆ n↑ between Zeeman splitting of Landau levels with the same orbital index n and opposite spin and the orbital spacing between nearby Landau levels with spin down orientation. Indeed keeping in mind the symmetry of γ s and the semiclassical quantization rule these ratios can directly be expressed as function of γ s :

∆ n↑ = (E n↑ -E n↓ ) (E (n+1)↓ -E n↓ ) = 2γ s , (3)
we remind that n is the orbital Landau level index. Equation (3) allows us to compare the Landau-levels obtained in the experiment and in our 6-band Kane model simulations with the semiclassical values for γ s from Eqs. (1,2).

For spin down we define ∆ n↓ = (E n↓ -E n↑ )/(E (n-1)↑ -E n↑ ), with this choice both ratios ∆ n↑ and ∆ n↓ collapse semi-classically on the same dependence which simplifies future discussions.

DISCUSSION

Fig. 4 (a) compares the theoretical semi-classical dependence 2γ s (E) with the ratios ∆ n↑ and ∆ n↓ obtained for Landau levels from numerical diagonalisation of the isotropic 6-band Kane model -these results show an excellent agreement between the semi-classical theory and numerical Landau levels already for n ≥ 4. This figure also shows that the energy dependence of the total phase Γ(C) + Γ N G qualitatively follows the Berry phase and is approximately equal to Γ(C)/2 + const. This approximate relation can viewed as a generalization to topological quantum wells of the exact cancellation (up to a topological winding number) between the two phases which occurs for electron-hole symmetric two band model [34,35]. Since for HgTe this compensation is only approximate, the anomalous shift of the Landau levels still follows the Berry phase providing a magnetotransport probe for the vertical spin projection of the quantum well wavefunctions.

The experimental data for ∆ ν = ∆ n↑,↓ are presented in Fig. 4 (b) for different filling factors ν, the position of the energy levels was extracted from the position of ρ xx maxima in Fig. 2 (a). It is seen that while the 6-band Kane model provides a good qualitative description of the experimental energy shifts (for e.g. in Fig. 3) the agreement is not quantitative. The agreement improves at high ν 15 however for small ν ≤ 6, ∆ ν is independent on energy. It is possible that in the low filling factor regime transport across incompressible stripes becomes important and make transport more sensitive to density gradients at the edges of the sample [START_REF] Dahlem | [END_REF]. This would make our assumption to extract the bulk energy levels position from ρ xx maxima invalid. On the theory side inter-Landau level exchange interactions [50,51] that were not included in the Kane model can maybe also give a substantial ν dependent correction to the energy shifts. In summary we have observed non-trivial Landau levels shifts with magnetic field in wide ∼20 nm HgTe quantum wells where Landau levels split in magnetic field and become degenerate again when magnetic field increases. Theoretical calculations within a 6-band Kane model allowed us to reproduce this picture qualitatively. We proposed a microscopical explanation based on the semiclassical quantization of the Landau levels: this shift is caused by a mixing of states from two quantum well subbands with different spin: the conduction band with spin 3/2 and the second well subband with spin 1/2. This mixing results in a change of the average electron spin with energy, and thus a change of the Berry phase which is directly connected to the spin in an isotropic quantum well. Taking into account the Berry phase and nongeometrical phase related to the electron orbital magnetization, we reproduced theoretically the ratio between the Zeeman and orbital Landau level splittings obtained = E|u 0 (x) e iS0(x)/ , where p x = ∂ x S 0 . We find that |u 0 (x) is an eigenvector of Ĥ(x, p x ); x and p x follow Hamilton equations with an Hamiltonian given by the corresponding eigenvalue Ĥs (x, p x )|u 0 (x) = (x, p x )|u 0 (x) .

We would like to go to higher order and seek the solution in the form |ψ = A|u 0 (x) e iS0(x)/ + |u 1 (x) e iS0(x)/ (|u 0 (x) and |u 1 (x) are normalized to unity, u 1 (x)|u 0 (x) = 0).

Note: we denote by triangular brackets | the matrix multiplying only, without integration. For further calculations it is useful to write the following identities:

e -iS0(x)/ (-i ∂ x ) A(x)e iS0(x)/ |u 0 (x) = S 0 (x)A(x)|u 0 (x) -i ∂ x [A(x)|u 0 (x) ]
e -iS0(x)/ (-i ∂ x ) 2 A(x)e iS0(x)/ |u 0 (x) = ( 10)

A(x)|u 0 (x) S 0 (x) 2 -i S 0 (x) -2i ∂ x [A(x)|u 0 (x) ] S 0 (x) + O( 2 ).
Now using (10) we compute the matrix element:

e -iS0/ u 0 | Ĥ|ψ(x) = e -iS0/ u 0 | V0 (x) + V1 (-i ∂ x ) + V2 (-i ∂ x ) 2 A(x)e iS0/ |u 0 = u 0 | Ĥs |u 0 A(x) -i V2,00 A(x)S 0 (x) -i V1,00 ∂ x A(x) -i2 V2,00 ∂ x A(x)S 0 (x) -i A(x) u 0 | V1 ∂ x |u 0 -2i A(x)S 0 (x) u 0 | V2 ∂ x |u 0 = (x, p x )A(x) -i V2,00 A(x)S 0 (x) -i u 0 | ∂ Ĥs ∂p x ∂ x [A|u 0 ] ,
where we introduced the operator:

∂ Ĥs ∂p x = V1 + 2 V2 p x
and the real (since V 1,2 are self adjoint) matrix elements V1,00 = u 0 | V1 |u 0 , V2,00 = u 0 | V2 |u 0 .

On the other hand since |ψ(x) is an eigenvector the following property takes place:

e -iS0/ u 0 | Ĥ|ψ(x) = EA.
On semiclassiclal trajectories (x, p x ) = E. Replacing the differentiation with respect to x to the differentiation with respect to time t we are thus led to the equation on the slowly varying complex amplitude A:

V2,00 A(x) dp x dt + u 0 | ∂ Ĥs ∂p x d dt [A|u 0 ] = 0. ( 11
)
We now want to separate the equations on the phase φ and modulus A of the complex quantity A = Ae iφ , for this purpose we rewrite Eq. ( 11) in the form:

V2,00 A(x) dp x dt + ∂ ∂p x dA dt + u 0 | ∂ Ĥs ∂p x d dt |u 0 A = 0, ( 12
)
where we used u 0 | ∂ Ĥs ∂px |u 0 = ∂ ∂px (that can be obtained by differentiating over p x the identity u 0 | Ĥs |u 0 = ). In Eq. ( 12) the coefficient before dA dt is real, what allows us to write an equation on the phase φ:

V2,00 A(x) dp x dt + ∂ ∂p x (A + iφ A) + u 0 | ∂ Ĥs ∂p x d dt |u 0 A = 0.
Further separating real and imaginary parts we find:

dφ dt = - ∂ ∂p x -1 Im u 0 | ∂ Ĥs ∂p x d dt |u 0 . (13)
To transform this equation on the phase into the usual form, we notice the relation:

∂ Ĥs ∂p x |u 0 = ∂ ∂p x Ĥs |u 0 -Ĥs ∂|u 0 ∂p x = ∂ ∂p x |u 0 + (E -Ĥs )| ∂u 0 ∂p x .
Since ∂ Ĥs ∂px is self adjoint we use this relation to transform Eq. ( 13) into:

dφ dt = -Im u 0 | d dt |u 0 - ∂ ∂p x -1 Im ∂u 0 ∂p x |(E -Ĥs ) d dt |u 0 . (14)
Eq. ( 14) can be simplified further by noticing that:

d dt |u 0 = ∂|u 0 ∂x dx dt + ∂|u 0 ∂p x dp x dt = ∂|u 0 ∂x ∂ ∂p x - ∂|u 0 ∂p x ∂ ∂x .
The last term gives a real contribution which does not contribute to the phase evolution, we are thus finally lead to the equation on the phase from Littlejohn [START_REF] Littlejohn | [END_REF]:

dφ dt = -Im u 0 | d dt |u 0 -Im ∂u 0 ∂p x |(E -Ĥs )| ∂u 0 ∂x = i u 0 | d dt |u 0 + i 2 ∂u 0 ∂p x |( Ĥs -E)| ∂u 0 ∂x - ∂u 0 ∂x |( Ĥs -E)| ∂u 0 ∂p x . (15)
Here the first term is the well known Berry phase, while the second one does not have the name. It can be integrated along orbits just as well as the Berry's-phase term, but the result in general cannot be represented as a line integral of a differential form, thus it is not "geometrical" in the same sense as Berry's phase term and thus can be called as "non-geometric" phase. Further we will use this name to denote the second term.

APPENDIX C: SEMICLASSICAL QUANTIZATION FOR HGTE QUANTUM WELL Without magnetic field, the Hamiltonian of the quantum well takes the form:

Ĥqw = Ĥqw (k x , k y ).
In magnetic field in Landau gauge the Hamiltonian becomes:

Ĥ = Ĥqw (-i ∂ x , eBx).
Rotation invariance gives the property (which can be checked by the direct calculation):

e -iθ Ĵz Ĥqw (k)e iθ Ĵz = Ĥqw (R z (θ)(k)),
where J z is z-component of the total angular momentum, R z (θ) is the rotation operator around z-axis to the angle θ.

We can thus write the wave function as:

|u(k, θ) = e -3iθ/2 e -iθ Ĵz |u(k, 0) , (16)
where factor phase e -3iθ/2 is introduced to make wavefunctions single valued upon full rotation θ → θ + 2π (it is necessary because J z is direct sum of spin 1/2 and 3/2 and Γ(C) in range 0 ≤ Γ(C) < 2π. The Berry phase Γ B (the first term in Eq. ( 15)), is thus simply connected with the spin of the wavefunctions:

Γ B = i dθ u(k, θ)| d dθ |u(k, θ) = dθ u 0 | Ĵz + 3 2 |u 0 = 2π J z + 3π,
where we introduced |u 0 = |u(k, 0) . The non-geometric contribution (the second term in Eq. ( 15)) can be also simplified taking into account the rotation invariance (16):

Γ N G = dtγ N G .
We first write γ N G in cylindrical coordinates:

γ N G =
γ N G = eB 2p u 0 | Ĵz + 1 2 ( Ĥs -E)| ∂u 0 ∂p + ∂u 0 ∂p |( Ĥs -E) Ĵz + 1 2 |u 0 ,
where the fact that Ĵz , e -iθ Ĵz and Ĥs all commute with each other was used. We thus find that γ N G does not depend on θ but only on the momentum amplitude. Since γ N G is constant along the semiclassical trajectory the total non-geometric phase is:

Γ N G = 2π ω c γ N G ,
where ω c = eB p d dp is the (semiclassical) cyclotron frequency at energy E = (p). Finally we find: (c) Landau levels position as the function of the perpendicular magnetic field calculated with (orange dots) and without (cyan dots) taking into account the well asymmetry due to gate voltage. The blue lines show the Landau levels with non-trivial behavior that were highlighted in Fig. 3 in the main text. The red lines correspond to the same Landau levels, but calculated with considering the well asymmetry due to gate voltage. The qualitative behavior of red lines including well anisotropy is similar to the blue ones, but a bit further from experimental Landau levels behavior. Panel (d) displays the ratio between Zeeman and orbital level splittings ∆ n↑/↓ (dots) as a function of electron energy taking into account the well asymmetry induced by the perpendicular electric field. Since the conduction band is spin-split (see Panel b) the universal scaling that observed in Fig. 4 (a) from the main text is no longer strictly followed.

Γ N G = 2π Re u 0 | Ĵz + 1 2 ( Ĥs -E)| ∂u 0 ∂ = 2π Im ∂u ∂θ | Ĥs -E| ∂u ∂ .
The electron density in the quantum wells was changed applying a top gate voltage. It is known that this results in an asymmetrical well deformation, which leads to the formation of spin-polarized bands due to the Rashba effect, which in turn leads to a shift of the Landau levels [27,52]. Here we show that the Rashba effect does not change much the qualitative picture of Landau levels behavior, which was discussed in the main text, and does not help to improve the quantitative agreement with the experiment. A rigorous approach to treat the gate voltage effect requires solving self-consistent Poisson and Schrödinger equations, in our case we used the simple approximation and modeled the effect of the top gate as a linear potential gradient inside the well (in other words a constant electric field). The gradient is zero at energies corresponding to the bottom of the conduction band and increases linearly with electron energy to 25 mV at energy E = 60 meV corresponding to electric field of 12.5 kV/cm. The slope has been chosen in the way that calculated Rashba splitting at energy 56.3 meV (corresponds to gate voltage of 7 V) was the same as obtained experimentally from Fourier analysis of Shubnikov -de Haas oscillations (see Fig. 5 (a)). The energy spectrum of the conduction band calculated with and without taking into account the well asymmetry due to top gate voltage is shown in Fig. 5 (b).

The comparison of the numerical results for Landau levels of the first conduction band versus magnetic field with (cyan dots) and without (orange dots) taking into account the Rashba effect is shown in Fig. 5 (c). An ensemble of three Landau levels highlighted by blue lines shows the non-trivial behavior where Landau levels split and come closer to each other again when magnetic field increases, it corresponds to the energy levels highlighted by lines in Fig. 3 in the main text. Red lines correspond to the same Landau levels, but calculated with considering the Rashba spin-splitting. The qualitative behavior of "red" Landau levels is close to blue lines, however at high magnetic field two upper red lines are more separate than the blue ones. Rashba effect thus does not enhance the intrinsic effect due to the mixing of the conduction band with the second well subband, but makes it a bit weaker.

Also as it can be seen from Fig 5 (d) the "universal" scaling from Fig. 4 (a) of the main text, where the ratio between Zeeman and orbital Landau level splitting collapsed on a single curve given by the doubled electronic phase 2γ s , is no longer valid in presene of Rashba interaction. This occurs because the derivation of the scaling relation assumed the same band structure for both spin orientations which is no longer correct when Rashba interaction is present. We see that due to Rashba interaction the dependence of ∆ n↑/↓ on energy becomes weaker with Landau level number increase, which is at odds with the experimental results (see Fig. 4 (b) in the main text) where this dependence becomes instead stronger at higher Landau level numbers.

To conclude this part, we have shown that Rashba interaction does not change the qualitative picture presented in the main manuscript, quantitatively the induced correction seem to show the opposite trend compared to the experimental results suggesting that other physical effects need to be included.

APPENDIX E: COMPARISON OF NON-ISOTROPIC AND ISOTROPIC LANDAU LEVELS CALCULATIONS

In order to calculate the position of Landau levels as a function of the perpendicular magnetic field, we used an isotropic model for the quantum well. This assumption enabled us to find eigenstates using only a small number of nearby oscillator wavefunctions and thus to solve the problem without diagonalizing very large matrices. In Fig. 6 the comparison of the calculations with (violet points) and without (red crosses) isotropic approximation is presented. In non-isotropic case 40 Landau levels were included in the calculation what made the diagonalizing matrix size being 24 000. The deviation of the calculation with isotropic approximation from it in non-isotropic case is seen being small and decrease at high Landau level numbers.

FIG. 2 .

 2 FIG. 2. (a) 2D color map of longitudinal resistivity ρxx versus magnetic field B and gate voltage Vg for 22 nm undoped HgTe quantum well. The energy axis calculated from Vg by using 6-band Kane model is also showed. Dots highlight the example of Landau levels with non-trivial behavior: "red" (highlighted by red dots) and "black" levels being split at B = 1 T degenerate at B = 2.4 T. (b) Longitudinal resistivity ρxx and Hall resistance ρxy versus Vg at magnetic field from 1.0 to 2.4 T. Dots correspond to the same Landau levels as in 2D color map (a). In the ρxx dependence on Vg y-axis corresponds to B = 1.0 T, the each next curve is shifted up by 200 Ω.

FIG. 4 .

 4 FIG. 4. Comparison between the doubled deviation 2γs of electron phase from free electron value 1/2 and the dimensionless ratios ∆ν in theory and experiment. (a) Dots show ∆ n↑/↓ from 6-band Kane model in isotropic approximation; blue line -Berry phase itself (not doubled) Γ(C) = 2π < Jz > +3π/2 and the black line is 2γs -the sum of Berry phase ΓC and non-geometrical phase ΓNG. (b) Shows the values of ∆ν obtained from the experiment compared to theoretical 2γs (black line). Thick green line with square symbols gives ∆ν for the Landau level highlighted by red dots in Fig. 2.

FIG. 5 .

 5 FIG. 5. (a) Fourier spectrum of longitudinal conductivity at gate voltage of 7 V. Due to Rashba splitting two peaks at frequencies 14.1 and 16.8 T are seen. Inset: longitudinal conductivity versus inverse magnetic field. Arrows indicate beatings due to Rashba splitting. (b-d): Calculations for (001) 20 nm HgTe strained quantum well by isotropic 6-band Kane model. (b) Energy spectrum of the conduction band without (black dashed line) and with (violet lines) taking into account well deformation arising as a consequence of the applied top gate voltage. The gate voltage effect was modeled as a linear well bottom slope growing linearly from zero at conduction band minimum to 25 mV ( electric field 12.5 kV/cm) at energy 60 meV.(c) Landau levels position as the function of the perpendicular magnetic field calculated with (orange dots) and without (cyan dots) taking into account the well asymmetry due to gate voltage. The blue lines show the Landau levels with non-trivial behavior that were highlighted in Fig.3in the main text. The red lines correspond to the same Landau levels, but calculated with considering the well asymmetry due to gate voltage. The qualitative behavior of red lines including well anisotropy is similar to the blue ones, but a bit further from experimental Landau levels behavior. Panel (d) displays the ratio between Zeeman and orbital level splittings ∆ n↑/↓ (dots) as a function of electron energy taking into account the well asymmetry induced by the perpendicular electric field. Since the conduction band is spin-split (see Panel b) the universal scaling that observed in Fig.4(a) from the main text is no longer strictly followed.

FIG. 6 .

 6 FIG.6. Landau levels position as a function of the perpendicular magnetic field in (001) 20 nm HgTe quantum well calculated by 6-band Kane model with (violet points) and without (red crosses) isotropic approximation. The deviation of the calculation with isotropic approximation from it in non-isotropic case is seen being small and decrease at high Landau level numbers.

  FIG.3. Landau levels position as a function of the perpendicular magnetic field in (001) 20 nm HgTe quantum well calculated by 6-band Kane model in isotropic approximation. Orange, green and red lines show the example of Landau levels with non-trivial behavior as in the experiment: orange and green levels come closer to each other with magnetic field increase. n is orbital Landau level number, arrows symbolize level spin directions.

ACKNOWLEDGEMENTS

We thank O.E. Raichev for his assistance in the implementation of the 6-band Kane model and S. Gueron and J.N. Fuchs for helpful discussions. This work was supported by the Russian Science Foundation (Grant No. 16-12-10041) and ANR grant SPINEX.

APPENDIX A: HAMILTONIAN OF THE SYSTEM

We used the 6-band Kane model for HgTe whose matrix Ĥ(k) form can be found in [17,47], the matrix form however does not explicitly show that the Hamiltonian is rotation invariant in the case when γ 2 = γ 3 = 0.9 (). Rotation invariance can be seen from the identity e -iθ Ĵz Ĥ(k)e iθ Ĵz = Ĥ(R z (θ)(k)), where R z (θ) is the rotation operator around z-axis by an angle θ, but it requires tedious calculations to be established.

We have thus found it useful to write the HgTe Hamiltonian as function of spin operators which make rotation invariance manifest. For this purpose we introduced spin 1/2 operators Ŝ and spin 1 operators L, their sum Ĵ = Ŝ + L can then describe both spin 3/2 and 1/2 bands from the 6-band Kane model As a function of Ŝ and L operators the Hamiltonian (in the isotropic approximation, and here without strain for simplisity) then reads:

where

) ;

(5)

E g , P, B, γ 0 , γ 12 = (γ 1 + γ 2 ) / 2, γ 23 = (γ 2 + γ 3 ) / 2 -Kane parameters of the system taken from [47].

The rotation invariance of the Hamiltonian in the form Eq. ( 4) is manifest as all the terms are writen as scalar products of vectors. Most of the terms have a familiar form in the point of view of spin physics: 0 (k 2 ) is a global energy offset, J(k 2 ) is an exchange interaction and D is a fine structure term describing an anisotropy in the k direction. The unusual part (from the spin-point of view) comes from the last "Zeeman" term proportional to P which has contribution in the form of an exchange energy dependent g-factor k • L Ŝ • L + c.c.. It seems likely that those terms are created by the underlying second order perturbation theory expansion in the k • p theory.

APPENDIX B: SEMICLASSICAL QUANTIZATION OF THE MULTICOMPONENT SCHROEDINGER EQUATION

Semi-classical quantization rules have been derived in several previous publications, however we haven't found an elementary derivation in the WKB spirit that would directly apply to the multi-component wave equations that stem from k • p theory. For example the derivation in [START_REF] Littlejohn | [END_REF] uses the Wigner transforms formalism, this makes the derivation conceptually appealing, but also less transparent. We have chosen an approach which directly mirrors the standard derivation of the WKB approximation for the usual Schöedinger equation. This allows us to make the derivation self-contained, and we have chosen to present it here.

Let us assume the multicomponent 1D Schroedinger equation in the following form:

Here all the operators V i are self-adjoint; also for simplicity to avoid problems due to non-commutation between px = -i ∂ x and Vi we have chosen V1,2 to be independent on x.

For the next calculation it is useful to introduce:

Assuming the wave function amplitude change much slower than the phase, we seek the solutions in the form:

Substituting ( 9) into (8) and expanding it into powers of , to lowest order in (since in this case derivatives act only on the exponent e iS(x)/ ) we have: Ĥ|ψ(x) = Ĥs (x, p x )|u 0 (x) e iS0(x)/ + O( )


Mohamad Houda

Paul Lescot

TWO COMPUTATIONS CONCERNING THE ISOVECTORS OF THE BACKWARD HEAT EQUATION WITH QUADRATIC POTENTIAL

We determine the isovectors of the backward heat equation with quadratic potential term in the space variable. This generalize the calculations of ). These results first appeared in the first author's PhD thesis (Rouen, 2013).

Backward Equation

The fundamental equation of Euclidean Quantum Mechanics of Zambrini is the backward heat equation with potential V (t, q):

     θ 2 ∂η V u ∂t = -θ 4 2 ∂ 2 η V u ∂q 2 + V (t, q)η V u (C (V ) 1 )
η V u (0, q) = u(q) where t represents the time variable, q the space variable, and θ is a real parameter, strictly positive (in physics, θ = √ ). This equation is not well-posed in general, but existence and uniqueness of a solution are insured whenever the initial condition u belongs to the set of analytic vectors for the operator appearing on the right-hand side of the equation (see e.g. [START_REF] Flato | Simple facts about analytic vectors and integrability[END_REF], Lemma 4, p. 429).

We shall denote by η u (t, q) = η 0 u (t, q) the solution of the backward heat equation with null potential :

   θ 2 ∂ηu ∂t = -θ 4 2 ∂ 2 ηu ∂q 2 (C (0) 1 )
η u (0, q) = η 0 u (0, q) = u(q). The authors wish to thank the referees for constructive comments.

1.1. Generalization of results due to Lescot-Zambrini (cf. [START_REF] Lescot | Probabilistic deformation of contact geometry, diffusion processes and their quadratures[END_REF]).

Theorem 1.1. Let a(t), b(t) and c(t) be continuous functions and

V (t, q) = a(t)q 2 + b(t)q + c(t).
For all initial conditions u, the solution η V u (t, q) of:

θ 2 ∂η V u ∂t = - θ 4 2 ∂ 2 η V u ∂q 2 + V (t, q)η V u (C (V ) 1 )
such that η V u (0, q) = u(q) is given by: η V u (t, q) = ϕ 1 (t, q) η u (ϕ 2 (t, q), ϕ 3 (t, q)) , (1.1) Date: 01 November 2019.

where: η u (t, q) = η 0 u (t, q) and ϕ 1 (t, q), ϕ 2 (t, q), ϕ 3 (t, q) only depend on a(t), b(t) and c(t) (via formulas (1.7), (1.8), (1.9), (1.11), (1.12), (1.13) and (1.14)).

Proof. We shall make formula 1.1 explicit with the following initial conditions:

   ϕ 1 (0, q) = 1 ϕ 2 (0, q) = 0 ϕ 3 (0, q) = q.
In this case : η V u (0, q) = η u (0, q) = u(q). We now differentiate formula 1.1 with respect to t and q; η u and its derivatives will be taken in (ϕ 2 (t, q), ϕ 3 (t, q)).

∂η V u ∂t = ∂ϕ 1 ∂t η u + ϕ 1 ( ∂η u ∂t ∂ϕ 2 ∂t + ∂η u ∂q ∂ϕ 3 ∂t ) , ∂η V u ∂q = ∂ϕ 1 ∂q η u + ϕ 1 ( ∂η u ∂t ∂ϕ 2 ∂q + ∂η u ∂q ∂ϕ 3 ∂q ) , ∂ 2 η V u ∂q 2 = ∂ 2 ϕ 1 ∂q 2 η u + 2 ∂ϕ 1 ∂q ∂η u ∂t ∂ϕ 2 ∂q + ∂η u ∂q ∂ϕ 3 ∂q + ϕ 1 ∂ 2 η u ∂t 2 ( ∂ϕ 2 ∂q ) 2 + ∂ 2 η u ∂q 2 ( ∂ϕ 3 ∂q ) 2 + 2 ∂ 2 η u ∂t∂q ∂ϕ 2 ∂q ∂ϕ 3 ∂q + ∂η u ∂t ∂ 2 ϕ 2 ∂q 2 + ∂η u ∂q ∂ 2 ϕ 3 ∂q 2 and V (t, q)η V u = a(t)q 2 ϕ 1 η u + b(t)q ϕ 1 η u + c(t) ϕ 1 η u .
So, it's enough to have:

θ 2 ϕ 1 ∂ϕ 2 ∂t -θ 2 ϕ 1 ( ∂ϕ 3 ∂q ) 2 + θ 4 2 (2 ∂ϕ 1 ∂q ∂ϕ 2 ∂q + ∂ 2 ϕ 2 ∂q 2 ) = 0 (1.2) θ 4 2 ϕ 1 ( ∂ϕ 2 ∂q ) 2 = 0 (1.3) θ 2 ϕ 1 ∂ϕ 3 ∂t + θ 4 2 (2 ∂ϕ 1 ∂q ∂ϕ 3 ∂q + ϕ 1 ∂ 2 ϕ 3 ∂q 2 ) = 0 (1.4) θ 4 ∂ϕ 2 ∂q ∂ϕ 3 ∂q = 0 (1.5) θ 2 ∂ϕ 1 ∂t + θ 4 2 ∂ 2 ϕ 1 ∂q 2 -ϕ 1 a(t)q 2 + b(t)q + c(t) = 0 . (1.6)
As ϕ 1 (0, q) = 1, the equation (1.3) gives us:

( ∂ϕ2 ∂q ) 2 = 0, then ϕ 2 = ϕ 2 (t)
and (1.5) is then automatically satisfied. So, the equation (1.2) implies:

ϕ 1 (t, q)( ∂ϕ 2 ∂t -( ∂ϕ 3 ∂q
) 2 ) = 0, then for all (t, q) :

∂ϕ 2 ∂t = ( ∂ϕ 3 ∂q ) 2 ∂ϕ 3 ∂q = A(t) ϕ 3 (t, q) = A(t)q + B(t) (1.7) and ϕ 2 (t) = t 0 A 2 (s) ds , (1.8)
then the equation (1.4) gives :

ϕ 1 ∂ϕ 3 ∂t + θ 2 ∂ϕ 1 ∂q ∂ϕ 3 ∂q = 0 ϕ 1 Ȧ(t)q + Ḃ(t) + θ 2 A(t) ∂ϕ 1 ∂q = 0 ϕ 1 (t, q) = k(t) e -1 θ 2 A(t) Ȧ(t)
2 q 2 + Ḃ(t)q , (1.9) therefore :

∂ϕ 1 ∂t = k(t) + k(t) θ 2 A(t) 1 2 ( Ȧ2 (t) A(t) -Ä(t))q 2 + ( Ȧ(t) A(t) Ḃ(t) -B(t))q e -1 θ 2 A(t) Ȧ(t) 2 q 2 + Ḃ(t)q , ∂ϕ 1 ∂q = -1 θ 2 ϕ 1 A(t) ( Ȧ(t)q + Ḃ(t))
and

∂ 2 ϕ 1 ∂q 2 = 1 θ 4 ϕ 1 A 2 (t) Ȧ2 (t)q 2 + 2 Ȧ(t) Ḃ(t)q -θ 2 Ȧ(t)A(t) + Ḃ2 (t) ,
then the equation (1.6) becomes :

k(t) - 1 2 Ä(t) A(t) + Ȧ2 (t) A 2 (t) -a(t) q 2 + k(t) - B(t) A(t) + 2 Ḃ(t) Ȧ(t) A 2 (t) -b(t) q + θ 2 k(t) -k(t)c(t) + k(t) 2 Ḃ2 (t) A 2 (t) -θ 2 k(t) 2 Ȧ(t) A(t) = 0.
As k(t) = 0, the equation (1.6) is equivalent to the following system :

               - 1 2 Ä(t) A(t) + Ȧ2 (t) A 2 (t) -a(t) = 0 - B(t) A(t) + 2 Ḃ(t) Ȧ(t) A 2 (t) -b(t) = 0 θ 2 k(t) -k(t)c(t) + k(t) 2 Ḃ2 (t) A 2 (t) -θ 2 k(t) 2
Ȧ(t) A(t) = 0
(1.10) knowing that :

           A(0) = 1 B(0) = 0 k(0) = 1 Ȧ(0) = 0 Ḃ(0) = 0
This is due to initial conditions at time t on the functions ϕ 1 (t, q), ϕ 2 (t, q), and ϕ 3 (t, q). Let

A(t) = e t 0 ρ(s)ds , (1.11) then ρ(t) = Ȧ(t) A(t) and ρ(t) + ρ 2 (t) = Ä(t) A(t) ,
therefore the first equation of (1.10) is equivalent to a Riccati equation :

ρ 2 (t) -ρ(t) = 2a(t) ρ(0) = Ȧ(0) A(0) = 0 . (1.12)
The second equation of (1.10) gives: 2 Ȧ(t)

A(t) Ḃ(t) -B(t) = A(t)b(t), then A -2 (t) Ḃ(t) = -2A -3 (t) Ȧ(t) Ḃ(t) + A -2 (t) B(t) = -A -2 (t) 2 Ȧ(t) A(t) Ḃ(t) -B(t) = -A -2 (t) (A(t)b(t)) = - b(t) A(t) , therefore B(t) is the solution of :    Ḃ(t) = -A 2 (t) t 0 b(s) A(s) ds B(0) = 0 . (1.13) Let ∆(t) = t 0 b(s)
A(s) ds and define k by: k

(t) = k(t) θ 2 c(t) -1 2 A 2 (t)∆ 2 (t) + θ 2 2 ρ(t) k(0) = 1.
(1.14) this makes the first equation of (1.7) satisfied.

Remark 1.2. If a(t), b(t) and c(t) were constants, then we could apply the method of Rosencrans (cf. [START_REF] Rosencrans | Perturbation Algebra of an Elleptic Operator[END_REF]).

Particular cases.

i) V(t, q) = a(t)q 2 (semiclassical case).

Then

B(t) = 0 , ∆(t) = 0 , k(t) = k(0) e 1 θ 2 t 0 c(s) - 1 2 A 2 (s)∆ 2 (s) + θ 2 2 ρ(s) ds = e 1 2 t 0 Ȧ(s) A(s) ds = A(t) and η V u (t, q) = A(t) e - q 2 2θ 2 ρ(t) η u (ϕ 2 (t), A(t)q)
is the solution of (C

(V )
1 ) such that η V u (0, q) = η u (0, q) = u(q). In this case, as in the next one, the formulas really involve only A(t). ia) V(t, q) = ω 2 2 q 2 . Then

ρ(t) = -ω tanh(ωt) , A(t) = 1 cosh(ωt) , ϕ 2 (t, q) = 1 ω tanh(ωt)
and according to i)

η V u (t, q) = 1 cosh(ωt) e ω 2θ 2 tanh(ωt)q 2 η u 1 ω tanh(ωt), q cosh(ωt)
is the solution of (C (V ) 1 ). This result is due to Zambrini (cf. [START_REF] Zambrini | From the geometry of parabolic PDE to the geometry of SDE[END_REF], p.227, (51)). The existence of such a formula was mentioned in [START_REF] Miller | Symmetry and separation of variables[END_REF], p. 96, and the similar result for the Schrödinger equation with quadratic potential (with cos(ωt) taking the place of cosh(ωt)) appeared in [START_REF] Miller | Symmetry and separation of variables[END_REF], p. 83. ii) V(t, q) = λ(t)q.

Then ρ(t) = 0,

A(t) = 1, Ḃ(t) = - t 0 λ(s) ds, B(t) = t 0 Ḃ(u) du = - t 0 ( u 0 λ(s) ds) du and η V u (t, q) = k(t) e - 1 θ 2 Ḃ(t)q η u (t, q + B(t))
is the solution of (C

(V ) 1 ), such that η V u (0, q) = η u (0, q) = u(q).
iia) V(t, q) = λq.

Then A(t) = 1, ∆(t) = λt, Ḃ(t) = -λt, B(t) = -λ 2 t 2 , k(t) = e - λ 2 6θ 2 t 3 and η V u (t, q) = e - 1 θ 2 ( λ 2 6 t 3 -λtq) η u (t, q - λt 2 2 )
is the solution of (C

(V ) 1
). We hereby recover a result of Lescot and Zambrini (cf. [START_REF] Lescot | Probabilistic deformation of contact geometry, diffusion processes and their quadratures[END_REF], p.219, (1)).

isovectors

Let G V defined by :

G V = {N | L N (I) ⊂ I}
and H V by :

H V = {N ∈ G V | ∂N S ∂S = 0} .
These definitions are taken from Lescot-Zambrini (cf. [START_REF] Lescot | Probabilistic deformation of contact geometry, diffusion processes and their quadratures[END_REF]).

We remind the reader that a G V is a Lie algebra for the usual Lie bracket of vectors fields, and H V is a Lie subalgebra of G V (cf. [START_REF] Lescot | Probabilistic deformation of contact geometry, diffusion processes and their quadratures[END_REF], p.210-211).

According to Lescot-Zambrini (cf. [START_REF] Lescot | Probabilistic deformation of contact geometry, diffusion processes and their quadratures[END_REF], p.214-215, 3.9'-3.29), for all N ∈ H V , we can write :

N q : = 1 2 ṪN (t)q + N (t) N t : = T N (t)
N S : = h(q, t, S) , for σ(t), l(t) and T N (t) satisfying :

-

...

T N q 2 -¨ q + σ + T N ∂V ∂t + 1 2 ṪN q ∂V ∂q + ∂V ∂q + ṪN V - θ 2 4 TN = 0 . (2.1)
Let us set l N = l. Proposition 2.1. Let K V defined by :

K V = {N ∈ H V / N = 0}, then K V is a Lie subalgebra of H V . Proof. Let N 1 , N 2 ∈ K V ;one has [N 1 , N 2 ] q = [N 1 , N 2 ](q) = N 1 (N q 2 ) -N 2 (N q 1 ) = N q 1 ∂N q 2 ∂q + N t 1 ∂N q 2 ∂t -(N q 2 ∂N q 1 ∂q + N t 2 ∂N q 1 ∂t ) = 1 2 ṪN2 ( 1 2 ṪN1 q + N1 (t)) + T N1 ( 1 2 TN2 q + ˙ N2 ) - 1 2 ṪN1 ( 1 2 ṪN2 q + N2 (t)) -T N2 (t)( 1 2 TN1 q + ˙ N1 ) = T N1 (t) ˙ N2 -T N2 (t) ˙ N1 + 1 2 ( N1 (t) ṪN2 -N2 (t) ṪN1 ) + 1 2 (T N1 (t) TN2 (t) -T N2 (t) TN1 (t))q and [N1,N2] = T N1 (t) ˙ N2 -T N2 (t) ˙ N1 + 1 2 N1 (t) ṪN2 -N2 (t) ṪN1 = 0. Then K V is stable by the Lie bracket, therefore K V is a Lie subalgebra of H V . Lemma 2.2. If V 1 (t, q) = V 2 (t, q) + C q 2 (for C a constant), then K V1 = K V2 . Proof. Replacing V 1 (t, q) by
its value in equation (2.1), we get:

-

...

T N q 2 -¨ q + σ + T N ∂V 2 ∂t + 1 2 ṪN q ∂V 2 ∂q -2 C q 3 + ∂V 2 ∂q -2 C q 3 + ṪN V 2 + C q 2 - θ 2 4 TN = 0 ⇔ - 1 4
...

T N q 2 -¨ q + σ + T N ∂V 2 ∂t + 1 2 ṪN q ∂V 2 ∂q + ∂V 2 ∂q - 2lC q 3 + ṪN V 2 - θ 2 4 TN = 0.
Therefore, if N ∈ K V2 , as l N = 0, we see that n satisfies (2.1) for V = V 1 . Then V 1 et V 2 give the same isovectors such as N = 0 .

Therefore : ...

K V1 = K V2 . Corollary 2.3. If V 1 = C q 2 (C =
T N q 2 -¨ q + σ + 1 2 ṪN q(-2 C q 3 ) + (-2

C q 3 ) + ṪN ( C q 2 ) - θ 2 4 TN = 0, that is - 1 4
...

T N q 2 -¨ q + σ -θ 2 4 TN -2 C q 3 = 0 . As T N , and σ depend only on t, the system is equivalent to:

       2C = 0 ¨ = 0 σ = θ 2 4 TN
...

T N = 0 .

(2.2)

But C = 0, therefore N = 0 and

H V1 = K V1 = (Lemma 2.2) K V2 ⊂ H V2 .
Remark 2.4. Thus, we recover another result of Lescot and Zambrini(cf. [START_REF] Lescot | Probabilistic deformation of contact geometry, diffusion processes and their quadratures[END_REF], p.220, (3)).

  0) and V 2 = 0, then H V1 = K V1 = K V2 and therefore: H V1 ⊂ H V2 .Proof. Let N ∈ H V1 , the equation (2.1) becomes :-1 4


Meng Sun

Yide Wang

Jingjing Pan

Direction of arrival estimation by a modified orthogonal propagator method with spline interpolation

Keywords: Direction of arrival (DoA), orthogonal propagator method (OPM), propagator method (PM), spline interpolation, denoising

Direction of arrival estimation by a modified orthogonal

propagator method with spline interpolation

I. INTRODUCTION

Direction of arrival (DoA) estimation is an important issue in the field of array signal processing and has been intensively studied for decades [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF], [START_REF] Tuncer | Narrowband and wideband DOA estimation for uniform and non-uniform linear arrays[END_REF]. It is important in many practical applications such as sonar, radar, source localization, and wireless communication [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF], [START_REF] Tuncer | Narrowband and wideband DOA estimation for uniform and non-uniform linear arrays[END_REF]. A lot of methods have been developed for DoA estimation in the literature. Among these methods, subspace-based algorithms or high-resolution methods such as MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] and ESPRIT [START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF] have attracted much attention due to their asymptotic infinite resolution and unbiased estimation performance. These methods are based on the assumption of the orthogonality between the signal subspace and noise subspace of the received signals. However, the above techniques require either the eigenvalue decomposition (EVD) or singular value decomposition (SVD) of the data covariance matrix. It may be too expensive and time-consuming for real time applications with a large amount of data, such as a massive Multiple-Input Multiple-Output system with a very large number of the antennas (hundreds or thousands). Besides, this kind of methods may have some limitation due to the high computational complexity, in some practical scenarios, for example, long time field work without guaranteed electronic energy.

To solve this problem, a concept named 'propagator' was introduced in array signal processing. Many methods have been proposed based on this concept, such as propagator method (PM) [START_REF] Munier | Spatial analysis using new properties of the cross-spectral matrix[END_REF], [START_REF] Marcos | Méthodes linéaires haute résolution pour l'estimation de directions d'arrivée de sources. performances asymptotiques et complexité[END_REF] and orthogonal propagator method (OPM) [START_REF] Marcos | The propagator method for source bearing estimation[END_REF]. These methods only use linear operations to find the propagator (which can be directly extracted through the data covariance matrix), without any EVD or SVD operation. The propagator has been proven effective in DoA estimation. Besides, the propagator-based methods are more suitable for dealing with large data, which can facilitate real time application. Nevertheless, the initial propagator-based algorithm does not take the background noise into consideration, as its basis relies on the noise free received signal model. Therefore, its performance suffers from degradation in low signal to noise ratio (SNR) scenarios. The least squares method is applicable in the estimation of the propagator, which reduces the impact of noise but does not eliminate it [START_REF] Munier | Spatial analysis using new properties of the cross-spectral matrix[END_REF].

Researchers attempt to reduce the noise impact by estimating the power of the noise [START_REF] Stoica | On estimating the noise power in array processing[END_REF], [START_REF] Morelli | Joint maximum likelihood estimation of CFO, noise power, and SNR in OFDM systems[END_REF], [START_REF] Liu | Noise power estimation in SC-FDMA systems[END_REF] or applying signal enhancement techniques [START_REF] Li | Enhanced GPR signal for layered media time-delay estimation in low-SNR scenario[END_REF], [START_REF] Volaric | A data driven compressive sensing approach for time-frequency signal enhancement[END_REF]. Stoica et al. propose an eigenvalue-based method [START_REF] Stoica | On estimating the noise power in array processing[END_REF] to estimate the noise power, which requires the EVD on the data covariance matrix. The noise power is estimated by calculating the average of the eigenvalues belonging to the noise subspace. The eigenvalue method is proven to be statistically efficient [START_REF] Stoica | On estimating the noise power in array processing[END_REF]. However, it is not the focus of this paper since the major merit of the propagator-based methods is to avoid any EVD or SVD operation. In addition, Stoica et al. [START_REF] Stoica | On estimating the noise power in array processing[END_REF], Marcos and Sanchez-Araujo [START_REF] Marcos | Méthodes linéaires haute résolution pour l'estimation de directions d'arrivée de sources. performances asymptotiques et complexité[END_REF] also propose two different linear methods. In their methods, the noise power can be extracted from the structure of the data covariance matrix without EVD or SVD operation. These two linear methods are more computationally efficient, yet less statistically efficient, compared with the eigenvalue based methods [START_REF] Stoica | On estimating the noise power in array processing[END_REF].

In this paper, a new way to reduce the influence of noise without estimating the noise power is proposed from another point of view. The noise is assumed as an additive temporally and spatially white Gaussian noise. The noise hence will only affect the main diagonal elements of the data covariance matrix. Therefore, we propose to restore the noise free diagonal elements of the data covariance matrix to reduce the influence of the noise. In a column (or line), the diagonal element can be estimated by the other elements of this column (or line), which could be done by an interpolation procedure. It is worth noting that interpolation techniques have been widely used in array signal processing. In [START_REF] Tuncer | Narrowband and wideband DOA estimation for uniform and non-uniform linear arrays[END_REF], array interpolation is used to map the sample covariance matrix of a non-uniform linear array to that of a uniform linear array, so that root-MUSIC and ESPRIT can be directly applied. Instead, the interpolation procedure in this paper is used to rebuild the noise free diagonal elements of the data covariance matrix. Moreover, when the incoming signals are uncorrelated or independent, the data covariance matrix is a A c c e p t e d M a n u s c r i p t Toeplitz matrix, which means that the diagonal elements of this matrix are all the same [START_REF] Horn | Topics in matrix analysis[END_REF]. Consequently, the interpolation needs to be performed only once to recover all the noise free diagonal elements of the matrix. In this paper, the cubic spline interpolation [START_REF] Schumaker | Spline functions:basic theory[END_REF] is adopted. The performance of the proposed modified OPM algorithm is evaluated on numerical data.

The remaining of this paper is organised as follows. The narrow band received signal model in far-field is presented in Section 2. In Section 3, the proposed modified OPM is introduced. Section 4 gives the simulation results to show the performance of the proposed method compared with the conventional OPM, OPM with elimination of noise called MPEN (method of propagator with elimination of noise) [START_REF] Marcos | Méthodes linéaires haute résolution pour l'estimation de directions d'arrivée de sources. performances asymptotiques et complexité[END_REF], OPM with EVD [START_REF] Stoica | On estimating the noise power in array processing[END_REF] and MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF]. Section 5 draws the conclusion.

II. SIGNAL MODEL

Consider a uniform linear array (ULA) with N = 2M + 1 isotropic antennas impinged by K far-field narrow band independent signals. The antennas are indexed by -M, . . . , 0, . . . , M . The received signals are corrupted by a temporally and spatially white Gaussian noise. In the following, the 0th element is taken as the reference point. Then, the complex envelope of the received signal at the mth antenna can be expressed as

r m (t) = K k=1 s k (t) exp(-j2mπ d λ sin θ k ) + n m (t) (1)
where s k (t) is the complex envelope of the kth signal received at the reference point with corresponding DoA θ k (θ k ∈ [-π/2, π/2]); n m (t) is the spatially and temporally additive white Gaussian noise at the mth antenna with zero mean and variance σ 2 ; d denotes the distance between two adjacent antenna and λ is the wavelength of the received signals. d is equal to half of the wavelength of the incoming signals.

Equation ( 1) can be written in the following vector form:

r(t) = As(t) + n(t) (2)
with the following notational definitions:

• r(t) = [r -M (t) r -M +1 (t) • • • r 0 (t) • • • r M (t)]
T is the (N × 1) received signal vector or the observation vector; the superscript T denotes the transpose operation;

• A = [a(θ 1 ) a(θ 2 ) . . . a(θ K )] is the (N × K) directional matrix; • a(θ k ) = [exp(jM π sin θ k ) . . . 1 . . . exp(-jM π sin θ k )] T is the (N × 1) directional vector; • s(t) = [s 1 (t) s 2 (t) • • • s K (t)] T is the (K × 1) source vector; • n(t) = [n -M (t) n -M +1 (t) • • • n 0 (t) • • • n M (t)]
T is the (N × 1) noise vector, with zero mean and covariance matrix σ 2 I, where I is the (N × N ) identity matrix.

Assuming that the noise is independent of the source signals, the covariance matrix R can be written as:

R = E r(t)r H (t) = AR s A H + σ 2 I (3)
where E(:) denotes the ensemble average; R s is the (K × K) dimensional covariance matrix of the source vector s(t) and the superscript H denotes the complex conjugate transpose operation.

III. MODIFIED OPM FOR DOA ESTIMATION

In order to reduce the noise impact, OPM-MPEN is proposed [START_REF] Marcos | Méthodes linéaires haute résolution pour l'estimation de directions d'arrivée de sources. performances asymptotiques et complexité[END_REF]. However, this method can only localize a maximum of N/2 -1 sources. In this section, we will present a new modified OPM algorithm which avoids such limitation. The signals are supposed to be totally uncorrelated to each other and independent from the noise. The elements in R can then be expressed as follows:

R(n -m) = R(m, n) =              K k=1 Ps k + σ 2 m = n K k=1 Ps k exp(jπ(n -m) sin θ k ) m = n (4)
where m, n ∈ [1, . . . N ], P s k is the power of the kth signal. R(n -m) is a function of n -m, which represents the expression of the element of the mth row and nth column of matrix R. From (4), the noise is presented only in the diagonal elements of the covariance matrix R. However, with finite snapshots, the non-diagonal elements contain residual from noise. In order to show the noise impact on both the diagonal and non-diagonal elements of the data covariance matrix with different snapshots, we define the noise uncertainty ratios (NURs) as the ratio between the elements y mn of the noise free data covariance matrix AR s A H and the elements z mn of matrix R -AR s A H (for non-diagonal elements, NURs= { Nt Nt i=1 r(t i )r H (t i ), N t is the number of snapshots. NURs on both diagonal and non-diagonal elements of the data covariance matrix are calculated in terms of the number of snapshots with 1000 independent Monte-Carlo trials and SNR =0 dB.

It can be seen from Fig. 1, for both the diagonal and nondiagonal elements, the influence of the noise decreases (NURs increases) when the number of snapshots increases. Moreover, the influence of the noise is greater (smaller NURs) to the diagonal elements than to the non-diagonal elements, especially when the number of snapshots is large. Therefore, we focus on the diagonal elements and the influence of the noise could be reduced by rebuilding the noise free diagonal elements from the non-diagonal elements of the data covariance matrix.

The noise free diagonal elements can be estimated by the other non-diagonal elements. Define

f (x) = R(x) = K k=1 P s k exp(jπx sin θ k ) with x = n -m, therefore x ∈ [-N +1, N -1]. When x = 0, f (0) = K k=1 P s k correspond-
ing to the value of the noise free diagonal element (unknown). f (0) can be estimated from its adjacent elements of f (x)

(or off-diagonal elements R(x)), x ∈ [-N + 1, -1] [1, N -1]
by the classical spline interpolation procedure [START_REF] Schumaker | Spline functions:basic theory[END_REF]. In addition, as R is a Toeplitz matrix, the diagonal elements of this matrix are all equal. Therefore, it is enough to estimate only one noise free diagonal element. As the variation of the norm of non-diagonal elements is slower than that of the elements themselves, in this paper, the spline interpolation is applied on the norm of the non-diagonal elements in the first column and first line of the data covariance matrix for the calculation of the corresponding noise free diagonal element f (0). It means by using |f

(-N + 1)|, |f (-N + 2)|, . . . , |f (-1)|, |f (1)
| . . . , |f (N -1)|, we can estimate the noise free diagonal element f (0) through spline interpolation. The detail is given as follows: to estimate the expression of |f (x)|, the spline function S(x) is built with 2N -2 sample points (-N + 1, . . . , -1, 1, . . . , N -1) and 2N -3 intervals. S i (x) is the spline function inside the ith interval. S(x) should obey the following constrains:

• In each interval, S(x) = S i (x) is a three order polynomial. • S(x) and its derivatives S (x), S (x) should be continuous for x ∈ (-N + 1, N -1).

• S(x i ) = |f (x i )|, x i = -N + 1, . . . , -1, 1, . . . , N -1. • S (-N +1) = S (-N +2) and S (N -1) = S (N - 2).
The principle of spline interpolation is to take a shape that minimizes the bending (under the constraint of passing through all sample points). Both S (x) and S (x) should be continuous for x ∈ (-N + 1, N -1). The above can be achieved with polynomials of degree 3 or higher. The classical spline interpolation is to use polynomials of degree 3 [START_REF] Schumaker | Spline functions:basic theory[END_REF]. Therefore, S i (x) can be expressed as S i (x) = a i x 3 + b i x 2 + c i x + g i , with a i , b i , c i , g i unknown coefficients. By using the principle of the spline interpolation, the 8N -12 unknown coefficients inside S(x) can be calculated. Then, the spline function S(x) can be rebuilt. Consequently, the noise free diagonal elements can be estimated by S(0), which is used to restore all the diagonal elements of R. The elements of the denoised data covariance matrix R nf can then be written as

R nf (m, n) = K k=1 P s k exp(jπ(n -m) sin θ k ). (5)
Since the noise is removed, OPM can be directly applied on the denoised data covariance matrix for DoA estimation. OPM [START_REF] Marcos | The propagator method for source bearing estimation[END_REF] partitions matrix A into two sub-matrices as follows:

A = A 1 A 2 (6)
where A 1 and A 2 are (K×K) and ((N -K)×K) dimensional matrices (N > K), respectively. Since θ k ∈ [-π/2, π/2] (for any θ i , θ j ∈ [-π/2, π/2], if θ i = θ j , sin θ i = sin θ j ), A is of full rank and A 1 is non-singular. Under the above condition, the following (K × (N -K)) dimensional propagator P can be defined

P H A 1 = A 2 (7)
or

[P H , -I N -K ]A = Q H A = 0 ( 8
)
where

I N -K is a ((N -K) × (N -K)) dimensional identity matrix.
Equation [START_REF] Stoica | On estimating the noise power in array processing[END_REF] shows that the directional vectors of sources are orthogonal to the columns of Q. Therefore, the DoA of each source can be estimated by searching the peak positions of the following OPM pseudo-spectrum:

P OPM (θ) = 1 a H (θ)Q(Q H Q) -1 Q H a(θ) . (9)
However, propagator P is unknown, which should be retrieved from data. Firstly, the data covariance matrix is partitioned into the following expression:

R nf = [G, H] (10)
where G and H are the (N × K) and (N × (N -K)) dimensional matrices, respectively. The relationship between propagator P and these two matrices (G and H) is given by

H = GP. (11)
Therefore, propagator P can be directly extracted through [START_REF] Li | Enhanced GPR signal for layered media time-delay estimation in low-SNR scenario[END_REF] by the least squares solution:

P = G + H ( 12
)
where the superscript + denotes the Moore-Penrose inverse operator.

Finally, the general steps of the proposed modified OPM are summarized as follows:

1) Estimate the covariance matrix R (3).

2) Make the estimated covariance matrix R Toeplitz, denoted by R T . The Hermitian Toeplitz matrix R T can be built by calculating the first line of this matrix,

r t (m) = 1 N -m N -m n=1 R(n, m+n), with m = 1, . . . , N -1.
3) Apply the spline interpolation on the norm of the 4M non-diagonal elements (the norm of the elements) of the first column and the first line of R T to estimate its diagonal element. 4) Construct the new data covariance matrix R nf by replacing the diagonal elements of R T with the above estimated diagonal element. 5) Apply the OPM on the new data covariance matrix R nf to estimate the DoA.

A c c e p t e d M a n u s c r i p t

IV. SIMULATIONS AND DISCUSSION

The performance of the proposed modified OPM is evaluated on numerical data with 4 simulations. Considering a ULA with 11 isotropic elements and 2 far-field narrow band uncorrelated signals with equal power, that is, N = 11, M = 5 and K = 2.

In the first simulation, the pseudo-spectrums of the conventional OPM, OPM-MPEN and proposed method are compared. The following two cases are considered:

• Case 1. θ 1 = -6 • and θ 2 = 7 • . • Case 2. θ 1 = -1 • and θ 2 = 3 • .
The data covariance matrix is estimated with 400 independent snapshots and SNR= 0 dB. The pseudo-spectrum search is performed over [-15 • , 15 • ] with step size 0.001 • . The angles corresponding to the two highest peak positions in the pseudospectrum allow estimating the DoAs of signals. Figs. 23show the pseudo-spectrums of the conventional OPM, OPM-MPEN and proposed method for the two different cases. For both cases, the proposed method is able to detect the true DoAs of signals. As shown in Figs. 23, two peaks corresponding to the DoAs of signals are clearly detected by the proposed method. For the conventional OPM and OPM-MPEN, they can detect the true DoAs of signals in large angle separation situation (in Case 1) and their pseudospectrums are similar to that of the proposed method. However, when the angle separation becomes small, like in Case 2, the conventional OPM fails in DoA estimation and the OPM-MPEN estimates the DoA with bias. In the second simulation, the performance of the proposed method is tested with respect to SNR. A Monte-Carlo process is constructed with 1000 independent runs of the method. The performance is assessed in terms of the RMSE of the estimated DoA defined as

RMSE = 1 KJ K k=1 J j=1 θkj -θ k 2 ( 13
)
where θkj denotes the estimated DoA of the kth incoming signal for the jth run; θ k denotes the true DoA of the kth incoming signal. The angles in Case 2 are used here. SNR ∈ [-6, 10] dB. The proposed method is compared with OPM-MPEN, OPM-EVD and MUSIC; moreover, the Cramér-Rao bound (CRB) results are also provided [START_REF] Stoica | MUSIC, maximum likelihood, and cramér-rao bound[END_REF], [START_REF] Stoica | The stochastic crb for array processing: A textbook derivation[END_REF]. Fig. 4 plots the RMSE of DoA estimation against SNR with the number of snapshots equals to 400. As expected, the RMSEs of the above mentioned methods continuously decrease when SNR increases. According to the principle of OPM [START_REF] Marcos | The propagator method for source bearing estimation[END_REF], the propagator based methods are efficient in terms of computational burdens, but their statistical performance is worse than that of the subspace based methods. Therefore, the gap (estimation performance) between the CRB and propagator based methods can be large. Same results can be found in Fig. 4 that the RMSE of MUSIC versus SNR is smaller than that of the propagator based methods at the beginning of SNR. However, with higher SNR (-2 dB), these methods tend to work similarly. In addition, MUSIC applies EVD operation with computational load O(N 3 ), which is greater than that of OPM (O(N 2 K + (N -K) 2 K), DoA spectrum search is not taken into account in computational complexity).

Compared with the other denoising procedures, the proposed method has a more significant decrease of RMSE than that of the OPM-MPEN; its performance is similar to that of OPM-EVD. The computational loads of the denosing procedures of OPM-MPEN, OPM-EVD, and the proposed method are shown in Table I. From Table I, the proposed method has the smallest computational complexity among the 3 compared methods.

Method

Computational loads OPM In the third simulation, the statistical performance of the modified OPM versus the angle separation between two incoming signals is assessed with a Monte-Carlo process of 1000 independent runs of the method. One of the incoming signals is fixed at θ 1 = 0 • while the other comes from θ 2 = θ 1 + ∆θ with same power. SNR = 0 dB, ∆θ varies from 0.5 • to 8 • . The number of snapshots is 400. From Fig. 5, it is clear that the mentioned methods fail to detect the DoAs when the angle separation is small (0.5 • to 2 • ). The performance of the proposed methods is similar to that of the second simulation. Moreover, the RMSEs of the mentioned methods continuously decrease as the angle separation increases. When the angle separation becomes greater, these methods tend to have similar performance. Moreover, the performance of the proposed modified OPM versus the number of snapshots is evaluated with a numerical example of 1000 independent Monte-Carlo trials. Case 2 is studied with SNR = 0 dB. Fig. 6 plots the RMSEs of the modified OPM, OPM-MPEN, OPM-EVD and MUSIC as function of the number of snapshots. In Fig. 6, the above mentioned methods fail when the number of snapshots is small (less than 10). With the increasing number of snapshots, the RMSEs decrease. Different from MUSIC, we can observe that the propagator based methods require good estimation of the data covariance matrix with a sufficient number of snapshots (> 50).

-MPEN O((N -K)K 2 + K 3 /3) OPM-EVD O(N 3 ) Modified OPM O(N + log N )

V. CONCLUSION

In this paper, we propose a modified OPM method for DoA estimation. By using the spline interpolation to rebuild the noise free diagonal elements of the data covariance matrix, the proposed method offers a new way to reduce the noise impact. The performance of the proposed method is tested on numerical data. The proposed modified OPM shows good performance in DoA estimation.

  mn |}; for diagonal elements, NURs= { N m=1 |y mm |}/{ N m=1 |z mm |}), where R = 1

Fig. 1 .

 1 Fig.1. NURs on both the diagonal and non-diagonal elements of the data covariance matrix versus the number of snapshots with 1000 independent Monte-Carlo trials, SNR of the received signal equals to 0 dB.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Case 1, pseudo-spectrums of the conventional OPM, OPM-MPEN and proposed modified OPM for DoA estimation. The two DoAs are θ 1 = -6 • and θ 2 = 7 • .

Fig. 4 .

 4 Fig. 4. Case 2, RMSE of estimation versus SNR with 400 snapshots.

Fig. 5 .

 5 Fig. 5. RMSE of DoA estimation as function of angle separation.

Fig. 6 .

 6 Fig. 6. RMSE of DoA estimation as function of the number of snapshots, SNR = 0 dB.

ACKNOWLEDGEMENT

The authors would like to thank the Shanghai Sailing Program (No.19YF1419100) and the National Natural Science Foundation of China (No.61701297) for funding part of this work.