text:
# Boundary contributions to specific heat and susceptibility in the spin-1/2 XXZ chain

(October 22, 2003)

## Abstract

Exact low-temperature asymptotic behavior of boundary contribution to specific heat and susceptibility in the one-dimensional spin-1/2 XXZ model with exchange anisotropy $`1/2<\mathrm{\Delta }1`$ is analytically obtained using the Abelian bosonization method. The boundary spin susceptibility is divergent in the low-temperature limit. This singular behavior is caused by the first-order contribution of a bulk leading irrelevant operator to boundary free energy. The result is confirmed by numerical simulations of finite-size systems. The anomalous boundary contributions in the spin isotropic case are universal.

Irrelevant operators in the renormalization-group (RG) theory are considered to give only subleading contributions to any physical quantity by definition. A subtle exception to this general rule can be found, however, when a bulk irrelevant operator has not so large boundary scaling dimension. Consider, for example, a Gaussian theory in (1+1)-dimension with a bulk perturbing operator $`O`$ having scaling dimension, $`x_o`$, larger than 2. First, let us impose periodic boundary condition to a one-dimensional (1D) system of length $`L`$. In second-order perturbation the irrelevant operator is found to give a contribution of order $`LT^{2x_o3}`$ to the specific heat of the periodic system at temperature $`T`$. This is small compared with the leading term of order $`LT`$. Next, let us suppose that the system has open boundaries. In general the specific heat has boundary contribution of order $`L^0`$. Indeed, first-order perturbation leads to a boundary contribution of order $`L^0T^{x_o^{}2}`$ to the specific heat, where $`x_o^{}`$ is boundary scaling dimension of the perturbing operator. If $`2<x_o^{}<3`$, then the boundary specific heat gives a dominant contribution at low enough temperatures, even though the operator is irrelevant.

This interesting possibility is in fact realized in the 1D spin-1/2 XXZ model, as recently pointed out by Fujimoto.Fujimoto  The model with open boundaries was studied earlier by de Sa and Tsvelik using Bethe ansatz,deSa  who found anomalous boundary contributions to the specific heat and spin susceptibility. Similar and more extensive calculationsEssler ; Asakawa ; Zvyagin ; FujimotoPRB  for the spin isotropic case have shown that the boundary contribution to the spin susceptibility at $`T=0`$ exhibits a Curie-like behavior with a logarithmic correction as a function of a weak external field. These results seemed somewhat mysterious, as irrelevant operators localized at a boundary can give at most a temperature-independent contribution to the susceptibility.EggertA ; FuruNaga ; Frojdh  The origin of these anomalous boundary contributions was not understood, until Fujimoto and Eggert recently examined boundary contributions from a bulk leading irrelevant operator to the specific heat and the susceptibility.Fujimoto ; FujimotoEggert  Their theory seems rather complicated, however, as it uses properties of boundary states in the boundary conformal field theory. The aim of this paper is to present an alternative, simpler, and straight-forward derivation of the exact results using the Abelian bosonization method. We also verify our analytical results with numerical calculations for the anisotropic XXZ model.

Our results are relevant to many 1D correlated electron systems with impurities. An impurity potential in a 1D system of repulsively interacting electrons is strongly renormalized and becomes perfectly reflecting in the low-energy limit.Kane ; FuruNaga93 ; Matveev  A long wire with a finite density of impurities is thus effectively cut into many pieces of finite length. Measurements of the specific heat and the susceptibility of such a wire will exhibit the anomalously large boundary contributions.

We analyze the spin-1/2 XXZ model in a weak magnetic field,

$`H`$ $`=`$ $`J{\displaystyle \underset{j=1}{\overset{L1}{}}}\left(S_j^xS_{j+1}^x+S_j^yS_{j+1}^y+\mathrm{\Delta }S_j^zS_{j+1}^z\right)h{\displaystyle \underset{j=1}{\overset{L}{}}}S_j^z`$ (1)
$`=`$ $`H_0h{\displaystyle \underset{j=1}{\overset{L}{}}}S_j^z,`$
where $`J>0`$ and $`L`$ is the number of sites. We are interested in the case where the exchange anisotropy $`\mathrm{\Delta }`$ is in the parameter range

$$\frac{1}{2}<\mathrm{\Delta }1,$$
(2)
for which we shall find an anomalous contribution to the boundary free energy.

First we consider a semi-infinite spin chain, taking the limit $`L\mathrm{}`$. Its free energy has, in addition to a bulk term, a boundary term of order 1 due to the presence of the boundary at $`j=1`$. We calculate the boundary free energy at low temperatures using the low-energy effective theoryAffleck  for the Hamiltonian (1), which is the Gaussian model perturbed by a nonlinear operator. The Hamiltonian density of the effective theory is

$$=\frac{v}{2}\left[\left(\frac{d\varphi }{dx}\right)^2+\mathrm{\Pi }^2\right]+\lambda \mathrm{cos}\left(\frac{2\varphi }{R}\right)\frac{h}{2\pi R}\frac{d\varphi }{dx}.$$
(3)
The lattice constant is unity. The bosonic field $`\varphi (x)`$ and its conjugate operator $`\mathrm{\Pi }(x)`$, defined on the positive half line, satisfy $`[\varphi (x),\mathrm{\Pi }(y)]=i\delta (xy)`$. The parameter $`R`$ and the velocity $`v`$ are functions of $`\mathrm{\Delta }`$ and $`h`$. At zero field $`h=0`$, the parameter $`\eta `$ defined by $`\eta =2\pi R^2`$ is

$$\eta =1\frac{1}{\pi }\mathrm{cos}^1(\mathrm{\Delta }),$$
(4)
with which the spin velocity $`v`$ at $`h=0`$ is written,

$$v=J\frac{\mathrm{sin}(\pi \eta )}{2(1\eta )}.$$
(5)
We note that in the parameter range (2) of our interest $`\eta `$ satisfies the inequality

$$\frac{2}{3}<\eta 1$$
(6)
at $`h=0`$. The original spin operator can be expressed with the fields $`\varphi `$ and $`\mathrm{\Pi }`$. In particular, $`S_j^z`$ is written as

$$S_j^z=\frac{1}{2\pi R}\frac{d\varphi }{dx}+(1)^{j1}a\mathrm{sin}\frac{\varphi (x)}{R},$$
(7)
where $`x=j`$ and $`a`$ is a constant calculated in Refs. Lukyanov99,  and corAmp, . For the semi-infinite wire the bosonic field obeys the boundary condition at $`x=0`$,EggertA ; corAmp

$$\varphi (x=0)=0.$$
(8)
Once Eqs. (7) and (8) are given, sign of the coupling constant $`\lambda `$ of the nonlinear operator is fixed to be positive. For the gapless regime ($`1<\mathrm{\Delta }1`$) of the XXZ chain the exact value of $`\lambda `$ is known to beLukyanov98 ; LukyanovT

$$\frac{\lambda }{v}=2\mathrm{sin}\left(\frac{\pi }{\eta }\right)\left[\frac{\mathrm{\Gamma }(\frac{1}{\eta })}{\pi }\right]^2\left[\frac{\mathrm{\Gamma }(1+\frac{\eta }{22\eta })}{2\sqrt{\pi }\mathrm{\Gamma }(1+\frac{1}{22\eta })}\right]^{\frac{2}{\eta }2}$$
(9)
under the regularization condition on the zero-temperature correlator in the Gaussian model,

$$e^{i\mu \varphi (x)}e^{i\mu \varphi (y)}=|xy|^{\mu ^2/2\pi },x,y,|xy|1.$$
(10)
We are now ready to calculate the leading contribution to the boundary free energy. We begin with the path integral representation of the partition function $`Z`$,

$$Z=𝒟\varphi \mathrm{exp}\left(_0^{1/T}𝑑\tau _0^{\mathrm{}}𝑑x\right),$$
(11)
where the Lagrangian density in the imaginary time is

$$=\frac{v}{2}\left[\left(\frac{\varphi }{x}\right)^2+\frac{1}{v^2}\left(\frac{\varphi }{\tau }\right)^2\right]+\lambda \mathrm{cos}\left(\frac{2\varphi }{R}\right)\frac{h}{2\pi R}\frac{\varphi }{x}.$$
(12)
To eliminate the last term ($`h`$), we change the field as

$$\varphi (x,\tau )\varphi (x,\tau )+\frac{hx}{2\pi Rv},$$
(13)
which yields a constant contribution $`\chi _0h^2/2`$ to the bulk free energy with the bulk susceptibility given by $`\chi _0=L/(2\pi \eta v)`$. This term is not of our main interest and will be discarded in the following analysis.

Since the operator $`\lambda \mathrm{cos}(2\varphi /R)`$ is an irrelevant perturbation in the gapless phase, the low-temperature expansion of the free energy can be obtained by expanding the partition function in powers of $`\lambda `$. The leading boundary contribution to the free energy then comes from the first-order term $`\lambda `$ in the expansion. After making the shift (13), the first-order term reads

$$Z_1=\lambda 𝒟\varphi _0^{1/T}𝑑\tau _0^{\mathrm{}}𝑑x\mathrm{cos}\left(\frac{2\varphi }{R}+\frac{2hx}{\eta v}\right)\mathrm{exp}\left\{\frac{v}{2}_0^{1/T}𝑑\tau _0^{\mathrm{}}𝑑x\left[\left(\frac{\varphi }{x}\right)^2+\frac{1}{v^2}\left(\frac{\varphi }{\tau }\right)^2\right]\right\}.$$
(14)
Under the boundary condition (8) the bosonic field can be expanded as

$$\varphi (x,\tau )=\sqrt{\frac{2T}{\pi }}\underset{n}{}_0^{\mathrm{}}𝑑ke^{i\omega _n\tau }\mathrm{sin}(kx)\stackrel{~}{\varphi }(k,i\omega _n),$$
(15)
where $`\omega _n=2\pi nT`$. Inserting the mode expansion (15) into Eq. (14) and performing the Gaussian integral of $`\stackrel{~}{\varphi }`$, we obtain

$$Z_1=Z_0\frac{\lambda }{T}_0^{\mathrm{}}𝑑x\mathrm{cos}\left(\frac{2hx}{\eta v}\right)\mathrm{exp}\left(\frac{8T}{\eta v}\underset{n}{}_0^{\mathrm{}}𝑑k\frac{e^{|\omega _n|/\mathrm{\Lambda }}\mathrm{sin}^2(kx)}{k^2+(\omega _n/v)^2}\right),$$
(16)
where $`Z_0=Z|_{\lambda =0}`$ and we have introduced the high-energy cutoff $`e^{|\omega _n|/\mathrm{\Lambda }}`$ with $`\mathrm{\Lambda }=v`$ to meet the condition (10). The integral and summation in the exponent can be easily performed, leading to the first-order contribution to the free energy $`F=T\mathrm{ln}Z`$, which is nothing but the leading term in the boundary free energy,

$$F_1=\lambda \left(\frac{\pi T}{v}\right)^{2/\eta }_\alpha ^{\mathrm{}}𝑑x\frac{\mathrm{cos}\left(\frac{2hx}{\eta v}\right)}{\left[\mathrm{sinh}\left(\frac{2\pi Tx}{v}\right)\right]^{2/\eta }},$$
(17)
where $`\alpha `$ is a short-distance cutoff of the order of half a lattice spacing. Note that $`F_1`$ has a contribution only from the boundary region within the distance $`v/2\pi T`$ from the edge of the spin chain. The exponent $`2/\eta `$ is the boundary scaling dimension of the irrelevant operator $`\mathrm{cos}(2\varphi /R)`$.

Our next task is to evaluate the boundary free energy $`F_1`$ at low temperatures. Let us first consider the simplest case $`h=0`$. To obtain the low-temperature expansion of the integral in Eq. (17), we change the variable as

$$I_1=_\alpha ^{\mathrm{}}\frac{dx}{\left[\mathrm{sinh}\left(\frac{2\pi Tx}{v}\right)\right]^{2/\eta }}=\frac{v}{2\pi T}_{u_0}^1𝑑uu^{\frac{2}{\eta }}(1u^2)^{\frac{1}{\eta }1},$$
(18)
where $`u_0=\mathrm{tanh}(2\pi T\alpha /v)2\pi T\alpha /v`$. For the parameter range (6) of our interest, the integral is divergent in the limit $`u_00`$. We integrate $`I_1`$ by parts to separate out the diverging contribution and obtain

$$\frac{2\pi TI_1}{v}=\frac{\eta }{2\eta }u_0^{1\frac{2}{\eta }}(1u_0^2)^{\frac{1}{\eta }1}\frac{1\eta }{2\eta }B(\frac{3}{2}\frac{1}{\eta },\frac{1}{\eta }1),$$
(19)
where we have set $`u_0=0`$ in the second term. The first term gives regular $`T^{2n}`$ terms ($`n=0,1,2,\mathrm{}`$) in the free energy, whereas the cutoff-independent second term leads to the anomalous contribution to the boundary free energy. From Eqs. (17) and (19) we find the leading term in the boundary free energy

$$F_1=E_1\frac{\lambda \eta }{2\sqrt{\pi }(2\eta )}\mathrm{\Gamma }(\frac{1}{\eta })\mathrm{\Gamma }(\frac{3}{2}\frac{1}{\eta })\left(\frac{\pi T}{v}\right)^{\frac{2}{\eta }1},$$
(20)
where $`E_1=F_1(T=0)`$ is the first-order boundary energy and subleading terms of order $`T^2`$ are ignored. Then the boundary contribution to the entropy is

$$S_b=\frac{F_1}{T}=\frac{\sqrt{\pi }}{2}\mathrm{\Gamma }(\frac{1}{\eta })\mathrm{\Gamma }(\frac{3}{2}\frac{1}{\eta })\frac{\lambda }{v}\left(\frac{\pi T}{v}\right)^{\frac{2}{\eta }2},$$
(21)
and the boundary contribution to the specific heat is

$$C_b=T\frac{S_b}{T}=\sqrt{\pi }(\frac{1}{\eta }1)\mathrm{\Gamma }(\frac{1}{\eta })\mathrm{\Gamma }(\frac{3}{2}\frac{1}{\eta })\frac{\lambda }{v}\left(\frac{\pi T}{v}\right)^{\frac{2}{\eta }2}.$$
(22)
For the XXZ chain with $`1/2<\mathrm{\Delta }<1`$, Eqs. (22) and (9) give the exact leading boundary contribution to the specific heat. For $`\mathrm{\Delta }<1/2`$ the leading term in the boundary specific heat is proportional to $`T`$. At $`\mathrm{\Delta }=1/2`$ there appears a logarithmic correction to the $`T`$-linear specific heat.

For a small but finite $`h`$ we may expand $`\mathrm{cos}(2hx/\eta v)`$ in Eq. (17) as

$$_\alpha ^{\mathrm{}}𝑑x\frac{\mathrm{cos}\left(\frac{2hx}{\eta v}\right)1}{\left[\mathrm{sinh}\left(\frac{2\pi Tx}{v}\right)\right]^{2/\eta }}\frac{2h^2}{\eta ^2v^2}_0^{\mathrm{}}\frac{x^2dx}{\left[\mathrm{sinh}\left(\frac{2\pi Tx}{v}\right)\right]^{2/\eta }},$$
(23)
where we have taken the limit $`\alpha 0`$ in the convergent integral. The integral in the right-hand side is elementary. After some algebra we find

$$_0^{\mathrm{}}\frac{x^2dx}{(\mathrm{sinh}x)^{2/\eta }}=\frac{2^{\frac{2}{\eta }4}\eta \mathrm{\Gamma }(\frac{1}{\eta })\mathrm{\Gamma }(3\frac{2}{\eta })}{(2\eta )\mathrm{\Gamma }(2\frac{1}{\eta })}[\pi ^22\psi ^{}(\frac{1}{\eta })],$$
(24)
where $`\psi ^{}(x)=d^2\mathrm{ln}\mathrm{\Gamma }(x)/dx^2`$. The boundary contribution to the susceptibility is then obtain as

$`\chi _b`$ $`=`$ $`{\displaystyle \frac{^2F_1}{h^2}}|_{h=0}`$ (25)
$`=`$ $`{\displaystyle \frac{\lambda \mathrm{\Gamma }(\frac{1}{\eta })\mathrm{\Gamma }(3\frac{2}{\eta })[\pi ^22\psi ^{}(\frac{1}{\eta })]}{4v^2\eta (2\eta )\mathrm{\Gamma }(2\frac{1}{\eta })}}\left({\displaystyle \frac{2\pi T}{v}}\right)^{\frac{2}{\eta }3}.`$
This is the exact leading contribution to the boundary susceptibility in the XXZ chain with $`1/2<\mathrm{\Delta }<1`$. Note that $`\chi `$ is diverging as $`T0`$. For $`\mathrm{\Delta }<1/2`$ the leading term in the susceptibility is independent of temperature. A logarithmic correction shows up at $`\mathrm{\Delta }=1/2`$.

At $`T=0`$ the boundary energy $`E_1`$ becomes

$`E_1`$ $`=`$ $`\lambda {\displaystyle _\alpha ^{\mathrm{}}}𝑑x{\displaystyle \frac{\mathrm{cos}\left(\frac{2hx}{\eta v}\right)}{(2x)^{2/\eta }}}`$ (26)
$`=`$ $`{\displaystyle \frac{\lambda (2\alpha )^{1\frac{2}{\eta }}}{2(\frac{2}{\eta }1)}}\mathrm{cos}\left({\displaystyle \frac{2h\alpha }{\eta v}}\right)`$
$`{\displaystyle \frac{\lambda \mathrm{\Gamma }(3\frac{2}{\eta })\mathrm{sin}[\pi (\frac{1}{\eta }1)]}{4(\frac{2}{\eta }1)(\frac{1}{\eta }1)}}\left({\displaystyle \frac{h}{\eta v}}\right)^{\frac{2}{\eta }1},`$
from which we obtain the boundary contribution to the total magnetization

$$M_b=\frac{E_1}{h}=\frac{\pi \lambda }{4\eta v}\frac{\mathrm{\Gamma }(3\frac{2}{\eta })}{\mathrm{\Gamma }(\frac{1}{\eta })\mathrm{\Gamma }(2\frac{1}{\eta })}\left(\frac{h}{\eta v}\right)^{\frac{2}{\eta }2},$$
(27)
where higher-order terms $`𝒪(h)`$ are ignored. This is the exact form of the leading term in the XXZ chain for $`1/2<\mathrm{\Delta }<1`$. The leading terms in $`C_b(T)`$, $`\chi _b(T)`$, and $`M_b(h)`$ we calculated above are found in agreement with those obtained independently by Fujimoto and Eggert,FujimotoEggert  after some analytical transformations are made in their results.

At the isotropic antiferromagnetic Heisenberg point ($`\mathrm{\Delta }=1`$), we cannot substitute $`\lambda `$ of Eq. (9) into Eqs. (22), (25), and (27), as it vanishes at $`\eta =1`$. To circumvent this difficulty, we adapt a RG improved perturbation theory in which we replace a constant $`\lambda `$ with a running coupling constant $`\lambda (E)`$ at energy scale $`E`$. At $`h=0`$, the one-loop scaling equations for the sine-Gordon Hamiltonian (3) are given by

$$\frac{dg}{dl}=\left(2\frac{1}{\pi R^2}\right)g,\frac{dR}{dl}=\frac{\pi g^2}{2R},$$
(28)
where $`g=\lambda /v`$ and $`dl=d\mathrm{ln}E`$. The coupling $`g`$ is marginally irrelevant at $`R=1/\sqrt{2\pi }`$ (i.e., at $`\mathrm{\Delta }=1`$) and decreases to zero logarithmically. To see this explicitly, we may write $`\sqrt{2\pi }R=1ϵ`$ and expand the RG equations (28) in lowest order in $`ϵ`$. We solve the simplified RG equations to find the solution corresponding to the antiferromagnetic Heisenberg chain,

$$\frac{\lambda (l)}{v}=\frac{1}{2\pi l},R(l)=\frac{1}{\sqrt{2\pi }}\left(1\frac{1}{4l}\right)$$
(29)
for $`l1`$. To obtain the leading term in $`C_b(T)`$, $`\chi _b(T)`$, and $`M_b(h)`$, we may substitute $`\eta =1`$ and

$$\lambda (T,h)=\frac{v}{2\pi \mathrm{ln}\left({\displaystyle \frac{1}{\mathrm{max}(T,h)}}\right)}$$
(30)
into the analytic formulas of the boundary free energy we obtained for $`1/2<\mathrm{\Delta }<1`$. For example, the boundary entropy $`S_b`$ is found from Eq. (21) to be

$$S_b=\frac{\pi \lambda (T)}{2v}=\frac{1}{4\mathrm{ln}(T_0/T)},$$
(31)
from which the boundary specific heat becomes

$$C_b=\frac{1}{4[\mathrm{ln}(T_0/T)]^2},$$
(32)
in agreement with Ref. FujimotoEggert, . Here $`T_0`$ is a constant of order $`J`$. The boundary contribution to the susceptibility is obtained from Eq. (25) as

$$\chi _b=\frac{\lambda (T)}{8\pi vT}[\pi ^22\psi ^{}(1)]=\frac{1}{24T\mathrm{ln}(T_0/T)}.$$
(33)
This result was confirmed later by the corrected calculation in Ref. FujimotoEggert, . The boundary contribution to the magnetization at zero temperature is found from Eq. (27) to be

$$M_b=\frac{\pi \lambda (h)}{4v}=\frac{1}{8\mathrm{ln}(h_0/h)},$$
(34)
in agreement with the exact result from the Bethe ansatz.Essler ; Asakawa ; Zvyagin  Here $`h_0`$ is another constant of order $`J`$. We note that there exist subleading terms to these formulas which are reduced only by a factor of $`1/\mathrm{ln}(T_0/T)`$ or $`1/\mathrm{ln}(h_0/h)`$. Even though they are important for quantitative analysis, we do not try to evaluate it here, as it is beyond the scope of this paper.

Finally, we present comparison between the analytic result of Eq. (26) and numerical data of energy eigenvalues of the XXZ spin chains with open boundaries. Using the density-matrix renormalization-group (DMRG) method,White1 ; White2  we have computed the lowest energy $`E(L,m)`$ of $`H_0`$ in the subspace of the magnetization per site $`m=(1/L)_{j=1}^LS_j^z`$. The maximum system size we calculated is $`L=800`$. The number of block states kept in the DMRG calculation is up to 200. The numerical errors due to the DMRG truncation, estimated from the difference between the data with 150 and 200 states kept, are typically less than $`10^5`$. The numerical data are therefore accurate enough for the discussion below.

The energy $`E(L,m)`$ is known to depend on the system size $`L`$ as

$$E(L,m)=L\epsilon _0(m)+\epsilon _1(m)+\frac{\epsilon _2(m)}{L+1}.$$
(35)
The bulk energy $`\epsilon _0(m)`$ can be obtained exactly by solving Bethe ansatz integral equations,BogoIK ; QinFYOA ; CabraHP  while the boundary energy $`\epsilon _1(m)`$ is related to Eq. (26). Here, we fit the numerical data to the formula (35) taking $`\epsilon _0`$, $`\epsilon _1`$, and $`\epsilon _2`$ as free parameters. The data of $`E(L,m)`$ at $`\mathrm{\Delta }=0.8`$ are shown in Fig. 1 for several values of $`m`$. We have found that the fitting works pretty well for all $`\mathrm{\Delta }`$s and $`m`$s calculated. In fact, the estimated values of $`\epsilon _0`$ obtained from the fitting are in excellent agreement with the exact values calculated from the Bethe ansatz equations (see Fig. 2). The bulk energy $`\epsilon _0(m)`$ has smooth quadratic dependence on $`m`$ since the irrelevant operator $`\mathrm{cos}(2\pi /R)`$ yields only subleading contributions to the bulk quantity.

According to Eq. (26) and the relation between the bulk magnetization and the external field ($`m=\chi _0h/L`$), the boundary energy $`\epsilon _1(m)`$ should depend on $`m`$ as

$$\epsilon _1(m)=\epsilon _1(0)h_1(2\pi m)^{\frac{2}{\eta }1}\frac{(2\pi m)^2}{2\chi _2}$$
(36)
for small $`m`$, where $`\chi _2`$ is a constant and the exact form of $`h_1`$ is

$$h_1=\frac{\lambda \mathrm{\Gamma }(3\frac{2}{\eta })\mathrm{sin}[\pi (\frac{1}{\eta }1)]}{2(\frac{2}{\eta }1)(\frac{1}{\eta }1)}.$$
(37)
Note that the boundary energy of the open spin chains should be twice that of the semi-infinite wire that we have considered in the analytic calculation. To estimate the coefficient $`h_1`$ numerically, we fit the data of $`\epsilon _1(m)`$ to Eq. (36) taking $`\epsilon _1(0)`$, $`h_1`$, and $`\chi _2`$ as free parameters. We show the $`m`$ dependence of $`\epsilon _1(m)`$ and the fitting curves in Fig. 3. Clearly, the plot of $`\epsilon _1(m)`$ versus $`(2\pi m)^2`$ shows large deviation from linear dependence on $`(2\pi m)^2`$, confirming the existence of the leading power-law term with the exponent $`2/\eta 1`$, which emerges from the irrelevant operator $`\lambda \mathrm{cos}(2\varphi /R)`$. The values of $`h_1`$ estimated from the fitting procedure are also shown in Fig. 3. We see that the numerical estimates agree extremely well with the analytic result (37). This agreement can be regarded as a numerical proof of our analytic result, Eq. (26).

In summary, we have obtained the exact leading boundary contributions to the specific heat and the susceptibility using the low-energy effective theory for the spin-1/2 XXZ model. Our method can be easily extended to other 1D systems (e.g., the extended Hubbard model and spin ladders in high magnetic fields) as long as their low-energy effective theory is given by Eq. (3), which contains nonuniversal constants $`v`$ and $`\lambda `$ to be determined with some other method, however. Our formulas (32), (33), and (34) for the spin isotropic case ($`\mathrm{\Delta }=1`$) are universal, in the sense that if we make any moderate deformation of the isotropic antiferromagnetic Heisenberg model that respects the SU(2) and translational symmetry and that does not generate a spin gap, then these formulas apply to some values of $`T_0`$ and $`h_0`$, except for cases where the coupling $`\lambda `$ is fine tuned to vanish. This implies that, for example, the 1D (extended) Hubbard model and other systems of interacting electrons with the spin SU(2) symmetry should universally exhibit the boundary contributions (32), (33), and (34).

###### Acknowledgements.

We are grateful to S. Fujimoto for sharing his results prior to publication. This work was supported in part by a Grant-in-Aid for Scientific Research on Priority Areas (Grant No. 12046238) and by NAREGI Nanoscience Project, Ministry of Education, Culture, Sports, Science and Technology, Japan.


format: text

text:
# Excitonic Effects and Optical Spectra of Single-Walled Carbon Nanotubes

## Abstract

Many-electron effects often dramatically modify the properties of reduced dimensional systems. We report calculations, based on an ab initio many-electron Green’s function approach, of electron-hole interaction effects on the optical spectra of small-diameter single-walled carbon nanotubes. Excitonic effects qualitatively alter the optical spectra of both semiconducting and metallic tubes. Excitons are bound by $`1`$ eV in the semiconducting (8,0) tube and by $`100`$ meV in the metallic (3,3) tube. These large many-electron effects explain the discrepancies between previous theories and experiments.

Synthesis and observation of single-walled carbon nanotubes (SWCNT) have advanced greatly in recent years, making possible the experimental study of the optical properties of individual SWCNTs ZMLi ; OConnell . If well understood, the optical response of SWCNTs may be used to characterize these nanotubes, to monitor and guide their separation by type Bachilo , and can be employed in device applications Misewich . However, measured optical transition frequencies deviate substantially from theoretical predictions based on one-particle interband theories. This deviation is not unexpected since many-body interactions should play a vital role in reduced dimensions Ando . Our ab initio results show that, indeed, many-electron effects can change qualitatively the optical spectra of SWCNTs. Strongly bound exictons are predicted in small diameter semiconducting nanotubes and even in some metallic tubes, and they dominate the optical response.

Below, motivated by recent experiments ZMLi ; Bachilo , we compute the optical absorption spectra of the three small-diameter SWCNTs: (3,3), (5,0), and (8,0). We use a recently developed approach in which electron-hole excitations and optical spectra of real materials are calculated from first principles in three stages Rohlfing : (i) we treat the electronic ground-state with ab initio pseudopotential density-functional theory (DFT) Sham , (ii) we obtain the quasiparticle energies $`E_{n𝐤}`$ within the GW approximation for the electron self-energy $`\mathrm{\Sigma }`$ Hybertsen  by solving the Dyson equation:

$$\left[\frac{^2}{2}+V_{ion}+V_{Hartree}+\mathrm{\Sigma }(E_{n𝐤})\right]\psi _{n𝐤}=E_{n𝐤}\psi _{n𝐤},$$

and (iii) we calculate the coupled electron-hole excitation energies $`\mathrm{\Omega }^S`$ and spectrum by solving the Bethe-Salpeter equation of the two-particle Green’s function Rohlfing ; Strinati :

$$\left(E_{c𝐤}E_{v𝐤}\right)A_{vc𝐤}^S+\underset{𝐤^{}v^{}c^{}}{}vc𝐤|K^{eh}|v^{}c^{}𝐤^{}A_{v^{}c^{}𝐤^{}}^S=\mathrm{\Omega }^SA_{vc𝐤}^S,$$

where $`A_{vc𝐤}^S`$ is the exciton amplitude, $`K^{eh}`$ is the electron-hole interaction kernel, and $`|c𝐤`$ and $`|v𝐤`$ are the quasielectron and quasihole states, respectively. We obtain the DFT wavefunctions and eigenvalues by solving the Kohn-Sham equations within the local density approximation (LDA) Sham  using a plane-wave basis with an energy cutoff of 60 Ry. We use ab initio Troullier-Martins pseudopotentials Martins  in the Kleinmann-Bylander form KB  ($`r_c=1.4`$ a.u.). To compare with experiments in which 4 Å diameter SWCNTs are grown inside zeolites ZMLi , we study the (3,3) and (5,0) tubes in the experimental geometry with a dielectric background of $`AlPO_4`$ dielectric . For the (8,0) tube, we work in a supercell with an intertube separation of at least $`9.7`$ Å to mimic experiments on isolated tubes OConnell ; Bachilo . In supercells, due to the long range of the screened Coulomb interaction in semiconduting tubes, unphysical interactions between periodic images can lead to deviations from the isolated case. Hence, we truncate the Coulomb interaction in a cylindrical geometry for the semiconducting tubes (we find negligible tube-tube interactions or image effects for metallic tubes where screening is complete). Due to depolarization effects in nanotubes Ajiki , strong optical response is observed only for light polarized along the tube axis ($`\widehat{z}`$). We only consider this polarization below.

For the metallic tubes (3,3) and (5,0), we find quasiparticle GW corrections to the LDA band energies similar to those in graphite, namely a $`15\%`$ stretching of the LDA eigenvalues away from the Fermi level ($`E_F`$) Fong . While this result is expected for large diameter metallic nanotubes which resemble a graphene sheet, we find it also holds for these small diameter metallic tubes where curvature effects lead to strong $`\sigma \pi `$ hybridization Blase .

Fig. 1a shows the quasiparticle density of states (DOS) for the metallic (3,3) tube featuring a number of prominent one-dimensional (1D) van Hove singularities (vHs) near $`E_F`$. Unlike predictions from simple tight-binding models Dress , these vHs are asymmetric about $`E_F`$ due to strong curvature effects. The arrow in the figure indicates optically allowed low-energy transitions. For the (3,3) metallic tube, the bands forming the first vHs below $`E_F`$ and the second vHs above $`E_F`$ meet at the Fermi level, but optical transitions between them are symmetry-forbidden. We calculate $`ϵ_2(\omega )`$ in two ways (see Fig. 1b). First, we neglect electron-hole interactions and find the existence of a symmetry gap, i.e. no electron-hole transition with energy below the prominent peak at $`\mathrm{}\omega =3.25`$ eV. Second, we include electron-hole interactions and solve the Bethe-Salpeter equation. In general, the electron-hole interaction kernel has two terms: an attractive direct term involving the screened Coulomb interaction and a repulsive exchange term involving the bare Coulomb interaction Rohlfing . Fig. 1b shows that for the (3,3) tube the direct term dominates: one bound exciton is formed with a binding energy of $`86`$ meV and an extent of $`50`$ Å along $`\widehat{z}`$. The surprising result of having a bound exciton in a metal stems from having the symmetry gap, which is possible for a quasi 1-D system where all k-states have well-defined symmetry. Also, the existence of a sole bound exciton is due to the metallic screening (the screening length is $``$ 3.2 Å): the effective electron-hole interaction along $`\widehat{z}`$ resembles an attractive $`\delta (z)`$ function, and in 1D, the Hamiltonian $`H=\frac{1}{2m^{}}\frac{d^2}{dz^2}|V_0|\delta (z)`$ has a single bound eigenstate.

Fig. 2a shows the quasiparticle band structure for the metallic (5,0) tube. According to the band-folding scheme Dress ; Hamada , this tube should be semiconducting. However, curvature effects Blase  lead to strong $`\sigma \pi `$ hybridization, forcing a band to cross $`E_F`$. Arrows in the figure indicate optically allowed interband transitions which give rise to the two peaks, labeled A and B in the optical spectrum in Fig. 2b. When neglecting electron-hole interactions, peak B has a lower intensity than A because the transitions contributing to B do not originate from a band extremum (vHs) but from the crossing at $`E_F`$. For this tube, electron-hole interactions do not bind excitons: while the screening length in the (5,0) tube is similar to that of the (3,3) tube, the symmetry of the bands in the (5,0) tube prohibits direct attraction between the electron-hole pairs contributing to peaks A and B. Thus the electron-hole interaction is governed by the repulsive exchange term. This effect, again, is peculiar to nanotubes: in traditional semiconductors, the attractive direct term dominates over the repulsive exchange term. Moreover, when electron-hole interactions are included, the exchange term has a larger effect on peak B and suppresses it greatly.

We now compare our results for the (3,3) and (5,0) tubes to experiments. In the work of Li et al. ZMLi , 4 Å diameter SWCNTs were grown inside zeolite channels, and three main peaks were found in the measured absorption spectra (see Table 1). While 4 Å diameter SWCNTs come in only three chiralities, (3,3), (5,0) and (4,2), it was not possible to assign directly the specific peaks to specific tubes experimentally. As shown in Table 1, our results for the (3,3) and (5,0) tubes are in excellent quantitative agreement with experiment and provide a concrete identification for two of the observed peaks. We conclude that the remaining peak at 2.1 eV is due to the (4,2) tube (other calculations ignoring electron-hole interactions point to the same conclusion CTChan ; Ordejon ; 42 ). Moreover, the many-electron suppression of peak B in the (5,0) spectrum explains the absence of any observed feature in the measured spectra at $`2.8`$ eV.

We now consider the 6.3 Å diameter semiconducting tube (8,0), in which we expect even larger excitonic effects. The (8,0) tube has a calculated LDA minimum band gap of $`0.60`$ eV. Quasiparticle corrections dramatically open the gap to $`1.75`$ eV. This correction is significantly larger than those in bulk semiconductors with similar LDA gaps: we again attribute this to the 1D nature of the SWCNTs which enhances Coulomb effects (as shown in model calculations Ando ).

Fig. 3a shows the calculated absorption spectrum for an isolated (8,0) tube. There are three distinct low-energy peaks (labeled A, B, and C) in the non-interacting spectrum. When electron-hole interactions are included, we find far more dramatic excitonic effects than in the metallic cases: each non-interacting peak gives rise to a series of visible exciton lines with large binding energies of $`0.99`$ eV, $`0.86`$ eV and $`1.00`$ eV for the lowest-energy excitons ($`A_1^{}`$, $`B_1^{}`$ and $`C_1^{}`$ respectively). These binding energies are more than ten times larger than those in bulk semiconductors with similar gaps, and excitonic effects qualitatively change the spectral profile. Again, these effects stem from the long range of the screened Coulomb interaction and the 1D nature of the SWCNTs: e.g., the binding energy of a 1D hydrogenic system is infinite due to the long-range Coulomb interaction Loudon . We note that the electron-hole interaction reverses the relative intensity of the first and second prominent optical peaks.

Theory predicts that there are two varieties of excitons in the (8,0) tube: bound excitons with energies below the non-interacting optical gap (A’ and B’ series) and resonant excitons with energies above the non-interacting optical gap (C’ series). Fig. 3b shows the real-space, electron-hole pair probability distribution $`|\mathrm{\Phi }(𝐫_e,𝐫_h)|^2`$ as a function of the electron position $`𝐫_e`$ for the photoexcited $`A_1^{}`$ bound exciton obtained by fixing the position of the hole $`𝐫_h`$ (the black star in the figure) on a carbon $`\pi `$ orbital. Fig. 3c and 3d show this correlation more quantitatively for the bound $`A_1^{}`$ and the resonant $`C_1^{}`$ exciton: $`|\mathrm{\Phi }|^2`$ is plotted along $`\widehat{z}`$ after integrating out the electron coordinates in the perpendicular plane. The extent of the bound part of both excitons is $``$ 25 Å.

Our results for the (8,0) tube are in excellent agreement with the experimental findings of Bachilo et al. Bachilo ; Weisman : by performing spectrofluorimetric measurements on various semiconducting SWCNTs, with diameters ranging from 0.62 to 1.31 nm and chiral angle from 3 to 28 degrees, optical transitions were assigned to specific individual (n,m) tubes. While their SWCNT samples did not contain the (8,0) tube, they obtained results for tubes with similar diameter and chirality and provided fits for the first and second optical transition energies ($`\nu _{11}`$ and $`\nu _{22}`$) as a function of diameter and chiral angle, which they demonstrate to work well for a wide range of (n,m) values. Their fits provide a ratio of $`\nu _{22}/\nu _{11}=1.17`$ for the (8,0) tube. The traditional non-interacting $`\pi `$-orbital only tight-binding model gives a ratio of 2, and the experimental deviation from 2 has been a puzzle Ichida ; Liu ; Mele . However, as shown in Table 2, our results for $`\nu _{11}`$ and $`\nu _{22}`$ in the (8,0) tube (peaks $`A_1^{}`$ and $`B_1^{}`$) and their ratio are in excellent agreement with the deduced experimental values. The deviation of $`\nu _{22}/\nu _{11}`$ from 2 is a consequence of both band-structure and many-electron effects: one needs to include both for a basic and quantitative understanding.

In conclusion, we study the optical absorption spectra of metallic and semiconducting small diameter SWCNTs and obtain excellent agreement with available experimental data. We show that electron-hole interactions (which can be either attractive or repulsive) play a crucial role, especially for semiconducting tubes, in understanding experimental results. Large excitonic features for both semiconducting and metallic tubes are seen to be due to the quasi-1D nature of SWCNTs, and the manner in which they effect the spectra depends strongly on the rotational symmetries of the tubes.

This work was supported by the NSF under Grant #DMR0087088, and the Office of Energy Research, Office of Basic Energy Sciences, Materials Sciences Division of the U.S. Department of Energy (DOE) under Contract #DE-AC03-76SF00098. Computer time was provided at the DOE Lawrence Berkeley National Laboratory (LBNL)’s NERSC center. Portions of this work were performed under the auspices of the DOE by the University of California Lawrence Livermore National Laboratory (LLNL) under contract No. W-7405-Eng-48. Collaborations between LLNL and LBNL were facilitated by the DOE Computational Materials Sciences Network.


text:
# Properties of galaxy dark matter halos from weak lensing

## 1. Introduction

The existence of massive dark matter halos around galaxies is widely accepted, based on different lines of evidence, such as flat rotation curves of spiral galaxies and strong lensing systems. However, relatively little is known about the properties of dark matter halos, such as their extent and shapes.

A promising approach to study the galaxy dark matter halos is weak gravitational lensing. The tidal gravitational field of the dark matter halo introduces small coherent distortions in the images of distant background galaxies. The weak lensing signal can be measured out to large projected distances from the lens, and hence provides a unique probe of the gravitational potential on large scales. The applications of this approach are numerous: one can infer masses of galaxies and compare the results to their luminosities (e.g., McKay et al. 2001), or one can attempt to constrain the halo mass profile (e.g., Fischer et al. 2000; McKay et al. 2001). Also, weak lensing can be used to constrain the shapes of halos by measuring the azimuthal variation of the lensing signal.

Unfortunately, one can only study ensemble averaged properties, because the weak lensing signal induced by an individual galaxy is too low to be detected. Nevertheless, as we will show here, weak lensing is a useful probe of the matter distribution in galaxies, and we expect it to make significant contributions to our understanding of galaxy formation in the coming years.

In these proceedings we highlight the main results from our weak lensing analysis. A detailed description of our results can be found in Hoekstra, Yee, & Gladders (2003).

## 2. Data & Analysis

We use 45.5 deg2 of $`R_C`$-band imaging data from the Red-Sequence Cluster Survey (e.g., Yee & Gladders 2001), which were taken with the CFH12k camera on the CFHT. A detailed description of the data reduction and object analysis can be found in Hoekstra et al. (2002), to which we refer for technical details.

For the analysis presented here, we select a sample of “lenses” and “sources” on the basis of their apparent $`R_C`$ magnitude. We define galaxies with $`19.5<R_C<21`$ as lenses, and galaxies with $`21.5<R_C<24`$ as sources which are used to measure the lensing signal. This selection yields a sample of $`1.2\times 10^5`$ lenses and $`1.5\times 10^6`$ sources.

The interpretation of the weak lensing signal (e.g., inferring sizes and masses for the galaxy halos), requires knowledge of the redshift distributions of both lenses and sources. The redshift distribution of the sample of lenses has been determined spectroscopically by the CNOC2 Field Galaxy Redshift Survey (e.g., Yee et al. 2000). The derived redshift distribution gives a median redshift $`z=0.35`$ for the lens sample. We use the redshifts and the colors of the galaxies observed in the CNOC2 survey to compute their rest-frame $`B`$ luminosity. For the source galaxies the situation is more complicated. These galaxies are generally too faint for spectroscopic surveys. Instead we use the photometric redshift distributions derived from both Hubble Deep Fields, which yields a median redshift of $`z=0.53`$ for the source galaxies.

Weak lensing measures the convolution of the galaxy distribution and the galaxy dark matter profiles. To examine the ensemble average properties of the dark matter halos properly, we need to account for the clustering of the lenses. This is done naturally in a maximum likelihood analysis, where a model for the mass distribution of individual galaxies is compared to the observations. We use the profile predicted by cold dark matter simulations (e.g., Navarro, Frenk, & White 1995, NFW hereafter). The NFW density profile is characterized by two parameters, a density contrast $`\delta _c`$ and a scale $`r_s`$

$$\rho (r)=\frac{\delta _c\rho _c}{(r/r_s)(1+r/r_s)^2},$$
(1)
where $`\rho _c`$ is the critical surface density at the redshift of the halo. The “virial” radius $`r_{200}`$ is defined as the radius where the mass density of the halo is equal to $`200\rho _c`$. The virial mass $`M_{200}`$ is defined as the mass enclosed within $`r_{200}`$, with a corresponding rotation velocity $`V_{200}=V_c(r_{200})`$.

In these proceedings we compare the NFW profile to the observations, with $`V_{200}`$ (or equivalently $`M_{200}`$) and $`r_s`$ as free parameters. Hoekstra et al. (2003) also consider a truncated isothermal sphere (TIS) model.

To infer the best estimates for the model parameters, one formally has to perform a maximum likelihood analysis in which the redshift of each individual galaxy is a free parameter, which has to be chosen such that it maximizes the likelihood. This approach is computationally not feasible, and instead we create mock redshift catalogs, using the observed redshift distributions from the CNOC2 survey (see Hoekstra et al. 2003 for details), which allows us to find estimates for the model parameters, close to the true maximum likelihood values.

## 3. Scaling relations

We first examine how the model parameters $`\delta _c`$ and $`r_s`$ scale with luminosity. We assume that the scaling relations have a power law dependence on the restframe $`B`$-band luminosity:

$$\delta _cL_\mathrm{B}^\beta ,\mathrm{and}r_sL_\mathrm{B}^{\eta /2},$$
(2)
Currently we use the apparent $`R_C`$-band magnitude as a crude measure of the luminosity. With upcoming multi-color data, we expect to significantly improve on our measurements. To derive joint constraints on $`\beta `$ and $`\eta `$, we marginalise over the values for $`\delta _c`$ and $`r_s`$. It is more convenient to constrain $`\beta +\eta `$ and $`\eta `$, and the results are presented in Figure 1.

Figure 1 shows that the sum $`\beta +\eta `$ is fairly well constrained. We obtain a best fit value of $`\beta +\eta =0.40_{0.13}^{+0.14}`$ (68% confidence). For $`\eta `$ we find $`\eta =0.77_{0.67}^{+0.63}`$. Also the relation between the maximum rotation velocity $`V_{\mathrm{max}}`$ and the luminosity is well described by a power law, with an exponent $`0.20_{0.14}^{+0.13}`$ (68% confidence) in good agreement with the observed slope of the Tully-Fisher relation (e.g., Verheijen 2001).

## 4. Mass and extent of halos

Many different lines of evidence suggest that galaxies are surrounded by massive dark matter halos, but observationally it is difficult to place constraints on the mass and extent of these halos, because of the lack of visible tracers that can be used to infer the gravitational potential. Weak gravitational lensing does not suffer from the latter requirement and consequently provides one of the most powerful measures of the mass distribution at large radii. With the current data it is not possible to distinguish between different mass profiles (as mentioned above we use the NFW profile), but for a given model, we can examine the masses and sizes of the halos.

Figure 2 shows the joint constraints on $`V_{200}`$ (the corresponding values for $`M_{200}`$ are indicated on the right) and $`r_s`$ for a galaxy with a luminosity of $`L_B=10^{10}h^2\mathrm{L}_\mathrm{B}`$(under the assumption that the luminosity does not evolve with redshift). The mass is well constrained and we find a best fit value of $`V_{200}=166\pm 5`$ km/s or $`M_{200}=(8.8\pm 0.7)\times 10^{11}h^1M_{}`$ (68% confidence). For the scale $`r_s`$ we find $`r_s=16.7_{3.0}^{+3.7}h^1`$ kpc (68% confidence).

In our maximum likelihood analysis we considered $`r_s`$ and $`V_{200}`$ free parameters. Numerical simulations, however, show that the parameters in the NFW model are correlated, albeit with some scatter. As a result, the NFW model can be considered as a one-parameter model. The dotted line in Figure 2 indicates this prediction. If the simulations provide a good description of dark matter halos, the dotted line should intersect our confidence region, which it does.

This result provides important support for the CDM paradigm, as the latter predicts a “size” of dark matter halos which is in excellent agreement with our observations. It is important to note that this analysis is a direct test of CDM (albeit not conclusive), because the weak lensing results are inferred from the gravitational potential at large distances from the galaxy center, where dark matter dominates. Most other attempts to test CDM are confined to the inner regions, where baryons are, or might be, important.

## 5. Shapes of halos

The average shape of dark matter halos can provide important information about the nature of dark matter. Numerical simulations of cold dark matter yield triaxial halos, with a typical ellipticity of $`0.3`$ (e.g., Dubinski & Carlberg 1991). Hence, in the context of collisionless cold dark matter, the theoretical evidence for flattened halos is quite strong. If the dark matter is interacting, it tends to produce halos that are more spherical (compared to cold dark matter). This difference is more pronounced in the central parts of the halo, where the density is high. On the large scales probed by weak lensing, the different types of dark matter (for reasonable interaction cross-sections) produce halos with similar shapes.

Weak gravitational lensing is potentially the most powerful way to derive constraints on the shapes of dark matter halos. The amount of data required for such a measurement, however, is large: the galaxy-galaxy lensing signal is small, and now one needs to measure an even smaller azimuthal variation. We also have to assume that halo is aligned with the galaxy. An imperfect alignment between light and halo will reduce the amplitude of the azimuthal variation detectable in the weak lensing analysis. Hence, weak lensing formally provides a lower limit to the average halo ellipticity.

To maximize the signal-to-noise ratio of the measurement one has to assign proper weights to the lenses: edge-on galaxies have maximal weight, whereas the lensing signal around face-on galaxies contains no information about the shape of the halo. We adopt a simple approach, and assume that the (projected) ellipticity of the dark matter halo is proportional to the shape of the galaxy: $`e_{\mathrm{halo}}=fe_{\mathrm{lens}}`$.

As before, we compute the model shear field, and compare this to the data. Figure 3 shows the resulting $`\mathrm{\Delta }\chi ^2`$ as a function of $`f`$. We find a best fit value of $`f=0.77_{0.21}^{+0.18}`$ (68% confidence). This suggests that, on average, the dark matter distribution is rounder than the light distribution. As discussed above, our analysis formally provides only a lower limit on the halo ellipticity, and the true ellipticity might be higher if some of the halos are misaligned with the light. Nevertheless, the fact that we detect a significant flattening implies that the halos are well aligned with the light distribution. Also note that the lensing signal is caused by a range of different galaxy types, for which our simple relation between the halo ellipticity and light distribution might not be valid.

A simple interpretation of the results is difficult, but a simple approach actually yields sensible results. For instance, the average ellipticity of the lens galaxies is $`e_{\mathrm{lens}}=0.261`$. Hence, the measured value of $`f`$ implies an average projected halo ellipticity of $`e_{\mathrm{halo}}=0.20_{0.05}^{+0.04}`$ (68% confidence), which corresponds to an projected axis ratio of $`c/a=0.66_{0.06}^{+0.07}`$ (68% confidence). Although the weak lensing yields a projected axis ratio, the result is in fair agreement with the results from numerical simulations.

A robust outcome of our analysis is that spherical halos ($`f=0`$) are excluded with 99.5% confidence. This result poses a serious problem for alternative theories of gravity, which attempt to explain the observations without the need of dark matter. In such theories the lensing signal we measure is effectively caused by point masses (the visible matter is confined to much smaller scales). As a result, in such theories one would expect an almost isotropic lensing signal around the galaxies (even if the quadrupole decreases as slow as $`r^\sqrt{3}`$, e.g., Binney this proceedings), which is not observed. Note that this test does not require knowledge of the actual deflection angles, which typically are unknown for these theories. Hence, our results provide strong support for the existence of dark matter.

## References

Dubinski, J. & Carlberg, R.G. 1991, ApJ, 378, 496

Fischer, P., et al. 2000, AJ, 120, 1198

Hoekstra, H., Yee, H.K.C., Gladders, M.D., Barrientos, L.F., Hall, P.B., & Infante, L. 2002, ApJ, 572, 55

Hoekstra, H., Yee, H.K.C, & Gladders 2003, ApJ, submitted, astro-ph/0306515

McKay, T.A., et al. 2001, ApJ, submitted, astro-ph/0108013

Navarro, J.F., Frenk, C.S., & White, S.D.M. 1995, MNRAS, 275, 56

Verheijen, M.A.W. 2001, ApJ, 563, 694

Yee, H.K.C., et al. 2000, ApJS, 129, 475

Yee, H.K.C., & Gladders, M.D. 2001, in “AMiBA 2001: High-z Clusters, Missing Baryons, and CMB Polarization”, ASP Conference Series, Eds. L.-W. Chen et al., astro-ph/0111431


text:
# A new sample of young Planetary Nebulae

## 1. Introduction

During the last few years many observational programs have been devoted to recognize new planetary and proto-planetary nebulae. Such studies were aimed to understand the process of formation of PNs by discovering new transition objects in the very short phase between the end of the AGB and the onset of the ionization.

A small sub-group of B-type stars, called BQ\[\] stars, defined as $`B_\mathrm{e}`$ with forbidden emission lines, were recognized as potential candidates to be new transition objects (Parthasarathy & Pottasch, 1989) on the basis of their IR excess. BQ\[\] stars, however, are not a well defined group, and there is still a controversy on their evolutionary stage.

## 2. Observations and Results

In order to find new very young Planetary Nebulae we have selected a sample of 16 hot post-AGB stars (BQ\[\]) from the most recent compilations, namely Parthasarathy, (1993); Conlon et al., (1993) and Parthasarathy et al., (2000). All the selected candidates have high galactic latitude, infrared excess, spectral type B1 I-II and the presence of nebular emission lines in their spectrum. In particular the last two requirements maximize the possibility to detect a radio nebula. The detection of radio nebulae associated to the selected targets would definitively assess the evolutionary status of this kind of objects, producing a unique sample of very young PNs.

The observations were carried out by using the VLA in two different runs, on June 8 and 10, 2001, at 8.4 GHz in CnB configuration.
We detected a total of 10 sources over 16 been observed. Our results are summarized in Table 1, where the radio flux density, with associated error ($`\sigma `$) and the rms of the map, are reported. The error associated to the flux density estimation ($`S`$) is a combination of the rms in the map plus the amplitude calibration $`error,\sigma _{cal}`$, $`\sigma =\sqrt{(rms)^2+(\sigma _{cal}S)^2}`$.

## 3. Discussion and Conclusions

In order to calculate distances to all the investigated sources, we have chosen to apply, among different methods, that developed by Tajitsu & Tamura (1998) which takes advantage of one characteristic common to all the sources in our sample, namely the presence of a strong far-infrared flux as detected by IRAS satellite and gives as by-product the total far-IR flux ($`F_{IR}`$).
No systematic difference in distances between detected and non-detected sources have been found. Therefore the non-detections are not related to larger distances but should be related to intrinsic characteristic of the source, such as the evolution of the ionization structure.

In order to check if our sample consists of young PNs we calculate some physical quantities, whose values can help in understanding the evolutionary stage of the nebula. Those are (Table 2):

* the brightness temperature, that for young nebulae should be of the order of $`T_B10^3K`$
* the Emission measure, that for young nebulae should be of the order of $`EM10^8cm^6pc`$
* the Infrared Excess, that for young nebulae should be $`IRE1`$

Our results have been compared to those obtained by Aaquist and Kwok (1991, AL91), who observed with the VLA at 15 GHz a sample of yPNs.
For the AK91 sample $`T_B`$ and $`EM`$ are systematically higher than those of our sample; this could imply that our sample consist of more evolved PNs. On the contrary, IREs for our sample are systematically higher than those reported by AK91, implying that our sample is formed by PNs particularly young.
This apparent contradiction is further complicated by the fact the infrared properties of both samples are quite similar: sources belonging to differenr samples occupy the same region in the IRAS color-color diagram and dust temperatures of both samples have quite similar distributions, implying analogous dust characteristics.

A possible cause of lower $`T_B`$ and $`EM`$ of our sample can be a systematic effect due to the different spatial resolution used in the two surveys, as both $`T_B`$ and $`EM`$ are function of the source angular size ($`\theta ^2`$). However, only 4 out of 10 detected sources were not resolved.

The apparent contradiction can be explained if we assume that sources of both samples are in the ionization bounded phase of radio nebula evolution, but our sample is less evolved and is characterized, on the average, by a lower radio luminosity when compared to the AK91 sample. Consistently, ionized masses for our sample are in the range $`3\times 10^51.6\times 10^3M_{}`$, much lower than typical values for evolved nebulae (Pottasch, 1984).

The detection of free-free radio emission in 10 of the observed sources indicates that ionization is already started in their circumstellar shells. The detected sources are in the very early stage of PNs evolution and constitute a unique sample to be studied to shed light on this quite poorly understood phase of stellar evolution.

Successive multi-frequency and high-resolution radio observations will allow to fully characterize the radio properties of these new objects.

### Acknowledgments.

The Very Large Array is a facility of the National Radio Astronomy Observatory which is operated by Associated Universities, Inc. under cooperative agreement with the National Science Foundation

## References

Aaquist, O.B. and Kwok, S. 1991, ApJ, 378, 599 (AK91)

Conlon, E.S., Dufton, P.L., McCausland, R.J.H. and Keenan, F.P. 1993, ApJ, 408, 599

Parthasarathy, M. and Pottasch, S.R. 1989 A&A, 225, 521

Parthasarathy, M. 1993, Luminous High-Latitude Stars, ASP Conferences Series, Vol.45, 173

Parthasarathy, M., Vijapurkar, J. and Drilling, J.S. 2000, A&AS, 145 269

Pottasch, S.R., 1984, Planetary nebulae - a study of late stages of stellar evolution, Dordrecht, D. Reidel Publishiing Co.

Tajitsu, A. and Tamura, S. 1998, AJ, 115, 198

Terzian, Y. and Dickey, J. 1973, AJ, 78, 875


text:
# Does tiny-scale atomic structure exist in the interstellar medium ?

## 1 Introduction

For many years, both observations and theory have provided extensive support for the existence of structure in the interstellar medium (ISM) on scales from $``$ 1 kpc down to $``$ 1 pc (cf. Dickey & Lockman (1990)). However, it has long been expected that structure on smaller scales ($`<1`$ pc) would not be prominent in the ISM (Heiles, 2000). Indeed, using median values for thermal pressure and temperature of the cold neutral medium (CNM) of $`P_{\mathrm{th}}2250`$ cm-3 K (Jenkins & Tripp, 2001) and T $``$70 K (Heiles & Troland, 2003), the expected volume density for the CNM clouds is about 30 cm-3. The median measured column density of 5$`\times 10^{19}`$ cm-2 (Heiles & Troland, 2003), would then indicate that the typical expected scale length for the CNM features is $`1`$ pc, in conformance with the standard theory and observations.

Consequently, it was quite surprising when observers began to find structure on $`10^{12}`$ AU scales in many different directions in the ISM. The apparent detections of the AU-sized structure in the cold neutral atomic hydrogen (HI) medium, called tiny–scale atomic structure \[“TSAS”, Heiles (1997)\], were obtained with the following techniques: (1) spatial mapping of HI absorption line profiles across extended background sources (Dieter et al., 1976; Diamond et al., 1989; Davis et al., 1996; Faison & Goss, 2001); (2) temporal and spatial variations of optical interstellar absorption lines against binary stars and globular clusters (Meyer & Blades, 1996; Meyer & Lauroesch, 1999); and (3) time variability of HI absorption profiles against pulsars (Deshpande et al., 1992; Frail et al., 1994; Johnston et al., 2003). Since the AU-sized structure was detected very frequently, it was thought that it is likely to be a general property of the ISM.

TSAS observations yield measurements of optical depth variation ($`\mathrm{\Delta }\tau `$) and the particular temporal baseline over which the variation is measured. The baseline can also be converted into a transverse angular size ($`L_{}`$) if the background source speed or distance is known, respectively. A straightforward interpretation is that the optical depth variations are due to a blob, whose transverse dimension is equal to $`L_{}`$, moving into and out of the line-of-sight. Assuming a simple spherical geometry, the HI volume density of these blobs can be estimated and is typically of order of $`10^4`$ cm-3, far above the canonically expected value. The inferred thermal pressure is then of the order of $`10^6`$ cm-3 K (assuming T $``$40–70 K), which is also much higher than the hydrostatic equilibrium pressure of the ISM or the standard thermal pressure of the CNM. These problems have been known for a long time and have caused much controversy. It is expected that such over-dense and over-pressured features should dissipate on a time scale of about 100 yrs and should therefore not be common in the ISM, yet the TSAS observations have indicated that they are apparently ubiquitous. Furthermore, Jenkins & Tripp (2001) found evidence for the existence of over-pressured gas ($`P/k>10^5`$ cm-3 K) using CI fine-structure excitation, indicating that pressure enhancements could be widespread in the ISM and may even be related to the TSAS.

Several explanations were proposed to reconcile the observations of extensive TSAS with theory. Heiles (1997) proposed that TSAS features are curved filaments and/or sheets that happen to be aligned along our line-of-sight. Deshpande (2000) suggested that TSAS blobs correspond to the tail of a hierarchical structure organization that exists on larger scales, rather than discrete structures with longitudinal dimension $`L_{}`$ and extraordinarily high volume densities. Gwinn (2001) proposed that optical depth fluctuations seen in multi-epoch pulsar observations are a scintillation phenomenon combined with the velocity gradient across the absorbing HI. Different explanations predict a different level of optical depth variations at a particular scale size. For example, Deshpande (2000) expects that optical depth variations would increase with the size of structure, while Gwinn (2001) predicts maximum variations on the very small spatial scales probed by interstellar scintillation. All suggested explanations, however, call for more observational data.

Motivated by the recent theoretical efforts in understanding the nature and origin of the TSAS we have undertaken new multi-epoch observations of HI absorption against a set of bright pulsars. We chose to observe the same sources as Frail et al. (1994) in order to enhance the number of available time baselines for comparison. This paper summarizes our first results from this project for three pulsars, B0823+26, B1133+16 and B2016+28. In Section 2 we summarize our observing and data processing strategies. Results on individual objects are presented in Section 3 and discussed in Section 4.

## 2 Observations and Data Processing

We have used the Arecibo telescope111The Arecibo Observatory is part of the National Astronomy and Ionosphere Center, operated by Cornell University under a cooperative agreement with the National Science Foundation. to obtain new multi-epoch HI absorption measurements against six pulsars previously studied by Frail et. al. (1994). The observing and data processing procedures were very similar to Stanimirović et al. (2003), who measured pulsar OH absorption profiles at Arecibo with the same backend, the Caltech Baseband Recorder. We had four observing sessions: August 2000, December 2000, September 2001 and November 2001, searching for HI absorption profile variations over time intervals from less than a day to 1.25 years. Data from Frail et. al. (1994) extend the time baselines to a decade.

Our basic data product was a cube of pulsar intensity as a function of pulsar rotational phase and radio frequency. ‘Pulsar-on’ and ‘pulsar-off’ spectra were then accumulated by finding the pulsar pulse in software, and extracting spectra during the pulse and between pulses, respectively. The pulsar absorption spectrum is created by generating the ‘pulsar-on’ – ‘pulsar-off’ spectrum for each scan, doing frequency switching to flatten the baseline, and accumulating all such spectra with a weight proportional to $`T_{\mathrm{PSR}}^2`$, where $`T_{\mathrm{PSR}}`$ is the antenna temperature of the pulsar. Final absorption and emission spectra have velocity resolution of 0.5 km s-1.

## 3 Results on individual sources

For each pulsar, we display a single absorption spectrum on top, and then one or more absorption spectrum differences between two epochs. Overlaid atop each difference spectrum is a $`\pm `$2-$`\sigma `$ significance envelope. To calculate this expected noise level, we have taken into account the following contributions. (1) The sky background contribution was estimated from the all sky survey at 408 MHz by Haslam et al. (1982) and scaled to 1.4 GHz using a spectral index of $`2.6`$. (2) The contribution from HI emission was estimated from the ‘pulsar-off’ spectra and scaled appropriately to match previously published observations. The effect of this contribution can be very significant. For example the rms noise on–line is three times higher than the rms noise off-line for the case of B2016+28. (3) The pulsar continuum emission itself is also sufficiently strong to contribute an additional $`1030\%`$ to the noise temperature. Information on individual objects and an upper limit on $`\mathrm{\Delta }\tau `$ are given in Table 1 (note that these are 2-$`\sigma `$ limits at the HI line). Our measurements and the Frail et al. (1994) results are discussed for each source below.

### 3.1 PSR B0823+26

Fig. 1 (left) compares HI absorption profiles obtained toward B0823+26 in 2000.6 and 2000.9. We find no significant change in absorption over this time span down to a $`\mathrm{\Delta }\tau `$ level of about 0.04. The time interval of 0.3 yr translates to a transverse distance of about 10 AU. Frail et al. (1994) also found almost no variations over a period of 0.6 yr, but did report variations of about 0.07 over a period of 1.1 yr. See the Discussion section below for an examination of the discrepancies between the two groups’ results.

### 3.2 PSR B1133+16

Fig. 1 (right) displays HI absorption profiles obtained toward B1133+16 in 2000.6 and 2000.9. We detect no variations down to a $`\mathrm{\Delta }\tau `$ level of about 0.02. During this period, the pulsar traveled 20 AU. Frail et al. (1994) also detected no significant variations on their 0.6 yr baseline on this pulsar, but did see variations in $`\tau `$ of about 0.04 over a period of 1.1 yr. We investigate the discrepancy below.

### 3.3 PSR B2016+28–Short and long timescales

For this pulsar, we have compared HI absorption profiles on our usual short ($`0.3`$ yr) timescale, but also on decade–long scales. On the left side of Fig. 2, in the middle panel we show the usual difference of absorption spectra from 2000.6 and 2000.9, while the bottom left panel investigates an 8.7 year interval by differencing our 2000.6 spectrum and the Frail et al. (1994) 1991.9 data. The 0.3 yr baseline, corresponding to a $`3`$ AU scale, shows only a marginal 2.6-$`\sigma `$ change in absorption, while the 8.7 yr, $`70`$ AU baseline, does not exhibit any significant variations at all. This result is very different from Frail et al. (1994), who found very large optical depth variations of $`\mathrm{\Delta }\tau >1`$ over periods of 0.6 and 1.7 yr.

The Frail et al. (1994) epoch 1991.9 absorption spectrum that we use for the above comparison is particularly different from the previous three Frail et al. epochs (their first one actually being from Clifton et al. (1988)). While it agrees well with our epoch $`2000`$ results, their 1991.9 profile exhibits significantly shallower absorption in all four principal absorption features than do their earlier epoch results. Frail et al. (1994) remarked that the apparent change could be caused by an incorrect normalization but were confident that was not the case.

We further investigate the apparent Frail et al. variations by comparing our 2000.6 spectrum with the pre-1991.9 Frail et al. profiles. As expected, there is a large difference in the depth of all four absorption features, as shown in the top plot on the right side of Fig. 2 where the 1990.9 Frail et al. spectrum is plotted with a dot-dashed line and our 2000.6 result is plotted as a thin solid line. Since all four lines exhibit the same trend, it seems reasonable that the differences may result from a slight calibration problem in the Frail et al. data. To test this hypothesis, we performed a least–squares fit for a single scale factor (plus an offset) that would minimize the difference between these two spectra.

After applying this fitted scale factor ($`0.900\pm 0.005`$) and a constant ($`0.097\pm 0.004`$) to the 1990.9 spectrum, we arrive at the “scaled” 1990.9 spectrum, shown as a thick solid line in the top right plot of Fig. 2. Note that this process brings the previously discrepant 1990.9 profile into excellent agreement with the 2000.6 result. The bottom right plot of Fig. 2 emphasizes this result by exhibiting the difference between the 2000.6 and scaled 1990.9 spectra. A single 3-$`\sigma `$ deviation is all that remains. Note that without scaling, this change would be much larger. We conclude that most of the apparent variations among the Frail et al. B2016+28 spectra, and hence between our epoch 2000.6 data and the early Frail et al. epochs, could result from a systematic calibration issue in the Frail et al. data. The measurement of pulsar absorption spectra is a very delicate process fraught with subtle issues related to the strongly varying pulsar signal. We speculate that the three-level correlation spectrometer used by Frail et al. was more prone to calibration problems than our current four-level spectrometer.

## 4 Discussion

Our observations probe the existence of temporal variations in HI absorption profiles which would be indicative of the existence of the TSAS on spatial scales of several tens of AU in several different directions. The canonical TSAS has $`L_{}=30`$ AU, extracted by Heiles (1997) from previous observations, meaning that TSAS should be common at the scales that we have studied. Yet in all of the short and long baselines we have investigated, we find only two marginal detections of changing absorption.

Our sensitivity is better than any level of variations reported previously (Frail et al., 1994; Johnston et al., 2003), at velocity resolution of 0.5 km s-1, meaning that our non-detections are a significant result. In addition, the upper limits on optical depth fluctuations, set by B0823+26 and B1133+16 at scales of 10–20 AU, $`\mathrm{\Delta }\tau =0.030.04`$, are strikingly low. A $`\mathrm{\Delta }\tau `$ of 0.03 corresponds to $`5\times 10^{18}`$ cm-2 for the column density fluctuations of the CNM, for assumed spin temperature of 50 K and the full-width half-maxima of absorption features of 2 km s-1. On slightly larger scales of 70–80 AUs probed by B2016+28 over periods of almost a decade we find a single 3-$`\sigma `$ spike of $`\mathrm{\Delta }\tau =0.15`$, a $`2.5`$-$`\sigma `$ feature of $`\mathrm{\Delta }\tau 0.18`$, and non-detections in all other cases.

Deshpande (2000) predicts opacity variations as an extension of HI opacity irregularities observed on larger scales using a single power law spectrum. Using the power spectrum of opacity distribution in the direction of Cas A with the power-law index of 2.75 and extrapolating down to AU scales, they predict $`\mathrm{\Delta }\tau <0.1`$ at scales $`<20`$ AU (from their Fig. 2), while at 50–100 AU $`\mathrm{\Delta }\tau `$0.2–0.4 is expected. This model generally expects that $`\mathrm{\Delta }\tau `$ increases with spatial scales. Our results for all three pulsars seem consistent with this picture so far.

Why are the (primarily) non–detections presented here and by Johnston et al. (2003) so different from the results of Frail et al. (1994)? We showed above that in the case of B2016+28 where multiple absorption features are visible it is possible to simply scale absorption spectra from different epochs and get rid of most of the variations, indicating that slight calibration problems may be responsible. A similar scaling technique could significantly reduce variations seen in some of the other Frail et al.’s pulsars (e.g. B0823+16, B1133+6, B1737+13). However, it is not obvious that this scaling could be justified since only single principal absorption features are visible in the spectra. In addition, Johnston et al. (2003) have already noted that many of the differences noted by Frail et al. may not be as significant as they stated.

This large number of TSAS non detections, together with recent results by Johnston et al. (2003), is disturbing and unexpected! While previous observations were frequently detecting TSAS, in six comparisons presented here, and in many comparisons in Johnston et al. (2003), there are only two marginal detections for the case of B2016+28 and one detection for B1641-45 (Johnston et al. 2003) over a period of almost twenty years. This raises the question of the existence of TSAS as traced at least by multi-epoch pulsars observations. Contrary to the previous belief that the TSAS is ubiquitous in the ISM, our results indicate that TSAS is rather rare and could be related to some kind of isolated, sporadic events such as superbubble or supernova expansion, stellar mass ejecta, or a bow-shock propagation.

## 5 Conclusions

We have compared multi-epoch HI absorption observations toward B0823+26, B1133+16 and B2016+28 over short and long periods. Except for two marginal changes in the case of B2016+28, we do not find significant changes in absorption spectra. This is very different from previous observations by Frail et al. (1994) who saw significant optical depth variations for the same pulsars over time periods of 1.1 and 1.7 yrs. We have shown that in the case of B2016+28 most of the variations in absorption profiles of Frail et al. (1994) could be due to a systematic calibration problem. A large number of non detections of the TSAS presented here, together with recent results by Johnston et al. (2003), suggests that the TSAS is not ubiquitous in the ISM and could be a rare phenomenon.

We are grateful to Caltech’s Center for Advance Computation and Research for the use of their facilities for data storage and processing. We express our thanks to Stuart Anderson and Rick Jenet for help with observations and data processing. We thank Carl Heiles for stimulating discussions, and an anonymous referee for insightful suggestions. SS acknowledges support by NSF grants AST-0097417 and AST-9981308. KED, AH, JTG, and JMW were supported by NSF Grant AST-0098540.