text:
# Transport Properties Calculation for a Quasi-Bidimensional System using T-Matrix Approximation

## 1 Introduction

In recent years one of the significant advances in the research of high-Tc superconductors is the detailed experiments on the pseudogap behavior in underdoped high-Tc superconductors. The pseudogap phenomena have been recognized as one of the most important issues. There are enormous studies from both experimental and theoretical point of view. Influence of the pseudogap formation shows up in different experimental probes, such as NMR relaxation time, Knight shift, neutron scattering, tunneling, photoemission, specific heat, optical conductivity and d.c. resistivity.

The pseudogap phenomena means the suppression of the low frequency spectral weight without any long range order. For cuprate superconductors the low frequency spectral weight begins to be strongly suppressed below some characteristic temperature T much higher than Tc. Typically T is several times larger that Tc in underdoped cuprates and their doping dependence is qualitatively different. While Tc deceases with underdoping, T increases in contrast. This suggest that T does not follow the real Tc line but instead some kind of mean-field critical temperature TMF.

In reality the origin of pairing interaction should be considered to be the antiferromagnetic spin fluctuations. There are studies dealing with pairing correlations obtained by the spin fluctuations on the basis of the fluctuation-exchange (FLEX) approximation . To explain the pseudogap phenomena many theories have been proposed, such as the argument based on BCS-BE cross-over , spin-charge separation induced by stripe order , RVB theory , strong superconducting phase fluctuations , magnetic scenario near the antiferromagnetic instability  and formation of charge density waves (CDW) or spin density waves (SDW) . A feature of the last scenario is that the pseudogap can persist in the superconducting state . The effect of the magnetic field on the pseudogap in high-Tc cuprates has been theoretically investigated by Yanase et.al.  considering the Landau quantization for the superconducting fluctuations as the main effect of the magnetic field. The pseudogap phenomena was studied in the context of non-Fermi Anderson model in  considering the electron-fluctuation interaction. Recent dimensional crossover study by Preosti et.al.  shows that the pseudogap effect is basically absent in three dimensions.

In this paper we describe the pairing scenario based on the strong coupling superconductivity. Generally the strong coupling superconductivity indicates the existence of the incoherent Cooper pairs as the Nozieres and Schmitt-Rink formalism . We think of the pseudogap as the phenomena brought about by the strong superconducting fluctuations. Maly et.al.  have introduced the idea of ”resonances” which means the scattering mechanism by the presence of metastable preformed pairs.

Through several years of extensive experimental work general consensus regarding the superconducting gap symmetry seems to be reach, that the superconducting gap has mainly $`d`$-wave character with possibility of a small mixture of other angular momentum state . A remarkable point of the pseudogap is that its structure in momentum space is the same as the superconducting $`d_{x^2y^2}`$ symmetry with continuous evolution through Tc. This implies that the pseudogap phenomena has close connection to the superconducting fluctuations . On the other hand using Intrinsic Tunneling Spectroscopy, Krasnov et.al.  found that the pseudogap is coexisting with the superconducting gap which imply that the two phenomena are different in nature. The precursor superconductivity scenario can not explain these experimental data and a possible explication might be the formation of CDW or SDW . Another striking characteristic of cuprate superconductors is the anisotropy of their structure. The quantity which most evidently displays the anisotropy of a particular material is the optical conductivity $`\sigma (\omega ,T)`$ and the corresponding d.c. resistivity $`\rho \left(T\right)`$ .

The opening of the pseudogap has drastic effect on the physical properties of the high-Tc cuprates . It is found that associated with the pseudogap the in-plane resistivity deviates from the $`T`$-linear behavior  and the T coefficient of the $`c`$-axis resistivity changes sign, showing semiconducting behavior at low temperature .

This paper is constructed as follow. In section 2 we give the model Hamiltonian and explain the theoretical framework. In section 3 we present the results obtained for the spectral function $`A(𝐤,\omega )`$ and the results for the transport properties studied in this paper. In section 4 we summarize the obtained results and give conclusions.

## 2 Theoretical Framework

### 2.1 Model Hamiltonian

In this section we describe the theoretical framework of this paper. We consider a microscopic model which incorporates both strong electron-fluctuation interaction in the $`CuO`$ planes and a weak interlayer coupling. Within each layer we consider the following two-dimensional model Hamiltonian which has $`d_{x^2y^2}`$-wave superconducting ground state:

$$H_l=\underset{𝐤,\sigma }{}\epsilon _𝐤c_{l𝐤\sigma }^+c_{l𝐤\sigma }+\underset{𝐪,𝐤,𝐤^{}}{}V_{\mathrm{𝐤𝐤}^{}}c_{l𝐤}^+c_{l𝐤+𝐪}c_{l𝐤^{}}^+c_{l𝐤^{}𝐪}$$
(1)
where $`c_{l𝐤\sigma }`$ is the annihilation operator for an electron with momentum $`𝐤`$ and spin $`\sigma `$ in the layer $`l`$. The electron dispersion relation is expressed as:

$$\epsilon _𝐤=2t\left(\mathrm{cos}k_x+\mathrm{cos}k_y\right)4t^{}\mathrm{cos}k_x\mathrm{cos}k_y\mu $$
(2)
where $`t`$ and $`t^{}`$ are the nearest and next-nearest neighbors hopping amplitudes and $`\mu `$ is the chemical potential. From now on we will consider $`t`$ equal to unity and $`t^{}=0.5t`$. The pair interaction is written as:

$$V_{\mathrm{𝐤𝐤}^{}}=Vf_𝐤f_𝐤^{}$$
(3)
where

$$f_𝐤=\mathrm{cos}k_x\mathrm{cos}k_y$$
(4)
is the $`d_{x^2y^2}`$-wave form factor. In the previous relation V is negative. The hopping between layers is included in:

$$H=\underset{l}{}H_lt_c\underset{l}{}\left(c_{li\sigma }^+c_{l+1i\sigma }+c_{l+1i\sigma }^+c_{li\sigma }\right)$$
(5)
where $`t_c`$ is the hopping amplitude between layers and $`c_{li\sigma }`$ is the annihilation operator for an electron within the planar site $`i`$, with spin $`\sigma `$ in the layer $`l`$. The important character of high-Tc cuprates is the momentum dependence of the interlayer hopping matrix element $`t_c\left(𝐤\right)`$. The band calculations  of these materials have shown that the c-axis tunneling matrix element is larger at $`(\pi ,0)`$ and symmetry related points. The functional form is:

$$t_c\left(𝐤\right)=\left(\mathrm{cos}k_x\mathrm{cos}k_y\right)^2$$
(6)
On the other side the coupling between layers has little effect on the pseudogap . The out of plane resistance is determined mainly by the properties of individuals $`CuO`$ layers. The large magnitude of the resistivity anisotropy reflects that the $`c`$-axis mean free path is shorter than the interlayer distance, and the carriers are tightly confined to the $`CuO`$ planes and also is the evidence of the incoherent charge dynamics in the $`c`$-axis direction.

The interlayer part of the Hamiltonian (5) will be considered for the calculation of $`c`$-axis charge dynamics and $`\rho _c\left(T\right)`$.

### 2.2 Self-Consistent T-Matrix Approximation

In this paper we focus on the calculation of the Green function directly on the real-frequency axis  in order to avoid the difficulties of controlling the accuracy of the calculations used in the numerical analytical continuations from the imaginary to real axis by Pade algorithm .

We carry out a self-consistent calculation for the spectral function:

$$A(𝐤,\omega )=\frac{1}{\pi }ImG(𝐤,\omega )$$
(7)
where the retarded Green function $`G(𝐤,\omega )`$ is given by:

$$G(𝐤,\omega )=\frac{1}{\omega \epsilon _𝐤\mathrm{\Sigma }(𝐤,\omega )}$$
(8)
where $`\mathrm{\Sigma }(𝐤,\omega )`$ is the retarded self-energy. The imaginary part of the self-energy can be expressed as:

$$Im\mathrm{\Sigma }(𝐤,\omega )=f_𝐤^2\underset{𝐤^{}}{}𝑑\omega ^{}\left[f\left(\omega ^{}\right)+n\left(\omega +\omega ^{}\right)\right]A(𝐤^{},\omega ^{})ImT(𝐤+𝐤^{},\omega +\omega ^{})$$
(9)
where $`f\left(\omega \right)`$ and $`n\left(\omega \right)`$ represent Fermi-Dirac and Bose-Einstein distributions. The real part of the retarded self-energy can be calculated directly or using Kramers-Kronig relation:

$$Re\mathrm{\Sigma }(𝐤,\omega )=p.v.\frac{d\omega ^{}}{\pi }\frac{Im\mathrm{\Sigma }(𝐤,\omega ^{})}{\omega \omega ^{}}$$
(10)
where $`p.v.`$ represents the principal value of the integral. The Hartree term is neglected in the self-energy. The T-Matrix is given by the following relation:

$$T(𝐤,\omega )=\frac{V^2\mathrm{\Pi }(𝐤,\omega )}{1+V\mathrm{\Pi }(𝐤,\omega )}$$
(11)
where $`\mathrm{\Pi }(𝐤,\omega )`$ can be evaluated directly from the spectral function as:

$$Im\mathrm{\Pi }(𝐪,\omega )=\pi \underset{𝐤}{}𝑑\omega ^{}f_𝐤^2\mathrm{tanh}\frac{\omega }{2T}A(𝐤,\omega ^{})A(𝐪𝐤,\omega \omega ^{})$$
(12)
The chemical potential is determined by fixing the density $`n`$ through the relation:

$$\frac{n}{2}=\underset{𝐤}{}𝑑\omega A(𝐤,\omega )f\left(\omega \right)$$
(13)
We can calculate the hole doping $`\delta `$ using the relation $`\delta =1n`$. Beside relation (13) we must take care of the sum-rule satisfied by the spectral function:

$$1=\underset{𝐤}{}𝑑\omega A(𝐤,\omega )$$
(14)
which must be evaluated at each step in the self-consistent calculation. In the present calculation we divided first Brilouinne zone into $`32\times 32`$ or $`64\times 64`$ lattice and the frequency integration was replaced by a discrete sum of $`1024`$ or $`512`$ points. We did not found differences in the behavior of spectral function using these two discretizations.

Performing the self-consistent calculation given by Eqs. (7)-(14) we obtain the spectral function $`A(𝐤,\omega )`$. We also calculated the density of states:

$$N\left(\omega \right)=\underset{𝐤}{}A(𝐤,\omega )$$
(15)
Having obtained the spectral function we can calculate different physical properties of the cuprates. The real part of the $`ab`$-plane optical conductivity is given by:

$$\sigma _{ab}\left(\nu \right)=\sigma _0𝑑\omega \frac{f\left(\omega \right)f\left(\omega +\nu \right)}{\nu }\underset{𝐤}{}\left[\left(\frac{\epsilon _𝐤}{k_x}\right)^2+\left(\frac{\epsilon _𝐤}{k_y}\right)^2\right]A(𝐤,\omega )A(𝐤,\omega +\nu )$$
(16)
where $`\sigma _0=e^2\pi /2h`$.

The in-plane conductivity $`\sigma _{ab}\left(\omega \right)`$ given by eq. (16) neglects the vertex correction. The in-plane conductivity is coherent in character and shows a Drude peak at low frequencies.

For the calculation of $`c`$-axis conductivity $`\sigma _c\left(\omega \right)`$ we consider the incoherent case. Incoherent conductivity corresponds to diffusive $`c`$-axis transmission and amounts in taking the averages of the spectral function $`N\left(\omega \right)`$ over all momenta.

$$\sigma _c\left(\nu \right)=\sigma _0𝑑\omega \frac{f\left(\omega \right)f\left(\omega +\nu \right)}{\nu }\underset{𝐤}{}t_c^2\left(𝐤\right)N\left(\omega \right)N\left(\omega +\nu \right)$$
(17)
The same relation appear in the calculation of optical conductivity in disordered systems . We calculate in the limit $`\omega 0`$ the in-plane resistivity $`\rho _{ab}=lim_{\omega 0}(1/\sigma _{ab}\left(\omega \right))`$ and $`c`$-axis resistivity $`\rho _c=lim_{\omega 0}(1/\sigma _c\left(\omega \right))`$.

In order to evaluate the thermopower we have to calculate :

$$L_{ij}=𝑑\omega \left(\frac{f\left(\omega \right)}{\omega }\right)\left[\underset{𝐤}{}\left(\frac{\epsilon _𝐤}{k_x}\right)^2A^2(𝐤,\omega )\right]^i\omega ^{j1}$$
(18)
and the thermopower can be calculated using:

$$S=\frac{k_B}{\left|e\right|T}\frac{L_{12}}{L_{11}}$$
(19)
## 3 Results

### 3.1 Spectral Function and Phase Diagram

In this section we present the obtained results for spectral function, using the self-consistent T-Matrix approximation described in previous section. In order to make the results more transparent we present the lowest order calculation results and the self-consistent calculations. We also present the phase diagram obtained for this model. As we can see both in Fig. 1 and Fig. 2 in the lowest order calculation we obtain a pseudogap along the direction $`(0,\pi )(\pi ,\pi )`$ and no pseudogap is obtained along the $`(0,0)(\pi ,\pi )`$ direction. This is in agreement with ARPES measurements . Janko et.al. claims that this anomaly is caused by nearly stable cooper pairs in the excited states . The same behavior is observed in the self-consistent calculation as we can see in Fig. 3 and Fig. 4. Qualitatively we found small differences in behavior between lowest and self-consistent calculation for the spectral function.

We also calculated the critical temperature Tc of the system using Thouless criterion . When $`1+V\mathrm{\Pi }(𝐤,\omega )=0`$ the superconductivity occurs. Thouless criterion is equivalent to BCS theory in the weak coupling limit. In our approximation $`T(𝐤,\omega )`$ can be regarded as the propagator of the Cooper pairs. Thouless criterion correspond to the situation in which $`T(𝐤,\omega )`$ has its pole at $`𝐤=0`$ and $`\omega =0`$. We are interested in the behavior of the system in the normal state, near the superconducting critical point, where the superconducting fluctuations are inhaced.

In Fig. 5 we present the phase diagram of the model where $`T_c`$ and $`T_{MF}`$ are the critical temperature calculate in the self-consistent and lowest-order calculated using Thouless criterion. We also calculated the critical temperature in the first order of perturbation theory and found that in the limit of strong coupling the self-consistent calculation of $`T_C`$ reach the limit of first order calculation of $`T_C`$. $`T_{MF}`$ can be viewed as the BCS critical temperature. The phase diagram corresponds to the underdoped region with the hole concentration $`\delta =0.1`$. Roughly speaking the pseudogap phenomena occurs in the region between Tc and TMF. This region is extremely wide in the strong coupling case. As we can see from the shape of spectral function $`A(𝐤,\omega )`$ the main effect of the electron-fluctuation interaction is the change of shape of $`A(𝐤,\omega )`$ which is seen near $`𝐤=(0,\pi )`$. This suggests that the electron-fluctuation interaction shifts the Fermi surface to the place far from $`𝐤=(0,\pi )`$. Similar with other reported self-consistent treatments of the T-matrix  the pseudogap is not seen clearly in the integrated density of states but only in the $`𝐤`$-dependent spectral function at extreme parameters ranges. In Fig. 6a we present the density of states for different couplings at the temperature $`T=0.35t`$ and in Fig. 6b the temperature dependence of the density of states for fixed coupling strength $`V=3t`$. As we can see from Fig. 6a, increasing the coupling strength, the density of states near $`\omega =0`$ slowly decreases but the central peak still remain. The central peak is due to the integration over the direction $`(0,0)(\pi ,\pi )`$ in the spectral function where no pseudogap appears. For the same temperature the central peak remain at the same frequency $`\omega 0`$. Decreasing the temperature the central peak start to shifts to lower frequencies and a hump start to appear in the density of states. The two humps appearing in Fig. 6a are due to the pseudogap formation in the spectral function near the $`(0,\pi )`$ and symmetry related points and their intensity increase with lowering the temperature.

### 3.2 Optical Conductivity and Resistivity

We have calculated the in-plane and $`c`$-axis optical conductivity using eqs. (16) and (17). We also have calculated the resistivity as function of temperature for both in-plane and $`c`$-axis direction. Our approximation assumes the independent electron propagation in each layer and is justified for small interlayer hopping. The interlayer hopping hopping matrix element $`t_c`$ is weak resulting in an incoherent c-axis transport. Many approaches have been suggested  to describe the c-axis transport properties. The optical conductivity was calculated combining well controlled analytical and numerical methods in the case of one-dimensional Mott-Hubbard insulator at zero temperature in . The $`c`$-axis optical conductivity and d.c. resistivity was calculated within $`tJ`$ model in .

In Fig. 7 we present the obtained results for the in-plane conductivity as function of frequency for two temperatures. The in-plane spectrum of $`\sigma _{ab}\left(\omega \right)`$ is dominated by Drude peak at $`\omega =0`$. The intensity of the peak grows by increasing the temperature. In Fig. 8 we present the $`c`$-axis optical conductivity for three temperatures at the doping level $`\delta =0.1`$. The electronic contribution to $`\sigma _c\left(\omega \right)`$ in the normal state of the underdoped cuprates is nearly $`\omega `$ independent and does nor form a Drude peak at $`\omega =0`$, indicating incoherent nature of the charge dynamics along the $`c`$-axis . It should be emphasize that the pseudogap appears only in $`\sigma _c\left(\omega \right)`$. Actually the development of the pseudogap in $`\sigma _c\left(\omega \right)`$ is correlated with a rapid narrowing of Drude width in $`\sigma _c\left(\omega \right)`$. We observe an anomalous broad peak in $`\sigma _c\left(\omega \right)`$ at finite frequencies which grows in intensity as the temperature is lowered with no discontinuity at the superconducting transition similar to the experimental results . In Fig. 9 we present the temperature dependence of the in-plane resistivity $`\rho _{ab}`$. We found a linear behavior for the in-plane resistivity at high temperature and a deviation from linearity below a characteristic temperature $`T^{}.`$ $`T^{}`$ may be defined by different methods . In this paper we consider the temperature $`T^{}`$ the characteristic temperature where the resistivity present a downturn, with a linear behavior above it. The arrows in Fig. 9 indicate $`T^{}`$. Increasing the doping concentration $`\delta `$, $`T^{}`$ decrease.

In Fig. 10 we present the $`c`$-axis resistivity. We found in this case the same linear behavior of resistivity as function of temperatures, at high temperatures and a upturn at lower temperatures $`\left(d\rho _c/dT<0\right)`$ signaling a semiconductor -like behavior. The semiconducting $`T`$ dependence of the $`c`$-axis resistivity is associated with the deepening of the pseudogap. Therefore the contribution to the conductivity from the quasiparticles near the $`(\pi ,0)`$ is suppressed by the pseudogap. The in-plane conductivity is determined by the contributions from ”cold-spots” and the pseudogap is not important for the in-plane conductivity, but has more effect on the $`c`$-axis conductivity. The $`c`$-axis conductivity is determined from the contributions from the ”hot-spots” where the pseudogap is large. Therefore the $`c`$-axis conductivity is remarkably suppressed by the pseudogap. This is in agreement with the experimental results for cuprates . Our results are qualitatively in agreement with the experimental results obtained for $`\left(R_{1x}Ca_x\right)Ba_2Cu_3O_{6+\delta }`$ where with $`Ca`$ doping the number of holes in $`CuO_2`$ planes increases , where $`R`$ stands for $`Y`$ and $`Lu`$. $`\rho _{ab}`$ is linear in $`T`$ above $`T^{}`$ and decreases more steeply below $`T^{}`$ .

### 3.3 Thermopower

Thermopower has become a simple and ideal probe for the studies of high $`T_C`$ superconductors. The thermopower is one of the several physical quantities which distinctly reveals the unusual normal state properties. The results obtained for the thermopower are calculated using eq. (19). The thermopower of cuprates  has the propriety that is not a monotonic function of temperature. Several authors have tried to explain the temperature behavior of thermopower in term of phonon drag, metallic diffusion, semiconducting model and others , but they cannot explain the linear behavior at high temperature  without considering the cylindrical Fermi surface .

In Fig. 11 we present the temperature dependence of the thermopower. It can be seen that the normal state thermopower decrease with increasing hole doping $`\delta `$ to a lower positive value, at a fixed temperature. The temperature dependence of thermopower is linear with a negative slope at high temperature and present a broad peak, shifted toward the low-temperature as the doping increase . The hump of thermopower in the normal state was reported as a kind of universal behavior of $`HTSC`$, corelated with the opening of the pseudogap in the normal state of a superconductor .

In Fig. 12 we present the doping dependence of the thermopower at the interaction strength $`V=5.5t`$ and at the temperature $`T=2t`$. At a fixed temperature much larger than the critical temperature, $`S`$ continuously decreases from underdoped to overdoped regions, changing its sign near optimal doping. The theoretical work presented in  shows that the doping dependence of the thermopower can be explained by the common band dispersion relation of high-Tc cuprates. Therefore the thermopower can be a criterium for measuring the hole doping in $`CuO`$ layers .

## 4 Conclusions

In this paper we have shown that the d-wave pairing scenario, based on the strong coupling superconductivity well explains some of the physical properties of the layered cuprates above the critical temperature $`T_c`$. We have considered that the strong coupling superconductivity takes place in underdoped cuprates and the pseudogap phenomena are its precursor.

We considered the T-Matrix approximation which is the dominant scattering vertex in the vicinity of the superconducting critical point. In this analysis we didn’t considered the vertex corrections. Anyway our theory appropriately explains the pseudogap phenomena in high-Tc cuprates. We did not take care into account the antiferromagnetic spin fluctuations which are important for explaining other properties of cuprates such as the NMR spin lattice relaxation $`1/T_1`$.

Based on this scenario we studied charge dynamics, in-plane and $`c`$-axis resistivity and thermopower and compared the results with the experimental results. As can be seen from previous section our theory fit well the experimental results from $`YBCO`$ class of layered superconductors. For the calculation of in-plane charge dynamics we obtained a spectrum dominated by a Drude term, while for $`c`$-axis dynamics we found an almost frequency independent behavior. For the $`c`$-axis resistivity we found, due to the opening of the pseudogap in the density of states, a semiconductor-like behavior. For the temperature dependence of thermopower we found a linearly dependence at high temperatures with a negative slope. The negative slope depends weakly on the doping level.

## 5 Acknowledgments

The research was supported by the CNCSIS contract no. 3447/1999.

Figure Captions

Figure 1

Spectral function in the lowest order approximation at $`T=0.6t`$ in the $`(0,\pi )(\pi ,\pi )`$ direction at the doping concentration $`\delta =0.1`$ and with the coupling strength $`V=5t`$

Figure 2

Spectral function in the lowest order approximation at $`T=0.6t`$ in the $`(0,0)(\pi ,\pi )`$ direction at the doping concentration $`\delta =0.1`$ and with the coupling strength $`V=5t`$

Figure 3

Spectral function in the self consistent approximation at $`T=0.6t`$ in the $`(0,\pi )(\pi ,\pi )`$ direction at the doping concentration $`\delta =0.1`$ and with the coupling strength $`V=5t`$

Figure 4

Spectral function in the self consistent approximation at $`T=0.6t`$ in the $`(0,0)(\pi ,\pi )`$ direction at the doping concentration $`\delta =0.1`$ and with the coupling strength $`V=5t`$

Figure 5

Phase diagram of the model at hole concentration $`\delta =0.1`$

Figure 6 a

Density of states for temperature $`T=0.35t`$ and for different coupling strength. The solid line correspond to $`V=t`$, the dotted line to $`V=3t`$ and the dashed line to $`V=6t`$. The hole concentration is $`\delta =0.1`$.

Figure 6 b

Density of states for different temperatures and for coupling strength $`V=3t`$. The solid line correspond to $`T=3.5t`$ and the dashed line to $`T=0.35t`$ at hole concentration is $`\delta =0.1`$.

Figure 7

In-plane optical conductivity at $`T=0.4t`$ (solid line) and $`T=3.4t`$ (dashed line) at the doping concentration $`\delta =0.1`$ and coupling strength $`V=5.5t`$

Figure 8

$`c`$-axis optical conductivity at different temperatures. $`T=0.6t`$ (solid line) $`T=1.0t`$ (dash line) $`T=1.6t`$ (dot line) in the limit of small frequencies at the coupling strength $`V=5.5t`$ and doping level $`\delta =0.1`$

Figure 9

Temperature dependence of the in-plane resistivity at different densities, $`\delta =0.1`$ (solid line) ($`T_C=0.21t`$), $`\delta =0.15`$ ( dashed line) ($`T_C=0.3t`$) and $`\delta =0.2`$ (dashed-dotted line) ($`T_C=0.36t`$), at the same coupling strength $`V=5.5t`$. The arrows indicate the $`T^{}`$ .

Figure 10

Temperature dependence of the $`c`$-axis resistivity at different densities, $`\delta =0.1`$ (solid line) ($`T_C=0.21t`$), $`\delta =0.15`$ ( dashed line) ($`T_C=0.3t`$) and $`\delta =0.2`$ (dotted line) ($`T_C=0.36t`$), at the same coupling strength $`V=5.5t.`$

Figure 11

Temperature dependence of thermopower at different densities, $`\delta =0.1`$ (solid line) ($`T_C=0.21t`$), $`\delta =0.15`$ ( dashed line) ($`T_C=0.3t`$) and $`\delta =0.2`$ (dotted line) ($`T_C=0.36t`$), at the same coupling strength $`V=5.5t.`$

Figure 12

Termopower as function of the doping concentration $`\delta `$ at the coupling strength $`V=5.5t`$ and $`T=2t`$


format: text

text:
# Problem Setup

## Problem Setup

Traditional observations of binary systems are carried out using photometric and spectroscopic methodics. The former gives the time dependence of brightness in a specific band $`I(t)`$ and the latter can give the time dependence of wavelength of some Doppler-shifted line $`\lambda (t)`$. Given ephemeris $`\phi (t)`$ is known, the dependencies $`I(t)`$ and $`\lambda (t)`$ can be converted with the help of Doppler formula to the light curve $`I(\phi )`$ and phase dependency of radial velocity $`V_R(\phi )`$.

During last ten years the observations of binary systems in the form of trailed spectrograms for some emission line $`I(\lambda ,t)`$ or in other terms $`I(V_R,\phi )`$ become widely used. A method of Doppler tomography  is suited to analyze the trailed spectrograms. This method provides obtaining a map of luminosity in the 2D velocity space from the orbital variability of emission lines intensity. The Doppler tomogram is constructed as a conversion of time resolved (i.e. phase-folded) line profiles into a map on ($`V_x`$, $`V_y`$) plane. To convert the distribution $`I(V_R,\phi )`$ to Doppler map $`I(V_x,V_y)`$ we should use the expression for radial velocity as a projection of velocity vector on the line of sight, i.e. $`V_R=V_x\mathrm{cos}(2\pi \phi )+V_y\mathrm{sin}(2\pi \phi )`$ (here we assume that $`V_z0`$, the minus sign before $`V_x`$ is for consistency with the coordinate system), and solve an inverse problem that is described by integral equation (see Appendix A of ):

$$I(V_R,\phi )=I(V_x,V_y)g(V_R+V_x\mathrm{cos}(2\pi \phi )V_y\mathrm{sin}(2\pi \phi ))𝑑V_x𝑑V_y,$$

where $`g(V)`$ is normalized local line profile shape (e.g., a Dirac $`\delta `$-function), and the limits of integration are from $`\mathrm{}`$ to $`+\mathrm{}`$. This inverse problem is ill-posed and special regularization is necessary for its solving (e.g. Maximum Entropy Method , Fourier Filtered Back Projection , Fast Maximum Entropy Method , etc., see also ). As a result we obtain a map of distribution of specific line intensity in velocity space. This map is easier to interpret than original line profiles, moreover, the tomogram can show (or at least gives a hint to) some features of flow structure. In particular, the double-peaked line profiles corresponding to circular motion of the gas (e.g. in accretion disk) become a diffuse ring-shaped region in this map. Resuming, we can say that components of binary system can be resolved in velocity space while they can not be spatially resolved through direct observations, so the Doppler tomography technique is a rather power tool for studying of binary systems.

Unfortunately the reconstruction of spatial distribution of intensity on the basis of Doppler map is an ambiguous problem since points located far from each other may have equal radial velocities and deposit to the same pixel on the Doppler map. So the transformation $`I(V_x,V_y)I(x,y)`$ is impossible without some a priori assumptions on the velocity field.

The situation changes drastically when one uses gas dynamical calculations alongside with Doppler tomography technique. In this case we are not need to cope with the inverse problem since the task is solved directly: $`\rho (x,y)\&T(x,y)I(x,y)`$ and $`I(x,y)\&V_x(x,y)\&V_y(x,y)I(V_x,V_y)I(V_R,\phi )`$. Difficulties can arise when converting the spatial distributions of density and temperature $`\rho (x,y)`$, $`T(x,y)`$ into the distribution of luminosity of specific emission line $`I(x,y)`$. For thick lines the formation of line profile should be described by radiation transfer equations (see, e.g., ), therefore to make our preliminary synthetic Doppler maps we assume that the matter is optically thin. As a first approximation we adopt the line luminosity as $`I=\rho `$ and as $`I=\rho ^2T^{1/2}`$ .

We also should emphasize the complexity of Doppler maps analysis for eclipsing binaries (see, e.g., ). The synthetic Doppler map produces the line emissions from any site of binary and suggests that they are visible at all orbital phases, in other words there are no eclipses and occultations of emission regions. Clearly, eclipsing systems violates that assumption. Usually, when dealing with observations, the eclipsed parts of trailed spectrograms are naturally excluded from input data for construction of Doppler tomograms. But conversion of gas dynamical simulation results into Doppler maps suggests using of full set of data. Thus when analyzing of synthetic Doppler maps for eclipsing binaries we should take in mind that some features can be occulted on some phases even for optically thin case.

## 3D gas dynamical simulations: the model and results

The full description of the used 3D gas dynamical model can be found in . Here we pay attention only to the main features of the model.

Let us consider the semidetached binary system with mass of accretor $`M_1`$, mass of donor-star $`M_2`$, separation $`A`$, and velocity of orbital rotation $`\mathrm{\Omega }`$. These parameters are connected by third Kepler’s low $`A^3\mathrm{\Omega }^2=G(M_1+M_2)`$. To describe the gas flow in this binary system we used the 3D system of Euler equations. The calculations were carried out in the non-inertial Cartesian coordinate system rotating with the binary system. To close the system of equations, we used the equation of state of ideal gas $`P=(\gamma 1)\rho \epsilon `$, where $`\gamma `$ is the ratio of heat capacities. To mimic the system with radiative losses, we accept in the model the value of adiabatic index close to unit: $`\gamma =1.01`$, that corresponds to the case close to the isothermal one .

To obtain numerical solution of the system of equations we used the Roe–Osher TVD scheme of a high approximation order  with Einfeldt modification . The original system of equations was written in a dimensionless form. To do this, the spatial variables were normalized to the distance between the components $`A`$, the time variables were normalized to the reciprocal angular velocity of the system $`\mathrm{\Omega }^1`$, and the density was normalized to its value111Because the system of equations can be scaled to density and pressure, the density scale was chosen simply for the sake of convenience. in the inner Lagrangian point $`L_1`$. We adopted the computational domain as a parallepipedon $`[{}_{}{}^{1}/_{2}^{}A\mathrm{}{}_{}{}^{1}/_{2}^{}A]\times [{}_{}{}^{1}/_{2}^{}A\mathrm{}{}_{}{}^{1}/_{2}^{}A]\times [0\mathrm{}{}_{}{}^{1}/_{4}^{}A]`$ (due to the symmetry of the problem calculations were conducted only in the top half-space). A sphere with a radius of $`{}_{}{}^{1}/_{100}^{}A`$ representing the accretor was cut out of the calculation domain. The boundary conditions were taken as ‘free outflow’ on the accretor star and on the outer edges of computational domain. In gridpoint corresponding to $`L_1`$ we injected the matter with parameters $`\rho =\rho (L_1)`$, $`V_x=c(L_1)`$, $`V_y=V_z=0`$, where $`c(L_1)`$ is a gas speed of sound in $`L_1`$ point. As the initial conditions we used rarefied gas with the following parameters $`\rho _0=10^5\rho (L_1)`$, $`P_0=10^4\rho (L_1)c^2(L_1)/\gamma `$, $`𝑽_0=0`$.

Analysis of considered problem shows that the gas dynamical solution is defined by three dimensionless parameters \[16–18\]: mass ratio $`q=M_2/M_1`$, Lubow-Shu parameter $`ϵ=c(L_1)/A\mathrm{\Omega }`$ , and adiabatic index $`\gamma `$. The value of adiabatic index was discussed above and we used the value $`\gamma =1.01`$. Analysis of our previous results  shows that the main characteristic features of 3D gas dynamical flow structure are qualitatively the same in wide range of parameters $`q`$ and $`ϵ`$. Therefore for model simulation we chose them as follows: $`q=1`$, $`ϵ={}_{}{}^{1}/_{10}^{}`$.

The results of calculations will be presented in coordinate system (see Fig. 1) which is widely used in Doppler mapping. The origin of coordinates is located in the center of the accretor, ‘$`x`$’-axis is directed along the line connecting the centers of stars, from accretor to the mass-losing component, ‘$`z`$’-axis is directed along the axis of rotation, and ‘$`y`$’-axis is determined so that we obtain a right-hand coordinate system (i.e. ‘$`y`$’-axis points in the direction of orbital movement of the donor-star). In Fig. 1 we put digits showing the phase angles of the observer in binary system. We also show a critical Roche lobe with shadowed donor-star and ballistic trajectory of a particle moving from $`L_1`$ point to accretor.

The morphology of gaseous flows in considered binary system can be evaluated from Figs 2, 3 and 4. In Fig. 2 the distribution of density over the equatorial plane and velocity vectors are presented. In this Figure we also put a gas dynamical trajectory of a particle moving from $`L_1`$ to accretor (a white line with circles) and a gas dynamical trajectory passing through the shock wave along the stream edge (a red line with squares, see also Fig. 6). In Fig. 3 we present the so called texture figure in equatorial plane which is visualization of velocity vectors field using Line Integral Convolution Method . In Fig. 4 the similar texture is presented for $`YZ`$ slice along a blue line in Fig. 3. Empty region in Fig. 4 corresponds to cross-section of the stream.

Analysis of presented results as well as our previous studies  show the significant influence of rarefied gas of circumbinary envelope on the flow patterns in semidetached binaries. The gas of circumbinary envelope interacts with the stream of matter and deflects it. This leads, in particular, to the shock-free (tangential) interaction between the stream and the outer edge of forming accretion disc, and, as the consequence, to the absence of ‘hot spot’ in the disc. At the same time it is seen, that the interaction of the gas of circumbinary envelope with the stream results in the formation of an extended shock wave located along the stream edge (‘hot line’). From Figs 2 and 3 it also seen that spiral shock tidally induced by donor-star appears (dotted line in Fig. 2). Appearance of tidally induced two-armed spiral shock was discovered by Matsuda . Here we see only one-armed spiral shock. In the place where the second arm should be we see the strong gas dynamical interaction of circumbinary envelope with the stream causing the formation of rather intensive shock along the outer edge of the stream. This gas dynamically induced shock probably prevents the formation of second arm of tidally induced spiral shock.

An analysis of the flow structure in $`YZ`$ plane (Fig. 4) shows that a part of the circumstellar envelope interacts with the (denser) original gas stream and is deflected away from the orbital plane. This naturally leads to the formation of ‘halo’. Following to , one can define ‘halo’ as that matter which: i) encircles the accretor being gravitationally captured; ii) does not belong to the accretion disc; iii) interacts with the stream (collides with it and/or overflows it); iv) after the interaction either becomes a part of the accretion disc or leaves the system.

Described gas dynamical features of flow structure are good candidates to be observed using Doppler mapping technique.

## Synthetic Doppler maps: methodics

The Doppler maps show the distribution of luminosity in the velocity space. Each point of flow has a three-dimensional vector of velocity $`𝑼=(U_x,U_y,U_z)`$ in observer’s (inertial) frame. In the case when observer is located in the orbital plane of the binary the Doppler map’s coordinate $`(V_x,V_y)`$ will coincide with $`U_x`$ and $`U_y`$. To define these coordinates for the case of inclined system we have to find a projection of vector $`𝑼`$ on the plane constituted by vectors $`𝒏`$ and $`𝒏\times 𝛀`$, where $`𝒏`$ is a direction from the observer to binary.

The line emissivity in the velocity space can be written as:

$$I(V_x,V_y)\underset{𝒪}{}\underset{U_x}{}\underset{U_y}{}I(x,y,z)\delta (U_x(x,y,z)\mathrm{sin}i+U_z(x,y,z)\mathrm{cos}iV_x)$$

$$\delta (U_y(x,y,z)\mathrm{sin}i+U_z(x,y,z)\mathrm{cos}iV_y)d𝒪dU_xdU_y,$$

where $`d𝒪=dxdydz`$, $`i`$ – inclination angle. Usually $`U_zU_x,U_y`$ in the most dense parts of the flow therefore we can neglect the third component of velocity $`U_z`$. This assumption allows us to take off ‘$`\mathrm{sin}i`$’ beyond the integral and to simplify the expression. Using $`V_x/\mathrm{sin}i`$ and $`V_y/\mathrm{sin}i`$ as coordinates for Doppler map will hide the dependency on the inclination angle.

The adopted dimensionless coordinate system for Doppler maps is shown in Fig. 5. We also put in Fig. 5 digits (the same as in Fig. 1) showing the phase angles of observer in binary system. The transformation of the donor-star from spatial to velocity coordinate system is very simple as it is fixed in the corotating frame. Every point $`𝒓`$ fixed in the binary frame has a velocity $`𝛀\times 𝒓`$ in corotation frame. This is linear in the perpendicular distance from the rotation axis and therefore the shape of the donor-star projected on the orbital plane is preserved (see Fig. 5). Since the velocity of each point of donor-star is perpendicular to the radius vector, all points of donor-star are rotated by 90 counter-clockwise between the spatial (see Fig. 1) and velocity coordinate diagrams .

Figure 5 depicts in velocity coordinates a critical Roche lobe with shadowed donor-star and ballistic trajectory of a particle moving from $`L_1`$ to accretor. On the velocity plane the accretor has coordinates (0,$`K_1`$), where $`K_1=A\mathrm{\Omega }M_2/(M_1+M_2)`$ or in dimensionless form (for adopted value of $`q=1`$) $`K_1={}_{}{}^{1}/_{2}^{}`$.

## Synthetic Doppler maps: analysis

As it was mentioned earlier we used both $`I=\rho `$ and $`I=\rho ^2T^{1/2}`$ to build synthetic Doppler maps. The synthetic Doppler map calculated with $`I=\rho `$ is presented in Fig. 6. This map shows the intensity integrated over the $`z`$-coordinate, i.e. taking into account all $`z`$-layers of computational grid. Figure 7 presents the Doppler map corresponding to $`I=\rho ^2T^{1/2}`$. The comparison of results presented in Figs 6 and 7 shows that all characteristic features are the same for different prescriptions of $`I`$, therefore for further analysis we will use results for $`I=\rho `$.

For the sake of comparison we have marked different features of flow structure and put these marks both in spatial and on velocity coordinates. Let us consider results presented in Figs 2 and 6. As it follows from these figures stream from $`L_1`$ (first part of curve marked by white line with circles) transforms to a spiral arm in the 2-nd and 3-d quadrants of the Doppler map. Shock wave caused by interaction of circumbinary envelope with the stream is located along the stream edge in spatial coordinates. Three last points (marked by larger symbols) of curves with circles and squares in Fig. 2 are the examples of two flowlines passing through the shock. The position of this shock on Doppler map is an spiral arm starting approximately from the center of donor-star in the direction of negative $`V_x`$. This arm lies above the arm caused by stream from $`L_1`$.

To analyze other features we have marked the corresponding sites both in spatial and velocity coordinates (see Fig. 8a,b). This comparison is a straightforward problem since we simply choose from all sites of gas dynamical flow structure those that fall in specific places on the Doppler map and mark them. It is seen that tidally induced spiral shock (dashed line in Fig. 2) or more precisely the dense post-shock zone (blue squares in Fig. 8a) produces bright arm in the first quadrant of Doppler map. It is important to note that gas of circumbinary envelope overflowing the stream (cyan diamonds and green up triangles) significantly increase the luminosity of the zone of Doppler map located below the zone corresponding to the stream. Analysis of other features can be conducted by comparison of Figs 8a and 8b as well (magenta circles, yellow stars and red asterisks).

## Doppler maps: observations

Doppler tomography of binaries is widely and actively used now. Since pioneer work by Marsh and Horne  a lot of Doppler map observations were made (see, e.g., \[25–28\]). Just to give to reader the idea on observable Doppler maps we present here in Fig. 9 two Doppler maps of HeII $`\lambda 4686`$Å for IP Peg obtained by Morales-Rueda, Marsh and Billington .

## Conclusions

Three-dimensional gas dynamical simulations of flow structure in semidetached binaries allow to build synthetic Doppler maps. In turn it gives a possibility to identify main features of the flow on the Doppler map without solution of ill-posed inverse problem. Comparison of synthetic tomograms with observations makes possible both to refine the gas dynamical model and to interpret the observational data.

## Acknowledgments

This work was supported by the Russian Foundation for Basic Research (grants 99-02-17619, 00-01-00392, 00-02-17253) and by grants of President of Russia (99-15-96022, 00-15-96722).


text:
# Phase Diagram of the Hubbard Model: Beyond the Dynamical Mean Field

## Introduction

The rich phenomenology of high-$`T_c`$ superconductors has stimulated strong experimental and theoretical interest in the field of strongly correlated electron systems. Common to all high-$`T_c`$ systems is the presence of antiferromagnetic ordering in undoped samples in proximity to a superconducting phase with a $`d`$-wave order parameter and the normal state pseudogap dominating the physics in underdoped samples. A successful theory must describe all these fundamental features at the same time.

The 2D Hubbard model in the intermediate coupling regime or closely related models like the t-J model are believed to capture the essential physics of the high-$`T_c`$ cuprates . The antiferromagnetic phase of the cuprates is well understood. In the strong coupling limit $`UW`$, where $`U`$ is the Coulomb repulsion and $`W`$ the bare bandwidth, the undoped Hubbard model reduces to the Heisenberg model, which has been proven to describe the low energy spin fluctuations of the cuprate parent compounds. However, off half-filling there is no complete understanding of the superconducting phase or the normal state pseudogap in the intermediate coupling 2D Hubbard model.

Finite size quantum Monte Carlo (QMC) calculations for the doped 2D Hubbard model in the intermediate coupling regime with $`UW`$, support the idea of a spin fluctuation driven interaction mediating $`d`$-wave superconductivity . However, the fermion sign problem and the fact that the number of degrees of freedom grows rapidly with the lattice size, limits these calculations to temperatures too high to study a possible transition. These calculations are also restricted to relatively small system sizes making statements for the thermodynamic limit problematic, and inhibiting studies of the low energy physics.

These shortcomings do not apply to the Dynamical Mean Field Approximation (DMFA), which is by construction in the thermodynamic limit. Unfortunately, the lack of non-local dynamics in the DMFA inhibits a possible transition to a state with a non-local ($`d`$-wave) order parameter.

However, it is well known within phenomenological theories that short-ranged antiferromagnetic spin fluctuations mediate pairing with $`d`$-wave symmetry and cause a pseudogap in underdoped samples . There is experimental evidence that the correlation length of dynamical spin fluctuations in the optimally doped cuprates is very small. Therefore, microscopic studies, which account for short-ranged dynamical correlations in addition to the local correlations of the DMFA are relevant and might succeed in describing the physics of the cuprates.

The recently developed Dynamical Cluster Approximation (DCA) is a fully causal approach which systematically incorporates these non-local corrections to the DMFA by mapping the lattice problem onto an embedded periodic cluster of size $`N_c`$. For $`N_c=1`$ the DCA is equivalent to the DMFA and by increasing $`N_c`$ the dynamic correlation length can be gradually increased while the calculation remains in the thermodynamic limit. Previous DCA calculations have indicated the presence of an extended d-wave superconducting phase at moderate doping.

In this manuscript we present calculations of the full phase diagram of the 2D Hubbard model studied with the DCA. To solve the cluster problem we use QMC. We choose $`N_c=4`$, the smallest cluster that includes non-local corrections while preserving the full translational and point group symmetries of the lattice. The results are compared to DMFA calculations, $`N_c=1`$, to illustrate the effect of the initial non-local corrections.

## Formalism

A detailed discussion of the DCA formalism was given in previous publications . The compact part of the free energy is coarse-grained in reciprocal space, projecting it onto a finite-sized cluster of $`N_c`$ points embedded in a self-consistently determined host. The cluster problem is solved using the Hirsch-Fye impurity algorithm modified to simulate an embedded cluster, and the spectra are analytically continued with the maximum entropy method. Once the cluster problem has been solved, lattice susceptibilities may be calculated. All calculations are done in the normal, paramagnetic state. We search for continuous phase transitions indicated by the divergence of the corresponding susceptibilities.

The Hubbard model is characterized by a near-neighbor hopping $`t`$ and an one-site repulsion $`U`$. We choose $`t=1/4`$ to establish a unit of energy and choose $`U=W=2`$ which is sufficiently large that for $`N_c4`$, a Mott gap is present in the half-filled model. The phase diagram of the Hubbard model, calculated with DMFA, displays a range of behaviors including Mott insulating, antiferromagnetic and Fermi liquid regimes. The inclusion of non-local corrections yields significant changes to the phase diagram, including the enhancement of the Mott phase at half filling, the suppression of antiferromagnetism, the introduction of a d-wave superconducting phase and non-Fermi liquid behavior. These different regimes are delineated by calculating the antiferromagnetic (Fig. 1), bulk magnetic (Fig. 2) and d-wave superconducting (Fig. 4) susceptibilities, and the single-particle self energy (Fig. 3).

## Antiferromagnetism

In Fig. 1, the inverse antiferromagnetic susceptibility is plotted versus temperature for dopings $`\delta =0`$ and $`0.05`$ and $`N_c=4`$ and $`1`$. We plot the inverse susceptibility to show that in contrast to finite-sized system calculations, the susceptibility diverges for low doping at $`T=T_N`$. This indicates a transition to an antiferromagnetic phase. In the DMFA for the 2D model, or as we found previously for the infinite-dimensional model, the antiferromagnetism persists to relatively high temperatures and dopings. The non-local dynamical fluctuations, included in the DCA for $`N_c>1`$, strongly suppress the antiferromagnetism. Their effect becomes pronounced for low temperatures and dopings. For example, when $`\delta =0`$, the $`N_c=1`$ and $`N_c=4`$ AF susceptibilities are identical at high temperatures due to the lack of non-local correlations, but separate as the temperature is lowered. The $`N_c=1`$ susceptibility diverges with mean-field exponent of about one; whereas the $`N_c=4`$ result diverges at a much lower temperature with a larger exponent. Consistent with the Mermin Wagner theorem, $`T_N`$ continues to fall for large values of $`N_c`$ (not shown).

## The Pseudogap and Non-Fermi Liquid Behavior

The bulk ($`𝐤=0`$) magnetic susceptibility and single-particle density of states (DOS) display evidence of a pseudogap for low doping $`\delta <0.2`$. We show the bulk magnetic susceptibility in the inset to Fig. 2 for three different dopings. For low to intermediate doping, it develops a peak at low temperatures, defining a temperature $`T^{}`$. $`T^{}T_N(N_c=1)`$, the mean field transition temperature (see Figs. 5 and 1). For $`N_c>1`$, it defines the temperature where short-ranged spin correlations first emerge. The underdoped bulk susceptibility data, $`\delta 0.075`$, may be scaled onto one curve by plotting versus $`T/T^{}`$ (not shown). A similar peak or downturn and scaling is seen in the Knight-shift data of the cuprates. The downturn of the susceptibility is accompanied by a loss of states near the Fermi energy. For temperatures $`T<T^{}`$, a pseudogap begins to develop in the DOS, as shown in Fig. 2. The pseudogap is widest, measured from peak to peak, at low doping, and vanishes for $`\delta 0.2`$. The depth of the pseudogap is greatest when $`\delta 0.05`$, and it vanishes as $`\delta 0`$, where it is replaced by a Mott gap of width $`U`$.

The slow fall of $`T^{}`$ with doping, indicates that the short-ranged spin correlations diminish slowly upon doping. Thus, even in the region $`T^{}>T>T_N`$ the antiferromagnetic correlations still have a significant effect. This is supported by the behavior of the self energy for temperatures $`T<T^{}`$. Once the DCA algorithm is converged, the lattice self energy is calculated by interpolating the cluster result on to the full lattice Brillouin zone. Thus the lattice self energy at any $`𝐤`$ is dominated by the cluster self energy at the nearest cluster momentum. For a Fermi liquid, the self energy $`\mathrm{\Sigma }(𝐤,\omega )(11/Z)\omega ib\omega ^2`$ where $`b>0`$ and $`1/Z>1`$. Our results show that near half filling, the self energy displays non-Fermi liquid behavior. This is illustrated in Fig 3 where we plot the low-frequency self energy at the DCA cluster momenta for $`\delta =0.05`$ and $`T=0.023`$. For momentum points near $`𝐤=(\pi ,0)`$, the imaginary part of the self energy crosses the Fermi energy almost linearly. Concomitant with this behavior is a pseudogap of width $`|J|4t^2/U`$ in the single-particle spectra $`A(𝐤,\omega )`$ for momenta near $`𝐤=(\pi ,0)`$ (not shown).

The pseudogap and the anomalies in the self energy vanish when $`T^{}`$ falls to zero. Here, as shown on the right of Fig. 3, the self energy becomes Fermi-liquid-like with quasiparticle weight $`Z1/2`$. A systematic study of the evolution of the single-particle spectra and the Fermi surface will be presented elsewhere.

It is important to stress that the pseudogap, the downturn of the bulk magnetic susceptibility and the non-Fermi liquid behavior in the self energy are absent when $`N_c=1`$ due to the lack of non-local fluctuations.

## Superconductivity

We searched for many different types of superconductivity, including s, extended-s, p and d-wave, of both odd and even frequency and we looked for pairing at both the zone center and corner. Of these, only the odd-frequency s-wave and even-frequency d-wave pair field susceptibilities at the zone center were strongly enhanced, and only the d-wave susceptibility diverged. This is illustrated in Fig. 4 where the odd-frequency s-wave and even frequency d-wave $`𝐪=0`$ susceptibilities are plotted versus temperature for $`U=2`$ and $`\delta =0.05`$. The s-wave and extended s-wave $`𝐪=0`$ even frequency susceptibilities are also plotted for comparison.

The behavior of the d-wave pair-field susceptibility as a function of temperature for $`N_c=1`$ and $`4`$ and $`\delta =0`$ and $`0.05`$ is shown in the inset to Fig. 4. For $`N_c=1`$ there is no tendency towards pairing. For the DMFA there is no pairing with symmetries lower than the lattice symmetry (i.e., p-, d-wave, etc.).

d-wave pairing is strongly enhanced for $`N_c=4`$ over the corresponding DMFA results. However, for $`\delta =0`$ the inverse susceptibility rises abruptly as the temperature is lowered and the Mott gap opens in the DOS. The Mott gap becomes more pronounced as $`N_c`$ increases, so that for larger clusters the gap prevents superconductivity even for $`U<W`$. If charge excitations are gapped, then pairing is suppressed. At half filling, for $`U=2`$ the gap is of order $`U`$, and is much larger than the magnetic exchange energy $`|J|4t^2/U=0.125`$, so that the opening of the Mott gap will suppress any magnetically mediated pairing. Away from half filling the width of the pseudogap in the charge excitation spectrum is much smaller, on the order of $`J`$, so magnetically mediated pairing is possible. For $`N_c=4`$ and $`\delta =0.05`$, the d-wave pair field susceptibility diverges at $`T_c0.021`$, with an exponent which is less than one, indicating that the fluctuations beyond DMFA which suppress the antiferromagnetism are also responsible for pairing.

The phase diagram of the system is shown in Fig. 5. In addition to the d-wave superconducting and antiferromagnetic phase boundaries, we include $`T^{}`$, the pseudogap temperature fixed by the the peak bulk susceptibility. At low temperatures, it serves as a boundary separating the observed Fermi liquid and non-Fermi liquid behavior. For $`T<T^{}`$ and $`\delta <0.2`$ the self energy shows non-Fermi liquid character for the parts of the Fermi surface closest to $`𝐤=(\pi ,0)`$; whereas, the low-temperature self energy is Fermi-liquid like for $`\delta 0.2`$. The d-wave transition temperature is maximum at $`\delta 0.05`$. The superconductivity persists to large doping, with $`T_c`$ dropping very slowly. In contrast to experimental findings, the pairing instability (preceded by an AF instability) persists down to very low doping. One possible reason for this is that the model remains very compressible down to very low doping $`\delta 0.025`$. This could be due to the lack of long-ranged dynamical spin correlations or stripe formation which could become more relevant as $`N_c`$ increases or when multiple Hubbard planes are coupled together. The effect of such additional non-local corrections ($`N_c>4`$) is presently unknown. However, we believe that a finite mean-field coupling between Hubbard planes will stabilize the character of the phase diagram presented here as $`N_c`$ increases.

## Conclusion

We have used QMC and DCA to study the effect of the initial non-local corrections on the phase diagram of the 2D Hubbard model at intermediate coupling $`U=W`$. The corrections make significant changes, including a strong suppression of the antiferromagnetism, the emergence of non-Fermi liquid (pseudogap) and d-wave superconducting regimes. The critical exponent of the pair-field susceptibility is smaller than one, whereas the antiferromagnetic susceptibility diverges with a critical exponent larger than one. This indicates that the same fluctuations that suppress antiferromagnetism upon doping, mediate pairing. At half-filling the formation of the Mott gap of width $`J`$ suppresses pairing.

## Acknowledgements

It is a pleasure to acknowledge useful discussions with P.G.J. van Dongen, B. Gyorffy, D. Hess, C. Huscroft, J. Keller, H.R. Krishnamurthy, S. Moukouri, Th. Pruschke, J. Zaanen . This work was supported by NSF grants DMR-0073308 and PHY94-07194. Computer support was provided by the Ohio Supercomputer Center.


text:
# The Development of Gas/Star Offsets in Tidal Tails

## 1. Introduction

Since the original work of Toomre & Toomre (1972) and Wright (1972), it has been recognized that the development of tidal tails during galaxy encounters is predominantly a kinematic effect – stars are tidally torn away from the parent galaxy and orbit as “test particles” in the evolving potential. While self-gravity and hydrodynamic forces can play a significant role in the development of substructure within tidal debris (see, e.g., Wallin 1990; Barnes & Hernquist 1992; Elmegreen, Kaufman, & Thomasson 1993), the global evolution of tidal tails is much simpler. Because both gas and stars respond to the gravitational potential in a similar fashion, they should have similar kinematics and spatial morphology on large scales. To be sure, differences exist between the mass distribution of stars and gas in tidal tails – gaseous tails are often much more extended than their stellar counterparts – but this is due largely to the differing radial distribution of stars and gas in the parent galaxy. Because H I disks extend further out than stellar disks (e.g., Broeils & van Woerden 1994; Broeils & Rhee 1997), there is more loosely bound gas to be ejected to large distances in the tidal debris. Nonetheless, the kinematic nature of the process predicts that when material is tidally ejected, stars and gas should remain largely cospatial in the resulting tidal tails.

Observationally, this situation holds true in many cases. Examples of cospatial gas/star tails include nearly all the nearby mergers such as the Antennae, NGC 7252, or NGC 4676 (Hibbard & van Gorkom 1996). This coevolutionary status of tidal gas and stars permits detailed matching between the theoretic predictions of stellar-dynamical models and the observed morphology and kinematics of the gas. Several investigators have taken advantage of this fact to produce dynamical models of specific merging systems (e.g., Stanford & Balcells 1991; Mihos et al. 1993; Hibbard & Mihos 1995; Yun 1999; Barnes & Hibbard 2000).

However, in a significant number of systems the H I and stellar tidal debris are not purely cospatial. In many of these cases, the observed gas/star offsets can be explained, at least qualitatively, by interpenetrating events during galaxy collisions. For example, the offset in the tidal bridge of NGC 7714/5 (Smith et al. 1997) is reasonably expected to occur since tidal bridges form from material on the interface between two interacting galaxies; collisional shocks will act to strip gas material out of the developing tidal bridge. In retrograde encounters, short tails and plumes form from material which passes entirely through the companion galaxy, leading to a strong decoupling of the stars and gas (compare, for example, Figures 11 and 12 of Mihos & Hernquist 1996). In yet other cases, offsets can be attributed to the fact that tidal material has fallen back through the the merging galaxies and may be subject to additional dissipation, as perhaps seen in the offset stellar and gaseous “shells” in Centaurus A (Schiminovich et al. 1994).

However, these explanations cannot account for a small number of cases in which offsets exist between the gas and stars in extended tidal tails. Notable examples include the tidal tails of NGC 520 (Hibbard, Vacca, & Yun 2000), 7714/5 (Smith et al. 1997), and, most dramatically, NGC 3690 (Hibbard & Yun 1999), in which the gas and stars show an offset of $``$ 20 kpc over much of the $``$ 180 kpc extent of the tail. In these examples, none of the previously mentioned mechanisms easily explains the observed offsets. These are all examples of long tidal tails, which develop from the side of the galaxy away from the companion and do not therefore come in physical contact with the companion. Unlike the tidal shells and loops in older mergers, the long tails here have not yet fallen back through their central objects, and, furthermore, the length of the tails precludes a retrograde geometry for the encounter.

At face value, these observations call into question the idea that the evolution of tidal tails is purely kinematic, and suggest that hydrodynamic forces may also play a role in the global evolution of the tidal tails. Hibbard et al. (2000) suggest that the offsets may arise from a variety of effects, including ionization of part of the tail from a starburst wind or from the stars in the tail itself, partial obscuration by dust in the tail, or by an interaction between the tidal gas and a hot ISM component, either a hot gaseous halo or a starburst wind. However, many of these solutions seem somewhat contrived, particularly with respect to the case of NGC 3690. In this case, the gas appears to lead the stellar tail, whereas a gas tail expanding into a hot halo would tend to trail its stellar counterpart. A starburst wind solution is also difficult to envision, as the wind would need to be extremely collimated over very large scales to result in such an extended, linear offset. While ROSAT observations of NGC 3690 (Heckman et al. 1999) do show evidence for a hot outflow, but the X-ray isophotes exist on scales much smaller than needed to explain the offset tidal gas ($``$ 45 kpc compared to the 180 kpc extent of the tidal tail). While the starburst wind may extend to even larger radii at lower temperatures and densities, there is no direct evidence for such an extended outflow.

In this paper we describe a simpler solution for the observed gas/star offsets in some tidal tails. In this model, the offsets arise from the differing radial distributions of gas and stars, coupled with interaction-driven dissipation and inflow in the tidal gas. This scenario explains the observed offsets without resorting to esoteric processes such as hot ISM interactions, starburst ionization, or dust obscuration.

## 2. Interaction Model

To explore the physics which might lead to gas/star segregation in tidal debris, we use numerical modeling of galaxy mergers to follow the kinematic evolution of both the stellar and gaseous components of tidal features. The models are calculated using an N-body treecode coupled with smoothed particle hydrodynamics to follow the evolution of the gas (TREESPH; Hernquist & Katz 1989). The galaxy models consist of a composite disk/bulge/halo system, in which the disk is exponential with mass $`M_d=1`$ and radial scale length $`h=1`$, the bulge is a flattened Hernquist model with mass $`M_b=1/3`$, scale radius $`a=0.4`$ and axial ratio $`b/a=0.5`$, and the halo is a truncated isothermal sphere of mass $`M_h=5.8`$, core radius $`r_c=1.0`$, and tidal radius $`r_t=10.0`$ (see Hernquist 1993 and MH96 for more details). In these units, one half mass rotation period for the disk is 12 units of time, and scaling the model to the Milky Way gives length, mass, and time scale factors of 3.5 kpc, $`5.6\times 10^{10}M_{\mathrm{}}`$, and 13 Myr respectively.

The model we present here is identical to the fiducial model of Mihos & Hernquist (1994, 1996), with one significant difference. In those previous models, the gas was distributed exponentially in the disk, with a radial scale length equal to that of the stellar disk and a total mass of $`M_g=0.1M_d`$. Those models were designed to study inflow and nuclear starbursts, and the gas distribution was chosen to mimic the molecular gas distribution in spiral galaxies. However, the H I in galaxies is typically much more extended than the starlight, so that the gas distribution used in those models was a poor representation of the gas likely to be ejected in tidal debris. To rectify this shortcoming, the models presented here use a different gas mass distribution: the gas is exponential out to a radius of $`2.5h`$, after which it follows a $`r^1`$ surface density distribution out to a truncation radius of $`8h`$. Because of the more slowly declining gas distribution at large radius, Model B has a higher gas mass: $`M_g=0.25M_d`$. Conceptually, Model B represents a more realistic gas distribution; the inner exponential distribution mimics the molecular gas content in galaxies, while the extended HI gas is represented by the slower drop off in density at large radius. The relatively high gas fraction in this model $`M_g=0.25M_d=0.05M_{dyn}`$ is more typical of late type Sb or Sc galaxies than luminous Sa’s (see compilations by Young & Scoville 1991; Roberts & Haynes 1994; McGaugh & de Blok 1997). However, given the large gas content and high $`M_g/L_B`$ ratios of systems like NGC 3690 which show gas/star offsets in their tidal debris, late-type spirals are in fact the likely progenitors of such systems (e.g., Hibbard et al. 2000).

The merger we simulate is again identical to the fiducial merger of MH96 – a merger between two equal mass disk galaxies, placed on an initially parabolic orbit with pericentric distance of $`r_{peri}=2.5h`$. One disk is exactly prograde ($`i=0^{}`$), while the second is highly inclined to the orbital plane ($`i=71^{},\omega =30^{}`$). Figure 1 shows the evolution of the gas and stars in the simulation. As the galaxies first pass one another, the tidal tails are launched and the galaxies move apart, reaching a maximum separation of $`15h`$ before turning around on their orbit and merging together. In terms of the global dynamics (morphology, kinematics, and timescales) this model evolves nearly identically to that of MH96; more details on the global evolution and the triggering of nuclear gas flows can be found in Mihos & Hernquist (1994) and MH96.

However, while the global evolution of this model is similar to that of MH96, on smaller scales striking differences exist between the tidal features of the two models. In MH96, the gas and stars in the two tails are cospatial, while here they exhibit strong differences. In the progade tail the gas is significantly offset from the stars – and in fact leads the stellar tail – while the inclined tail shows a clear bifurcation in the gaseous component. Neither of these features are seen in the MH96 tails (see Figures 3 and 4 of MH96). The gas-leading prograde tail in the current model captures many of the features of the tail in Arp 299, including both the gas/star offset and the fact that the gas and stars become cospatial again at the end of the tail.

Figure 1 alone is enough to weaken models for gas/star offsets which propose interactions between tail gas and a hot ISM or ICM; since no such hot gaseous phase is present in these models, it is not a necessary condition for gas/star segregation. Similarly, starburst outflows are not necessary, nor is in-situ ionization of the gas from the stellar component of the tails. In fact, given the restricted set of physics modeled in the simulation, the offsets must be explainable either by the simple differences in spatial distribution of gas and stars in the progenitor disks, or by the hydrodynamic effects which occur within the gaseous disks. We now examine the development of the tidal tails in detail to separate these two effects.

## 3. Kinematics of Extended Mass Distributions

To decouple hydrodynamic effects from those due to an extended disk of material, a second model was run, identical to the first except that all the gas particles were converted into collisionless stellar particles. The evolution of these collisionless “pseudo-gas” particles is not meant to mimic the evolution of gas, but simply to ask whether an extended mass distribution could result in spatial offsets from a purely kinematic evolution. We leave the question of whether or not it could evolve that way when hydrodynamic effects are included to the following section, which studies the full hydrodynamic evolution of the extended gas disks. For the current exercise, Figure 2 compares the distribution of stars and collisionless pseudo-gas in the galaxies late in the interaction where the offset had developed in the full hydrodynamic model (cf. the final frame of Figure 1).

Although differences can be seen in the distribution of stars and pseudo-gas, the two populations clearly coexist throughout the tails. Pseudo-gas both leads and trails the stellar tail, forming a sheath around the stellar tails. This distribution is similar to that seen in other models of interacting systems, and in fact can be seen in the original collision models of Toomre & Toomre (1972), where rings of particles from subsequentially larger initial radii formed tails which were both longer and broader. This result was further explored in Mihos et al. (1998), who examined the evolution of extended test particle distributions in tidal debris, and showed in more detail how material from large radius forms broader, more extended tidal sheath. The difference between those models and the collisionless model shown in Figure 2 is that in the latter model the pseudo-gas is distributed at all radii, so that the sheath around the stellar tails is filled in by pseudo-gas from the inner regions of the disk. The bottom panels of Figure 2 show this more clearly, comparing the distributions of the inner and outer pseudo-gas, where “inner” and “outer” are defined by gas initially interior or exterior to the radius at which the initial gas distribution changed from an exponential to a $`r^1`$ distribution. Again, the outer pseudo-gas does not lead the stellar tails, but rather is distributed throughout.

While the inclusion of an extended component in the galaxy results in material which leads the stellar tails, it also produces a trailing component – again, the extended pseudo-gas is found in a broader sheath surrounding the stellar component. Comparing such a configuration to observed H I tails is a bit difficult, as there are no hydrodynamical effects included in this collisionless calculation. Certainly one expects that on small scales H I should be more concentrated along the tails, as dissipation acts to keep the internal velocity dispersion of the gas low. On the other hand, on larger scales this calculation argues that a more extended gas distribution could in principle lead to a more diffuse, extended distribution of tidal gas than is suggested by the stellar tidal tails (see, for example, H I imaging of the Antennae by Hibbard & von Gorkom 1996). Any deeper inference must appeal to models which include hydrodynamic effects, which we address in the next section. However, what is clear is that, by itself, an extended spatial distribution of gas does not result in the strong gas/star offset observed in the full hydrodynamic model, or in the tidal tail of Arp 299. Hydrodynamical effects internal to the galaxies must play a role in the formation of these offsets.

## 4. Hydrodynamical Effects

The collisionless model demonstrates that extended mass distributions alone cannot be the entire explanation for offset tidal debris. We now examine the role that hydrodynamical effects play in the evolution of the tidal debris. Naively, one might expect a simple comparison of the models with and without gasdynamics (i.e. the models of Figures 1 and 2) would reveal how gas physics have shaped the evolution of the tail particles. In practice, however, it is not that simple – due to dissipation in the extended and relatively massive gas disks drives further orbital decay, making the hydrodynamic model merge slightly faster than the collisionless one. As a result, the tails in each model evolve slightly differently, making a direct comparison difficult.

To alleviate this problem and study the evolution of the tails in more detail, we introduce the concept of “phase space partners.” We focus on the prograde tail, and isolate gas and star particles within a small piece of the tail where the gas/star offset is most dramatic (Figure 3a). We then examine the model at a time one half rotation period before the initial collision and locate the identified particles in their preencounter position (Figures 3b,c). For each particle, we find the particle from the other component nearest to it in phase space, i.e. for each gas particle $`i`$ we find the star particle $`j`$ with the minimum phase space distance $`\mathrm{\Delta }_{ij}`$ where

$$\mathrm{\Delta }_{ij}^2=\frac{|\stackrel{}{x}_i\stackrel{}{x}_j|^2}{h^2}+\frac{|\stackrel{}{v}_i\stackrel{}{v}_j|^2}{v_c^2}$$

where $`h`$ is the disk scale length and $`v_c`$ is the circular velocity. Because of the large particle number in the model (a total of 65,536 disk particles and 24,576 gas particles in each disk), the phase space separation between the chosen particles and their partners is very small. This ensures that comparing the evolution of the gas particles and their stellar partners are tracking nearly identical phase space origins. By comparing the evolution of the particles and their partners we can thus separate the stellar dynamical evolution of the gas particles from their hydrodynamical evolution.

We begin by showing the evolution of the selected gas and star particles which eventually form the scrutinized portion of the tail (Figure 3d). Before the encounter they are wound into the disk; as the galaxies approach one another the isolated particles unwind until the moment of perigalacticon. At this point, gas and stars are both spatially and kinematically coincident; however, shortly thereafter the gas particles overtake and lead the stellar particles, forming the offset tails. The thickness of the tails reflects their random velocities; the compression of the gas at perigalacticon effectively dissipates the random kinetic energy injected by the collision, so that the gas tail is thinner than the stellar tail.

We next compare the evolution of the gas particles with their stellar phase space partners (Figure 3e). Again, this comparison shows how the gaseous component would have evolved in the absence of hydrodynamic forces. Two things are apparent from this exercise. First, the distribution of gas particles is more compact than that of their stellar partners; again, dissipation has acted to damp out random motion and cause the gas to evolve more cohesively. Second, the gas particles and their stellar partners are essentially cospatial – no offset exists between the components. Hydrodynamic forces act to compress the gas but do not accelerate it; there is no net force applied to the gas (e.g., from other gas particles in the simulation). The gas in this portion of the tidal tail is moving on an essentially collisionless trajectory.

If the tidal gas is acting collisionlessly on large scales, then why is it offset from the collisionless stellar particles? In other words, why is the stellar tail gas poor – what happened to the gas that should have tracked with the stellar tidal tails? Figure 3f shows the evolution of the stars and their gaseous phase space partners. Shortly after collision, the gas partners begin to fall back in towards the inner regions of their host galaxy, unlike the stars which continue to expand outwards. Because of the exponential distribution of stars, the partner gas comes from the inner regions of the disk, where it is more subject to the dissipation which drives radial inflows (e.g., MH96, Barnes & Hernquist 1996). The gas which would have been cospatial with the stellar tidal tails has lost energy and angular momentum and has decoupled from the stellar tidal debris.

Figure 4 directly compares the energy and spin angular momentum content (with respect to the center of the parent galaxy) for the different subsamples of tail particles. Comparing the gaseous and stellar components, we see that, as expected due to its initially larger radii, the gas starts out less tightly bound and with higher angular momentum than the stars. After the galaxies initially collide ($`T30`$), the tidal tails are ejected and both gas and stars end up with a net increase in angular momentum and decrease in binding energy. The path each component takes is markedly different, however. The stars become transiently more bound as they fall into the potential well of the companion, after which they gain energy through the spin-orbit coupling of the interaction and form the tidal tail. In contrast, the gas gains energy monotonically and moves smoothly into the tidal tail. (Recall here that the gas and star subsets sample different regions of initial phase space, so their evolutionary history can be quite different even though they end up near each other in the resulting tails.)

Comparing the energy and angular momentum for the gas and its stellar partners, the values track very tightly together, reflecting the largely collisionless global evolution of the tidal gas. However, we again see significant differences between the stars and their gas partners, in that as the the tails begin to form, the gas partners decouple from the stellar component, losing energy and angular momentum and flowing inwards. Note that the decoupling occurs not at the moment of closest approach – which might be expected if the decoupling was due to shocks from the interpenetrating ISMs – but after the collision, when the galaxy develops a strong self-gravitating response to the collisional perturbation. This result is consistent with the idea that gas loses angular momentum as a result of the internal response of the disk, driven inwards by gravitational torques from the stellar bar and spiral features (see, e.g., Barnes & Hernquist 1991, 1996; Mihos & Hernquist 1996). Once the galaxies finally merge (at $`T65`$) further decoupling occurs as gas is again driven inwards due to the merging process.

In summary, then, the gas/star offset observed in the models arises due to a combination of the different spatial distribution of gas and stars in the disk and the collisional evolution of the gaseous component. However, it is not the hydrodynamics of the gas present in the tails which drives the decoupling; rather, it is the evolution of the gas which would have been in the tails, but was driven inwards by the collision which shapes the offset. In this sense, the formation of the offset is again an internal response of the galaxies involved, rather than due to subsequent evolution of the tidal debris. Certainly the offset is not due to any hydrodynamic interaction between the tail gas and any other gaseous component such as a hot halo or IGM.

## 5. Geometrical Effects

The previous discussion has focussed on the prograde tail, coming from the disk whose rotational plane matched the orbital plane of the encounter. Since for this galaxy the encounter is exactly prograde ($`i=0^{}`$), the response of the disk is maximized, both in the ejection of tidal tails and the onset of dissipational inflows. Studies of developing caustics in tidal tails show their strong dependence on disk inclination (Wallin 1990). In encounters where the perturbation acts entirely in the disk plane, as in a purely prograde encounter, the developing caustics lead to strong shocks and density enhancements in the gas. If the encounter is more inclined, the caustics become more three dimensional, reducing the effects of the caustic-driven shocks. It is natural therefore to ask how sensitive our result is to varied encounter geometries, and we assess this by examining the structure of the tail in the more inclined ($`i=71^{}`$) disk.

Figures 5a-c shows the inclined tail projected onto the disk plane of its host galaxy. Unlike the morphology of the prograde tail, the inclined tail shows gas distributed across the tidal tail, with a spatial spread similar to that of the stellar component (as seen, for example, in Figure 2). Further inspection, however, reveals an interesting twist to tidal tail formation – the inclined tail actually shows a somewhat bifurcated structure. This bifurcation is due in part to material from the prograde companion which has been accreted by the inclined host galaxy and has formed a secondary, slightly trailing tidal tail. Material in this secondary tail is a mixture of original and accreted material, while the leading tail consists solely of material from the original disk. This mixing of tidal debris from the two galaxies into a single tail could in principle lead to peculiar abundance patterns in tidal tails, which could host material with a range of abundances from different progenitors.

Although a significant fraction of gas in the interior portion of the inclined tail has been accreted from the prograde companion, the majority of the gas in the tail comes from the host galaxy. In other words, the lack of an observed gas/star offset is not because the gas tail has been filled in via accretion, but rather because the offset never developed in the inclined tail. An edge-on view of the inclined tail (Figure 5d-e) shows how the tail has warped out of its original disk plane, reducing the effects of collisional shocks at driving inflow along the tail. The lower angular momentum gas which was easily driven inwards in the prograde tail now remains largely in place in the inclined tail. Geometry does thus play a role in the establishment of offsets by altering the efficiency at which the inner gas decouples from the developing tidal tail. This argues that tails which show a marked offset – such as that in Arp 299 – probably arise from encounters in which the host disk was very prograde.

## 6. Discussion

The fact that we can account for the offset between the gas and stars in tidal tails without the need for a separate hot ISM component argues against these offsets necessarily arising from interactions with gaseous halos or starburst winds. We emphasize here that the initial conditions are far from being contrived; on the contrary, the radial distribution of H I in galaxies is typically much more extended than the stellar component (e.g., Broeils & van Woerden 1994; Broeils & Rhee 1997). In truth, the extended gas disk in our model is a more realistic initial galaxy model than those used by MH96 or Barnes & Hernquist 96, where purely exponential gas distributions were used. In essence, better initial conditions have yielded a more accurate description of the formation and evolution of the tidal tails, showing the natural development of the gas/star offset.

These models indicate that the offset is produced when the tail comes from a very prograde encounter; if the encounter is fairly inclined, the dissipation in the gas component is reduced and the offset inhibited. How does this prediction relate to the observed H I properties of mergers? Most nearby mergers do not show strong offsets, consistent with the idea that the geometry must be fairly specific to yield offsets. In two of the more extreme examples of tidal offsets – NGC 3690 and NGC 520 – the orbital geometry inferred from the measured kinematics of the system suggests a prograde encounter (Hibbard et al. 2000), although with significant uncertainties in the exact geometry. If our explanation of the offset holds true, our constraints on the geometry are even more severe than those inferred from the global kinematics of these systems.

Our model also explains the fact that the offset gas tail in NGC 3690 appears to curve back onto the optical tail at large radius (the “hook” described by Hibbard & Yun 1999). Material in the tail at large radius comes from the most loosely bound material at initially large radius in the progenitor disk. Coming from the outskirts of the disk, this material is less susceptible to the inflow driven by the self-gravitating response of the host stellar disk. This material develops in an almost purely collisionless evolution so that the stars and gas remain cospatial. It is only further inwards along the tail that the offset develops, since it is here where the tidal material would have come from deeper in the progenitor galaxy, and be more prone to dissipative effects.

Given the different spatial distributions of gas and stars in the tidal tails, it is also interesting to ask if we see kinematic segregation between the components as well. This has particularly important ramifications for the frequent exercise of modeling specific galaxies based on H I or H$`\alpha `$ velocity fields. Because of the increased computational cost of running hydrodynamic simulations, collisionless N-body models are typically used to match the observed (gas-phase) velocities (e.g., Stanford & Balcells 1991; Mihos et al. 1993; Hibbard & Mihos 1995; Yun 1999; Barnes & Hibbard 2000). If severe kinematic segregation exists in the tidal debris, models derived from matching collisionless kinematics to gas velocities may be in serious error. To examine this possibility, Figure 6 shows a simulated position-velocity plot for both the gas and stars, observed in the plane of the prograde tail such that the portion of prograde tail which shows the offset moves directly towards the observer. This viewing geometry results in the maximum observed kinematic offset between the stars and gas in the prograde tail. These position-velocity plots show that the kinematic offsets are observable, with offset amplitudes of $``$ 50–70km s<sup>-1</sup>, but that they do not radically affect the inferred global kinematics of the system. While detailed model-data matching may hint at differences between the stellar and gas kinematics (see, for example, the discussion in Hibbard & Mihos 1995), the derived models should not be seriously misled by these effects.

Finally, we emphasize that our explanation of gas/star offsets in tidal debris is not the only possible mechanism which may produce offsets. Certainly many examples exist of kinematic decoupling which occurs when the encounter draws tidal material through the inner portions of galaxies, such as retrograde encounters or late infall from tidal tails. The anti-correlation noted by Hibbard et al. (2000) between gas and stars in the tidal debris around Arp 220 are one possible example of this situation, where the stubby tidal features seen in the optical are most likely not due to an extremely prograde encounter. Also, as Hibbard et al. (2000) point out in their discussion of gas/star offsets, many interacting galaxies possess strong starburst winds (Heckman, Armus, & Miley 1987, 1990) and models suggest that collisionally heated halo gas may exist (Barnes & Hernquist 1996); the interaction between the tidal H I and these hot gas components may yet be significant in some cases. What our model shows is that the observed offsets in long tidal tails are easily explained without the need for a hot component, but they do not rule out the presence – or even the effects – of such a hot component.

I thank Lars Hernquist and John Hibbard for many useful discussions. This work is supported through a grant of computing time from the San Diego Supercomputer Center, and by the National Science Foundation through a CAREER Fellowship grant #9876143.