text:
# Dissociation of Diatomic Molecule by Energy-Feedback Control

## 1 Introduction

During last decade a growing interest has been observed in the control problems for molecular systems in classical and quantum formulation . One of the simplest problems of that type is dissociation problem for diatomic molecule . It often serves as a benchmark for comparison of classical and quantum calculations. In the paper  possibilities of dissociation a molecule by monochromatic (single frequency) laser field have been explored for the case of hydrogen fluoride (HF) molecule using Chirikov’s resonance overlap criterion. In  the case of two-frequency (two-laser) control field was investigated. It was shown that intensity of two-frequency field required for dissociation can be reduced compared to one-frequency case. In  the further reduction of the control field intensity has been demonstrated by means of chirping (frequency modulation) the laser frequency with constant chirping rate.

New possibilities for changing of physical and chemical properties are provided by using feedback. In  methods based on geometric control theory (inverse control) were proposed for molecular systems. In the papers  a general method for investigation of physical systems by means of feedback controlling actions was proposed. It was shown in  by example of the problem of escape from a potential well that feedback allows to reduce control intensity required for overcoming a potential barrier by several orders of magnitude. It appears that using speed-gradient method  for design of feedback control action allows to achieve desired level of energy by excitation of intensity inversely proportional to the system dissipation degree. Particularly, for conservative systems it allows to reach any energy level by (ideally) arbitrarily small control action.

In this paper the algorithm of the feedback frequency chirping of the controlling field frequency is proposed based on dependence of the natural frequency of the nonlinear molecule oscillations on its energy, further referred to as frequency-energy characteristics or frequency-energy curve (FE-curve). On-line implementation of the algorithm requires measurement of molecule energy only. This can be easily generalized for an ensemble of molecules, replacing exact energy of the separate molecule in the control criterion by the average energy per single molecule. However, for control algorithm design the knowledge of the frequency-energy relation is required. In this paper a classical Morse potential is considered, which allows obtaining this relation in a simple analytical form. For other potentials, where the analytical solution is too complicated or impossible, the frequency-energy relation can be obtained numerically. Efficiency of the proposed algorithm for the case of HF molecule is illustrated by computer simulation. It is shown that the time of dissociation can be reduced significantly when using proper feedback excitation.

## 2 Basic equations

Let us consider a diatomic molecular system under the action of the external laser field. Dynamics of such a system can be described by the following controlled Hamiltonian

$$H=\frac{p^2}{2m}+\mathrm{\Pi }(r)\mu (r)u(t),$$
(1)
where $`m`$ is reduced mass, $`p`$ is momentum, $`\mathrm{\Pi }(r)`$ is potential of interatomic interaction, $`\mu (r)`$ is dipole moment of the molecule, $`u(t)`$ is intensity of external field. The value $`u(t)`$ serves as control variable. E.g. for monochromatic control field $`u(t)=E\mathrm{cos}(\omega t)`$, where $`E`$ and $`\omega `$ are strength and frequency of the external laser field. Substitution of (1) into Hamilton equations

$$\frac{r}{t}=\frac{H}{p},\frac{p}{t}=\frac{H}{r}$$
(2)
yields the following equation of molecular motion

$$m\ddot{r}=\mathrm{\Pi }^{}(r)+\mu ^{}(r)u(t).$$
(3)
For the case of harmonic excitation equations of motion is as follows:

$$m\ddot{r}=\mathrm{\Pi }^{}(r)+E\mu ^{}(r)\mathrm{cos}(\omega t).$$
(4)
Let us use Morse form of the interatomic potential

$$\mathrm{\Pi }(r)=D\left(1e^{\alpha (ra)}\right)^2D=D\left(e^{2\alpha (ra)}2e^{\alpha (ra)}\right),$$
(5)
where $`D`$ is the bond energy, $`a`$ is the equilibrium interatomic distance. The corresponding force is

$$f(r)=\mathrm{\Pi }^{}(r)=2\alpha D\left(e^{2\alpha (ra)}e^{\alpha (ra)}\right).$$
(6)
The dipole moment can be represented in the form

$$\mu (r)=Are^{\xi r^4}\mu ^{}(r)=A\left(14\xi r^4\right)e^{\xi r^4},$$
(7)
where $`A`$ and $`\xi `$ are the known constants. Thus the equation of motion 4 reads

$$m\ddot{r}=2\alpha D\left(e^{2\alpha (ra)}e^{\alpha (ra)}\right)+EA\left(14\xi r^4\right)e^{\xi r^4}u(t).$$
(8)
For sufficiently small $`\xi `$ the simplified expression for dipol potential can be used:

$$\mu (r)=Ar,\mu ^{}(r)=A,$$
(9)
In the vicinity of equilibrium $`ra`$ equation (8) with harmonic excitation can be reduced to the equation of linear forced oscillations

$$m\ddot{r}=Cr+E\mu ^{}(a)\mathrm{cos}(\omega t),C\stackrel{\text{def}}{=}\mathrm{\Pi }^{\prime \prime }(a)=2\alpha ^2D,$$
(10)
where $`C`$ is linear stiffness of the bond. The natural frequency of the linear system is

$$\mathrm{\Omega }_0\stackrel{\text{def}}{=}\sqrt{C/m}=\alpha \sqrt{2D/m}.$$
(11)
## 3 Control of the dissociation process

In a nonlinear system the resonance frequency is function of the amplitude of oscillations, or, in other words, it depends on the energy of the oscillator. Thus the faster dissociation can be achieved if the excitation frequency is changing (decreasing in our case) while the energy increases. The simplest control algorithm of this kind can be done by linear chirping:

$$\omega =\mathrm{\Omega }_0\epsilon t$$
(12)
where $`\epsilon `$ is a constant (chirping rate), which characterizes the speed of the frequency decrease. In this case equation of motion (4) takes the form

$$m\ddot{r}=\mathrm{\Pi }^{}(r)+E\mu ^{}(r)\mathrm{cos}\phi ,\dot{\phi }=\omega (t)=\mathrm{\Omega }_0\epsilon t.$$
(13)
(Note that in this case the representation $`\phi =\omega t`$ is not valid since the frequency $`\omega `$ is not constant.

More effective dissociation can be achieved by taking advantage of the feedback control. In this case the frequency of excitation can be changed according to the oscillator energy so that at any amplitude of oscillations the excitation will act at the resonance frequency. Let $`\mathrm{\Omega }(W)`$ be the natural frequency of the diatomic molecular oscillator at the specified energy $`W`$. We call the function $`\mathrm{\Omega }(W)`$ frequency-energy curve of FE-curve. It contains important information about dynamics of the molecular motion. The key idea of our approach is that if the curve $`\mathrm{\Omega }(W)`$ is known and the energy $`W=W(t)`$ is observable then control algorithm

$$m\ddot{r}=\mathrm{\Pi }^{}(r)+E\mu ^{}(r)\mathrm{cos}\phi ,\dot{\phi }=\mathrm{\Omega }(W)$$
(14)
will give the desired excitation at the resonance frequency at any time instant.

The relation $`\mathrm{\Omega }(W)`$ in principle can be calculated for any known potential $`\mathrm{\Pi }(r)`$ using the integral of energy, which is realized when the external force is absent

$$\frac{1}{2}m\dot{r}^2+\mathrm{\Pi }(r)=WT=\sqrt{2m}\underset{r_1}{\overset{r_2}{}}\frac{dr}{\sqrt{W\mathrm{\Pi }(r)}},$$
(15)
where $`W`$ is a constant value of energy, $`T`$ is period of the oscillations, $`r_1`$ and $`r_2`$ are the minimum and the maximum possible values of radius $`r`$ for the specified potential energy (the solutions of the equation $`\mathrm{\Pi }(r)=W`$). If period $`T`$ is known then frequency can be obtained as $`\mathrm{\Omega }=2\pi /T`$. Obviously integral (15) can not be calculated analytically for arbitrary potential and the desired relation $`\mathrm{\Omega }(W)`$ should be evaluated by numeric integration.

Luckily for the Morse potential (5) the integral can be evaluated analytically in the following way. Substitution $`\zeta =\zeta (r)=e^{\alpha (ra)}`$ transforms the integral (15) in the case of Morse potential to the simple form which can be easily integrated

$$T=\frac{1}{\alpha }\sqrt{\frac{2m}{D}}\underset{\zeta _1}{\overset{\zeta _2}{}}\frac{d\zeta }{\sqrt{\frac{W}{D}\zeta ^2+2\zeta 1}}=\frac{\pi }{\alpha }\sqrt{\frac{2m}{W}}.$$
(16)
The above formula is valid if $`W<0`$. Otherwise, the motion of the system is not periodic, since the molecule dissociates. Equation (16) provides the following simple relation between energy and natural frequency of Morse oscillator

$$\frac{W}{D}=\left(\frac{T_0}{T}\right)^2=\left(\frac{\mathrm{\Omega }}{\mathrm{\Omega }_0}\right)^2,$$
(17)
where $`\mathrm{\Omega }_0`$ is defined by (11), $`T_0=2\pi /\mathrm{\Omega }_0`$. The control law (14) therefore takes the form

$$\dot{\phi }=\mathrm{\Omega }(W)=\mathrm{\Omega }_0\sqrt{\frac{W(t)}{D}}.$$
(18)
where $`W(t)`$ is the current value of the molecule energy.

## 4 Calculations

Efficiency of the proposed control algorithm was confirmed by computer simulations for dissociation of the hydrogen fluoride (HF) molecule that has become a benchmark problem for controlled dissociation . The parameters of the model 8, corresponding to HF molecule  are as follows: $`m=1732`$, $`D=0.2101`$, $`\alpha =1.22`$, $`a=1.75`$, $`A=0.4541`$, $`\xi =0.0064`$, $`E=0.1`$. All quantities are given in atomic units. The initial condition for excitation is molecule equilibrium at the bottom of the potential well: $`r=a`$, $`\dot{r}=0`$. The equation was integrated numerically by central differences method with time step $`dt=T_0/160`$ (integration with smaller time steps yields practically the same results.

The result of excitation at the resonance frequency is shown in Fig. 1a. The graph shows that the motion is quasiperiodic and maximum value of energy, which can be achieved by this method, is $`0.68D`$. Since the natural frequency of the molecule is decreasing when the energy is growing, more efficient excitation can be achieved at lower frequencies. It is even possible to dissociate the molecule if the excitation frequency is close to $`0.841\mathrm{\Omega }_0`$ — see Fig. 2a. The figure shows that at $`t13.5T_0`$ the energy crosses the zero level, and at $`t=14.5T_0`$ the energy becomes a positive constant. This practically means that the potential energy is negligible and atoms are moving in the opposite directions with constant velocities (total dissociation). However this result is not structurally stable: small variation in the excitation frequency can noticeably increase the dissociation time and even break the dissociation.

Fig. 2a,b correspond to linear frequency chirping (12) at two different chirping rates. The dissociation is also achieved, but it is slower then in the previous case and the results are very sensitive to the variations of $`\epsilon `$. Additional flexibility to the chirping process can added by using the relation $`\omega =\omega _0\epsilon t`$, where $`\omega _0`$ is a constant less then $`\mathrm{\Omega }_0`$. Optimal selection of both parameters is a pretty complicated task. In more detail this problem is investigated in .

In this paper we will concentrate on the feedback chirping, based on the frequency-energy curve (17). Fig. 3a shows the energy as a function of time for the feedback control (18). At the initial stage the energy is growing faster then in other cases, but after the level about $`W=0.5D`$ is reached, the efficiency of the excitation decreases so that dissociation is not realized. Probably this phenomenon is caused by wrong phasing which prevents the molecule dissociation although the resonance condition is all the time fulfilled. Indeed, the sign of the multiplier $`\mu ^{}(r)`$ in equation (4) due to formula (7) is changing for high energies, so that the excitation starts acting in the opposite phase. To overcome this difficulty the following algorithm can be used: when the energy of the molecule reaches some certain value $`W_{}`$ then the sign of the excitation force is permanently changing to the opposite. The result of such an excitation is shown in Fig. 3b. The dissociation is achieved at $`t=10T_0`$, which is faster then in all previous cases.

Thus for the selected level of control ($`E=0.1`$ a.u.) the best result is obtained for feedback excitation. However excitation with shifted constant frequency gives the competitive result. The situation changes if we consider lower level of excitation. If this level is sufficiently small then excitation at any constant frequency can not provide the dissociation. On the contrary, excitation with changing frequency can in principle give dissociation for any arbitrarily small level of control. For illustration let us consider the 5 times smaller control then in previous calculations: $`E=0.02`$ a.u. In this case dissociation at a constant frequency is not possible. Results of linear chirping and feedback control excitation are given in Fig. 4. For simplicity these calculations use approximation $`\mu ^{}(r)\mu ^{}(a)`$. The plots show that both methods provide dissociation. However, linear chirping achieves dissociation at $`t=89T_0`$, while energy-feedback control achieves it at $`t=34T_0`$. Thus the feedback control provides more then two times faster dissociation then the best variant of the linear chirping.

Although the proposed algorithm requires measurement of molecule energy only, its on-line implementation requires changing controlling action at the time intervals less then $`10^{12}`$ s. This can be achieved with modern ultrafast lasers. The approach can be easily generalized to control of an ensemble of molecules, replacing exact energy of the separate molecule in the control objective by the average energy per single molecule.


format: text

text:
# Constraints on reionization from the thermal history of the intergalactic mediumBased on observations made at the W.M. Keck Observatory which is operated as a scientific partnership between the California Institute of Technology and the University of California; it was made possible by the generous support of the W.M. Keck Foundation. Based on public data released from the VLT/UVES Commissioning and Science Verification and from the OPC program 65.O-296A (PI S.-D’Odorico) at the VLT/Kueyen telescope, ESO, Paranal, Chile.

## 1. Introduction

Neutral hydrogen in the intergalactic medium (IGM) along the line of sight to quasars at redshifts $`z6`$ produces hundreds of Ly$`\alpha `$ absorption lines. The fact that not all flux is absorbed (i.e., the absence of a ‘Gunn-Peterson’ trough, Gunn & Peterson 1965) requires that the universe be ionized to a level far higher than can be attributed to residual ionization from recombination. At lower redshifts, $`z3`$, observed stars and quasars produce enough ionizing photons to explain the high levels of ionization, but the nature of the sources responsible for converting most of the IGM from neutral to ionized remain uncertain, as does the epoch of reionization (e.g., Barkana & Loeb 2001).

The observed mean flux decrement $`D_A`$ blueward of the quasar’s Ly$`\alpha `$ emission line increases with the redshift of the quasar, both because the intensity of the ionizing background radiation decreases above $`z=4`$ (e.g., McDonald & Miralda-Escudé 2001) and because the mean density $`\overline{\rho }`$ of the universe – and hence the neutral fraction $`x\overline{\rho }T^{0.7}/\mathrm{\Gamma }_{\mathrm{H}\mathrm{I}}`$ for fixed values of the photo-ionization rate $`\mathrm{\Gamma }_{\mathrm{H}\mathrm{I}}`$ and temperature $`T`$ – increases. Recently, Becker et al. (2001) and Djorgovski et al. (2001) observed a sudden increase in $`D_A`$ in the spectra of redshift $`z6`$ quasars discovered by the SLOAN digital sky survey. Such a sharp rise has been predicted to mark the transition associated with a sudden epoch of reionization (e.g., Cen & Ostriker 1993; Gnedin 2000). Similarly, a sudden increase in the He II opacity has been detected around $`z3`$ (Reimers et al. 1997; Heap et al. 2000; Kriss et al. 2001), associated with helium reionization.

Another way to study the reionization history of the IGM is to investigate its thermal evolution. Because its cooling time is long, the low-density IGM retains some memory of when and how it was reionized (e.g., Miralda-Escudé & Rees 1994; Hui & Gnedin 1997; Haehnelt & Steinmetz 1998). The combined effects of photo-ionization heating and adiabatic expansion introduce a tight temperature-density relation in the unshocked IGM, which can be approximated by a power-law $`T=T_0(\rho /\overline{\rho })^{\gamma 1}`$ for densities around the cosmic mean (Hui & Gnedin 1997). A change in these parameters influences the shapes of the Ly$`\alpha `$ lines, because thermal broadening and Jeans smoothing determine the line widths (e.g., Theuns, Schaye & Haehnelt 2000). Schaye et al. (1999) used hydrodynamical simulations to demonstrate that one can accurately calibrate the relation between the minimum line width ($`b`$) as a function of column density ($`N_{\mathrm{H}\mathrm{I}}`$) on the one hand, and the underlying $`T\rho `$ relation on the other. Schaye et al. (2000) applied the method to observations in the redshift range 2.0-4.5 and found that $`T_0`$ peaks at $`z3`$, which they interpreted as evidence for the reionization of He II. Ricotti, Gnedin & Shull (2000) used pseudo hydrodynamical simulations and found a similar temperature increase, albeit only at the 0.5$`\sigma `$ level. McDonald et al. (2001) found no evidence for temperature evolution, but their analysis neglected the important temperature dependence of Jeans smoothing. Finally, the analysis of Zaldarriaga et al. (2001) neglected hydrodynamical effects all together.

Here, we provide new evidence for a relatively sudden increase in $`T_0`$ between redshifts $`z3.5`$ and 3.0, using a new method based on a wavelet decomposition of the absorption spectrum. We then use the measured values of $`T_0`$ at higher redshift to constrain the epoch of hydrogen reionization $`z_\mathrm{H}`$ and find that the data require $`z_\mathrm{H}9`$ for any reasonable value ($`T_06\times 10^4\mathrm{K}`$) of the hydrogen reionization temperature.

## 2. Helium reionization

### 2.1. Wavelet analysis

Our analysis uses wavelets to characterize line-widths. By picking an appropriate wavelet scale (we used $`15`$ km s-1), we find that the amplitude $`A`$ of the wavelet anti-correlates with the widths of the lines – and hence the temperature $`T_0`$ of the absorbing gas, $`AT_0^1`$. By examining the statistics of $`A`$ along the spectrum, we can look for changes in temperature in an objective way, given that the wavelet decomposition is unique. Full details can be found in Theuns & Zaroubi (2000) and Theuns et. al (2002), here we give only a brief summary of the underlying reasoning.

To investigate whether a region of given size $`V`$ of the QSO spectrum has an unusual temperature, we compare the cumulative probability distribution $`𝒞_V(A)`$ of the wavelet amplitudes in that region, with $`𝒞(A)`$ for the spectrum as a whole. If the region is unusually hot, it will tend to have very few large wavelet amplitudes, and hence the maximum difference $`\mathrm{\Delta }(V)=𝒞_V𝒞`$ will be large. (Note that by construction, $`|\mathrm{\Delta }|1`$, and it is defined for a given region). Conversely, cold regions will have large, negative $`\mathrm{\Delta }(V)`$. In panel (a) of Fig. 1, we have plotted $`\mathrm{\Delta }`$ for regions of size 5000 km s-1 along a mock spectrum of a simulation111The cosmological parameters for this vacuum-energy dominated, flat, cold dark matter, smoothed particle hydrodynamics (SPH) simulation are $`(\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda },\mathrm{\Omega }_bh^2,h,\sigma _8)=(0.3,0.7,0.019,0.65,0.9)`$. The simulation box is $`12h^1`$ Mpc on a side, and gas and dark matter are represented with 2563 particles each.. The mock spectrum has a jump in $`T_0`$ from $`T_0=2.2\times 10^4`$ K below redshift $`z=3.3`$ to $`T_0=1.5\times 10^4`$ K above it. As expected, there is a corresponding jump in $`\mathrm{\Delta }`$ at $`z=3.3`$, from positive values at low $`z`$ to negative ones at higher $`z`$.

To determine whether a high value of $`|\mathrm{\Delta }|`$ is statistically significant – and hence whether we have identified a region with an unusual temperature – we proceed as follows. We repeat the procedure with spectra where the absorption lines are randomly scrambled (with replacement, see Theuns et al. 2002), in order to wash-out any correlations in the wavelet amplitudes, resulting from intrinsic temperature fluctuations. Using these randomized spectra, we can construct the statistical probability $`P(\mathrm{\Delta })`$ from the fraction of regions in the randomized spectra, that have a given value of $`\mathrm{\Delta }`$. Given $`P(\mathrm{\Delta })`$, we can determine how likely a value of $`\mathrm{\Delta }`$ – and hence of a temperature deviation – in the original spectrum is. Performing this analysis, we find that the detected change in $`\mathrm{\Delta }`$ in the mock spectrum has a statistical significance of more than 99.5 per cent, shown as the full line in panel (a) (positive values refer to regions hotter than average, and vice versa for negative values). This means that in only one of out 200 random realisation, do we, by chance, get values of $`|\mathrm{\Delta }|0.2`$. Note that we only use the simulation to generate the spectrum, and not to assign the statistical significance of a change in $`\mathrm{\Delta }`$. In the preparation of the mock spectrum, we have imposed the same biases as are present in the real data, by adding noise and instrumental broadening to the simulated lines, and by scaling the mean absorption to the observed value. We therefore believe that the method can be applied to real data with confidence.

We have applied the same wavelet analysis to high-resolution echelle spectra of quasars Q0055-2169 (emission redshift $`z_{\mathrm{em}}=3.6`$, Kim et al 2001), the combined spectrum of Q0302-003 ($`z=[3,3.27]`$, Kim et al 2001) and APM 0827+5255 ($`z=[3.27,3.7]`$, Ellison et al 1999) and Q1422+231 ($`z_{\mathrm{em}}=3.6`$, Rauch et al 1997) (panels b, c and d of Fig. 1 respectively). In each spectrum, we find a cold region at high redshift, and a hot region at lower redshift, each significant at the more than 99 percent level when compared to randomized spectra. Panels (b) and (c) appear very similar to the mock spectrum of panel (a), where we had imposed a sudden temperature increase below $`z=3.3`$. Note that the implied temperature evolution is exactly opposite of what one would expect from photo-heating in the optically thin limit, in which case the IGM will gradually cool down. We take this as strong evidence that a large fraction of He II is reionized around redshift $`z3.3`$.

The mean wavelet amplitude scales approximately inversely with temperature: $`AT_0^1`$ (Theuns et al. 2002). We have used hydrodynamical simulations to calibrate the proportionality constant, and investigate its dependence on $`\gamma `$. Applying the calibration to the above QSO sample, we obtain values for $`T_0`$ in good agreement with those obtained by Schaye et al (2000) from the cutoff in the line width – column density relation, i.e., $`T_010^{4.1}`$ for $`z3.4`$, and $`T_010^{4.3}`$ at $`z3.0`$. This assumes the value of $`\gamma `$ as determined by Schaye et al., but reasonable changes of $`\gamma `$ do not change $`T_0`$ by more than $`10`$ per cent (Theuns et al. 2002). Simulations that include radiative transfer are required to investigate whether such a temperature change is consistent with He II reionization.

Other heating mechanisms, for example shock heating by galactic winds, do not have a major influence on the value of $`T_0`$ deduced from fitting the cutoff in the $`bN_{\mathrm{H}\mathrm{I}}`$ diagram, at least as long as the volume fraction of shocked gas remains small. This is because the method is based on identifying the narrowest lines in a region, irrespective of whether there is also a set of much broader lines. In contrast, the wavelet method used here examines all lines in a stretch of spectrum. So the fact that both methods find similar values for $`T_0`$ suggests that photo-heating is indeed the dominant heating mechanism, and that the volume        Temperature evolution of the IGM after $`z=3.4`$. The points with error bars are determined from fitting the cutoff in the $`bN_{\mathrm{H}\mathrm{I}}`$ relation (Schaye et al 2000). The curves emanating from the star all assume that $`T_0=10^{4.4}\mathrm{K}`$ at $`z=3.4`$ as a result of He II reionization, and differ in the assumed UV-background at lower redshifts. The solid line includes photo-heating by H I, He I and He II. The short dashed curve neglects heating from He II, the long dashed curve neglects heating from both H I and He I, and finally the dot-dashed curve ignores photo-heating all together. All models fit the inferred temperature evolution extremely well.

filling factor of gas that has been shocked by winds is small.

### 2.2. Thermal evolution

After reionization, the evolution of $`T_0`$ is given by

$$\frac{1}{T_0}\frac{dT_0}{dt}\frac{1}{\mu }\frac{d\mu }{dt}=2H+\frac{\mu \mathrm{\Delta }_ϵ}{\frac{3}{2}k_BT_0},$$
(1)
where $`H`$ is the Hubble parameter, $`k_B`$ Boltzmann’s constant, $`\mu `$ the mean molecular weight, and $`\mathrm{\Delta }_ϵ`$ is the effective radiative cooling rate \[in $`\mathrm{erg}\mathrm{g}^1\mathrm{s}^1`$\]. $`\mathrm{\Delta }_ϵ`$ is negative (positive) for net cooling (heating) and includes photo-electric heating and cooling via recombination, excitation, inverse Compton scattering, collisional ionization, and bremsstrahlung (we use the rates listed in Appendix B of Theuns et al. 1998). For gas around the mean density, the dominant cooling process is the adiabatic expansion of the universe (the first term on the right-hand side of equation 1), except at $`z>7`$ where inverse Compton cooling off the cosmic microwave background is more efficient. The radiative heating and cooling rates depend on the ionization balance of the gas, which itself depends on the temperature. By coupling equation 1 to the differential equations for the ionization balance (listed in Appendix B of Theuns et al. 1998), we can thus solve for the evolution of $`T_0`$ given a model for the ionizing background, an initial value for $`T_0`$, and the initial ionization state. We have tested this procedure using full hydrodynamic simulations and find that the evolution of $`T_0`$ is reproduced very well.

We assume the universe to be permeated by a uniform UV-background with a power-law spectral shape,

$$J=\{\begin{array}{cc}J_{\mathrm{H}\mathrm{I}}\left(\frac{\nu }{\nu _{\mathrm{H}\mathrm{I}}}\right)^{1.8}\hfill & \nu <\nu _{\mathrm{He}\mathrm{II}}\hfill \\ J_{\mathrm{He}\mathrm{II}}\left(\frac{\nu }{\nu _{\mathrm{H}\mathrm{I}}}\right)^{1.8}\hfill & \nu \nu _{\mathrm{He}\mathrm{II}}\hfill \end{array}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2\mathrm{sr}^1\mathrm{Hz}^1,$$

where $`\nu _{\mathrm{H}\mathrm{I}}`$ is the H I ionization threshold. (Note that we allow for a step between the intensities of H I and He II ionizing photons, but normalize $`J_{\mathrm{He}\mathrm{II}}`$ at the hydrogen ionization edge.) We set $`J_{\mathrm{H}\mathrm{I}}=4\times 10^{23}`$, which yields a photo-ionization rate of $`\mathrm{\Gamma }_{\mathrm{H}\mathrm{I}}=10^{13}\mathrm{s}^1`$ and vary $`J_{\mathrm{He}\mathrm{II}}`$ (and the corresponding ionization rate $`\mathrm{\Gamma }_{\mathrm{He}\mathrm{II}}`$) as described below. Because the photo-heating rate is independent of the amplitude of the UV-background as long as the gas is highly ionized, the exact value of $`J_{\mathrm{H}\mathrm{I}}`$ after hydrogen reionization is unimportant. Since He II is not always highly ionized, the thermal evolution does depend on $`J_{\mathrm{He}\mathrm{II}}`$.

The predicted evolution of $`T_0`$ is compared with the data from Schaye et al (2000) in Figure 2. We assume the IGM to be in ionization equilibrium at temperature $`T_0=10^{4.4}\mathrm{K}`$ at redshift $`z=3.4`$ as a result of He II reionization. The curves emanating from the big star then show the subsequent evolution of $`T_0`$, for the following imposed UV-backgrounds. The solid curve is for a constant ionizing rate of $`\mathrm{\Gamma }_{\mathrm{H}\mathrm{I}}=\mathrm{\Gamma }_{\mathrm{He}\mathrm{II}}=10^{13}\mathrm{s}^1`$, and the corresponding photo-heating from all species. Increasing or decreasing $`\mathrm{\Gamma }_{\mathrm{H}\mathrm{I}}`$ by a factor 10 does not change the evolution appreciably. The short and long dashed curves ignore photo-electric heating of He II and of both H I and He I respectively. Finally, the dot-dashed line ignores photo-heating all together. The measurements of Schaye et al. (2000) below $`z3`$ are clearly consistent with a thermal evolution dominated by adiabatic cooling.

Prior to He II reionization, some regions will already be ionized in He II by local sources. If most of the universe is reionized significantly later, then such differences in reionization epoch will lead to spatial variations in $`T_0`$. The wavelet analysis by Theuns et al. (2002) can detect variations in $`T_0`$ of order 50 per cent, over a region of size 5000 km s-1, yet no such fluctuations were found in the data. The wavelet analysis by Zaldarriage (2001) also failed to detect any such temperature fluctuations. This suggests that the temperature increase is the result of the overall increase in the far UV-background following the percolation of He III regions, which prior to reionization are too small to be detected by current methods.

## 3. Hydrogen reionization

Given the success in reproducing the $`z<3.4`$ temperature evolution, we now turn to higher redshifts. Abel & Haehnelt (1999) performed radiative transfer calculations of the expansion of an ionization bubble around a QSO and found typical post-reionization temperatures of $`4\times 10^4`$ K, including He II reionization. To put a conservative upper limit on the hydrogen reionization redshift $`z_{\mathrm{H}\mathrm{I}}`$, we assume a higher value of $`T_0=6\times 10^4\mathrm{K}`$. Assuming the gas to be in ionization equilibrium just after reionization, we compute $`T_0(z)`$ for a given $`z_{\mathrm{H}\mathrm{I}}`$ as before. In Fig. 3 the predicted thermal evolution for various values of $`z_{\mathrm{H}\mathrm{I}}`$ is compared with the data.

The temperature decreases rapidly with decreasing redshift because at these temperatures photo-heating cannot compensate for the rapid expansion cooling. Consequently, the temperature drops quickly below the values measured at $`z4`$, unless hydrogen reionization occurred relatively recently. The solid curves are for $`\mathrm{\Gamma }_{\mathrm{H}\mathrm{I}}=10^{13}`$ s-1 and a range of reionization redshifts, as indicated in the figure. For $`z_\mathrm{H}=6.2`$, we show $`\mathrm{\Gamma }_{\mathrm{H}\mathrm{I}}=10^{14}\mathrm{s}^1`$ for comparison. The models with late H I reionization, $`z_\mathrm{H}8`$, fit the data best. Models with $`z_\mathrm{H}9`$ have $`z4`$ temperatures that are significantly below the measured values. For example, the $`z_\mathrm{H}=9.2`$ curve has a reduced $`\chi ^2`$ of 4.5 for the 5 $`z>3.4`$ data points.

The contribution of He II photo-heating at $`z5`$ is uncertain. It is likely that stars dominate the H I ionizing background at such high redshifts (e.g., Madau, Haardt, & Rees 1999), and stars emit very few He II ionizing photons. However, very massive and extremely metal poor stars could lead to a non-negligible He II ionizing background (Tumlinson & Shull 2000). The QSO contribution is also uncertain, although the paucity of faint point sources in the Hubble Deep Field does provide some constraints (Haiman, Madau & Loeb 1999). If He II reionized at $`z3`$, as argued in the previous sections, then the He III abundance, $`x_{\mathrm{He}\mathrm{III}}`$, should be small at higher $`z`$. We therefore set $`\mathrm{\Gamma }_{\mathrm{He}\mathrm{II}}=10^{15}\mathrm{s}^1`$ so that $`x_{\mathrm{He}\mathrm{III}}0.1`$ at redshift $`z=3.5`$. This limits the plausible contribution from He II photo-heating, and allows us to put a conservative upper limit on the redshift of H I reionization, $`z_\mathrm{H}<9`$.

The importance of He II heating is illustrated by the solid curve connecting the filled squares. This model has $`z_\mathrm{H}=10.2`$, but the He II ionizing background is increased such that $`x_{\mathrm{He}\mathrm{III}}0.6`$ at $`z=3.5`$. This increases $`T_0(z=4)`$ significantly, although it still falls below the measured values. Finally, the long dashed line is for a model with $`z_\mathrm{H}=20`$ and $`J_{\mathrm{He}\mathrm{II}}`$ further increased so that $`x_{\mathrm{He}\mathrm{III}}0.96`$ at $`z=3.5`$. The temperature of this model is consistent with the data, yet He II is ionized at the more than 90 percent level as early as $`z=5`$. Such a high level of ionization conflicts with the evidence that He II reionizes at $`z3`$ from the observed He II opacities and the associated increase in $`T_0`$.

All this leads us to the following conclusions. Two independent methods consistently find a rather sudden increase in the temperature of the IGM over the range $`z`$ 3.5–3.0, which we associate with He II reionization. If this interpretation is correct, then the He III fraction must be low at higher redshifts. Therefore, above redshifts 3.5, He II photo-heating cannot be significant and the IGM cools rapidly following H I reionization. The high values of the IGM temperature at $`z4`$ then require that H I reionization occurred late as well, $`z_\mathrm{H}<9`$, for any reasonable value $`T_06\times 10^4`$ K for the H I reionization temperature. More plausible reionization temperatures of $`T_04\times 10^4`$ K and $`2\times 10^4`$ K would constrain the hydrogen reionization redshift further to $`z_\mathrm{H}8`$ and $`z_\mathrm{H}7`$, respectively.

Acknowledgments TT thanks PPARC for the award of an Advanced Fellowship. JS is supported by a grant from the W.M. Keck Foundation. We acknowledge support from the ‘Physics of the Intergalactic Medium’ network set up by the European Commission. Research was conducted in cooperation with Silicon Graphics/Cray Research utilizing the Origin 2000 super computer at DAMTP, Cambridge.


text:
# Boundary Conditions in an Electric Current Contact

## Abstract

In most electronic devices, electric current of both types (electrons and holes) flows through a junction. Usually the boundary conditions have been formulated exclusively for open circuit. The boundary conditions proposed here bypass this limitation by the first time, as far as we are aware. Besides, these new boundary conditions correctly describe current flow in a circuit, i.e., closed circuit conditions, which are the usual operation conditions for electronic devices and for the measurement of many transport properties. We also have generalized the case (as much as it is possible in a classical treatment), so self-consistent boundary conditions to describe current-flow through a contact between two arbitrary conducting media are developed in the present work. These boundary conditions take into account a recently developed theory: influence of temperature space inhomogeneity due to the interfaces and quasi-particles temperature-mismatch on thermo-generation and recombination. They also take into account surface resistance, surface recombination rates and possible temperature discontinuities at the interface due to finite surface thermo-conductivity. The temperature difference between current-carriers and phonon subsystems is also included in this approach.

This work addresses the problem of defining appropriate boundary conditions (BCs) in a closed electronic circuit of any type in the simplest general case. It is a continuation and generalization of previous studies published elsewhere.g1 ; g2 ; g3 ; g4 ; g5  They have been undertaken to clarify different processes arising in connection with real working conditions of electronic devices in general,Singh ; Sze  and in the search for applications to design new optoelectronic and thermoelectric devices.Aguado ; DiSalvo ; Fujiwara

In spite of great advances in the field of contemporary electronics, the effects of any boundary between two arbitrary materials have been studied, in practice, only for open circuits. The presence of both charge carriers (like electrons and holes in semiconducting materials; Cooper pairs in superconductors; etc…) in circuits containing heterojunctions, makes a big difference when the circuit is closed. Additionally, an intrinsic nonlinear effect due to the temperature difference between electrons and holes, from one side, and phonons, from the other (hot electrons), is usually forgotten (see, for example Ref. Tautc, ; Anselm, ; Ioffe, ). This effect has been studied in Ref. g5, .

Actually, one can mention very few works handling closed circuits and non-equilibrium carriers (see Ref. Balmush,  and reference therein). In more recent worksg1 ; g2 ; g4  it has been shown that it is compulsory to take into account non-equilibrium carriers in bipolar and $`p`$-type semiconductors even in a linear approximation by electric field; furthermore, recombination starts to play an important role even in the linear regime and cannot be neglected.g3 ; g4 ; g5

In practice, people use zero-current BCs to define the functionality of solid state electronic devices in general.Singh ; Sze  This is an ideal assumption that works well in many cases, particularly when there is no different kind of charge and heat carriers and interfaces involved. For a proper description of real performance conditions in heterojunctions, it should not be ignored that all devices are working in a mode such that electric current flows through the ends of a semiconducting structure and in a whole closed circuit, completed with another semiconductor, a normal metal or a superconducting material. Typically the semiconductors have both charge carriers, electrons and holes, and important recombination processes take place at the ends. All these facts demand a more stringent choice of BCs, and in particular of non-vanishing current BCs at the borders between different elements of a circuit to setup the correct working conditions of solid state electronic devices in general. These BCs become more important nowdays when accurate measurements of current to the level of holes and electrons are feasable.Fujiwara

In what follows, general BCs at the contact between two arbitrary conducting media will be obtained. Next we exemplify the usual case of a metal-semiconductor boundary. Finally we comment our results.

BCs (see, for example, Ref. Landau, ) can be obtained from the transport equations for non-equilibrium carriers.Bonch  In the static case they look like:

$$\mathrm{div}j_n=eR_n,\mathrm{div}j_p=eR_p.$$
(1)
Here $`R_{n,p}`$ are electron and hole recombination rates;

$$j_n=\sigma _n\left[\frac{d\stackrel{~}{\phi }_n}{dx}\alpha _n\frac{dT_n}{dx}\right],j_p=\sigma _p\left[\frac{d\stackrel{~}{\phi }_p}{dx}\alpha _p\frac{dT_p}{dx}\right],$$
(2)
are electron and hole electric currents; $`\sigma _{n,p}`$—electron and hole conductivity;

$$\stackrel{~}{\phi }_{n,p}=\phi \frac{\mu _{n,p}}{e}$$
(3)
are the electrochemical potentials (for exhaustive discussion see Ref. gm, ), where $`\phi `$ is the electric potential and $`\mu _{n,p}`$ are electron and hole chemical potentials; the latters are related by $`\mu _p=ϵ_g\mu _n`$ ($`ϵ_g`$ is the semiconductor band gap) only under thermodynamic equilibrium conditions (for exhaustive discussion see Ref. g3, ); ($`e`$) is the electron charge ($`e>0`$); $`\alpha _n`$ and $`\alpha _p`$, electron and hole power coefficients; $`T_n`$ and $`T_p`$—electron and hole temperatures (for exhaustive discussion see Ref. gm, ). The recombination rates are usually written out in the following form: Neamen

$$R_n=\frac{\delta n}{\tau _n},R_p=\frac{\delta p}{\tau _p}.$$
(4)
Here $`\tau _{n,p}`$ are electron and hole lifetimes; $`\delta n=nn_0`$ and $`\delta p=pp_0`$ are electron and hole non-equilibrium concentrations; $`n`$ and $`p`$ are the full concentrations; $`n_0`$ and $`p_0`$ are the corresponding equilibrium values; $`\delta n`$ and $`\delta p`$ are related to fluctuations in the chemical potentials $`\delta \mu _{n,p}`$.g4  Notice that these and the Fermi quasilevels $`\delta \stackrel{~}{\phi }_{n,p}(x)=\delta \phi (x)(1/e)\delta \mu _{n,p}(x)`$ are inhomogeneous across a boundary.

It follows from Maxwell’s equations that $`\mathrm{div}(j_n+j_p)=0`$ in the steady state, so that taking into account Eq. (4) one obtains from Eq. (1) the unphysical condition: $`(\delta n/\tau _n)=(\delta p/\tau _p)`$. Gurevich et al. (in Ref. g4, ) have shown that Eq. (4) is thermodynamically incorrect; moreover, Eqs. (4) should always be replaced by the following relationship (according to Ref. g3, ):

$$R_n=R_p=R,R=\frac{\delta n}{\tau _n}+\frac{\delta p}{\tau _p}.$$
(5)
The system of Eqs. (1) is incomplete as we have three unknown functions $`\delta n(x)`$, $`\delta p(x)`$, and $`\phi (x)`$. Poisson equation might be used to complete it. For simplicity, we will assume that all characteristic lengths are much bigger than the Debye’s radius,Silin  so that:

$$\delta n(x)=\delta p(x)$$
(6)
and Poisson equation becomes unnecessary.

Let us assume that the boundary between media 1 and 2 lies at $`x=0`$. Integrating Eqs. (1) with $`x`$ in a short range from $`x=\delta `$ to $`x=+\delta `$ and taking the limit $`\delta 0`$ one obtains:

$$j_n(+0)j_n(0)=+eR_s,j_p(+0)j_p(0)=eR_s.$$
(7)
Here

$$R_s=\underset{\delta 0}{lim}_\delta ^{+\delta }R(x)𝑑x$$
(8)
is the surface recombination rate. Making use of Eq. (5) it can be rewritten as:g3

$$R_s=S_n\delta n+S_p\delta p,$$
(9)
or, with the help of Eq. (6),

$$R_s=S\delta n,S=S_n+S_p.$$
(10)
Notice that the recombination rate $`R`$ in Eq. (8) is $`x`$ dependent. Each one of the contacting media changes its properties in a distance close to the Debye’s radius and becomes inhomogeneous in that region. This explains the $`x`$ dependence of $`R`$ in Eq. (8).

Let us now integrate Eq. (1) with $`dx`$ from $`\epsilon `$ to $`\delta `$, and with $`d\epsilon `$ from $`\delta `$ to $`+\delta `$. One obtains:

$`j_n(+0)=`$ $`\sigma _n^s\left\{\stackrel{~}{\phi }_n(0)\stackrel{~}{\phi }_n(+0)\alpha _n^s\left[T_n(+0)T_n(0)\right]\right\}+eR_s^{n+}`$ (11)
$`j_p(+0)=`$ $`\sigma _p^s\left\{\stackrel{~}{\phi }_p(0)\stackrel{~}{\phi }_p(+0)\alpha _p^s\left[T_p(+0)T_p(0)\right]\right\}eR_s^{p+}.`$
Here $`\sigma _n^s`$ and $`\sigma _p^s`$ are electron and hole surface conductivity:

$$(\sigma _{n,p}^s)^1=\underset{\delta 0}{lim}_\delta ^{+\delta }(\sigma _{n,p})^1(\xi )𝑑\xi ,$$
(12)
$`\alpha _{n,p}^s`$ are thermopower surface coefficients:

$$\alpha _{n,p}^s=\frac{1}{T(\delta )T(\delta )}\underset{\delta 0}{lim}_{T(\delta )}^{T(\delta )}\alpha _{n,p}(T)𝑑T,$$
(13)
$`R_s^{n+,p+}`$ are surface recombination rates for electrons and holes to the right of the boundary $`x=0`$:

$$R_s^{n+}=\sigma _n^s\underset{\delta 0}{lim}_\delta ^{+\delta }𝑑\epsilon \sigma _{n}^{}{}_{}{}^{1}(\epsilon )_\epsilon ^\delta 𝑑xR(x)R_s^{p+}=\sigma _p^s\underset{\delta 0}{lim}_\delta ^{+\delta }𝑑\epsilon \sigma _{p}^{}{}_{}{}^{1}(\epsilon )_\epsilon ^\delta 𝑑xR(x).$$
(14)
Using the same procedure, the following BCs are obtained for the left:

$`j_n(0)=`$ $`\sigma _n^s\left\{\stackrel{~}{\phi }_n(0)\stackrel{~}{\phi }_n(+0)\alpha _n^s\left[T_n(+0)T_n(0)\right]\right\}eR_s^n`$ (15)
$`j_p(0)=`$ $`\sigma _p^s\left\{\stackrel{~}{\phi }_p(0)\stackrel{~}{\phi }_p(+0)\alpha _p^s\left[T_p(+0)T_p(0)\right]\right\}+eR_s^p.`$
$`R_s^{n,p}`$ are surface recombination rates for electrons and holes to the left of the boundary $`x=0`$:

$$R_s^n=\sigma _n^s\underset{\delta 0}{lim}_\delta ^{+\delta }𝑑\epsilon \sigma _{n}^{}{}_{}{}^{1}(\epsilon )_\delta ^\epsilon R(x)𝑑xR_s^p=\sigma _p^s\underset{\delta 0}{lim}_\delta ^{+\delta }𝑑\epsilon \sigma _{p}^{}{}_{}{}^{1}(\epsilon )_\delta ^\epsilon R(x)𝑑x.$$
(16)
From Eqs. (8), (14), and (16) it becomes evident that

$$R_s^{n+}+R_s^n=R_s^{p+}+R_s^p=R_s$$
(17)
therefore BCs (7), (11), and (15) are not independent and only two of them should be used.

One should base particular decision of which BCs to use on experimental setup or physical sense of the problem to solve. Besides, this decision depends a lot on particular properties of the contact (Schottky barrier, $`p`$-$`n`$-junction, $`n^+`$-$`n`$-contact, etc…). We will demonstrate below how to choose the correct BCs for a metal-semiconductor junction.

We would like to emphasize two important facts. Firstly, it is not enough to define general surface recombination rates for the correct definition of transport effects on contact (this follows from Eqs. (11) and (15)). One has to use particular surface recombination rates for every type of current carriers (electrons, holes, etc.) defined on both sides of the contact. As far as we are aware, this has not been taken into account previously by anybody.

Secondly, we have used electron’s and hole’s Fermi quasilevels in Eqs. (11) and (15). One assumes that the electron’s Fermi quasilevel is measured from the bottom of the conduction band, and that the hole’s quasilevel is measured from the top of the valence band. Usually the position of the top of the valence band and the bottom of the conduction band are different in heterocontacts. That is why the reference (zero) points for the Fermi quasilevels are different at $`x=+0`$ and $`x=0`$.

Next we will describe how to take into account the above mentioned differences in a metal–semiconductor boundary.

It is very common to inject electric current into semiconductor samples through metal contacts. These contacts have an interesting peculiarity: semiconductors can be characterized by the presence of both electrons and holes, while metals have only electrons as current carriers.

Let us take into account this fact in Eqs. (7), (11), and (15). We assume that we have a metal to the left of $`x=0`$ and a semiconductor to the right. Therefore we have that $`j_p^m(0)=0`$, so that Eq. (7) leads to $`j_p^s(+0)=eR_s`$ (the subscript indicates the corresponding media; here “s” stands for semiconductor and “m” for metal).

We can simplify Eq. (15) due to the absence of surface recombination in a metal (it has only electrons as charge carriers): $`R_s^n=R_s^p=0`$. Therefore, from Eq. (17) it follows that:

$$R_s=R_s^{n+}=R_s^{p+}.$$
(18)
There are no positive charge carriers in metals, which means that electron current in the metal, $`j_n^m(0)`$, should be equal to the whole current $`j_0`$. In this case, the BCs (15) reduce to:

$$j_0=\sigma _n^s\left\{\phi _m\phi _s(0)\frac{\mu _m}{e}+\frac{\mu _s(0)}{e}+\frac{\mathrm{\Delta }\epsilon _c}{e}\alpha _n^s\left[T_n^s(0)T_n^m(0)\right]\right\}.$$
(19)
Notice that the electrical and the chemical potential of the metal have not changed. Let us remind that $`\mu _m`$ and $`\mu _s`$ are calculated from the bottom of the conduction band of each material. This leads to the additional term in Eq. (19) $`(\sigma _n^s/e)\mathrm{\Delta }\epsilon _c`$, where $`\mathrm{\Delta }\epsilon _c`$ is the distance between the metal and the semiconductor conduction bands.

It follows from the second equation of system (15) that $`\sigma _p^s=0`$. This is to be expected since the absence of holes in a metal means the absence of holes surface conductivity.

Summarizing all the above mentioned facts, we can write Eq. (11) as

$$j_n^s(0)=\sigma _n^s\left\{\phi _m\phi _s(0)\frac{\mu _m}{e}+\frac{\mu _s(0)}{e}+\frac{\mathrm{\Delta }\epsilon _c}{e}\alpha _n^s\left[T_n^s(0)T_n^m(0)\right]\right\}+eR_s,j_p^s(0)=eR_s.$$
(20)
The last condition in (20) is equal to the second condition in (7). The first equation in (7) can be rewritten as

$$j_n^s(0)j_0=eR_s.$$
(21)
Notice that the first equation (20) follows from Eqs.(19) and (21). All this yields the new BCs for a metal-semiconductor junction which are now formed by Eq. (19), Eq. (21), and the second of equations (20).

In conclusion, we have formulated general BCs, corresponding to current flow through the boundary between two conducting media. These conditions consider possible jumps of the electron’s and hole’s Fermi quasilevels and of the electric potentials at the boundary; they also take into account the surface recombination rates. The general procedure has been applied to a usual contact, namely the metal-semiconductor boundary, but the method is valid for all types of contacts ($`n^+`$-$`n`$, $`p^+`$-$`p`$, $`p`$-$`n`$, Schottky barrier, etc…) between different conducting materials.

This work was partially supported by CONACyT, Mexico. JG has received support from DINAIN (Universidad Nacional de Colombia), Colombia.


text:
# One-body Properties of Nuclear Matter with Off-shell Propagation

## I Introduction

Single-particle properties of nucleons are modified inside the strongly interacting nuclear matterMahaux and Sartor (1991). The optical potential describes the average interaction of a nucleon with the medium. At negative energies, the single-particle potential reflects the modifications of the quasi-particles in the nuclear matter. The spectral function of a nucleon is a measure of the energy distribution of a plane wave state in the system, and can be observed in electron scattering experiments.

Single-particle properties can be calculated from the Brueckner-Hartree-Fock (BHF) theory Köhler (1992); Baldo et al. (1992). Hole-hole diagrams must be added to the lowest order results to obtain a meaningful imaginary part of the self-energy. Because of the violation of the Hugenholz-Van Hove property by the BHF approach, the overall normalization of the single-particle energies is sometimes a problem. BHF calculations use a quasi-particle approximation for the propagators; therefore, a finite order diagram has kinematical limits on the total energy in the scattering. A finite total energy in the interaction could introduce distortions in the far energy tails of the calculated spectral functions.

Recently the in-medium $`T`$-matrix approximation has been applied to the nuclear matter Vonderfecht et al. (1993); Alm et al. (1996); Dickhoff (1998); Bożek (1999a); Dickhoff et al. (1999); Dewulf et al. (2001); Bozek and Czerski (2001). One obtains nontrivial self-energies by summing a ladder of hole-hole and particle-particle propagators. The self-energies obtained are $`\mathrm{\Phi }`$-derivable Baym (1962); Bozek and Czerski (2001); within this approximation, the single-particle properties are consistent with global thermodynamical properties of the system. A necessary requirement for the in-medium $`T`$-matrix approach is to use fully self-consistent self-energies and propagators. The propagators are dressed with the imaginary and real parts of the self-energy, and the full spectral function for nucleons must be taken. It requires a serious numerical effort in the calculations, but using a quasi-particle approximation in the $`T`$-matrix ladder gives too strong pairing Bożek (1999a); Bozek (2000a); Dickhoff et al. (1999), too much scattering Dickhoff et al. (1999); Bożek (1999a), and incorrect single-particle energies Bozek and Czerski (2001).

In the present work we calculate single-particle self-energies and spectral functions using fully self-consistent propagators. We take the complete energy dependence of the self-energies and spectral functions into account in order access reliably the high energy tails of the spectral functions. We present results on the single-particle potential, scattering width and spectral function, obtained with a separable Paris interaction with several partial waves (App. B). In Appendix A we describe in details the numerical methods that allow us to tackle the problem.

We find that normal nuclear density is at the limit of the region of the $`{}_{}{}^{3}S_{1}^{}^3D_1`$ pairing instability. Owing to the dressing of propagators in the gap equation the superfluid gap is very small, if not vanishing (depending on the details of the interaction in the deuteron channel). We leave the detailed analysis of the pairing in symmetric nuclear matter with realistic interactions to a different work, and study in the following the normal state of nuclear matter.

## II In-medium $`T`$-matrix

In this Section we recall the basic equations of the approach used, more details and discussion can be found in Kadanoff and Baym (1962); Danielewicz (1984); Kraeft et al. (1986); Alm et al. (1996); Vonderfecht et al. (1993); Dickhoff (1998); Bożek (1999a). The $`T`$-matrix approximation sums ladder diagrams with dressed particle-particle and hole-hole propagators for the in-medium two-particle propagator. The retarded $`T`$-matrix is

$`𝐩|T(𝐏,\mathrm{\Omega })|𝐩^{^{}}=V(𝐩,𝐩^{^{}})`$
$`+{\displaystyle \frac{d\omega _1d\omega _2d^3q}{32\pi ^5}V(𝐩,𝐪)A(p_1,\omega _1)A(p_2,\omega _2)}`$
$`{\displaystyle \frac{\left(1f(\omega _1)f(\omega _2)\right)}{\mathrm{\Omega }\omega _1\omega _2+iϵ}}𝐪|T(𝐏,\mathrm{\Omega })|𝐩^{^{}}`$ (1)
where $`𝐩_{\mathrm{𝟏},\mathrm{𝟐}}=𝐏/2\pm 𝐪`$ and $`f(\omega )`$ is the Fermi distribution. A partial wave expansion of the in-medium $`T`$-matrix is performed

$`p|T_{ll^{^{}}}^{(JST)}(P,\mathrm{\Omega })|p^{^{}}=V_{ll^{^{}}}^{(JST)}(p,p^{^{}})`$
$`+{\displaystyle \underset{l^{^{\prime \prime }}}{}}{\displaystyle \frac{d\omega _1d\omega _2q^2dq}{8\pi ^4}V_{ll^{^{\prime \prime }}}^{(JST)}(p,q)}`$
$`A(p_1,\omega _1)A(p_2,\omega _2)_{P,q}`$
$`{\displaystyle \frac{\left(1f(\omega _1)f(\omega _2)\right)}{\mathrm{\Omega }\omega _1\omega _2+iϵ}}q|T_{l^{^{\prime \prime }}l^{^{}}}^{(JST)}(P,\mathrm{\Omega })|p^{^{}}`$ (2)
after angle averaging the intermediate two-particle propagator ($`\mathrm{}=\frac{d\mathrm{\Omega }}{4\pi }\mathrm{}`$)

The imaginary part of the retarded self-energy is obtained by closing a pair of external vertices in the $`T`$-matrix with a fermion propagator

$`\mathrm{Im}\mathrm{\Sigma }(p,\omega )={\displaystyle \frac{d\omega _1d^3k}{16\pi ^4}A(k,\omega _1)}`$
$`(𝐩𝐤)/2|\mathrm{Im}T(𝐩+𝐤,\omega +\omega _1)|(𝐩𝐤)/2_A`$
$`\left(f(\omega _1)+b(\omega +\omega _1)\right)`$ (3)
where

$`A(p,\omega )=`$
$`{\displaystyle \frac{2\mathrm{I}\mathrm{m}\mathrm{\Sigma }(p,\omega )}{\left(\omega p^2/2m\mathrm{Re}\mathrm{\Sigma }(p,\omega )\right)^2+\mathrm{Im}\mathrm{\Sigma }(p,\omega )^2}}`$ (4)
is the self-consistent spectral function of the nucleon and $`b(\omega )`$ is the Bose distribution. The real part of the self-energy is related to $`\mathrm{Im}\mathrm{\Sigma }`$ by a dispersion relation

$$\mathrm{Re}\mathrm{\Sigma }(p,\omega )=\mathrm{\Sigma }_{HF}(p)+𝒫\frac{d\omega ^{^{}}}{\pi }\frac{\mathrm{Im}\mathrm{\Sigma }(p,\omega ^{^{}})}{\omega ^{^{}}\omega }$$
(5)
with $`\mathrm{\Sigma }_{HF}(p)`$ being the Hartree-Fock self-energy. Eqs. (II), (II), (5) and (II) are to be solved iteratively and at each iteration the chemical potential $`\mu `$ is adjusted to fulfill the condition on the density

$$\rho =\frac{d\omega d^3p}{16\pi ^4}A(p,\omega )f(\omega )=.16\text{fm}^3.$$
(6)
The $`T`$-matrix approximation sums ladder diagrams contributing to the ground state energy. In this way it regularizes the short range core in the nucleon interaction similarly as the BHF approach. At the same time, the self-consistent $`T`$-matrix approximation is $`\mathrm{\Phi }`$-derivable Baym (1962). It was shown for a model interaction that if a fully self-consistent $`T`$-matrix calculation is performed the thermodynamic consistency relations are fulfilled Bozek and Czerski (2001). In the $`T`$-matrix approximation a nonzero imaginary self-energy appears and leads to a nontrivial spectral function. It is crucial to keep the full off-shell dependence of the propagators in the calculations. It is possible to treat the resulting energy integrals using numerical methods described in the Appendix A. All the calculations are done using the Fermi energy as the origin of the energy scale. The figures are also plotted using this convention. The calculations are done using a separable parameterization of the Paris potential Haidenbauer and Plessas (1984, 1985). The Fermi energy obtained is $`21.4\text{MeV}`$ and the binding energy $`15.1\text{MeV}`$ at normal density. This is not the saturation point for this interaction, but the values quoted above are quite reasonable and give some confidence in the single-particle properties we want to study.

## III Optical potential

The real-part of the self-energy defines the single-particle pole in the propagator

$$\omega _p=p^2/2m+\text{Re}\mathrm{\Sigma }(p,\omega _p).$$

The free dispersion relation is modified due to interactions with the medium. The resulting effective potential $`\text{Re}\mathrm{\Sigma }(p,\omega _p)`$ is attractive. It leads in particular to a reduction of the effective mass $`m^{}=\frac{pdp}{d\omega _p}`$. We find $`m^{}.85m`$ at the Fermi momentum, this value depends somewhat on the chosen interaction. The effective mass can be written Mahaux and Sartor (1991) as the product of the $`k`$-mass $`m_k/m=\left(1+\frac{m\text{Im}\mathrm{\Sigma }(k,\omega )}{kk}|_{\omega =\omega _k}\right)^1`$ and the $`\omega `$-mass $`m_\omega /m=\left(1\frac{\text{Im}\mathrm{\Sigma }(k,\omega )}{\omega }|_{\omega =\omega _k}\right)`$. We find $`m_k/m.62`$ and $`m_\omega /m1.37`$ at the Fermi momentum.

The real part of the self-energy is the sum of the Hartree-Fock and the dispersive contribution (5). The real part of the dispersive self-energy is negative around the Fermi energy and around the quasi-particle pole (Fig. 1). The value of the real part of the self-energy on-shell is the effective potential felt by the quasi-particle in the medium. For positive energies of the particle it corresponds to the optical potential.

The depth of the single particle potential is around $`70\text{MeV}`$ and it is decreasing with momentum. In the range of positive single particle energies ($`280\text{MeV}<p<500\text{MeV}`$), the single-particle potential can be fitted with the form Perey and Buck (1962)

$$\text{Re}\mathrm{\Sigma }(p,\omega _p)H_0\mathrm{exp}\left(\beta ^2k^2/4\right)$$

with $`H_0=70\text{MeV}`$ and $`\beta =.58\text{fm}`$ (dashed-dotted line in Fig. 1). The range of nonlocality $`\beta `$ is small reflecting the fact that the effective mass $`m^{}`$ is not very different from the free one. The depth of the potential for positive energies is consistent with values obtained from phenomenological analyzes of the optical potential Varner et al. (1991). The single particle potential is weakly dependent on the temperature; at $`T=4\text{MeV}`$ $`\text{Re}\mathrm{\Sigma }(p,\omega _p)`$ is shifted up by less than $`1\text{MeV}`$ from its value at zero temperature.

In Fig. 2 is presented the energy dependence of the real part of the self-energy. The energy dependence of the single particle potential is responsible for the time nonlocality of the optical potential. In the whole range of frequencies in Fig. 2 some energy dependence of the potential can be seen. This energy dependence is decreasing with momentum. For energies close to the Fermi energy the energy dependence of the single-particle potential is similar to the one presented in Ref. Baldo et al. (1992). The energy dependence determines the quasi-particle strength

$$Z_p=\left(1\frac{\text{Re}\mathrm{\Sigma }(p,\omega )}{\omega }|_{\omega =\omega _p}\right).$$
(7)
At the Fermi momentum we find $`Z_{p_F}.73`$. On the other hand, in the range $`50\text{MeV}<\omega <200\text{MeV}`$ $`\text{Re}\mathrm{\Sigma }(p,\omega )`$ is relatively flat as function of energy.

## IV Scattering width

The off-shell scattering width is nonzero on both sides of the Fermi energy but larger at positive energies (Fig. 3). At negative energies the tail of the scattering width extends very far. It is a general feature of self-consistent calculations. Dressed propagators in the $`T`$-matrix ladder make the imaginary part of the $`T`$-matrix nonzero even at very negative energies. On the other hand, quasi-particle approximations have a kinematical limit on the lowest energy in the scattering. The scattering width decreases with momentum. It can be understood by the fact that at low momenta the interaction is dominated by the strongest $`S`$ waves. Other calculations Benhar et al. (1992); Köhler (1992); Baldo et al. (1992); Bozek (2000b); Lehr et al. (2000); Bożek (1999a); Dickhoff et al. (1999); Dewulf et al. (2001) of the variational, BHF, Born, or $`T`$-matrix type give qualitatively but not quantitatively similar results. The differences are partly due to different short range properties of the interactions used.

Consistently with general features of Fermi liquids at zero temperature, we find that $`\text{Im}\mathrm{\Sigma }(p,\omega =0)=0`$ (Fig. 3). The same relation, coming form the restricted phase-space, is fulfilled by other approximations Benhar et al. (1992); Baldo et al. (1992); Lehr et al. (2000). At higher temperatures a nonzero scattering width at the Fermi energy appears.

As can be seen from Fig. 4 the scattering width at the Fermi energy is increasing as the square of the temperature, as expected from phase space arguments Abrikosov et al. (1963).

The large value of the scattering width above the Fermi energy is due to the short range part of the interaction potential. Splitting the total width into contributions from different partial waves shows that the $`{}_{}{}^{3}S_{1}^{}^3D_1`$ partial wave is dominant (Fig. 5). For momenta up to twice the Fermi momentum and energies around the Fermi energy, the deuteron partial wave gives by far the most important contribution to the scattering width. It is also in this kinematical region that the overall scattering width is the largest. We have observed a similar large contribution due to this partial wave for the Monagan Mongan (1969) and Yamaguchi Yamaguchi (1954) separable potentials. Consistently, for all these interactions we obtain similar values for $`Z_{p_F}.7`$. Large values of the scattering width at large energies lead to long tails in the spectral functions.

The imaginary part of the self-energy at the quasi-particle pole is the scattering width of the quasi-particle.

From Fig. 6 we see that around the Fermi momentum, narrow quasi-particle excitation exist. The scattering width behaves as $`\text{Im}\mathrm{\Sigma }(p,\omega _p)(pp_F)^2`$, as expected from restricted phase space for scattering Abrikosov et al. (1963). The on-shell scattering width does not reflect the very large values of the off-shell scattering width (Fig. 3). For momenta close to the Fermi surface it is small; and particles with large energies have a scattering probability proportional to the total density and cross section. At finite temperature quasi-particles at the Fermi surface aquire a finite life-time (Fig. 6).

## V Spectral function

The role of correlation induced by the medium can be judged by the modifications of the spectral function. Nontrivial, dispersive self-energies lead to broad spectral function with non-Lorentzian shapes and long tails. High energy parts of the spectral functions could be revealed by electron scattering experiments. It is important to calculate the form of the spectral functions including contributions form short range correlations. In Fig. 7 are presented the spectral functions for three representative momenta. For zero momentum the spectral function has a broad peak below and a background part above the Fermi energy. It is the reverse

for $`p=355\text{MeV}`$. Following the behavior of the scattering width, the spectral function for any momentum goes to zero at the Fermi energy. For momenta close to the Fermi energy the spectral function has a very sharp peak at the single-particle energy. Its has also a significant background part which cannot be ignored in the sum rules or in the calculation of effective interactions between quasi-particles. Similarly as for the scattering width the self-consistent calculation gives rise to a long tail in the spectral function at negative energies.

The spectral function can be used to obtain the momentum distribution

$$n(p)=\frac{d\omega }{2\pi }A(p,\omega )f(\omega ).$$
(8)
Short range correlations modifying the spectral function are reflected in the nucleon momentum distribution. The free Fermi distribution is depleted below $`p_F`$ and a high momentum tail in $`n(p)`$ is formed.

At the Fermi momentum the discontinuity in the occupancy is reduced from $`1`$ to $`Z_{p_F}`$ (7). The Fermi momentum in the interacting system should be the same as in the free one Luttinger (1960); Baym (1962). This relation is approximately fulfilled in the present calculation, but not as well as in our previous work using only $`S`$ wave interactions Bozek and Czerski (2001). A possible source of the discrepancy could be the use of the partial-wave expansion which spoils the $`\mathrm{\Phi }`$-derivability of the self-energies. At finite temperature the Fermi surface is washed out, but the depletion of the Fermi sea and the high momentum tail are the same. It means that the short range correlations stay essentially the same at $`T=4\text{MeV}`$. The same can also be seen by inspecting the self-energies. At finite temperature the scattering width is increased only in the vicinity of the Fermi energy.

Long tails in the spectral function have implications for saturation properties of nuclear matter. Increased kinetic energy due to the high momentum tail in $`n(p)`$ is compensated by increased removal energies due to the negative energy tail of the spectral function. A detailed discussion at different densities will be presented elsewhere.

## VI Summary

We study the single-particle properties of nucleons in nuclear matter using a conserving, in-medium $`T`$-matrix approximation. The calculations are done in a self-consistent way with dressed propagators in the ladder diagrams. To our knowledge, it is the first such calculation in the literature using realistic interactions and several partial waves. The full discretization of the spectral function and self-energies allows to discuss the details of their energy dependence. We find that the basic features of a consistent approximation to a Fermi liquid system are fulfilled. The scattering width is zero at the Fermi energy, and the quasi-particles at the Fermi surface have infinite lifetime. The momentum occupancy has a discontinuity of $`Z_{p_F}`$ at the Fermi momentum. The off-shell scattering width is very large at energies $`300400\text{MeV}`$ and momenta below $`400\text{MeV}`$. In this region the main contribution comes from the $`{}_{}{}^{3}S_{1}^{}^3D_1`$ partial wave. The $`{}_{}{}^{3}P_{1}^{}`$ and $`{}_{}{}^{1}S_{0}^{}`$ partial waves give also important contributions to the on-shell self-energies. The scattering width on-shell is not very large, maximally $`26\text{MeV}`$. Finite temperature effects are concentrated around the Fermi surface. At finite temperature, the scattering width gets finite around the Fermi energy without modifying the real single-particle potential and the short range correlations. The quasi-particle renormalization factor is $`Z_{p_F}.73`$ and is largely independent on the interaction used. The effective mass is $`m^{}.85m`$.

###### Acknowledgements.

This work was partly supported by the KBN under Grant 2P03B02019.

## Appendix A  Numerical methods

Calculations using self-consistent spectral functions require the evaluation of energy integrals in Eqs. (II) and (II). This is the main numerical difficulty for any self-consistent approximation; only in the last years self-consistent approaches have been applied to the study of high $`T_c`$ superconductors Haussmann (1994); Pedersen et al. (1997); Kyung et al. (1998); Letz and Marsiglio (1999) and nuclear matter Dickhoff (1998); Dickhoff et al. (1999); Bożek (1999a); Bozek and Czerski (2001); Dewulf et al. (2001). Some calculations are performed in the imaginary time formalism Haussmann (1994); Pedersen et al. (1997) which requires a numerical procedure for the analytical continuation to calculate the spectral function. A simpler way is to use the real-time formalism operating with the retarded $`T`$-matrix, the retarded self-energy, and the spectral function Kyung et al. (1998); Bożek (1999a). The real-time formalism was also used for other calculations of the nuclear matter Dickhoff et al. (1999); Dewulf et al. (2001) performed at zero temperatures. To deal with the off-shell propagation a numerical parameterization of the energy dependence of the spectral function $`A(p,\omega )`$ by a set of Gaussians has been used Dickhoff (1998); Dickhoff et al. (1999). Alternatively, the spectral function can be represented as a sum of $`\delta `$ functions Letz and Marsiglio (1999); Dewulf et al. (2001). The above two methods can be easily applied at zero temperature, where a narrow Gausian or one of the $`\delta `$ functions describes the quasi-particle peak for momenta close to the Fermi surface.

The alternative approach, here employed, uses a direct discretization of the spectral function and self-energies as functions of the energy. If the discretization in $`\omega `$ is uniform the energy integrals in Eqs. (II) and (II) can be performed efficiently by Fourier transforms. This method can be used directly at finite temperature Bożek (1999a); Kyung et al. (1998), but not for low temperatures where the quasi-particles are very narrow around the Fermi surface (Fig. 6). At zero temperature we have

$$A(p_F,\omega )=2\pi Z_{p_F}\delta (\omega E_F)+B(p,\omega )$$
(9)
where the background part $`B`$ is smooth and can be discretized. Generally, for any momentum close to the Fermi momentum the spectral function is rapidly changing in the vicinity of the quasi-particle peak; therefore, it is has to be split into a quasi-particle peak and a smooth background. The sharp peak is approximated by a $`\delta `$ function and the smooth part is discretized. We have found numerically that this separation in the spectral function is needed if the width of the quasi-particle peak $`2\text{Im}\mathrm{\Sigma }(p,\omega _p)`$ is smaller than $`3`$ times the spacing in the energy discretization $`\mathrm{\Delta }\omega `$. We use the following representation of the spectral function

$$A(p,\omega )=2\pi 𝒵_p\delta (\omega \omega _p)+B(p,\omega ).$$
(10)
The background part is defined as

$$B(p,\omega )=\{\begin{array}{c}\frac{2\mathrm{I}\mathrm{m}\mathrm{\Sigma }(p,\omega )}{(\omega \omega _p)^2+\mathrm{Im}\mathrm{\Sigma }(p,\omega )^2};\\ \text{for}(\omega \omega _p)^2+\mathrm{Im}\mathrm{\Sigma }(p,\omega )^2>\eta \\ \\ 2\mathrm{I}\mathrm{m}\mathrm{\Sigma }(p,\omega )/\eta ;\\ \text{for}(\omega \omega _p)^2+\mathrm{Im}\mathrm{\Sigma }(p,\omega )^2<\eta \end{array}.$$
(11)
The parameter $`\eta `$ is set to cut off the rapidly changing part of the spectral function. We use $`\eta 4\mathrm{\Delta }\omega ^2`$. The strength of the quasi-particle component $`𝒵_p`$ in Eq. (10) is adjusted to conserve the sum rule $`\frac{d\omega }{2\pi }A(p,\omega )=1`$. The weight of the singular part $`𝒵_p`$ is not the same as the renormalization factor $`Z_p`$ (Eq. 7). In Fig. 9 is shown an example of the separation of the background and singular parts of the spectral function $`A(p=255\text{MeV},\omega )`$.

For a separable interaction $`V^\alpha (p,p^{^{}})=_{i,j}\lambda _{i,j}^\alpha g_i^\alpha (p)g_j^\alpha (p^{^{}})`$ the $`T^\alpha `$-matrix for the partial wave $`\alpha =(JST)`$ is $`p|T^\alpha (\omega ,P)|p^{^{}}=_{i,j}T_{ij}^\alpha (\omega ,P)g_i^\alpha (p)g_j^\alpha (p^{^{}})`$. Where the matrix $`T_{ij}^\alpha =(\lambda _{ij}^\alpha J_{ij}^\alpha )^1`$ with the coefficients $`J_{ij}^\alpha `$ given by the integrals

$`J_{ij}^\alpha (\mathrm{\Omega },P)={\displaystyle \frac{d\omega _1d\omega _2d^3k}{32\pi ^5}g_i^\alpha (k)g_j^\alpha (k)}`$
$`{\displaystyle \frac{A(p_1,\omega _1)A(p_2,\omega _2)}{\mathrm{\Omega }\omega _1\omega _2+iϵ}}\left(1f(\omega _1)f(\omega _2)\right).`$ (12)
Substituting (10) for the spectral functions we have

$`\mathrm{Im}J_{ij}^\alpha (\mathrm{\Omega },P)={\displaystyle \frac{d\omega d^3k}{32\pi ^4}g_i^\alpha (k)g_j^\alpha (k)B(p_1,\omega )\left(12f(\omega )\right)B(p_2,\mathrm{\Omega }\omega )}`$
$`{\displaystyle \frac{d^3k}{8\pi ^3}g_i^\alpha (k)g_j^\alpha (k)B(p_1,\mathrm{\Omega }\omega _{p_2})𝒵_{p_2}\left(1f(\mathrm{\Omega }\omega _{p_2})f(\omega _{p_2})\right)}`$
$`+\mathrm{Im}{\displaystyle \frac{d^3k}{8\pi ^3}g_i^\alpha (k)g_j^\alpha (k)\frac{𝒵_{p_1}𝒵_{p_2}\left(1f(\omega _{p_1})f(\omega _{p_2})\right)}{\mathrm{\Omega }\omega _{p_1}\omega _{p_2}+iϵ}}.`$ (13)
The energy integral in the first term of Eq. (A) is a convolution of the functions $`B(p_1,\omega )\left(12f(\omega )\right)`$ and $`B(p_2,\omega )`$. It can be evaluated by Fourier transformations using FFT transform algorithms Press et al. (1986). Moreover, the Fourier transform and its inverse can be performed outside of the momentum integral Bożek (1999b). The second term of (A) is a standard two-dimensional integral without singularities. The integral in the last term of (A) is of the type often occurring in quasi-particle calculations, such as the $`G`$-matrix or the quasi-particle $`T`$-matrix approximations. Angle averaging separately the numerator and the denominator under the integral we have

$`\mathrm{Im}{\displaystyle \frac{d^3k}{8\pi ^3}g_i^\alpha (k)g_j^\alpha (k)\frac{𝒵_{p_1}𝒵_{p_2}\left(1f(\omega _{p_1})f(\omega _{p_2})\right)}{\mathrm{\Omega }\omega _{p_1}\omega _{p_2}+iϵ}}`$ $``$
$`\mathrm{Im}{\displaystyle \frac{k^2dk}{2\pi ^2}g_i^\alpha (k)g_j^\alpha (k)\frac{<𝒵_{p_1}𝒵_{p_2}\left(1f(\omega _{p_1})f(\omega _{p_2})\right)>_{P,k}}{\mathrm{\Omega }<\omega _{p_1}+\omega _{p_2}>_{P,k}+iϵ}}`$ $`=`$
$`{\displaystyle \frac{k_0^2}{2\pi }}{\displaystyle \frac{<𝒵_{p_1}𝒵_{p_2}\left(1f(\omega _{p_1})f(\omega _{p_2})\right)>_{P,k_0}}{<\omega _{p_1}+\omega _{p_2}>_{P,k}/k|_{k=k_0}}}`$ (14)
where $`k_0`$ is the solution of

$$\mathrm{\Omega }=<\omega _{p_1}+\omega _{p_2}>_{P,k_0}.$$

The angle average in the denominator $`<\omega _{p_1}+\omega _{p_2}>_{P,k}`$ is a function of the total and relative momentum ($`P`$ and $`k`$). It can be represented using a one-dimensional function

$$<\omega _{p_1}+\omega _{p_2}>_{P,k}=\frac{2}{Pk}\left(F(P/2+k)F(P/2k)\right)$$
(15)
with

$$F(x)=_0^xp𝑑p\omega _p.$$

This effectively one-dimensional parameterization allows to perform the integral corresponding to the quasi-particle part of the spectral function very efficiently without using any parabolic approximation for $`\omega _p`$. The same method can also be used in $`G`$-matrix calculations. The real part of the integrals $`J_{ij}^\alpha `$ is obtained using the dispersion relation

$$\text{Re}J_{ij}^\alpha (\mathrm{\Omega },P)=\frac{\mathrm{\Omega }^{^{}}}{\pi }\frac{\text{Im}J_{ij}^\alpha (\mathrm{\Omega }^{^{}},P)}{\mathrm{\Omega }^{^{}}\mathrm{\Omega }}.$$

The calculation of the energy integral in the equation for the self-energy (II) proceeds very similarly. We have

$`\mathrm{Im}\mathrm{\Sigma }(p,\omega )={\displaystyle \frac{d\omega _1d^3k}{16\pi ^4}B(k,\omega _1)}`$
$`(𝐩𝐤)/2|\mathrm{Im}T(𝐩+𝐤,\omega +\omega _1)|(𝐩𝐤)/2_A`$
$`\left(f(\omega _1)+b(\omega +\omega _1)\right)`$
$`+{\displaystyle \frac{d^3k}{8\pi ^3}}`$
$`(𝐩𝐤)/2|\mathrm{Im}T(𝐩+𝐤,\omega +\omega _k)|(𝐩𝐤)/2_A`$
$`\left(f(\omega _k)+b(\omega +\omega _k)\right).`$ (16)
The energy integral in the first term is again of the form of a convolution and is calculated using Fourier transforms. The integral in the second term is a standard two-dimensional integral. Unlike in the $`G`$-matrix approximation, the calculation of the self-energy requires the knowledge of the full off-shell $`T`$-matrix. Restricting oneself to contributions from on-shell $`T`$-matrix leads to erroneous results Schnell (1996).

In the iteration of the self-consistent set of equations (II, II, II, 5 and 6) it is advantageous to parameterize the energy dependence of the off-shell quantities with respect to the Fermi energy. The iterations are much faster, since a modification of the chemical potential between iterations does not change the most important features of the $`\omega `$ dependence of $`\mathrm{\Sigma }(p,\omega )`$ and $`T(\mathrm{\Omega },P)`$, e.g. $`\text{Im}\mathrm{\Sigma }(p,\omega =0)=0`$ and the pairing singularity for $`T(\mathrm{\Omega }=0,P=0)`$. The absolute energy scale is recovered only in the calculation of the total density (6) and final observables. This way of proceeding with the iterations can make use of self-energies calculated for other densities or temperatures to start a new iteration.

## Appendix B N-N interaction and numerical parameters

We use a separable parameterization of the Paris potential Haidenbauer and Plessas (1984, 1985). It contains all the $`J=0`$, $`J=1`$, $`J=2`$, and the $`{}_{}{}^{3}D_{3}^{}^3G_3`$ partial waves. For the most important $`{}_{}{}^{1}S_{0}^{}`$ and $`{}_{}{}^{3}S_{1}^{}^3D_1`$ partial waves we choose the rank $`3`$ and $`4`$ respectively. We have observed that the use of the separable parameterization of the Paris potential for the $`{}_{}{}^{3}P_{1}^{}`$ partial wave leads to very high values of the effective mass $`m^{}.95m`$. Since the $`{}_{}{}^{3}P_{1}^{}`$ phase shifts are not well reproduced by the parameterization of Ref. Haidenbauer and Plessas (1984), we choose the Mongan $`I`$ parameterization Mongan (1969) for this partial wave. It is important to check that no unphysical bound states occur in the off-shell $`T`$-matrix for the given choice of $`NN`$ potentials. The contribution from such unphysical bound states, even far from the considered on-shell energies, would spoil the calculated self-energies.

The use of the Fourier transform algorithm for the calculation of energy integrals requires a fixed spacing in the energy grid. It means that we have to set finite ranges for the kinematical variables in the model. The single particle momenta are in the range $`p<1700\text{MeV}`$. The total momentum of a nucleon pair is limited by $`3400\text{MeV}`$. The energy dependence of the $`T`$-matrix $`T(\mathrm{\Omega },P)`$ is calculated for $`2500\text{MeV}+P^2/4m<\mathrm{\Omega }<3500\text{MeV}`$. The energy range for the single particle self-energy is taken from $`1700\text{MeV}`$ to $`3500\text{MeV}`$. The energy dependent functions are discretized with a grid spacing of $`1.28\text{MeV}`$.

To obtain results with stability between iterations better than $`1\%`$ we need typically 7 iterations at zero temperature. The convergence of the iterations is faster and more stable at finite temperature.