Pseudo-scalar {\boldm $ q\bar q $} bound states at finite temperatures within a Dyson-Schwinger--Bethe-Salpeter approach
S. M. Dorkin L. P. Kaptari B. Kämpfer

   










 

-10mm
=170mm
=237mm



















P
^2
^2
⟨
⟩
  k
 q 
 p 
P

    Bogoliubov Lab. Theor. Phys., 141980, JINR, Dubna,  RussiaInternational University Dubna, Dubna, Russia Bogoliubov Lab. Theor. Phys., 141980, JINR, Dubna, RussiaHelmholtz-Zentrum Dresden-Rossendorf, PF 510119, 01314 Dresden, GermanyHelmholtz-Zentrum Dresden-Rossendorf, PF 510119, 01314 Dresden, GermanyInstitut für Theoretische Physik, TU Dresden, 01062 Dresden, Germany
     The combined Dyson-Schwinger–Bethe-Salpeter equations are employed at non-zero temperature. The truncations refer to a rainbow-ladder approximation augmented with an interaction kernel which facilitates a special temperature dependence. At low temperatures, T → 0, we recover a quark propagator from the Dyson-Schwinger (gap) equation which delivers, e.g. mass functions B, quark renormalization wave function A,
    and two-quark condensate q q̅ smoothly interpolating to the T = 0 results, despite the broken O(4) symmetry in the heat bath and discrete Matsubara frequencies. Besides the Matsubara frequency difference entering the interaction kernel, often a Debye screening mass term is introduced when extending the T = 0 kernel to  non-zero temperatures. At larger temperatures, however, we are forced to drop this Debye mass  in the infra-red part of the longitudinal interaction kernel to keep the melting of the two-quark condensate in a range consistent with lattice QCD results. Utilizing that quark propagator for the first few hundred fermion Matsubara frequencies we evaluate the Bethe-Salpeter vertex function in the pseudo-scalar q q̅ channel for the lowest boson Matsubara frequencies and find a competition of  q q̅ bound states and quasi-free two-quark
    states at T =  O (100 MeV). This indication of pseudo-scalar meson dissociation below the anticipated QCD deconfinement temperature calls for an improvement of the approach, which is based on an interaction adjusted to the meson spectrum at T = 0.
    Pseudo-scalar qq̅ bound states at finite temperatures within a Dyson-Schwinger–Bethe-Salpeter approach 
         B. Kämpfer
        
    =======================================================================================================
    
    
    
    § INTRODUCTION
    
     The description of mesons as quark-antiquark bound states within the framework of the Bethe-Salpeter (BS) equation
     with momentum dependent quark mass functions, determined by the Dyson-Schwinger (DS equation, is able
     to explain successfully many spectroscopic data, such as
     meson masses <cit.>,
     electromagnetic properties
     of pseudoscalar mesons and their radial excitations <cit.>
      and  other   observables <cit.>.
     Contrary to purely  phenomenological models, like  the quark bag  model,
     such a formalism  maintains   important features of QCD, such as dynamical chiral
     symmetry breaking, dynamical quark dressing, requirements of the renormalization group theory etc., cf. Ref. <cit.>.
     The main ingredients here are  the full quark-gluon vertex function and  the dressed gluon
     propagator,  which are entirely determined  by  the running coupling
     and the bare quark mass parameters. In principle, if one were able to solve the Dyson-Schwinger equations
     for the quark propagator, the gluon propagator, the quark-gluon and inter-gluon vertices ,
     supplemented by ghosts if necessary, then the approach would  not depend on  any additional parameters.
    In practice,  one restricts oneself to calculations within effective models which
      specify the dressed vertex function Γ_μ (for quark-gluon coupling) and interaction kernel D_μν
      (for the gluon propagator).
      The rainbow-ladder approximation <cit.> is a model with rainbow truncation of the
      vertex function Γ_μ→γ_μ in the quark DS equation  and
      a specification of the dressed quark-quark interaction kernel (for gluon exchanges)
       as g^2 D_μν(k) → G(k^2) D_μν^free(k) within the Landau gauge. (Here, γ_μ
       is a Dirac gamma matrix and D_μν stands for the gluon
       propagator; g is the coupling strength and k denotes a momentum.)
    
       The model is completely specified once a form
    is chosen for the “effective coupling” G(k^2). The ultraviolet
    behavior is chosen to be that of the QCD running coupling α_s(k^2)=g^2(k^2)/4π computed within one-loop
    approximation; the ladder-rainbow truncation then generates the correct perturbative
    QCD structure of the DS and BS
    equations. Moreover, the  ladder-rainbow truncation preserves such an  important feature of the
    theory as the maintenance of the Nambu-Goldstone theorem in the chiral limit,
    according to which   the spontaneous chiral symmetry breaking results in an appearance of a
     (otherwise absent) scalar term in the quark propagator of the DS  equation. As a consequence, in the BS
     equation a massless pseudoscalar bound state should appear. By using the Ward  identities, it has been proven (see,
     e.g.  Refs. <cit.>) that, in the chiral limit, the DS equation for the quark propagator and
     the BS equation for a massless pseudo-scalar in ladder approximation are completely  equivalent.
     It  implies  that such a  massless bound state (pion) can be interpreted as a Goldstone boson.
     This results in a straightforward understanding    of the pion as  a Goldstone boson and  quark-antiquark
     bound state at the same time.
    
        Another important property of the DS and BS equations is that, at zero temperatures, or equivalently - in vacuum -
        they are explicitly Poincaré  invariant.    This  frame-independency of the approach provides
        a useful tool in studying processes where  a rest frame for mesons cannot or  needs not be defined.
        A rather different situation is encountered at finite temperatures,   where a  particle is embedded in
        a heat bath. In this case,  the propagator   of the particle is defined by  two four-vectors,
        the four-momentum p  of the particle and a unit four-vector u
        characterizing the heat bath itself <cit.>.
        The scalar parts of the propagator then depend on p^2 and (up). An existence of the
        statistical thermal factor  exp(-  H/T), H being the   QCD hamiltonian, implies the choice of the unit vector
        u  as u=( 0,1),   i.e. the heat bath at rest <cit.>.
     As a consequence, the O(4) symmetry of equations is lost and only the spatial translation invariance
         is maintained. This requires a separate treatment of the transversal and longitudinal, with respect to the heat bath,
       parts of the kernel with the need of additional functions in parametrizing  the quark propagators and the   BS vertex functions.
    
        In quantum field theory, a system at finite temperature can be
     described within the  imaginary-time formalism,
    which uses the Matsubara frequancies <cit.>. Due to finiteness of the heat bath
     temperature T and the requirements of periodicity (bosons) or antiperiodicity (fermions) of the
     considered  fields <cit.>,
      the Fourier transform  to  Euclidean momentum space becomes discrete in the time direction, resulting
     in a discrete spectrum of the energy, known as the Matsubara frequencies. Consequently, the
     interaction kernel and the DS and BS solutions must become also discrete with respect to these frequencies.
    It is worth mentioning that there is no straight forward connection of the obtained solutions
    in the discrete Euclidean space to  the corresponding continuum quantities in  Minkowski space <cit.>.
    In particular, within the Matsubara formalism one can define the mass of the bound state
    in different ways in dependence of the choice of the total momentum P=(,iΩ_N);
    at N=0 one defines the so-called screening mass, which has no direct relation with the usual
    inertial mass. At || 0 and N 0 the corresponding pole mass P^2=-M^2 is also discrete
    and hitherto it is note clearly established  how to relate it to the inertial mass. In principle,
    at each particular choice  of P
    one can define a quark correlator and try to define the inertial mass by computing the
    momentum Fourier transforms and the resulting spatial integral  of the correlator (analogue of the 2-points  Schwinger
    functions, see e.g. <cit.>) and   by taking the time derivative of the logarithm
     of the obtained expression, see also <cit.>.
     Evidently, such a procedure requires knowledge of the pole mass for all, or at least for a large enough number
     of Matsubara frequencies. It is worth mentioning that a  growing interest nowadays is in approaches based on
     real-time calculations of quark and gluon propagators and their spectral functions within the functional
     renormalization group <cit.>. These approaches also require analytical continuation of
     the propagators from Euclidean to Minkowski spaces.
    
    
     In the present paper, we calculate the pole masses of pseudo-scalar
     quark-antiquark bound-states only for the few
     first Matsubara frequencies and investigate their dependence on temperature in a large interval, including
     the "dissociation region". We anticipate that the main peculiarities
     of  masses computed  at some particular Matsubara frequencies will
     reflect the general T dependence of the inertial masses at high temperatures.
    We treat the bound states within the BS formalism within the
    same approach as the one used in solving the DS equation, i.e. with the rainbow truncation and
    AWW interaction kernel <cit.>.
    Recall that the merit of the approach is that, once the effective parameters are fixed
     (usually the effective parameters of the
     kernel are chosen, cf. Refs. <cit.>,  to reproduce the known
     "experimental" data from lattice calculations, such as the quark mass function and/or  quark condensate,
     the whole spectrum of known mesons is supposed to  be described,   on the same footing,
     including also excited states. The achieved  amazingly good description  of
     the mass spectrum with only  few effective parameters encourages one to employ the same approximations
     to the  DS and BS equations also at finite temperatures with the hope that, once an adequate
     description of the quark propagators at non-zero temperature    is accomplished, the corresponding solution
     can be implemented into the BS equation for mesons  to investigate the meson properties
     in hot matter.
    
       At low temperatures, the properties of hadrons   in nuclear matter
       are expected to change in comparison with the vacuum ones, however, the main  quantum
       numbers, such as  spin and orbital momenta, space and inner parities etc. are maintained.
       The hot environment may modify the hadron masses, life times  (decay constants) etc. Contrarily,
       at sufficiently large temperature in hot and dense strongly interacting matter,
         (phase) transitions may occur, related to quark deconfinement phenomena, as e.g.
       dissociation of hadrons into quark-gluon degrees of freedom. Therefore,
       this temperature region is of great interest,
       both from a theoretical and experimental point of view. So far,
       the truncated DS and BS formalism has  been mostly used
        at large   temperatures  to investigate the critical phenomena
        near and above the  pseudo-critical and  (phase) transition values predicted by lattice simulation data
       (cf. Refs. <cit.>).
      It has been found that, in order to achieve an agreement of   the  model results with lattice data,
      a modification of the vacuum interaction kernel is required.
      First of all it has been observed that modifications  of the
      (pure phenomenological,  nonperturbative) infra-red (IR) term by including the T dependence of the
      perturbative Debye mass results in an essential suppression of the longitudinal interaction kernel
      at large temperatures <cit.>.
      As a consequence, the predicted critical temperatures T_c of chiral symmetry restoration
      in the chiral limit and also the critical temperatures for  the maximum of chiral susceptibility, as well as of
       inflection point of the mass function or of the quark condensate,
       are much smaller, by about 50%, than that expected from the lattice calculations.
       It became evident that  modifications of the IR term at large temperatures are necessary.
       Namely, the IR  term has to vanish abruptly in this region and to be replaced by another
       phenomenological kernel.  For instance, it has been  suggested <cit.>  to employ
        a kernel with a  Heavyside step-like behaviour  in the vicinity
       of the (pseudo-) critical temperature  T_c. Then,  it becomes possible to  achieve a
      rather reliable description of such quantities as the quark spectral function,
      plasmino modes, thermal masses etc., see also Ref. <cit.>. However, a use   of  such a
       discontinuously  modified interaction in
      the BS equation in the whole temperature range becomes hindered.
    
      Another strategy of solving the DS equation in a
      larger interval of temperatures  is to  utilize directly   the available  results of nonperturbative lattice QCD (LQCD)
      calculations  to fit, point by point, the interaction kernel
      at given  temperatures. In such a way one achieves a good description
      of the quark mass function and condensate  for  different   temperatures,  including the region
      beyond  T_c <cit.>. The success of such approaches demonstrates that the
      rainbow approximation to the DS equation with a proper choice of the interaction kernel  is quite
      adequate in  understanding the properties  of quarks in  a hot environment. Nevertheless, for
       systematic studies of quarks
      and hadrons within the BS equation, on needs a  smooth parametrization of the
       kernel  in the whole interval of the considered
       temperatures. In view of still scarce LQCD data, such a direct parametrization from "experimental" data is problematic.
       An alternative method  is to solve simultaneously   a (truncated) set of Dyson-Schwinger equations
      for the quark and gluon propagators within some additional approximations <cit.>.
      This approach also provides
      good description of the quark mass function and condensate   in vicinity of T_c,
      however it becomes too cumbersome in attempts
      to solve the BS equation,  since in this case one should solve a too large system of equations.
    It should also be  noted that  there are other  investigations of the quark propagator
    within the rainbow truncated   DS equation,
    which   employ solely the vacuum parameters in  calculations of  T dependencies of quarks <cit.>
       without further attempts to accommodate    to LQCD results in the kernel. As a result    one
      finds that  the  critical behaviour of the propagators (e.g.  chiral symmetry restoration)  starts at temperatures
       much smaller than the ones  expected from  LQCD.
    
     In the present paper we are interested in a qualitative study of the masses of quark-antiquark bound states at moderate and
     large temperatures. Particular attention is paid to the region where the bound state mass becomes larger
     than the sum of two masses of free quarks. This temperature range probably indicates a dissociation of the
     ground bound state and existence of only highly excited correlations.
     For such a qualitative analysis we employ the simplest
     version of the kernel which consists of only the infra-red term of the Maris-Tandy <cit.> kernel, known also as the
     Alkofer-Watson-Weigel <cit.>  kernel, referred  to as the  AWW model.
     For further kernels cf. <cit.>,
     in particular to symmetry preserving issues and going beyond the rainbow-ladder approximation.
      As mentioned above, including of the Debye mass in this term results in  too low values of the
      critical temperatures, and a proper modification  is needed. For this sake, in the present
      paper we ignore  the Debye mass which allows to obtain
      solutions of the DS equation with characteristics close to that of lattice calculations.
       However, for completeness we also
      present some results of calculations with accounting of the Debye mass too.
      We start with the AWW interaction kernel, known at T=0  to provide  a good description  of
     the mass spectrum of scalar, pseudo-scalar, vector and tensor mesons.
    Then we extend the DS and BS equations  to finite temperatures  and solve the corresponding equations
    in a large interval of T and Matsubara frequencies. Results of calculations are compared with masses of free quarks.
    
    
     Our paper is organized as follows.
     In Sec. <ref>,  we recall  the truncated Dyson-Schwinger
     (tDS) and truncated Bethe-Salpeter (tBS) equations in vacuum and at finite temperatures.
      The rainbow approximation for the DS equation kernel in vacuum is  introduced
      and  the  system of equations for the quark propagator is solved.
    
     Numerical solutions for the quark propagator at finite bare quark masses
     are discussed in Sec. <ref>, where
     the inflection points of the quark condensate  and the mass function are considered as a
     definition of the pseudo-critical
     temperature T_c. It is shown that, for finite quark masses,
     the inflection method determines the pseudo critical temperatures by ∼ 50% smaller than the ones
     obtained by other approaches, e.g. by LQCD calculations.
      It is found that, by merely ignoring  the Debye mass in the infra-red term, it becomes possible to
     reconcile the model and lattice QCD results; the such obtained values of
     T_c∼ 135 MeV is  only by about  10-15% smaller then the lattice results (reported, e.g. in
     Refs. <cit.>).
     In Sec. <ref> we generalize the tBS equation to  finite temperature by exploring   the same procedure  as
     that used in  generalizing the tDS equation to finite T. As already mentioned, due to the breaking of the O(4) symmetry,
     the tBS equation cannot be longer considered for mesons at rest. The total momentum P^2 = -M^2_qq is replaced  as
     P( 0,iM_qq) → P( P,iE), where now the energy E becomes discrete and depends on the Matsubara
     frequency Ω_N = 2π N T.  In dependence of the chosen components of the total momentum P, one can
     define at least three different ground state  masses:
     (i) the thermal mass at P=0, M^2= (2π N)^2T^2,
     (ii) the imaginary screening mass at zero Matsubara frequency n=0, P^2 = -M^2 and
     (iii) the pole mass M^2 = (2π N)^2T^2 -  P^2. In principle, these masses are
     not obligatorily the same at given temperature T, as discussed in Refs. <cit.>, for instance.
     In the present paper, we consider the pole masses which at
      given temperature and Matsubara frequency are continuous functions of the three momentum P^2.
    
    
       Summary and conclusions are collected in Sec. <ref>.
    
    
    § BASIC FORMULAE
    
    
     §.§ Truncated Dyson-Schwinger  and Bethe-Salpeter    equations in vacuum
    
    To determine the bound-state mass of a quark-antiquark pair one needs to solve the
    DS and the homogeneous   BS equations, which
    in the rainbow ladder approximation    and in Euclidean space  read
    
        S^-1(p)= S_0^-1(p) + 4/3∫d^4 k /(2π)^4[g^2 D_μν(p-k) ]γ_μ S(k)
        γ_ν  ,
        
              Γ(P,p) =   -4/3∫d^4k/(2π)^4γ_μ S(η_1 P+k) Γ(P,k) S(-η_2 P+k))γ_ν [g^2     D_μν(p-k) ]   ,
    
    where η_1 and η_2=1-η_1 are the partitioning parameters  defining the quark momenta
    as p_1,2=k ±η_1,2 P with P, p and k denoting the total and relative momenta
    of the bound system, respectively;[Usually,  for quarks of masses m_1,2 the partitioning
    parameters  are chosen as  η_1,2=m_1,2/(m_1+m_2). However, in general, the tBS solution
    is independent of the choice of η_1,2.    ]Γ(P,k) stands for  the tBS vertex function being a 4× 4 matrix,
    S_0(p)= (iγ· p +m)^-1  and
    S(p)= ( iγ· p A(p)  +B(p))^-1 are the
       propagators of bare and dressed  quarks, respectively, with mass parameter m and
        the dressing functions  A(p) and B(p).  At zero temperature the above equations are O(4) invariant and
        the propagator functions A(p) and B(p)  depend solely on p^2= p^2+p_4^2. The
        total momentum P = ( 0, iM_qq̅), for a particle at rest,  is   an external parameter
        for (<ref>); the momenta of individual quarks are
        p_1,2^2 = -M^2_qq̅/4 +k^2 ± iM_qq k_4, where the Euclidean components
        k and k_4 of the relative momentum are defined as
         k^2 = ^2 + k_4^2 which, for the BS equation,  are the integration variables. In order
         to reduce the  four-dimensional integral  to a two-dimensional one, usually
         the phase space volume is parametrized  by the Euclidean modulus k,
         a hyperangle  χ_k defined  as | k| = ksinχ_k, k_4=kcosχ_k and by the
         standard angular space volume dΩ_ k. Then the tBS vertex functions
         Γ(P,k) are decomposed over the  hyperspherical harmonics basis which allow to carry out  the spacial
         angular integration  analytically (see for details Refs. <cit.>).
          The resulting tBS equation represents a system of two-dimensional integral
         equations, where the vectors p_1,2 vary in  a restricted complex domain located within
         parabolas p_1,2^2 = -M^2_qq̅/4 +k^2 ± iM_qq̅ kcosχ_k.
        In  Euclidean space  the Dirac matrices    γ_4=γ_0,_E=-i_M are  Hermitian and  obey
        the anti-commutation relation {γ_μ,γ_ν}=2δ_μ,ν; for
        the four-product one has (a· b)= a b +a_4b_4 .
              The masses M_qq̅  of mesons as bound states  of a m_1-quark and m_2-antiquark follow
       from the  solution of the tBS equation at P^2=-M_qq̅^2  in specific J^PC channels,
       with the solution of the tDS equation (<ref>) as input into the calculations in
       Eq. (<ref>). The interaction between  quarks in the  pair is encoded in  g^2 D_μν,
        imagined     as gluon exchange. For consistency, the same interaction is to be employed
         in the tDS equation (<ref>) for the inverse dressed quark propagator.
    
    Often, the coupled equations of the quark propagator S, the gluon propagator D_μν and the
    quark-gluon vertex function Γ_μ (not to be confused with the BS vertex function
    Γ(P,p)), all with full dressing (and, if needed, supplemented by ghosts and their
    respective vertices), are  considered as an
    integral formulation being equivalent to QCD.
      In practice, due to  numerical problems, the finding of the  exact solution of the
    system of coupled equations for  S-D_μν-Γ_μ can hardly be accomplished  and
    therefore some approximations <cit.> are appropriate.
     Being interested in dealing with mesons as quark-antiquark bound states
     of the BS equation (<ref>), one has to provide the quark
      propagator which depends on the gluon propagator and vertex as well, which in turn
      depend on the quark propagator. Leaving  a detailed discussion of the variety of approaches in dressing
     of the gluon propagator and vertex function in DS equations
     (see e.g. Refs. <cit.> and references therein
     quoted) we mention only that in solving the DS equation for the quark propagator
     one usually employs truncations of the  exact interactions and replaces the gluon propagator combined with the vertex
     function by an effective interaction kernel [g^2 D_μν].
     This leads to the tDS  equation  for
     the quark propagator which may be referred   to as the gap equation.
    In explorative  calculations, the choice of the form of the effective interaction is inspired by
    results from calculations of Feynman diagrams within  pQCD maintaining
    requirements of symmetry and asymptotic behaviour  already implemented,
    cf. Refs. <cit.>.
     The  results of such calculations, even in the simple case of accounting only for one-loop
     diagrams with proper regularization
     and renormalization procedures,  are rather cumbersome for further use in
     numerical calculations, e.g. in the framework of BS or Faddeev equations.
     Consequently, for practical purposes,  the wanted exact results are
     replaced by  suitable parametrizations of the vertex and the gluon propagator.
     Often, one employs  an  approximation which corresponds to one-loop calculations of diagrams with
     the  full  vertex function Γ_ν,
      substituted by the free one, Γ_ν(p,k)→γ_ν (we suppress the color structure and
      account cumulatively for the strong coupling later on).
       To emphasize the replacement of combined gluon propagator and vertex (in the Landau gauge) we use the notation
    [g^2D_μν], where an additional power of g from the second undressed vertex is
    included.
    
    
    
     §.§ Choosing an interaction kernel
    
    
     Note  that the nonperturbative behaviour of the kernel [g^2 D_μν]
     at small momenta, i.e. in the infra-red (IR) region, nowadays is
     not uniquely determined and, consequently, suitable models are needed.
    There are several ansätze in the literature for the IR  kernel, which can be
    formally classified in the two groups: (i) the IR part is parametrized  by  two
    terms - a delta distribution at zero momenta and an exponential, i.e. Gaussian term,
    and (ii) only the Gaussian term is considered.
     In principle, the IR term must be supplemented
    by a ultraviolet (UV) one, which assures the correct asymptotics at large momenta.
    A detailed investigation <cit.> of
     the interplay of these two terms has shown that, for bound states,
     the IR part  is dominant for light (u,  d and s)  quarks
     with a decreasing role for heavier (c and b) quark  masses  for which the
     UV part may be quite important in forming  mesons with masses M_qq̅> 3-4 GeV as bound states.
     In the vacuum, if one is interested in an analysis of light mesons  with M_qq̅≤ 3-4 GeV, the UV term can be omitted.
    
    
    Following examples in the literature <cit.>
    the interaction kernel in the rainbow approximation in the Landau gauge is chosen
    as
    
        g^2(k^2)  D_μν (k^2) =(
         D_IR(k^2)+ D_UV(k^2)) ( δ_μν-k_μ k_ν/k^2), 
              
        D_IR(k^2)=
                4π^2 D k^2/ω^6 e^-k^2/ω^2,
           
         D_UV(k^2) =
                 8π^2 γ_m F(k^2)/ln[τ+(1+k^2/Λ_QCD^2)^2] ,
    
     where the first term originates from  the effective IR part of the interaction
     determined by soft, non-perturbative effects, while the second one ensures the correct
    UV asymptotic behaviour of the QCD running coupling, F(k^2) = { 1 - exp(-k^2/[4m_t^2])}/k^2 with m_t
    as an adjustable parameter, m_t ≃ 0.5 GeV, and  τ = e^2 -1,
    Λ_QCD =0.234 GeV, and γ_m = 12/(33 - 2N_f ) for N_f as active flavors.
     In what follows we restrict ourselves to  the simplest version of the model, namely with the interaction kernel
    where the UV term (the effect of which is found <cit.> to be negligibly small)
    is ignored at all. As said above, this interaction
     is the  Alkofer-Watson-Weigel <cit.> kernel. A prominent
    feature of such an interaction is that, with  only a few adjustable parameters -  D, ω  in the IR term and
     quark mass parameter m in the bare quark propagator S_0 - it
    provides a  good description of the pseudoscalar, vector and tensor meson
     mass spectra <cit.> at zero temperature. It should be noted that, despite
    the negligible contribution, the inclusion  of the UV term   in the numerical calculations causes
    additional uncertainties, related to the divergency of the integrals and with requirements of regularization and
     subtraction procedures <cit.>.
     Therefore, at finite temperatures a tempting choice of the interaction  is to keep the AWW kernel the same as in vacuum.
      Notice also that, even at zero temperatures, the tBS equation
    becomes rather complicate for numerical solutions since it involves the quark propagator functions
    in the complex Euclidean space, where they can (actually they do) exhibit pole-like singularities. A rather detailed analysis
    of solving the tBS equation in vacuum in presence of poles has been reported in Ref. <cit.>.
    Reiterate that, within the rainbow approximation, the Euclidean P_4 = iM_qq̅ is an external parameter
     in the tBS equation, while k_4 is an integration variable.
    
     
    
     §.§ Finite temperatures
    
     The theoretical treatment of  systems at non-zero temperatures differs from the
     case of zero temperatures.   In this case,  a preferred  frame is determined
     by the local rest system of the  thermal bath. This means that the O(4) symmetry
     is broken and,  consequently, the dependence of the quark propagator on p and p_4  requires a separate treatment.
     To describe the propagator in this case a third function C is needed,
     besides the functions A and B introduced above for vacuum.
      Yet, the theoretical formulation of the field theory
    at finite temperatures can be performed in at least two, quite different, frameworks which
    treat  fields  either with ordinary time variable
    t    ( - ∞  < t < ∞), e.g.  the termo-field dynamics
    (cf. <cit.>) and path-integral formalism (cf. <cit.>),
    or with imaginary time it= τ (0< τ< 1/T) which is known as the Matsubara
    formalism <cit.>.
      In this paper we utilize the imaginary-time formalism within which the
    partition function is defined and  all calculations are performed in Euclidean space.
     Since, due to the  Kubo–-Martin-–Schwinger condition <cit.>
     for equilibrium systems, the (imaginary) time evolution is restricted to the interval [0… 1/T],
     the quark  fields become   anti-periodic  in time with the period 1/T.
     In such a case, the Fourier transform is not longer continuous and the energies
     p_4 of particles become discrete <cit.> which are known as the Matsubara frequencies, i.e.
      p_4=ω_n=π T (2n+1) for Fermions (n∈ℤ).
     The inverse quark propagator is now parametrized as
    
        S^-1( p,ω_n)=i A(^2,ω_n^2)+iγ_4 p_4 C(^2,ω_n^2)+
        B(^2,ω_n^2).
    
    Accordingly, the interaction kernel  is decomposed in to a transversal and longitudinal part
    
        [g^2 D_μν (,Ω_mn)] =
        ^T_μν D^T(,Ω_mn,0)+
        ^L_μν D^L(,Ω_mn,m_g),
    
    where Ω_mn=ω_m-ω_n
    and the gluon Debye screening    mass m_g is introduced in the longitudinal part of the
    propagator, where k^2=^2+Ω_mn^2+m_g^2 enters. The Debye mass
    describes  perturbatively the screening of chromoelectric fields at large temperatures, therefore
    it is relevant for the perturbative UV term in the limit of quark-gluon plasma where the light-quark bound states
    do not longer exists, and the system ie to be described by quark and gluon degrees of greedom.
     As for the nonperturbative  pure phenomenological  IR term,  it is not a priory clear
     whether the Debye mass  has to be implemented  in the IR term or not.
     We consider both possibilities, with and without the Debye mass in the AWW model.
    
    The scalar coefficients D^L,T are defined below.
    The  projection operators _μν^L,T can be written as
    
        ^T_μν={[              0,         μ, ν=4,; δ_αβ-k_α k_β^2;  μ,ν=α,β=1,2,3, ].
                                                 
               ^L_μν =δ_μν-k_μ k_νk^2- ^T_μν.
    
    By using the Feynman rules for finite temperatures, see e.g. Ref. <cit.>, it can be shown that
    the gap equation has  the same form as in case of T=0, Eq. (<ref>), except that within the Matsubara formalism
    the integration over k_4 is replaced by the summation over the corresponding  frequencies, formally
    
        ∫d^4 k(2π)^4⟶ T∑_n=-∞^∞∫d^3 k(2π)^3.
    
    The explicit form of the system of equations for A,B and C to be solved can be found in
    Refs. <cit.>.
    
    The form of the  interaction kernel  is taken  the same as at T=0,  except that at finite T
    the Debye mass can modify the definitions of the four momenta. The information on these
     kernels is even more sparse than in the case of T=0. While the effective parameters
     of the kernel in vacuum can be adjusted to some known experimental data,
     e.g. the meson mass spectrum from the BS equation, at finite temperature one
     can rely  on results of  QCD calculations, e.g. by using results of the nonperturbative lattice calculations.
     There are some indications, cf. <cit.>, that at low temperatures the gluon propagator is insensitive
      to the temperature impact, and the interaction can be chosen  as at T=0 with D^T=D^L <cit.>.
     However, in a hot and/or dense medium  the gluon is also subject to medium
     effects and thereby acquires modes with with
    finite transversal (known also as the Meissner mass) and longitudinal (Debye or electric)
     masses. Generally, these masses  appear as independent   parameters with
     contributions depending on the considered  process <cit.>.  In practice, the perturbative expressions are employed. The
    role of the Meissner masses   in the tDS equation at zero chemical potential  is not yet well
    established and requires separate investigations. This is beyond the goal of the present paper where
    only  the role of Debye mass in  the solution of the tDS equation, m_g, is analysed
    for the IR term. In most approaches based on the tDS equation within the rainbow approximation,
      it is also  common practice to ignore the effects of Meissner masses.
     This is inspired by the results  of a tDS equation analysis in the high temperature and density
      region <cit.> which indicate that the Meissner mass  is of no importance in tDS equation.
       At this level, the Matsubara frequencies and the Debye mass are the only T depending  part of the kernel.
     The Debye mass is well defined in the weak-coupling regime. In <cit.>
     it was found that in the leading order
    
        m_g^2=α_s
        π3[ 2 N_c+N_f]T^2,
    
    where N_c=3 and N_f denote the number of active color and flavor degrees of freedom, respectively;
      the running coupling  α_s in the one-loop approximation  is
    
        α_s(E)≡g^2(E)4π=f(E) 12π11 N_c-2N_f
    
    with E being the energy scale. For the temperature range considered in the present
    paper we adopt f(E)→ 2, which is an often employed
     choice for the Debye mass in the tDS equation <cit.>.
      For N_f=4   equation (<ref>) results in
     the  Debye mass m_g^2 = 16π^2T^2/5, which is commonly used in
      literature <cit.>. It should be noted, however,  that such a choice
     of f(E) is not unique. It may vary in some interval, in   dependence on
      the employed method of  infrared regularization <cit.>.
      Since the Debye mass enters  as an additional energy parameter in
      k^2=^2+Ω_mn^2+m_g^2, which determines the Gaussian form of the longitudinal part of the
       AWW kernel (<ref>),  an increase of   m_g^2 results in a shift of the tDS solution towards lower temperatures
       leaving, at the same time, the shape  of the solution practically  unchanged.
    
    The transversal and longitudinal parts of the interaction kernel (<ref>)
    can be  cast in the form
    
        D^T(,Ω_mn,0)=D_IR(^2+Ω_mn^2), 
              
        D^L(,Ω_mn,m_g)=D_IR(^2+Ω_mn^2+m_g^2).
    
    In the present paper we use the parameters of the AWW model for the interaction kernel
    
      ω=0.5 GeV, D=1 GeV^2, m_l=m_u,d=5 MeV, m_s=  115 MeV.
     
      Recall that the AWW model  provides  values for the vacuum quark condensate in a narrow corridor,
     - qq̅_0 = [(0.244 - 0.251)  GeV]^3, and
     the correct mass spectra for pseudo-scalar, vector and tensor  mesons
      as quark-antiquark bound states <cit.>.
    
    
    
    
    
    
    § SOLUTION OF THE TDS EQUATION
    
     As mentioned above, the AWW model retains  only the IR term in the interaction kernel.
     Thus, it allows to avoid
     the logarithmic divergency  of the integrals originating from the  UV term and additional
     uncertainties related to different schemes adopted  for regularization and subtraction
     procedures,  see discussions in Refs. <cit.>.
    Despite the fact that,   without the UV term within the AWW model all the relevant integrals in tDS and tBS equations are convergent,  further calculations involving the solution A,B,C of tDS
      can encounter  divergences (not directly connected with the interaction kernel) which have to be  properly
      regularized and  subtracted. Such a situation occurs
      in calculations of the quark condensate which, by definition,  is quadratically divergent.
    For instance, at T=0 and m_q 0,
    
        qq̅_0 =
        -1/8π^2∫  p^3  B(p^2) / p^2 A^2( p^2)+B^2( p^2) d p,
     where p^2=^2 +p_4^2. At large  p^2
    the asymptotic solution of the tDS equation becomes  A( p^2)→ 1, B( p^2)→ m_q
    and the  integral (<ref>) is manifestly quadratically divergent. The same situation occurs at finite temperatures.
    Usually the quark condensate is regularized by
    subtracting at large momenta the asymptotic quark mass m_q and defining the regularized  (subtracted)
    condensate as (see,  Ref. <cit.>)
    
        qq̅_0^l-m_l/m_h qq̅_0^h =  qq̅_ren.^l-m_l/m_h qq̅_ren.^h,
    
     where m_l and m_h denote the  mass of  light, e.g. u,d, and heavy, e.g. s, quarks, respectively.
      Exactly the same procedure is applied to determine  the quark condensate at finite T,
      see also Ref. <cit.>.
      The remaining multiplicative divergences can be removed by normalizing to
       quark condensate  at zero temperature.
    
    We solve numerically the tDS equation for the functions A, B and C  by an iteration
    procedure. The summation over ω_n is truncated  at a sufficiently
    large value of n=N_max, where in our calculations N_max∼ 200 for
     low temperatures and N_max∼ 64 for temperatures T > 80 -100 MeV are utilized.
     The integration over the momentum || is replaced by a  Gaussian integral sum,
     with the Gaussian mesh of N_G=96 points. In Fig. <ref>
    we present the solutions for A and B at the lowest Matsubara drequency n=0 as functions
    of the three-momentum modulus | k| with and without taking into account the Debye mass.
    
    
      From this figure one infers that  the effect of the Debye mass is rather large at low and moderate values of the
      momentum | k|. Since for the AWW kernel (<ref>) the Debye mass enters as a
      Gaussian exponential factor in the longitudinal part of the integral,
       i.e. ∫ d||...∼exp(-m_g^2/ω^2) ∫ d||...,
       the effect increases with increasing  temperature, resulting in a considerable  suppression of the
      solution at large T. To analyse the role of the Debye mass we solved the tDS equation also in a
      large interval of T
      with and without the Debye mass. Notice that, since the solution of the tDS is not an observable,
     one does not have experimental quantities to be  compared with.
       Instead, one can try to reconcile   the model calculations with results of "exact" LQCD
      data. Suitable quantities of the tDS equation which can be
      related to corresponding quantities from LQCD are, e.g. the dependencies on temperature
      of the mass function B and the quark condensate ⟨ qq̅⟩.
    
      Results of such calculations are presented in Fig. <ref>, where we exhibit
      the mass function B, which can be considered as mass parameter according to the
      decomposition (<ref>), at zeroth momentum ||=0 and zero Matsubara frequency ω_0 =π T
      and the regularized condensate ⟨ qq̅⟩ (<ref>, <ref>)
       normalized at T=0.
      
    It is seen from the figure that the solutions of the mass function  and the quark condensate, calculated
    with and without the Debye mass taken into account,
     are essentially different at large temperatures. As expected,
     including the Debye mass shifts the results to lower temperatures.
    To compare the model results with LQCD results
     we use the method of the maximum of the chiral susceptibility, i.e. the maxima
    of the derivatives of  B and/or  qq̅ with respect to the quark bare mass, as well as
    the inflection point of the mass function or of the condensate, i.e. the maxima of the
    corresponding derivatives with respect to the temperature <cit.>:
    
        χ_B(T)=d^2 B(0,ω_0)dT^2;   χ_qq̅(T)=d^2  qq̅dT^2.
    
    The (pseudo-) critical temperature
    T_c is fixed by the condition  . χ_B(T) |_T=T_c=0 and/or .χ_qq̅(T) |_T=T_c =0, see Fig. <ref>.
    
     
    
    
    
    Figures <ref> and <ref> clearly demonstrate that the inflection points at finite quark bare masses
    provide much smaller (pseudo-) critical temperatures T_c if the Debye mass is taken into
    account.  For the AWW model,  one has T_c∼ O(100) MeV, cf. also Ref. <cit.>.
    Without  the Debye mass in  the IR term,  the positions of the inflection points occur
      at T_c∼ O(130-135) MeV, for both, the solution B and  the regularized quark condensate. This
      pseudo-critical temperature  is quite close to that obtained in LQCD calculations <cit.>
      which report T_c∼ O(145 - 155) MeV.
      This implies that, for a better agreement with the lattice results,
    the model interaction kernel must not contain the Debye mass. Recall that, for the AWW model the Debye mass
    enters in the longitudinal  interaction kernel exponentially, ∝exp(-16π^2 T^2/ω^2), hence
    sizably   suppresses the solution B at large T.  This is a hint that, generally in models based on the rainbow
     approximation, the Debye mass has to be taken into account  only in the perturbative part of
     the interaction, e.g. in the UV term. Note that in the literature there are also other attempts to modify the
     interaction kernel  and to introduce a suitable tuned additional T dependence at large T <cit.>.
     In what follows in solving the tBS equation for qq̅ bound states in
     most calculations we ignore the Debye mass. However, for illustration of
      the effects of m_g, we present some results where the solution of the tDS contains m_g as well.
    
    
     
    
    § TBS EQUATION AT FINITE TEMPERATURES
    
    
     §.§ Partial decomposition of the BS vertex function
    
    In the present paper we focus on the tBS equation for pseudo-scalar states.
     Prior an analysis  of  the tBS equation at finite temperature, we briefly  recall the calculations at
     T=0. The O(4) symmetric
     solution of Eq. (<ref>) in Minkowski space for the pseudo-scalar mesons can be written in the form
     (cf. Refs. <cit.>
        Γ(P,p)=[ F_1 +(pP) F_2 p̂ +  F_3P̂ +  F_4(p̂P̂ -P̂p̂)]γ_5,
    
     where the scalar vertex functions F_i  are functions on p^2 solely. In (<ref>) the
     notation p̂ = p^μγ_μ and P̂ = P^μγ_μ is adopted.
     It has been found that the contribution of the last term in Eq. (<ref>), proportional to the tensor
     matrix σ_μν=i[γ_μ,γ_ν]/2,
     is negligibly small <cit.> and, for the qualitative analysis we are interested in the present paper,
     can be safely omitted. Then, at T=0, there remain only three components in the BS vertex.
    
    
     To extend the  equation to finite temperatures one has to take into account  the
     broken O(4) symmetry which  requires separate consideration of  space-like and time-like products of four vectors.
     Consequently, the first three terms in Eq. (<ref>)  transform into six components:
      Γ(P,q)≡Γ̃(P,q) γ_5=1/2 [
    A_1 + A_2 P_0γ_0 - A_3 (q̃p̃) q_0 γ_0 -
    A_4 q_0 P_0 q̂̃̂  + A_5 (q̃p̃) q̂̃̂- A_6p̂̃̂ ] γ_5,
    
    where the scalar functions A_i depend separately on (P_0,q_0) and (| |,| |).
    The unit vectors q̃ and p̃ in (<ref>) are purely spatial, i.e. q̃ =(0,/||)
     and p̃=(0,/||). If one considers only the zeroth Matsubara frequency for the four-vector P,
     P=(0,),  the decomposition (<ref>) reduces to the one used before     for calculations of
     screening masses <cit.>.
      To pass to Euclidean space, recall that for a meson at rest at T=0 the four product (Pq)  transforms as
      (P_0q_0)_M → -(P_4q_4)_E=-iM_qq̅q_4, where q_4 is the integration variable in the BS equation.
      Within the rainbow approximation, the interaction kernel does not depend on the total momentum P.
       Therefore,  in the tBS equation, P_4 plays the role of an external parameter
       which defines the pole of the  two-particle Green function, P_4^2 = -M^2_qq̅. If the two-particle system
        is not at rest  then  P_4^2 +  P^2= (iE)^2 +  P^2=-M^2_qq̅,
       where E is the total energy of the meson.
      At finite temperatures the Feynman rules in Euclidean space <cit.> result in the
      same formal procedure for the tBS equation
     as the one used in deriving the tDS equation, i.e. formally, the relative momentum becomes
     discrete, q=(, q_4)→ q_n=( , ω_n)
     and the integration over q_4 is replaced by summation over the Matsubara frequencies
      ω_n, cf. Eq. (<ref>). The total energy of the meson becomes also discrete,
      iE→ iE_N=i Ω_N=2π i N T, where Ω_N is the Matsubara frequency for bosons,
       with N ∈ℤ.
     Then  the BS equation in Euclidean space reads
      Γ̃(P_N,p_n)   =4/3 T∑_m ∫d^3 q/(2π)^3γ_μ S^(+)(1) Γ̃(P_N,q_m)
     S̃^(-)(2) γ_ν D_μν (κ_mn),
     
    where Γ̃ is defined by Eq. (<ref>),
    κ_mn=(-,ω_n-ω_m) and S̃^(-)(2) =γ_5   S^(-)(2) γ_5.
    Correspondingly, the quark propagators S^(+)(1) and S̃^(-)(2) are
          
    S^(+)(1)=iγ⃗·  (/2  +  )σ_V(1) -iγ_4·  (   P_4/2  +q_4 )σ_C(1)+σ_S(1), 
          S̃^(-)(2)=iγ⃗· ( /2  -  )σ_V(2)+
     iγ_4·  (   P_4/2  -q_4) σ_C(2)+ σ_S(2),
     
    where σ_V,C,S are defined by the solution of the tDS equation at the same temperature T,
    
        σ_V,(C,S)(1,2)=A(_1,2,ω_m_1,2), (C(_1,2,ω_m_1,2),B(_1,2,
        ω_m_1,2))^2 A^2(_1,2,ω_m_1,2) + ω_m_1,2^2 C^2(_1,2,ω_m_1,2)+B^2(_1,2,ω_m_1,2).
    
    
     §.§ Angular integration
    
    As   seen in Eqs. (<ref>)-(<ref>)   the tBS equation
    implicitly  depends on three spatial solid angles,
     (θ_,φ_), (θ_,φ_) and (θ_,φ_).
    The dependence on (θ_,φ_) and (θ_,φ_) follows from
     the interaction kernel D(κ^2=(-)^2+(ω_m-ω_n)^2 ), which consists of two parts,
    D_1(κ^2)∼ e^-κ^2/ω^2 and D_2=κ^2/ω^2 D_1(κ^2))≡ακ^2exp(-ακ^2).
    The  dependence on (θ_,φ_) and  also on (θ_,φ_) comes from the
    propagator functions σ_F_1(1) and
    σ_F_2(2).
     The angular parts of D_1 and D_2 of the kernel can be handled by
     decomposing them over the  spherical harmonics
     Y_lm() and Y^*_lm() as
    
    D_1(κ^2)≡ e^-κ^2/ω^2=
     e^-[(ω_n - ω_m)^2 +(|| - ||)^2]/ω^2
     4π∑_lm f_l^(s)(2|| ||/ω^2)  Y_lm() Y^*_lm()
     
    where the coefficients f_l^(s)(2|| ||/ω^2) (here ω^2 denotes the slope
    parameter of the AWW kernel, not to be confused with a Matsubara frequency ω_n) are
    proportional to the scaled spherical  Bessel functions of the first kind, I_n(x),
    
    f_l^(s)(x)=√(π/2x) e^-xI_l+1/2(x),
    
    with x≡ (2|| ||/ω^2).
    The second part of the kernel can be inferred from  Eq. (<ref>) by observing that
    D_2(ακ^2)=-αd  D_1(ακ^2)/dα, where α=1/ω^2.
    
    The angular dependence on (θ_,φ_) and  (θ_,φ_), which comes from the
    propagator functions, is more involved. Since   σ_F(1,2) are numerical solutions of
    the tDS equation, an analytical expression for the decomposition over the corresponding spherical
    harmonics is lacking. Nevertheless, we decompose the product of two propagator functions σ_F_1(1) and
    σ_F_2(2) as
    σ_F_1(p_1) σ_F_2(p_2)=∑_L σ_L^F_1F_2(ω_n,Ω_N, ||,|| )
    ∑_M Y_LM()  Y_LM^*() ,
    
        where p_1,2=q ±1/2 P,  and the coefficients
        σ_L^F_1F_2(ω_n,Ω_N, ||,|| ) have to be computed numerically
        from the solution of the tDS equation
    σ_L^F_1F_2(ω_n,Ω_N, ||,|| )=2π∫ d(cosθ_Pq)
    σ_F_1(p_1)σ_F_2(p_2)    P_L(cosθ_Pq).
    
    In the  expressions above, P_L(cosθ_Pq) are the Legendre polynoms;
      the superscripts  F_1F_2 denote the decomposition of  diagonal F_1F_2 = VV, SS, CC
    and   non-diagonal products VS,  SV,  VC,  CV,  CS,  SC of the propagator functions, respectively.
    Now we are in a position to reduce the tBS equation to a system of linear algebraic equation and   solve it
    numerically.
    
    
    
    § NUMERICAL METHODS
    
    Equation (<ref>) is a 4× 4 matrix in the Dirac spinor space.
    To find the corresponding equations for the partial components A_i (i=1… 6) we consecutively multiply the
    l.h.s. and r.h.s. of the equation by the matrices
          [       G̃_1=1/2Î,             oooo    G̃_2 =1/2γ_4,; G̃_3 =i1/2γ⃗/||,  G̃_4= i1/2γ⃗/|| ]
    and compute the corresponding traces and find the integral equations of the partial
    components A_i as a  system of three-dimensional integral equations
    with an infinite  summation over the Matsubara frequencies. Prior to the analytical   angular
    integration, we solve the tDS numerically to find the coefficients
    σ_L^F_1F_2(ω_n,Ω_N, ||,|| ) also numerically.
     Then we use the decomposition (<ref>) for the interaction kernel and integrate
     analytically  the resulting expressions. Note  that additional angular dependencies emerge
     also from the traces which, together with (<ref>) and (<ref>), make the final expression
     rather cumbersome. It may contain products of up to four spherical harmonics for each solid
     angle of the three-vectors ,   and  , e.g. products of the form
     ∑_i_i,m_i(⋯)  Y_l_1m_1() Y_l_2m_2() Y^*_l_3m_3() Y^*_l_4m_4(),
     cf. Eqs. (<ref>)-(<ref>). The analytical integration results
     in a  more cumbersome expression which contains a series of products of 3j-, 6j- and 9j-symbols to be
      summed up over all l_i m_i for each vector ,   and  .
    Calculations can be essentially simplified if one uses packages for analytical calculation.
     We employ the Maple packages to calculate traces and to manipulate the  Racah
    symbols <cit.> to perform explicitly the summation over all quantum numbers.
    The remaining integration over the momentum || is calculated by
    discretizing the integral with a Gaussian mesh. The interval for
     ||=[0,∞]  is truncated by a sufficiently large value of ||_max= O(10-15) GeV.
      To re-arrange the Gaussian nodes closer
    to the origin ||∼ 0 we apply an appropriate  mapping of the  mesh by changing the variables as
    
    y=x_0 1+x/1-x
    with x_0 chosen to assure ||=||_max at the last Gaussian node. In our calculations
    we use  Gaussian meshes with N_G=36 and N_G=48 points, respectively.
     Eventually,  the resulting system   of linear
    equations reads schematically
    
        X=   SX,
    
    where  the vector
    
        X^T= (
        [{ A_1^n(q_i)}_i=1^N_G]_n=1^M_max,
         [{A_2^n( q_i)}_i=1^N_G]_n=1^M_max,
         …,
         [{A_6^n( q_i)}_i=1^N_G]_n=1^M_max ),
    
     for a given Matsubara frequency Ω_N and given value of | |,
    represents the sought solution in the form of a group of sets of
    partial wave components A_α^n(q_i),
      specified on the integration mesh of the order N_G and the maximum M_max
      for the Matsubara frequencies ω_n. In our calculations we use M_max= 25 - 30, which is
      quite sufficient for summations over ω_n within the AWW model with only the Gaussian-like IR term.
     Then the matrix S is of   dimension N_S× N_S, where N_S=α×M_max× N_G with
     α=6.
    Since the system  (<ref>) is homogeneous,
    the  eigenvalue  solution  is obtained from the condition Δ=(S-𝕀)=0.
    At fixed T and Matsubara frequency
    N,  only the free external parameter remains as the three-vector of the
    two-quark system ||. The zeros of the determinant of S-1
    are sought by changing   || from a minimal value ||∼ 0 to a maximum value ||=2π N T.
    More details about the numerical algorithm of solving the BS equation have been reported in
    Refs. <cit.>. Note  that, with the adopted truncations for the Matsubarra  summation and
    the employed   Gaussian  meshes, the resulting dimension of the determinant is rather large, N_S∼ 5000 - 8000.
    Consequently, calculations of the  elements and the determinant itself are
    lengthy and require some computer time. To increase the computer efficiency
    in  solving the tDS equation and in computing the matrix elements of the matrix S we use parallel computations
    on a PC farm within the MPI ( Message Passing Interface) standard package with a
     large enough amount of cores  (∼ 250-300). The efficiency  of such calculations is directly
    proportional to the number of cores used in parallel.
     A more involved  situation occurs for the  determinant Δ. A straightforward  use of parallel
     computation  is hindered by the fact that, if one employs  the pivoting in the Gaussian elimination method,
     the communication time among processors
    increases  rapidly with an increase of number of  processors  and could be even larger than the calculations with a single
    core. We developed a method of computing determinants in parallel with an optimization of the
    computing and communication time in such a way that, in  our case of large dimensions,
     the total   time  can be reduced by a factor of 10-20. The corresponding algorithm and
    the employed code will be reported elsewhere.
    
    
    
    § RESULTS
    
    We solve the tBS equation for pseudo-scalar ground states for two values of the Matsubara frequency N=1 and N=2.
    At   given temperature and Matsubara frequency the possible value of the pole mass is restricted
    to the interval 4π^2 N^2 T^2 < M_qq̅^2 ≤ 0 which corresponds to 0<||^2 < 4π^2 N^2 T^2.
    The maximum value of the pole mass corresponds to the limit of thermal mass, i.e. ||→ 0.
    For the lowest Matsubara frequency, N=1, this limit can occur already at   T∼ 100 MeV,
     which implies that the solution of tBS at large temperatures approaches the thermal limit ||→ 0.
    In most of our calculations, the Debye mass has been omitted. Results of calculations for the ground state
    pseudo-scalar pole masses are presented in Fig. <ref> for the Matsubara frequencies  N=1 (solid curve)
    and N=2 (dashed curve) as functions of temperature. In both cases, the mass rapidly increases with increasing
    temperature, and already at T≥ 100 MeV (for N=1) and T ≥ 80 MeV (for N=2)
    becomes larger than the maximum value of the  mass of two quasi-free quarks at the same value of T.
    This is demonstrated in Fig. <ref>, where the mass of two quarks, defined as the solution of
    tDS equation, M_2q=2B( 0,ω_0)/A( 0,ω_0), is presented by the dotted curve.
    
    
    
     
    
    For temperatures, where M_qq̅ > M_2q,  the bound ground  state, in the "canonical" sense, does not
    exists or can be considered as unstable against the dissociation into two correlated quasi-free quarks.
    
     Now we focus on the case of large temperatures for the solution with  Matsubara frequency N=1.
    As mentioned above, the maximum value of the pole mass is limited by the
    thermal mass at =0,  i.e. M_qq̅→ 2π N T.
    In  Fig. <ref>, we present an illustration of such a limit, where the real and imaginary parts
     of the determinant Δ are depicted as a function  of M_qq̅ at two values of the temperature, T=110 MeV and T=120 MeV. It can be seen that
    the determinant converges to zero, i.e. to the solution of the tBS equation, exactly in the region of the
    corresponding thermal mass. In   Fig. <ref>,
      this is depicted by  two vertical arrows. Note that, at large values of
    T  and for N=1, we have not found  solutions of the  tBS equation  with smaller pole masses.
    
    
    
    
     We investigate now also the influence  the Debye mass on the bound states for N=1. For this
    we consider the tBS equation with the solution of tDS equation with Debye mass taken into account.
     Results are exhibited in Fig. <ref>, where the solid curve denotes the pole mass without the Debye mass,
     cf. Fig. <ref>, and the dashed curve illustrates the effect of the Debye mass.
     As expected, at low temperatures,
     where m_g is small, the two solutions are practically the same. At intermediate values of T,  the interaction
     kernel is more compact, hence  it leads to larger values of the masses. At larger temperatures
     both solutions approach the thermal mass limit, see Fig. <ref>.
     For completeness, in Fig. <ref> we also present the values of two quasi-free quark masses
     M_2q≡ 2B( 0,ω_0)/A( 0,ω_0), as solutions of tDS equation with
     Debye mass (dot-dashed curve) and without Debye mass (dotted curve).
      It is seen that with Debye mass taken into account the
     pole mass becomes larger than the two free masses at already T∼ 90 MeV, i.e.
     the dissociation occurs at lower temperatures. This is in
     agreement with the behaviour of the inflection points χ_B
      of the mass solution B   and χ_qq̅ of the quark condensate qq̅, cf.  Fig. <ref>.
    
    
    
    
    
    § SUMMARY
    
     We have investigated   the solution of the truncated Dyson-Schwinger  (tDS)
     equation  at finite temperature within the rainbow approximation by employing
     the Alkofer-Watson-Weigel (AWW) model which consists  only
     of the infra-red term in a more complex interaction kernel. The  solution
     of the tDS equation is a prerequisite for a
     consistent solution of the truncated Bethe-Salpeter (tBS) equation for
     quark-antiquark bound states at finite temperature within the
     same approximation.   The ultimate goal is to use  the tBS equations,
      for an analysis of the behaviour of hadrons in hot matter, including possible (phase) transitions
     and  dissociation effects. For this goal we investigate to what extent the model, which provides a fairly good
     description of ground-state mesons at zero temperatures, can be applied to the truncated tDS equation at finite temperatures.  We find that a direct inclusion of the Debye mass  in the interaction kernel results
      in a too low value of the pseudo-critical temperature T_c, defined as the inflection point of the mass function
      B( 0,ω_0) or/and of the normalized quark condensate qq̅, being by O(50%) smaller than  the ones found in lattice QCD calculations. Since in the considered model at a given temperature   T,  the  Debye mass  m_g enters into the longitudinal part of the interaction
      as a (constant) Gaussian exponential factor ∝exp(-16π^2 T^2/5ω^2),
       it makes problematic the attempts of obtaining larger values  of the QCD deconfinement temperature
       T_c close to  the lattice values.
       We argue that, since the Debye mass m_g is essentially a perturbative quantity, it
      has to be taken into account in the perturbative ultra-violet term to maintain a correct behaviour
     at large momenta, while in the nonperturbative infra-red term it can be omitted at all. We find that by
      omitting the Debye mass m_g  one is able to obtain critical temperatures closer
      (smaller by only 10-15%) to that from lattice 2+1 flavour QCD.   To achieve a better agreement with lattice data the  IR term requires further meticulous analyse of the tDS equation at finite
      temperatures  which will be done elsewhere.
    
      As in the case of zero temperature, T=0, we consider the tBS equation within the rainbow approximation with the quark propagators
      entering the tBS equation found preliminarily as solutions of the tDS equation at the same
      temperature.  In our calculations, we restrict the boson Matsubara frequency Ω_N=2π N T to N=1 and N=2. Larger values
      of N provide  large values of M_q q̅ (N>2), far above masses of known lightest pseudo-scalar mesons.  For each Matsubara frequency Ω_N,   the ground state  mass M_qq̅^2 = 4π^2 N^2 T^2-^2
       is defined as the first (lowest) zero of the corresponding determinant as a function of ||. Such  a mass is referred to as  the Matsubara pole mass.
    
      For both frequencies  Ω_N, N=1 and N=2, we find that the pole
      masses rapidly increase with increase of  temperature. This is in agreement with the behaviour
      of the screening masses at large temperatures reported in literature, see e.g.
      Refs. <cit.>.
    It has been found that, for the lowest Matsubara frequency, where  N=1,  a solution of the tBS
    equation at large temperatures, T>110 MeV, can be found only in  the limit of thermal masses, i.e. at
      || → 0.  We analyse  the behaviour of the pole masses as a function of temperature
      with respect to the T dependence of the masses of free quarks obtained from the
      tDS equation.
      At large values of T, T>100 MeV,  the  pole masses become larger than the sum of two  quarks.
      This implies that at large T the ground state of two quark does not
      occur in the sense   as commonly adopted  in quantum mechanics where the binding energy is negative.
      This can be interpreted as dissociation instability of the state
      against a fragmentation  into a state  of two quasi-free quarks.
       The Matsubara pole masses do not have a direct relation to the inertial masses, however, one can
       expect that the obtained behaviour of the pole masses follow the same tendency as the
       inertial masses at the same temperature.
    
    
    
    
    
    § ACKNOWLEDGMENTS
    
    This work was supported in part by the Heisenberg - Landau program
    of the JINR - FRG collaboration.  LPK appreciates the warm hospitality at the
    Helmholtz Centre Dresden-Rossendorf. The authors acknowledge discussions with Dr. J. Vorberger.
    
    99fishJPA
                               C. S. Fischer, S. Kubrak,  R. Williams,
                               Spectra of heavy mesons in the Bethe-Salpeter approach,
                               Eur. Phys. J.  A51  (2015) 10.
    Thomas
                               T. Hilger,  M. Gómez-Rocha, A. Krassnigg,
                               Light-quarkonium spectra and orbital-angular-momentum decomposition in 
                               a Bethe-Salpeter equation approach,
                               Eur. Phys. J.  C77  (2017)  625.
     rob-1
                   P. Maris, C. D. Roberts, Pi- and K meson Bethe-Salpeter amplitudes,
                   Phys. Rev. C  56 (1997)  3369.
    MT
                           P. Maris, P. C. Tandy,  Bethe-Salpeter study of vector meson masses and decay constants,
                           Phys. Rev.    C 60  (1999)  055214.
    ourLast
                        S. M. Dorkin, L. P. Kaptari, B. Kämpfer,   	Accounting for the analytical properties of the quark propagator from the Dyson-Schwinger equation,
                        Phys. Rev.  C 91 (2015)  055201.
    DorkinByer
                    S. M. Dorkin, M. Beyer, S.S. Semikh, L.P. Kaptari,
                    Two-Fermion Bound States within the Bethe-Salpeter Approach,
                     Few Body Syst.  42 (2008) 1.
    
     Hilger
                     T. Hilger, C. Popovici, M. Gomez-Rocha, A. Krassnigg,
                     Spectra of heavy quarkonia in a Bethe-Salpeter-equation approach,
                      Phys. Rev.  D 91  (2015)  034013.
    Maris:2003vk
                            P. Maris, C. D. Roberts,  Dyson-Schwinger equations: A tool for hadron physics,
                             Int. J. Mod. Phys.  E 12 (2003) 297.
    Holl:2004fr
                         A. Holl, A. Krassnigg, C. D. Roberts,  Pseudoscalar meson radial excitations,
                            Phys. Rev.  C 70  (2004) 042203.
    Blank:2011ha
                       M. Blank, A. Krassnigg,   Bottomonium in a Bethe-Salpeter-equation study,
                       Phys. Rev.  D 84 (2011) 096014;
    
                       M. Blank, A. Krassnigg, A. Maas,
                       Rho-meson, Bethe-Salpeter equation, and the far infrared,
                       Phys. Rev.  D 83 (2011) 034020.
    Jarecke:2002xd
                             D. Jarecke, P. Maris, P. C. Tandy,    Strong decays of light vector mesons,
                             Phys. Rev.  C 67 (2003) 035202.
    Krassnigg:2004if
                              A. Krassnigg, P. Maris,  Pseudoscalar and vector mesons as q anti-q bound states,
                                J. Phys. Conf. Ser.  9  (2005) 153.
    Roberts:2007jh
                        C. D. Roberts, M. S. Bhagwat, A. Hol, S. V. Wright,   Aspects of hadron physics,
                        Eur. Phys. J. ST  140 (2007) 53.
    ourFB
                  S. M. Dorkin, T. Hilger, L. P. Kaptari, B. Kämpfer,   Heavy pseudoscalar mesons in a Dyson-Schwinger–Bethe-Salpeter approach,
                  Few Body Syst.  49  (2011) 247;
    
                   S. M. Dorkin, L. P. Kaptari, C. Ciofi degli Atti, B. Kämpfer,   Solving the Bethe-Salpeter Equation in Euclidean Space,
                    Few Body Syst.  49  (2011) 233.
    wilson
                    S. -x. Qin, L. Chang, Y. -x. Liu, C. D. Roberts, D. J. Wilson,   Interaction model for the gap equation,
                    Phys. Rev. C  84  (2011) 042202; 
    
                     S. -x. Qin, L. Chang, Y. -x. Liu, C. D. Roberts,  D. J. Wilson,  Investigation of rainbow-ladder truncation for excite`d and exotic mesons,
                    Phys. Rev.  C 85 (2012) 035202.
    
    JM-1
                   P. Jain, H. J. Munczek,
                   Calculation of the pion decay constant in the framework of the Bethe-Salpeter equation,
                   Phys. Rev.  D 44  (1991) 1873.
    JM-2
                   H. J. Munczek, P. Jain,
                   Relativistic pseudoscalar q anti-q bound states: Results on Bethe-Salpeter wave functions and decay constants,
                   Phys. Rev    D 46 (1992) 438.
    rob-2
                   M. R. Frank, C. D. Roberts,
                   Model gluon propagator and pion and rho meson observables
                    Phys. Rev.     C 53  (1996) 390.
    Alkofer
                     R. Alkofer, P. Watson, H. Weigel,
                      Mesons in a Poicare Bethe Salpeter approach,
                     Phys. Rev.   D 65 (2002) 094026.
    fisher
                      C. S. Fischer, P. Watson, W. Cassing,
                      Probing unquenching effects in the gluon polarisation in light mesons,
                      Phys. Rev.    D 72  (2005) 094025.
      physRep
                   R. Alkofer,  L. von Smekal,
                   The Infrared behavior of QCD Green's functions: Confinement dynamical symmetry breaking, and hadrons as relativistic bound states,
                   Phys. Rept.  353  (2001) 281.
    smekal
              A. Bender, C. D. Roberts, L. von Smekal,
              Goldstone theorem and diquark confinement beyond rainbow ladder approximation,
              Phys. Lett.   B 380 (1996) 7.
    scadron
                R. Delbourgo, M. D. Scadron,
                Proof of the Nambu-goldstone Realization for Vector Gluon Quark Theories,
                J. Phys. G: Nucl. Phys.  5  (1979) 1621.
    munczek
               H. J. Munczek,
               Dynamical chiral symmetry breaking, Goldstone's theorem and the consistency of the Schwinger-Dyson and Bethe-Salpeter Equations,
               Phys.  Rev.  D 52 (1995) 4736.
    hatsuda
                               T. Hatsuda,Y. Koike, S. H. Lee,
                               Finite temperature QCD sum rules reexamined: rho, omega and A1 mesons,
                               Nucl. Phys.  B394 (1993) 221.
    Ayale
                            A. Ayala,  A. Bashir,
                            Longitudinal and transverse fermion boson vertex in QED at finite temperature in the HTL approximation,
                            Phys. Rev.  D 64  (2001) 025015.
    RobProg
                          C. D. Roberts, S. M. Schmidt,
                          Dyson-Schwinger equations: Density, temperature and continuum strong QCD,
                          Prog. Part. Nucl. Phys. 45 (2000) S1-S103.
    matzubara
                      T. Matsubara,
                      A New approach to quantum statistical mechanics,
                      Prog. Theor. Phys.,  14  (1955) 351.
    periodic
                 R. Haag, N. Hugenholtz, M. Winnink,  On the Equilibrium States in Quantum Statistical Mechanics,
                 Commun. Math. Phys.  5 (1967) 215.
    KMS
              J. Derezinski, C. -A. Pillet, KMS states,
              pillet.univ-tln.fr/data/pdf/KMS-states.pdf
     kapusta
                         J. I. Kapusta, "Finite-Temperature Field Theory", Cambridge University Press, N.Y., 1989.
    abrikosov
                     A. A. Abrikosov, L. P. Gorkov, I. E. Dzyaloshinski,
                     "Methods of Quantum Field Theory in Statistical Physics", (editor) R. A. Silverman,
                     PRENTICE-HALL, INC., Englewood Cliffs, New  Jersey, 1963.
    fisher1
              C. S. Fischer, J. M. Pawlowski, A. Rothkopf, C. A. Welzbacher,          . 	
              Bayesian analysis of quark spectral properties from the Dyson-Schwinger equation,
              e-Print: arXiv: 1705.03207 [hep-ph].
    roberts
                     Kun-lun Wang,Yu-xin Liu, Lei Chang, C. D. Roberts, S .M. Schmidt,
                     Baryon and meson screening masses,
                      Phys.Rev.  D 87 (2013) 074038.      
        yaponetz
                     M. Ishii, H. Kouno,  M. Yahiro,  Model prediction for temp. dep. of meson pole masses from LQCD on meson screening masses,
                     Phys. Rev.  D 95 (2017) 114022.
    robertsSchwinger
                     C. D. Roberts,
                     Hadron Properties and Dyson-Schwinger Equations,
                    Prog. Part. Nucl. Phys.  61 (2008) 50.
    physrepMorley
                    P. D. Morley, M. B. Kislinger,
                    Relativistic Many Body Theory, Quantum Chromodynamics and Neutron Stars/Supernova,
                     Phys. Rept.  51 (1979) 63.
    blaske
                   D. Blaschke, A. Dubinin, A. Radzhabov, A. Wergieluk,
                   Mott dissociation of pions and kaons in hot, dense quark matter,
                    Phys.Rev.  D96 (2017) 094008.
    frg1
                          J. M. Pawlowski, N. Strodthoff, N. Wink,
                         	Finite temperature spectral functions in the O(N)-model,
                          E-print: arXiv 1711.07444 [hep-th].
    frg2
                                  N. Strodthoff,
                                  Self-consistent spectral functions in the O(N) model from the functional renormalization group,
                                  Phys. Rev.  D95 (2017) 076002.
    frg3
                              J.M. Pawlowski, N. Strodthoff,
                              Real time correlation functions and the functional renormalization group,
                              Phys. Rev.  D92 (2015) 094009.
    Parameterslattice
           P. Maris, C. D. Roberts,
           QCD modeling of hadron physics,
          Nucl. Phys. Proc. Suppl.  161 (2006) 136-152.
    williams
                                 C. D. Roberts, A. Williams,
                                 Dyson-Schwinger equations and their application to hadronic physics,
                                  Prog. Part. Nucl. Phys.  33  (1994) 477.
    kitaizyRpberts
                            S-X Qin, L.  Chang, Y-X Liu, C. D. Roberts,
                            Quark spectral density and a strongly-coupled QGP,
                            Phys. Rev.  D 84  (2011) 014017.
    MTrenorm
                         P. Maris, C. D. Roberts, S. M. Schmidt, P. C. Tandy,
                         T - dependence of pseudoscalar and scalar correlations,
                         Phys. Rev.  C 63  (2001) 025202.
    rischke
                            Si-xue Qin, Dirk H. Rischke,
                            Quark Spectral Function and Deconfinement at Nonzero Temperature,
                            Phys. Rev.  D88  (2013) 056007.
    ourModernPhys
                         S. M.Dorkin,  M. Viebach, L. P. Kaptari, B. Kämpfer,
                         Extending the truncated Dyson-Schwinger equation to finite temperatures,
                         J. Mod.Phys. 7 (2016) 2071; e-Print: arXiv:1512.06596 [nucl-th].
    BlankKrass
                      M. Blank, A. Krassnigg,
                      The QCD chiral transition temperature in a Dyson-Schwinger-equation context,
                       Phys. Rev.  D 82  (2010) 034006.
    robLast
                     Fei Gao, Jing Chen, Yu-Xin Liu, Si-Xue Qin, Craig D. Roberts, Sebastian M. Schmidt,
                     Phase diagram and thermal properties of strong-interaction matter,
                      Phys. Rev.  D 93  (2016), 094019;
    
                      F. Gao, S.-X. Qin, Y.-X. Liu, C. D. Roberts, S. M. Schmidt,
                      Zero mode in a strongly coupled quark gluon plasma,
                     Phys. Rev.  D 89  (2014) 076009.
     kitazavaPRD80
                              F. Karsch, M. Kitazawa,
                              Quark propagator at finite temperature and finite momentum in quenched lattice QCD,
                              Phys. Rev.  D 80  (2009) 056001.
    fischerPRD90
                  C. S. Fischer, J. Luecker, C. A. Welzbacher,
                  Phase structure of three and four flavor QCD,
                       Phys. Rev.  D 90  (2014) 034022.
       FischerRenorm
                          C. S. Fischer, J. A. Mueller,
                          Chiral and deconfinement transition from Dyson-Schwinger equations,
                           Phys. Rev.  D 80  (2009) 074029.
    fischerGluon
                        G. Eichmann, C. S. Fischer,  C. A. Welzbacher,
                        Baryon effects on the location of QCD’s critical end point,
                        Phys. Rev.   D 93  (2016)  034013;
    
                        C. S. Fischer, J. Lueckerb, C. A. Welzbacher,
                        Locating the critical end point of QCD,
                        Nucl. Phys.  A 931  (2014) 774.
    otherKernels1
                        Hui-Feng Fu, Qing Wang,
                        Quark propagator in a truncation scheme beyond the rainbow approximation,
                        Phys. Rev.  D 93 (2016) 014013.
    otherKernels2
                        D. Binosi, Lei Chang, J. Papavassiliou, Si-Xue Qin, C. D. Roberts,
                        Symmetry preserving truncations of the gap and Bethe-Salpeter equations,
                        Phys. Rev.  D 93 (2016) 096010.
    otherKernels3
                         B. El-Bennich, G. Krein,  E Rojas, F E. Serna,
                         Excited hadrons and the analytical structure of bound-state interaction kernels
                          Few Body Syst.  57 (2016) 955.
    yoki
              Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg, K. K. Szabo,
              The QCD transition temperature: results with physical masses in the continuum limit II,
              JHEP 0906 (2009) 88.
    Fischer:2008uz
                           C. S. Fischer, A. Maas,  J. M. Pawlowski,  On the infrared behavior of Landau gauge Yang-Mills theory,
                          Ann. Phys.  324  (2009)  2408.
    PenningtonUV
                              M. R. Pennington, D. J. Wilson,
                               Are the Dressed Gluon and Ghost Propagators in the Landau Gauge presently determined in the confinement regime of QCD?,
                             Phys. Rev.  D 84 (2011)  094028;
                                 [ibid.  D 84 (2011) 119901, Erratum]. 
    souglasInfrared
                              N. Souchlas,
                               	Infrared and Ultraviolet QCD dynamics with quark mass for J=0,1 mesons,
                               arXiv:nucl-th/1006.0942, June 2010.
      umezawa
                H. Umezawa, H. Matsumoto, M. Tachiki, "Termo Field Dynamics and Condensed States"
                 (North-Holland, Amsterdam, 1982).
    nemi
           A. J. Niemi,  G. W. Semenoff,
           Finite Temperature Quantum Field Theory in Minkowski Space,
            Ann. Phys. 152 (1984) 105.
    landshof
                Michel le Bellac,
               "Thermal Filed Theory",     Cambridge Monographs on Mathematical Physics,
               General Editors:  P. V. Landshoff, D. R. Nelson, D. W. Sciama, S. Weinberg,
               Cambridge University Press 1996 .
    Sch
                 P. C.  Martin, J. Schwinger,
                 Theory of many particle systems. 1.
                 Phys. Rev., 115 (1959) 1342.
    landsman
                      N. P. Landsman,  Ch. G. van Weert,
                      Real and Imaginary Time Field Theory at Finite Temperature and Density,
                      Phys. Rept.   145  (1987) 141.
    tereza
                    A. Cucchieri, A. Maas, T. Mendes,
                    Infrared properties of propagators in Landau-gauge pure Yang-Mills theory at finite temperature,
                     Phys. Rev.  D 75  (2007) 076003.
    Tagaki1
                    M. Harada, S. Tagaki,
                    Phase transition in two flavor dense QCD from the Schwinger-Dyson equation,
                           Prog.  Theor. Phys.  107  (2002) 561.
    meiss
                     D. K. Hong, V. A. Miransky, I. A. Shovkovy, L. C. R. Wijewardhana,
                     Schwinger-Dyson approach to color superconductivity in dense QCD,
                      Phys. Rev.  D61 (2000), 056001, [ibid. D 62,  (2000),059903, Erratum].
    AlkoferDebye
                       R. Alkofer, P. A. Amundsen,   K. Langfeld,
                       Chiral Symmetry Breaking and Pion Properties at Finite Temperatures,
                       Z. Phys.  C 42  (1989) 199.
    Tagaki
                    S. Tagaki,
                    Phase structure of hot and / or dense QCD with the Schwinger-Dyson equation,
                     Prog. Theor. Phys.  109  (2003) 233.
    kalinovsky
                                 D. Blaschke, G. Burau, Yu. L. Kalinovsky, P. Maris, P. C. Tandy,
                                 Finite T meson correlations and quark deconfinement,
                                Int. J . Mod. Phys.  A16 (2001) 2267-2291
               A. Bender, D. Blaschke, Yu. Kalinovsky, C. D. Roberts,
               Continuum study of deconfinement at finite temperature,
               Phys. Rev. Lett.   77 (1996) 3724.
    our2013
                   S. M. Dorkin, L. P. Kaptari, T. Hilger, B. Kämpfer,
                   Analytical properties of the quark propagator from a truncated Dyson-Schwinger equation in complex Euclidean space,
                    Phys. Rev.  C 89 (2014) 034005.
    Llewellyn
                     S. Llewellyn,
                     A relativistic formulation for the quark model for mesons,
                     Ann. of  Phys.  53 (1969) 521.
    Fischerrenorm1
                          C. S. Fischer, J. Luecker,
                          Propagators and phase structure of Nf=2 and Nf=2+1 QCD,
                          Phys. Lett.  B 718  (2013) 1036.
    racah
                   S. Fritzsche,
                   Maple procedures for the coupling of angular momenta. An up-date of the Racah module,
                    Computer Physics Communication  180 (2009) 2021;
                    Program summary URL:
