format: text

text:
# Non-emission line young stars of intermediate mass

## 1 Introduction

Pre-main sequence stars of intermediate mass ($`2M/M_{}8`$) which show emission lines in their spectra are called Herbig Ae/Be (HAEBE) stars. They were first discussed as a group by Herbig (1960). He identified 26 stars which are of spectral type A or B, located in obscured starforming regions and illuminating a bright nebulosity in its immediate vicinity. Additions to the original list were made by Finkenzeller & Mundt (1984) and Herbig & Bell (1988). Recently, a more extended catalogue of HAEBE stars and related objects was published by Thé et al. (1994). In this catalogue 287 HAEBE stars and related objects are listed in five tables, which include stars with later spectral type (G0 or earlier) and those found relatively isolated from star forming clouds. Out of this, only 109 stars which are listed in Table 1 of the catalogue, are recognized as either HAEBE stars or potential candidate members. Other stars listed in the catalogue are either of very uncertain or unknown spectral type or have not been identified to belong to any specific group.

The following set of properties are often taken as a working definition of HAEBE stars (Waters & Waelkens 1998). (a) Spectral type B or A with emission lines. (b) Infrared excess due to hot or cool circumstellar dust or both. (c) Luminosity class III to V. The emission lines are believed to be formed in a stellar wind originating from a hot and extended chromosphere around the star (Bouret & Catala 1998), and is closely connected with the accretion related disc activity. The near and far infrared excesses which characterize the spectral energy distribution (SED) of these stars is attributed to the presence of significant amounts of circumstellar dust around them with a wide range in temperature (e.g. Hillenbrand et al. 1992). Submillimeter and millimeter observations have clearly established the existence of dust. Dust masses estimated from these studies range from $`10^5`$ to $`10^1M_{}`$ (Waters & Waelkens 1998 and references therein). Scattering of central star light by the circumstellar dust manifests in large values of intrinsic polarization measured for HAEBE stars ( e.g. Grinin 1994). However, the geometry of the circumstellar environment of HAEBE stars is still a matter of debate. Evidence for the presence of discs as well as envelopes has been found. Recently Natta et al. (2001) have argued that irradiated discs with a puffed-up inner wall of optically thick dust provide a good fit to the observations over the entire range of wavelengths.

Though the pre-main sequence (PMS) nature of the HAEBE stars is now well established, several questions concerning their PMS evolution remain to be answered. Do all PMS stars of intermediate mass go through the Herbig Ae/Be phase? What happens to the attendant circumstellar material around a HAEBE star by the time it evolves into a main sequence star? There have been suggestions in literature that HAEBE stars evolve into Vega-like stars ( eg. Malfait et al. 1998, Waters & Waelkens 1998). Vega-like stars are characterized by substantial far infrared excess due to cool dust, relatively low near infrared excess, low polarization and lack of emission lines in their spectra. The dust masses found around them are a few orders of magnitudes lower than that of HAEBEs. Also, Vega-like discs in general are gas depleted (Lagrange et al. 2000). Do all HAEBE stars pass through a Vega-like phase with gas depleted discs? These questions are critical to our understanding of the nature of the pre-main sequence evolution of intermediate mass stars.

A study of non-emission line young stellar candidate objects listed in Table 5 of the catalogue by Thé et al. (1994) may shed some light on the issues raised above. There are fourteen stars listed in this table. Typical PMS properties are less clearly seen in these stars. One of them, $`\beta `$ Pic, is a bonafide Vega-like star. The evolutionary status of other stars is not very clear. They are believed to be, as the authors suggest, transition objects between PMS and MS phase.

In this paper we present the results of a study of four non-emission line stars listed in the afore-mentioned catalogue. Spectroscopic and polarimetric observations of these stars were carried out. In Section 2 of this paper we present our observations. Together with the information available from literature in different wavelength ranges, we discuss the structure of circumstellar environment and the evolutionary status of these stars in Section 3. Summary of our study is presented in Section 4.

## 2 Observations

Medium resolution ($`\lambda /\mathrm{\Delta }\lambda 3000`$) optical CCD spectra were obtained for stars HD 36917 (V372 Ori), HD 36982 (LP Ori) and HD 37062 (V361 Ori) with the Optometrics Research (OMR) spectrograph on the 2.3 meter Vainu Bappu Telescope (VBT) and for HD 58647 with the Universal Astronomical Grating Spectrograph (UAGS) on the 1 meter telescope at the Vainu Bappu Observatory, Kavalur, India. Log of spectroscopic observations is given in Table 1. The prototype Herbig Ae/Be star $`ABAur`$ was also observed with the UAGS on the 1 meter telescope and is included in Table 1. All spectra were bias subtracted, flat-field corrected, extracted and wavelength calibrated in the standard manner using the IRAF111 IRAF is distributed by National Optical Astronomy Observatories, USA. reduction package. In view of the presence of surrounding diffuse H$`II`$ region nebulosity in the direction of three of the programme stars, the background sky subtraction was performed in the following way. The stellar spectra were extracted by summing up over ten pixels (plate scale = $`0.2\mathrm{}/pixel`$) perpendicular to the dispersion axis. The sky background subtracted was obtained by summing up ten pixels, ten pixels ($`2\mathrm{}`$) away from the star on either side. The spectra were corrected for instrumental response and brought to a relative flux scale using the spectrophotometric standard observed on the same night. Each spectrum spans a wavelength range of $`2400\AA `$, centered roughly at $`H\alpha (\lambda 6562\AA )`$. Reduced spectra of HD 36917, HD 36982 and HD 37062 along with that of $`ABAur`$ are presented in Figure 1(a) and that of HD 58647 in Figure 1(b). Optical linear polarization measurements were made with a fast star-and-sky chopping polarimeter (Jain & Srinivasulu 1991) coupled at the $`f/13`$ Cassegrain focus of the 1 meter telescope at the Vainu Bappu Observatory, Kavalur of the Indian Institute of Astrophysics. A dry-ice cooled R943-02 Hamamatsu photomultiplier tube was used as the detector. All measurements were made in the $`V`$ band with an aperture of $`15\mathrm{}`$ . The instrumental polarization was determined by observing unpolarized standard stars from Serkowski (1974). It was found to be $`0.1\%`$, and has been subtracted vectorially from the observed polarization of the programme stars. The zero of the polarization position angle was determined by observing the polarized standard stars from Hsu & Breger (1982). The position angle is measured from the celestial north, increasing eastward. Observed polarization and position angle values are presented in Table 2. Object name and date of observations are presented in column 1 and 2 of the table. Percentage polarization in V band and the probable error in the measurement are given in columns 3 and 4 and position angle and probable error in columns 5 and 6.

It can be seen from Figure 1(a) and 1(b) that emission lines are not present in the spectra of these stars except HD 58647. In Figure 1(a) we have also included a spectrum of $`ABAur`$, a prototype Herbig Ae star, in the same wavelength range and of similar resolution as of the other spectra for comparison. Typical HAEBE emission features such as $`H\alpha `$, HeI ($`\lambda 5875`$ & $`\lambda 6678`$) and OI ($`\lambda 7774`$)which are prominent in $`ABAur`$ are not seen in the spectra of HD 36917, HD 36982 and HD 37062. In contrast, HD 58647 has a strong $`H\alpha `$ line in emission but does not show any other characteristic HAEBE emission features. These stars cannot be unequivocally classified as HAEBE stars.

We note that three of the programme stars, viz., HD 36917, HD 36982, HD 37062 are towards the direction of the Orion Nebula and the diffuse H$`II`$ region present there is projected on to the line of sight to these stars. Therefore a careful subtraction of the surrounding nebular emission from the observed spectra is very important. This is illustrated in Figure 2 which presents the raw spectrum with superposed nebular lines and the reduced spectrum with surrounding nebular emission subtracted for HD 37062.

In Table 3 we present polarization data for these stars already existing in literature. The polarization values measured are typical of that of pre-main sequence stars except, again, for HD 58647 which shows a relatively low value of polarization. Also, LP Ori for which we have more than one polarization measurements shows variability which is again an indicator of the youth of the star.

## 3 Discussion

All the stars in our sample have earlier been treated as HAEBE stars by different authors (e.g. Malfait et al. 1998, Yudin 2000, Valenti et al. 2000) though in the catalogue by Thé et al. (1994) they are identified as non-emission line stars. Spectra presented in Figure 1 show that the characteristic HAEBE emission features ($`HeI`$ $`\lambda 5875`$ & $`\lambda 6678`$ and $`OI`$ $`\lambda 7774`$) are conspicuously absent in these stars. Their classification as HAEBE stars is, therefore, doubtful. A strong $`H\alpha `$ emission is present in HD 58647, but as argued later in this paper this star is likely to be a classical Be star. Weak $`H\alpha `$ emission, filling in the absorption core is present in HD 36917 even upon the subtraction of the possible contribution of the surrounding nebulosity. By comparing with synthetic spectra (Pickles 1985) of normal spectrophotometric standard stars of similar spectral type, we determine the equivalent width of the $`H\alpha `$ emission in this star to be $`2.5`$ $`\AA `$. This upper limit (for any residual nebular emission) to the $`H\alpha `$ equivalent width for HD 36917 is much smaller than those for typical Herbig Ae/Be stars (for the Herbig Ae/Be prototype star AB Aur we find $`H\alpha `$ emission equivalent width to be $`44`$ $`\AA `$). Interpreted as a measure of accretion rate (e.g. Hillenbrand et al. 1992, Muzerolle et al. 2001) this $`H\alpha `$ emission line flux implies mass accretion rate ($`\dot{M}`$) of order $`10^7`$ which is much smaller ( by a factor of 10 to 100) than those in typical Herbig Ae/Be stars ($`\dot{M}10^6M_{}/yr10^4M_{}/yr`$ depending on the spectral type; Hillenbrand et al. 1992).

Three of the stars that we have observed viz., HD 36917, HD 36982 and HD 37062 are in the direction of the known H$`II`$ region (Sharpless No. 281) towards the the Orion Nebula. From proper motion and radial velocity studies it has been established that these stars are kinematically connected with the Orion Nebula Cluster (ONC) which is at the northwestern end of Orion A cloud (van Altena et al. 1988, McNamara & Huels 1983, McNamara 1976, Tian et al. 1996, Hillenbrand 1997). They belong to the stellar population in the inner $`2.5`$ $`pc`$ of the ONC which is at a distance of $`470\pm 70`$ $`pc`$ (Genzel et al. 1981; Walker 1969; Tian et al. 1996; Hillenbrand 1997). ONC as a whole is characterized by a mean age of $`<1Myr`$ and an age spread which is probably less than 2 Myr with an internal velocity dispersion of $`2`$ $`kms^1`$ (Hillenbrand 1997; Tian et al. 1996). The kinematic association with the Orion Nebula Cluster strongly constrains the ages of these stars to within a few $`10^6`$ years. They are young stars of intermediate mass. In the following we use the average distance to the ONC $`(470`$ $`pc)`$ as the distance to its member stars.

To complement our observations, we have collated available information from literature on our sample stars. We present them in Table 4. The spectral types are taken from SIMBAD except that of HD 36917 for which it is from Lavato & Abt (1976). $`vsini`$ values are taken form Yudin (2000) and van den Ancker et al. (1998). The optical B and V magnitudes are derived from Tycho magnitudes except for HD 37062 for which B and V magnitudes are from SIMBAD. Tycho photometric errors are less than 0.02 mag. The J, H, K magnitudes are from 2MASS catalogue. Errors in J,H,K magnitudes are less than 0.07 mag. Reddening,$`E(BV)`$ estimated from spectral types and photometric magnitudes is listed in column 9.

We have constructed near infrared (NIR) colour-colour diagram from the 2MASS magnitudes for our programme stars which is shown in Figure 3. Along with the four programme stars, HAEBE stars and main sequence stars are also plotted in the diagram. The colours for the main sequence stars are from Koornneef (1983). HAEBE stars are taken from Thé et al. (1994) and their colours are derived from 2MASS magnitudes. The two parallel dotted lines form reddening band for normal stellar photospheres. These lines are parallel to the reddening vector and bound the range in the colour-colour diagram within which stars with purely reddened normal stellar photospheres can fall (Lada & Adams 1992). It can be seen that all the four stars are distinctly separated from the region occupied by HAEBE stars. Their near infrared charcterestics are different from that of HAEBE stars. The near infrared excesses of all the four stars are considerably lower than that of HAEBE stars. HD 36982 and HD 37062 have very little near infrared excess, if any. The near infrared excess in HAEBE stars is attributed to reradiation from hot dust less than $`1AU`$ from the star. The low near infrared excess shown by our programme stars, then, would strongly suggest the absence of submicron sized dust grains close to the star.

In Table 5 we present far infrared data for the stars and the quantities estimated from them except for HD 37062 which does not have an IRAS entry. HD number of the stars, their IRAS source names and the IRAS flux densities at 12, 25, 60 and 100$`\mu `$ are listed in the first six columns. In column 7 we list the fractional infrared luminosity ($`L_{ir}/L_{}`$) estimated from IRAS fluxes. $`L_{ir}=4\pi d^2F_{ir}`$ with
$`F_{ir}=[20.653f_{12}+7.53f_{25}+4.578f_{60}+1.762f_{100}]\times 10^{14}Wm^2`$
and $`d`$ being the distance to the star (Cox 2000). The IRAS flux densities (Jy) at 12, 25, 60 and 100 microns are given by $`f_{12}`$, $`f_{25}`$, $`f_{60}`$, and $`f_{100}`$ respectively. $`L_{}`$ is computed from $`M_v`$ using standard bolometric corrections (Cox 2000) where $`M_v`$ is evaluated from dereddened $`V`$ magnitude and distance to the star. Since the star HD 36982 is located below the ZAMS in the colour - magnitude diagram (cf. Fig. 4) $`L_{}`$ for this star is taken to be the ZAMS luminosity expected for it’s spectral type. The black body colour temperatures derived from the ratio of fluxes at $`25\mu `$ and $`60\mu `$ are considered to be the dust temperatures and are listed in column 8. The dust temperature derived for HD 36917 is a lower limit since the IRAS $`60\mu `$ flux density is only an upper limit. Dust masses (in units of $`M_{}`$, the mass of earth = $`6\times 10^{27}g`$) listed in column 9 are computed using the relation,

$`M_d=4\pi a\rho _dd^2F_{ir}/3Q_\nu \sigma T_d^4`$

assuming a grain size of $`a=1\mu `$, dust grain material density $`\rho _d=2g/cc`$, absorption efficiency $`Q_\nu =0.5`$, and with $`F_{ir}`$ computed from the flux densities at four IRAS bands. Since IRAS flux densities at $`100\mu `$ are upper limits for all the stars we have used flux densities expected at $`100\mu `$ for the derived colour temperatures in estimating $`F_{ir}`$.

The fractional infrared luminosities $`L_{ir}/L_{}`$ estimated for HD 36917 and HD 36982 are quite significant ($`L_{ir}/L_{}`$ $`0.15`$). Here it is assumed that the IRAS flux densities quoted do represent emission from these sources and are not due to other sources in the IRAS beam. We note that IRAS point sources 05323-0536, 05327-0529 and 07236-1404 have positional coincidences with HD 36917 , HD 36982 and HD 58647 to within $`4\mathrm{}`$, $`2\mathrm{}`$ and $`2\mathrm{}`$ respectively. The surface density of IRAS point sources in this region is $`3\times 10^3`$ to $`10^2sources/(arcmin)^2`$. Therefore, for a given star, in the 12 and 25 $`\mu `$ bands, the probability that an unrelated IRAS point source is in the IRAS beam ($`0.75\mathrm{}\times 4.5\mathrm{}`$) is only $`13\%`$. If the IRAS fluxes observed for these stars are dominated by thermal emission from the circumstellar dust present around these stars, which is most likely the case in view of the relatively low surface density of IRAS point sources and low probability of chance projections, then for the derived dust temperatures listed in Table 5 the emitting dust is at distances of $`10AU`$ (HD 36917) and $`100AU`$ (HD 36982) from the central star, assuming the dust to be distributed in an optically thick disc as in some models of HAEBE stars (Hillenbrand et al. 1992). On the other hand if the dust is distributed in optically thin shells as in Vega-like stars, then the dust is located at $`100AU`$ (HD 36917) and $`1000AU`$ (HD 36982) from the stars. These dimensions of dust shells or rings are quite similar to those of Vega-like stars and are compatible with the dust responsible for emission being part of the circumstellar environment. It is clear from Table 5 that these two stars have significant amounts of cold circumstellar dust present around them. The dust masses estimated are systematically lower than those for HAEBE stars but higher than that found around Vega-like stars.

In Figure 4 we present the far infrared spectral energy distribution for the three stars. Their IRAS flux densities are plotted against the wavelength together with the expected photospheric flux densities at IRAS wavelengths (Song 2001) for a star of that spectral type. It is clearly seen that HD 36917 and HD 36982 have considerable excess at far infrared wavelengths whereas HD 58647 has a low excess and an energy distribution that is decreasing with wavelength.

It is clear from the above discussion that the three member stars of Orion Nebula Cluster (ONC) are extremely young and lack emission lines in their spectra. This would imply the absence of a hot and extended chromosphere around them where the emission lines are thought to be formed. In the case of HD 36982 and HD 37062, which are early B type stars it is possible that this is an evolutionary effect. They are possibly at the end of their pre-main sequence phase, which is roughly $`12`$ $`Myr`$. The absence of $`HeI`$ line in emission, which is an indicator of accretion, also suggests that these stars are beyond their accretion phase. Almost complete absence of near infrared excess further supports the absence of an inner accretion disc. This evolutionary picture is further supported by the fact that HD 36917, which is of a later spectral type and thus has a longer PMS life time, shows near infrared excess though at a much lower level than seen in HAEBEs. Binarity of the star cannot account for the near infrared excess. The $`JH`$ and the $`HK`$ colours computed by adding up the individual fluxes expected for each component of the binary ($`B9.5+A0.5`$) and estimating the combined magnitudes, are only $`0.006`$ and $`0.002`$ whereas the observed colours are 0.257 and 0.397 respectively. Thus it could be concluded that the inner disc has not been completely disrupted in this star. .

However, we do not rule out the possiblity of the stars not passing through a HAEBE phase with emission lines and near infrared excess. A different formation mechanism or a very destructive cluster environment can drastically alter the PMS properties and evolutionary sequence that a young star passes through.

In any case these stars fit well with the definition of Vega-like stars though the far infrared excesses, fractional infrared luminosities $`(L_{ir}/L_{})`$, and the dust masses computed for HD 36982 and HD 36917 are much higher than that for the prototype Vega-like stars and for ”old PMS” (OPMS) and ”young main sequence” (YMS) systems discussed recently by Lagrange et al. (2000). Their circumstellar dust may not be the debris product but rather what is left over from their PMS phase. These stars, then, are the youngest Vega-like stars hitherto known.

Unlike the stars which are members of the Orion Nebula Cluster, HD 58647 may not be a young star. This object is about 30 degree away from the Orion complex and is not associated with any star forming cloud. Though it shows an excess at the near infrared wavelengths, it’s far infrared excess is small. $`L_{ir}/L_{}`$ indicates the scarcity of dust present around the star. Low polarization observed reinforces this fact. Colour temperature derived from the IRAS fluxes is very high and is much larger than the dust sublimation temperature. The dust mass estimated formally from IRAS fluxes is negligibly low and is unphysical as the colour temperature does not represent the dust temperature. Clearly, the excess shown by the star in the IRAS wavelengths cannot be attributed to circumstellar dust. The infrared excess, both near and far, most likely is due to free-free emission from an ionized region around the star. We suggest that HD 58647 is a classical Be star. The relatively large value of $`vsini280`$ $`kms^1`$ supports this claim. The presence of $`H\alpha `$ line in emission observed in the spectra is consistent with the star being a classical Be. Moreover, its IRAS flux density distribution varies nearly as $`F_\nu `$ $`\alpha `$ $`\nu ^{0.6}`$ which is very similar to that for classical Be stars (Taylor 1990).

A colour-absolute magnitude diagram constructed for all the four programme stars is shown in Figure 5. The zero-age main sequence data is taken from Schmidt-Kaler (1965). For the ONC member stars a distance of $`470`$ $`pc`$ is assumed and a Hipparcos distance of $`277`$ $`pc`$ is used for HD 58647 in computing the absolute magnitude. Major contribution to the error bars shown in the figure results from uncertainities in distance. Further, the average interstellar value of 3.1 is used for the ratio $`R=A_v/E(BV)`$ . Significance of their position in the colour-magnitude diagram (CMD) is discussed below.

HD 36917 which is a spectroscopic binary (B9.5 + A0.5) (Lavato & Abt 1976) is found to be way above the main sequence in the colour - magnitude diagram. Increase in the brightness caused by binarity would only account for $`0.7mag`$. Thus it’s location in the colour - magnitude diagram is consistent with the star being a pre-main sequence and the age indicated by its kinematic association with the ONC.

The other two stars which are members of the ONC, viz., HD 36982 and HD 37062 fall below the zero-age main sequence (in particular HD 36982) in the CMD. A possible explanation for their anomalous position is that they have an anomalous circumstellar extinction component. If they have already reached the zero age main sequence then the extinction towards HD 36982 and HD 37062, implied by their location in the CMD and computed from the observed and absolute V magnitudes for a distance of $`470`$ $`pc`$ are $`2.12`$ $`mag`$ and $`0.7`$ $`mag`$ respectively. Such large extinction and a rather low $`E(BV)`$ cannot be produced by interstellar grains of submicron size. This would imply the presence of large grains around these stars. The neutral extinction produced by these grains can be even larger than that estimated assuming the stars to be zero age main sequence stars, without causing any change in the colour excess. The possiblity of these stars being intrinsically brighter than a ZAMS star of similar spectral type cannot be ruled out. They could be pre-main sequence stars with very high extinction.

## 4 Summary

We have obtained optical spectra of four non-emission line stars listed in the catalogue by Thé et al. (1994). Emission lines which characterize the typical HAEBE spectrum are not seen in any of the stars. Also, their near infrared properties are very unlike HAEBE stars. The amount of circumstellar dust present is also systematically lower than that found around HAEBE stars. We argue that these stars cannot be unambigu

ously classified as HAEBE stars. Nevertheless, three of the stars which are kinematic members of Orion Nebula Cluster are very young. This association constrains their ages to be less than a few $`Myr`$. Absence of emission lines, low near infrared excesses and presence of far infrared excesses in these objects make them somewhat similar to Vega-like stars though the dust masses estimated for them are much higher than that of prototype Vega-like stars. These stars may then represent the youngest examples of Vega phenomenon and may well be the intermediate mass counter parts of weak line T Tauri stars. The observed reddening and estimated extinction for these stars indicate presence of larger than submicron sized dust grains around these stars. One of the stars in our sample, HD 58647, is more likely to be a classical Be star as evident from the low $`L_{ir}/L_{}`$, scarcity of circumstellar dust, low polarization, presence of $`H\alpha `$ emission and appreciable near infrared excess and far IR spectral energy distribution consistent with free-free emission similar to other well known classical Be stars.

ACKNOWLEDGMENTS
We thank the referee Dr. M. van den Ancker for critical comments and useful suggestions.


text:
# Mean field dynamics of fermions and the time-dependent Hartree-Fock equation

## 1 Introduction

In this article we consider the Hamiltonian dynamics of systems of fermions and derive the time-dependent Hartree-Fock equation in the mean field limit. We follow the approach of Spohn, who derived a mean field dynamical equation (the time-dependent Hartree equation) for mean field systems of distinguishable particles, remarking at the time that “the convergence of the mean field limit with statistics included is an open problem”.

In Spohn’s theory the initial $`N`$-body density operator $`D_N`$ is assumed to be a product state $`D^N`$, i.e., the particles are statistically independent and identically distributed. The mean field limit is investigated in the Schrödinger picture, where $`D_N(t)`$ obeys the von Neumann equation

$$i\mathrm{}\frac{d}{dt}D_N(t)=\underset{1jN}{}[L_j,D_N(t)]+\frac{1}{N}\underset{1i<jN}{}[V_{ij},D_N(t)]$$
(1)
with $`V_{ij}`$ denoting the two-body potential $`V`$ acting between the $`i^{th}`$ and $`j^{th}`$ particles and $`[,]`$ denoting the commutator. The limit as $`N\mathrm{}`$ of the $`n`$-body density operator $`D_{N:n}(t)`$ is shown to converge to $`D(t)^n`$, where $`D(t)`$ obeys a time-dependent Hartree equation. (The subscript :n appearing in $`D_{N:n}`$ is our notation for the $`n^{th}`$ partial trace, defined in equation (4) below.) Spohn’s ideas have been generalized to open systems in . There are other theories of quantum mean field dynamics, e.g., the algebraic theory of , but to our knowledge the problem of including quantum statistical effects remains unsolved.

The problem is that Fermi-Dirac or Bose-Einstein statistics constrain the possible initial condition of (1) to have the appropriate symmetry, which is typically inconsistent with the product form $`D^N`$. An $`N`$-body density operator with Fermi-Dirac symmetry can never have the form $`D^N`$ and a Bose-Einstein density operator can only have the form $`D^N`$ if $`D`$ is a pure state (i.e., if the system of bosons is in a condensed state). The remedy for this problem, for fermions, is to replace the hypothesis that the initial state be a product state with a hypothesis that is consistent with Fermi-Dirac statistics, e.g., that the initial states are Slater determinants.

The role of the factorization hypothesis $`D_N(0)=D^N`$ is to permit the closure of the BBGKY hierarchy by setting the two-body state $`D_{N:2}`$ equal to $`DD`$. Closing the hierarchy this way results in the time-dependent Hartree equation. This kind of closure hypothesis is implicit in the Stosszahlansatz that leads to Boltzmann’s kinetic equation for gases . Kac noted that, for Boltzmann’s equation, the factorization $`f_{N:2}=ff`$ is only realized in the limit $`N\mathrm{}`$, and he called this behavior the Boltzmann property . Later authors  developed Kac’s ideas; what is now called the propagation of chaos is an important tool in rigorous kinetic theory . We have noted that Boltzmann’s closure Ansatz is inconsistent with the Pauli Exclusion Principle, and needs to be replaced by another closure Ansatz when the particles are fermions. The novelty of our approach consists in replacing the condition of asymptotic independence of the particles by a condition that describes the correlations of Slater determinants. This condition, called Slater closure is defined in Definition 2.1 below.

Assuming that $`\{D_N(0)\}`$ is a sequence of initial states for (1) that has Slater closure, we can prove that $`\{D_N(t)\}`$ has Slater closure for all $`t>0`$. This phenomenon could be called the propagation of Slater closure because it is like the “propagation of chaos” mentioned above. Since $`\{D_N(t)\}`$ has Slater closure, the two-body density operator $`D_{N:2}(t)`$ is approximately equal to $`(D_{N:1}(t)D_{N:1}(t))\mathrm{\Sigma }_2`$ when $`N`$ is large, where $`\mathrm{\Sigma }_2`$ is the two-body operator defined by

$$\mathrm{\Sigma }_2(xy)=xyyx.$$

Substituting $`(D_{N:1}(t)D_{N:1}(t))\mathrm{\Sigma }_2`$ for $`D_{N:2}(t)`$ in the BBGKY hierarchy leads one to conjecture that, when $`N`$ is large, the single-body density operator should nearly obey the time-dependent Hartree-Fock (TDHF) equation

$`i\mathrm{}{\displaystyle \frac{d}{dt}}F(t)`$ $`=`$ $`[L,F(t)]+[V,(F(t)F(t))\mathrm{\Sigma }_2]_{:1}`$
$`F(0)`$ $`=`$ $`D_{N:1}(0).`$
Theorem 3.1 confirms this conjecture.

Theorem 3.1 states that the distance in the trace norm between $`D_{N:1}(t)`$ and the corresponding solution $`F(t)`$ of the TDHF equation tends to $`0`$ as $`N`$ tends to infinity. The trace norms of $`D_{N:1}(t)`$ and $`F(t)`$ are separately equal to $`1`$, so it is significant that their difference $`D_{N:1}(t)F(t)`$ converges to $`0`$ in the trace norm. A crucial detail of the proof is Lemma 5.1, which states that the operator norm of $`D_{N:1}`$ tends to $`0`$ if $`\{D_N\}`$ has Slater closure. Much of the rest of the proof lies in bounding the trace norm of $`D_{N:1}(t)F(t)`$ by an expression involving the operator norm of $`D_{N:1}(0)`$.

The use of the trace norm to measure the distance between two density operators is quite natural. A density operator $`D`$ corresponds to a quantum state through the assignment $`B\text{Tr}(DB)`$ of expectation values to bounded observables $`B`$. Thus, two density operators $`D`$ and $`D^{}`$ are within $`\epsilon `$ of one another in the trace norm if and only if they correspond to quantum states that give expectations differing by no more than $`\epsilon `$ for all observables $`B`$ with $`B1`$.

In this article, we assume that the two-body potential $`V`$ is a bounded operator. Recent work on the time-dependent Schrödinger-Poisson equation  suggests that it may be possible to prove a theorem similar to our Theorem 3.1 when $`V`$ is the Coulomb potential. This work shall be published in a separate paper.

The rest of this article is organized as follows: the next section discusses fermionic density operators and defines Slater closure. The N-particle Hamiltonian and the associated time-dependendent Hartree-Fock equation are described in Section 3. That section concludes with the statement of our main result, Theorem 3.1, whose proof spans Sections 4, 5, and 6. Section 7 is an appendix relating the operator form of the TDHF equation, used throughout this paper, to the formulation of the TDHF equation as a coupled system of wave equations, which may be more familiar to some readers.

## 2 Fermionic density operators and Slater closure

Let $`IH`$ be a Hilbert space, supposed to be the space of wavefunctions for a certain type of quantum system (a “component” or “particle”). Then the Hilbert space of wavefunctions for a system consisting of $`N`$ distinguishable components or particles of that type is $`IH_N=IH^N`$. If the components are not distinguishable, but obey Fermi-Dirac statistics, then the appropriate Hilbert space of wavefunctions is the antisymmetric subspace $`𝒜_NIH_N`$. To define this subspace, it is convenient first to define unitary transposition and permutation operators on $`IH_N`$. The transposition operator $`U_{(ij)}`$ is defined by extending the following isometry defined on simple tensors

$$U_{(ij)}(x_1x_2\mathrm{}x_i\mathrm{}x_j\mathrm{}x_N)=x_1x_2\mathrm{}x_j\mathrm{}x_i\mathrm{}x_N$$

to all of $`IH_N`$. For any $`\pi `$ in the group $`\mathrm{\Pi }_N`$ of permutations of $`\{1,2,\mathrm{},N\}`$, one may define the permutation operator $`U_\pi `$ as $`U_{(i_kj_k)}\mathrm{}U_{(i_2j_2)}U_{(i_1j_1)}`$, where $`(i_kj_k)\mathrm{}(i_2j_2)(i_1j_1)`$ is any product of transpositions that equals $`\pi `$.

The antisymmetric subspace may now be defined as

$$𝒜_N=\{\psi IH_N:U_\pi \psi =\text{sgn}(\pi )\psi \pi \mathrm{\Pi }_N\}.$$

One may verify that

$$P_{𝒜_N}=\frac{1}{N!}\underset{\pi \mathrm{\Pi }_N}{}\text{sgn}(\pi )U_\pi $$

is the orthogonal projector whose range is $`𝒜_N`$.

The pure states of an $`N`$-fermion system correspond to the orthogonal projectors $`P_\psi `$ onto one-dimensional subspaces of $`𝒜_N`$. That is, a pure state is given by

$$P_\psi (\varphi )=\psi ,\varphi \psi $$

for some $`\psi 𝒜_N`$ of unit length. The statistical states of the $`N`$-fermion system are the positive trace class operators or density operators $`D`$ on $`𝒜_N`$ of trace $`1`$. These can be identified with density operators $`D`$ on all of $`IH_N`$ whose eigenvectors lie in $`𝒜_N`$, i.e., such that

$$D=\underset{i=1}{\overset{\mathrm{}}{}}\lambda _iP_{\psi _i}$$

for some orthonormal system $`\{\psi _i\}`$ in $`𝒜_N`$ and a family of positive numbers $`\lambda _i`$ that sum to $`1`$. It follows that these fermionic densities are those density operators that satisfy

$$DU_\pi =U_\pi D=\text{sgn}(\pi )D\pi \mathrm{\Pi }_N.$$
(2)
If a density operator $`D`$ on $`IH_N`$ commutes with every permutation operator $`U_\pi `$ then it is symmetric. In particular, fermionic densities are symmetric by (2).

If $`\{e_j\}_{jJ}`$ is an orthonormal basis of $`IH`$ then

$$\{e_{j_1}e_{j_2}\mathrm{}e_{j_N}:j_1,j_2,\mathrm{},j_NJ\}$$

is an orthonormal basis of $`IH_N`$. Since $`𝒜_N`$ is the range of $`P_{𝒜_N}`$ and since

$$P_{𝒜_N}(e_{j_1}e_{j_2}\mathrm{}e_{j_N})=0$$

unless all of the indices $`j_i`$ are distinct, the set

$$\{P_{𝒜_N}(e_{j_1}e_{j_2}\mathrm{}e_{j_N}):j_1,j_2,\mathrm{},j_n\text{all distinct}\}$$

is a spanning set for $`𝒜_N`$. In fact it is an orthogonal basis for $`𝒜_N`$, each vector having norm $`1/\sqrt{N!}`$. Vectors of the form $`\sqrt{N!}P_{𝒜_N}(e_{j_1}e_{j_2}\mathrm{}e_{j_N})`$ are known as Slater determinants.

The trace class operators on a Hilbert space $`IH`$ form a Banach space $`𝒯(IH)`$ with the norm $`T_{tr}=\text{Tr}(|T|)`$. The important inequality

$$TB_{tr}T_{tr}B$$
(3)
holds whenever $`B`$ is a bounded operator of norm $`B`$ and $`T𝒯(IH)`$. It is this basic inequality that will produce our key estimates.

For $`nN`$, the $`n^{th}`$ partial trace is a contraction from $`𝒯(IH^N)`$ onto $`𝒯(IH^n)`$. The $`n^{th}`$ partial trace of $`T`$ will be denoted $`T_{:n}`$, and may be defined as follows: Let $`𝒪`$ be any orthonormal basis of $`IH`$. If $`T𝒯(IH^N)`$ and $`n<N`$ then

$`T_{:n}(w),x=`$
$`{\displaystyle \underset{z_1,\mathrm{},z_{Nn}𝒪}{}}`$ $`<T(wz_1\mathrm{}z_{Nn}),xz_1\mathrm{}z_{Nn}>`$ (4)
for any $`w,xIH^n`$. If a trace class operator $`T𝒯(IH^N)`$ satisfies (2) then so does $`T_{:n}`$, i.e., the partial trace defines a positive contraction from $`𝒯(IH^N)`$ to $`𝒯(IH^n)`$ that carries fermionic densities to fermionic densities.

In the following definition, and throughout this article, we use the superscript ⊗n to denote the $`n^{th}`$ tensor power of an operator, and we use the notation $`\mathrm{\Sigma }_n`$ for $`n!P_{𝒜_n}`$, i.e.,

$$\mathrm{\Sigma }_n=\underset{\pi \mathrm{\Pi }_n}{}\text{sgn}(\pi )U_\pi .$$

The $`n^{th}`$ tensor power of an operator $`A`$ on $`IH`$ is the operator $`A^n`$ on $`IH_n`$ defined on simple tensors by

$$A^n(x_1x_2\mathrm{}x_n)=Ax_1Ax_2\mathrm{}Ax_n.$$

###### Definition 2.1

For each $`N`$, let $`D_N`$ be a symmetric density operator on $`IH_N`$. The sequence $`\left\{D_N\right\}`$ has Slater closure if, for each fixed $`n`$,

$$\underset{N\mathrm{}}{lim}D_{N:n}D_{N:1}^n\mathrm{\Sigma }_n_{tr}=0.$$

This terminology is motivated by the observation that, if $`\mathrm{\Psi }_N`$ is a Slater determinant in $`𝒜_N`$ and $`P_{\mathrm{\Psi }_N}`$ denotes the orthogonal projector onto the span of $`\mathrm{\Psi }_N`$ then

$$\left(P_{\mathrm{\Psi }_N}\right)_{:n}=\frac{N^n(Nn)!}{N!}\left(P_{\mathrm{\Psi }_N}\right)_{:1}^n\underset{\pi \mathrm{\Pi }_n}{}\text{sgn}(\pi )U_\pi ,$$
(5)
for this implies the following:

###### Proposition 2.1

For each $`N`$ let $`\mathrm{\Psi }_N`$ be a Slater determinant in $`𝒜_N`$, and let $`P_{\mathrm{\Psi }_N}`$ denote the orthoprojector onto the span of $`\mathrm{\Psi }_N`$. Then $`\left\{P_{\mathrm{\Psi }_N}\right\}`$ has Slater closure.

## 3 The time-dependent Hartree-Fock equation

We are going to prove that, in the mean field limit, the time-dependent Hartree-Fock equation describes the time-evolution of the single-particle state in systems of fermions. We state our theorem in this section and go on to prove it in the three subsequent sections.

First we describe the $`N`$-particle Hamiltonian. Let $`L^{(N)}`$ be a self-adjoint operator on $`IH`$, where $`L^{(N)}`$ may depend on $`N`$ in an arbitrary manner. The free motion of the $`j^{th}`$ particle is governed by

$$L_j^{(N)}=I^{j1}L^{(N)}I^{Nj},$$

where $`I`$ denotes the identity operator on $`IH`$. The interaction between the particles has the form $`1/N`$ times the sum over pairs of distinct particles of a two-body potential $`V`$. Let $`V`$ be a bounded Hermitian operator on $`IHIH`$ that commutes with the transposition operator $`U_{(12)}`$. Define the operator $`V_{12}`$ on $`IH_N`$ by

$$V_{12}\left(x_1x_2\mathrm{}x_N\right)=V(x_1x_2)x_3\mathrm{}x_N$$

and for each $`1i<jN`$ define $`V_{ij}=U_\pi ^{}V_{12}U_\pi `$ where $`\pi `$ is any permutation with $`\pi (i)=1`$ and $`\pi (j)=2`$. Let

$$H_N=\underset{1jN}{}L_j^{(N)}+\frac{1}{N}\underset{1i<jN}{}V_{ij}$$
(6)
be the N-particle Hamiltonian operator on $`IH_N`$. The von Neumann equation for the $`N`$-particle density operator $`D_N(t)`$ is

$$i\mathrm{}\frac{d}{dt}D_N(t)=\underset{1jN}{}[L_j^{(N)},D_N(t)]+\frac{1}{N}\underset{1i<jN}{}[V_{ij},D_N(t)]$$
(7)
and has the solution

$$e^{iH_Nt/\mathrm{}}D_N(0)e^{iH_Nt/\mathrm{}}.$$
(8)
Next we define the time-dependent Hartree-Fock equation. Let $`L^{(N)}`$ and $`V`$ be as above, and let $`U_{(12)}`$ denote the transposition operator on $`IHIH`$. The time-dependent Hartree-Fock (TDHF) equation for a density operator $`F(t)`$ on $`IH`$ is

$`i\mathrm{}{\displaystyle \frac{d}{dt}}F(t)`$ $`=`$ $`[L^{(N)},F(t)]+[V,F_2^{}(t)]_{:1}`$
$`F_2^{}(t)`$ $`=`$ $`(F(t)F(t))\left(IU_{(12)}\right)=F(t)^2\mathrm{\Sigma }_2`$ (9)
(the subscript :1 on the last commutator denotes partial contraction). Following , we define a mild solution of equation (9) to be a continuous function $`tF(t)`$ from $`[0,\mathrm{})`$ to the real Banach space of Hermitian trace class operators such that

$$F(t)=e^{i\frac{t}{\mathrm{}}L^{(N)}}F(0)e^{+i\frac{t}{\mathrm{}}L^{(N)}}\frac{i}{\mathrm{}}_0^te^{i\frac{ts}{\mathrm{}}L^{(N)}}[V,F_2^{}(s)]_{:1}e^{+i\frac{ts}{\mathrm{}}L^{(N)}}𝑑s$$

for all $`t0`$. The results proved in  show that (9) has a global mild solution555The solution obtained in Theorem 4.2 of  is indeed defined for all positive times because the nonlinearity of TDHF satisfies condition (4.1) of  — see Proposition 3.5 there. for each Hermitian trace class initial data $`F(0)`$. Furthermore,

$$F(t)=U^{}F(0)U$$
(10)
for some unitary operator depending on $`t`$ and $`F(0)`$. In particular, the operator norm of $`F(t)`$ is constant.

The relationship between the $`N`$-particle system and the TDHF equation is the subject of our main theorem. Recall the definition (2.1) of Slater closure.

###### Theorem 3.1

For each $`N`$, let $`D_N(t)`$ the solution to (7) whose initial value $`D_N(0)`$ is a symmetric density. Let $`F^{(N)}(t)`$ be the mild solution of the TDHF equation (9) whose initial value is $`F^{(N)}(0)=D_{N:1}(0)`$.

If $`\left\{D_N(0)\right\}`$ has Slater closure then $`\left\{D_N(t)\right\}`$ has Slater closure and

$$\underset{N\mathrm{}}{lim}D_{N:1}(t)F^{(N)}(t)_{tr}=0$$

for all $`t>0`$.

## 4 Two hierarchies, and their difference

Consider the $`N`$-particle von Neumann equation (7). From now on we will suppose that the initial $`N`$-particle density operator $`D_N(0)`$ is symmetric, i.e., that

$$U_\pi ^{}D_N(0)U_\pi =D_N(0)$$

for all $`\pi \mathrm{\Pi }_N`$. (Recall that, in particular, fermionic densities are symmetric.) The symmetry of the Hamiltonian (6) ensures that $`D_N(t)`$ remains symmetric for all $`t`$. From (7) and the symmetry of $`D_N(t)`$, it follows that the partial trace $`D_{N:n}(t)`$ satisfies666In the proof of Theorem 3.1, we proceed as if $`D_N`$ was a classical solution of equation (7); the case of general initial data as in the statement of the Theorem follows from a density argument and the unitarity of the group generated by the $`N`$-body Schrödinger operator.

$`i\mathrm{}{\displaystyle \frac{d}{dt}}D_{N:n}(t)`$ $`=`$ $`{\displaystyle \underset{1jn}{}}[L_j^{(N)},D_{N:n}(t)]+{\displaystyle \frac{1}{N}}{\displaystyle \underset{1i<jn}{}}[V_{ij},D_{N:n}(t)]`$ (11)
$`+{\displaystyle \frac{Nn}{N}}{\displaystyle \underset{1in}{}}[V_{i,n+1},D_{N:n+1}(t)]_{:n}`$
The system of equations (11) for $`D_{N:1},D_{N:2},\mathrm{},D_{N:N1}`$ together with the equation (7) for $`D_N`$ is called the N-particle hierarchy. For our estimates later on, it is convenient to rewrite the equations (11) of the hierarchy as

$`i\mathrm{}{\displaystyle \frac{d}{dt}}D_{N:n}(t)`$ $`=`$ $`_n^{(N)}(D_{N:n}(t))+{\displaystyle \underset{1in}{}}[V_{i,n+1},D_{N:n+1}(t)]_{:n}`$ (12)
$`+_n(t,N,D_N(0))`$
with

$`_n^{(N)}()`$ $`=`$ $`{\displaystyle \underset{1jn}{}}[L_j^{(N)},]`$
$`_n(t,N,D_N(0))`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{1i<jn}{}}[V_{ij},D_{N:n}(t)]`$ (13)
$`{\displaystyle \frac{n}{N}}{\displaystyle \underset{1in}{}}[V_{i,n+1},D_{N:n+1}(t)]_{:n}.`$
Next we describe another hierarchy, built from “the bottom up” out of solutions to the TDHF equation, in contrast to the hierarchy we have just considered, which is built from “the top down” starting with solutions to (7). If $`F`$ is a trace class operator, define $`F_1^{}=F`$ and

$$F_n^{}=F^n\mathrm{\Sigma }_n$$

for $`n>1`$. When $`F`$ depends on $`t`$ we write $`F_n^{}(t)`$ instead of $`F(t)_n^{}`$. The notation $`F_2^{}(t)`$ has already been used in the TDHF equation (9).

###### Proposition 4.1

If $`F`$ is a classical solution 777Consistently with the previous footnote, we are proceeding as if $`F_N`$ was a classical solution of (9); the case of general initial data as in the statement of the Theorem is recovered by a density argument and the continuous dependence on initial data of the mild solution to (9): see Theorem 4.1 of . of the TDHF equation (9) then

$$i\mathrm{}\frac{d}{dt}F_n^{}(t)=\underset{j=1}{\overset{n}{}}[L_j^{(N)},F_n^{}(t)]+\underset{j=1}{\overset{n}{}}[V_{j,n+1},F_{n+1}^{}(t)]_{:n}+_n(F(t))$$

where $`_n`$ is defined on trace class operators by $`_1(X)=\mathrm{𝟎}`$ (the zero operator) and

$$_n(X)=\underset{j=1}{\overset{n}{}}[V_{j,n+1},X^{n+1}\underset{kj}{}U_{(k,n+1)}]_{:n}\mathrm{\Sigma }_n$$
(14)
for $`n>1`$.

Proof: For any trace class operator $`X`$,

$`{\displaystyle \underset{j=1}{\overset{n}{}}}[V_{j,n+1},X_{n+1}^{}]_{:n}`$
$`=`$ $`{\displaystyle \underset{j=1}{\overset{n}{}}}[V_{j,n+1},X^{n+1}\left(I{\displaystyle \underset{k=1}{\overset{n}{}}}U_{(k,n+1)}\right)\mathrm{\Sigma }_nI_{(IH)}]_{:n}`$
$`=`$ $`{\displaystyle \underset{j=1}{\overset{n}{}}}[V_{j,n+1},X^{n+1}\left(I{\displaystyle \underset{k=1}{\overset{n}{}}}U_{(k,n+1)}\right)]_{:n}\mathrm{\Sigma }_n.`$
The first equality in (4) holds because

$$\mathrm{\Sigma }_{n+1}=\left(I\underset{k=1}{\overset{n}{}}U_{(k,n+1)}\right)\mathrm{\Sigma }_nI_{(IH)},$$

and the second equality in (4) holds because $`\mathrm{\Sigma }_nI_{(IH)}`$ commutes with $`\underset{j=1}{\overset{n}{}}V_{j,n+1}`$. From the TDHF equation (9) we calculate that

$`i\mathrm{}{\displaystyle \frac{d}{dt}}F_n^{}(t)=i\mathrm{}{\displaystyle \frac{d}{dt}}F(t)^n\mathrm{\Sigma }_n`$
$`=`$ $`i\mathrm{}\left\{{\displaystyle \underset{j=1}{\overset{n}{}}}F(t)^{j1}{\displaystyle \frac{d}{dt}}F(t)F(t)^{nj}\right\}\mathrm{\Sigma }_n`$
$`=`$ $`{\displaystyle \underset{j=1}{\overset{n}{}}}\left\{[L_j^{(N)},F(t)^n]+[V_{j,n+1},F(t)^{n+1}\left(IU_{(j,n+1)}\right)]_{:n}\right\}\mathrm{\Sigma }_n`$
$`=`$ $`{\displaystyle \underset{j=1}{\overset{n}{}}}[L_j^{(N)},F_n^{}(t)]+_n(F(t))`$
$`+{\displaystyle \underset{j=1}{\overset{n}{}}}[V_{j,n+1},F(t)^{n+1}\left(I{\displaystyle \underset{k=1}{\overset{n}{}}}U_{(k,n+1)}\right)]_{:n}\mathrm{\Sigma }_n.`$
By the identity (4), the last sum in (4) equals $`\underset{j=1}{\overset{n}{}}[V_{j,n+1},F_{n+1}^{}(t)]_{:n}`$, proving the proposition. $`\mathrm{}`$

Now let $`D_N(t)`$ be a solution of the N-particle von Neumann equation (7) and let $`F(t)`$ be a solution of the TDHF equation (9). For $`1nN`$ define the $`n^{th}`$ difference

$$E_{N,n}(t)=D_{N:n}(t)F_n^{}(t).$$
(17)
From the N-particle hierarchy equations (12) and (13) and Proposition 4.1, it follows that

$`i\mathrm{}{\displaystyle \frac{d}{dt}}E_{N,n}(t)`$ $`=`$ $`_n^{(N)}(E_{N,n}(t))+{\displaystyle \underset{j=1}{\overset{n}{}}}[V_{j,n+1},E_{N,n+1}(t)]_{:n}`$ (18)
$`+_n(t,N,D_N(0))_n(F(t))`$
for $`n=1,2,\mathrm{},N1`$. The characters $``$ and $``$ were chosen to evoke the words “error” and “remainder.” Indeed, in the next section we find bounds on these error terms under conditions on $`D_N(0)`$ and $`F(0)`$. The rest of this section is devoted to showing how such bounds lead to an upper bound on the differences $`E_{N,n}(t)`$.

To this end, let us define

$$\text{Err}(t,N,n)=_n(t,N,D_N(0))_n(F(t)).$$
(19)
Let $`U_{n,t}^{(N)}`$ denote the unitary operator $`\mathrm{exp}\left(\frac{it}{\mathrm{}}_{j=1}^nL_j^{(N)}\right)`$ on $`IH_n`$ and define isometries $`𝒰_{n,t}^{(N)}`$ on the trace class operators by

$$𝒰_{n,t}^{(N)}()=e^{\frac{it}{\mathrm{}}_n^{(N)}}()=U_{n,t}^{(N)}()U_{n,t}^{(N)}.$$

Then $`Z_{N,n}(t)=𝒰_{n,t}^{(N)}(E_{N,n}(t))`$ has the same trace norm as $`E_{N,n}(t)`$ and satisfies

$$\frac{d}{dt}Z_{N,n}(t)=\frac{i}{\mathrm{}}\underset{j=1}{\overset{n}{}}[V_{j,n+1},Z_{N,n+1}(t)]_{:n}\frac{i}{\mathrm{}}𝒰_{n,t}^{(N)}\text{Err}(t,N,n)$$
(20)
for $`n=1,2,\mathrm{},N1`$. From (20) it follows that

$`E_{N,n}(t)_{tr}`$ $`=`$ $`Z_{N,n}(t)_{tr}`$
$``$ $`Z_{N,n}(0)_{tr}+{\displaystyle \frac{2Vn}{\mathrm{}}}{\displaystyle _0^t}Z_{N,n+1}(s)_{tr}𝑑s`$
$`+`$ $`{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^t}𝒰_{n,t}^{(N)}(\text{Err}(s,N,n))_{tr}𝑑s`$
for $`n=1,2,\mathrm{},N1`$. Recalling that $`Z_{N,n+1}(s)_{tr}=E_{N,n+1}(s)_{tr}`$ and that $`𝒰_{n,t}^{(N)}(\text{Err}(s,N,n))_{tr}=\text{Err}(s,N,n)_{tr}`$, the preceding inequality becomes

$$E_{N,n}(t)_{tr}\epsilon (N,n,t)+\frac{2Vn}{\mathrm{}}_0^tE_{N,n+1}(s)_{tr}𝑑s$$
(21)
if we define

$$\epsilon (N,n,t)=E_{N,n}(0)_{tr}+\frac{1}{\mathrm{}}_0^t\text{Err}(s,N,n)_{tr}𝑑s.$$
(22)
Beginning from (21) and iterating the inequality $`m`$ times (for some $`mNn1`$) we obtain our desired bound on the trace norm of $`E_{N,n}(t)`$:

$`E_{N,n}(t)_{tr}`$ $``$ $`{\displaystyle \underset{k=0}{\overset{m}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n+k1}{n1}}\right)\left({\displaystyle \frac{2Vt}{\mathrm{}}}\right)^k\epsilon (N,n+k,t)`$
$`+`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{n+m1}{n1}}\right)\left({\displaystyle \frac{2Vt}{\mathrm{}}}\right)^m\underset{s[0,t]}{sup}\left\{E_{N,n+m+1}(s)_{tr}\right\}.`$
## 5 Error estimates

In this section we collect the error estimates that will be used to prove Theorem 3.1.

If $`D_N(0)`$ is a density operator then the solution $`D_N(t)`$ of the N-particle von Neumann equation (7) is a density operator for all $`t>0`$, and it is clear from (13) that

$$_n(t,N,D_N(0))_{tr}3\frac{n^2}{N}V$$
(24)
for all $`t`$.

###### Lemma 5.1

If $`\left\{D_N\right\}`$ has Slater closure then

$$\underset{N\mathrm{}}{lim}D_{N:1}=0.$$

Proof: The trace norm of $`D_{N:1}^2`$ equals the sum of the squares of the eigenvalues of $`D_{N:1}`$. Since the operator norm of $`D_{N:1}`$ equals its largest eigenvalue, it follows that $`D_{N:1}D_{N:1}^2_{tr}^{1/2}`$. But $`D_{N:1}^2=\left\{D_{N:1}^2U_{(12)}\right\}_{:1}`$, whence

$$D_{N:1}^2_{tr}=\left\{D_{N:2}D_{N:1}^2(IU_{(12)})\right\}_{:1}_{tr}D_{N:2}D_{N:1}^2\mathrm{\Sigma }_2_{tr}.$$

The Slater closure of $`\left\{D_N\right\}`$ implies that the right-hand side of the preceding inequality tends to $`0`$ as $`N\mathrm{}`$. $`\mathrm{}`$

###### Lemma 5.2

If $`F`$ is a density operator then $`F_n^{}_{tr}1`$ for all $`n`$.

Proof: Since $`\mathrm{\Sigma }_n=\left(\mathrm{\Sigma }_n\right)^{}=\frac{1}{n!}\left(\mathrm{\Sigma }_n\right)^2`$ commutes with $`F^n`$, it follows that $`F_n^{}=F^n\mathrm{\Sigma }_n=\frac{1}{n!}\mathrm{\Sigma }_n\left(F^n\right)\mathrm{\Sigma }_n`$ is a nonnegative operator. Thus, the trace norm of $`F_n^{}`$ equals its trace. This trace is

$$\underset{j_1,\mathrm{},j_nJ}{}<e_{j_1}\mathrm{}e_{j_n},F^n\mathrm{\Sigma }_n(e_{j_1}\mathrm{}e_{j_n})>$$

where $`\{e_j\}_{jJ}`$ is an orthonormal basis for $`IH`$ consisting of eigenvectors of $`F`$. This sum may be taken over distinct indices $`j_1,\mathrm{},j_nJ`$, since $`\mathrm{\Sigma }_n`$ annihilates all tensor products $`e_{j_1}\mathrm{}e_{j_n}`$ with repeating factors, so that

$`\text{Tr}\left(F_n^{}\right)`$ $`=`$ $`{\displaystyle \underset{\stackrel{distinct}{j_1,\mathrm{},j_nJ}}{}}<e_{j_1}\mathrm{}e_{j_n},F^n\mathrm{\Sigma }_n(e_{j_1}\mathrm{}e_{j_n})>`$
$`=`$ $`{\displaystyle \underset{\stackrel{distinct}{j_1,\mathrm{},j_nJ}}{}}e_{j_1}\mathrm{}e_{j_n},F^n(e_{j_1}\mathrm{}e_{j_n})`$
$``$ $`\text{Tr}\left(F^n\right)=1`$
as asserted. $`\mathrm{}`$

The next lemma provides a bound on the trace norm of the remainder term $`_n(F)`$ when $`F`$ is a density operator. The bound is proportional to the operator norm of $`F`$.

###### Lemma 5.3

Let $`_n`$ be as in (14) and let $`F`$ be a density operator. Then

$$_n(F)_{tr}2n(n1)VF.$$
(25)
Proof: From (14) we see that $`_n(F)`$ equals

$$\left\{\underset{\stackrel{j,k=1}{jk}}{\overset{n}{}}\left(V_{j,n+1}F^{n+1}U_{(k,n+1)}F^{n+1}U_{(k,n+1)}V_{j,n+1}\right)\left(\mathrm{\Sigma }_nI\right)\right\}_{:n}.$$

Since $`U_{(k,n+1)}`$ commutes with $`F^{n+1}`$ and since $`\mathrm{\Sigma }_nI`$ commutes with $`_{j,k:jk}U_{(k,n+1)}V_{j,n+1}`$, it follows that $`_n(F)`$ equals

$$\left\{\underset{\stackrel{j,k=1}{jk}}{\overset{n}{}}V_{j,n+1}U_{(k,n+1)}\left(F_n^{}F\right)\right\}_{:n}\left\{\left(F_n^{}F\right)\underset{\stackrel{j,k=1}{jk}}{\overset{n}{}}U_{(k,n+1)}V_{j,n+1}\right\}_{:n}.$$

Since the trace norm of a trace class operator equal the trace norm of its adjoint, it follows that

$`_n(F)_{tr}`$ $``$ $`2{\displaystyle \underset{\stackrel{j,k=1}{jk}}{\overset{n}{}}}\left\{V_{j,n+1}U_{(k,n+1)}\left(F_n^{}F\right)\right\}_{:n}_{tr}`$ (26)
$``$ $`2n(n1)\left\{V_{n1,n+1}U_{(n,n+1)}\left(F_n^{}F\right)\right\}_{:n}_{tr}.`$
But one may verify directly that

$$\left\{V_{n1,n+1}U_{(n,n+1)}\left(F_n^{}F\right)\right\}_{:n}=(I^{n1}F)V_{n1,n}F_n^{},$$
(27)
so that, by (3) and Lemma 5.2,

$$\left\{V_{n1,n+1}U_{(n,n+1)}\left(F_n^{}F\right)\right\}_{:n}_{tr}FVF_n^{}_{tr}FV.$$

Substituting this in (26) yields (5.3).

To verify (27), choose an orthonormal basis $`\{e_j\}_{jJ}`$ for $`IH`$ and check that the operators on both sides of (27) have the same matrix elements relative to the basis $`\{e_ie_j:i,jJ\}`$. $`\mathrm{}`$

Let $`F(t)`$ be a solution of the TDHF equation (9). Since the (operator) norm of $`F(t)`$ is constant, it follows from Lemma 5.3 that

$$_n(F(t))_{tr}2n(n1)VF(0)$$

for all $`t0`$. With this estimate for $`_n`$ and estimate (24) for $`_n`$, it follows that $`\text{Err}(t,N,n)`$ of equation (19) satisfies

$$\text{Err}(t,N,n)_{tr}2n^2V\left(\frac{2}{N}+F(0)\right)$$

and $`\epsilon (N,n,t)`$ of equation (22) satisfies

$`\epsilon (N,n,t)`$ $`=`$ $`2n^2V{\displaystyle \frac{t}{\mathrm{}}}\left({\displaystyle \frac{2}{N}}+F(0)\right)`$ (28)
$`+E_{N,n}(0)_{tr}.`$
## 6 Proof of the theorem

Equipped with the estimates of the preceding sections, we proceed to the proof of Theorem 3.1.

So, let us assume that $`D_N(0)`$ is a symmetric density for each $`N`$ and that the sequence $`\{D_N(0)\}`$ has Slater closure. Let $`D_N(t)`$ be the solution of (7) with initial value $`D_N(0)`$, and let $`F^{(N)}(t)`$ be the solution of the TDHF equation (9) whose initial value is $`F^{(N)}(0)=D_{N:1}(0)`$. Let $`\left\{F^{(N)}\right\}_n^{}(t)`$ denote $`\left\{F^{(N)}(t)\right\}^n\mathrm{\Sigma }_n`$ and let $`E_{N,n}(t)`$ denote the difference between $`D_{N:n}(t)`$ and $`\left\{F^{(N)}\right\}_n^{}(t)`$.

We have the upper bound (LABEL:Bardos\_inequality) for the trace norm of $`E_{N,n}(t)`$, into which we now substitute the estimates (28). In the same stroke, we will replace the binomial coefficients $`\left(\genfrac{}{}{0pt}{}{n+k1}{n1}\right)`$ with the larger quantities $`(n+k)^n/n!`$ and we will use the fact that $`\underset{s[0,t]}{sup}\left\{E_{N,n+m+1}(s)_{tr}\right\}2`$ by Lemma 5.2. Also, let us set $`T=2Vt/\mathrm{}.`$ We obtain

$`E_{N,n}(t)_{tr}`$ $``$ $`{\displaystyle \frac{1}{n!}}{\displaystyle \underset{k=0}{\overset{m}{}}}(n+k)^nE_{N,n+k}(0)_{tr}T^k`$ (29)
$`+`$ $`{\displaystyle \frac{1}{n!}}{\displaystyle \underset{k=0}{\overset{m}{}}}(n+k)^{n+2}\left({\displaystyle \frac{2}{N}}+F^{(N)}(0)\right)T^{k+1}`$
$`+`$ $`{\displaystyle \frac{2}{n!}}(n+m)^nT^m`$
for $`mNn1`$. Fix $`T`$ to be less than $`1`$, i.e., fix $`t<\mathrm{}/2V`$. For fixed $`n`$, consider the limit of the right-hand-side of (29) as $`N`$ and $`m`$ tend to infinity. The individual terms (fixed $`k`$) tend to $`0`$, for $`F^{(N)}(0)`$ tends to $`0`$ by Lemma 5.1 and $`E_{N,n+k}(0)_{tr}`$ tends to $`0`$ thanks to the hypothesis that $`\{D_N(0)\}`$ has Slater closure (recall that $`F^{(N)}(0)=D_{N:1}(0)`$). On the other hand, the series on the right-hand-side of (29) are dominated, uniformly with respect to $`m`$, by a series that converges absolutely for $`T<1`$. It follows that

$$\underset{N\mathrm{}}{lim}E_{N,n}(t)_{tr}=0$$
(30)
if $`t<\mathrm{}/2V`$. When $`n=1`$, this shows that $`\underset{N\mathrm{}}{lim}D_{N:1}(t)F^{(N)}(t)_{tr}=0`$ and consequently

$$\underset{N\mathrm{}}{lim}D_{N:1}^n(t)\mathrm{\Sigma }_n\left\{F^{(N)}\right\}_n^{}(t)_{tr}=0$$

for $`n>1`$ and $`t<\mathrm{}/2V`$. From (30) again it follows that, for any $`n`$ and any $`t<\mathrm{}/2V`$,

$$\underset{N\mathrm{}}{lim}D_{N:n}(t)D_{N:1}^n(t)\mathrm{\Sigma }_n_{tr}=0,$$

i.e., $`\{D_N(t)\}`$ has Slater closure. This proves the theorem up to $`t=\mathrm{}/2V`$.

Let $`\tau =\mathrm{}/3V`$; the previous argument shows that the theorem holds for $`t[0,\tau ]`$. At time $`\tau `$, it is no longer the case that $`D_{N:1}(\tau )=F^{(N)}(\tau )`$. However, $`E_{N,n+k}(\tau )_{tr}0`$ and $`F^{(N)}(\tau )0`$ as $`N`$ tends to infinity — to see this, use (10) and the fact that $`F^{(N)}(0)0`$. An argument nearly identical to the one above shows that the theorem holds for $`t[\tau ,2\tau ]`$. This argument may be repeated to establish the conclusion of the theorem on each interval of the form $`[k\tau ,(k+1)\tau ]`$ for each nonnegative integer $`k`$, and hence for all $`t>0`$. $`\mathrm{}`$

## 7 Appendix: TDHF equations for wavefunctions

The time-dependent Hartree-Fock equations were invented by Dirac . He wrote the TDHF equations as a system of coupled Schrödinger equations for $`N`$ time-dependent electronic orbitals, and thence obtained the TDHF equation for the density operator. In this article we have studied the TDHF in the latter form: as a von Neumann equation for the density operator. This appendix explains how to recast the wavefunction formulation of the TDHF equations into the language of density operators used in this paper.

The starting point in this discussion is the linear $`N`$-body Schrödinger

$$i\mathrm{}\frac{}{t}\mathrm{\Psi }_N=\frac{\mathrm{}^2}{2}\underset{k=1}{\overset{N}{}}\mathrm{\Delta }_{x_k}\mathrm{\Psi }_N+\frac{1}{N}\underset{1k<lN}{}V(x_kx_l)\mathrm{\Psi }_N$$
(31)
where $`\mathrm{\Psi }_N\mathrm{\Psi }_N(t,x_1,\mathrm{},x_N)`$ is the $`N`$-particle wavefunction. (Note that the interaction term has been multiplied by $`1/(N)`$. This scaling has been introduced so that $`N\mathrm{}`$ may yield a mean-field equation for the single-particle density, namely, the TDHF equation.) The dynamics defined by (31) is unitary on $`L^2((IR^3)^N)`$. Therefore,

$$\left|\mathrm{\Psi }_N(t,x_1,\mathrm{},x_N)\right|^2𝑑x_1\mathrm{}𝑑x_N=1$$

for all $`t0`$ if the same equality holds at $`t=0`$ (as is the case if $`|\mathrm{\Psi }_N|^2`$ is meant to be interpreted as the probability density of the system of $`N`$ particles in its configuration space). In the language of (1), $`L_k=\frac{\mathrm{}^2}{2}\mathrm{\Delta }_{x_k}`$ while $`V_{kl}`$ denotes the multiplication by $`V(x_kx_l)`$. The TDHF equations corresponding to (31) may be written as a system of $`N`$ coupled Schrödinger equations for orthonormal orbitals $`\psi _1(t,x),\psi _2(t,x),\mathrm{},\psi _N(t,x)`$:

$`i\mathrm{}{\displaystyle \frac{}{t}}\psi _k(t,x)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2}}\mathrm{\Delta }_x\psi _k(t,x)+\psi _k(t,x){\displaystyle \frac{1}{N}}{\displaystyle \underset{l=1}{\overset{N}{}}}{\displaystyle V(xz)|\psi _l(t,z)|^2𝑑z}`$ (32)
$`{\displaystyle \frac{1}{N}}{\displaystyle \underset{l=1}{\overset{N}{}}}\psi _l(t,x){\displaystyle V(xz)\psi _k(t,z)\overline{\psi _l(t,z)}𝑑z}`$
The $`N`$ orbitals remain orthonormal at all times; if $`\psi _1(t,x),\psi _2(t,x),\mathrm{},`$ $`\psi _N(t,x)`$ is a solution of (32) and

$`{\displaystyle \psi _k(0,x)\overline{\psi _l(0,x)}𝑑x}`$ $`=`$ $`\delta _{kl}\text{ then}`$
$`{\displaystyle \psi _k(t,x)\overline{\psi _l(t,x)}𝑑x}`$ $`=`$ $`\delta _{kl}\text{ for all }t0.`$
One way to obtain the TDHF equations (32) from the linear $`N`$-particle Schrödinger equation (31) is to solve a variational problem which would lead to (31) if unconstrained, but with the constraint that the $`N`$-particle wave equation remain a Slater determinant at all times . This constraint is imposed for the sake of obtaining an computationally amenable approximation to (31), and it is not justified on physical grounds. In effect, this paper proves that the constraint maintaining Slater determinants at all times is rigorously justified in the mean-field limit.

To see how the orbital form (32) of the TDHF equations relates to the TDHF equation (9) discussed in this paper, we shall first rewrite (9) as an equation for the integral kernel of a time-dependent density operator. To do this, we need to know how to translate the partial trace into the language of integral operators, for equation (9) involves a partial trace. Let $`T`$ be a trace class operator on $`L^2(IR^m\times IR^n)`$ having an integral kernel $`\rho (x,\xi ,y,\eta )`$ with $`x,yIR^m`$ and $`\xi ,\eta IR^n`$. The partial trace

$$T_{:m}\text{ is the operator with integral kernel}\rho (x,z,y,z)𝑑z.$$

We may now convert the TDHF equation (9) into an integro-differential equation for a time-dependent integral kernel: Let $`\rho \rho (t,x,y)`$ be the integral kernel of the operator $`F(t)`$ that appears in (9). Then $`F_2^{}(t)`$ has integral kernel

$$\rho (t,x_1,y_1)\rho (t,x_2,y_2)\rho (t,x_1,y_2)\rho (t,x_2,y_1),$$

while $`[V,F_2^{}(t)]_{:1}`$ has integral kernel

$$(V(x_1z)V(y_1z))(\rho (t,x_1,y_1)\rho (t,z,z)\rho (t,x_1,z)\rho (t,z,y_1))𝑑z,$$

and the TDHF equation (9) in the language of integral kernels is

$`i\mathrm{}{\displaystyle \frac{}{t}}\rho (t,x,y)={\displaystyle \frac{\mathrm{}^2}{2}}(\mathrm{\Delta }_x\mathrm{\Delta }_y)\rho (t,x,y)`$ (33)
$`+{\displaystyle (V(xz)V(yz))(\rho (t,x,y)\rho (t,z,z)\rho (t,x,z)\rho (t,z,y))𝑑z}.`$
It may be verified that a solution $`\psi _1(t,x),\psi _2(t,x),\mathrm{},\psi _N(t,x)`$ to the orbital form of the TDHF equations (32) yields a solution $`\rho (t,x,y)`$ to the integro-differential equation (33) via

$$\rho (t,x,y)=\frac{1}{N}\underset{k=1}{\overset{N}{}}\psi _k(t,x)\overline{\psi _k(t,y)}.$$

The rest of this appendix is meant to serve as a key for reading this paper with the Schrödinger wave equation (31) in mind.

To the wavefunction $`\mathrm{\Psi }_N(t,)`$ solution of the $`N`$-particle Schrödinger equation (31) one associates the operator $`D_N(t)`$ with integral kernel

$$\rho _N(t,x_1,\mathrm{},x_N,y_1,\mathrm{},y_N)=\mathrm{\Psi }_N(t,x_1,\mathrm{},x_N)\overline{\mathrm{\Psi }_N(t,y_1,\mathrm{},y_N)}.$$
(34)
The natural Hilbert space $`IH`$ in this context is $`IH=L^2(IR^3)`$, and $`IH^N`$ is isomorphic to $`L^2((IR^3)^N)`$ through the identification

$$\psi _1\mathrm{}\psi _N\underset{k=1}{\overset{N}{}}\psi _k(x_k).$$

The corresponding representation of the permutation group $`\mathrm{\Pi }_N`$ is given by the formula

$$(U_\pi \mathrm{\Psi })(x_1,\mathrm{},x_N)=\mathrm{\Psi }(x_{\pi ^1(1)},\mathrm{},x_{\pi ^1(N)})$$

for $`\pi \mathrm{\Pi }_N`$ and $`\mathrm{\Psi }_NL^2((IR^3)^N)`$. Hence the projection $`P_{𝒜_N}`$ is given by the formula

$$(P_{𝒜_N}\mathrm{\Psi }_N)(x_1,\mathrm{},x_N)=\frac{1}{N!}\underset{\pi \mathrm{\Pi }_N}{}\text{sgn}(\pi )\mathrm{\Psi }_N(x_{\pi (1)},\mathrm{},x_{\pi (N)}).$$

A wavefunction $`\mathrm{\Psi }`$ is antisymmetric if it is in the image of $`P_{𝒜_N}`$, or equivalently, if

$$\mathrm{\Psi }(x_{\pi (1)},\mathrm{},x_{\pi (N)})=\text{sgn}(\pi )\mathrm{\Psi }(x_1,\mathrm{},x_N)$$

for all $`\pi \mathrm{\Pi }_N`$. If $`\mathrm{\Psi }_N`$ is antisymmetric then the rank $`1`$ orthogonal projector $`P_{\mathrm{\Psi }_N}`$ with integral kernel as in (34) is fermionic in the sense of (2). Property (2) extends to all convex combinations of such projectors $`P_{\mathrm{\Psi }_N}`$, the fermionic density operators discussed in this article. Property (2) implies that fermionic density operators commute with all of the operators $`U_\pi `$, so that the integral kernel of a fermionic density operator is symmetric in the sense that

$$\rho _N(x_1,\mathrm{},x_N,y_1,\mathrm{},y_N)=\rho _N(x_{\pi (1)},\mathrm{},x_{\pi (N)},y_{\pi (1)},\mathrm{},y_{\pi (N)}))$$

for all $`\pi \mathrm{\Pi }_N`$.

In the case where $`\mathrm{\Psi }_N=\psi _1\mathrm{}\psi _N`$ with $`\psi _1,\mathrm{},\psi _N`$ orthonormal, one finds that the Slater determinant $`\sqrt{N!}P_{𝒜_N}(\mathrm{\Psi }_N)`$ truly is a determinant:

$$\sqrt{N!}(P_{𝒜_N}\mathrm{\Psi }_N)(x_1,\mathrm{},x_N)=\frac{1}{\sqrt{N!}}\left|\begin{array}{cccc}\psi _1(x_1)& \psi _1(x_2)& \mathrm{}& \psi _1(x_N)\\ \psi _2(x_1)& \psi _2(x_2)& \mathrm{}& \psi _2(x_N)\\ \mathrm{}& \mathrm{}& & \mathrm{}\\ \psi _N(x_1)& \psi _N(x_2)& \mathrm{}& \psi _N(x_N)\end{array}\right|.$$

Acknowledgement.

This research was supported by the French-Austrian “Amadeus” program (ÖAD V6) and by the Austrian START project “Nonlinear Schrödinger and quantum Boltzmann equations” of N.J.M. Also, F.G. was supported by the Institut Universitaire de France. The first three authors thank the ESI in Vienna for its hospitality. We also express our gratitude to Xavier Blanc, Eric Cancès and Claude Le Bris for numerous and helpful discussions on this paper.

Résumé

On montre dans ce travail que les équations d’évolution de Hartree-Fock décrivent la limite de l’équation de Schrödinger à $`N`$ corps pour $`N`$ tendant vers l’infini et une constante de couplage en $`O(1/N)`$ et pour des données initiales proches de déterminants de Slater. On ne considère ici que le cas de potentiels d’interaction binaires, symétriques et bornés. Lorsque $`N\mathrm{}`$, on montre que la suite des traces partielles “à un corps” de l’opérateur densité à $`N`$ corps converge, au sens des opérateurs à trace, vers la solution de l’équation de Hartree-Fock sous forme opératorielle.

Mots-clés: Equations de Hartree-Fock; Problème à $`N`$ corps quantique; Approximation de champ moyen


text:
# Hadro-Chemistry and Evolution of (Anti-) Baryon Densities at RHIC

\[

## Abstract

The consequences of hadro-chemical freezeout for the subsequent hadron gas evolution in central heavy-ion collisions at RHIC and LHC energies are discussed with special emphasis on effects due to antibaryons. Contrary to naive expectations, their individual conservation, as implied by experimental data, has significant impact on the chemical off-equilibrium composition of hadronic matter at collider energies. This may reflect on a variety of observables including source sizes and dilepton spectra.

\]

The systematic study of hadron production in central collisions of heavy nuclei at (ultra-) relativistic energies has revealed strong evidence for a high degree of equilibration achieved in the course of these reactions. In particular, it has been shown that, over a large range of collision energies (from SIS to SPS) , the measured final state hadron abundances can be well described by essentially two parameters, the temperature $`T`$ and baryon chemical potential $`\mu _B`$ (which bear little sensitivity on centrality ). Importantly, the inferred collision-energy dependence reinforces confidence in the notion that one is actually probing different regions of the QCD phase diagram in these experiments. The analyses are in line with the naive expectation that with increasing CM-energy $`\sqrt{s}`$ of the colliding nuclei, the central rapidity regions are increasingly governed by interactions between low-$`x`$ partons of the incoming nucleons, thus leading to a decrease in the net baryon content. This is reflected in a systematic decrease of the baryon chemical potential with $`\sqrt{s}`$, accompanied by an increase in temperature. In fact, for the top SPS energies ($`\sqrt{s}=17`$ AGeV) and for the first RHIC data at $`\sqrt{s}=130`$ AGeV, the such extracted $`(\mu ,T)`$ values are very close to the expected phase boundary to the Quark-Gluon Plasma (QGP). Under RHIC conditions this implies a copious production of baryon-antibaryon pairs. Whether those are solely of thermal origin, or what exactly the underlying microscopic (partonic) production mechanism is, remains a matter of debate at present. It is clear, however, that, as with any hadronic observable, the impact of the later hadronic stages has to be well understood before any firm conclusions can be drawn.

At SPS and RHIC/LHC energies the total hadron densities implied by the hadro-chemical analyses are not small, and one cannot expect the system to decouple at this point. It has been realized, however, that hadronic cross sections for elastic scattering, especially when involving intermediate resonances (e.g., $`\pi \pi \rho `$, $`\pi N\mathrm{\Delta }`$, $`\pi KK^{}`$), are typically much larger than those for inelastic (number-changing) reactions. The associated timescales for thermal equilibration are thus much smaller than for the chemical one. At SPS energies and above, thermal equilibrium in the expanding hadron gas can be maintained for about 10 fm/c. This is much below typical chemical relaxation times and has led to the picture of a chemical freezeout at $`(\mu _{ch},T_{ch})`$ (as inferred from hadron ratios) with $`T_{ch}`$ 160-180 MeV, followed by a thermal decoupling at significantly lower temperatures, $`T_{th}`$ 100-120 MeV. The existence of this intermediate hadronic phase is well supported as it importantly figures into both hadronic  and electromagnetic observables .

An immediate consequence of this picture is that all hadron species which are not subject to strong decays, have to be effectively conserved in number subsequent to chemical freezeout. In statistical mechanics language, this can be described by the build-up of (effective) chemical potentials, $`\mu _\pi `$, $`\mu _K`$, $`\mu _\eta `$, etc., in the hadronic evolution towards thermal freezeout. Corresponding thermodynamic trajectories at SPS energies have been constructed . In this context, antibaryons play a special role. Although their production from AGS to RHIC energies follows chemical-freezeout systematics, this cannot be expected a priori, since their annihilation cross section on baryons is large (e.g., $`\sigma _{p\overline{p}n\pi }50`$ mb with an average $`n`$=5-6 at the typical thermal energies). Under SPS conditions, this issue has been resolved in ref. , where it is shown that the backward reaction of multi-pion fusion, in connection with the rather large oversaturation of the pion phase space towards thermal freezeout (entailing $`\mu _\pi ^{th}70MeV`$), is capable of approximately maintaining the chemical-freezeout number of antibaryons. On the other hand, with a $`\overline{p}/p`$ ratio of 5-7%  and a final $`\pi /B`$ ratio (after strong decays) of around 5, the detailed treatment of antibaryons does not exert a large influence on the bulk thermodynamics of the hadronic evolution. The main point of this note is to show that this is no longer true at RHIC energies. Based on the experimental fact that also at RHIC the chemical freezeout abundances survive the hadron gas phase, we will argue that the explicit conservation of antibaryon number significantly affects the composition (and possibly evolution) of the hadronic expansion at collider energies. In particular, we will construct an explicit isentropic trajectory in the $`\mu _B`$-$`T`$ plane including effective number-conservation constraints.

Before we proceed with a more quantitative assessment, let us further motivate our objective by the following estimate: RHIC data at $`\sqrt{s}=130`$ AGeV for central Au-Au collisions  imply a midrapidity density of charged pions plus kaons of about $`dN_{ch}/dy600`$, so that $`dN_{\pi +K}/dy900`$. On the other hand, the total baryon number has been estimated to be $`dN_B/dy100`$  so that $`dN_{B+\overline{B}}/dy165`$ (using $`\overline{B}/B0.65`$ ). Since a (anti-) baryon carries about twice the entropy than that of pion, it follows that the former account for at least 25% of the total entropy at midrapidity.

Along the lines of standard hadro-chemical analyses we base the following calculations for the thermodynamic state variables on an ideal hadron gas***for a discussion of the influence of hadronic in-medium effects, see, e.g., refs. . including a set of resonance states comprised of 37 meson and 37 baryon species (up to masses $`m_M=1.7`$ GeV and $`m_B=2GeV`$, respectively). We will focus on conditions representative for midrapidities in central Au+Au collisions at $`\sqrt{s}=200`$ AGeV (very similar results emerge for a net-baryon free environment appropriate for LHC). For definiteness, we employ a specific entropy (i.e., entropy per net baryon) of $`S/N_B^{net}=250`$ which, together with a chemical freezeout temperature of $`T_{ch}=180`$ MeV, results in a baryon chemical potential of $`\mu _B^{ch}=24`$ MeV. This agrees well with (moderate) extrapolations of very similar thermal model analyses  based on RHIC data at $`\sqrt{s}=130`$ AGeV. By definition, all meson-chemical potentials are zero at $`T_{ch}`$, and we also neglect (small) isospin ($`\mu _I`$) and strangeness ($`\mu _s`$) chemical potentials. As usual, the effective pion number, which is subject to conservation, is defined by including all strongly decaying resonances with lifetimes shorter than the duration of the interacting hadronic phase, i.e.,

$$N_\pi ^{eff}=V_{FB}\underset{i}{}N_\pi ^{(i)}\varrho _i(T,\mu _i),$$
(1)
and likewise for $`K`$, $`\eta `$ and $`\eta ^{}`$. In eq. (1), $`V_{FB}`$ denotes the fireball 3-volume (which is related to the centrality of the collision but does not affect thermodynamic quantities) and $`\varrho _i`$ the number density of hadron $`i`$. The effective pion-number $`N_\pi ^{(i)}`$ of a given resonance is determined by the number of pions in its decay modes including the branching ratios, e.g., $`N_\pi ^{(\rho )}=2`$, $`N_\pi ^{(\mathrm{\Delta })}=N_\pi ^{(K^{})}=1`$, $`N_\pi ^{(N(1520))}=0.551+0.452`$. The same weighting applies to the corresponding effective chemical potential of each resonance, e.g., $`\mu _\rho =2\mu _\pi `$, $`\mu _\mathrm{\Delta }=\mu _N+\mu _\pi `$, $`\mu _K^{}=\mu _\pi +\mu _K`$, $`\mu _{N(1520)}=\mu _N+1.45\mu _\pi `$, which means that the strong decays and formations are in relative chemical equilibrium. For a few states this distinction is not very sharp, most notably $`\omega `$ and $`\varphi `$ mesons. Their free lifetimes are well above hadronic fireball lifetimes, but in-medium effects are expected to alter this situation . In the former case, some of the $`\omega `$’s and $`\varphi `$’s decay after chemical freezeout without being regenerated, i.e., their number decreases; in the latter case, their abundances stay (at least for a while) in chemical equilibrium which, in the case of large pion- and kaon-chemical potentials, could even increase their numbers. In our calculation we employ an intermediate scenario keeping the individual $`\omega `$ and $`\varphi `$ numbers approximately constant (which is numerically close to the chemical equilibrium assumption). An important point now is the treatment of antibaryons as first pointed out in ref. . We also conserve their number by introducing another effective chemical potential $`\mu _{\overline{B}}^{eff}`$; this implies, e.g., $`\mu _{\overline{N}}=\mu _N+\mu _{\overline{B}}^{eff}`$, $`\mu _{\overline{\mathrm{\Delta }}}=\mu _N+\mu _{\overline{B}}^{eff}+\mu _\pi `$.

A thermodynamic trajectory in the $`\mu _N`$-$`T`$ plane is now constructed starting from chemical freezeout under the additional assumption of total entropy conservation (as applied in ideal fluid dynamics). Together with net baryon number conservation (which is, of course, exact) this amounts to keeping the specific entropy per net baryon fixed, $`S/N_B^{net}=s/\varrho _B^{net}=250`$ for RHIC with

$`s`$ $`=`$ $`{\displaystyle \underset{i}{}}d_i{\displaystyle }{\displaystyle \frac{d^3k}{(2\pi )^3}}[\pm f(\omega _i)\mathrm{ln}f(\omega _i)`$ (3)
$`+(1f(\omega _i))\mathrm{ln}(1f(\omega _i)],`$
$`\varrho _B^{net}`$ $`=`$ $`{\displaystyle \underset{i}{}}d_{B_i}{\displaystyle \frac{d^3k}{(2\pi )^3}\left[f^{B_i}(\mu _{B_i},T)f^{\overline{B}_i}(\mu _{\overline{B}_i},T)\right]}`$ (4)
denoting entropy- and net baryon number-density, respectively (in eq. (3), upper (lower) signs refer to fermions (bosons), and $`f`$ are the pertinent thermal distribution functions). The resulting trajectory is shown by the full line in Fig. 1. For comparison, we also display a trajectory without antibaryon-number conservation (dashed line). The difference is rather dramatic, indicating a strongly increased baryon chemical potential towards the regions of expected thermal freezeout at around $`T_{th}120`$ MeV. This is induced by an approximately linear increase of the antinucleon chemical potential, which is displayed in Fig. 2 along with nucleon, pion, kaon, and eta chemical potentials. At chemical freezeout the sum of baryon and antibaryon density is close to nuclear saturation density $`\rho _0=0.16`$ fm-3, $`\varrho _{B+\overline{B}}1.1\varrho _0`$. On the other hand, at $`T=120`$ MeV, $`\varrho _{B+\overline{B}}0.3\varrho _0`$ for the solid-line trajectory, but $`\varrho _{B+\overline{B}}0.02\varrho _0`$ for the dashed-line trajectory.

One can readily verify that the solid-line trajectory in Fig. 1 is consistent with the measured $`\overline{p}/p`$ ratio, which is approximately given by $`\mathrm{exp}[\mathrm{\Delta }\mu _N/T]`$ with $`\mathrm{\Delta }\mu _N=\mu _N\mu _{\overline{N}}`$. It is also important to note that in the $`\mu _{\overline{B}}^{eff}0`$ case, the pion-chemical potentials remain close to zero. Thus, $`\overline{B}`$ conservation at collider energies is at the origin of large meson chemical potentials in the late hadronic stagesThis is opposite to SPS conditions, where, following the arguments of ref. , pion-oversaturation is responsible for regenerating antibaryons.. The reason can be traced back to the substantial amount of entropy stored in $`B\overline{B}`$ excitations. At a given temperature, the smaller available entropy per pion at fixed pion number can only be realized by a finite $`\mu _\pi `$, as also noted in ref. . To illustrate this fact, we display in Fig. 3 the entropy per pion in a $`\pi `$-$`\rho `$ gas. One finds that for $`T=120`$ MeV, $`S/N_\pi `$ is reduced by about 20% when raising $`\mu _\pi `$ from zero to its actual value of $``$80 MeV corresponding to the solid-line RHIC-trajectory of Fig. 2. This is in line with the rough estimate of the entropy carried by $`B\overline{B}`$ pairs at RHIC above.

It furthermore follows that at given temperature the fireball volume $`V_{FB}`$ of the system with finite $`\mu _{\overline{B}}^{eff}`$ is substantially reduced, i.e., the expanding hadron gas cools faster. E.g., at $`T=120`$ MeV the ratio of the 3-D fireball volumes for the $`\mu _{\overline{B}}^{eff}0`$ and $`\mu _{\overline{B}}^{eff}>0`$ solution for the trajectories amounts to a factor of $``$3. Whether this will translate into smaller source size radii as extracted via HBT correlation measurement (and thus contribute to an explanation of the current ’HBT-puzzle’ ) remains to be investigated within a (hydro-) dynamical simulation. Other observable consequences of the increased densities in the hadronic phases at RHIC concern low-mass dilepton emission, where in-medium effects are rather sensitive to (the sum of) anti-/baryon densities . Accordingly improved estimates have been presented in ref. .

In our analysis we did not allude to underlying mechanisms that could provide a constant antibaryon abundance in the hadronic phase at RHIC. As mentioned above, at SPS energies this issue could be resolved  within a thermal framework by accounting for multi-pion fusion reactions in the presence of a large $`\mu _\pi `$. Although it is conceivable that this mechanism is also operative at RHIC, the situation is more involved due to the feedback of the $`\overline{B}`$ abundance on $`\mu _\pi `$. This problem is not easily addressed within transport or cascade simulations either due to the difficulties of treating the backward reaction of multi-pion fusion (see ref.  for recent progress). Corresponding trajectories that have been extracted from, e.g., UrQMD transport calculations  are therefore resembling the dashed lines in Fig. 1.

In summary, based on the notion of hadro-chemical freezeout, we have constructed (isentropic) thermodynamic trajectories for the hadronic phase in ultra-relativistic heavy-ion collisions. Consistency with the observed particle ratios has been enforced by effective conservation laws, implemented via finite chemical potentials for stable hadrons (w.r.t. strong interactions). The main emphasis was on consequences of conserved antibaryon-number. At SPS energies and below, these are small, but at collider energies (RHIC and LHC) we have found appreciable modifications in the composition of the expanding hadronic system. In particular, the entropy carried by $`B\overline{B}`$ pairs is not small and has been identified as the origin of large meson-chemical potentials, which in turn leads to smaller and thus denser systems at temperatures prior to thermal freezeout. This could have important consequences for hadronic and electromagnetic observables, such as source radii and medium effects in low-mass dilepton production. The former are rather sensitive to the detailed expansion dynamics of a heavy-ion reaction, which we did not address here. For the same reason, the question where on the given trajectories thermal freezeout in the hadronic phase occurs has not been answered, but needs to be investigated in hydrodynamical and/or microscopic transport simulations.

This work was supported by the U.S. Department of Energy under Grant No. DE-FG02-88ER40388.