Floquet spectrum for anisotropic and tilted Dirac materials under linearly polarized light at all field intensities
J. C. Sandoval-Santana V. G. Ibarra-Sierra A. Kunold Gerardo G. Naumis
cond-mat.mes-hall






AIP/123-QED

Sample title]Floquet spectrum for anisotropic  and tilted Dirac materials
under linearly polarized light at all field intensities




jcarlosss@fisica.unam.mx


vickkun@fisica.unam.mx


akb@azc.uam.mx


naumis@fisica.unam.mx




             


The Floquet spectrum in an anisotropic tilted Dirac semimetal 
modulated by linearly polarized light is addressed through the
solution of the time-dependent Schrödinger equation for the
two-dimensional Dirac Hamiltonian via the Floquet theorem.
The time-dependent wave functions and the quasienergy spectrum
of the two-dimensional Dirac Hamiltonian under the normal incidence
of linearly polarized waves are obtained for an arbitrarily intense
electromagnetic radiation. We applied a set of unitary transformations
to reduce the Schrödinger equation to an ordinary second-order
differential Hill equation with complex coefficients.
Through the stability analysis of this differential equation,
the weak and strong field regimes are
clearly distinguished in the quasi-spectrum. In the
weak electric field regime, above a certain threshold given by the field parameters, the spectrum mostly resembles
that of free electrons in graphene.
Below this threshold, in the strong electric field regime,
the spectrum abruptly becomes highly anisotropic and a gap opens up.
As an example,  we apply the results to the particular case of borophene.



[
    Gerardo G. Naumis
    December 30, 2023
=====================





§ INTRODUCTION


The superior physical, mechanical an chemical properties
of two-dimensional (2D) materials makes them
an ideal playground to study new and exciting kinds
of quantum phases
<cit.>.
Their remarkable electronic and optical features
have intensely driven the development
of novel and innovative
optoelectronic devices <cit.>
as broad band optical modulators
<cit.>,
solar cells
<cit.>,
infrared photodetectors
<cit.>
and hybrid plasmonic devices
<cit.>.

In the last decade persistent efforts have been made
to harness the unique features of graphene's so-called
dressed electrons
<cit.>
to design different kinds of optoelectronic devices.
Electromagnetic dressing, attained when
electrons strongly couple to electromagnetic fields,
substantially renormalize the energy and the velocities
in graphene<cit.>.
In turn, renormalized parameters highly
depend on the light polarization: in graphene circularly polarized fields
open a dynamical gap in the Dirac point while linearly polarized fields
live it intact. Particularly, dressed electrons under linearly polarized
light induce an anisotropy of the electron dispersion relation <cit.>.
Electromagnetic dressing could, therefore, be used to tune the electronic
and optical
properties of graphene, including band gap and carrier
velocities<cit.> which are clearly manifest in
several measurable physical properties as the photocurrent
<cit.>.
A similar approach has been adopted in ultrafast material science
promising optical and
mechanical control of the physical properties of 
2D materials<cit.>.
In this case, it is time-dependent
strain that acts as  pseudo-electromagnetic field
<cit.>.

In graphene, the weak-field regime, where light-matter coupling
is perturbative<cit.>,
is well understood.
It may be pictured as a quantized-photon field interacting
with massless Dirac fermions having a conic dispersion relation.
However, in the intense-field regime, this perturbative expansion cannot
converge. The Dirac dispersion relation is highly distorted
as a result of the electronic dressing rendering the quantized-photon
picture invalid.
It is clear then that this emerging field requires new
theoretical tools that go beyond the conventional techniques
<cit.>.
There are numerous proposals to approach
this problem
<cit.>,
but, even for graphene, some aspects of the interaction
between carriers and strong time-driven fields remain elusive.

Another essential question that remains to be fully answered
is how such an intense electromagnetic radiation would affect
other Dirac materials<cit.>
as borophene or black-phosphorus
<cit.>,
for example.
Particularly borophene is a remarkable anisotropic material
<cit.>.
After being theoretically
predicted thirty years ago<cit.>,
it was until 2015 that it was synthetized  <cit.>.
Borophene turns out to be stronger and even more flexible than graphene.
It is a good conductor of both electricity and heat and it
is expected to be a superconductor <cit.>
with relatively high transition temperatures.

In a series of recent papers, we have investigated time-driven
anisotropic Dirac Hamiltonians as the one that describes borophene 
<cit.>. 
However, the mathematical complexity of the problem required several
approximations that are only valid for very intense 
fields<cit.>.
Therefore, the critical link between
weak fields, treatable with perturbation theory,
and strong fields is still missing.
In particular, the cases of linearly and elliptically polarized light were
studied in the intense field regime <cit.>,
where the rather convoluted Hill equation with time dependent coefficients
was simplified into the very well known Mathieu equation with constant
coefficients.
Nevertheless, the complete study of the Hill equation entails
the interplay between the strong and weak field intensity regimes,
key to comprehend the formation of the quasienergy spectrum.
This effect can be thoroughly studied by analyzing
the quasi-energy spectrum using the Floquet theory
<cit.>.

In this paper we address the general problem of a particle that
obeys the  anisotropic Dirac Hamiltonian subject to linearly
polarized light with an arbitrarily large field intensity.
The calculation of the quasienergy spectrum and the wave function
is achieved through a set of unitary transformations
that enables to reduce the matrix differential equation
into a scalar differential equation.
This, in turn, is readily solved via the Floquet theorem
and a Fourier spectral decomposition of the periodical part
of the solution.
It is shown that in the intense field regime the quasienergy
spectrum abruptly develops an anisotropic structure, absent
in the weak field regime.

The paper is organized as follows.
In Sec. <ref> we introduce the low-energy effective 
two-dimensional anisotropic Dirac Hamiltonian and we study
the case of borophene subject to an arbitrarily intense
linearly polarized electromagnetic field.
Subsequently, in Sec. <ref>,
we analyze the quasienergy spectrum that emerges from the
Hill equation by means of the Floquet approach.
Also, we find the time-dependent wave functions.
Finally, we summarize and conclude in Sec. <ref>.





§ TWO-DIMENSIONAL ELECTRONS IN A TILTED DIRAC CONE
SUBJECT TO ELECTROMAGNETIC FIELDS





 §.§ The anisotropic Dirac Hamiltonian


We start by considering a low-energy anisotropic Dirac Hamiltonian
close to one of the Dirac points.
In the particular case of 8-Pmmn borophene  it is given by 
<cit.>

    Ĥ= ħ v_t k_yσ̂_0
    +ħ[v_x k_xσ̂_x +v_y k_yσ̂_y],

where k_x and k_y are the components of the two-dimensional
momentum vector k,  σ̂_x
and σ̂_y are the Pauli matrices,
and σ̂_0 is the 2× 2 identity matrix.
The three velocities in the anisotropic 8-Pmmn borophene Dirac
Hamiltonian (<ref>)
are given by v_x=0.86v_F, v_y=0.69v_F and v_t=0.32v_F where
v_F=10^6 m/s <cit.>
is the Fermi velocity.
In Eq. (<ref>), the
last two terms give rise to the familiar form of the
kinetic energy leading to the Dirac cone and
the first one tilts the
Dirac cone in the k_y direction.
These two features are contained in the
energy
dispersion relation<cit.>

    E_η,k=(v_t/v_y) k̃_y +νϵ,

where

    ϵ=√(k̃_x^2+k̃_y^2),

and ν= ± 1 is the band index.
In Eq. (<ref>), we used the set of
renormalized moments
k̃_x=ħ v_x k_x, k̃_y=ħ v_y k_y.
The corresponding free electron wave function is,

    ψ_ν(k)
      = 1/√(2)[
      [            1; νexp(i θ_k ) ]]

where θ_k= tan^-1 (k̃_y/k̃_x).




 §.§ Linearly polarized waves
 

Now we consider a charge carrier, described by the two-dimensional
anisotropic Dirac Hamiltonian, subject to an electromagnetic wave that 
propagates along a direction perpendicular to the surface of the crystal.
The effects of the electromagnetic field
are introduced in the Hamiltonian
(<ref>) through the Peierls 
substitution<cit.> 
ħk→ħk-eA
where A=(A_x,A_y) is the vector potential
of the electromagnetic wave.
A considerable simplification
can be achieved by adopting
a gauge in which A only
depends on time. 
The Hamiltonian (<ref>)
is thus transformed into

    Ĥ= v_t/v_y(k̃_y-e v_y A_y)σ̂_0
    
     +(k̃_x-e v_x A_x)σ̂_x 
     +(k̃_y-e v_y A_y)σ̂_y.

Assuming that the electromagnetic wave is
linearly polarized, the vector potential
can be written as

    A=E_0/Ωcos(Ω t)r̂,

where r̂=(1,0) is the polarization vector,
E_0 is the uniform amplitude of the electric field
and Ω is the angular frequency of the electromagnetic wave.
Observe that here the field A is not quantized and
is treated clasically. Thus, our results are valid for a field
with a large number of photons, which can be represented
by a quantum coherent field. 
In the Schrödinger equation corresponding to
(<ref>),

    i ħd/d tΨ(t)
      =ĤΨ(t),

the two dimensional spinor can be expressed as
Ψ(t)
 =(Ψ_A(t),Ψ_B(t))^⊤, where
A and B label the two sublattices.

The main difficulty in
deducing the wave function's explicit form
resides in that the Hamiltonian
(<ref>)
couples the differential equations for
the Ψ_A(t) and Ψ_B(t)
spinor components due to
the terms that are proportional
to σ̂_x and σ̂_y.
To uncouple the spinor components we proceed
as follows.
First, applying a 45^∘ rotation around the k_y axis
of the form

    Ψ(t)
    =exp[-i/ħ(π/4)σ̂_y]Φ(t),

conveniently transforms the non-diagonal σ̂_x
matrix into σ̂_z.
Indeed, substituting (<ref>) into 
Eq. (<ref>) yields

    i d/dϕΦ(ϕ)
    =2/ħΩ[(v_t/v_y) 
    k̃_yσ̂_0+ Π̃_xσ̂_z
    + k̃_yσ̂_y ]Φ(ϕ)   ,

where the only non-diagonal remaining term
is the one proportional to σ̂_y.
In the foregoing equation, ϕ= Ω t/2,
Π̃_x=k̃_x-ζ_x cos(2 ϕ) and
ζ_x=ev_xE_x/Ω.
The spinor components of
Φ(ϕ)=(Φ_+(ϕ),Φ_-(ϕ))^⊤
are given by
Φ_+(ϕ)=[Ψ_A(ϕ)+Ψ_B(ϕ)]/√(2) and 
Φ_-(ϕ)=[Ψ_A(ϕ)-Ψ_B(ϕ)]/√(2).
Second, the term proportional to σ̂_0
in  Eq. (<ref>)
is removed by adding a time-dependent
phase to the wave function

    Φ(ϕ)=exp[-2i(v_t/v_y)k̃_y/ħΩ ϕ σ̂_0]χ(ϕ),

where χ(ϕ)=(χ_+1(ϕ),χ_-1(ϕ))^⊤. 
Finally, after inserting Eq. (<ref>) into Eq.
(<ref>), we follow the
procedure shown in Appendix <ref>.
The resulting diferential equation
takes on the form
of the Hill equation <cit.>

    χ”(ϕ)+F(ϕ)χ(ϕ)=0   ,

where the matrix F(ϕ) is defined as

    F(ϕ)=[a+q_1cos(2ϕ)
    +q_2cos(4ϕ)]σ̂_0
    +q_3sin(2ϕ)σ̂_z.

The Hill equation parameters

    a   =   (2/ħΩ)^2(ϵ^2+ζ^2_x/2),
    
    q_1   =   -8(k̃_x/ħΩ)(ζ_x/ħΩ),
    
    q_2   =   2(ζ_x/ħΩ)^2,
    
    q_3   =    4 i(ζ_x/ħΩ).

are expressed in terms of the ratios of the characteristic energies
of the system.
Thereby, ϵ/ħΩ
is the ratio of the electron kinetic energy to the photon energy,
ζ_x/ħΩ  is the ratio of
the work done on the charged carries by the electromagnetic wave
to the photon energy and k̃_x/ħΩ
is the ratio of the x part of the electron kinetic energy
to the photon energy.

Expressing (<ref>) as
a second order differential equation
is quite advantageous for the calculations that follow.
First, the evolution operator that propagates 
the state χ in time must be diagonal
since F(ϕ) is solely composed of the diagonal
matrices σ̂_0 and σ̂_z.
As a consequence of this,
the scalar differential
equations for the χ_+1(ϕ) and χ_-1(ϕ)
spinor components decouple. Moreover, the differential
equation for the χ_-1(ϕ) component turns out to be
the complex conjugate of the one for χ_+1(ϕ).
Both differential equations may be summarized by

    χ_η”(ϕ)+[a+ q_1cos(2ϕ)+q_2cos(4ϕ).
    
    +. η q_3 sin(2ϕ)]χ_η(ϕ)=0,

where η=± 1.

This alternative form of the Schrödinger equation considerably
simplifies the computation and analysis of the stability spectrum. 








§ QUASIENERGY SPECTRUM AND WAVE 
FUNCTION





 §.§ Quasienergy spectrum


In this section, we analyse the quasienergy spectrum
produced by the Hill equation (<ref>),
and its relations with the stability of its solutions.
The determination of the stability regions of this differential equation 
is quite challenging mainly due to the
imaginary coefficient q_3.
The real coefficients a, q_1 and q_2 alone give rise
to the Whittaker-Hill equation <cit.>, widely discussed
in the spectral theory of periodic differential equations.
The imaginary coefficient q_3, however,
introduces additional difficulties that
are rarely addressed in the literature
<cit.>.
Among other things, it yields complex characteristic values.
Despite the added complexity, Eq. (<ref>)
may be approached by Whitakker's original assumption
<cit.> that the solution
should take the Floquet normal form

    χ_η(ϕ)=exp[i μ_ηϕ] u_η(ϕ),

due to the periodicity of the Hamiltonian. Also,
according to the Floquet theorem,
u_η(ϕ) must be a function with period π.
The function
μ_η(a) is termed the characteristic exponent.
As detailed in the following sections, the
characteristic exponent and the quasienergy
spectrum are closely connected.
Thanks to the periodicity of the u_η(ϕ) function,
χ_μ(ϕ) might be expressed as
a Fourier series expansion. In this manner,
the second-order differential
equation with time-variable periodic coefficients 
is traded for
a matrix eigenvalue problem. The eigenvalues
that stem from it are
precisely the characteristic exponents μ_η(a).
These have the form

    μ_η(a)
    =1/πcos^-1[1+Δ_η(0)
    ( cos(√(a)π) -1 )],

where a and Δ_η(0)  are given by Eqs. 
(<ref>) and (<ref>),
respectively. For a detailed calculation of
μ_η(a) refer to Appendix <ref>.

In general, the solutions of the Hill equation 
(<ref>)
fall either on stable (Im[μ_η(a)] =0)
or unstable regions (Im[μ_η(a)]  0)
depending on the values taken by the coefficients a, q_1,
q_2 and q_3
<cit.>. 
Surprisingly, the solutions mostly
fall on the stable regions when these coefficients are
restricted by the parametrization of
Eqs. (<ref>)-(<ref>)
for seemingly arbitrary domain spaces of
k̃_x, k̃_y and ζ_x.
Fig. <ref>
shows a density plot the characteristic exponent 
spectrum for μ_+1(a) as a function of
the normalized momenta k̃_x/ħΩ
and k̃_y/ħΩ. Over the entire parameter
space covered in this figure the solutions
of the Hill equation are stable.

One of most striking features of the spectrum
emerges when one compares
the quasienergies at low and high electric field amplitudes.
The contrast between the effects of
weak (E_x= 10^-4V/m, ζ_x/ħΩ=10^-3)
and strong  (E_x=9.7V/m, ζ_x/ħΩ=5)
electric field amplitudes is shown in
Figs. <ref> (a)
and (b) respectively.
In the weak electric field regime
(see Fig. <ref> (a))
the isolines of the characteristic exponent μ_+1(a)
as a function of the normalized momenta
form a pattern of concentric circles,
In the strong electric field regime
(see Fig. <ref> (b))
a spectrum of vertical lines, perpendicular
to the electric field direction, emerges
close to the Dirac point (k̃_x=k̃_y=0).
The vertical grill is surrounded by an elliptical
outline that is approximately delimited by the
white circle of radius ϵ=ζ_x.
In the vicinity of this contour (ϵ>ζ_x)
the spectrum takes on
an elliptical form but further out (ϵ≫ζ_x)
the spectrum recovers the circular shape
observed in the low electric field regime.

The appearance of the vertical grill in the
quasienergy spectrum is due to the interplay between
the terms cos(2ϕ), sin(2ϕ) and cos(4ϕ)
in Eq. (<ref>).
For example, ϵ<ζ_x implies that
q_2>| q_3| and | q_3|>q_1,
therefore the most significant term is cos (4ϕ).
By contrast, ϵ>ζ_x implies that
q_1>| q_3| and | q_3|>q_1,
and the dominant element is  cos (2ϕ).  

Further details of the spectrum are appreciated in 
Fig. <ref> where we show a 3D plot
of the quasienergy μ_+1(a) as a function
of k̃_x and k̃_y in the high electric field
regime (ζ_x/ħΩ=5).





 §.§ Wave function and Floquet spectrum

As mentioned above,
the wave function in (<ref>) must satisfy the
Floquet theorem as a result of the time periodicity of
the Hamiltonian (<ref>).
Hence, following the Floquet theorem, Eqs. (<ref>),
(<ref>)
and (<ref>) can be combined into the wave function

    Ψ(ϕ)
     =𝒩exp[-i/ħ(π/4)σ̂_y] 𝒰(ϕ)χ(0),

where 𝒩 is a normalization constant, 
χ(0)=(χ_+1(0),χ_-1(0))^⊤ is the
initial state vector in
(<ref>) and 𝒰(ϕ) denotes the time
evolution operator
such that χ(ϕ)=𝒰(ϕ)χ(0).
As we pointed out before, because F(ϕ) is a diagonal matrix,
the evolution operator must also be
diagonal. Thus, it can be expressed
quite generally as

    𝒰(ϕ)=
    [ exp[-i ε_+1ϕ]u_+1(ϕ)                    0;                    0 exp[-i ε_-1ϕ]u_-1(ϕ) ],

where u_-1(ϕ)=u_+1^*(ϕ) have period π and

    ε_η= 2v_t/v_yk̃_y/ħΩ
    + μ_η(a).

It is easy to verify that
the wave function (<ref>) reduces to the
free-particle wave function (<ref>)
when the electric field vanishes.
From the Floquet theory
<cit.>,
the time evolution operator is periodic
𝒰(ϕ)=𝒰(ϕ+lπ)
and quasienergy can be expressed as

    ℰ_η,l=ε_η+l
     =2v_t/v_yk̃_y/ħΩ+ μ_η,l(a),

where

    μ_η,l(a)=μ_η(a)+l,     l∈Z,

is the characteristic exponent
for different Brillouin zones which are tagged by
the integer subscript l.

In Fig. <ref>, we show the quasienergy for 
ℰ_+1,0.
At the bottom of this figure, a density plot of μ_+1,0(a)
is shown for reference. One notes that the spectrum consists of
the quasienergies μ_+1(a), plotted in Fig.
<ref>, on to the tilt
that comes from the first term of Eq. (<ref>).
It arises from the anisotropic character of the
Hamiltonian (<ref>).
Even though it strongly distorts the symmetry
of the Dirac cone, it has been shown that
interband transitions are not affected by it
in the zero-temperature limit <cit.>.

Different cross sections of the quasienergy spectrum μ_η,l(a) 
are shown in Fig. <ref>
for the Brillouin zones l=-1,0,1.
Panels (a) and (b) show the spectrum section planes
k̃_y=0 and k̃_x=0, respectively,
in the  weak electric field regime
( ζ_x/ħΩ= 10^-3, E_x=10^-3V/m,
and Ω=50×10^9Hz) .
In this case, both k̃_y=0 and k̃_x=0
cross sections are almost identical since
the electric field is not intense enough to provoke
any distortion to the free particle spectrum.
The quasienergies for the many Brillouin zones
as functions of k̃_x [panel (a)]
or k̃_y [panel (b)] have the form of a triangular function.
In the strong electric field regime
(ζ_x/ħΩ=5, E_x=9.7V/m
and Ω=50×10^9Hz)
both cuts are radically different.
In panel (c) we note that the k̃_y=0 cross section
of the spectrum resembles those corresponding
to the weak electric field regime also exhibiting
a triangular shaped function of k̃_x.
However, in the k̃_x=0 plane [panel (d)]
the presence of the electromagnetic field
becomes evident as the spectrum is warped
approximately in the domain -2 < k̃_y < 2.
At both ends of this range the quasienergy abruptly
recovers the triangular feature that characterizes
the spectrum in the absence of electromagnetic radiation.
This distortion is roughly bounded by the two black
vertical lines that correspond to the black
dots in Fig. <ref> (b),
where the strength of the electromagnetic field is comparable
to the energy of the unperturbed system, i.e. ϵ=ζ_x.

A quite robust feature of the spectrum is
the preservation of the gapless Dirac cone
in the vicinity of the Dirac point:
the quasienergies μ_+1,l and μ_-1,l as
well as ℰ_+1,l and ℰ_-1,l touch
at the tip of the Dirac point despite the 
intensity of the linearly polarized electromagnetic field.
Nevertheless, the conic dispersion relation is stretched
along the k_y direction as a result of the
renormalization of the v_y component of the velocity
due to electronic dressing.
Despite the absolute absence of a gap
a the tip of the Dirac cone,
it is possible  to open up gaps in other
zones of the spectrum.
In Fig. <ref> (d),
the appearance of a small gap
between
μ_+1,l and μ_-1,l can be appreciated 
at k̃_y/ħΩ≈± 2.








§ CONCLUSIONS


We investigated the quasienergy spectrum of
an anisotropic tilted Dirac material subject to an arbitrarily
intense linearly polarized electromagnetic field.
To this end,  we studied the behavior of 8-Pmmn borophene under
the normal incidence of a linearly polarized field.
We worked out the time-dependent wave function
and the quasienergy spectrum from the Schöendiger equation
via the Floquet theory.
The quasienergy spectrum exhibits a sharp difference between the
weak and strong electromagnetic field regimes. While in the first
the quasienergy as a function of the quasimomenta
(k̃_x and k̃_y)
is highly isotropic the latter presents an anisotropic
pattern in the low energy region that resembles a grid
aligned perpendicularly to the direction of the
radiation's electric field.
This pattern abruptly disappears beyond the
threshold where the free kinetic energy
of the carriers is larger than the energy associated
to the electric field.
Near this threshold, a gap opens up.
Probably the most astonishing feature of this
spectrum is that, eventhough the gapless
Dirac cone is preserved, the dispersion relation
is stretched along the direction perpendicular
to the field's polarization. This is an outcome
of the electronic dressing and the consequent 
rescaling of the velocity.
This mechanism could be exploited to tune
the electronic properties of Dirac materials
through the field parameters.
A fundamental aspect, yet to be addressed, is
the capability of circularly polarized fields
to adjust these properties.




§ AUTHOR'S CONTRIBUTIONS


All authors contributed equally to this work.



§ ACKNOWLEDGEMENTS


This work was supported by DCB UAM-A grant numbers
2232214 and 2232215, and UNAM DGAPA PAPIIT IN102620.
 J.C.S.S. and V.G.I.S acknowledge the total support from 
DGAPA-UNAM fellowship. 



§ DATA AVAILABILITY


The data that support the findings of this study are available from the corresponding author
upon reasonable request





§ 


This Appendix covers the method developed by  Whittaker
<cit.>
to determine the characteristic exponent
(<ref>).
Eq. (<ref>) is the starting point.
This equation is periodic
and therefore its solution must comply with the
Floquet theorem. Thus, the solution is given by

    χ_η(ϕ)=e^iμ_ηϕu_η(ϕ)   ,

where η=± 1, u_η(ϕ) is a function with period π
and μ_η denotes the characteristic exponent.
Thanks to the periodicity of u_η(ϕ), the
wave function can
be expanded as a Fourier series as

    χ_η(ϕ)=e^iμ_ηϕ∑_r=-∞^∞C^(η)_2re^i 2r ϕ  .

Substituting the preceding equation into (<ref>)
and rearranging  the coefficients,
we obtain the following recurrence relation

    γ_2rC^(η)_2(r-2)+α_2rC^(η)_2(r-1)+C^(η)_2r
    
    +β_2rC^(η)_2(r+1)+γ_2rC^(η)_2(r+2)=0,

where 

    α_2r=1/2q_1-iη q_3 /a-(μ+2r)^2,


    β_2r=1/2q_1+iη q_3 /a-(μ+2r)^2,


    γ_2r=1/2q_2/a-(μ+2r)^2.

The equation parameters a, q_1, q_2 and q_3 are defined in  Eqs.
(<ref>)-(<ref>).
The recurrence relation (<ref>) can be put in the form
of a linear equation as

    𝒜_r(μ,η,a,q_1,q_2,q_3)C^(η)=0

where
C^(η)
=(C^(η)_2,C^(η)_4,C^(η)_6,...)^⊤
and

    𝒜_r(0,η,a,q_1,q_2,q_3)=
          [       1    α_2r    γ_2r       0       0       0       0       0       0;  β_2r-2       .       .       .       .       .       .       .       0;       0       .       1     α_4     γ_4       0       0       .       0;       0       .     β_2       1     α_2     γ_2       0       .       0;       0       .     γ_0     β_0       1     α_0     γ_0       .       0;       0       .       0    γ_-2    β_-2       1    α_-2       .       0;       0       .       0       0    γ_-4    β_-4       1       .       0;       0       .       .       .       .       .       .       . α_-2r+2;       0       0       0       0       0       0   γ_-2r   β_-2r       1 ].

To avoid the trivial solution we demand that the determinant of
the precedent matrix vanishes:

    Δ_η(0)
    =det[𝒜_r(0,η,a,q_1,q_2,q_3)]
    =0.

It can be proven that this determinant may be written in the
compact form<cit.>

    sin^2(μ_η(a) π/2)
    =Δ_η(0)sin^2(√(a)π/2)  .

Solving the above equation for μ_η(a)
we obtain

    μ_η(a)= 1/πcos^-1[1+Δ_η(0)(cos(√(a)π)-1) ] .

This expression presents the advantage
of being efficiently evaluated by
using numerical methods. Strictly speaking
𝒜_r
is infinite-dimensional, however good numerical
convergence of μ_η(a) is achieved by
cutting it down to a 400× 400 matrix.



 







