{'fraction_non_alphanumeric': 0.03488049966000957, 'fraction_numerical': 0.01337295690936107, 'mean_word_length': 4.802718106093819, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nuclear superfluidity at finite temperature

The equation of motion for the two-fermion two-time correlation function in the pairing channel is considered at finite temperature. Within the Matsubara formalism, the Dyson-type Bethe-Salpeter equation (Dyson-BSE) with the frequency-dependent interaction kernel is obtained. Similarly to the case of zero temperature, it is decomposed into the static and dynamical components, where the former is given by the contraction of the bare interaction with the two-fermion density and the latter is represented by the double contraction of the four-fermion two-time correlation function, or propagator, with two interaction matrix elements. The dynamical kernel with the four-body propagator, being formally exact, requires approximations to avoid generating prohibitively complicated hierarchy of equations. We focus on the approximation where the dynamical interaction kernel is truncated on the level of two-body correlation functions, neglecting the irreducible three-body and higher-rank correlations. Such a truncation leads to the dynamical kernel with the coupling between correlated fermionic pairs, which can be interpreted as emergent bosonic quasibound states, or phonons, of normal and superfluid nature. The latter ones are, thus, the mediators of the dynamical superfluid pairing. In this framework, we obtained the closed system of equations for the fermionic particle-hole and particle-particle propagators. This allows us to study the temperature dependence of the pairing gap beyond the Bardeen-Cooper-Schrieffer approximation, that is implemented for medium-heavy nuclear systems. The cases of 68Ni and 44,46Ca are discussed in detail.


I. INTRODUCTION
Superfluidity and superconductivity phenomena in nuclear systems, after their recognition in late 1950's [1], attracted a tremendous amount of theoretical effort since then. Although the theory of Bardeen, Cooper and Schrieffer (BCS) [2] appeared to be very successful when applied to superconductivity in metals, building a consistent theory of nuclear pairing correlations was complicated by the nature of the nuclear forces. In particular, the presence of the repulsive core in the nucleon-nucleon interaction and strong in-medium correlations made the direct applicability of the BCS theory to nuclear matter and finite nuclei problematic. The development of powerful many-body methods, such as numerous variants of perturbative and cluster expansions, the correlated basis function method, the Monte Carlo approach and others, together with the advancements of the nucleon-nucleon potentials, has helped significantly to clarify microscopic mechanisms of nuclear superfluidity, eventually going far beyond the BCS theory [3].
While the observations, such as the rotational anomalies in pulsar periods and measurements of their surface temperatures evidence unambiguously about superfluidity of neutron stars, theoretical models still vary considerably in the description of its characteristics, for instance, the pairing gaps. The use of different nucleon-nucleon interactions and regularization techniques as well as different treatments of polarization effects may cause substantial differences in model predictions. Refining the models of superfluidity in both symmetric and asymmetric nuclear matters, in particular, clarifying the role of induced pairing in screening and antiscreening is a topic of active research [4][5][6][7][8][9][10].
Investigation of pairing correlations in finite nuclei seems to be less intense. In most of the applications to nuclear structure calculations, rather simplistic concepts of pairing like BCS, Hartree-Fock-Bogoliubov or Gor'kov Green functions are employed, which is unavoidable to make otherwise sophisticated many-body calculations feasible. The accuracy of such simplified treatments of pairing is comparable with the errors introduced by other model approximations, such as neglecting high-rank many-body correlations, multiparticle interactions and coupling to the continuum, to name a few. However, with the progress of those aspects also more accurate treatment of pairing correlations should be considered.
One of the most intriguing issues in strongly-coupled many-body systems is the emergence of collective phenomena. Understanding their significant role in the formation of the nuclear ground and excited states has been progressed impressively over the decades, since Bohr and Mottelson [11,12]. The impact of collective effects on nuclear pairing was investigated in various phenomenological frameworks [13][14][15][16][17][18], which revealed that coupling between nucleons and collective surface vibrations (particle-vibration coupling, or PVC) can be responsible for a large fraction of the pairing gap. The PVC effects identified with the major contribution to the induced pairing are widely known to be of prime importance in electronic condensed matter systems, where they induce superconductivity by reversing the sign of the repulsive Coulomb interaction [3,19,20].
Atomic nuclei embedded in stellar environments are of special interest. Their response to various changes of those environments and the associated nuclear reactions are complicated by the thermodynamical conditions, first of all, by the non-zero temperature. It is widely recognized that the modifications of nuclear properties by finite temperature can noticeably influence the star evolution scenarios. In particular, the dependence of nuclear superfluidity on temperature may play a non-negligible role for the electron capture in collapsing stars and r-process nucleosynthesis in neutron star mergers. Whereas the temperature-dependent BCS is well understood and known for the superfluid to non-superfluid phase transition at the critical temperature T c ∼ 0.6∆(T = 0), the temperature dependence of the induced pairing is more complicated as well as the induced pairing itself.
This kind of pairing at finite temperature was investigated more intensely for nuclear matter, although abundant shell-model Monte Carlo studies for finite nuclei are available, see Refs. [21,22] and references therein. In recent Ref. [10] finite-temperature calculations of the singlet pairing gap in dilute neutron matter were performed. The authors investigated the pairing gaps and the critical temperature of the superfluid phase transition. The V low-k interaction derived from the Argonne potential AV18 was employed for the static kernel of the pairing gap equation and the effective interaction from the Skyrme family was used for the RPA vertices, which determine the dynamical kernel with the induced pairing. At higher densities the full RPA lead to stronger screening than the reference Landau approximation. As previously, for instance, in the studies of the BEC-BCS crossover and the liquid-gas phase transition in hot and dense nuclear matter [23], it was pointed out that the pairing gap and the phase transition temperature are sensitive to the approximation used to describe the medium polarization effects responsible for the induced pairing.
In the present work we aim at investigating the temperature dependence of the induced pairing in finite nuclei. Technically, we consider nuclear correlation functions in the equation of motion (EOM) framework, which is one of the most universal methods known across the many areas of quantum physics from condensed matter to quantum chemistry [24][25][26][27][28][29][30]. Following our previous developments reported in Ref. [31] for the zero-temperature case, in Section II we generate the EOM for the fermionic pair, or particleparticle, propagator, but now evolving in the domain of imaginary time introduced by Matsubara [32]. We show that the four-fermion correlation function in the dynamical kernel of the resulting EOM, which is responsible for the induced pairing, can be approximated with various degrees of accuracy, in analogy with the zero-temperature case. At finite temperature, however, this kernel carries a non-trivial temperature dependence, which is different from that of the static kernel implied in the BCS theory. The impact of the latter temperature dependence on nuclear pairing gaps is studied numerically in Section III. The conclusions are drawn in the summary Section IV.

II. FERMIONIC PAIR PROPAGATOR IN A HEATED CORRELATED MEDIUM
In analogy with Ref. [31], we stay within the formalism of correlation functions, such as the Green functions, or propagators. As the propagators are directly related to observed excitation spectra and ground state properties of the manybody systems, this formalism is one of the most convenient and powerful ones in the description of phenomena that occur in strongly-coupled media. Following Matsubara [32], the temperature-dependent propagator of a fermionic pair in a heated correlated medium can be defined as a thermal average [33] with the chronological ordering T τ of the one-fermion fields in the imaginary time domain of the Wick rotated picture: The operator H is given by H = H − λN , where H is the many-body Hamiltonian with the antisymmetrized matrix elementsv 1234 = v 1234 − v 1243 , λ is the chemical potential, and N is the particle number operator. Here and in the following the number subscript denotes the full set of the single-particle quantum numbers in a given representation and the imaginary time variables τ are related to the real times t as τ = it. The fermionic fields ψ 1 and ψ † 1 satisfy the usual anticommutation relations, and the angular brackets in Eq. (1) stand for the thermal average [33,34] with the summation over the expectation values in the eigenstates of the Hamiltonian |ν weighted with the probabilities w ν of finding the system in those states within the grand canonical ensemble. In the present work the Hamiltonian of Eq. (3) is confined by the two-body interaction, however, as in the zero-temperature case and as it follows from the discussion below, the generalization to multiparticle forces is straightforward. The first equation of motion probing the evolution of the correlation function of a fermionic pair defined by Eq. (1) with the imaginary time can be generated by the differentiation of this function with respect to the first time variable τ : where we adopted the notation: for an arbitrary operator A. After the evaluation of the commutators the first EOM reads: where the single-particle energies ε 1 are ε 1 = t 11 − λ and we assumed that the working basis diagonalizes the onebody part of the Hamiltonian. The norm matrix in the pp-channel N 121 2 is the thermal average: where the one-body density matrix ρ 11 obeys ρ 11 = ψ † 1 ψ 1 = δ 11 n 1 with n 1 being, in general, the correlated fermionic occupancies at the given temperature T . In Eq. (8) we adopt the antisymmetrized Kronecker symbol δ 121 2 = δ 11 δ 22 − δ 21 δ 12 and the notation n 12 = 1 − n 1 − n 2 .
Differentiating the last term on the right hand side of the first EOM (7) F 121 2 (τ −τ ) = T τ [V, ψ 1 ψ 2 ](τ )(ψ 2 ψ 1 )(τ ) with respect to the second time argument τ generates the second EOM: Applying the operator (∂ τ − ε 1 − ε 2 ) to the first EOM (7) and combining Eqs. (7) and (9), we perform the Fourier transformation to the domain of the Matsubara's discrete energy variable ω n = 2πnT . The spectral image in this domain is defined by the relation: In this way, we obtain: where the free particle-particle propagator is introduced as: The interaction kernel of Eq. (11) has the meaning of Tmatrix and reads: The components T (0) 343 4 and T (r) 343 4 (ω n ) are, formally, the Fourier images of the two last terms on the right hand side of Eq. (9), i.e., so that T (0) is the instantaneous, or static, part of the Tmatrix and T (r) is its dynamical part. In analogy with the case of the particle-hole response [35] and the zerotemperature particle-particle response [31], Eq. (11) can be transformed to an equation of the Dyson type by introducing the kernel K(ω) irreducible with respect to the uncorrelated pp-propagator (12): or K(ω n ) = T (irr) (ω n ). Thus, at finite temperature the EOM for the propagator of a fermionic pair also acquires the form of the Dyson Bethe-Salpeter equation (Dyson-BSE) [36]. To further specify the interaction kernel of the latter equation, one has to evaluate the commutators of Eq. (14). For the static part, we find, similarly to the case of T = 0 [24,31]: where σ (2) ijkl is the correlated part of the two-body density andΣ 11 is the mean-field part of the single-particle self-energyΣ The latter ones can be absorbed in the uncorrelated propagator, so that the Dyson-BSE takes the form: 12,1 2 (ω n ) + whereG (0) 12,1 2 (ω n ) is the uncorrelated particle-particle propagator in the mean field while the kernel K does not contain the mean-field term in its static part. The full interaction kernel of Eq. (21) can then be written as: K =Ñ −1 (K (0) + K (r) )Ñ −1 , and the symbol "˜" indicates the mean-field character of the quantity. Remarkably, the static kernel has the same form as at T = 0 because of its instantaneous character, however, it depends implicitly on temperature via the fermionic densities.
The presence of the static kernel is a direct consequence of the instantaneous nature of the bare interactionv, that was our initial assumption about the Hamiltonian (3). In general, the fermionic bare interaction does not have to be instantaneous, for instance, it can be mediated by a boson, whose exchange between fermions must have retardation. In the case of nuclear forces, that is the meson exchange between two nucleons in the vacuum. We will see in the following that the in-medium analog of this type of interaction can be generated in the present EOM framework, if the dynamical kernel is taken into account. Certain approximations, such as cluster decompositions of the dynamical kernel, bring its structure to the boson-exchange form, where the emergent bosons are correlated fermionic pairs, and the intermediate propagators of these bosons are associated with retardation effects of the in-medium induced interation. As we showed in detail in Refs. [31,35], the in-medium dynamical kernel in the form of the phononexchange interaction is completely analogous to the mesonexchange interaction between the nucleons in the vacuum. In turn, the latter interaction should be, in principle, derivable from the Hamiltonian of quarks and gluons. However, if one wants to rely on the scale separation and consider nucleons as elementary degrees of freedom, the consistent framework implies neglecting the time-dependence of the meson-exchange interaction. The latter approximation is widely used in the low-energy nuclear physics in the socalled "ab initio" calculations.
The most common practice for various applications of the many-body theory is to treat the EOM (21) in the simplest approximation, which retains only the static part K (0) of the kernel and neglects the correlations originating from K (r) . Such an approach forms the content of the self-consistent particle-particle random phase approximation (RPA) [37], the analogous particle-hole RPA and the self-consistent quasiparticle random phase approximation, or SCQRPA, which combines both of them and demonstrates great success in applications to exactly-solvable models [38]. In nuclear physics, moreover, the common practice is not to compute the static kernel according to Eq. (18), but rather use effective interactions adjusted to finite nuclei or the Gmatrix of the Brückner's type.
More and more applications of such approaches to nuclear systems, however, indicate that confining by only the static part of the interaction kernel can not lead to satisfactory results. The most spectacular examples are nuclear excitation spectra and associated decay properties, where the part K (r) of the kernel associated with dynamical processes induced by the medium plays a decisive role [35,39,40]. In the description of superfluid nuclear matter, this part of the kernel produces an interplay of screening and antiscreening effects which can be revealed, for instance, in the calculations of pairing gaps [5,9].
The evaluation of the commutators determining the dynamical kernel leads to the following result: which, in complete analogy to the case of the particleparticle propagator at zero temperature [31], is determined by the two-time four-fermion correlation functions contracted with two interaction matrix elements in all possible ways, which lead to a four-leg interaction kernel. Each term of Eq. (23) contains a propagator of three particles and one hole (3p − 1h). Rather than generating new EOM's for such higher-rank propagators, we will follow the approach of cluster decomposition of Eq. (23) including up to twofermion correlation functions: where we implied that [R i2,q2 G lk,nm ](τ − τ ) ≡ R i2,q2 (τ − τ )G lk,nm (τ − τ ) and the finite-temperature particle-hole response function R is introduced as Eqs. (24 -27) are, again, completely analogous to the zerotemperature case [31]. The Fourier transformation of the latter kernel to the domain of the Matsubara frequencies requires calculation of the following generic integral which yields: In Eq. (30) we have introduced the matrix elements of the normal ρ νν 12 and pairing α µν 12 , β κν 12 transition densities: where the former connects the states |ν of the given Nparticle system and the latter connect the states |ν with the states |µ , |κ of the N ±2-particle systems, respectively. The frequencies in the denominators correspond to the associated energy differences. Similarly to Refs. [31,35], it is convenient to also introduce the vertices of the emergent normal and pairing phonons as follows: where the presence of the two upper indices indicates that these vertices characterize transitions between excited states, in contrast to the case of zero temperature, where only transitions between the ground and excited states were considered. Then, the first component of the dynamical kernel takes the following form: with Θ µν;ν ν (+) 121 2 while the second component reads: where Σ µν;ν ν (+) 121 2 The two remaining components K (r;21) 121 2 (ω n ) and K (r;22) 121 2 (ω n ) can be found from Eqs. (33 -36) with the help of the symmetry relations of Eqs. (26,27). Finally, all the expressions can be analytically continued to the domain of real energies.
It is easy to see that in the approximation of Eq. (24 -27) to the dynamical kernel the many-body problem can be formulated as a closed scheme. For that, one would need to supplement Eq. (21) with an analogous EOM for the particle-hole response function (28) and for the lower-rank single-fermion propagator. With the cluster decomposition confined by the two-fermion propagators, all these propagators can be, in principle, found by a self-consistent iterative procedure. Possibilities to implement such a program for nuclear systems will be investigated elsewhere, and for the rest of this work we will focus on calculations of the quantity which is commonly referred to as pairing gap.
The equation for the pairing gap can be obtained from Eq. (21) if, for instance, the energy argument of the pair propagator is close to the transition frequency from the ground state of the N -particle system to the ground state of the (N + 2)-particle system. Then, the equation for the pairing transition density α µν = α s in the vicinity of this frequency ω s reads: If we assume that the ground state of the reference nucleus is approximated by the finite-temperature BCS-like variational ansatz, where the finite-temperature S-wave pairing gap ∆ 1 can be related to the pairing transition density as where the bar denotes the conjugate or the time-reversed state [37]. In this approximation, at the frequency ω s = 2λ Eq. (37) takes the form of the well-known pairing gap equation: which has formally the same structure as the finitetemperature BCS equation, but with a more complicated interaction kernel whose both the static and dynamical components include the Ñ −1 ...Ñ −1 factors. The dynamical part, although taken in the static limit, carries the retardation effects and the additional temperature dependence. Notice that Eq. (42) has the same form regardless the approximations made for its static K (0) and dynamical K (r) parts. As the fully self-consistent treatment of those kernels is difficult even in the approximation made for K (r) in Eqs. (24 -27), further approximations can be made. Besides the most strong BCS-like one neglecting correlations, such as the complete dynamical part K (r) and the terms with σ (2) in K (0) , one can make weaker approximations. For the static kernel K (0) this could be the G-matrix, various kinds of preprocessing of the bare interactions, such as the renormalization group or low-k, and, eventually, effective interactions. For the dynamical kernel K (r) the R, G or both correlation functions appearing in Eqs. (24 -27) can be approximated by their uncorrelated mean-field analogs or, alternatively, correlations in these propagators can be only partly relaxed. This type of approaches were applied, for instance, for the nuclear matter calculations of Refs. [5,41,42], to name a few.

III. DETAILS OF CALCULATIONS, RESULTS AND DISCUSSION
The numerical implementation of the approach of Eq. (42) for the pairing gap with the kernel of Eq. (43) aimed at the investigation of the temperature dependence of the induced pairing, i.e., essentially of the role of the second term in Eq. (43). Therefore, at this point we kept the static kernel as simple as in Ref. [31], namely described by the effective monopole-monopole force with adjustable strength to avoid complications like taming the bare interaction with the hard core. The latter will be investigated elsewhere.
As in the previous implementations of the relativistic finite-temperature approaches with PVC, first we solve the closed set of the relativistic mean field (RMF) equations using the non-linear sigma-model and the NL3 parametrization [43], where the Fermi-Dirac thermal fermionic occupation numbers self-consistently modify the classical meson fields. The procedure generates a set of temperaturedependent single-particle Dirac spinors and the corresponding single-nucleon energies, which serve as the working basis {1,ε 1 }. Then, the finite-temperature relativistic random phase approximation (FT-RRPA) equations are solved to obtain the phonon vertices g m ≡ g νν and their frequencies ω m ≡ ω νν . In this implementation we relaxed correlations in the particle-particle propagator of Eqs. (24)(25)(26)(27), because they are known to be less important than the correlations in the particle-hole propagator. This means, technically, (i) neglecting the terms with γ-vertices in the Eqs. (34,36) and (ii) replacing the pairing transition densities α and β with their uncorrelated analogs and, simultaneously, the pairing frequencies ω (±±) with the sums of the single-particle energies. Thus, the model space for the dynamical kernel is formed by the set of the obtained FT-RRPA phonons and the thermal RMF single-particle states coupled in the pp⊗phonon, hh⊗phonon and ph⊗phonon configurations. Particle-hole configurations with the energies ε ph ≤ 100 MeV and the antiparticle-hole (αh) ones with ε αh ≥ −1800 MeV, with respect to the positive-energy continuum, limited the particle-hole basis for the FT-RRPA calculations of the phonons. The set of phonons included vibrations with spins and parities J π = 2 + , 3 − , 4 + , 5 − , 6 + below 20 MeV. The phonon modes with the reduced transition probabilities B(EL) equal or more than 5% of the maximal one (for each J π ) were included in the model space, and the single-particle intermediate states entering the matrix elements g νν nk in the summations of Eqs. (34,36) with the energy differences |ε k − ε n | ≤ 50 MeV were included in the summation. The same truncation criteria were applied for all temperature regimes, that is justified by our previous calculations [44][45][46]. Note that we did not take into account the phenomenological static pairing in the mean field calculations and in the calculations of the PVC vertices as it was done in the first application of the approach at T = 0 in Ref. [31], therefore, the results at zero temperature are slightly different (see a more detailed discussion below). The main features of the solutions, however, remain intact.
The solutions obtained for the neutron pairing gap in  Fig. 1. At this temperature the result is nearly identical to that at T = 0. In this calculation we illustrate (i) the role of the PVC effects in the formation of the pairing gap in the present calculation scheme and (ii) the sensitivity of the pairing gap values to the parameter G, which defines the strength of the static part of the interaction kernel K (0) . As in the previous work of Ref. [31], the latter kernel was taken in the form of the monopole-monopole interaction, which is given in detail, for instance, in Ref. [47], so that the present study is fully focused on the features of the dynamical kernel K (r) . The parameter G is the only free parameter used in solving the pairing gap equation (42) and, ideally, will be eliminated in the ab-initio calculations, where the static kernel is determined explicitly according to Eq. (18). Here we adjust this parameter to reproduce on average the experimental value of the pairing gap obtained in the BCS+PVC calculation, where the averaging is weighted with the orbital degeneracy 2j + 1. The experimental value of the neutron pairing gap in 68 Ni was extracted from the mass tables of Ref. [48] with the aid of the commonly used three-point formula [49]. Such calculation scheme allows for understanding the role of the dynamical PVC effects in the formation of the pairing gap. One can see from Fig. 1 that the dynamical PVC is responsible for more than 50 % of the pairing gap value. Its contribution is slightly above 50 % in the peripheral energy regions with respect to the Fermi energy (FE) and increases to 60-70% for the states close to the FE, where the pairing gap values exhibit a smooth maximum. The presence of such maximum is attributed to the functional form of the dynamical kernel K (r) , namely its propagator structure. This result is qualitatively consistent with our previous work [31], where a different calculations scheme was employed, as pointed out above, and with the results of Ref. [16,50]. As in the case of nuclear matter [4][5][6][7][8][9][10], the dynamical kernel is sensitive to details of the approximation made and to the calculation scheme. In the present implementation the experimental value of the pairing gap at T = 0 and low temperatures T ≤ 0.4 MeV is best reproduced in BCS+PVC approach at G = 12.6 MeV, so that with this parameter value ∆ n ≈ 1.6 MeV. The latter value of G is adopted for the BCS+PVC calculations at higher temperatures, which are presented in Figs. 2 and 3 and discussed below. Fig. 2 illustrates the temperature dependence of the neutron pairing gap in 68 Ni in the BCS+PVC approach in comparison with the pure BCS model. For the latter case the parameter G was increased accordingly to reproduce the experimental pairing gap, and the calculated pairing gap is almost state-independent except for the energy window border, where it smoothly decreases to the zero value because of the "soft pairing window" implied in the monopole forces [47]. In this way, it is possible to illuminate the difference in the temperature evolution between the descriptions with only the static and with both static and dynamical kernels. In the conventional BCS, where only the static kernel is taken into account, the temperature dependence is fully determined by the factor 1 − 2f 2 (T ) in Eq. (42), while the static kernel itself has no explicit temperature dependence. Implicitly, this kernel depends on temperature, if its matrix elements are computed in a self-consistent cycle, but this dependence is relatively weak. Some additional temperature dependence may originate from the two-body density if the static kernel is calculated microscopically via Eq. (18). Otherwise, from the transformation (29) it follows that the explicit temperature dependence of the interaction kernel is the consequence of the time dependence, i.e., of the retardation effects present in the dynamical kernel components (24)(25)(26)(27). For this study, the pairing gaps were calculated at temperatures 0 ≤ T ≤ 1. BCS case the pairing gap decreases quickly with the temperature increase: at T = 1.2 MeV it already disappears, while in the BCS+PVC calculation the gap values are still quite sizable grouping around 1 MeV. Another observation is that the pairing gap retains its peaked character in the BCS+PVC approach even when its average value decreases with the temperature increase.
The critical temperature in the theory of superfluidity is known as the temperature, at which the pairing gap vanishes. The canonical BCS relationship between the critical temperature and the value of the pairing gap at T = 0 ∆ 0 is T c ≈ 0.6∆ 0 . In our BCS calculation, the pairing gap disappears at the temperature below ≈ 1.1 MeV, so that the coefficient between ∆ 0 and T c is close to the canonical value. In the BSC+PVC approach with the additional temperature dependence of the dynamical kernel one could expect a different ratio between the ∆ 0 and T c values and also a variation of this ratio from state to state. These trends are illustrated in Fig. 3, where we display the pairing gap as a function of temperature for selected neutron states in 68 Ni, which was obtained in the BCS+PVC calculations. Namely, we show this function for the examples of states near and far from the Fermi energy. The average pairing gap is also shown, and these results are compared to the average BCS pairing gap. This representation helps determine more accurately the values of the critical temperatures for all the cases. The first observation from Fig. 3 is that in the BCS+PVC calculations the pairing gaps for all the states collapse at the same critical temperature. In the approximation described above we get the value T c ≈ 1.4 MeV for all the states, independently on the pairing gap values for these states at T = 0. Another observation is that the critical temperature in the approach with the dynamical kernel is higher than the BCS critical temperature at the same values of the T = 0 average pairing gaps. This indicates that the retardation effects in the dynamical component of the in-medium nucleon-nucleon interaction can  help superfluidity survive at higher temperatures than it is expected in simpler models with only the static kernels.
Figs. 4 and 5 display the analogous calculations for two calcium isotopes, 44 Ca and 46 Ca. The upper panels of Fig.  4 show the results for the pairing gaps in these two nuclei obtained within BCS+PVC and pure BCS approximation at T = 0.2 MeV, which are nearly identical to T = 0 results. Both calculations are performed with the same static pairing strength (only slightly different for 44 Ca and 46 Ca), in order to isolate the PVC effects. In contrast to the case of 68 Ni, in calcium isotopes the PVC produces peaks in the pairing gaps mainly around the Fermi energy while affecting very little the pairing gaps of the peripheral states. As in Ref. [31], this may be due to a stronger cancellation between the self-energy and exchange PVC terms for the peripheral states in these nuclei. One can also notice some irregularities in the BCS+PVC pairing gaps in 46 Ca, namely the small gaps at 1d 3/2 and 2s 1/2 . The remaining variance with Ref. [31] pertains to the differences in the calculation schemes, that is discussed below. The lower panels show the temperature dependence of the BCS+PVC pairing gaps for the states closest to the Fermi energy and their weighted averages in comparison to the behavior of the pure static pairing gaps. The latter are computed with readjusted strength to reproduce the experimental pairing gaps at T = 0. These calculations illustrate the drastic difference in the superfluid phase transition temperature in the two models. Similarly to the case of 68 Ni, in the BCS+PVC calculations the critical temperature is considerably higher than that in the BCS, with approximately the same ratio between the BCS+PVC and BCS values. Thus, this feature of the dynamical kernel remains robust also in calcium isotopes.
In Fig. 5 one can see a more global description of the temperature dependencies of the pairing gaps in 44 Ca and 46 Ca. As in the previous case, the BCS+PVC and pure BCS calculations were performed with different values of the static pairing strength, so that the weighted average pairing gap values reproduce the T = 0 data in both approaches. The general result here is a noticeably slower temperature evolution of the BCS+PVC pairing gaps, as compared to the BCS ones, because of the presence of the non-trivial temperature dependence in the dynamical kernel in the former case. This is observed for all states under study. Remarkably, the shape of the pairing gap as function of energy evolves with temperature becoming less and less peaked in the Fermi energy region with the temperature growth. As briefly mentioned above, the calculation scheme employed in this work is somewhat different from the one of Ref. [31] for T = 0 calculations. There are technical reasons why we do not adopt the same calculation scheme here. In the latter scheme, on the first step, we initiated the procedure by running the RMF+BCS calculations. The BCS equation was solved self-consistently with the RMF in the usual static-kernel approximation with the monopole forces and pairing strength adjusted to the odd-even mass differences, since the phonon vertices are not yet available at this step. On the second step, the obtained pairing gap was subsequently used in the relativistic QRPA (RQRPA), and the phonon vertices and frequencies were extracted. These vertices and frequencies were used on the final step for solving the gap equation in the BCS+PVC approximation with both static and dynamical kernels. Notice that in this scheme, using the RQRPA phonons is only partially self-consistent, because in the final BCS+PVC calculation the strength constant of the static pairing kernel has to be refitted to reproduce the final pairing gap. In addition, the denominators of the dynamical kernel contain the pure single-particle energies (no quasiparticle energies) in combination with RQRPA phonon frequencies, which have the information about superfluid pairing.
To adopt an analogous calculation scheme at finite temperature, we would need to (i) generalize the RQRPA to finite temperature, which is a non-trivial task on its own, (ii) take into account a more complex structure of the phonon vertices, which means (iii) re-deriving the entire approach for the pairing propagator and, thus, for the pairing gap in the thermal RMF+BCS (or Hartree-(Fock)-Bogoliubov) quasiparticle basis. This would be essentially a more complicated and advanced solution at finite temperature, which goes beyond the scope of the present article, but which will be considered in the future.
To avoid such complications, in this work we did not include pairing correlations on the first two steps, and, instead, ran thermal RMF and finite-temperature relativistic RPA to obtain the mean-field and phonon characteristics. On one hand, this reduces the accuracy of the phonon calculations at T = 0, but, on the other hand, allows us to avoid the inconsistency between the denominators of the dynamical kernel, which would contain the energies and frequencies obtained in the quasiparticle picture, and exponential factors with the temperature dependence, which are not yet adopted to superfluid pairing.
These differences in the calculation schemes, which are both approximate and may be replaced by a more accurate one in the future, are responsible for somewhat different behavior of the resulting pairing gaps obtained in this work at T = 0, as compared to Ref. [31]. However, their enhancements for the states surrounding the Fermi energy and the general trends in nickel and calcium nuclei remain similar, although in the Ref. [31] the peaks of the pairing gaps are not always centered at the Fermi surface. This may be because of only partial self-consistency described above, which modifies the location of the poles of the dynamical kernel. Slightly weaker PVC effects in Ref. [31] are observed mainly due to the reduction of the PVC vertices, as they are multiplied with the combinations of the occupation factors in RQRPA [47].

IV. SUMMARY AND OUTLOOK
The equation of motion for the two-time two-fermion correlation function in a strongly coupled many-body system at finite temperature is considered. We show that, as in the case of zero temperature, the EOM for this propagator takes the form of the Dyson-Bethe-Salpeter equation with the interaction kernel, which is split into the static and dynamical components. This kernel includes, in principle, all the in-medium physics derived from the underlying bare two-fermion interaction. While the static component of the kernel depends on the correlated two-fermionic density, the dynamical component contains a higher-rank fermionic propagator. The latter, in the case of the symmetric form of the kernel, is represented by the propagator of four fermions. Factorization of this propagator allows for the truncation of the many-body problem at the level of two-body correlation functions whose EOM's, together with those for the one-fermion correlation function discussed in Ref. [35] form a closed system of integral equations. The equation for the temperature-dependent pairing gap, which is related to the residue of the two-time particle-particle propagator, is formulated as a static limit of the EOM for this propagator. To extract the pairing gap, which is one of the most common characteristics of superfluidity, we adopt the BCS-like variational ansatz for the ground state wave function. The resulting equation allows for an extension of the BCS approximation to correlations of higher complexity, which introduce an additional non-trivial temperature dependence of the pairing gap.
The effects of the dynamical kernel at finite temperature are illustrated in the calculations of the neutron pairing gaps for 68 Ni and two calcium isotopes, 44 Ca and 46 Ca. The pole character of this kernel gives rise to the peak of the pairing gap around the Fermi surface at all temperatures when the pairing gap has non-zero values. We find that the time dependence, mostly the retardation, present in this kernel translates to a different temperature dependence of the pairing gap than the one of the BCS approximation. In particular, the presence of the dynamical term leads to noticeably higher values of the critical temperature. This finding may be important for numerous applications. For instance, in applications to nuclear astrophysics, such as the r-process nucleosynthesis in the neutron star mergers and the supernovae explosion, the nuclear input for temperatures below 1-2 MeV is involved. Crossing the critical temperature, i.e., the superfluid phase transition can affect considerably the excitation spectra snd, thus, various reaction rates. An example of the electron capture is discussed in Ref. [51]. Therefore, the form of the dynamical kernel and its temperature dependence has to be computed as accurate as possible. Self-consistent calculations of this kernel as well as the static one, desirably in an ab-initio framework, thus, remain an important topic for future research.