text:
# Thermodynamics and Crossover Phenomena in the Correlation Lengths of the One-Dimensional t-J Model

## I Introduction

The t-J model is one of the most fundamental systems of strongly correlated electrons. The two-dimensional version has attracted much attention because it is believed that it describes the basic interactions in the copper-oxygen planes of high-$`T_c`$ superconductors. For the one-dimensional (1D) t-J model much progress has been achieved using various analytical and numerical techniques.Ogata et al. (1991); Hellberg and Mele (1991); Bares and Blatter (1990); Kawakami and Yang (1990); Jüttner et al. (1996); Chen and Lee (1996); Kobayashi et al. (1996); Nakamura et al. (1997) At the supersymmetric point $`J/t=2`$ the model is solvable by the Bethe ansatz and ground state properties as well as the excitation spectra have been obtained exactly.Bares and Blatter (1990) Because the two critical excitations of spin and charge type are separated, the properties can be described by two independent $`c=1`$ Virasoro algebras. By a combination of finite-size results from the Bethe ansatz and conformal field theory (CFT) it is therefore also possible to calculate the critical exponents of algebraically decaying correlation functions.Kawakami and Yang (1990) This explicitly shows that the t-J model at the supersymmetric point behaves as a Tomonaga-Luttinger liquid (TLL) for all electron densities. Thermodynamic quantities at this special point have been obtained by thermodynamic Bethe ansatz Okiji and Suga (1997) as well as by a combination of a Trotter-Suzuki mapping leading to a quantum transfer matrix (QTM) and the Bethe ansatz.Jüttner et al. (1996) Exact results are also available in the limit $`J/t0`$ where the t-J model is equivalent to the Hubbard model with $`U/t\mathrm{}`$, showing again TLL behavior.Frahm and Korepin (1990) It is therefore believed that the t-J model shows TLL properties for all $`0J/t2`$, what is supported by various numerical calculations.Ogata et al. (1991); Hellberg and Mele (1991); Chen and Lee (1996); Kobayashi et al. (1996) In these numerical works there is also general agreement that the t-J model phase separates for $`J/t=2.8`$ to $`3.5`$ depending on the electron density.

Already Ogata et al. Ogata et al. (1991) conjectured a third phase with a spin gap in the low-density region for $`J/t>2`$. However, by calculating the spin susceptibility for an electron density $`n=1/3`$ on small chains no evidence for a spin gap was found. Also the variational quantum Monte Carlo (QMC) calculations of Hellberg and Mele Hellberg and Mele (1991) could not confirm the appearance of a spin gap. By using the same method but other trial wave functions a phase with Luther-Emery (LE) properties was found by Chen and Lee Chen and Lee (1996) and Kobayashi et al. Kobayashi et al. (1996) at low densities and $`2<J/t<3.1`$. However, the obtained result strongly depends on the trial wave functions used in the calculations. Completely different phase boundaries with the spin-gap phase extending into the high density region have been obtained by Nakamura et al. Nakamura et al. (1997) using a renormalization group treatment of the Tomonaga-Luttinger model under the assumption that the spin gap is caused by an attractive backward scattering process. They argue that the spin gap was underestimated in the numerical calculations, because it is the result of a marginal operator leading to an exponentially small gap.

In this paper we want to examine thermodynamics of the one-dimensional t-J model in the whole $`J/t`$-parameter region. In particular, we are interested in the temperature dependence of correlation lengths. We use a grandcanonical description of the model, where the Hamiltonian is given by

$`H`$ $`=`$ $`t{\displaystyle \underset{i,\sigma }{}}P(c_{i,\sigma }^{}c_{i+1,\sigma }+c_{i+1,\sigma }^{}c_{i,\sigma })P`$
$`+`$ $`J{\displaystyle \underset{i}{}}\left(𝑺_i𝑺_{i+1}{\displaystyle \frac{n_in_{i+1}}{4}}\right)h{\displaystyle \underset{i}{}}S_i^z\mu {\displaystyle \underset{i}{}}n_i`$
with a magnetic field $`h`$, a chemical potential $`\mu `$ and $`P`$ being the projection operator onto the Hilbert-subspace without double occupancy. To study thermodynamic properties also away from the supersymmetric point the TMRG provides a powerful numerical tool. This method is particularly suited, because the thermodynamic limit is performed exactly and it has been applied successfully to various one-dimensional systems before.Bursill et al. (1996); Wang and Xiang (1997); Shibata (1997); Klümper et al. (1999) Compared with QMC methods it has further the advantage of never suffering under the fermion sign problem restricting QMC often to relatively high temperatures.

In Sec. II we give a brief introduction into the calculation of thermodynamic quantities by the TMRG algorithm and point to a novel algorithm with several advantages. Our results for the supersymmetric point are shown in Sec. III where we also compare the numerics with the Bethe ansatz results by Jüttner et al. Jüttner et al. (1996) In particular we calculate for this special point the temperature dependence of various correlation lengths and static correlation functions and discuss the low-temperature behavior in comparison to CFT predictions by Kawakami and Yang Kawakami and Yang (1990) in detail. Because the t-J model is realized by the Hubbard model in the limit $`Ut`$ the case $`J=2t^2/|U|<t`$ is physically very relevant. Results for $`J/t=0.35`$, a value often used in literature, are given in Sec. IV. As already mentioned the t-J model phase separates for $`J/t`$ large. The meaning of this in the grand-canonical ensemble is explained in Sec. V. Without using any assumption about the ground state or the low energy effective theory, the existence of a LE phase is proven in Sec. VI by calculating directly spin susceptibilities as well as spin-spin and density-density correlation lengths. In Sec. VII we study the t-J model with an additional Ising-like anisotropy. The final section is devoted to our conclusions and in appendix A we compare numerical and exact results for the density-density correlation function of free spinless fermions.

## II TMRG

After a decomposition of the Hamiltonian $`H`$ into even ($`H_e`$) and odd parts ($`H_o`$) the partition function is expressed by means of the Trotter formula

$$Z=\text{Tr}\text{e}^{\beta H}=\underset{M\mathrm{}}{lim}\text{Tr}\left\{\left[\text{e}^{ϵH_e}\text{e}^{ϵH_o}\right]^M\right\}$$
(2)
with $`ϵ=\beta /M`$ and $`\beta `$ being the inverse temperature leading to a classical model on a lattice with checkerboard structure. The column-to-column transfer matrix (QTM),$`T_M`$, is a non-symmetric matrix describing the evolution along the spatial direction. The thermodynamic limit with fixed Trotter number M is performed exactly, because the free energy of the infinite chain is given solely by the largest eigenvalue $`\mathrm{\Lambda }_0`$ of this QTM

$$f_{\mathrm{},M}=\frac{T}{2}\mathrm{ln}\mathrm{\Lambda }_0,$$
(3)
where $`\mathrm{\Lambda }_0`$ is unique and a real, positive number for all temperatures. Since a vanishing gap between leading and next-leading eigenvalue indicates a phase transition, such a degeneracy is not possible for a 1D quantum system at finite temperature. The calculations are simplified by the conservation laws for spin and particle number or equivalently for the number of particles with spin up ($`N_{}`$) and spin down ($`N_{}`$) leading to a block structure of $`T_M`$. Expectation values of local operators can be expressed in terms of the left and right eigenvectors belonging to the largest eigenvalue and calculated directly. This is used for the inner energy $`u`$ and the magnetization $`m`$. The behavior of a two-point correlation function is given by

$$A_0A_r=\underset{\alpha }{}M_\alpha \text{e}^{r/\xi _\alpha }\text{e}^{\text{i}k_\alpha r}$$
(4)
with $`M_\alpha `$ being matrixelements. The correlation lengths $`\xi _\alpha `$ can be calculated from the QTM by

$$\xi _\alpha ^1=\frac{1}{2}\mathrm{ln}\left|\frac{\mathrm{\Lambda }_0}{\mathrm{\Lambda }_\alpha }\right|$$
(5)
and the wavevectors $`k_\alpha `$ by

$$k_\alpha =\frac{1}{2}\mathrm{arg}\left(\frac{\mathrm{\Lambda }_\alpha }{\mathrm{\Lambda }_0}\right)+n\pi (n=0\text{or}n=1),$$
(6)
where $`\mathrm{\Lambda }_\alpha `$ are eigenvalues in the block of the QTM with non-vanishing matrixelements $`M_\alpha `$. This condition is controlled by a direct evaluation of

$$M_\alpha =\frac{\mathrm{\Psi }_0^L|T_M(A_0)|\mathrm{\Psi }_\alpha ^R\mathrm{\Psi }_\alpha ^L|T_M(A_r)|\mathrm{\Psi }_0^R}{\mathrm{\Lambda }_0\mathrm{\Lambda }_\alpha },$$
(7)
where $`\mathrm{\Psi }_0^L|`$, $`|\mathrm{\Psi }_0^R`$ are the left and right eigenvectors belonging to the largest eigenvalue $`\mathrm{\Lambda }_0`$ and $`\mathrm{\Psi }_\alpha ^L|`$, $`|\mathrm{\Psi }_\alpha ^R`$ the left and right eigenvectors belonging to $`\mathrm{\Lambda }_\alpha `$. The value for the wavevector $`k_\alpha `$ is not unique, because in the checkerboard decomposition of the partition function a local transfer matrix covers 2 sites. To overcome this ambiguity we have also applied a different Trotter-Suzuki mapping leading to a classical lattice with alternating rows. Because the transfer matrix in this formalism can be formulated for a single column, the wavevector in Eq. (6) is determined unambiguously. Details of this modified TMRG algorithm will be published elsewhere.Sirker and Klümper (2002) We have convinced ourselves that the accuracy of both algorithms is of the same order and use both of them in the following. During all calculations we use the infinite DMRG algorithm to increase the Trotter number $`M`$ being equivalent to a decrease of the temperature $`T`$. The number of states kept in the DMRG varies between 64 and 240 and the fixed parameter $`ϵ`$ between $`0.025`$ and $`0.05`$.

## III The supersymmetric point

To check the numerics we have calculated several thermodynamic quantities for the supersymmetric point $`J/t=2`$ and compared with the exact results by Bethe ansatz.Jüttner et al. (1996) We have chosen three chemical potentials corresponding in the low-temperature limit to a high, medium and low electron density (see Fig. 1).

Obviously the particle density $`n`$ at a given chemical potential depends on temperature. Counting the degrees of freedom per lattice site immediately implies $`n=2/3`$ for any finite $`\mu `$ in the limit $`T\mathrm{}`$, whereas the density at finite temperature is given by $`n=n(\mu ,T)`$.

For the same chemical potentials with $`N=100`$ states retained in the DMRG algorithm and $`ϵ=0.05`$ the free energies and their accuracy are shown in Fig. 2 in comparison to the Bethe ansatz results by Jüttner et al. Jüttner et al. (1996) In all three cases the accuracy is of the order $`10^4`$ if $`T/t>0.1`$. For lower temperatures $`0.01T/t0.1`$ the errors for the medium and high electron density grow up to $`10^3`$, whereas the error in the low density case remains of the same order. More accurate results can be obtained if more states are retained or an extrapolation in the number of states $`N`$ and the small parameter $`ϵ`$ is made. However, the obtained accuracy is sufficient for our purposes. From CFT it is known that the low-temperature asymptotics is given by

$$f=e_0\frac{\pi }{6}\left(\frac{c_s}{v_s}+\frac{c_c}{v_c}\right)T^2$$
(8)
where $`v_{s,c}`$ are the velocities of the spinon and holon excitations, respectively. Here the central charges $`c_{s,c}`$ are equal to 1. According to Eq. (8) we have fitted the numerical data and tried to determine the errors of the fit parameters by a variation of the fit region from $`T/t[0.01:0.02]`$ to $`T/t[0.01:0.05]`$. The estimates for the ground-state energies $`e_0`$ coincide with the exact results with deviations that are slightly larger than the errors calculated by the fit procedure (see table 1). The remaining deviation is due to systematic errors. The parameter $`1/v_s+1/v_c`$ is difficult to obtain from such a fit, because the extent of the low-temperature region in which Eq. (8) is valid is very small (see Bethe ansatz results in Fig. 3) making it necessary to restrict the fit to the temperature region defined above. In such a small interval, however, a change of the parameter $`1/v_s+1/v_c`$ by a factor of 2 corresponds to deviations in the fitted free energy of the order $`10^3`$ only, so that the TMRG data are not accurate enough to determine this parameter. However, it is possible to calculate the velocities from TMRG by using the results for the susceptibilities and correlation lengths to be discussed later on. We do not want to pursue this further, because the estimation of $`v_s`$ and $`v_c`$ is not our main goal.

To calculate the specific heat $`c_n`$, we first have to calculate the entropy $`S`$. This can be done by using the relation $`S=(uf)/T`$, where the inner energy $`u`$ is directly calculated in the TMRG algorithm as a local expectation value or alternatively, by the numerical derivative $`S=f/T`$. At first sight it seems better to calculate the entropy without using numerical derivatives, however, a local expectation value directly involves the eigenvectors of the QTM which are generally less accurate than the corresponding eigenvalues. By comparing both results with the exact ones we have convinced ourselves that the numerical derivative leads to more accurate results. Because in the calculations only the chemical potential $`\mu `$ but not the particle density $`n`$ can be fixed, we calculate the specific heat by using the thermodynamic relation

$$c_n=T\left(\frac{S}{T}\right)_n=T\left[\left(\frac{S}{T}\right)_\mu \left(\frac{n}{T}\right)_\mu ^2\left(\frac{n}{\mu }\right)_T^1\right].$$
(9)
The second term is again evaluated using numerical derivatives.

This procedure enhances the numerical errors, but we are able to reproduce the exact results within errors of the order $`10^3`$ down to temperatures of $`T/t0.2`$ as shown in Fig. 3. The CFT result for the free energy, Eq. (8), leads to a linear temperature dependence of the specific heat at low temperatures with a coefficient given by $`\pi (1/v_c+1/v_s)/3`$. The Bethe ansatz results clearly show this linear behavior at very low temperatures but within the TMRG this temperature region is not accessible because of the errors caused by the numerical derivatives used in Eq. (9).

To calculate the spin susceptibility $`\chi _s`$, a small magnetic field $`h=10^2`$ is applied. From the resulting magnetization $`m`$ the susceptibility at vanishing magnetic field is evaluated by

$$\chi _s|_{h=0}=\frac{m}{h}.$$
(10)
Similarly the charge susceptibility $`\chi _c`$ or compressibility is given by

$$\chi _c=\frac{n}{\mu },$$
(11)
where again numerical derivatives are used and the variation of $`\mu `$ is typically of the order $`10^2`$.

Here we are able to reproduce the exact results down to temperatures of $`T/t=0.01`$ as shown in Fig. 4. The high temperature asymptotics of the spin susceptibility is given by the pure paramagnetic part $`\chi _s(T\mathrm{})2s(s+1)/9T`$ what is easily understood when $`\chi _s`$ is expressed as a sum of two-point spin correlations. For $`T=0`$ the spin susceptibility has a finite value as expected from the linear dispersion of the spinon excitations. Within CFT the zero temperature value is given explicitly by $`\chi _s=1/(2\pi v_s)`$ in good agreement with the numerical results. The maxima at finite temperatures are determined by the band structure. The charge susceptibility is given in the same way as sum over two-point density-density correlations leading to the high-temperature asymptotics $`\chi _c2/9T`$. For $`T=0`$ the charge susceptibility is again given by CFT as $`\chi _c=\xi _c^2(Q)/(\pi v_c)`$ with $`\xi _c(Q)`$ being the dressed charge. Note, that the charge susceptibility is diverging for $`T0`$ in the two limiting cases $`n0`$ and $`n1`$, because $`v_c0`$ in both cases. The spin susceptibility $`\chi _s`$ shows different behavior with divergence for $`T0`$ only in the limiting case $`n0`$ as only here $`v_s0`$.

To summarize our numerical findings for $`c`$, $`\chi _s`$ and $`\chi _c`$ in a qualitative manner we may use a picture of the ground state as a “liquid” consisting of bound singlet pairs of electrons. Above this ground state there are two different elementary excitations. First, a momentum transfer onto any individual pair is possible by keeping the bound pair intact (holon excitation). Second, a breaking of any bound pair into its components (electrons) with the two spin-1/2 objects coupling either to $`S=0`$ or 1 (spinon excitation). Once the “free” electrons are produced by spinon excitations, they may acquire momentum and energy individually which we refer to as incoherent single particle motion. Basically, any structure observed in the $`T`$-dependence of $`c`$, $`\chi _s`$ and $`\chi _c`$ can be attributed to the saturation of a particular type of excitation. In $`\chi _s`$ ($`\chi _c`$) we see a finite temperature maximum at a temperature of the order of $`v_s`$ ($`v_c`$). Note that the characteristic holon temperature is always lower than or equal to the spinon temperature. This is rather natural as the saturation of spinons is marked by the practical absence of pairs such that saturation of holon excitations must have occurred at lower or equal temperature. In principle, these structures are also visible in the specific heat. In most cases, however, only the spin-maximum is clearly noticeable. For generic densities the small holon-maximum disappears in the low temperature regime of $`c(T)`$ with rather steep slope due to the spinon-excitations. At lower densities the holon and spinon structures are located at similar temperatures. In all cases, the higher temperature region is dominated by a broad maximum due to incoherent single particle motion. The corresponding energies are of the order of spinon excitations at particle densities close to 1 in which case the structures merge.

### III.1 Correlation lengths

For the correlation lengths of the supersymmetric t-J model no exact results are available. To obtain nevertheless a measure for the accuracy of the numerical data we have calculated the leading correlation length for a system of free spinless fermions (see appendix A). The comparison with the exact result shows that TMRG yields reliable results down to temperatures $`T/t0.1`$. Because the free fermion model is harder to tackle within the DMRG scheme due to the rapidly decaying spectrum of the reduced density matrix, we believe that for the t-J model a similar or even better accuracy is obtained.

We concentrate now on the calculation of correlation lengths for the supersymmetric t-J model to study the crossover from the high T lattice into the quantum critical regime determined by $`Tt`$. In the high T lattice region we expect non-universal properties dependent on the microscopic Hamiltonian $`H`$, whereas the quantum critical regime should show universal TLL properties. Respecting the selection rules, the density-density (d-d CL’s) and longitudinal spin-spin correlation lengths (s-s CL’s) are in the block of the QTM with unchanged quantum numbers ($`\mathrm{\Delta }N_{}=\mathrm{\Delta }N_{}=0`$). To distinguish between them, the matrixelement $`M_\alpha `$ in Eq. (7) has to be calculated explicitly. With zero magnetic field and isotropic spin interactions, it turns out that $`M_\alpha `$ is either non-zero for the density or for the $`S^z`$ operator so that all d-d and s-s CL’s are different from each other. This is a consequence of the $`SU(2)`$ spin symmetry. If a magnetic field is applied this symmetry is broken and less stringent selection rules apply such that all eigenvalues discussed above contribute to the s-s as well as to the d-d correlation function. On the other hand the singlet pair operator $`P_i^s=c_{,i+1}c_{,i}`$ changes the quantum numbers by $`\mathrm{\Delta }N_{}=\mathrm{\Delta }N_{}=\pm 1`$ and the eigenvalues corresponding to the singlet pair CL’s are found in the block with these quantum numbers. The leading d-d and s-s CL’s (times temperature) and the corresponding wavevectors for $`J/t=2.0`$, $`\mu =1.4`$ and zero magnetic field are shown in Fig. 5.

At high temperatures the largest s-s CL is given by a real negative eigenvalue leading to $`\pi `$-oscillations. However, at a well defined crossover temperature $`T_c0.8`$ a real, positive eigenvalue becomes largest, which is associated with $`k=0`$. Regarding only the largest spin CL this means that there is a non-analyticity at $`T_c`$. We want to point out that any thermodynamic quantity derived from the free energy is an analytic function at finite $`T`$. Phase transitions and corresponding singularities only occur at $`T=0`$. However, quantities describing the asymptotics of correlation functions (i.e. CL) may show non-analyticities even at finite temperature. Such crossovers between CL’s are characteristic for the non-universal high T lattice regime, whereas no crossovers are expected in the universal quantum critical regime described by CFT. The second largest spin CL is given at low temperatures by a pair of complex conjugated eigenvalues. Complex eigenvalues always appear in pairs, because the QTM is nonsymmetric but real and lead to incommensurate oscillations, i.e. $`k0,\pi `$. An incommensurate wavevector $`k`$ is depending on temperature as seen in the inset of Fig. 5. The leading d-d CL shows no oscillations at high temperatures but is crossed at lower temperatures by a real, negative eigenvalue leading to $`\pi `$-oscillations before again the non-oscillating part dominates. Below $`T_c0.5`$ the largest d-d CL is given by a complex eigenvalue pair. The wavevector for this incommensurate part is again shown in the inset of Fig. 5.
It is illuminating to regard first the d-d correlation for a system of free spinless fermions. For such a system the CL’s can be easily calculated (see appendix A) and it turns out that a non-oscillating part and a part with $`k=\pi `$ arise, both with a CL given by $`\xi ^1=2\text{arcsinh}(\pi T)`$. Therefore the CL’s at high temperature behave as $`\xi ^1\mathrm{ln}(T)`$ consistent with the numerical results. In the low temperature limit it follows that $`\xi ^11/(2\pi T)`$ what is expected from CFT when the velocity $`v`$ and the scaling dimension $`x`$ are equal to $`1`$. The matrixelement for the CL is given by $`M=2/(\pi ^2+\beta ^2)`$ indicating that $`M0`$ when $`T0`$. This is reasonable, because at zero temperature the CL’s diverge and the matrixelements have to approach zero so that the series in Eq. (4) is not diverging but sums up to an algebraic function. The numerical results also show that the matrixelements $`M_\alpha 0`$ if $`T0`$.

To understand the behavior of the CL’s at low temperatures in detail it is useful to regard the results from the Bethe ansatz and the finite-size scaling technique in CFT at zero temperature.Kawakami and Yang (1990) In general, a two-point correlation function of the scaling fields $`\varphi _{\mathrm{\Delta }^\pm }(r,\tau )`$ with conformal weights $`\mathrm{\Delta }^\pm `$ is given by

$`\varphi _{\mathrm{\Delta }^\pm }(r,\tau )\varphi _{\mathrm{\Delta }^\pm }(0,0)`$
$`=`$ $`{\displaystyle \frac{\text{e}^{\text{i}(2\pi 2k_F2k_F)D_cx}\text{e}^{\text{i}(2\pi 2k_F)D_sx}}{(r\text{i}v_c\tau )^{2\mathrm{\Delta }_c^+}(r+\text{i}v_c\tau )^{2\mathrm{\Delta }_c^{}}(r\text{i}v_s\tau )^{2\mathrm{\Delta }_s^+}(r+\text{i}v_s\tau )^{2\mathrm{\Delta }_s^{}}}}`$
where the conformal weights are determined by quantum numbers of the elementary excitations. Considering here only the case of vanishing magnetic field and using the notation of Ref. Kawakami and Yang, 1990, the conformal weights reduce to a simple form
$$\mathrm{\Delta }_c^\pm (𝑰,𝑫)=\frac{1}{2}\left(\frac{I_c}{2\xi _c(Q)}\pm \xi _c(Q)\left(D_c+\frac{D_s}{2}\right)\right)^2+N_c^\pm $$
(12b)
$$\mathrm{\Delta }_s^\pm (𝑰,𝑫)=\frac{1}{4}\left(I_s\frac{I_c}{2}D_s\right)^2+N_s^\pm .$$
(12c)
$`\xi _c(Q)`$ is the dressed charge which can be calculated explicitly by solving an integral equation and $`(I_c,I_s,D_c,D_s,N_c,N_s)`$ are quantum numbers of the excitation. To be specific, $`I_c`$ counts the total number of holes, $`I_s`$ the number of holes with respect to the up spins, $`D_c`$ ($`D_s`$) the number of holons (spinons) being transfered from one Fermi point to the other and $`N_c`$ ($`N_s`$) the number of charge (spin) particle-hole excitations. The Fermi momentum $`k_{F()}`$ for up (down) spin electrons is given by
$$k_{F()}=\frac{\pi }{2}(n\pm 2m)$$
(12d)
where $`m`$ is the magnetization. If no magnetic field is applied then $`k_F=k_F=k_F`$. As shown in Fig. 1 the density in the zero temperature limit for $`\mu =1.4`$ is approximately $`n_{T0}0.524`$ and therefore $`k_F0.823`$. In the inset of Fig. 5 the circles denote $`k=\pm 2k_F`$ at zero temperature indicating that the leading d-d and s-s CL’s with incommensurate oscillations correspond to the $`2k_F`$-oscillating part. From Eq. (12) we derive

$`n(r)n(0)`$ $`=`$ $`\text{const.}+A_0r^2+A_2r^{\alpha _c}\mathrm{cos}(2k_Fr)`$ (13a)
$`+A_4r^{\beta _c}\mathrm{cos}(4k_Fr)`$
for the equal-time d-d correlation, where $`A_i`$ are matrixelements. The non-oscillating part is due to the lowest particle-hole excitation $`(0,0,0,0,1,0)`$, the $`2k_F`$ part due to a $`(0,0,\pm 1,1,0,0)`$ excitation, while the $`4k_F`$ part arises from the excitation $`(0,0,\pm 1,0,0,0)`$. Thus, the critical exponents can be calculated from Eq. (12b,12c) leading to
$$\alpha _c=1+\frac{\xi _c^2(Q)}{2},\beta _c=2\xi _c^2(Q).$$
(13b)
The equal-time s-s correlation has the same form as Eq. (13a), but the constant as well as the $`4k_F`$ part are absent and the matrixelements are different. The critical exponent of the $`2k_F`$ part is the same but the corresponding excitation is now $`(0,0,0,\pm 1,0,0)`$.<sup>1</sup><sup>1</sup>1The excitations $`(0,0,\pm 1,1,0,0)`$ and $`(0,0,0,\pm 1,0,0)`$ lead to the same critical exponent but they can be distinguished by applying a magnetic field because the first shows $`2k_F`$-oscillations whereas the second oscillates with $`2k_F`$. As shown in Ref. Kawakami and Yang, 1990 the dressed charge $`\xi _c(Q)`$ varies between $`\sqrt{2}`$ for $`n=0`$ and $`1`$ at half-filling. Thus, $`\alpha _c`$ is always the smallest critical exponent and therefore the $`2k_F`$ parts of the correlation functions dominate at zero temperature.

The results for the correlation functions shown here are absolutely consistent with TLL theory when the identification

$$K_\rho =\frac{\xi _c^2(Q)}{2}$$
(14)
is used. For the Tomonaga-Luttinger model it is also possible to calculate multiplicative logarithmic corrections to the algebraic terms of Eq. (13a) which are not directly accessible within Bethe ansatz and CFT. It turns out that the $`2k_F`$ parts of the d-d and the s-s correlation have different logarithmic corrections given by $`\mathrm{ln}^{3/2}r`$ ($`\mathrm{ln}^{1/2}r`$) for the d-d (s-s) correlation.Giamarchi and Schulz (1989)

The long-distance asymptotics of correlation functions at small finite temperature can still be obtained from conformal invariance by the usual mapping of the complex plane onto a strip of width $`1/T`$ with periodic boundary conditions. Eq. (12) is then replaced by

$`\text{e}^{\text{i}(2\pi 2k_F2k_F)D_cx}\text{e}^{\text{i}(2\pi 2k_F)D_sx}`$ $`\left({\displaystyle \frac{\pi T}{v_c\mathrm{sinh}(\pi T(x\text{i}v_c\tau )/v_c)}}\right)^{2\mathrm{\Delta }_c^+}\left({\displaystyle \frac{\pi T}{v_c\mathrm{sinh}(\pi T(x+\text{i}v_c\tau )/v_c)}}\right)^{2\mathrm{\Delta }_c^{}}`$ (15)
$`\times `$ $`\left({\displaystyle \frac{\pi T}{v_s\mathrm{sinh}(\pi T(x\text{i}v_s\tau )/v_s)}}\right)^{2\mathrm{\Delta }_s^+}\left({\displaystyle \frac{\pi T}{v_s\mathrm{sinh}(\pi T(x+\text{i}v_s\tau )/v_s)}}\right)^{2\mathrm{\Delta }_s^{}}.`$
Thus, all CL’s diverge in the low temperature limit as

$$\xi =\frac{1}{2\pi T\left(\frac{x_c}{v_c}+\frac{x_s}{v_s}\right)}=\frac{\gamma }{T}$$
(16)
where the coefficient $`\gamma =v_s/2\pi `$ ($`\gamma =v_c/2\pi `$) for the non-oscillating part of the s-s (d-d) CL,

$$\gamma =\frac{2v_c}{\pi \left(2\frac{v_c}{v_s}+\xi _c^2(Q)\right)}$$
(17a)
for both $`2k_F`$ parts whereas
$$\gamma =\frac{v_c}{2\pi \xi _c^2(Q)}$$
(17b)
for the $`4k_F`$ part of the d-d CL. Contrary to zero temperature the answer to the question which part dominates depends not only on the scaling dimensions $`x_{s,c}`$ but also on the velocities $`v_{s,c}`$! In Fig. 5 the values for these $`\gamma `$’s are denoted by symbols at $`T=0`$ showing good agreement for the non-oscillating parts between the numerics and the CFT results. The $`4k_F`$-oscillating d-d CL is so small that it is difficult to obtain numerically. From Eq. (17a) the $`2k_F`$ s-s and d-d CL’s are expected to be equal in the low-temperature limit. However, the numerical result shows that they are well separated even at the lowest accessible temperatures. This is a consequence of the different logarithmic corrections stated above. Quantitatively we can regard a multiplicative logarithmic term $`\mathrm{ln}^\alpha r`$ as an effective, distance dependent correction of the scaling dimension:

$$x^{}=x\frac{1}{2}\frac{\mathrm{ln}(\mathrm{ln}^\alpha r)}{\mathrm{ln}r}$$
(18)
The relevant length scale is given by the correlation length at the considered temperature, $`r\xi (T)`$. Using this correction together with Eq. (16) leads to splitting of the $`2k_F`$ s-s and d-d CL at finite temperature and to an excellent agreement with the numerical results (see triangle downs in Fig 5).

Next, we regard the singlet pair CL’s shown in Fig. 6.

The CL leading over the entire temperature range is non-oscillating, whereas the next-leading shows incommensurate oscillations. Again the zero temperature result from Eq. (12d) identifies this contribution as the $`2k_F`$ part. The form of the algebraically decaying singlet pair correlation $`P_s(r,0)`$ at zero temperature is given by CFT

$$P_s(r,0)=C_0r^{\delta _p}+C_2r^{ϵ_p}\mathrm{cos}(2k_Fr).$$
(19a)
The non-oscillating term is due to the excitation $`(2,1,\pm 1/2,1,0,0)`$ leading to an exponent
$$\delta _p=1+\frac{2}{\xi _c^2(Q)}$$
(19b)
and the $`2k_F`$ part due to the excitation $`(2,1,\pm 1/2,0,0,0)`$ leading to
$$ϵ_p=\frac{2}{\xi _c^2(Q)}+\frac{\xi _c^2(Q)}{2}.$$
(19c)
Thus, the coefficient $`\gamma `$ in Eq. (16) is now

$$\gamma =\frac{v_c}{\pi \left(\frac{2}{\xi _c^2(Q)}+\frac{v_c}{v_s}\right)}$$
(20a)
for the non-oscillating part, whereas
$$\gamma =\frac{v_c}{\pi \left(\frac{2}{\xi _c^2(Q)}+\frac{\xi _c^2(Q)}{2}\right)}$$
(20b)
for the $`2k_F`$ part. The value for the non-oscillating part is again in good agreement with the numerical results whereas the numerical accessible temperature range is not sufficient to compare the $`2k_F`$ part with CFT (see symbols in Fig. 6).

To illustrate the dependence of the leading s-s and singlet-pair CL on electron density, there are plots for three different chemical potentials in Fig. 7.

The s-s CL is largest for the Heisenberg case ($`n=1`$) and gets suppressed, as expected, with an increasing concentration of holes. On the other hand the singlet-pair CL is nearly independent of the particle density. The s-s CL shows a non-analyticity at a temperature $`T_c`$, where the $`k=\pi `$ part is crossed by the non-oscillating part. With increasing particle density this non-analyticity is shifted to lower temperatures and disappears for chemical potentials $`\mu >0.47`$. Therfore no singularity is visible in the Heisenberg limit of the t-J model. $`T_c`$ as a function of the chemical potential shows an algebraic behavior. A very similar phenomena is described in Ref. Klümper et al., 2001 where the s-s CL of the $`1/2`$-XXZ chain in a magnetic field is investigated. In this case the magnetic field plays the same role as the chemical potential here. However, the singularity there is associated with a crossover from commensurate oscillations at $`T>T_c`$ to incommensurate oscillations at $`T<T_c`$. Recently, temperature dependent CL’s have also been studied analytically in a generalized Hubbard model related to the supersymmetric t-J model.Sakai and Klümper (2001)

### III.2 Static correlation functions

Within the TMRG imaginary time correlation functions are easy to calculate.I.Peschel et al. (1999) The fatal point is that the analytical continuation of the imaginary time results to real times is an ill-posed problem. We therefore want to restrict ourselves to the calculation of static correlation functions defined by

$$G(r,z=0)=_0^\beta 𝑑\tau G(r,\tau )$$
(21)
where $`G(r,\tau )`$ is the correlation function for distance $`r`$ and imaginary time $`\tau `$. Note, that the static correlation function is calculated to higher accuracy than the equal-time correlation. Nevertheless the spatial CL’s are identical in both cases except for the trivial correlations caused by particle-hole excitations which have vanishing matrixelements in the static case (see Eq. 15). The static correlation function can be expressed through the largest eigenvalue and corresponding eigenvectors of the QTM:

$$G(r,z=0)=\frac{ϵ}{M\mathrm{\Lambda }_0^{r+1}}\mathrm{\Psi }_0^L|\stackrel{~}{T}_MT_M^{r1}\stackrel{~}{T}_M|\mathrm{\Psi }_0^R$$
(22)
for distances $`r1`$ with the definition

$$\stackrel{~}{T}_M=\underset{k=0}{\overset{M}{}}T_M(A_{ϵk})$$
(23)
where $`T_M(A_{ϵk})`$ denotes the usual transfer matrix $`T_M`$ where the considered operator $`A`$ is added at imaginary time position $`\tau =ϵk`$ and

$$G(r=0,z=0)=\frac{ϵ}{\mathrm{\Lambda }_0}\mathrm{\Psi }_0^L|\widehat{T}_M|\mathrm{\Psi }_0^R$$
(24)
with

$$\widehat{T}_M=\underset{k=0}{\overset{M}{}}T_M(A_0,A_{ϵk}).$$
(25)
These equations are valid within the alternative mapping to the classical model with alternating rowsSirker and Klümper (2002) where a local transfer matrix covers only one site (see Sec. I). Note, that these formulas become more complicated if the usual decomposition of the Hamiltonian into even and odd parts is applied, because in this case even and odd distances as well as distance $`r=1`$ have to be treated separately. To check the accuracy of the TMRG algorithm we have calculated the static density-density correlation function for free spinless fermions (see appendix A) showing that the TMRG results are reliable even at short distances $`r`$.

In Fig. 8 (Fig. 9) the static longitudinal s-s correlation function for $`J/t=2.0`$, $`\mu =1.4`$ and $`T=2.0`$ ($`T=0.1`$) is shown. From Fig. 5 we expect $`\pi `$-oscillations of the correlation function for $`T=2.0`$ and $`2k_F`$-oscillations for $`T=0.1`$ because the non-oscillating part does not show up in a static correlation function as explained before.

By using $`S_0^zS_r^z=A\mathrm{exp}(r/\xi )\mathrm{cos}(kr+\delta )`$ as a fit function we have estimated the CL’s as well as the wavevectors for both temperatures directly from the correlation functions (see legends in Fig. 8 and Fig. 9). In both cases the results coincide within errors of the order $`10^4`$ with that obtained by using directly the eigenvalues of the QTM. The perfect coincidence of the CL’s estimated directly from the eigenvalues and that from a fit of the correlation function is not astonishing. In both cases the same QTM is used in the numerics where several eigenvalues are calculated directly using diagonalization routines in the first case whereas only the largest eigenvalue and the corresponding eigenvectors are needed to calculate the correlation function by some time-consuming matrix-vector multiplications in the second case. However, it provides a good consistency check of the numerics.

## IV The physically relevant region

As mentioned in the introduction the t-J model is a simplification of the Hubbard model obtained in the limit $`Ut`$. Thus, one is especially interested in the case where the spin exchange coupling $`J=2t^2/|U|`$ is smaller than the hopping $`t`$. The TMRG is applicable for any value of $`J/t`$ with almost the same accuracy and our algorithm remains stable even in the limiting cases $`J/t1`$ or $`J/t1`$. Problems as reported by Ammon et al. Ammon et al. (1999) never occurred in our algorithm. To be specific, we never observed that the algorithm breaks down after performing only a few or more RG-steps although we have retained up to 240 states to calculate the correlation lengths accurately. In the following calculations we set $`J/t=0.35`$, a value often used in the literature.

When regarding the specific heat shown in Fig. 10 a two-peak structure is again obvious in the high and medium density case and there seems to be a linear regime consistent with TLL theory, although we have noticed in the supersymmetric case that numerical errors grow up for low temperatures. The structure of the spin and charge susceptibilities (see Fig. 11) looks also rather similar to the supersymmetric case.

The arguments given for the high and low temperature regime of the susceptibilities are also applicable here. It is easy to understand that the absolute values of the spin susceptibility are now larger, because the weaker antiferromagnetic coupling $`J`$ suppresses less efficiently a ferromagnetic alignment of spins. In principle, the qualitative discussion of $`c`$, $`\chi _s`$ and $`\chi _c`$ is very similar to that of the integrable case. There are, however, two modifications. First, in the charge susceptibility for low densities we observe two maxima. We do not have any simple explanation of this fact. Note that a similar observation has been made for the Hubbard chain. Second, due to the different energy scales of spinon and incoherent single particle motion at any density there are two maxima in $`c(T)`$ even at densities close to 1.

We also calculated again several s-s and d-d CL’s as shown in Fig. 12. The CL’s are now often difficult to distinguish and the crossing of CL’s is “smoother” and shifted to higher temperatures when compared with the supersymmetric case. This clearly shows that the spin-exchange interaction is responsible for the various crossovers of CL’s and in Appendix A we show that all such crossovers vanish for free spinless fermions. It is also striking that incommensurate oscillations are now present even at high temperatures. The leading s-s CL shows oscillations which are identified by Eq. (12d) as $`2k_F`$ whereas one of the d-d CL shows $`4k_F`$ oscillations. The other d-d CL from Fig. (12) is non-oscillating and leading in the low-temperature limit. These kind of oscillations are expected from TLL theory (see Eq. (13a)) and we also note that the CL’s seem to diverge again as $`\xi 1/T`$ for low temperatures. This provides additional evidence that the t-J model at $`J/t=0.35`$ belongs to the same universality class (TLL) as the supersymmetric model.

## V Phase separation

At $`J/t`$ large, the attractive Heisenberg interaction in Eq. (I) dominates the kinetic energy term. In a canonical ensemble the model therefore phase separates into a high density and a low density region in order to optimize the Heisenberg energy. In our grand-canonical description of the model, however, there remains a competition between the chemical potential term and the Heisenberg exchange energy. If we ignore the kinetic energy completely - what is exact if total phase separation occurs - a simple picture evolves. Because the ground state energy of a $`s=1/2`$-Heisenberg chain per particle is given by $`J\mathrm{ln}2`$, we expect in the limit $`T0`$ an empty state if $`\mu <J\mathrm{ln}2`$ and a state with $`n=1`$ if $`\mu >J\mathrm{ln}2`$. Therefore phase separation can only be present at one special point characterized by $`\mu J\mathrm{ln}2`$.

In the Figures 13(a), 13(b), 13(c) the density at a constant temperature is shown as a function of the chemical potential for three different parameters $`J/t`$. Fig. 13(a) shows that at $`J/t=3.1`$ the compressibility $`n/\mu `$ in the limit $`T0`$ is not diverging, indicating that the phase separated region is not yet reached. At $`J/t=3.3`$ (Fig. 13(b)) the density jumps from $`n=0`$ to $`n0.8`$ at $`\mu J\mathrm{ln}2`$ in the limit $`T0`$ as expected for the phase separated region. The phase separation is here between the empty and an electron rich state. At $`J/t=3.5`$ we find full phase separation, indicated by the jump of the density from 0 to 1 for $`T0`$ (see Fig. 13(c)). Our calculations confirm that the fully phase-separated state is destroyed by introducing holes into the Heisenberg chain island as stated by Ogata et al. Ogata et al. (1991)

## VI Luther-Emery Phase

As already mentioned in the introduction a phase with a spin gap is expected at least for low densities and values of $`J/t`$ close to phase separation. Because we expect the spin gap to be caused by a marginal operator leading to an exponentially small gap, the TMRG is not suited to determine the phase boundaries. However, in regions where the spin gap is larger than the lowest accessible temperature the TMRG can show the existence of such a phase without using any additional assumptions. A spin gap is directly visible in the spin susceptibility going to 0 for $`T0`$.

Therefore, we have calculated the spin susceptibility for two different values of $`J/t`$ close to phase separation and have chosen chemical potentials so that the density is given by $`n0.2`$ for low temperatures (see Fig. 14). In both cases a very small spin gap appears. The quadratic dispersion of a gapped 1D system leads to $`\chi \mathrm{exp}(\mathrm{\Delta }/T)/\sqrt{T}`$ for the low-temperature asymptotics. Using this function for a fit of the numerical data, we find $`\mathrm{\Delta }=0.05\pm 0.01`$ in both cases.

Another proof of LE-properties of the t-J model is given by the calculation of s-s and d-d CL’s. In Fig. 15 these CL’s are shown as usual as temperature times CL versus temperature. A spin gap is connected with finite s-s CL’s $`\xi _s`$ and we therefore expect that $`T\xi _s0`$ for $`T0`$. On the other hand the d-d CL’s $`\xi _n`$ are not affected and should still diverge as $`\xi _n1/T`$ for low temperatures. This picture seems to be consistent with the numerical results. However, we are not able to present numerical data for lower temperatures, which could support this scenario further.

## VII t-J Model with Ising anisotropy

In this section we want to consider the t-J model with modification of Hamiltonian (I) by replacing

$$𝑺_i𝑺_{i+1}S_i^xS_{i+1}^x+S_i^yS_{i+1}^y+\mathrm{\Delta }S_i^zS_{i+1}^z\frac{\mathrm{\Delta }}{4}n_in_{i+1}$$
(26)
where $`\mathrm{\Delta }>1`$. We expect that such an anisotropy enhances superconducting correlations relative to d-d correlations as has been explicitly shown in an exactly solvable anisotropic t-J model by Bariev et al. Bariev et al. (1995) However, their version includes some unphysical parity breaking terms (keeping $`PT`$-symmetry) making it interesting to investigate if the same is true for the model defined here. For the Heisenberg chain it is known that an Ising-like anisotropy promotes long range spin order and causes a spin gap. The situation is more complicated in the t-J model: The charge sector is unaffected and the charge excitations remain critical (i.e. d-d correlations decay algebraically in the ground state). The expected long range spin order could be hidden when the density $`n1`$, because the spin operators act on the physical lattice whereas the spins are coupled to the “electron lattice”. The long range order would then be visible only in a string order parameter $`\stackrel{~}{S}_0^z\stackrel{~}{S}_r^z`$ where all empty sites between $`0`$ and $`r`$ are omitted. This has been emphasized also by Pruschke and Shiba Pruschke and Shiba (1991) who have studied the limit $`J/t0`$. In any case there has to be a spin gap also in the t-J model if $`\mathrm{\Delta }>1`$. However, in the Heisenberg chain the gap is given by $`\mathrm{\Delta }E\mathrm{exp}\left\{J/(\mathrm{\Delta }1)\right\}`$ leading to an exponential small gap and this gap is further reduced in the t-J model for densities $`n<1`$ making it often undetectable in the thermodynamical data.

In Fig. 16 a spin gap $`\mathrm{\Delta }_s=0.054\pm 0.01`$ is visible in the susceptibility data for $`\mu =1.5`$ corresponding to $`n0.98`$ in the low-temperature limit. This is in agreement with the gap $`\mathrm{\Delta }^{HB}=0.043`$ for the Heisenberg chain with the same anisotropy. For $`\mu =1.8`$ ($`n_{T0}0.54`$) and $`\mu =2.0`$ ($`n_{T0}0.27`$) the spin gap is so small that it is not visible in the accessible temperature range.

The specific heat (Fig. 17) looks qualitatively similar to the isotropic case. Again two peaks (shoulders) are visible corresponding to spinon and holon excitations.

In Fig. 18 several crossovers in the CL’s are visible. In the low-temperature regime the leading d-d and longitudinal s-s CL’s diverge as $`\xi 1/T`$ showing that the charge excitations remain critical and that the long-range spin order is hidden. As expected the singlet pair CL dominates at low temperatures over the d-d CL supporting that an Ising-like anisotropy could be a mechanism leading to superconductivity in higher dimensions.Bariev et al. (1995)

## VIII Conclusions

We have described briefly two variants of the TMRG algorithm, where the novel one allows to determine wavevectors unambiguously solely from the eigenvalues of the QTM and makes the calculation of static correlation functions much easier. At the supersymmetric point of the one-dimensional t-J model an excellent agreement between the numerical results and Bethe ansatz was found. Additionally, the numerical results for correlation lengths and static correlation functions have been successfully tested for free spinless fermions. This shows that the TMRG is well suited even for fermionic systems and yields accurate results over a wide temperature range what was sometimes questioned before.Ammon et al. (1999) In particular, we have concentrated on the calculation of correlation lengths at the supersymmetric point to study the crossover from the non-universal high T lattice into the quantum critical regime ($`Tt`$). The non-universal regime is characterized by various crossovers between CL’s with different wavevectors whereas the CL’s are non-crossing and diverging as $`1/T`$ in the universal TLL regime. A good coincidence between predictions for the low-temperature asymptotics of CL’s by CFT and the numerical results was shown, but it was important to take also the logarithmic corrections into account. For $`J/t=0.35`$, a value often considered physically relevant, we observed properties that were rather similar to the supersymmetric case. In particular, the t-J model at this parameter point belongs also to the TLL universality class. In our grand-canonical description phase separation has the meaning of phase coexistence at one chemical potential $`\mu `$. This is due to the remaining competition between the Heisenberg and the chemical potential term. If $`J/t`$ is large enough so that total phase separation occurs this competition leads to an empty state for $`T=0`$ if $`\mu <J\mathrm{ln}2`$ and to a state with $`n=1`$ if $`\mu >J\mathrm{ln}2`$. Our data also show that the fully phase separated state is destroyed by introducing holes into the Heisenberg island. By calculating directly the spin susceptibility for small densities and values of $`J/t`$ near phase separation, we have proven the existence of a spin-gap phase without making additional assumptions. This was also supported by the calculation of s-s and d-d CL’s, indicating that all spin CL’s are finite at zero temperature, whereas the d-d CL’s are unaffected and diverging still as $`1/T`$. Finally we have studied the t-J model with an additional Ising-like anisotropy and have shown that singlet pair correlations are enhanced so that a tendency towards superconductivity is expected in higher dimensions. Further on, we have settled that the expected long range spin order is hidden away from half-filling, because the spins are coupled to the “electron lattice” whereas the operators act on the physical lattice.

###### Acknowledgements.

J.S. thanks Ch. Scheeren for communicating his results from Bethe ansatz and is also grateful for valuable discussions with A. Kemper. This work is supported by the DFG in SP1073.

## Appendix A Free spinless fermions

Here we calculate the density-density correlation function for a system of free spinless fermions on a lattice. The Hamiltonian is given by

$$H=t\underset{q}{}\mathrm{cos}(q)c_q^{}c_q$$
(27)
and the density operator by

$$\rho _q(t)=\frac{1}{V}\underset{k}{}c_{k+q/2}^{}(t)c_{kq/2}(t).$$
(28)
By applying Wick’s theorem the d-d correlation function separates into two-point functions

$`c_r^{}(t)c_r(t)c_0^{}(0)c_0(0)={\displaystyle \frac{1}{V^2}}\left({\displaystyle \underset{k}{}}n(ϵ_k)\right)^2`$
$`+`$ $`{\displaystyle \frac{1}{V^2}}{\displaystyle \underset{k,q}{}}\text{e}^{\text{i}qr}c_{k\frac{q}{2}}^{}(t)c_{k\frac{q}{2}}(0)c_{k+\frac{q}{2}}(t)c_{k+\frac{q}{2}}^{}(0)`$
where $`ϵ_k=t\mathrm{cos}(k)`$ and $`n(ϵ_k)`$ is the Fermi function. By using the Matsubara formalism, inserting the known results for the electrons Green’s function and then going back to real time, this is easily transformed into

$`c_r^{}(t)c_r(t)c_0^{}(0)c_0(0)=n^2`$
$`+`$ $`{\displaystyle \frac{1}{V^2}}{\displaystyle \underset{q,k}{}}\text{e}^{\text{i}qr}{\displaystyle \frac{\text{e}^{\text{i}(ϵ_{k\frac{q}{2}}ϵ_{k+\frac{q}{2}})t}}{(\text{e}^{\beta ϵ_{k\frac{q}{2}}}+1)(\text{e}^{\beta ϵ_{k+\frac{q}{2}}}+1)}}.`$
Regarding only the connected part and transforming to Matsubara frequencies we get in the free fermion case

$$G(r,z=0)=\frac{1}{(2\pi )^2}𝑑q𝑑k\frac{n(ϵ_{k+\frac{q}{2}})n(ϵ_{k\frac{q}{2}})}{ϵ_{k\frac{q}{2}}ϵ_{k+\frac{q}{2}}}\mathrm{cos}(qr)$$
(31)
for the static density-density correlation function as defined in Eq. (21).

In Fig. 19 we show results for two characteristic temperatures obtained by TMRG using Eqs. (22) and (24) in comparison to the exact results from a numerical evaluation of Eq. (31). The absolute error is largest for short distances and decreases exponentially with distance. However, even at $`T/t=0.1`$ the error always remains smaller than $`10^2`$ showing that TMRG yields useful results not only for the asymptotic behavior of static correlation functions but also at short distances.

Returning to Eq. (A) and setting $`t=0`$ we find

$`c_r^{}c_rc_0^{}c_0`$ $`=`$ $`\left[{\displaystyle \frac{1}{2\pi }}{\displaystyle _\pi ^\pi }𝑑q{\displaystyle \frac{\text{e}^{\text{i}qr}}{\text{e}^{\beta \mathrm{cos}(q)}+1}}\right]`$ (32)
$`\times `$ $`\left[{\displaystyle \frac{1}{2\pi }}{\displaystyle _\pi ^\pi }𝑑k{\displaystyle \frac{\text{e}^{\text{i}kr}}{\text{e}^{\beta \mathrm{cos}(k)}+1}}\right].`$
Now we can in principle evaluate these integrals by closing the integration path in the upper (lower) half plane and using Cauchy’s formula. The poles are at the points

$$q,k=\pm \frac{\pi }{2}\pm \text{i}\text{arcsinh}\left((2n+1)\frac{\pi }{\beta }\right)$$
(33)
and thus, the dominant contribution is

$$c_r^{}c_rc_0^{}c_0M\text{e}^{r/\xi }(1\mathrm{cos}(\pi r))$$
(34a)
where the matrixelement is given by
$$M=\frac{2}{\pi ^2+\beta ^2}$$
(34b)
and the CL by
$$\xi =\frac{1}{2\text{arcsinh}(\pi /\beta )}.$$
(34c)
In Fig. 20 the leading CL received by Eq. (34c) is shown together with TMRG results calculated by Eq. (5). A good agreement is obtained with errors of the numerical data remaining smaller than $`10^3`$ for temperatures down to $`T/t0.1`$.


format: text

text:
# Deuterium Toward Two Milky Way Disk Stars: Probing Extended Sight Lines with the Far Ultraviolet Spectroscopic Explorer

## 1 Introduction

The abundance of deuterium in the interstellar medium (ISM) is a strong constraint on chemical evolution scenarios. Unlike most elements, the net cosmic abundance of deuterium is not enhanced by stellar processes or supernovae, so its abundance has been declining from the primordial value due to destruction by stellar processing. The primordial abundance of deuterium is itself an important quantity, as it can be used to measure $`\eta `$, the universal ratio of baryons to photons, which can be compared to the predictions of Big Bang nucleosynthesis (Reeves et al., 1973; Boesgaard & Steigman, 1985; Burles, Nollett, & Turner, 2001).

One approach to finding the initial deuterium abundance is to measure it in high-redshift QSO absorption-line systems, assuming that the clouds which give rise to these systems have undergone very little stellar processing (Webb et al., 1997; Tytler et al., 1999; O’Meara et al., 2001; Pettini & Bowen, 2001; Levshakov et al., 2002). However, the large range in D/H values measured in these systems suggests that the gas in some of these clouds may in fact have undergone some stellar processing (Prantzos & Ishimaru, 2001; Fields et al., 2001). It is therefore crucial to understand how chemical evolution affects the abundance of deuterium relative to other elements. Several groups are measuring the abundance of deuterium and other elements along many sight lines in the Galactic ISM in order to establish the mean D/H ratio in the Galactic disk, and to look for variations in the D/H ratio. These variations can shed light on the effects of star-formation history and mixing processes on the D/H ratio. The local D/H value can then be compared to regions that have undergone varying amounts of stellar processing (e.g., high velocity clouds, or intergalactic H I clouds), which will illuminate the effects of stellar processing on the D/H ratio in these systems.

The Copernicus satellite was used to measure the D I and H I Lyman series absorption toward $`20`$ stars, and the results suggested a D/H ratio of 1.5$`\times `$10-5, with a dispersion that may reflect real variations in the ISM (see Vidal-Madjar & Gry 1984 for a review). Several measurements of D/H in the local ISM have been made with the Hubble Space Telescope (HST) (Linsky et al., 1995; Linsky & Wood, 1996; Linsky, 2003; Lemoine et al., 1996; Vidal-Madjar et al., 1998; Hébrard et al., 1999; Sahu et al., 1999; Vennes et al., 2000), though only the Ly$`\alpha `$ transitions of D I and H I could be analyzed, limiting the number of available sight lines to those with low enough H I columns that the D I absorption is not blended with the H I line but large enough H I columns that D I can be detected. Three sight lines observed with the Interstellar Medium Absorption Profile Spectrograph (IMAPS) spanned larger distances than the HST sight lines (250 to 500 pc versus $``$ 100 pc), and like the Copernicus results they showed evidence of variation in D/H with location (Jenkins et al., 1999; Sonneborn et al., 2000).

One of the science goals of the Far Ultraviolet Spectroscopic Explorer (FUSE) is to measure D/H in the Milky Way (Moos et al., 2000). FUSE is uniquely suited to do this because it observes the spectral region from 905 to 1187 Å which includes all of the D I Lyman transitions except Ly$`\alpha `$. The ability to use multiple D I transitions means that FUSE is able to measure $`N`$(D I) along higher column density sight lines than missions that only observe Ly$`\alpha `$ (e.g., HST), and for many more sight lines than short term missions (e.g., IMAPS). Moos et al. (2002, and references therein) reported results for the first set of seven sight lines ranging from 37 to 179 pc, five of which had both D I measurements (from FUSE) and H I (from HST). The results were consistent with a single value of D/H$`=(1.52\pm 0.08)\times 10^5`$ (the uncertainty is $`1\sigma `$ in the mean). However, none of these sight lines extend as far as the three sight lines observed by IMAPS. Analysis of more extended sight lines is crucial for the verification and understanding of variations in the D/H ratio in the ISM.

In this paper we present FUSE measurements of D/H toward two distant stars in the Galactic disk: HD 195965 and HD 191877. These stars are at distances greater than the IMAPS targets, extending the distance and range of $`N`$(H I) for D/H measurements in the ISM. Information about these two sight lines is listed in Table 1. The FUSE, STIS, and ground-based observations and data reduction are described in § 2. Section 3 describes some properties of the sight lines. In § 4 we describe both the curve of growth and the absorption profile fitting analyses of the D I column density toward both stars. In § 5 we discuss the determination of the O I and N I abundances along these sight lines, and in § 6 we re-examine the H I column densities. Section 7 contains a discussion of the results and conclusions.

## 2 The Data

### 2.1 FUSE Observations

Information for the FUSE observations is given in Table 2. The FUSE mission and instrument are described by Moos et al. (2000) and Sahnow et al. (2000). The spectra were obtained through the large (30′′$`\times `$30′′) LWRS aperture. FUSE consists of 4 co-aligned optical channels, two optimized for longer wavelength (LiF1 and LiF2; 1000 – 1187 Å) and two optimized for shorter wavelengths (SiC1 and SiC2; 905 – 1100 Å). All of the deuterium lines used in this paper are located in the SiC channels. The Ly$`\beta `$ D I line was not detected in either spectrum because it was blended with the strong H I line, but some lines of other elements used in our analysis are located in the LiF channels. Between exposures the focal plane assemblies were moved so that the spectra sampled different portions of the detector. This shifting lessens the effects of detector defects when the individual spectra are combined.

The raw spectra were processed through the FUSE calibration pipeline (CALFUSE v2.1.6) at the Johns Hopkins University. The pipeline screens data for passage through the South Atlantic Anomaly and low Earth limb angle pointings and corrects for thermal drift of the gratings, thermally-induced changes in the detector read-out circuitry, and Doppler shifts due to the orbital motion of the satellite. Finally, the pipeline subtracts a constant detector background and applies wavelength and flux calibration. After the pipeline reduction the individual spectra were co-added to produce the final calibrated spectrum. The data from the four channels were analyzed separately, because there are slight differences in the spectral resolution between channels, which can cause problems if the channels are co-added. This has the benefit of allowing at least two measurements of each line, a further safeguard against detector defects.

The presence of bright airglow lines at specific regions of the spectrum has resulted in a long-term depletion of electrons in the detector in those regions, and this can cause localized shifts in the spectrum in the dispersion direction. This effect, called “x-walk,” may change the shapes and areas of absorption lines. The effect is more pronounced in certain FUSE detectors. We have avoided using the equivalent widths of absorption lines in regions where the x-walk effect is pronounced, so these detector artifacts do not affect the analysis in this paper. The x-walk distortion is described in the FUSE Observer’s Guide available at the FUSE website (http://fuse.pha.jhu.edu/).

Figure 1 shows a portion of the FUSE spectrum of HD 195965, and Figure 2 shows the same spectral region for HD 191877. At $`\lambda 920`$ Å, the shortest wavelengths used in this paper, the signal to noise ratio is $``$30:1 per $``$20 km s-1 resolution element in both spectra, and is considerably higher at longer wavelengths.

### 2.2 STIS Observations

We observed HD 195965 with the Space Telescope Imaging Spectrograph (STIS) on the HST on 2001 October 9 as part of observing program GTO-8487 (P.I.: Moos). Brief information about the observation is given in Table 2. The spectrum was obtained in ACCUM mode using the FUV MAMA, the E140H grating ($`\lambda _c=1307`$ Å), and the $`0.1\mathrm{}\times 0.03\mathrm{}`$ aperture. We reduced the data using CALSTIS (v2.12a) with all of the normal reduction procedures enabled in the standard data pipeline, including the two-dimensional scattered light removal algorithms (see Landsman & Bowers 1997; Lindler & Bowers 2000; see also Howk & Sembach 2000). A description of the STIS pipeline processing steps can be found in the $`HST`$ Data Handbook (Brown et al., 2002).

The extracted STIS spectrum has $`S/N2530`$ per spectral resolution element. The spectral resolution of $`2.6`$ km s-1 (R = $`\lambda /\mathrm{\Delta }\lambda 200,000`$) is $`11`$ times better than that available with $`FUSE`$. The absolute wavelength scale is accurate to $`0.51.0`$ pixel, or about 3-6 mÅ ($`0.61.3`$ km s-1; see Leitherer et al. 2001). The post-processed STIS data have very low residual backgrounds. We find that the strongly saturated cores of the O I 1302.169 Å and C II 1334.532 Å lines have residual intensities of $`1`$%.

### 2.3 Ground-based Observations

High-resolution echelle spectra of HD 195965 and HD 191877 were taken with the 2.7 m telescope of the University of Texas McDonald Observatory on 2001 December 29. The “2dcoudé” spectrograph (Tull et al., 1995) was used with 1 camera setup and a Tektronix charge-coupled device (CCD) as the detector, which provided high spectral resolution (R = $`\lambda /\mathrm{\Delta }\lambda `$ 160,000; $`\mathrm{\Delta }`$v $``$ 2 km s-1) and covered the wavelength range 6000 Å to 8600 Å. This high resolution mode yielded spectra containing 12 disjoint echelle orders approximately 40 Å wide, separated by approximately 150 Å at 6000 Å. The inter-order separation increased with wavelength until there was approximately 280 Å between orders at near infrared wavelengths. This setup provided spectra on interstellar K I. The resolution was determined from the full width at half maximum (FWHM) of Thorium-Argon (Th-Ar) calibration lamp lines.

The data were reduced in a standard way utilizing the NOAO SUN/IRAF software (v2.11.3). Dark, bias, and flat field lamp exposures were taken to remove any instrumental effects due to the CCD detector. Comparison spectra were taken periodically throughout the night, typically every two hours. The average bias exposure was subtracted from all raw stellar, comparison (Th-Ar) and flat images. The scattered light was fitted by low order polynomials in both the dispersion and spatial directions, and removed. Pixel to pixel sensitivity variations were removed by dividing the normalized average flat field into the stellar spectra. The extracted spectra were placed on a heliocentric wavelength scale using the Th-Ar comparison spectra and corrections for the Earth’s motion. The spectra were co-added and normalized to unity yielding a final spectrum with high signal to noise, approximately 100:1. Figure 3 shows the spectral region near the K I 7698.964 Å line for both sight lines after correcting the velocities to the Local Standard of Rest (for HD 195965 $`V_{LSR}=V_{helio}+16.7`$ km s-1, and for HD 191877 $`V_{LSR}=V_{helio}+17.4`$ km s-1).

## 3 Properties of the Sight Lines

HD 195965 ($`l`$=85.71, $`b`$=5.00) was classified by Morgan, Code, & Whitford (1955) as type B0V. It is listed by Humphreys (1978) as belonging to the Cygnus OB7 association, which is at a distance of $``$830 pc in the Local spiral arm. Diplas & Savage (1994) give a distance of 794 pc from spectroscopic parallax (with an uncertainty of $`25\%`$). At this distance it would be 69 pc above the Galactic plane, still well within the thin disk, where most of the cool gas is located. The Hipparcos catalog lists an uncertain distance of the star at d$`=520\pm _{120}^{230}`$ pc (ESA, 1997). Diplas & Savage (1994) find log $`N`$(HI)=20.9 based on IUE data.

HD 191877 ($`l`$=61.57, $`b`$=-6.45) is type B1.0Ib (Morgan et al., 1955), and spectroscopic parallax gives a distance of $``$2200 pc (Diplas & Savage, 1994), again with an uncertainty of $`25\%`$. This places the star on the edge of the Sagittarius arm and 250 pc below the plane of the Milky Way, well out of the thin disk (scale height $`100`$ pc, Dickey & Lockman 1990). It is likely that most of the gas along this sight line lies within $`1000`$ pc of the Sun. This is supported by the fact that the H I column toward HD 191877 found by Diplas & Savage (1994) (log $`N`$(HI)=20.9) is similar to that toward HD 195965, even though the distance to HD 191877 is more than twice as great as the distance to HD 195965. The two sight lines therefore likely probe a similar distance range in the disk. The Hipparcos distance to HD 191877 is d$`>1160`$ pc (ESA, 1997).

The ground-based optical spectra of HD 195965 and HD 191877 in Figure 3 reveal that there are multiple clouds along the line of sight to each star. Along both sight lines there appear to be two major components of K I, as well as several weaker components. The STIS spectrum of HD 195965 (Figure 4) confirms the presence of multiple velocity components. Two components are visible in the O I 1355.598 Å profiles, and the S I 1425.030 Å profile shows that the strongest of these is actually two components. The structure in the S I line matches well the K I component structure. The Cl I line shows two additional components at more negative velocities, only one of which is strong enough to appear in the O I 1355.598 Å line. The Cl I line is usually a good tracer of H2 (Jura, 1974), so it appears that the molecular hydrogen velocity structure along the sight line is not complex. The strong O I 1302.168 Å line is completely saturated and shows absorption extending over $`\pm `$30 km s-1, along with damping wings at even higher velocities.

Both sight lines contain a significant amount of H2. Most of the unlabeled features in Figures 1 and 2 are H2 absorption lines. Lines from the lower rotational levels ($`J=0`$,1,2) are mostly saturated. These contain most of the H2 column density, so we did not attempt to determine the H2 abundance. However, based on the Cl I column toward HD 195965, we estimate roughly $`(510)\times 10^{19}`$ cm-2 on that sight line (Jura, 1974). The HD 191877 sight line probably contains a similar amount. There are clearly multiple components along the HD 195965 sight line (revealed by the Cl I line), but they fall within $`\pm `$15 km s-1, so at FUSE resolution the H2 distribution can be treated as a single component (see § 4.1).

The widths of the K I lines in the spectra of both stars (and the STIS lines in the HD 195965 spectrum) are larger than the resolution of the data. This suggests that there are no strong, narrow (unresolved) features on either sight line. We can use the STIS data for HD 195965 to show this more conclusively for that sight line. Figure 5 shows the apparent column density of S I derived from two different lines: 1295.653 Å and 1296.174 Å. Although the oscillator strengths differ by a factor of $``$3 (Morton, 1991), the derived column density is very similar for the two lines, supporting the conclusion that narrow, saturated features are not present toward HD 195965. (See Savage & Sembach 1991 for a discussion of apparent column density profile comparisons.) In other words, while cold clouds are clearly present, they do not contain a significant amount of the material along these sight lines. Most of the column arises in the diffuse ISM.

To test the effects of the lower resolution of FUSE on the analysis of blended lines, we smoothed the absorption profiles of the S I 1295.653 Å and 1296.174 Å lines in the STIS data by 15 km s-1 to approximate FUSE resolution. We then determined the column densities by summing the apparent optical depth over the line profile. The total column density derived from the 1295.653 Å line is $`N`$(S I)=$`(1.21\pm 0.02)\times 10^{13}`$ cm-2, while that derived from the 1296.174 Å line is $`N`$(S I)=$`(1.34\pm 0.05)\times 10^{13}`$ cm-2. Both values are well within the $`1\sigma `$ uncertainties of the column densities derived from the unsmoothed profiles, $`N`$(S I)=$`(1.26\pm 0.08)\times 10^{13}`$ cm-2 from the 1295.653 Å line, and $`N`$(S I)=$`(1.33\pm 0.18)\times 10^{13}`$ cm-2 from the 1296.174 Å line. The measured equivalent widths are $`17.2\pm 0.9`$ mÅ and $`7.4\pm 1.0`$ mÅ, similar to the weaker N I and O I lines along the HD 195965 sight line. We also carried out this test on the K I lines shown in Figure 3. Toward HD 191877, the column densities from the unsmoothed profile, $`(4.02\pm 0.06)\times 10^{11}`$ cm-2, agrees well with the column density derived from the smoothed profile, $`(3.97\pm 0.09)\times 10^{11}`$ cm-2. The test was harder to perform on the HD 195965 sight line because of residuals from poorly removed airglow lines, but the column density derived from the unsmoothed profile, $`(3.24\pm 0.16)\times 10^{11}`$ cm-2, is in marginal agreement with that derived from the smoothed profile, $`(2.91\pm 0.19)\times 10^{11}`$ cm-2. With measured equivalent widths of $`49.7\pm 2.1`$ mÅ and $`62.0\pm 1.0`$ mÅ toward HD 195965 and HD 191877, respectively, these lines are similar in strength to the stronger O I, N I, and D I lines along both sight lines. These tests provide more evidence that there are no unresolved, saturated components along either sight line.

## 4 The D I Column Density

The simplest approach to analyzing the chemical composition along a line of sight is to measure the equivalent widths of all lines of a species, assume the material lies in a single component, and construct a Doppler-broadened curve of growth (COG). Since the FUSE bandpass contains all of the transitions of D I except for Ly$`\alpha `$, in principle this is possible using the current data. In reality, many of the D I transitions are not usable due to interference from neighboring absorption lines. Table 3 lists some of the important interfering lines for the unusable D I transitions along these sight lines. The lower order D I lines are blended with strong H I lines, while higher order lines are often blended with strong H2 lines. However, in the spectra of each sight line there are four D I transitions that were either free of blending or were blended with weak H2 lines that could be removed: $`\lambda `$919.102, $`\lambda `$920.713, $`\lambda `$922.900, and $`\lambda `$925.974. A more complex but potentially more accurate method is to model the absorption lines with a profile fitting program. This technique can be especially advantageous when the velocity structure of the cloud components is complicated. In this section we describe both of these approaches to determining the D I abundances toward HD 195965 and HD 191877.

### 4.1 Molecular Hydrogen Removal

There are many H2 lines in the spectral range around the deuterium lines, and several can potentially affect the measurement of the D I equivalent widths. The lines that affect the four D I lines used in this analysis are listed in Table 4. We include only lines with rotational level J$``$5, because no higher $`J`$ level lines were found in the spectra. The lines listed in Table 4 are expected to be strong enough to affect the D I lines. The 920.803 Å line does not fall within the 920.713 Å D I line, but it could affect the continuum on the red side of the line.

To remove the contamination, a line of the same rotational level in a clean part of the spectrum was found to use as a template, following Jenkins et al. (1999). The template line was chosen from the same detector segment as the line being corrected, in order to avoid differences in the line spread function between segments. The parameters of the template lines are listed in Table 4. A Gaussian function was fitted to the template line, and the fit was scaled in optical depth space by the ratio of the $`f\lambda `$ values to have the strength appropriate for the contaminating line. The template line was chosen to have as close to the same value of $`f\lambda `$ as the contaminating line as possible in order to avoid saturation effects. The scaled fit was then shifted to the wavelength of the contaminating line and divided out of the observed spectrum around the D I lines. Another H2 line close in wavelength to the D I line was used to register the position of the scaled line before dividing it out, in order to compensate for any relative errors in the wavelength scale. Figures 6 and 7 show the corrected profiles (solid lines) along with the uncorrected profiles (dashed lines) for both sight lines as a function of LSR velocity.

### 4.2 Curve of Growth Analysis

Once the H2 lines were removed, it was necessary to determine the location of the continuum at the position of the D I absorption lines. This was done by fitting low order polynomials across the lines. The situation is complicated by the fact that the D I lines fall on the blue damping wings of the associated H I lines. We used the steeply-sloped blue sides of the H I lines to constrain the fits, and when possible used the red wings of the H I lines as a guide in determining the optimal continuum. Figures 6 and 7 show the absorption profiles of the four D I lines after correction for H2 for HD 195965 and HD 191877, respectively, along with the adopted continuum placement.

We measured the equivalent widths of the D I lines using the procedures described by Sembach & Savage (1992). The lines were integrated over a velocity range that depended upon the details of the local continuum and line blending. The equivalent widths are listed in Table 5 for HD 195965 and Table 6 for HD 191877. A theoretical single-component COG was fit to the equivalent widths by minimizing $`\chi ^2`$ while varying the column density $`N`$ and the Doppler parameter $`b`$. The final COGs are shown in Figures 8 and 9 for HD 195965 and HD 191877, respectively. The final uncertainties in the equivalent widths were determined from the dispersion of the points about the COG, assuming that a single component COG is appropriate for the sight lines. These error bars are a combination of statistical uncertainties and continuum fitting errors, as described by Sembach & Savage (1992), and include the effects of fixed-pattern noise.

Placing the continuum near the 925.974 Å line was problematic in the spectra of both stars, because of absorption from other elements on the blue side of the absorption line (see Figures 6 and 7). This results in large uncertainties in the equivalent widths. Because of this the 925.974 Å line is not given as much weight as the other lines in the COG fitting. Table 7 contains the COG D I column density results for both sight lines.

The inset in each COG figure shows a contour plot of the $`\mathrm{\Delta }\chi ^2`$ values in the $`N`$ and $`b`$ dimensions, with contours of the minimum $`\mathrm{\Delta }\chi ^2`$+1, 4, and 9 shown (corresponding to 1, 2, and 3$`\sigma `$). The quoted error bars for $`N`$ shown in Table 7 are the maximum extent of the projection of the 2$`\sigma `$ contour onto the $`N`$ axis, and likewise for $`b`$ (listed in the captions to Figures 8 and 9).

### 4.3 Absorption Profile Fitting

Column densities were also measured independently with the profile fitting procedure Owens.f, following the method described by Hébrard et al. (2002). Briefly, the software simultaneously fits Voigt profiles to all of the unsaturated interstellar lines detected in all FUSE segments, using a $`\chi ^2`$ minimization procedure with many free parameters. The free parameters characterize the physical properties of the interstellar clouds (radial velocities, column densities, temperatures, micro-turbulent velocities) as well as the continuum shapes (modeled with polynomials) and possible instrumental effects (line spread function widths, wavelength solution variations). The error bars are obtained using the standard $`\mathrm{\Delta }\chi ^2`$ method; they include statistical as well as systematic uncertainties (Hébrard et al., 2002).

The fits include the D I lines 916.6 Å, 919.1 Å, and 920.7 Å, the N I lines 951.1 Å, 951.3 Å, and 959.5 Å, as well as Fe II, and molecular H2, HD, and CO lines. No O I lines were included in the fits because they are all saturated. Since the H2 lines were included in the fits, the lines were not removed as for the COG analysis (§ 4.1). Samples of the fits are plotted on Figure 10, which presents two spectral windows for each target. The final fits include 30 and 22 spectral windows for HD195965 and HD191877, respectively. The results for $`N`$(D I) and $`N`$(N I) are reported in Table 7. We obtain $`\mathrm{log}N`$(Fe II)$`=14.81\pm 0.02`$ and $`14.95\pm 0.04`$ ($`2\sigma `$) for HD 195965 and HD 191877, respectively.

## 5 The O I and N I Column Densities

We measured the equivalent widths of O I and N I lines in the FUSE spectra using the same techniques as described for D I. The measured equivalent widths are listed in Table 5 for HD 195965 and Table 6 for HD 191877. Figure 11 shows single-component COGs for O I and N I toward HD 195965. The O I lines in the FUSE bandpass that are not affected by H2 absorption do not cover a large range of oscillator strengths, and they mostly fall on the flat part of the COG. The STIS spectrum of HD 195965 adds two O I lines: the 1355.598 Å line on the linear part of the COG, and the 1302.169 Å line, which exhibits damping wings and is on the square root portion of the COG. The STIS lines cover a larger range of oscillator strengths, so the O I abundance toward HD 195965 is well constrained by the FUSE+STIS COG. The N I lines in the FUSE spectra also cover a large range of oscillator strengths, so an accurate determination of the column density is possible. The derived column densities for O I and N I toward HD 195965 are listed in Table 7.

Figure 12 shows the COGs for O I and N I toward HD 191877. As was the case for HD 195965, the N I abundance is well constrained, but the O I abundance toward HD 191877 is less certain because of the lack of STIS data. A single-component COG fit to the FUSE points yields a column density of log $`N`$(O I)$`=17.24\pm _{0.22}^{0.42}`$ cm-2 and the Doppler broadening parameter $`b=8.2\pm _{0.9}^{0.5}`$ km s-1($`2\sigma `$). We also set a lower limit of 230 mÅ ($`3\sigma `$) on the 1302.169 Å equivalent width using an IUE spectrum of HD 191877 (see § 6). This point rules out the COG found for the FUSE points alone, and restricts the O I column density to log $`N`$(O I)$`17.42`$ cm-2, with $`b=7.85`$ km s-1. The $`b`$-values derived from the COGs for N I and O I toward HD 195965 are almost identical, suggesting that the velocity behavior of these two species is coupled. If we assume that this is also true for the HD 191877 sight line and force $`b=7.5`$ km s-1 for the O I COG (the solid line in the top panel of Figure 12), the derived column density of O I is log $`N`$(O I)=$`17.54\pm _{0.12}^{0.20}`$ cm-2 ($`2\sigma `$). Although the HD 195965 sight line suggests that O I and N I are coupled, we cannot know with certainty whether this is also true along the HD 191877 sight line. This exercise illustrates the fact that the O I column density toward HD 191877 is not well constrained by the lines in the FUSE spectrum. We adopt the larger $`N`$(O I) column density as our preferred value, but note that it should be treated with caution.

We also used the damping wings of the O I 1302.169 Å line in the STIS spectrum of HD 195965 to make an independent determination of the O I abundance toward the star. We fit a single component Voigt profile to the normalized O I line using an oscillator strength f = 0.04887 and a natural damping constant $`\gamma =5.75\times 10^8`$ s-1 (Morton 1991). The best fit model has a central velocity $`v_{helio}=+2.1`$ km s-1, b = $`6.3\pm 0.3`$ km s-1, and $`\mathrm{log}N`$(O I) = $`17.80\pm _{0.09}^{0.05}`$ ($`2\sigma `$). This column density is in excellent agreement with the COG result listed in Table 7. The single-component fit to the O I 1302.169 Å line is shown in Figure 13. The fit shown in the figure reproduces the observed absorption profile core and damping wings well, with the minor exception of a slight underestimate in the absorption strength between +21 and +40 km s-1. This additional absorption is likely due to very low column density gas along the sight line that does not contribute much to the total column density; a component with an O I column density of $`10^{13}`$ cm-2 could account for this minor discrepancy. Experimentation with additional components does not change the best-fit result substantially. Performing a fit with the O I column density divided among two components separated by $``$ 12 km s-1 as suggested by the O I 1355.598 Å line shown in Figure 4 yields a best fit nearly indistinguishable from the single-component fit shown in Figure 13.

The O I 1302.169 Å line is very strong, and can be used to test for the presence of high-velocity components that might contribute to the column density of H I but be outside of velocity integration range for other elements. The most obvious candidate is the component between +21 and +40 km s-1, which contributes only $`N`$(O I)$`10^{13}`$ cm-2, corresponding to $`N`$(H I)$`3\times 10^{16}`$ cm-2 using the O/H ratio of $`3.43\times 10^4`$ from Meyer (2001), or $`3\times 10^5`$ of the total H I column density (see § 6). No other components this strong are visible in Figure 13, so we conclude that high-velocity H I is not an issue for the HD 195965 sight line.

The STIS spectrum of HD 195965 contains both the optically thin 1355.598 Å line and the damped 1302.169 Å line. We can use this fact to test for consistency between the oscillator strengths of the two lines (see Sofia, Cardelli, & Savage 1994). The O I column densities derived from these lines are $`N`$(O I)$`=(5.7\pm 1.1)\times 10^{17}`$ cm-2 and $`(6.3\pm 0.5)\times 10^{17}`$ cm-2, respectively (uncertainties are $`1\sigma `$). These results are identical within the uncertainty, which suggests that any error the oscillator strength of either line is less than 20% (the uncertainty in the 1355.598 Å equivalent width). The oscillator strength for the 1355.598 Å line listed in Table 5 includes the correction from Welty et al. (1999).

## 6 The H I Column Density

We have re-measured the H I column densities for the HD 191877 and HD 195965 sight lines using existing STIS and IUE data. Single-component Voigt profiles were fit to the H I absorption toward each star. The fits were limited to the 1170 to 1260 Å spectral region, with the overall continuum levels determined at longer wavelengths, as needed. For both sight lines, the fit is insensitive to the assumption of a single absorbing component. Adding additional components or changing the $`b`$-value of the absorption (assumed to be in the 1 to 100 km s-1 range) did not alter the fits significantly since the quality of the fits is governed primarily by the steep walls of the line cores and the broad radiation-damping wings. The Voigt profile fitting method has been described in detail for H I determinations in previous FUSE D/H papers (see Moos et al. 2002, and references therein). Our analyses follow many of these same principles and methods.

For HD 195965, we used the STIS E140H data available in this wavelength region (see §2.2) to determine the H I column density. We paid special attention to making sure that the flux calibration mapping at the ends of the individual echelle orders was consistent between adjoining orders. Scattered light is not a significant source of uncertainty in the final measurement since we used the damping wings of the lines rather than the saturated portion of the line core to determine the H I column density. We find a best-fit value of log $`N`$(H I) = 20.95$`\pm `$0.05 ($`2\sigma `$) for the fit to the high quality STIS data shown in Figure 14. This value is slightly higher than, but within the errors of, the former value of 20.90+/-0.09 ($`1\sigma `$) found by Diplas & Savage (1994) using the ”continuum reconstruction” technique on existing IUE data. The reconstruction method and profile fitting methods are complimentary to each other. The fit is shown in the bottom panel of Figure 14, together with the fits appropriate for the $`\pm 2\sigma `$ column density ranges (dashed lines). For this sight line, the higher quality STIS data more tightly defines the extent of the H I absorption than the IUE data, and as can be seen in the top panel of Figure 14, the resulting fit is consistent with the IUE small aperture data analyzed previously.

The primary difficulty in fitting the STIS profile of HD 195965 arises in the wavelength region between 1222 and 1225 Å, where the fit slightly overestimates the amount of absorption present. We explored other models that fit this region better (at the expense of fit quality at other wavelengths) and found that in all cases the resulting column density was within the listed uncertainties. Therefore, while we are uncertain of the exact cause of the discrepancy in this wavelength region, we believe the resulting H I column density is a robust estimate. Note that this portion of the STIS spectrum contains the border between two echelle orders, and that the anomaly is not apparent in the IUE spectrum.

For HD 191877, we re-examined the IUE small aperture data (SWP02837) available for this sight line. A more recent exposure (SWP148425) taken through the small aperture clearly has problems with the flux levels as a result of poor centering; we did not use this low quality spectrum. We find log $`N`$(H I) = 21.05$`\pm `$0.10 ($`2\sigma `$), which again is slightly larger than, but consistent with, the previous estimate of 20.90+/-0.10 ($`1\sigma `$) found by Diplas & Savage (1994). The fit is shown in Figure 15, along with the $`\pm 2\sigma `$ error bars. As can be seen from this figure, a fit with a smaller column, like the one preferred by the previous investigators, is clearly below the true column density (see the $`2\sigma `$ \[log $`N=20.95`$\] model in the figure). Diplas & Savage note that there may be a small amount ($`0.01`$ dex) of contamination by stellar absorption in the H I Ly$`\alpha `$ profile. This contamination does not affect our $`N`$(H I) estimate significantly.

## 7 Discussion

The sight lines analyzed in this paper are the longest ones for which the D I column densities have been determined with FUSE. As such, they present the opportunity to look for variations in D I/H I, D I/O I, and D I/N I beyond the local ISM sight lines presented by Moos et al. (2002). Table 8 contains column density ratios of various elements toward both stars, along with the mean values for the local ISM from Moos et al. (2002). In this section we discuss the behavior of these column density ratios beyond the local ISM.

### 7.1 D/H

The deuterium abundances along seven sight lines have been analyzed with FUSE (Moos et al., 2002; Kruk et al., 2002; Friedman et al., 2002; Sonneborn et al., 2002; Lemoine et al., 2002; Lehner et al., 2002; Wood et al., 2002; Hébrard et al., 2002). Five of these had both D I and H I measurements, and the results were consistent with a constant value of D/H in the local ISM. Earlier Copernicus observations showed evidence for variations in D/H at larger distances (see Vidal-Madjar & Gry 1984 for a review), evidence which was strengthened by subsequent IMAPS investigations (Jenkins et al., 1999; Sonneborn et al., 2000). The main obstacle to verifying and understanding this variation is the paucity of suitable deuterium sight lines that probe distances $``$250 pc.

Figure 16 shows the D/H ratios for HD 195965 and HD 191877 plotted as a function of distance from the Sun, along with five of the seven original FUSE sight lines with reliable H I measurements, and several others from HST, IMAPS, and Copernicus (see Moos et al. 2002). The D/H ratios toward both stars are lower than the average value found for the original seven FUSE sight lines. This contrasts with the lack of variation seen in the local ISM. The new sight lines support the existence of variations in the D/H ratio at distances larger than a few hundred pc, as also suggested by IMAPS andCopernicus.

In the simplest model of chemical evolution, stellar processing (astration) leads to a decline in the D/H ratio because D is easily destroyed in stars. Thus, the D/H ratio is expected to vary among regions of the Galaxy that have undergone different amounts of star formation in the past, provided that the gas in these regions has not mixed (e.g. Tosi et al. 1998; Chiappini et al. 2002). In a small volume like the local ISM ($`r=100`$ pc), Moos et al. (2002) showed that the time between supernova $`t_{SN}`$ is roughly equal to the mixing time scale (assumed to be the sound crossing time $`t_s`$), suggesting that the gas in such a region is well-mixed. If a larger region is considered, say 1000 pc in radius to approximate the scales probed by HD 195965 and HD 191877, the crossing time increases by a factor of 10, while the time between supernovae decreases by a factor of 100, so that $`t_{SN}t_s`$. Supernovae are not distributed randomly over this region but are spatially correlated in regions of star formation, meaning that substantial inhomogeneities will arise between star-forming regions and more quiescent regions. For these reasons it is not surprising to see variations over the distances probed by HD 195965 and HD 191877.

A more precise treatment of mixing in the ISM was carried out by de Avillez & Mac Low (2002). Using numerical simulations, they found that the mixing time scales were a few times 108 years. Chemical inhomogeneities in the ISM are long-lived, so variations in D/H are expected over 1000 pc length scales. In fact, since the mixing time scales increase only slightly with length scale in their models, it is unlikely that the local ISM (e.g., the original seven FUSE sight lines) would have a constant D/H value if any sources of inhomogeneities (i.e., supernovae) have occurred within the past $`10^8`$ years.

We have investigated several alternative explanations for the low D/H values. One possibility is that some of the D column is in high-velocity components which are outside of the velocity integration range used to determine $`N`$(D I) (Jenkins et al., 1999). We have used the O I 1302.169 Å line in the STIS spectrum of HD 195965 to rule this out for that sight line (see § 5). Note that $`50\%`$ of the D I column would have to be missed if the D/H ratio for the HD 191877 sight line were as high as the Moos et al. (2002) mean local ISM value. Figure 7 shows that this is not the case, although we cannot rule out a much smaller systematic offset for HD 191877.

A second possibility is that a significant amount of D is locked up in HD molecules. To test this idea, we have measured the equivalent widths of several strong, isolated Lyman series HD lines in the FUSE spectra: (3-0) R(0) 1066.271 Å, (5-0) R(0) 1042.847 Å, and (4-0) R(0) 1054.288 Å. The column densities derived from the equivalent widths are $`2.4\times 10^{14}`$ cm-2 for the HD 195965 sight line, and $`3.0\times 10^{14}`$ cm-2 for the HD 191877 sight line. For both sight lines this is $`3.5`$% of the total D column density, much lower than the uncertainty in $`N`$(D I).

### 7.2 O/H

The D/O and O/H ratios toward HD 195965 in Table 8 are quite different from the average values in the local ISM found by Moos et al. (2002). This is a surprising result, since Meyer et al. (1998) found very little variation in O/H in their survey, which covered distances and H I column densities similar to those of the sight lines in this paper. Both the O I and H I column densities toward HD 195965 were determined through at least two independent methods, with very similar results. The H I was measured using both IUE and STIS data, and the O I was measured with a single-component COG which combined FUSE and STIS data, as well as independent estimates from the optically thin 1355.598 Å line and the damped 1302.169 Å line in the STIS spectrum. Both of these column densities are therefore quite firmly established, and we believe that it is very unlikely that there is a substantial systematic error in either quantity.

HD 195965 appears to lie on a very high metallicity sight line, with supersolar O/H (Sofia & Meyer, 2001). This may be a region where the signature of stellar processed gas is still visible (as is not the case for the local bubble). Humphreys (1978) found that this star is in the Cygnus OB7 association. If so, it may be in a peculiar region, perhaps in a metal-enriched cloud of ejecta from the star itself or a neighbor. Since gas with high metallicity is thought to have undergone stellar processing, D/H in such gas is expected to be low. This expected anticorrelation between O/H (or N/H) and D/H was not seen in the first seven FUSE sight lines, although it could have been masked by the small number of sight lines studied (Steigman, 2002). Figure 17 compares the O/H and D/H value for HD 195965 with those of the original seven FUSE sight lines. Unfortunately HD 195965 is the only point at high metallicity, so no firm conclusions can be drawn. However, it is possible that this sight line shows evidence of the expected anticorrelation in the ISM (see also Steigman 2002). Note that an error in the H I column density toward HD 195965 would tend to move that point along a diagonal line with a positive slope, and thus could not produce the expected trend of decreasing D/H with rising O/H.

The HD 191877 sight line lacks STIS data, so we could not constrain the O I and H I abundances with multiple estimates. For these reasons we are less confident in the O I column densities found for the HD 191877 sight line. As discussed in § 5, the derived log $`N`$(O I) is very uncertain, although our preferred value, found by coupling the $`b`$-values of N I and O I, is consistent with the FUSE local ISM ratio and the Meyer et al. (1998) ratio. Unfortunately, with the existing information, we cannot make any strong conclusions about O/H toward HD 191877.

### 7.3 N/H

The N/H ratios toward HD 195965 and HD 191877 in Table 8 are higher than the average values in the local ISM found by Moos et al. (2002). Figure 18 shows a plot of N/H versus D/H along both sight lines, compared to the four FUSE sight lines for which N/H was measured. The figure is very suggestive of an anticorrelation between D/H and N/H.

Complicating the interpretation of this plot is the fact N/H ratio is susceptible to ionization effects, especially at low column densities which do not provide adequate shielded from ionizing photons, such as the column densities along the local ISM sight lines (Jenkins et al., 2000). The dashed line in Figure 18 shows the ISM value of N/H found by Meyer et al. (1997). The sight lines used in that study had higher H I column densities than the FUSE local ISM sight lines, so the FUSE sight lines are very likely underestimating N/H, because a significant fraction of N is ionized. The apparent trend in Figure 18 may simply reflect the tendency for $`N`$(N II)/$`N`$(N I) to decrease with increasing $`N`$(H I).

However, an anticorrelation between D and N abundances would be expected if D is destroyed in stars while N is produced (e.g. Steigman et al. 2002). The overabundance of heavy elements accompanied by the underabundance of D seems to fit well within a simple chemical evolution scenario. While ionization effects probably play a role in the high N/H ratio on these sight lines, they may not be the only factor. A comprehensive N/H survey using FUSE spectra is currently underway, and may shed light on this issue.

Several sight lines do not follow the expected trend of low D/H accompanied by high N/H or O/H. Sonneborn et al. (2002) found a low O/H ratio toward BD+284211 with FUSE, while their D/H ratio is close to the Local ISM mean. Jenkins et al. (2000) found a low D/H toward $`\delta `$ Orionis using IMAPS, but also a somewhat low N/H ratio. These deviations suggest that the simple picture of increasing heavy elements abundances with decreasing D/H is not complete.

### 7.4 Future Prospects for D/H Measurements

Measuring D/H values in the Milky Way becomes increasingly difficult as $`N`$(H I) increases, perhaps reaching a fundamental confusion limit. The COGs in Figures 8 and 9 show that the high order D I Lyman transitions are close to saturation. At higher column densities, these transitions will fall on the flat part of the COG and will not give precise column density information. Even higher order transitions at shorter wavelengths, where the D I lines might not yet be saturated at higher column density, are not useful because the H I transitions begin to blend together. These facts make it unlikely that many Galactic disk sight lines more distant than those studied here will be useful for D/H measurements. However, the sight lines analyzed in this paper demonstrate that D/H properties over the distance range from 500 to 2500 pc can be explored with FUSE, even if these sight lines contain a significant amount of H2. This first glimpse at D/H in this distance range suggests that it will be very useful for understanding mixing and chemical evolution in the ISM.

The chemical evolution model of Chiappini et al. (2002) predicts an increase in D/H with galactocentric radius, accompanied by a decrease in N/H and O/H, leading to a pronounced gradient in D/O and N/O. The gas in the outer parts of the galaxy has undergone less stellar processing, and in fact Chiappini et al. (2002) predict that the D/H ratios in the outer disk will approach the primordial value. The sight lines studied here have roughly the same galactocentric radii as the local ISM sight lines of Moos et al. (2002), so we cannot test this model, but with FUSE observations of sight lines in the inner and outer Galactic disk this may prove possible.

We would like to thank Ed Jenkins, Gary Steigman, Jeffrey Linsky, Cristina Oliveira, and Paule Sonnentrucker for enlightening discussions. This work is based on data obtained for the Guaranteed Time Team by the NASA-CNES-CSA FUSE mission operated by the Johns Hopkins University. Financial support to U. S. participants has been provided by NASA contract NAS5-32985. French participants are supported by CNES. This work was partly supported by NASA Guaranteed Time Observer funding to the STIS Science Team under NASA contract NAS5-30403 and is based upon observations obtained with the NASA/ESA Hubble Space Telescope, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. This work has made use of the profile fitting procedure Owens.f developed by M. Lemoine and the FUSE French Team. GH would like to thank Sylvestre Lacour and Arielle Moullet. DK thanks the staff at McDonald Observatory, especially David Doss.