text:
# A density-functional approach to fermionization in the 1D Bose gas

## 1 Introduction

Density-functional theory (DFT) provides a unique framework for the treatment of quantum many-body systems beyond the realm of perturbation theory. While ground-state DFT, based on the theorems of Hohenberg, Kohn, and Sham  deals with the energy and the density profile of inhomogeneous ground states, the time-dependent generalization of these theorems by Runge and Gross  in principle allows the study of time-dependent processes and excited states. Further development of the time-dependent DFT has shown that it is advantageous to consider functionals simultaneously of the density and the current  in order to go beyond the simplistic adiabatic approximation and the microscopic Navier-Stokes equations for the electron fluid have been derived elegantly this way .

In this paper we will apply DFT to the study of fermionization in the system of 1D bosons with contact interaction. In the case of infinite interaction strength, the model is also known as the Tonks-Girardeau gas  and maps exactly to the system of non-interacting spinless fermions with identical energy spectrum and single-particle density. In fact, it has been shown by Cheon and Shigehara  that the exact mapping between the original Bose system and a Fermi system can be extended to arbitrary interaction strength at the expense of a highly singular interaction in the fermionic picture. For the homogeneous Bose gas, exact solutions for stationary states at arbitrary interaction strength can be found using the Bethe ansatz . Although the exact solutions show a continuous transition between the perturbative regime of weakly-interacting bosons and the strongly correlated, fermionized regime, they do not prove very useful for the study of time-dependent processes or inhomogeneous situations. Neither do they provide us with a simple, intuitive picture of how strong the degree of fermionization is in a given system. The crossover from the quasicondensate to the fermionized Tonks gas has also been studied in Refs.

Inspired by the works of Haldane , Sutherland , and others on the concept of exclusion statistics describing a crossover between fermionic and bosonic statistics in 1D systems, we study a model which explicitely allows this transition and provides an intuitive picture of the degree of fermionization. In order to devise a practical scheme for treating time-dependent and inhomogeneous systems we employ DFT and develop a time-dependent Kohn-Sham formalism based on the auxiliary model system of $`N`$ non-interacting spin-$`(\nu 1)/2`$ fermions. This model is chosen because the spin degeneracy may simulate the level attraction or the bunching of single-particle quasi momenta in the interacting bosonic system . The interaction energy of the Bose system is simulated by the kinetic energy of the spin-degenerate fermions. In this work we will study the model in the simplest and most generic approximation: the adiabatic local density approximation (ALDA). The limiting case of infinite interaction strength is obtained easily and is treated exactly with $`\nu =1`$. The opposite limit of weak interaction can also be treated accurately within the proposed formalism with $`\nu =N`$ where the perturbative Gross-Pitaevskii and Bogoliubov equations are recovered asymptotically. In the general case of arbitrary interaction strength, the spin degeneracy $`\nu `$ is fixed by requiring the correct low-energy asymptotics of the excitation spectrum, which is analyzed in the framework of linear-response theory. The resulting model is suited for the study of time-dependent processes in inhomogeneous 1D Bose gases close to the Tonks-gas limit. The properties of this approximate model are analyzed and systematic improvements are suggested.

Before developing the general theory we want to point out that DFT is not commonly used to discuss the theory of ultra-cold bosonic gases although many standard approximations can be derived in this framework. Only few papers therefore explicitely suggested to apply DFT in this context, see e.g. Refs. , where the time-dependent DFT of superfluids  is particularly far developed. The case of strongly interacting bosons in 1D has been considered by Kolomeisky et al.  who suggested a generalization of the Gross-Pitaevskii equation of the following form:

$$i\mathrm{}\frac{}{t}\psi (x,t)=\left\{\frac{\mathrm{}^2}{2m}\frac{^2}{x^2}+v_{\mathrm{ext}}(x,t)+\varphi (n(x,t))\right\}\psi (x,t)$$
(1)
Here, $`\psi (x,t)`$ is a time-dependent complex field and $`n(x,t)=|\psi (x,t)|^2`$ is the one-dimensional density. The function $`\varphi (n)`$ provides the nonlinear term. Kolomeisky et al. suggested to use the chemical potential of the Tonks-Girardeau gas

$$\varphi (n)=\varphi ^{(TG)}(n)\frac{\pi ^2\mathrm{}^2}{2m}n^2$$
(2)
whereas the Gross-Pitaevskii equation is recovered for

$$\varphi (n)=\varphi ^{(GP)}(n)g_{1D}n.$$
(3)
We will refer to Eq. (1) as the bosonic LDA because this equation may be derived as a Kohn-Sham equation in the ALDA using a Bose condensate as the auxiliary non-interacting system for the Tonks gas and the weakly interacting Bose gas, respectively. The same equation also applies for 1D bosons with arbitrary interaction strengths and was used by Öhberg and Santos . In this case, $`\varphi (n)`$ represents the exactly known chemical potential of the Lieb-Liniger model (at $`T=0`$). The bosonic LDA equation (1) is almost identical to a set of hydrodynamic equations derived from general hydrodynamic arguments assuming local equilibrium, the only difference being a kinetic energy pressure term (see e.g. ).

It may be questioned whether the choice of a Bose condensate as the auxiliary non-interacting system of the Kohn-Sham formalism is the ideal starting point for treating a strongly interacting Bose gas in 1D which is not condensed. Indeed it has been shown in the Tonks-Girardeau limit of impenetrable point bosons, that phase coherence is grossly overestimated by Eq. (1) leading to wrong predictions for interference properties . Phase and density are conjugate variables in a quantum field theory and may coexist in the weakly-interacting Bose condensate, which is close to a classical field. It is known already in the perturbative regime, that phase fluctuations are very important in 1D . In the ultimate limit of the Tonks-Girardeau gas, the concept of phase looses its meaning and the governing equations should be rather based on the density alone. For this reason the hydrodynamic approximations employed in Refs. , where the concept of a phase does not appear, are on a safer footing than the bosonic LDA equation of Kolomeisky where phase-coherence is built into the theory.

Moreover, the excitation spectrum and density of the the Tonks-Girardeau gas is identical to the corresponding properties of a gas of independent spinless fermions due to the Bose-Fermi mapping. Already in his original paper, Lieb pointed out, that the fermionic character of the excitation spectrum prevails at finite interaction strength and introduced two elementary branches of excitations which he called type I and type II. Type I excitations result from exciting a particle from the edge of the 1D Fermi sphere to an unoccupied orbital. Thus type I excitations are the particle-hole excitations of the highest possible energy for a given momentum. Type II excitations, on the contrary, are the excitations of the lowest possible energy for a given momentum. They can be seen as taking a particle from the inside of the Fermi sphere and putting it to the first free orbital on the outside.

Only type I excitations are found by standard bosonic methods like Bogoliubov perturbation theory or linear response in Eq. (1) whereas type II excitations and the fermionic character of the spectrum remain hidden. We want to mention at this point that a possible connection between type II excitations and dark solitons in variations of Eq. (1) has been pointed out in the literature .

## 2  An exotic Kohn-Sham scheme

### 2.1 Variational principles

The basic scheme of DFT is the Hohenberg-Kohn variational principle. It allows the ground-state energy $`E_0`$ and one-particle density $`n(x)`$ of a given many-particle system to be found by variationally minimizing the energy functional

$$E_v[n]=F_{\mathrm{HK}}[n]+v_{\mathrm{ext}}(x)n(x)𝑑x.$$
(4)
Here $`v_{\mathrm{ext}}(x)`$ is the external potential and the Hohenberg-Kohn functional $`F_{\mathrm{HK}}[n]`$ is universal in that it does not depend on the external potential. The scheme of Kohn and Sham proceeds by partitioning the unknown functional $`F_{\mathrm{HK}}[n]`$ as

$$F_{\mathrm{HK}}[n]T_s[n]+E_{\mathrm{KS}}[n],$$
(5)
where $`T_s[n]`$ is the kinetic energy of a ficticious, auxiliary system of particles interacting only with a single-particle potential $`v_s`$, which has the same density $`n(x)`$ as the original system. The external potential $`v_{\mathrm{ext}}`$ is contained in $`v_s=v_{\mathrm{ext}}+v_{\mathrm{KS}}`$ where $`v_{\mathrm{KS}}[n]=\delta E_{\mathrm{KS}}/\delta n`$ is a universal functional of the density depending only on the structure and internal interactions of the original and auxiliary many-body problem but not on the external potential. Once an (approximate) expression for $`v_{\mathrm{KS}}[n]`$ has been found, the density of the system in a given external potential $`v_{\mathrm{ext}}`$ is easily found by solving the Kohn-Sham equations, a set of nonlinear independent-particle equations. Kohn and Sham in their original publication used the auxiliary system of noninteracting spin-1/2 fermions while interested in systems of interacting electrons. However, it has been pointed out that the choice of the auxiliary system is somewhat arbitrary and rather a matter of convenience than dictated by laws of nature .

### 2.2 Hamiltonian

We start with the usual Hamiltonian for $`N`$ bosons in 1D with point interactions of Ref. , augmented with a possibly time-dependent external potential $`v_{\mathrm{ext}}(x,t)`$:

$$H=\underset{i=1}{\overset{N}{}}\left\{\frac{\mathrm{}^2}{2m}\frac{}{x_i}+v_{\mathrm{ext}}(x_i,t)\right\}+2g_{1\mathrm{D}}\underset{i<j}{}\delta (x_ix_j),$$
(6)
where $`g_{1\mathrm{D}}=2\mathrm{}^2/(m|a_{1\mathrm{D}}|)`$ is the interaction parameter and $`a_{1\mathrm{D}}`$ is the one-dimensional scattering length.

The exact many-body eigenstates and the complete excitation spectrum of the Hamiltonian (6) in the homogeneous system where $`v_{\mathrm{ext}}=0`$ with periodic boundary conditions have been found by Lieb and Liniger using the Bethe ansatz . The ground state energy is

$$E_0^{\mathrm{hom}}=N\frac{\mathrm{}^2n^2}{2m}e(\gamma )$$
(7)
where $`e(\gamma )`$ is the dimensionless energy per particle in the Lieb-Liniger model and $`\gamma =2/(n|a_{1\mathrm{D}}|)`$ is the single dimensionless parameter of the homogeneous gas with the one-dimensional single-particle density $`n=N/L`$ of $`N`$ particles in a box of length $`L`$. The function $`e(\gamma )`$ is defined as the solution of a Fredholm equation and can be obtained numerically to any desired accuracy. The chemical potential of the homogeneous Lieb-Liniger gas at density $`n`$ is given by

$$\mu _{\mathrm{LL}}(n)=\frac{d}{dN}E_0^{\mathrm{hom}}=\frac{2\mathrm{}^2}{m|a_{1\mathrm{D}}|^2}f(\gamma ),$$
(8)
where $`f(\gamma )=[3e(\gamma )e^{}(\gamma )]/\gamma ^2`$ is a dimensionless function. The asymptotic behaviour of $`e(\gamma )`$ and $`f(\gamma )`$ is known for large and small $`\gamma `$  and, furthermore, the functions are tabulated for intermediate values of $`\gamma `$ in Ref.  111The dimensionless chemical potential $`f(\gamma )`$ is known to have the expansion $`f(\gamma )=\pi ^2\gamma ^2\frac{16\pi ^2}{3}\gamma ^3+𝒪(\gamma ^4)`$ for $`\gamma \mathrm{}`$ and $`f(\gamma )=2\gamma ^1\frac{2}{\pi }\gamma ^{\frac{1}{2}}+𝒪(\gamma ^1)`$ for $`\gamma 0`$. We give the following convenient rational approximation $`f(\gamma )[\pi ^2(48+2\gamma ^2+3\gamma ^3)]/\{3\gamma [8\pi ^2+8\pi \sqrt{\gamma }+\gamma (2+\gamma )^3]\}`$, which deviates with a maximum relative error of $`0.12`$ from the numerical result at $`\gamma 2`$. It has the same expansion as given above for $`\gamma \mathrm{}`$. For small $`\gamma `$ it introduces an error of order $`𝒪(\gamma ^{1/2})`$. .

### 2.3 Functionals and Kohn-Sham equations

Consider a physical system of $`N`$ noninteracting and indistinguishable fermions with spin $`(\nu 1)/2`$. For simplicity let us consider the case where $`\nu `$ divides $`N`$ and assume that all orbitals are nondegenerate. In the ground state $`\stackrel{~}{N}`$ orbitals will be occupied with $`\nu `$ particles each when $`\nu \stackrel{~}{N}=N`$. The kinetic energy is

$$T_s[n]=\frac{\mathrm{}^2}{2m}\nu \underset{i=1}{\overset{\stackrel{~}{N}}{}}\phi _i^{}(x)\frac{^2}{x^2}\phi _i(x)𝑑x.$$

The Kohn-Sham equations are found by functional differentiation with respect to $`\phi _i^{}(x)`$:

$$\epsilon _i\phi _i(x)=\left\{\frac{\mathrm{}^2}{2m}\frac{^2}{x^2}+v_{\mathrm{ext}}(x)+v_{\mathrm{KS}}[n](x)\right\}\phi _i(x).$$
(9)
The potential $`v_{\mathrm{KS}}[n]=\frac{E_{\mathrm{KS}}[n]}{n}`$ is the mean-field potential of Kohn-Sham theory and represents the unknown rest of the functional. In regular electronic DFT, $`v_{\mathrm{KS}}[n]`$ would be written as a sum of the Hartree potential and the exchange-correlation potential. The density is found self-consistently by summing over all occupied ($``$) orbitals

$$n(x,t)=\nu \underset{i}{}|\phi _i(x,t)|^2$$

The ground state density $`n(x)`$ is obviously time-independent and only the $`\stackrel{~}{N}`$ levels with the lowest energies $`\epsilon _i`$ will be occupied. The generalization to the corresponding time-dependent Kohn-Sham equations along the lines of the Runge-Gross theorems  is obvious:

$$i\mathrm{}_t\phi _i(x,t)=\left\{\frac{\mathrm{}^2}{2m}\frac{^2}{x^2}+v_{\mathrm{ext}}(x,t)+v_{\mathrm{KS}}[n](x,t)\right\}\phi _i(x,t).$$
(10)
We note that Eq. (10) is still exact. We now proceed in approximating $`v_{\mathrm{KS}}[n]`$ in a local-density approximation.

To this end, the remaining unknown term $`E_{\mathrm{KS}}[n]`$ is approximated by $`E_{\mathrm{KS}}[n]E_{\mathrm{KS}}^{\mathrm{LDA}}(n)=E_0^{\mathrm{hom}}T_s^{\mathrm{hom}}`$, where $`T_s^{\mathrm{hom}}`$ is the noninteracting kinetic energy of the homogeneous gas. It is simply found by summing the contributions from the $`\stackrel{~}{N}`$ orbitals in the Fermi sphere and multiplying with the degeneracy $`\nu `$:

$`T_s^{\mathrm{hom}}`$ $`=`$ $`\nu {\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \underset{k=(\stackrel{~}{N}1)/2}{\overset{(\stackrel{~}{N}1)/2}{}}}\left({\displaystyle \frac{2\pi k}{L}}\right)^2`$ (11)
$`=`$ $`N{\displaystyle \frac{\mathrm{}^2\pi ^2n^2}{6m\nu ^2}}\left(1{\displaystyle \frac{1}{\stackrel{~}{N}^2}}\right)`$
where $`L`$ is the length of a box with periodic boundary conditions, $`n=N/L=n`$ is the line density and $`\stackrel{~}{N}`$ has been assumed odd.

Before taking the derivative $`\frac{F_{\mathrm{HK}}[n]}{n}`$ in order to find an expression for $`v_{\mathrm{KS}}`$ we have to consider what to do with the term $`1/\stackrel{~}{N}`$. If $`\stackrel{~}{N}=1`$ the contribution from the kinetic energy vanishes and we arrive at the bosonic LDA of Eq. (1). However, if $`\stackrel{~}{N}`$ is of the order of $`N`$, we may disregard the term $`1/\stackrel{~}{N}`$ in the limit of large $`N`$. In fact it is not consistent to keep this term if we approximate the Lieb-Liniger energy with the bulk value as this is only valid in the limit $`N\mathrm{}`$. In the following we will therfore drop this term and arrive at

$`v_{\mathrm{KS}}^{\mathrm{LDA}}(x,t)`$ $`=`$ $`{\displaystyle \frac{d}{dN}}\left(E_0^{\mathrm{hom}}E_{\mathrm{kin}}\right)|_{n=n(x,t)}`$ (12)
$`=`$ $`\mu _{\mathrm{LL}}(n){\displaystyle \frac{\mathrm{}^2\pi ^2n^2}{2m\nu ^2}},`$
where the chemical potential $`\mu _{\mathrm{LL}}(n)`$ of the homogeneous (Lieb-Liniger) gas at density $`n`$ is given by Eq. (8). The remaining parameter $`\nu `$ will be determined by a consistency condition on the excitation spectrum that will be derived in the following paragraph.

### 2.4 Excitation spectrum from linear-response theory

Lieb pointed out that the excitation spectrum of the 1D Bose gas resembles that of the 1D Fermi gas and likewise has a particle-hole type of structure — at any value of the interaction strength ! The spectrum of the homogeneous gas may be equally well described by two branches of elementary excitations $`ϵ^\mathrm{I}(k)`$ and $`ϵ^{\mathrm{II}}(k)`$. Both branches have the same slope for $`k0`$ and therefore lead to the same speed of sound $`v_s=dϵ^\mathrm{I}/dk|_{k=0}=dϵ^{\mathrm{II}}/dk|_{k=0}`$. Bogoliubov perturbation theory yields an approximation only for the branch $`ϵ^\mathrm{I}(k)`$. In the time-dependent Kohn-Sham model introduced above, a particle-hole type excitation spectrum arises naturally. The elementary excitations may be found by solving for the resonant frequencies of the time-dependent equation (10) in the linear-response limit. The resulting excitation energies are given by differences of the orbital energies of the stationary Kohn-Sham orbitals modified by interaction contributions from the nonlinear mean-field $`v_{\mathrm{KS}}^{\mathrm{LDA}}`$.

The derivation of the linear response equations is the same as for the LDA in the much studied electronic systems and very similar to the time-dependent Hartree-Fock or random-phase approximation (see e.g. ) and is summarized in the appendix. Here we want to note that the linear response correction to the spectrum of the Kohn-Sham system is proportional to the derivative $`dv_{\mathrm{KS}}/dn`$. From the explicit linear response equations for the homogeneous system it can be seen that type I and type II excitations have a different slope at $`k=0`$, except for the situation where

$$v_{\mathrm{KS}}^{}(n)=0.$$
(13)
This argument is detailed in the Appendix where the linear response equations in the small momentum limit are solved exactly. Using the expression (12) for $`v_{\mathrm{KS}}`$ we easily find the correct value for the number of modes in the fermionic Kohn-Sham equations (10):

$$\stackrel{~}{N}=N\sqrt{\frac{\gamma ^3}{2\pi ^2}f^{}(\gamma )}.$$
(14)
The value of $`\stackrel{~}{N}/N1/\nu `$ is a measure of the degree of fermionization in the system. When $`\stackrel{~}{N}/N=1`$ the system is completely fermionized whereas $`\stackrel{~}{N}/N0`$ indicates a purely bosonic system with a Bose condensate as the ground state. The functional dependence of $`\stackrel{~}{N}/N`$ on the interaction strength $`\gamma `$ is shown in Figure 1.

Qualitatively the picture agrees with the results of Girardeau and Wright  who modelled the BEC-Tonks crossover by a mixture of a pair-correlated Bose liquid with a non-interacting Fermi gas in a variational framework. We would like to point out, that $`\stackrel{~}{N}`$ is a global parameter whereas $`\gamma `$ is local and varies in an inhomogeneous system. Therefore the condition (14) should be used with care and we expect the present model to be most useful for nearly homogeneous systems.

Although it would be desirable to find an approximation for the condensate fraction from our model, it should be noted that Kohn-Sham DFT as used here generally does not model the many-body wavefunction but gives direct access only to energies and the diagonal part of the single-particle density matrix. The information about off-diagonal long range order and condensate fraction, however, is contained in the inaccessible off-diagonal part of the single-particle density matrix. Nevertheless it should be pointed out that the fractional occupation of the Kohn-Sham orbitals $`\stackrel{~}{N}^1`$ in our model, which is reminiscent of a condensate fraction, vanishes in the thermodynamic limit where $`N\mathrm{}`$ at $`\gamma =const.`$ consistent with the general result that there is no Bose condensate in the $`T=0`$ homogeneous interacting Bose gas .

An important question related to the condensate fraction is local phase coherence . From the construction of the Kohn-Sham model it appears difficult to define expectation values for the phase operator in this scheme. However, the time-dependent DFT does contain information about local phase coherence to the degree that is necessary to correctly describe the supression of interference patterns in the dynamics.

The approximation scheme defined by Eqs. (10), (12), and (14) is the main subject of this paper and shall be referred to as fermionic LDA. In the remainder of this paper we will analyze the predictions of this model and compare with the exact solutions of the Lieb-Liniger model and other approximations.

## 3 Properties

### 3.1 Excitations of the homogeneous gas

It is interesting to calculate the dispersion for the elementary excitations of type I and II, which can be done analytically since the linear-response corrections vanish due to Eq. (13). In the homogeneous gas, the mean field $`v_{\mathrm{KS}}`$ is constant and the Kohn-Sham single-particle energies become

$$ϵ_p=\frac{2\pi ^2\mathrm{}^2}{mL^2}p^2+v_{\mathrm{KS}},$$

where $`p`$ is any integer and $`\mathrm{}k_p=\mathrm{}p\mathrm{\hspace{0.17em}2}\pi /L`$ is the momentum. In the ground state all orbitals with $`p<F`$ are occupied where $`F=(\stackrel{~}{N}1)/2`$ is the index of the Fermi momentum $`\mathrm{}k_F`$. The energies of particle-hole excitations of momentum $`\mathrm{}k_q`$ are given by the difference

$`ϵ_pϵ_{pq}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}n^2{\displaystyle \frac{4\pi ^2}{N^2}}q(2pq).`$ (15)
Figure 2 illustrates the discrete excitation spectrum obtained for a finite number of particles in a box with periodic boundary conditions. Type I and II excitations are the upper and lower bounds of the elementary excitation spectrum, respectively. The fermionic DFT is seen to slightly overestimate the energies of type I excitations and underestimate the energies of type II excitations at large momentum while the correct asymptotes are obtained for small momenta.

From Eq. (15) we can construct explicit expressions for type I and type II excitations. Type I excitations are defined by exciting a particle with the Fermi momentum $`\mathrm{}k_F`$ into an unoccupied orbital with $`p=F+q`$ whereas type II excitations take a particle from inside the Fermi sphere to the lowest unoccupied orbital with $`p=F+1`$. In terms of the dimensionless momentum

$$\stackrel{~}{p}=\frac{\mathrm{}k_q}{\mathrm{}n}=\frac{2\pi q}{N}$$

of the excitations we obtain

$$ϵ^\mathrm{I}(\stackrel{~}{p})=\frac{\mathrm{}^2}{2m}n^2\stackrel{~}{p}\left[2\pi \left(\frac{1}{\nu }\frac{1}{N}\right)+\stackrel{~}{p}\right].$$
(16)
Correspondingly we find for type II

$$ϵ^{\mathrm{II}}(\stackrel{~}{p})=\frac{\mathrm{}^2}{2m}n^2\stackrel{~}{p}\left[2\pi \left(\frac{1}{\nu }+\frac{1}{N}\right)\stackrel{~}{p}\right].$$
(17)
Expressions for the elementary excitation spectrum of the interacting gas can be derived using equation (14). In the thermodynamic limit we obtain

$`ϵ^\mathrm{I}(\stackrel{~}{p})`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}n^2\stackrel{~}{p}\left[\sqrt{2\gamma ^3f^{}(\gamma )}+\stackrel{~}{p}\right],`$ (18)
$`ϵ^{\mathrm{II}}(\stackrel{~}{p})`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}n^2\stackrel{~}{p}\left[\sqrt{2\gamma ^3f^{}(\gamma )}\stackrel{~}{p}\right].`$ (19)
These dispersion relations ought to be compared with the exact ones and with earlier approximations. First of all we note that both dispersions relations indeed lead to the same speed of sound

$$\frac{dϵ^{\mathrm{I}/\mathrm{II}}}{dp}|_{p=0}=c_s\frac{2\mathrm{}}{m|a_{1\mathrm{D}}|}\sqrt{\frac{\gamma }{2}f^{}(\gamma )},$$
(20)
which is identical to the exact result of the Lieb-Liniger model.

The type II excitations, as discussed by Lieb, have the character of elementary excitations for $`\stackrel{~}{p}\pi `$ and at this point the slope $`dϵ^{\mathrm{II}}/d\stackrel{~}{p}`$ vanishes. This condition is strictly fulfilled by Eq. (19) only in the Tonks-gas limit $`\gamma \mathrm{}`$ where the equations become exact. In the general case, $`ϵ^{\mathrm{II}}(\stackrel{~}{p})`$ of Eq. (19) has a maximum at $`\stackrel{~}{p}_m<\pi `$ and $`\stackrel{~}{p}_m`$ approaches 0 as $`\gamma 0`$. This observation is related to the fact that the umklapp transition of taking one particle from one side of the Fermi-sphere to the other is not described correctly in the considered fermionic Kohn-Sham model. Therefore our approach is limited to describing type II excitations in the long-wavelength limit and works best for large $`\gamma `$.

### 3.2 Comparison with bosonic LDA

The bosonic LDA of Eq. (1) is simply the single-mode case with $`\stackrel{~}{N}=1`$ of the general fermionic Kohn-Sham formalism in this paper. The linear-response equations (31) for this case coincide with the Bogoliubov-de Gennes equations derived in Refs. . It is instructive to derive explicit expressions for the excitations of a homogeneous gas, which appear not to be available in the literature. Eqs. (31) now become

$`(ϵ_qϵ_0+{\displaystyle \frac{N\varphi ^{}}{L}})u_q+{\displaystyle \frac{\nu \varphi ^{}}{L}}v_q`$ $`=`$ $`\mathrm{}\omega u_q`$ (21)
$`(ϵ_qϵ_0+{\displaystyle \frac{N\varphi ^{}}{L}})v_q+{\displaystyle \frac{\nu \varphi ^{}}{L}}u_q`$ $`=`$ $`\mathrm{}\omega v_q`$ (22)
where $`\varphi ^{}=d\varphi /dn`$ and we have set $`X_{hp}=u_{ph}`$ and $`Y_{hp}=v_{hp}`$. With

$$ϵ_qϵ_0=\frac{2\pi ^2\mathrm{}^2q^2}{mL^2}$$

we find the simple solution

$$\mathrm{}\omega =\pm \sqrt{\frac{\mathrm{}^2k_i^2}{2m}(\frac{\mathrm{}^2k_i^2}{2m}+2n\varphi ^{})}.$$
(23)
Only the plus sign contributes here, the minus sign is a well understood artifact of linear response theory. In terms of the dimensionless momentum $`\stackrel{~}{p}=k/n`$ we find the following final result for the excitation spectrum of the bosonic LDA with $`\varphi (n)=\mu _{\mathrm{LL}}(n)`$:

$$\mathrm{}\omega =\frac{\mathrm{}^2}{2m}n^2\sqrt{\stackrel{~}{p}^2(\stackrel{~}{p}^22\gamma ^3f^{}(\gamma ))}.$$
(24)
In contrast to the excitation spectrum described by Lieb, there is only one excitation branch in the linear response of the bosonic LDA. Comparing with Eqs. (21) and (22) we see that the first terms in the expansions around $`k=0`$ and $`k=\mathrm{}`$ are identical but there are discrepancies in between. The speed of sound is again the exact result $`c_s`$ as a consequence of the compressibility sum rule which is obeyed by the LDA.

We can recover the well-known Bogoliubov approximation in the limit of small $`\gamma `$ by expanding $`f(\gamma )=2\gamma ^1+𝒪(\gamma ^{1/2})`$ to find the usual dispersion relation

$$(\mathrm{}\omega )_{\mathrm{Bog}}=\frac{\mathrm{}^2}{2m}n^2\sqrt{\stackrel{~}{p}^2(\stackrel{~}{p}^2+4\gamma )}.$$
(25)
The speed of sound becomes

$$v_s^{\mathrm{Bog}}=\frac{\mathrm{}}{m}n\sqrt{\gamma }=\frac{\mathrm{}}{m}\sqrt{\frac{2n}{|a_{1\mathrm{D}}|}}.$$

Figure 3 shows the dispersion relations derived in this paragraph for a particular value of the interaction parameter $`\gamma `$. It can be seen clearly that the Bogoliubov approximation for the speed of sound is wrong, which is no surprise as $`\gamma `$ is not small in this example. However, it is also interesting to see that the dispersion derived from the bosonic LDA deviates significantly from the fermionic LDA.

### 3.3 Limits of strong and weak interaction

In the strongly-interacting limit we find $`\nu =1`$ and $`\stackrel{~}{N}=N`$. In this case $`v_{\mathrm{KS}}=0`$ and the Kohn-Sham equations become the equations for non-interacting spinless fermions. Due to the Bose-Fermi mapping theorem of Girardeau this gives the correct description of both the full excitation spectrum as well as the (time-dependent) diagonal single-particle densities. In particular we obtain the following equations for type I and II excitations:

$`ϵ^\mathrm{I}(\stackrel{~}{p})`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}n^2\stackrel{~}{p}\left[2\pi +\stackrel{~}{p}\right]`$ (26)
$`ϵ^{\mathrm{II}}(\stackrel{~}{p})`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}n^2\stackrel{~}{p}\left[2\pi \stackrel{~}{p}\right]`$ (27)
It is interesting to note that the single-mode equation (1), which here is exactly the equation studied by Kolomeisky , does not have the correct limit of Eq. (26). Instead we find

$$\mathrm{}\omega =\frac{\mathrm{}^2}{2m}n^2\sqrt{\stackrel{~}{p}^2(\stackrel{~}{p}^22\pi ^2)}.$$
(28)
This limit underestimates the type I excitation energy. It has been suggested in Ref.  to link dark solitary waves in the single-mode equation to the type II excitations. The dispersion relation, however, gives only qualitative agreement and underestimates the energy at large momenta. We want to stress that our fermionic LDA, on the contrary, yields the exact excitation spectrum in the Tonks-Girardeau limit.

In the limit of weak interactions we should obtain the single-mode equation with

$$\varphi (n)=\mu _{\mathrm{LL}}(n),$$
(29)
as was discussed in the context of the derivation of the Kohn-Sham equations. It was also discussed that the LDA approach for finite interactions followed in the main part of the paper cannot describe the limit of the Bose condensate with a single occupied mode correctly. Indeed it can be seen that Eq. (18) does not reduce to the Bogoliubov form (25) for small $`\gamma `$, although it has the correct low- and high-energy asymptotics. In order to obtain this limit correctly, a finite size correction should be added to the density functional.

### 3.4 Static density profiles in the Thomas-Fermi approximation

A consistency check for the static density profiles found from the fermionic LDA can be derived from applying a Thomas-Fermi approximation. The idea is to treat the Kohn-Sham model system of Eq. (9) as non-interacting Fermions (equivalent to a Tonks gas) in the external potential $`v_{\mathrm{ext}}(x)+v_{\mathrm{KS}}^{\mathrm{LDA}}(x)`$. Approximating now the density by the Thomas-Fermi approximation or hydrodynamic LDA as in Ref.  leads to exactly the same equations for the density as the Thomas-Fermi approximation for the original system of interacting bosons discussed in Ref. . This result is independent of the specific choice of $`\stackrel{~}{N}`$ and provides a check for the consistency of the LDA between the different model systems.

## 4 Conclusions

In this paper we have devised a variational scheme based on Kohn-Sham DFT to calculate ground-state properties and time-dependent processes in the system of 1D bosons with short-range interaction. The scheme will be particularly useful for the regime of strong interaction where fermionization is important and bosonic perturbation theory ceases to be applicable. In contrast to the previously suggested nonlinear Schrödinger equation of Kolomeisky , our theory has the correct strong-interaction limit for the density and for ground-state and excitation energies. The degree of fermionization in the model is determined by requiring consistency of the low-energy, low-momentum excitation spectrum. Employing the local-density approximation yields a parameter-free model. Interesting applications and test cases may arise in situation where coherence properties and long-wavelength excitations are important. Due to the construction of our model, it is best suited for situations with a large number of bosons and almost homogeneous densities. A specific example may be the study of shock waves in the homogeneous gas as studied by Damski . In the decay of a shock wave, the usual bosonic LDA fails due to the fact that coherence is overestimated. Exact results are only available for the Tonks-Girardeau limit so far. For finite interaction strengths, our fermionic LDA can make predictions about the importance of coherences and numerical studies of this problem are under way. Another situation where the current theory can be applied is the 1D Bose gas in a weak optical lattice. Interesting questions to study are the nature of excitations and the phase diagram as well as a comparison to the predictions of bosonic LDA.

The DFT scheme studied in this paper is based on the careful choice of the kinetic energy functional determined by the Kohn-Sham auxiliary independent-particle system and the rest of the functional was approximated in the ALDA. The next logical step to improve performance of our model is to go beyond the ALDA and employ standard correction schemes based on current-density functional theory .

I would like to thank William Reinhardt, Kiril Tsemekhman, and Subhasis Sinha for encouragement and useful discussions. Support from the NSF, the Alexander von Humboldt Foundation, and ECT (Trento) is gratefully acknowledged.

## Appendix A Linear response equations

Following the standard scheme of linear response theory  we introduce a time-dependent perturbation

$$v_{\mathrm{ext}}(x,t)=v_{\mathrm{ext}}(x)+\lambda [f^+(x)e^{i\omega t}+f^{}(x)e^{i\omega t}]$$

in Eq. (10) where $`\lambda `$ is small. We look for solutions of Eq. (10) through first order in $`\lambda `$. To zeroth order we obtain the stationary equation (9). Up to first order we expect to find solutions of the form

$$\phi _h(x,t)=\left\{\phi _h(x)+\lambda \underset{p}{}(X_{hp}e^{i\omega t}+Y_{hp}^{}e^{i\omega t})\phi _p(x)\right\}e^{\frac{iϵ_ht}{\mathrm{}}}$$
(30)
Here, the linear response of the occupied orbitals can be expanded in terms of unoccupied ($``$) orbitals because the time-evolution according to the Kohn-Sham equations preserves the orthonormality of the occupied orbitals to every order in $`\lambda `$. Substituting the expression (30) into the time-dependent equation (10) and equating the first order terms with different time-dependence separately and finally setting $`f^\pm =0`$ yields the following equations for small-amplitude free oscillations around the stationary state:

$`{\displaystyle \underset{h^{},p^{}}{}}(_{hph^{}p^{}}X_{h^{}p^{}}+_{hph^{}p^{}}Y_{h^{}p^{}})`$ $`=`$ $`\mathrm{}\omega X_{hp}`$ (31)
$`{\displaystyle \underset{h^{},p^{}}{}}(_{hph^{}p^{}}^{}Y_{h^{}p^{}}+_{hph^{}p^{}}^{}X_{h^{}p^{}})`$ $`=`$ $`\mathrm{}\omega Y_{hp}`$ (32)
In the basis of stationary Kohn-Sham orbitals \[i.e. solutions of Eq. (9)\] the matrices read

$`_{hph^{}p^{}}`$ $`=`$ $`(ϵ_pϵ_h)\delta _{pp^{}}\delta _{hh^{}}+\nu v_{ph^{}p^{}h}^{}`$ (33)
$`_{hph^{}p^{}}`$ $`=`$ $`\nu v_{pp^{}h^{}h}^{},`$ (34)
where

$$v_{ijkl}^{}=\phi _i^{}(x)\phi _j^{}(x)\frac{dv_{\mathrm{KS}}}{dn}|_{n(x)}\phi _k(x)\phi _l(x)dx.$$
(35)
A simple basis transformation may be applied to rewrite Eqs. (33) to (35) in a different basis, e.g. the coordinate representation. For the specific case of the bosonic LDA with $`\stackrel{~}{N}=1`$ we find the Bogoliubov–de Gennes equations of Ref. . In the limit $`\gamma 1`$ the regular Bogoliubov–de Gennes equations of Gross-Pitaevskii theory are recovered.

For studying the homogeneous system where $`v_{\mathrm{ext}}=0`$ and $`n`$ is constant we introduce the usual box quantization with $`\varphi _m(x)=e^{ik_qx}/\sqrt{L}`$ and $`k_q=2\pi q/L`$. We find

$$ϵ_i=\frac{\mathrm{}^2k_i^2}{2m}+v_{\mathrm{KS}}$$

and $`v_{ijkl}=\delta _{ik,jl}v_{\mathrm{KS}}^{}/L`$ with $`v_{\mathrm{KS}}^{}dv_{\mathrm{KS}}/dn`$. Due to translational symmetry the matrix equations block and it makes sense to label the blocks according to the momentum transfer $`\mathrm{}k_p\mathrm{}k_h=2\pi q\mathrm{}/L`$. With defining $`_{p,p^{}}^q=_{(pq)p,(p^{}q)p^{}}`$ we obtain

$`_{p,p^{}}^q`$ $`=`$ $`(ϵ_pϵ_{pq})\delta _{p,p^{}}+{\displaystyle \frac{\nu v_{\mathrm{KS}}^{}}{L}}`$ (36)
$`_{p,p^{}}^q`$ $`=`$ $`{\displaystyle \frac{\nu v_{\mathrm{KS}}^{}}{L}}`$ (37)
Due to the presence of the Fermi sphere there are restrictions on the possible values of the indices as $`k_p`$ always has to be outside and $`k_{pq}`$ has to be inside.

We now consider the two cases $`q=1,2`$ explicitely in order to calculate the slopes of both the type I and type II excitation branches at zero momentum. For $`q=1`$ corresponding to the smallest possible momentum transfer, the matrices $``$ and $``$ have one entry each and Eqs. (31) and (32) can be recast as a 2 by 2 matrix eigenvalue equation. This equation can be solved easily to yield the eigenvalues $`\mathrm{}\omega _1`$ and $`\mathrm{}\omega _1`$, where $`(\mathrm{}\omega _1)^2=\alpha ^2+2\alpha \nu v_{\mathrm{KS}}^{}/L`$ with $`\alpha =\mathrm{}^2N2\pi ^2/(mL^2\nu )`$ and only the positive solution is physically relevant. Using Eqs. (12) and (8) we find that $`\mathrm{}\omega _1/p_1=c_s`$ where $`p_1=2\pi \mathrm{}/L`$ is the momentum transfer for this excitation and $`c_s`$ is the speed of sound in the Lieb-Liniger model given by Eq. (20). We would like to point out that this result is independent of $`\nu `$ and can be traced back to the compressibility sum rule which is granted by the LDA.

We are now going to check this result for the speed of sound by calculating the energies for the next larger momentum transfer with $`q=2`$ where 2 particle-hole excitations give the first possibility to have different slopes for the type I and type II excitation branches. If both branches have the same slope and yield the same speed of sound, the energy difference between these two excitations should vanish in the thermodynamic limit. We have solved the eigenvalue equations (31) and (32) which are 4 by 4 with a specific symmetry. We omit the lengthy expressions but note the following observations: One of the solutions is independent of $`\nu `$ and gives the speed of sound as in Eq. (20) up to order $`1/L`$. The other solution does depend on $`\nu `$ and for the difference between the squared slopes we find

$$\left[(\mathrm{}\omega _2^{(1)})^2(\mathrm{}\omega _2^{(2)})^2\right]/p_2^2=\frac{n}{2m}v_{\mathrm{KS}}^{}+𝒪(1/L).$$
(38)
We see that this term remains finite in the thermodynamic limit unless $`v_{\mathrm{KS}}^{}`$ vanishes, which gives a stringent criterion for the choice of $`\nu =N/\stackrel{~}{N}`$ as discussed in Section 2.4. If the criterion of vanishing $`v_{\mathrm{KS}}^{}`$ is violated, the model will not only develop phonon branches with different speed of sound but additionally mean-field instabilities may occur as indicated by complex eigenvalues of the linear response equations. In the thermodynamic limit of the homogeneous case this occurs whenever $`\stackrel{~}{N}`$ is larger than the value determined by Eq. (14).

## References


format: text

text:
# Probing the Redshift Desert Using the Gemini Deep Deep Survey: observing galaxy mass assembly at 𝑧>1

## 1. Introduction: the ‘redshift desert’

The ‘redshift desert’ is the well known gap between $`z=1`$ and $`z=2`$ where there is a dearth of spectroscopic redshifts. This is illustrated in Figure 1 which shows a color-redshift plot of some of the deepest Keck spectroscopic samples now available. Generally galaxies fall in too two groups: pure magnitude selected samples which now reach as faint as $`R_{AB}<24`$ and redshifts $`z<1.2`$ and at higher redshifts the ‘Lyman Break’ galaxies. These are color-selected via their Lyman break and are easily found for $`z>3`$ (Steidel et al. 2003) and are being extended down to $`z=2`$ (Erb et al. 2003).

These samples are selected optically and as a result the galaxies selected at high-redshift have turned out to be those with the strongest ultraviolet (UV) emission because for $`z>1`$ the observed optical probes the rest-frame UV. Of note are the Lyman break galaxies which all have spectra characteristic of starburst galaxies: blue UV continua, ISM absorption lines and evidence of outflows.

Other high-redshift populations, which emit less strongly in the UV, have been found. Examples include the sub-mm selected SCUBA sources (Chapman et al. 2003) and the class of objects known as ‘Extremely Red Objects’ selected via their strong near-infrared (NIR) emission and red optical-NIR colors (e.g. Smail et al. 2002). While some spectra of such objects have been achieved there has been a lack of success in securing large samples of both red and blue galaxies above $`z=1`$.

The reason for this is primarily that the K-correction above $`z=1`$ has such a strong effect. Most galaxies apart from strong starbursts are significantly less luminous in the near-UV than at optical wavelengths, this is especially true of galaxies with SEDs characteristic of old stellar populations. Furthermore there is a lack of strong emission lines, even in starforming galaxies, between \[OII\] (3727Å) and Lyman-$`\alpha `$ (1216Å). Thus spectroscopy must reach sufficient signal:noise to robustly detect the continuum and hence achieve absorption line redshifts.

## 2. Galaxy Spectroscopy with Nod & Shuffle

Our goal with the ‘Gemini Deep Deep Survey’ (GDDS) was to achieve a complete NIR magnitude limited sample of both red and blue galaxies pushing well above $`z=1`$ with the aim of reaching $`z=2`$. In order to do this we need to be able to secure absorption line spectra of the reddest galaxies — candidate elliptical galaxies with redshifts $`z>1.5`$ which are the hardest to get redshifts of as all the flux comes out in the red and there are no emission lines. Therefore one needs efficient red spectroscopy. In principle the best option would be near-IR 1–2$`\mu `$m spectroscopy but too date there are no known examples of $`z>1`$ normal galaxies with absorption line NIR redshifts. This is because the sky is extremely bright in the NIR and because the detectors are very noisy. Moreover the field-of-view of the instrumentation is limited by the expense of NIR detector technology.

For the GDDS the compromise approach we adopted was to push on the red-end of the optical passband (5000-10 000Å) using CCD technology. The readnoise and sky brightness are reduced with respect to the near-IR, but so are the fluxes of the target objects. A simple estimate shows that a $`L^{}/3`$ $`z=1.5`$ elliptical could be as faint as $`I=24.5`$ requiring a $`100`$ ksec spectroscopic exposure to reach $`S/N=5`$ in the continuum at low resolution. Such long exposures are subject to many possible systematic effects making sky removal difficult — objects are $`<1\%`$ of the sky. However if such long exposures were possible then the optical approach would have a grant advantage: because CCDs are relatively inexpensive instruments typically offer much wider fields-of-view and hence greater multiplex advantage.

We adopted an approach known as ‘nod & shuffle’ (Glazebrook & Hawthorn 2001) to facilitate long exposures with accurate sky-subtraction in the red. (Figure 2). Briefly part of the CCD acts as a storage area and the telescope is nodded between an ‘object’ and a ‘sky position’. (In our case these were two neighboring positions along a slit). The spectrum at the sky position is stored in the storage area by ‘charge shuffling’ which introduces no additional noise provided a CCD has good charge transfer efficiency.

The nodding & shuffling is typically done many times throughout the course of a long exposure before reading out. Sky subtraction simply consists of subtracting the image from itself after being shifted by the shuffle distance.

For the GDDS we cooperated with the Gemini Observatory to implement nod & shuffle on the optical GMOS spectrograph (Crampton et al. 2000). This is notable as the first implementation of nod & shuffle on an 8m telescope. Commissioning tests in July 2002 indicated that sky-subtraction performance better than 0.1% was being achieved with a 60s object$`+`$60s sky nod & shuffle cycle. For more details about nod & shuffle with GMOS and the modes used with the GDDS see Abraham et al. (2003).

## 3. The Gemini Deep Deep Survey

The scientific design of the GDDS required spectroscopy of a near-IR selected sample to $`z=2`$. Targets were selected from the Las Campanas Infrared Survey which has deep optical and NIR data over many square degrees and has good photometric redshift estimates (Chen et al. 2003). The advantages of near-IR selection are well known (Glazebrook et al. 1994). If one selects in the $`K`$-band then the K-corrections are SED independent to roughly $`z=2`$. Moreover the light is dominated by the old stellar population, selecting in $`K`$ is a very good proxy for selecting by stellar mass. Thus a near-IR selected sample can select massive galaxies at high-redshift if they exist even if they are not actively star-forming and have weak UV. (The K20 survey has similar motivations and goals but is less deep; see Cimatti et al., 2002.) The survey design consisted of:

1. A primary $`K<20.8`$ selected sample. The photo-z distribution of the LCIRS indicated we needed to go this deep in order to penetrate to $`z=2`$. For the reddest galaxies at this limit $`I=24.5`$.
2. A fill-in sample of low-priority targets, anything with $`I<24.5`$ where there was space on the mask.
3. Photometric redshifts were used to weed out objects with $`z<0.8`$. This has to be done or else the low redshift foreground dominates. At low redshifts photo-z’s work the best because of the vastly better calibration, so this robustly leaves us with a $`0.8<z<2`$ sample.
4. Four fields were picked which lay within $`5^{}\times 5^{}`$ GMOS tiles, typically 70 slits could be allocated within each field. With nod & shuffle the slits only needed to be 2 arcsec long.
5. Each field was observed for $`100`$ ksec total using a 60s$`+`$60s nod & shuffle cycle, reading out every 1800s. Low resolution ($`R500`$) 5000–10,000Å GMOS spectra were taken.

All four fields were observed on Gemini-North between August 2002 and July 2003 and have been completely reduced. Redshifts and spectral classes were assigned by eye (by committee!) — in particular for the early-type galaxies by comparing with low-redshift templates. The final catalogs are currently being prepared for publication.

Some example spectra are shown in Figure 3, absorption line redshifts are obtained for both blue star-forming galaxies (ISM absorptions in FeII and MgII) and red early-type spectra (from matching to low-redshift UV templates which exhibit the same broad F features).

## 4. GDDS early results

Figure 4 shows the GDDS color-redshift plot from the first two fields for the $`K`$-selected primary sample. The spectroscopic completeness is 77%. There are 36 galaxies with spectroscopically confirmed galaxies with secure redshifts $`1.2<z<2`$. It can be seen that we are picking up a complete range of SED types to $`z1.7`$ and identifying objects much redder than Lyman break galaxies at comparable redshifts.

It is of course natural to ask ‘how massive are these galaxies’ ? The typical Lyman Break objects have stellar masses $`<10^{11}`$M (Papovich et al. 2001). We might expect that our objects being much brighter in the $`K`$-band will be considerably more massive.

To answer this question we have fit SEDs to our multicolor data to determine $`K`$-band mass-to-light values for the stellar population and hence derive stellar masses. Our procedure will be described in detail elsewhere; however we note that to first order the $`K`$-band $`M/L`$ values vary much less than they do in the optical (typically only a factor of two) and are far less sensitive to extinction and metallicity uncertainties. We do use throughout a full calculation fitting over a grid of 95,864 model galaxy SEDs computed using the PEGASE.2 code (Fioc & Rocca-Volmerange 1997). These models include a full range of dust extinction, metallicity, multiple bursts and different long-term star-formation histories.

Our results are shown in Figure 5. This shows the masses of our galaxies as a function of our $`K`$-band limit. It can be seen that for $`K<20.8`$ we reach a mass limit of $`>10^{10}`$M at $`z=2`$, and that a considerable number of galaxies have masses $`>10^{11}`$M, greatly in excess of the LBG population. This is natural given the galaxies are red and $`K`$-luminous. Many of these are spectroscopically confirmed to have old stellar populations. Such an abundance of massive galaxies may present a considerable problem for $`\mathrm{\Lambda }`$CDM models of hierarchical galaxy formation which typically predict a factor of $`30`$ drop in the number of $`10^{11}`$M galaxies between $`z=1`$ and $`z=2`$. (Baugh et al. 2003).

We are currently working on finalizing the space density of these massive galaxies at $`z>1.5`$, this work will be published as Glazebrook et al. (2004). The beauty of the GDDS is that we have spectra as well. We can assess whether a galaxy is spectroscopically old or a reddened starburst. We can also determine abundances from the spectra (see Savaglio et al. 2004, in press).

## 5. Conclusions

The Gemini Deep Deep Survey has succeeded in obtaining spectroscopic redshifts of galaxies in the former redshift desert. We have closed the $`1<z<2`$ gap, but more importantly we are highly complete for both red and blue galaxies. A catalog paper (Abraham et al. 2004) is in preparation. Our results to date are:

1. An abundance of massive galaxies at $`z>1.5`$. (Glazebrook et al. 2004, in prep.) These may be problematic for $`\mathrm{\Lambda }`$CDM models of heircharchical growth of massive halos.
2. Many of the massive red galaxies at $`z>1`$ have spectra of pure old stellar populations (McCarthy et al. 2003, in prep.)
3. Massive starburst galaxies at $`1.3<z<2`$ have high metal columns and are close to solar in abundances (Savaglio et al. 2004 in press).

At the time of this meeting we had just received the Gemini data for the final two of our four fields. The complete sample will be more than double the size of the one presented here. To come we have been awarded Cycle 12 time to image our fields with the Advanced Camera for Surveys on the Hubble Space Telescope. Thus we hope to connect our results to physical galaxy morphology. It is clear we are now living in an age where the formation of galaxies and the Hubble Sequence can be observed directly.

## 6. Acknowledgments

This work is based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the Particle Physics and Astronomy Research Council (United Kingdom), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), CNPq (Brazil) and CONICET (Argentina). Karl Glazebrook acknowledges genereous funding from the David & Lucile Packard Foundation.

## References

Abraham R.G., Glazebrook K., McCarthy P., Carlberg R., Chen H-W, Crampton D., Hook I., Jørgensen I., Marzke R., Murowinski R., Roth K., Savaglio S., 2003, Gemini Observatory Newsletter, 2003, 26, 6.

Baugh C. M., Benson A. J., Cole S., Frenk C. S., Lacey C., 2003, In: n the proceedings ‘The Mass of Galaxies at Low and High Redshift’, eds. R. Bender, A. Renzini, 91

Chapman S. C., Blain A. W., Ivison R. J., Smail I. R., 2003, Natur, 422, 695

Chen H., Marzke R. O., McCarthy P. J., Martini P., Carlberg R. G., Persson S. E., Bunker A., Bridge C. R., Abraham R. G., 2003, ApJ, 586, 745

Cimatti A., Mignoli M., Daddi E., Pozzetti L., Fontana A., Saracco P., Poli F., Renzini A., Zamorani G., Broadhurst T., Cristiani S., D’Odorico S., Giallongo E., Gilmozzi R., Menci N., 2002, A&A, 392, 395

Crampton D., Fletcher J. M., Jean I., Murowinski R. G., Szeto K., Dickson C. G., Hook I., Laidlaw K., Purkins T., Allington-Smith J. R., Davies R. L., 2000, SPIE, 4008, 114

Cohen J. G., 2001, AJ, 121, 2895

Erb D. K., Shapley A. E., Steidel C. C., Pettini M., Adelberger K. L., Hunt M. P., Moorwood A. F. M., Cuby J., 2003, ApJ, 591, 101

Fioc M., Rocca-Volmerange B., 1997, A&A, 326, 950

Glazebrook K., Peacock J. A., Collins C. A., Miller L., 1994, MNRAS, 266, 65

Glazebrook K., Bland-Hawthorn J., 2001, PASP, 113, 197

Papovich C., Dickinson M., Ferguson H. C., 2001, ApJ, 559, 620

Smail I., Owen F. N., Morrison G. E., Keel W. C., Ivison R. J., Ledlow M. J., 2002, ApJ, 581, 844

Steidel C. C., Giavalisco M., Pettini M., Dickinson M., Adelberger K. L., 1996, ApJ, 462, L17

Steidel C. C., Adelberger K. L., Shapley A. E., Pettini M., Dickinson M., Giavalisco M., 2003, ApJ, 592, 728


text:
# Detection of cool dust around the G2V star HD 107146

## 1 Introduction

Protostellar disks dissipate over a period of several Myr as their constituent dust and gas either accretes onto the star, is dispersed by processes such as stellar winds and photoevaporation, or aggregates into planetesimals (Hollenbach, Yorke, & Johnstone 2000). Possible observational examples of disk systems in the process of dissipating are the “transitional disks”, a thus far small class of objects which includes sources such as V819 Tau (Skrutskie et al. 1990), HR4796A (Jura et al. 1993; Jayawardhana et al. 1998; Koerner et al. 1998; Schneider et al. 1999; Telesco et al. 2000) and HD141569A (Weinberger et al. 1999; Augereau et al. 1999; Fisher et al. 2000).

Dramatic evidence for disk dispersal through planetesimal formation is the detection of numerous planets, and some planetary systems, around nearby stars (Marcy, Cochran, & Mayor 2000). Planetesimals are believed to grow into planets through collisional agglomeration. But collisions between planetesimals are also expected to create a cascade of smaller particles, and in this way a relatively old circumstellar disk may regenerate its own dust.

Such second generation dust was first detected around main-sequence stars by IRAS (Aumann et al. 1984). The IRAS results indicated that $`15`$% of main-sequence stars possess dusty disks (Lagrange, Backman, & Artymowicz 2000), the brightest of which have been detected at sub-millimeter wavelengths (Zuckerman & Becklin 1993; Greaves et al. 1998; Sylvester, Dunkin, & Barlow 2001; see also Wyatt, Dent, & Greaves 2003). Sub-millimeter observations are useful because they place strong constraints on the dust mass of the disk, due to the low optical depth of dust grains in the sub-millimeter compared to their opacities at shorter wavelengths. Because of their small angular size and weak emission, only four debris disk systems have been spatially resolved at these wavelengths to date (Holland et al. 1998; Greaves et al. 1998). The maps of these systems show interesting asymmetries that have been interpreted as the dynamical signature of a planetary companion (Holland et al. 2003; Wilner et al. 2002; Greaves et al. 1998).

Although the majority of the well-studied debris (and transitional) disk systems surround early-type stars (Vega, Fomalhaut, $`\beta `$ Pic, HR4796A, and HD141569A are all A stars), the expectation is that the process of disk dissipation and the regeneration of dust through planetesimal collisions is also a part of the evolutionary history of lower mass stars and, indeed, our own solar system. Thus, it is of considerable interest to identify debris disks associated with young solar-mass stars, since the detailed study of these objects can provide direct insight into the evolutionary history of our solar system. The $`ϵ`$ Eri system ($`d=3.2`$ pc) is one exciting example (Greaves et al. 1998; Hatzes et al. 2000). In this Letter, we present sub-millimeter observations of a dusty disk around a close solar analogue, the nearby young G2V star HD 107146.

HD 107146 was first identified as an “excess dwarf” on the basis of IRAS colors by Silverstone (2000). We observed it as part of a program to provide ground-based support for the “Formation and Evolution of Planetary Systems” SIRTF Legacy Project777Information concerning this program can be found at http://feps.as.arizona.edu. HD 107146 was selected for inclusion in the Legacy project based on its distance (28.5 pc; Perryman et al. 1997) and its high Ca II H and K index $`\mathrm{log}R_{\mathrm{HK}}^{}=4.28`$, which is indicative of youth. An age range can be estimated from its lithium equivalent width, 125 mÅ (Wichmann, Schmitt, & Hubrig 2003), which places it at the lower envelope of the 125 Myr Pleiades distribution (Soderblom et al. 1993) but above the distribution for the 250 Myr old M34 (Jones et al. 1997) at the effective temperature of a G2V star. Wichmann et al. also show that its space motions are similar to the Pleiades moving group. Its x-ray luminosity, $`L_X=2\times 10^{29}`$ erg s-1 (Voges et al. 1999), is again similar to solar-type Pleiads (Micela et al. 1999) but weaker than that of the average solar-type stars in the 80 Myr old Alpha Per cluster (Prosser et al. 1996). Together, these indicators suggest an age for HD 107146 between 80 and 200 Myr. Its location in the HR diagram is also consistent with this age range but allows the possibility that it may be as young as 30 Myr if the star has just reached the main sequence.

Our observations, made at sub-millimeter and mid-infrared wavelengths, are detailed in §2. The results from the imaging and fits to the spectral energy distribution (SED) are presented in §3. The SED shows clear evidence for a large inner hole, although the hole is too small to be resolved in the sub-millimeter maps presented here. In §4, we compare the properties of the HD107146 star-disk system to other debris and transitional disk systems. We also discuss prospects for future observations that may be able to determine whether the inner hole is maintained by the dynamical effect of an unseen orbiting companion.

## 2 Observations

HD 107146 was mapped using the SCUBA bolometer array at the James Clerk Maxwell Telescope (JCMT) on Mauna Kea, Hawaii in dry, stable conditions during February 17-19, 2003. Precipitable water vapor levels were less than 0.5 to 1 mm during the observations and zenith optical depths ranged from 0.1 to 0.18 at 850 µm and 0.4 to 0.85 at 450 µm. Observations were made simultaneously at 450 and 850 µm using a 64 point “jiggle” pattern to produce Nyquist sampled images at each wavelength. Two maps were rejected in the reduction process due to an error in the telescope tracking software at high elevation (Coulson, personal communication). The final maps presented here are the median of 23 maps made at a range of image rotation angles on the array rebinned to a $`1^{\prime \prime }`$ rectangular grid in equatorial coordinates. Pointing was checked between individual maps using Virgo A and 3C273 and was accurate to an rms error of $`2\stackrel{}{\mathrm{.}}5`$. Calibration was carried out by observing the CRL 618 planetary nebula and bootstrapped to observations of Mars later in the night. Based on the gain variations from night to night, the calibration accuracy is estimated to be 20% at 850 µm and 30% at 450 µm. The total on source integration time was 4.0 hours and the noise in the final maps was 4 and 13 mJy beam-1 at 850 and 450 µm respectively. The resulting peak signal-to-noise ratio is $`7`$ in both maps but the calibration uncertainty dominates the error in the absolute flux measurements.

We obtained mid-IR photometry of HD 107146 from the Keck II Telescope on Mauna Kea, Hawaii using the facility instrument LWS. We observed the star on February 19 and 20, 2003 UT using filters centered at 11.7 µm (10.5–12.9 µm) and 17.8 µm (17.3–18.2 µm). Contemporaneously and at similar airmasses, we observed the bright standard stars $`\alpha `$ CrB and $`\mu `$ UMa from Tokunaga (1988) for photometric calibration. Both nights were characterized by dry conditions with an estimated 0.5–1.5 mm of precipitable water, as determined from sky dip measurements from the JCMT. Seeing conditions were poor and variable, leading to mid-IR images with 0.4″ to 0.6″ FWHM. Observations were conducted in the standard “chop/nod” mode, which involves switching between three sky positions using fast chopping of the secondary mirror and slower nodding of the telescope itself. This allows for effective subtraction of the very bright and variable thermal emission from the sky and telescope in the reduction process.

HD 107146 was well-detected in both filters, with a formal signal-to-noise ratio $`>`$15 for the 17.6 µm data and $`>`$50 for the 11.7 µm data. However, the seeing conditions limit the photometric precision, due to varying image quality from the science target to the standard stars. These conditions also prevented any useful constraints on any extended mid-IR emission from HD 107146. Photometric errors were first determined from the scatter in measurements done with apertures scaled by the FWHM of the images. We then added in quadrature an error term to account for seeing mismatch (10 mJy at 11.7 µm, 8 mJy at 17.6 µm), estimated from analyzing curve of growth photometry for standard stars observed during the course of the entire night. The absolute flux calibration is based on the flux of Vega ($`\alpha `$ Lyr) compiled by Tokunaga (1988). The LWS and SCUBA fluxes are listed in Table 1.

## 3 Results

The resulting maps of the emission at 450 and 850 µm are shown in Figure 1. There is a $`2\stackrel{}{\mathrm{.}}9`$ offset between the two maps, as measured by gaussian fits to the data clipped at the FWHM level. The centroid of the emission is consistent with the stellar position to within $`4\stackrel{}{\mathrm{.}}4`$ (450 µm) and $`1\stackrel{}{\mathrm{.}}6`$ (850 µm). It is not clear how significant are such small offsets given the relatively low signal-to-noise in the data and the possibility of systematic pointing errors as large as $`2\stackrel{}{\mathrm{.}}5`$. No other source within this angular range of the star was apparent in the Digital Sky Survey, 2MASS images, the Keck imaging at 11.7 and 17.8 µm, and the FIRST 20 cm survey. In addition, Blain et al. (1999) predict $`100300`$ sources deg-2 with fluxes greater than 20 mJy at 850 µm, implying a probability less than $`1.4\times 10^3`$ that an unrelated background object this bright would be found within $`4\stackrel{}{\mathrm{.}}4`$ of the star. Since the probability of unrelated submillimeter emission is very low, we assume that the SCUBA source is a disk associated with the star and discuss its properties in this context.

Allowing for pointing errors, the resolutions of the maps are $`14\stackrel{}{\mathrm{.}}5`$ at 850 µm and $`8^{\prime \prime }`$ at 450 µm. The gaussian fits give sizes of $`15\stackrel{}{\mathrm{.}}1\times 14\stackrel{}{\mathrm{.}}9`$ and $`13\stackrel{}{\mathrm{.}}2\times 10\stackrel{}{\mathrm{.}}9`$ (position angle $`35^{}`$) respectively. The elongation at 450 µm is apparent at the half-power level and is therefore not due to the JCMT beam pattern which shows significant non-circularity only at the 10% level. The slight extension in the 450 µm map is also seen in crosscuts through the image along the major axis. To increase the signal-to-noise ratio, an average crosscut over a range of position angles $`35^{}\pm 10^{}`$ centered on the peak of emission is shown for each wavelength in the lower panels of Figure 1. The same crosscut averaging, shown as dashed lines, was performed for a map of Mars, taken from the last night of observations, February $`19^{\mathrm{th}}`$, when the Martian diameter was $`5\stackrel{}{\mathrm{.}}7`$. We conclude that the disk around HD 107146 is marginally resolved at the $`8^{\prime \prime }`$ resolution of the 450 µm data. Subtracting the beam size in quadrature from the gaussian fit gives an angular size of $`10\stackrel{}{\mathrm{.}}5\times 7\stackrel{}{\mathrm{.}}4`$ for the disk corresponding to $`300\mathrm{AU}\times 210\mathrm{AU}`$.

The fluxes measured from the LWS and SCUBA observations are listed in Table 1. The SED of the source, from optical to sub-millimeter wavelengths, is plotted in Figure 2. The UBVRI photometry is from Landolt (1983), the near-infrared fluxes in the JHK bands are from 2MASS, and the 10 µm point is from Palomar observations (Metchev, Hillenbrand, & Meyer 2003). IRAS fluxes were determined by color correcting the quoted values in the Faint Source Catalog. The IRAS 12 µm flux is consistent with the Keck 11.7 µm observation and is not plotted. The stellar photosphere was fit by a Kurucz model ($`T_{\mathrm{eff}}=5750\mathrm{K},\mathrm{log}g=4.5`$, solar metallicity) with a power law extrapolation beyond 10 µm. Strong excess emission is apparent beyond 25 µm. In order to compare the properties of the excess with that detected from other debris disk systems, we fit the disk SED using a single temperature modified black body with emission efficiency, $`Q_\lambda =1\mathrm{exp}[(\lambda _0/\lambda )^\beta ]`$, which has asymptotic behavior, $`Q_\lambda =1`$ for $`\lambda \lambda _0`$ and $`Q_\lambda =(\lambda _0/\lambda )^\beta `$ for $`\lambda \lambda _0.`$ The critical wavelength, $`\lambda _0`$, was set to 100 µm for consistency with the assumptions made in previous analysis of debris disks SEDs (Dent et al. 2000; Wyatt, Dent & Greaves 2003).

The parameters of the modified black body fit were measured using a least squares fit to the data. Errors were estimated via fits to multiple simulations of the dust excess SED. Simulated data points were drawn from a gaussian distribution with the mean and standard deviation as determined by the observed data. The resulting distribution of best fit parameter values are $`T=51\pm 4`$ K, $`\beta =0.69\pm 0.15`$. Similarly low values of $`\beta `$ are found for other disks around main-sequence stars (Dent et al. 2000). The dust mass, based on the fitted flux at 850 µm and assuming a dust mass absorption coefficient, $`\kappa _{850}=1.7`$ cm2 g-1, is $`M_\mathrm{d}=0.10\pm 0.02M_{}`$. This value of $`\kappa _{850}`$ is chosen for consistency with Holland et al. (1998) and Greaves et al. (1998), but is on the high end of calculated values (Pollack et al. 1994). For the full range of $`\kappa _{850}=0.41.7`$ cm2 g-1 discussed in Pollack et al., the corresponding mass range is $`M_\mathrm{d}=0.100.43M_{}`$. As with all sub-millimeter observations, the inferred dust masses do not include a potentially dominant mass component that resides in larger bodies (grains and planetesimals) subtending a negligible solid angle. The mean parameter fit is shown in Figure 2.

The single temperature fit is a simplification that is warranted by the small number of data points. Nevertheless, the large dip in the SED at 25 µm imposes a strong limit on the mass of warmer dust that may be present in the system. To illustrate this constraint, a $`T=100`$ K, $`\beta =0.7`$ modified black body component was added to the SED fit and increased until the IRAS 25 µm upper limit was exceeded. This maximum allowable contribution, which corresponds to a mass limit of $`M_\mathrm{d}(T=100\mathrm{K})=7\times 10^4M_{},`$ is shown in Figure 2 but was not included in the overall fit that is shown. The Keck 18µm measurement tightly constrains the presence of still warmer dust to a $`3\sigma `$ limit of $`M_\mathrm{d}(T=200\mathrm{K})<3\times 10^5M_{}`$.

The limits on dust cooler than 50 K are less stringent since such dust emits less per unit mass and does so at longer wavelengths. The minimum grain temperature is 23 K for blackbody particles at the measured outer radius of the disk, 150 AU from the star. In practice, the effect of a range of cool dust temperatures is indistinguishable from changes in the wavelength dependence of the grain emissivity, parameterized by $`\beta `$. However, without changing the parameters of the fit in Figure 2, $`0.13M_{}`$ of 23 K dust can be added before the 850 µm $`3\sigma `$ upper limit is exceeded. Thus, a significant fraction of the total dust mass may reside in cooler dust.

## 4 Discussion and Summary

We have detected strong sub-millimeter excess emission from HD 107146, a G2V star with an age estimated to lie in the range $`80200`$ Myr. Based on a fit to the SED of the system, the mass of the emitting dust is estimated to be $`0.10M_{}`$ or larger. We also find that the disk is marginally resolved at 450 µm. The constraints placed by these observations on the mass, temperature, and physical extent of the dust in the system are compared in Table 2 with the properties of other well-studied debris disk systems.

Four of the six stars in Table 2 are A stars, a probable bias due to the relatively high luminosity of these systems at far-infrared and sub-millimeter wavelengths. Nevertheless, despite the range in central star masses and luminosities, the disks have similar properties and follow several trends. The HD 107146 disk is quite massive, comparable in mass to the $`\beta `$ Pic disk. Although it lies noticeably above the trend of decreasing mass with age that is defined by HR4796A, $`\beta `$ Pic, Fomalhaut, and Vega (Holland et al. 1998), it is within the scatter in the mass-age relation found for larger samples of dust disks (Wyatt et al. 2003). The fractional dust luminosity ($`L_d/L_{}`$) of HD107146 is also large; it is half that of $`\beta `$ Pic, $`>10`$ times that of Fomalhaut, and is consistent with the trend of decreasing $`L_d/L_{}`$ with age (Spangler et al. 2001). Table 2 also indicates an apparent trend of decreasing outer disk radius with age. Perhaps this trend, admittedly of low statistical significance at present, will be verified when a larger number of debris disk systems have been spatially resolved.

We also find that the local minimum in the SED at 25 µm places a strong limit on the amount of warm dust in the HD 107146 disk, $`M_\mathrm{d}(\mathrm{T}=51\mathrm{K})/M_\mathrm{d}(\mathrm{T}=100\mathrm{K})140`$. This lack of warm dust implies that the disk does not extend all the way to the star. The dust temperature, and therefore the size of this inner hole, depends on the grain size distribution and optical properties. For grains that emit as blackbodies at 850 µm, 51 K dust would lie at 31 AU (Wyatt et al. 1999). This is a lower limit to the inner radius as smaller grains could achieve this temperature at greater distances from the star. The other disks listed in Table 2 have similar size inner holes, as determined from spatially resolved images in thermal emission or scattered light, and are comparable in size to the Kuiper belt in our solar system.

Inner holes are not, by themselves, long lasting because Poynting-Roberston drag will cause dust to spiral in from the outer disk onto the central star on $`10`$ Myr timescales. Such inward migration would be evident in the SED by the presence of warm dust (Jura et al. 1998). Thus, the existence of inner holes has been explained as a consequence of either the dynamical sweeping of an orbiting companion (e.g., Vega; Wilner et al. 2002) or the sublimation of icy grains (e.g., HR 4796A; Jura et al. 1998). Since water ice sublimates at temperatures greater than 100 K (Pollack et al. 1994), this explanation can not apply to the HD 107146 disk where the shape of the SED places a strong limit on the mass of such warm dust. The alternative explanation, that the inner hole is dynamically maintained by a closely orbiting companion, would require a fairly low-mass companion. The best current limit on the existence of close companions is from Palomar adaptive optics imaging by Metchev & Hillenbrand (2002) who found no companions at detection limits of 11.2, 11.7, and 15.2 in absolute K band magnitude at angular separations of $`0\stackrel{}{\mathrm{.}}5,1^{\prime \prime }`$, and $`2^{\prime \prime }`$ respectively. The corresponding mass limits for an age of 100 Myr are approximately 30, 25, and $`10M_\mathrm{J}`$ (Burrows et al. 1997).

While it has been well-recognized that the dynamical sculpting of disks by orbiting companions can produce inner holes and ring-like structures, recent studies have shown that the migration of dust in the presence of a residual gas disk can also induce ring-like structures in the dust distribution. For example, Takeuchi & Artymowicz (2001) have shown that the dust structures seen in the HR4796A and HD141569A systems are qualitatively similar to those expected to result from the coupling between gas and dust in disks. However, orbiting companions may also induce significant departures from axisymmetry in the dust distribution (e.g., Liou & Zook 1999; Ozernoy et al. 2000; Quillen & Thorndike 2002; Moro-Martin & Malhotra 2002), whereas such non-axisymmetric structures cannot be produced by dust migration. Thus, demonstrating the existence of asymmetries in the dust distribution, either in thermal emission or scattered light, as well as measuring the gas content of the HD107146 disk, is needed to distinguish between these two possibilities.

At a distance of 28.5 pc, HD 107146 is relatively nearby. Compared to the other sources listed in Table 2, it is more distant than the debris disk systems that have been spatially resolved at submillimeter wavelengths ($`\beta `$ Pic, Vega, Fomalhaut, and $`ϵ`$ Eri), but it is at half the distance of HR4796A. Since the HD 107146 disk is just resolved at the shortest operating wavelength of SCUBA, and the submillimeter excess is relatively bright, the Submillimeter Array should be able to map the morphology of the emitting dust in greater detail. Phase referencing will provide a more accurate absolute position, and such observations could determine whether the dust is distributed axisymmetrically (e.g., in a ring) or in a more asymmetric distribution. For example, these observations may confirm the marginal evidence for an offset between the stellar position and the peak of the $`450\mathrm{\mu m}`$ emission found in the observations presented here (Figure 1). The offset, if real, may result from an asymmetric dust distribution, as has been found for several of the other debris disk systems in Table 2.

The dust disk asymmetry can also be addressed by imaging the disk in scattered light. The possibility of detecting scattered light from Vega-like stars is of great interest given the few such systems that have been detected thus far, despite extensive deep surveys (e.g., Kalas & Jewitt 1996). Using the quantities listed in Table 2, we can make a rough estimate of the relative strengths of the scattered light from each disk. If we ignore the disk inclination and assume a similar dust grain size distribution and albedo, the scattering area will be proprtional to $`M_d`$ and the flux of scattered light will be roughly proportional to $`M_dL_{}/R_{\mathrm{out}}^2d^2`$, where $`d`$ is the distance to the star and $`R_{\mathrm{out}}`$ is the physical extent of the disk. The factor of $`R_{\mathrm{out}}^2`$ arises in this expression because grains of a given size located farther from the star intercept less starlight. Hence the surface brightness of the disk will be proportional to $`M_dL_{}/R_{\mathrm{out}}^4`$. A potentially more relevant comparison is of the contrast between the expected surface brightness from the disk and the flux from the star, which is proportional to $`M_dd^2/R_{\mathrm{out}}^4`$ (cf. Jura et al. 1998). For the disk around HD107146, the expected surface brightness of the scattered light is one half that of the disk around $`\beta `$ Pic and the expected contrast is nine times larger, suggesting that the scattered light from HD107146 may not only be quite bright but may also stand out against the glare of the star. Given the low dust temperature deduced for the HD107146 disk ($`<100`$K), the grains are expected to be icy with a high albedo, which favors the detection of scattered light.

Nevertheless, there are significant uncertainties associated with this estimate. Detecting reflected light will be more challenging if the disk is face-on rather than edge-on. In addition, the sub-millimeter measurements from which $`M_d`$ is derived are primarily sensitive to large grains ($`100\mathrm{\mu m}`$), whereas the scattered light observations will be particularly sensitive to much smaller grain sizes. Thus, scattered light from the HD107146 disk will be weaker if the grain size distribution is significantly skewed to large grain sizes. Conversely, scattered light measurements can help to constrain the grain size distribution of the disk (e.g., Artymowicz, Burrows, & Paresce 1989).

Future spectroscopic observations of this disk with SIRTF, particularly in the $`3040`$ µm range, will place additional constraints on the grain composition and size distribution (e.g. Wolf & Hillenbrand 2003). SIRTF is also expected to discover many more disks over a large range of ages and stellar masses with consequent improvements for our understanding of the formation and evolution of planetary systems.

We thank Anneila Sargent for counseling during an extended review process, Hervé Aussel, Remo Tilanus, Iain Coulson, and Randy Campbell for advice on instrumentation, Eric Mamajek for helpful discussions on the age of HD 107146, and Alan Tokunaga for making the Keck time available. We acknowledge support from NSF grant AST-0324328 (JPW), the Beatrice Watson Parrent Fellowship at the University of Hawaii and NASA grant HST-HF-01152.01 (MCL), and NASA contract 1224768 administered through JPL (MRM, JMC, LAH). This research has made use of the SIMBAD database.

## 5 References

Artymowicz, P., Burrows, C., & Paresce, F. 1989, ApJ, 337, 494

Augereau, J. C., Lagrange, A. M., Moillet, D., & Menard, F. 1999, A&A, 350, L51

Aumann, H. H., Beichman, C. A., Gillett, F. C., de Jong, T., Houck, J. R., Low, F. J., Neugebauer, G., Walker, R. G., & Wesselius, P. R. 1984, ApJ, 278, L23

Blain, A. W., Kneib, J.-P., Ivison, R. J., & Smail, Ian 1999, ApJ, 512, L87

Burrows, A., Marley, M., Hubbard, W. B., Lunine, J. I., Guillot, T., Saumon, D., Freedman, R., Sudarsky, D., & Sharp, C. 1997, ApJ, 491, 856

Dent, W. R. F., Walker, H. J., Holland, W. S., & Greaves, J. S. 2000, MNRAS, 314, 702

Greaves, J. S. et al. 1998, ApJ, 506, L133

Greaves, J. S., Mannings, V., & Holland, W. S. 2000, Icarus, 143, 155

Fisher, R. S., Telesco, C. M., Pina, R. K., Knacke, R. F., & Wyatt, M. C. 2000, ApJ, 532, L141

Hatzes, A.P., et al. 2000, ApJ, 544, L145

Holland, W. S. et al. 1998, Nature, 392, 788

Hollenbach, D. J., Yorke, H. W., & Johnstone, D. 2000, in “Protostars and Planets IV”, ed. V. Mannings, A.P. Boss & S.S. Russell (Tucson: University of Arizona Press), 401

Jayawardhana, R. Fisher, S., Hartmann, L., Telesco, C., Pina, R., & Fazio, G. 1998, ApJ, 503, L79

Jones, B. F., Fischer, D., Shetrone, M., & Soderblom, D. R. 1997, AJ, 114, 352

Jura, M., Zuckerman, B., Becklin, E. E., Smith, R. C. 1993, ApJ, 418, L37

Jura, M., Malkan, M., White, R., Telesco, C., Pina, R., & Fisher, R. S. 1998, ApJ, 505, 897

Kalas, P., & Jewitt, D., 1996, AJ, 111, 1347

Koerner, D. W., Ressler, M. E., Werner, M. W., & Backman, D. E. 1998, ApJ, 503, L83

Lagage, P. O., & Pantin, E. 1994, Nature, 369, 628

Lagrange, A.-M., Backman, D. E., & Artymowicz, P. 2000, in “Protostars and Planets IV”, ed. V. Mannings, A.P. Boss & S.S. Russell (Tucson: University of Arizona Press), 639

Landolt, A. U. 1983, AJ, 88, 853

Liou, J.-C. & Zook, H. A. 1999, AJ, 118, 580

Marcy, G. W., Cochran, W. D., & Mayor, M. 2000, in “Protostars and Planets IV”, ed. V. Mannings, A.P. Boss & S.S. Russell (Tucson: University of Arizona Press), 1285

Metchev S. A. & Hillenbrand, L. A. 2003, in “Debris disks and the formation of planets”, eds. L. Caroff & D. Backman, ASP Conf. Series

Metchev S. A., Hillenbrand, L. A., & Meyer, M. R. 2003, ApJ, submitted

Micela, G., Sciortino, S., Harnden, F. R., Kashyap, V., Rosner, R., Prosser, C. F., Damiani, F., Stauffer, J., & Caillault, J.-P. 1999, A& A, 341, 751

Moro-Martin, A., & Malhotra, R. 2002, AJ, 124, 2305

Ozernoy, L. M., Gorkavyi, N. N., Mather, J. C., & Taidakova, T. T. 2000, ApJ, 537, L147

Pantin, E., Lagage, P. O., & Artymowicz, P. 1997, A&A, 327, 1123

Perryman, M. A. C. et al. 1997, A& A, 323, L49

Pollack, J. B., Hollenbach, D., Beckwith, S., Simonelli, D. P., Roush, T., & Fong, W. 1994, ApJ, 421, 615

Prosser, C. F., Randich, S., Stauffer, J. R., Schmitt, J. H. M. M., & Simon, T. 1996, AJ, 112, 1570

Quillen, A. C., & Thorndike, S. 2002, ApJ, 578, L149

Schneider, G., et al. 1999, ApJ, 513, L127

Silverstone, M. 2000, UCLA PhD thesis

Skrutskie, M. F., Dutkevich, D., Strom, S. E., Edwards, S., Strom, K. M., & Shure, M. A. 1990, AJ, 99, 1187

Soderblom, D. R., Jones, B. F., Balachandran, S., Stauffer, J. R., Duncan, D. K., Fedele, S. B., & Hudon, J. D 1993, 105 2299

Spangler, C., Sargent, A. I., Silverstone, M. D., Becklin, E. E., & Zuckerman, B. 2001, ApJ, 555, 932

Sylvester, R. J., Dunkin, S. K., & Barlow, M. J. 2001, MNRAS, 327, 133

Telesco, C. M., et al. 2000, ApJ, 530, 329

Tokunaga, A. T. 1988, Infrared Telescope Facility (IRTF) Photometry Manual

Voges, W. et al. 1999, A& A, 349, 389

Weinberger, A. J., et al. 1999, ApJ, 525, L53

Wichmann, R., Schmitt, J. H. M. M., & Hubrig, S., 2003, A& A., 399, 983

Wilner, D. J., Holman, M. J., Kuchner, M. J., & Ho, P. T. P. 2002, ApJ, 569, L115

Wolf, S., & Hillenbrand, L. A. 2003, ApJ, in press

Wyatt, M. C., Dent, W. R. F., & Greaves, J. S. 2003, MNRAS, 342, 876

Wyatt, M. C., Dermott, S. F., Telesco, C. M., Fisher, R. S., Grogan, K., Holmes, E. K., & Pina, R. K. 1999, ApJ, 527, 918

Zuckerman, B., & Becklin, E. E. 1993, ApJ, 414, 793

Figure 1: SCUBA images of the dust emission around HD 107146. Coordinates are offset from the J2000 position of the star epoch 2003.13, $`\alpha =12^\mathrm{h}19^\mathrm{m}6.46^\mathrm{s},\delta =16^{}32^{}53\stackrel{}{\mathrm{.}}4`$, indicated by the star symbol. The left panel shows the emission at 850 µm. Contour levels begin at $`2\sigma `$ and increment by $`\sigma `$, where $`\sigma =4`$ mJy beam-1 is the noise in the map. The greyscale runs from $`\sigma `$ to $`7\sigma `$. The right panel shows the emission at 450 µm; contours begin at $`3\sigma `$ and increment by $`\sigma `$, where $`\sigma =13`$ mJy beam-1. The greyscale runs from $`2\sigma `$ to $`8\sigma `$. The beam sizes, $`14\stackrel{}{\mathrm{.}}5`$ at 850 µm and $`8^{\prime \prime }`$ at 450 µm, are indicated by the hashed circles in the lower right corner of each figure. Inspection shows that the half power point of the 450 µm image ($`3.5\sigma `$) is slightly greater than the beam size. Averaged cuts of the source and Mars are shown in the lower panels for each wavelength. The cuts are centered on the peak of emission and averaged over positions angles $`35^{}\pm 10^{}`$, outlined by dotted lines in each image. The Mars profile is shown as the dashed line and is significantly narrower than HD107146 in the higher resolution 450 µm image.

Figure 2: Optical to sub-millimeter spectral energy distribution of HD 107146. The open circle at 25 µm is an upper limit from the IRAS Faint Source catalog. All other points are detections with error bars shown when the error exceeds the symbol size. The double peaked distribution is modeled as the sum of a Kurucz model of the stellar spectrum and a modified black body fit to the points longward of 25 µm ($`T=51\mathrm{K},M_\mathrm{d}=0.10M_{},\beta =0.7`$). The individual contribution of each component is shown by the dotted lines and the sum of the two by the solid line. The dashed line shows the maximum allowable $`T=100`$ K dust component that fits the constraint of the IRAS 25 µm upper limit. This component is not included in the overall fit shown by the solid line.