text:
# 1 Variation of the relative perturbation corrections 𝐸_𝐺^(𝑛)/𝐸_𝐺⁽⁰⁾ to the ground state energy as a function of the coupling strength (𝑔₋) for 𝑡/𝜔₀=1.0 for (a) MLF method and (b) LF method. 𝐸_𝐺^(𝑛) is the nth order perturbation correction to the ground state energy and 𝐸_𝐺⁽⁰⁾= 𝐸⁽⁰⁾_{+,0} is the unperturbed ground state energy. 2,..,8 denote the value of 𝑛 in 𝐸_𝐺^(𝑛).

Two-site two-electron Holstein model:
a perturbation study

Jayita Chatterjee111 e-mail: moon@cmp.saha.ernet.in and A. N. Das

Theoretical Condensed Matter Physics Group
Saha Institute of Nuclear Physics
1/AF Bidhannagar, Kolkata 700064, India

PACS No.71.38. +i, 63.20.kr

Abstract

The two-site two-electron Holstein model is studied within a perturbation method based on a variational phonon basis obtained through the modified Lang-Firsov (MLF) transformation. The ground-state wave function and the energy are found out considering up to the seventh and eighth order of perturbation, respectively. The convergence of the perturbation corrections to the ground state energy as well as to the correlation functions are investigated. The kinetic energy and the correlation functions involving charge and lattice deformations are studied as a function of electron-phonon ($`e`$-ph) coupling and electron-electron interactions for different values of adiabaticity parameter. The simultaneous effect of the $`e`$-ph coupling and Coulomb repulsion on the kinetic energy shows interesting features.

1. Introduction

The Holstein model  is one of the fundamental models describing the interactions of conduction electrons with lattice vibrations and has been studied widely for a long time using different analytical and numerical methods. However, conventional analytical methods based on variational and perturbation approach are found to be much less satisfactory than the numerical methods; particularly in a wide region of intermediate values of electron-phonon ($`e`$-ph) coupling and hopping. The importance of a reliable analytical approach is crucial due to the failure of standard analytical approximations on the Holstein model; the Migdal approximation  is valid for a weakly coupled adiabatic system, while the Lang-Firsov canonical transformation  could not be suitable far from the nonadiabatic limit and strong coupling limit. The failure of the standard analytical techniques in the intermediate range of coupling has been observed by many workers  while comparing with exact results. Recently we  developed a perturbative expansion using different phonon bases obtained through the LF, MLF and MLF with squeezing (MLFS) transformations for a two-site single electron Holstein model. We found that the perturbation corrections within the MLF and MLFS methods are much smaller than the corresponding LF values for weak and intermediate couplings. The MLF method shows a very satisfactory convergence in perturbation expansion and is able to produce almost exact results for the entire region of the $`e`$-ph coupling from weak to intermediate values of hopping.

Our aim is to extend the applicability of our approach, developed for a single electron system, to a many-electron system. To examine that whether this perturbation method based on MLF phonon basis works well or not with increasing number of electrons and in presence of electron-electron ($`e`$-$`e`$) interaction, we consider a two-site two-electron Holstein model which is the minimal model to investigate the effect of $`e`$-$`e`$ interaction on the polaronic properties and the combined effects of the electronic correlation and the $`e`$-ph coupling on the kinetic energy of the system. It may be mentioned that the polaron crossover has been studied in two-site and four-site systems in presence of $`e`$-$`e`$ interaction by different groups  with the Holstein model previously. But most of the studies were done at the zero phonon approximation level with MLF or MLFS transformed Hamiltonian. However, the combined effect of $`e`$-$`e`$ interaction and $`e`$-ph interaction on the kinetic energy has not been properly studied and this is the main motivation of the present study.

In our study, we have developed a perturbation series using the MLF phonon basis and examined the convergence by comparing contributions of different orders of perturbation to the energy and the correlation functions for the ground state. The variation of different correlation functions involving charge and lattice deformation and the kinetic energy of the system with $`e`$-ph coupling are presented and discussed for different values of adiabaticity parameters and Coulomb repulsion.

2. Model Hamiltonian and formalism

The two-site two-electron Holstein Hamiltonian is

$`H`$ $`=`$ $`{\displaystyle \underset{i,\sigma }{}}ϵn_{i\sigma }{\displaystyle \underset{\sigma }{}}t(c_{1\sigma }^{}c_{2\sigma }+c_{2\sigma }^{}c_{1\sigma })+U{\displaystyle \underset{i}{}}n_in_i+Vn_1n_2`$ (1)
$`+`$ $`g_1\omega _0{\displaystyle \underset{i,\sigma }{}}n_{i\sigma }(b_i+b_i^{})+g_2\omega _0{\displaystyle \underset{i,\sigma }{}}n_{i\sigma }(b_{i+\delta }+b_{i+\delta }^{})+\omega _0{\displaystyle \underset{i}{}}b_i^{}b_i`$
where $`i`$ =1 or 2, denotes the sites. $`c_{i\sigma }`$ ($`c_{i\sigma }^{}`$) is the annihilation (creation) operator for the electron with spin $`\sigma `$ at site $`i`$ and $`n_{i\sigma }`$ (=$`c_{i\sigma }^{}c_{i\sigma }`$) is the corresponding number operator, $`g_1`$ and $`g_2`$ denote, respectively, the on-site and inter-site $`e`$-ph coupling strengths. $`t`$ is the usual hopping integral. $`b_i`$ ($`b_i^{}`$) is the annihilation (creation) operator for the phonons corresponding to interatomic vibrations at site $`i`$, $`\omega _0`$ is the phonon frequency, $`U`$ and $`V`$ are, respectively, the on-site and inter-site Coulomb repulsion between the electrons.

Introducing new phonon operators $`𝐚=(b_1+b_2)/\sqrt{2}`$ and $`d=(b_1b_2)/\sqrt{2}`$ the Hamiltonian is separated into two parts ($`H=H_d+H_𝐚)`$ :

$`H_d`$ $`=`$ $`{\displaystyle \underset{i,\sigma }{}}ϵn_{i,\sigma }{\displaystyle \underset{\sigma }{}}t(c_{1\sigma }^{}c_{2\sigma }+c_{2\sigma }^{}c_{1\sigma })+U{\displaystyle \underset{i}{}}n_in_i+Vn_1n_2`$ (2)
$`+`$ $`\omega _0g_{}(n_1n_2)(d+d^{})+\omega _0d^{}d`$
and

$$H_𝐚=\omega _0\stackrel{~}{𝐚}^{}\stackrel{~}{𝐚}\omega _0n^2g_+^2$$
(3)
where $`g_+=(g_1+g_2)/\sqrt{2}`$, $`g_{}=(g_1g_2)/\sqrt{2}`$, $`\stackrel{~}{𝐚}=𝐚+ng_+`$, $`\stackrel{~}{𝐚}^{}=𝐚^{}+ng_+`$, $`n=n_1+n_2`$ and $`n_i=n_i+n_i`$.

$`H_𝐚`$ describes a shifted oscillator which couples only with the total number of electrons $`n`$, which is a constant of motion. $`H_d`$ represents an effective $`e`$-ph system where phonons directly couple with the electronic degrees of freedom. We now use the MLF transformation where the lattice deformations produced by the electron are treated as variational parameters . For the present system the relevant transformation is,

$$\stackrel{~}{H_d}=e^RH_de^R$$
(4)
where $`R=\lambda (n_1n_2)(d^{}d)`$, $`\lambda `$ is a variational parameter related to the displacement of the $`d`$ oscillator.

The transformed Hamiltonian is then obtained as

$`\stackrel{~}{H_d}`$ $`=`$ $`\omega _0d^{}d+{\displaystyle \underset{i,\sigma }{}}ϵ_pn_{i\sigma }{\displaystyle \underset{\sigma }{}}t[c_{1\sigma }^{}c_{2\sigma }\mathrm{exp}(2\lambda (\mathrm{d}^{}\mathrm{d}))+\mathrm{c}_{2\sigma }^{}\mathrm{c}_{1\sigma }\mathrm{exp}(2\lambda (\mathrm{d}^{}\mathrm{d}))]`$ (5)
$`+`$ $`\omega _0(g_{}\lambda )(n_1n_2)(d+d^{})+U_e{\displaystyle \underset{i}{}}n_in_i+V_en_1n_2`$
where

$`ϵ_p`$ $`=`$ $`ϵ\omega _0(2g_{}\lambda )\lambda `$ (6)
$`U_e`$ $`=`$ $`U2\omega _0(2g_{}\lambda )\lambda `$ (7)
$`V_e`$ $`=`$ $`V+2\omega _0(2g_{}\lambda )\lambda `$ (8)
For the two-site two-electron system there are six electronic states of which three belong to the triplet states as follows:

$`|\psi _{T1}`$ $`=`$ $`c_1^{}c_2^{}|0_e`$ (9)
$`|\psi _{T2}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(c_1^{}c_2^{}+c_1^{}c_2^{})|0_e`$ (10)
$`|\psi _{T3}`$ $`=`$ $`c_1^{}c_2^{}|0_e`$ (11)
and other three states are singlet as given below:

$`|\psi _1`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(c_1^{}c_1^{}c_2^{}c_2^{})|0_e`$ (12)
$`|\psi _2`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(c_1^{}c_1^{}+c_2^{}c_2^{})|0_e`$ (13)
$`|\psi _3`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(c_1^{}c_2^{}c_1^{}c_2^{})|0_e`$ (14)
The triplet states do not mix with each other or couple with the $`d`$-oscillator. In the following we will only consider the singlet states and the $`d`$-oscillator which are mutually coupled through the off-diagonal terms of the Hamiltonian (5). For a perturbation method it is desirable to use a basis where the major part of the Hamiltonian becomes diagonal. For this purpose we choose the basis set $`|\psi _+,N`$, $`|\psi _1,N`$ and $`|\psi _{},N`$ defined as

$`|\psi _+,N`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{a_N^2+b_N^2}}}(a_N|\psi _2+b_N|\psi _3)|N`$ (15)
$`|\psi _1,N`$ $`=`$ $`|\psi _1|N`$ (16)
$`|\psi _{},N`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{a_N^2+b_N^2}}}(a_N|\psi _3b_N|\psi _2)|N`$ (17)
where $`|\psi _+`$, $`|\psi _1`$ and $`|\psi _{}`$ represent the bonding, nonbonding and antibonding electronic states, respectively and $`|N`$ denotes the $`N`$th excited oscillator state in the MLF phonon basis where $`N`$ takes integer values from 0 to $`\mathrm{}`$.
Here $`a_N=(U_eV_e)+\sqrt{(U_eV_e)^2+16t_e^2(N,N)}`$ and $`b_N=4t_e(N,N)`$.
where $`t_e(N,N)=t_e_{p=1}^N(1)^p\frac{(2\lambda )^{2p}}{p!}N_{C_p}`$, $`N_{C_p}=\frac{N!}{p!(Np)!}`$,
and $`t_e=t\mathrm{exp}(2\lambda ^2)`$.

The hopping term in Eq. (5) has both diagonal and offdiagonal matrix elements in the chosen basis (Eqs. 15-17). The states $`|\psi _e,N`$ (for $`e=+,1,`$) are the eigenstates of the unperturbed Hamiltonian

$`\stackrel{~}{H_0}`$ $`=`$ $`\omega _0d^{}d+{\displaystyle \underset{i}{}}ϵ_pn_i{\displaystyle \underset{\sigma }{}}t_eF(d^{}d)[c_{1\sigma }^{}c_{2\sigma }+c_{2\sigma }^{}c_{1\sigma }]`$ (18)
$`+`$ $`U_e{\displaystyle \underset{i}{}}n_in_i+V_en_1n_2`$
It may be noted that $`t_eF(d^{}d)`$ in the above Eq. (18) is the diagonal part of the operator $`t\mathrm{exp}[\pm 2\lambda (d_i^{}d_i)]`$ in Eq. (5) in the phonon occupation number representation.

$`\mathrm{Since}\mathrm{exp}[2\lambda (d^{}d)]`$ $`=`$ $`\mathrm{exp}[2\lambda d^{}]\mathrm{exp}[2\lambda d]\mathrm{exp}[{\displaystyle \frac{4\lambda ^2}{2}}[d^{},d]_{}]`$
$`=`$ $`e^{2\lambda ^2}{\displaystyle \underset{p}{}}{\displaystyle \underset{q}{}}(1)^q{\displaystyle \frac{(2\lambda )^{p+q}}{p!q!}}(d^{})^p(d)^q,`$
$`\mathrm{so}F(d^{}d)`$ $`=`$ $`{\displaystyle \underset{p=1}{\overset{N}{}}}(1)^p(2\lambda )^{2p}(d^{}d)^p/(p!)^2.`$
For the choice of the unperturbed Hamiltonian in Eq. (18) the perturbation Hamiltonian ($`\stackrel{~}{H}\stackrel{~}{H_0}`$) has no diagonal matrix element in our chosen basis. This is not achieved when one replaces $`F(d^{}d)`$ by 1 and $`a_N`$, $`b_N`$ become $`N`$-independent. The convergence in perturbation expansion for the latter case is found to be weaker than that when $`F(d^{}d)`$ is included in $`H_0`$.

The unperturbed energies corresponding to the eigenstates $`|\psi _e,N`$ ($`e=+,1,`$) are

$`E_{+,N}^{(0)}`$ $`=`$ $`N\omega _0+2ϵ_p+{\displaystyle \frac{(U_e+V_e)}{2}}{\displaystyle \frac{1}{2}}\sqrt{(U_eV_e)^2+16t_e^2(N,N)},`$
$`E_{1,N}^0`$ $`=`$ $`2ϵ_p+U_e`$
$`\mathrm{and}\mathrm{E}_{,\mathrm{N}}^{(0)}`$ $`=`$ $`N\omega _0+2ϵ_p+{\displaystyle \frac{(U_e+V_e)}{2}}+{\displaystyle \frac{1}{2}}\sqrt{(U_eV_e)^2+16t_e^2(N,N)}`$ (19)
The part of the Hamiltonian, which is not included in $`\stackrel{~}{H_0}`$, is taken as the perturbation $`\stackrel{~}{H_1}`$. For the chosen basis, the state $`|\psi _+,0`$ has the lowest unperturbed energy $`E_{+,0}^{(0)}`$. The general off-diagonal matrix elements of $`\stackrel{~}{H_1}`$ between the states $`|\psi _\pm ,N`$ and $`|\psi _1,M`$ are calculated as (for $`(NM)>0`$)

$`N,\psi _\pm |\stackrel{~}{H_1}|\psi _\pm ,M`$ $`=`$ $`2t_e(N,M)(a_N^{^{}}b_M^{^{}}+a_M^{^{}}b_N^{^{}})\mathrm{for}\mathrm{even}(\mathrm{N}\mathrm{M}).`$ (20)
$`=`$ $`0\mathrm{for}\mathrm{odd}(\mathrm{N}\mathrm{M}).`$
$`N,\psi _\pm |\stackrel{~}{H_1}|\psi _{},M`$ $`=`$ $`2t_e(N,M)(a_N^{^{}}a_M^{^{}}b_M^{^{}}b_N^{^{}})\mathrm{for}\mathrm{even}(\mathrm{N}\mathrm{M}).`$ (21)
$`=`$ $`0\mathrm{for}\mathrm{odd}(\mathrm{N}\mathrm{M}).`$
$`N,\psi _1|\stackrel{~}{H_1}|\psi _+,M`$ $`=`$ $`2a_M^{^{}}P(N,M)2t_e(N,M)b_M^{^{}}\mathrm{for}\mathrm{odd}(\mathrm{N}\mathrm{M}).`$ (22)
$`=`$ $`0\mathrm{for}\mathrm{even}(\mathrm{N}\mathrm{M}).`$
$`N,\psi _1|\stackrel{~}{H_1}|\psi _{},M`$ $`=`$ $`2b_M^{^{}}P(N,M)2t_e(N,M)a_M^{^{}}\mathrm{for}\mathrm{odd}(\mathrm{N}\mathrm{M}).`$ (23)
$`=`$ $`0\mathrm{for}\mathrm{even}(\mathrm{N}\mathrm{M}).`$
$`\mathrm{where}t(N,M)`$ $`=`$ $`t_e(2\lambda )^{NM}\sqrt{{\displaystyle \frac{N!}{M!}}}[{\displaystyle \frac{1}{(NM)!}}+{\displaystyle \underset{R=1}{\overset{M}{}}}[(1)^R`$
$`{\displaystyle \frac{(2\lambda )^{2R}}{(NM+R)!R!}}M(M1)\mathrm{}(MR+1)]]`$
$`P(N,M)`$ $`=`$ $`\omega _0(g_{}\lambda )[\sqrt{M}\delta _{N,M1}+\sqrt{M+1}\delta _{N,M+1}]`$
$`a_N^{^{}}`$ $`=`$ $`{\displaystyle \frac{a_N}{\sqrt{(a_N^2+b_N^2)}}}\mathrm{and}b_N^{^{}}={\displaystyle \frac{b_N}{\sqrt{(a_N^2+b_N^2)}}}`$
Now one has to make a proper choice of $`\lambda `$ so that the perturbative expansion is satisfactorily convergent. In our previous work we obtained good convergence of the perturbation series for the value of $`\lambda `$ which minimizes the unperturbed ground state energy. Following the same procedure, i.e. minimizing $`E_{+,0}^{(0)}`$ with respect to $`\lambda `$ we obtain

$`\lambda =g_{}/[1+{\displaystyle \frac{4t_e}{\omega _0[\sqrt{r^2+4}r]}}]`$ (24)
where $`r=(U_eV_e)/2t_e`$.

The corrected ground state wave function may be written as,

$`|\psi _G`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N_G}}}\left[|\psi _+,0+{\displaystyle \underset{N=1,3,..}{}}C_N^1|\psi _1,N+{\displaystyle \underset{N=2,4,..}{}}C_N^\pm |\psi _\pm ,N\right]`$ (25)
$`\mathrm{where}N_G`$ $`=`$ $`1+{\displaystyle \underset{N=1,3,..}{}}(C_N^1)^2+{\displaystyle \underset{N=2,4,\mathrm{}}{}}((C_N^+)^2+(C_N^{})^2)\mathrm{is}\mathrm{the}\mathrm{normalization}`$
$`\mathrm{factor}.`$
The coefficients $`C_N^1`$, $`C_N^+`$ and $`C_N^{}`$ for different values of $`N`$ are determined using the matrix elements of $`\stackrel{~}{H_1}`$ connecting different unperturbed states and their energies following Ref. . These Coefficients are determined up to the seventh order of perturbation while the ground-state energy is found out up to the eighth order of perturbation.

We have also studied the static correlation functions involving charge and lattice deformations $`n_1u_1_0`$ and $`n_1u_2_0`$, where $`u_1`$ and $`u_2`$ represents the lattice deformations at sites 1 and 2 respectively. These correlation functions may be written in terms of relevant operators as

$$n_1u_{1,2}=n_1[2g_++\frac{(𝐚+𝐚^{})}{2}]\pm n_1[\lambda (n_1n_2)+\frac{(d+d^{})}{2}]$$
(26)
The final form of the correlation functions are obtained as

$`n_1u_{1,2}_0`$ $`=`$ $`2g_+{\displaystyle \frac{2\lambda }{N_G}}[{\displaystyle \frac{a_0^2}{a_0^2+b_0^2}}+{\displaystyle \underset{N=1,3,5,..}{}}(C_N^1)^2`$ (27)
$`+`$ $`{\displaystyle \underset{N=2,4,..}{}}{\displaystyle \frac{1}{a_N^2+b_N^2}}(a_NC_N^+b_NC_N^{})^2]\pm {\displaystyle \frac{1}{N_G}}[{\displaystyle \frac{a_0}{\sqrt{a_0^2+b_0^2}}}C_1^1`$
$`+`$ $`{\displaystyle \underset{N=2,4,..}{}}{\displaystyle \frac{\sqrt{N}C_{N1}^1}{\sqrt{a_N^2+b_N^2}}}(a_NC_N^+b_NC_N^{})`$
$`+`$ $`{\displaystyle \underset{N=3,5,..}{}}{\displaystyle \frac{\sqrt{N}C_N^1}{\sqrt{a_N^2+b_N^2}}}(a_{N1}C_{N1}^+b_{N1}C_{N1}^{})]`$
The kinetic energy of the system in the ground state $`E_{ke}=t<\psi _G|`$ $`_\sigma [c_{1\sigma }^{}c_{2\sigma }`$
$`\mathrm{exp}(2\lambda (d^{}d))`$ $`+c_{2\sigma }^{}c_{1\sigma }\mathrm{exp}(2\lambda (d^{}d))]`$ $`|\psi _G>`$, has also been calculated using the ground state wave function $`|\psi _G`$, evaluated up to the seventh order of perturbation.

3. Results and discussions

From the relevant analytical expressions, the quantities of our interest have been calculated within the LF and the MLF perturbation methods considering 35 phonon states (which is more than sufficient for $`g_{}2.5`$ in the transformed phonon basis). For all numerical calculations we take $`ϵ=0`$, $`g_2=0`$ (i.e., $`g_+`$ = $`g_{}`$) in this work.

In figure 1 we have shown the relative perturbation corrections to the ground-state energy, i.e., the ratios of the perturbation corrections of different orders to the unperturbed ground-state energy as a function of $`g_{}`$ for $`t/\omega _0=1`$. It appears from figure 1 that the convergence is more or less satisfactory in both the MLF and LF approaches. Apart from the second order correction, higher order energy corrections are very small within the MLF approach except in a small region of $`g_{}0.6`$ to 0.8 where the convergence is relatively weaker. The energy corrections of fifth to eighth order are, however, negligibly small even in this region. The energy correction of any order within the LF method is larger than the corresponding correction within the MLF approach. For smaller values of $`t/\omega _0`$ the corrections are much less and the convergence is better as expected from our previous study .

To examine the applicability of our method for higher values of $`t`$ and also in presence of $`e`$-$`e`$ interaction, we plot in figure 2 the relative perturbation corrections of different orders to the ground-state energy as a function of $`g_{}`$ for $`t/\omega _0=2`$ for both zero and nonzero values of $`e`$-$`e`$ interactions. Figure 2(a) shows that the convergence is satisfactory and except the second order correction, higher order energy corrections are small. In presence of electronic interactions (nonzero $`U`$ and $`V`$), magnitudes of the relative energy corrections increase but the convergence is still good (figure 2(b)). We calculate the ground-state energy within the MLF method up to the eighth order of perturbation. As pointed out previously, the fifth to eighth order energy corrections are so small that the energy calculated up to the eighth order of perturbation may be treated as an almost exact result. We find that our results exactly match with that obtained by the exact diagonalization study  for the same set of parameters.

The on-site and inter-site charge-lattice deformation correlations, $`n_1u_1_0`$ and $`n_1u_2_0`$, are determined for the ground state considering up to the seventh order of perturbation within the MLF method. In figure 3 we have plotted these quantities against $`g_{}`$. When we examine the convergence by comparing the values of the same correlation function, calculated up to different orders of perturbation we find that the values obtained considering up to the fourth to the seventh order of perturbation, remain almost same. This ensures a good convergence of our perturbation method for the correlation function and thus the results obtained are expected to match satisfactorily with the exact result. The unperturbed LF and the MLF results are also plotted in figure 3 for a comparison, which shows that the unperturbed MLF results for the correlation functions are much closer to the exact result than the corresponding LF result. The relative difference between the exact MLF result and the MLF (ZPA) result is much less for $`n_1u_1_0`$ than for $`n_1u_2_0`$. This is owing to the fact that $`n_1u_2_0`$ is a very sensitive function of the perturbation corrections to the wave function compared to $`n_1u_1_0`$. The MLF perturbation result predicts a long tail for $`n_1u_2_0`$ in the strongly coupled region where both the LF (ZPA) and the MLF (ZPA) give zero value for $`n_1u_2_0`$. This result is a signature of the retardation effect even in the strong coupling region where even the MLF (ZPA) do not predict the retardation effect.

To show the crossover from the delocalized (large) to localized (small) polaron we plot the quantity $`\chi _{12}`$ (= $`(n_1n_2)(u_1u_2)_0/g_+`$) against $`g_{}`$ in figure 4(a) where the MLF perturbation results for different values of $`t/\omega _0`$ are shown. The ZPA results within MLF and LF methods are also given for $`t/\omega _0=1.0`$ to compare them with the exact perturbation results. The MLF (ZPA) result is very close to the MLF perturbation result. The MLF perturbation result shows a smooth crossover from large to small polaron as $`e`$-ph coupling strength increases. If we identify the crossover point as the point of inflection of the curve obtained within the MLF (perturbation) approach we find that for $`t/\omega _0=1`$ the crossover occurs at a critical $`e`$-ph coupling $`g_{}`$ ($`g_{}0.8`$) when $`U=V=0`$. At this point the coefficient of one phonon state ($`C_{N=1}^1`$ in Eq. (25)) changes its sign. Similar behaviour has been also observed in case of our previous single polaron studies . The critical value of the $`e`$-ph coupling where polaron crossover occurs increases with increasing value of $`t/\omega _0`$, as expected.

In figure 4(b) we have shown the effect of $`e`$-$`e`$ interaction on the polaron crossover by plotting $`\chi _{12}`$ against $`g_{}`$ for several values of $`U`$ and $`V`$. For this two-site two-electron system the polaronic properties we have studied, depend on $`(UV)`$ rather than $`U`$ and $`V`$ independently. From figure 4(b) it is seen that the critical coupling, where the polaron crossover takes place, increases and the crossover becomes more abrupt with increasing $`(UV)`$. Increasing on-site Coulomb repulsion makes the crossover region more abrupt and extends the large polaron region. Similar results were also reported earlier . The inter-site Coulomb repulsion, on the other hand, makes the crossover smoother and the crossover takes place at a lower value of $`e`$-ph coupling.

The variation of the kinetic energy of the system with $`e`$-ph coupling strength $`g_{}`$ is shown in figure 5 for $`U=V=0`$ and for different values of $`t/\omega _0`$ (= 0.1, 0.5, 1.0 and 1.5). We find that our perturbation results for the kinetic energy match exactly with the results obtained by exact diagonalization study  for the corresponding set of parameters. It may be mentioned that for $`t/\omega _0=`$ 1.0, 1.5 the convergence is not good in a very small region of $`g_{}`$ around the points indicated by the arrows in the figure. Except these narrow regions the convergence of our perturbation method is excellent in the rest of the parameter space we studied.

It is well known that the kinetic energy of a tight-binding electronic system is suppressed by the introduction of either $`e`$-ph coupling or by Coulomb repulsion. An important question one may ask what is the simultaneous role of the Coulomb repulsion and the $`e`$-ph interaction in lowering the kinetic energy of the system. To address this question we have plotted the kinetic energy against $`e`$-ph coupling strength for different values of $`U`$ for $`t/\omega _0=1.0`$ in figure 6. It is seen from figure 6 that for large $`U`$ values there is a flat region where the kinetic energy remains unaffected by increase of $`e`$-ph coupling. Higher the value of $`U`$ larger is the flat region. Beyond this flat region the kinetic energy is suppressed almost exponentially with increasing $`g_{}`$. The occurance of this flat region and then a rapid fall in the kinetic energy for higher values of $`U`$ may be understood from the results of figure 4(b) where it is found that for larger values of $`U`$ the polaron crossover is abrupt but there is an extended large polaron region (in $`g_{}`$ space). The value of $`\chi _{12}`$, which is a measure of the difference in the lattice distortions produced at the charge residing site and the neighbouring site, is small and almost constant in this large polaron region. As a result a flat behaviour in the kinetic energy is obtained. Beyond a critical value of $`g_{}`$ an abrupt increase in $`\chi _{12}`$ would lead to a sharp fall in the kinetic energy, which is evident in figure 6.

In figure 7 the variation of the kinetic energy with $`U`$ for different values of $`g_{}`$ is shown. It is found that for $`g_{}=1.0`$ and $`1.5`$ there is a region where the kinetic energy increases with $`U`$ in contrast to the usual behaviour that the kinetic energy is suppressed by the Coulomb correlation. The latter behaviour is of course seen in major parameter space. To understand this unusual feature of increasing kinetic energy with $`U`$ for finite $`g_{}`$, we plot the variation of $`\chi _{12}`$ as well as the kinetic energy with $`U`$ in the inset of figure 7 for $`g_{}=1.0`$. It is readily seen that as $`U`$ increases the correlation function $`\chi _{12}`$ decreases which results in an increase in the kinetic energy. At the crossover region, $`\chi _{12}`$ shows a sharp fall which is reflected by the abrupt rise in the kinetic energy. For large values of $`U`$, $`\chi _{12}`$ is small and the polaron is a large one. In this case the lattice distortions uniformly spread throughout the whole lattice (hence the charge particle behaves more like a normal electron) and the kinetic energy is suppressed with increasing $`U`$ as commonly expected.

In figure 8 the variation of the kinetic energy with $`U`$ for different values of $`t/\omega _0`$ is shown for a fixed value of $`e`$-ph coupling ($`g_{}=1.0`$). The kinetic energy increases with $`U`$, shows a peak and then decreases. Similar behaviour is also observed in figure 7. The peak is broader for higher values of $`t`$. The initial enhancement of the kinetic energy with $`U`$ is due to rapid spread up of the polaron size or decrease in the value of $`\chi _{12}`$ whereas in the large polaron region $`U`$ plays its conventional role of suppressing the kinetic energy. With increasing $`t`$ the kinetic energy increases as well as the large polaron region is extended. We have also plotted $`\chi _{12}`$ in the same figure for $`t/\omega _0=1.5`$ to show the rise in the kinetic energy is related to rapid fall in $`\chi _{12}`$. For $`t/\omega _0=2.0`$, the polaron is a large one for $`g_{}=1`$ and so the small polaron region is not seen for this case.

In summary we apply the MLF perturbation method on a two-site two-electron Holstein model taking the diagonal part of the Hamiltonian (in the momentum space of the MLF basis) as the unperturbed Hamiltonian. A good convergence is achieved by our approach. The ground state energy and the kinetic energy obtained within this method match exactly with those obtained by exact diagonalization study. The behaviour of different correlation functions involving charge and the lattice distortions are reported as a function of $`e`$-ph coupling. Effect of the Coulomb repulsion and the adiabaticity parameter on the large to small polaron crossover is studied. On-site Coulomb repulsion makes the polaron crossover more abrupt. Our study on the kinetic energy in presence of both $`e`$-ph interaction and Coulomb repulsion shows that for large values of $`U`$, the kinetic energy is not suppressed by $`e`$-ph coupling in a range which is wider larger the value of $`U`$. In this range the polaron is a large one with almost an uniform value of $`\chi _{12}`$. This makes the kinetic energy insensitive to the e-ph coupling. When $`e`$-ph coupling crosses a critical value a rapid crossover to a small polaron occurs and in this region the $`e`$-ph coupling is effective in suppressing the kinetic energy. On the other hand for finite $`e`$-ph coupling the kinetic energy increases initially with $`U`$ and then decays. This is ascribed to the fact that the former region corresponds to a small polaron or a crossover region where the main role of $`U`$ is to spread the size of the polaron and make it more delocalized. This results in an increase in the kinetic energy whereas for large $`U`$ values the polaron is a large one and the on-site correlation suppresses the kinetic energy in its usual way.


format: text

text:
# Inner Structure of Protostellar Collapse Candidate B335 Derived from Millimeter-Wave Interferometry1footnote 11footnote 1Based on observations carried out with the IRAM Plateau de Bure Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain).

## 1 Introduction

The dense core in the B335 dark globule is generally recognized as the best protostellar collapse candidate. This dense core is nearby (250 pc, Tomita, Saito & Ohtani 1979), isolated, and nearly spherical. It contains a deeply embedded low luminosity young stellar object ($`3`$ L) discovered at far-infrared wavelengths by Keene et al. (1983) and detected by IRAS only at $`\lambda 60`$ $`\mu `$m. Detailed radiative transfer models based on the theory of inside-out collapse (Shu 1977) provide very good fits to spectral line profiles of the dense gas tracers CS and H2CO observed at $`10^{\prime \prime }`$$`30^{\prime \prime }`$ resolution (Zhou et al. 1993; Choi et al. 1995).

Recent studies have raised serious doubts about the inside-out collapse interpretation of the molecular line profiles. Observations of the CS(5-4) line at higher angular resolution, about 500 AU, show that the high velocity emission arises from the inner part of a bipolar outflow, not from the gravitational acceleration of dense gas close to the protostar (Wilner et al. 2000). In addition, studies of starless cores show that some molecular species, in particular CS, are frozen out onto grains at densities characteristic of the inner regions of B335 (Caselli et al. 2002, Tafalla et al. 2002). If protostellar heating does not promptly desorb CS molecules from icy grain mantles, then this species likely makes a poor probe of the gas kinematics in the presumed infall zone.

Since the density field is strongly coupled to the velocity field, observations of dust emission and absorption provide an alternative means to infer the dynamical state of the dense core. Harvey et al. (2001) describe deep near-infrared extinction measurements toward B335 that show a density profile consistent with the inside-out collapse theory, in particular an $`r^2`$ density profile for the outer envelope and an inner turnover towards free-fall. The best fit model has an infall radius of 6500 AU ($`26^{\prime \prime }`$), similar to that derived from molecular line studies (Zhou et al. 1993; Choi et al. 1995). The extinction data also show the dramatic effect of the bipolar outflow driven by the protostar. However, extinction measurements could not be made interior to $`3500`$ AU ($`14^{\prime \prime }`$), where no background stars are visible, even with NICMOS on HST. As a result, the extinction data could not distinguish the inside-out collapse model from, e.g. a highly unstable Bonnor-Ebert sphere, whose density profiles depart significantly from each other only at radii $`2500`$ AU.

Observations of long wavelength dust emission provide another method for probing density structure. This technique is nearly as direct as observations of dust extinction. The intensity of the emission provides an integral along the line-of-sight of the product of the density, temperature, and opacity of the dust. By modeling the dust temperature and specific mass opacity, the observed intensities may be used to constrain the density distribution. The analysis of dust emission complements extinction work because the dust emission method becomes most effective in the high column density inner regions of dense cores where dust extinction becomes large and difficult to penetrate.

Dust emission has been used very successfully to measure density structure of dense cores using bolometer cameras at millimeter and submillimeter wavelengths on large single dish telescopes, including the IRAM 30 meter (Ward-Thompson, Motte & Andre 1999, Motte & Andre 2001) and the JCMT (Ward-Thompson et al. 1994; Visser, Richer & Chandler 2001, Shirley et al. 2000). Unfortunately, the central regions of the nearest protostellar cores are generally comparable in size to the beamwidths of these telescopes and are poorly resolved. The density structure at smaller scales can only be probed with interferometers. With interferometer observations, the density structure may be strongly constrained through analysis of dust emission directly in the visibility domain. So far, detailed modeling of millimeter continuum visibility data in this context has been rarely applied (Keene & Masson 1990, Hogerheijde et al. 1999, Mundy, Looney & Welch 2000)

B335 has been imaged with millimeter interferometers in several previous studies (Chandler & Sargent 1993, Hirano et al. 1992, Saito et al. 1999). In this paper, we present observations of dust continuum emission from B335 at both 3 mm and 1.2 mm obtained with the IRAM Plateau de Bure Interferometer (PdBI) that obtain considerably better sensitivity. These observations sample well the density profile of B335 from size scales of $`5000`$ AU to $`500`$ AU, affording an order of magnitude better resolution than has been possible with single dish telescopes. The visibilities can discriminate the behaviour of the various physical models that all match the near infrared extinction in the outer regions of B335.

## 2 Observations

Continuum emission from B335 was observed with the IRAM PdBI simultaneously with the CS(5-4) line emission reported by Wilner et al. (2000). In brief, a single pointing was observed using 5 antennas in the D configuration giving baseline lengths from the shadowing limit of 15 meters to nearly 80 meters. The tuning was single (lower) sideband at 100 GHz (3.0 mm) and double sideband at 245 GHz (1.2 mm). The half power field of view is $`50^{\prime \prime }`$ at 3.0 mm and $`20^{\prime \prime }`$ at 1.2 mm. The absolute flux scale was determined through observations of the standard source MWC349, assumed to be 1.09 Jy at 3.0 mm at 1.82 Jy at 1.2 mm, with estimated uncertainties of 20% and 10%, respectively. Continuum visibility records at the two wavelengths were formed for each 60 second integration of the digital correlator ($`3\times 160`$ MHz at 1.2 mm, $`1\times 160`$ MHz at 3.0 mm), excluding a few channels contaminated by molecular line emission, which resulted in 2372 records at 3.0 mm and 4824 records at 1.2 mm (after flagging occasional bad records). In addition to amplitude and phase, each record also contains a variance measure, determined from the system temperatures and antenna gains.

## 3 Constructing Model Visibilities

We analyze the dense core structure using the interferometer measurements directly, in the visibility domain, without producing images. While this approach is computationally intensive, it recognizes the limitations of standard Fourier inversion and deconvolution techniques, and allows a much more direct comparison with models than analyzing images, avoiding the problems caused by synthesized beam characteristics.

The B335 visibilities are compared to theoretical models of protostellar envelope structure by constructing synthetic visibilities, taking account of (1) the dust continuum radiative transfer, and (2) the specifics of the observations, including the exact $`(u,v)`$ sampling and primary beam attenuation. The analysis of dust emission requires more complexity than dust extinction, because the models must include assumptions about the temperature distribution and specific mass opacity in addition to the model density field. Since we cannot consider an exhaustive list of protostellar envelope models, we instead fix attention on a few widely promoted density fields, including simple power law descriptions, inside-out collapse, and the Bonnor-Ebert sphere. These models all successfully match the near-infrared extinction data in the outer regions of B335, beyond $`5000`$ AU, but they predict substantial differences in the inner regions that are sampled by the interferometer.

### 3.1 Analytic Description for Power Law Distributions

The measured visibilities are given by the Fourier transform of the antenna weighted intensity distribution emergent from the globule. In our analysis, we must numerically model visibilities for comparison with the data. However, several simplifying assumptions exist that provide an analytic solution which affords some intuitive understanding of how the various physical quantities affect the visibilities.

For optically thin dust emission, the observed intensity at impact parameter $`b`$ from a spherically symmetric globule is given by

$$I_\nu (b)=2_b^{R_{\mathrm{out}}}B_\nu (T_d(r))\kappa _\nu (r)\rho (r)\frac{r}{\sqrt{r^2b^2}}𝑑r$$
(1)
where $`R_{\mathrm{out}}`$ is the outer radius of the globule, $`T_d(r)`$ the dust temperature, $`\rho (r)`$ the density, $`\kappa _\nu `$ the specific dust opacity, and $`B_\nu (T)`$ is the Planck function. If the density, and temperature follow simple radial power laws, $`\rho r^p`$, $`Tr^q`$, and the opacity is a power law in frequency that does not vary along the line-of-sight, $`\kappa _\nu \nu ^\beta `$, then if the emission is assumed to be in the Rayleigh-Jeans (R-J) regime ($`h\nu /kT<<1`$) from a globule with infinite outer boundary, the emergent intensity also has a simple power law form, i.e. $`I(\nu ,b)\nu ^{2+\beta }b^{(p+q1)}`$ (Adams 1991). If we neglect the attenuation resulting from the antenna pattern of the interferometer, then the visibility distribution is also a power law, i.e. $`V(\nu ,u)\nu ^{2+\beta }u^{(p+q3)}`$, where $`u`$ is the baseline length. With these approximations, the slope of the visibility profile — like the intensity profile — is determined solely by the sum of the power law indices assumed for the temperature and density. Note that a more steep density falloff with radius results in a less steep decrease in visibility with increasing baseline length. The relative brightness of the emission at different frequencies is determined by $`\beta `$, the dust opacity index.

### 3.2 Numerical Calculation of Model Visibilities

The synthetic visibility datasets used to compare with the observations are produced by a series of programs written in the Interactive Data Language (IDL) developed by Research Systems Inc.

The first step consists of simulating the emergent intensity distribution. Given inputs of model density and temperature distributions, and an adopted specific mass opacity, the program calculates the intensity image consisting of a $`512\times 512`$ pixel array at resolution $`0\stackrel{}{\mathrm{.}}5`$ per pixel, using the full Planck function for the emissivity, and integrates the radiative transfer equation through the model globule. These parameters determine the shape of the emergent intensity, and also the relative brightness at 1.2 mm and 3.0 mm (from the frequency dependence of the opacity). In general, each model is normalized by matching to the 1.3 mm flux measurement of Shirley et al. (2000), $`F=0.57\pm 0.09`$ Jy within a circular aperture (top-hat) of $`40^{\prime \prime }`$ diameter. The matching is achieved by iteratively scaling the product $`\kappa _\nu \rho `$ ($`\kappa _\nu `$ is the opacity, $`\rho `$ the density). We adopt an outer boundary of $`R_{\mathrm{out}}=0.15`$ pc in the models; this is well constrained both by the extinction measurements (Harvey et al. 2001) and by molecular line observations (e.g. Zhou et al. 1990). Modifying the outer boundary assumption by as much as 50% was also found to have almost no discernable effect on our results.

The next step is to simulate observations of the model intensity images. The program approximates the antenna pattern of the IRAM dishes as a Gaussian of FWHM $`50^{\prime \prime }(100`$ GHz$`/\nu )`$. Varying illumination, focus and pointing during the observations will cause the actual antenna pattern of the IRAM dishes to depart slightly from a pure Gaussian, but we expect this deviation to be small enough to have little impact on our analysis. We must account for the fact that the position of the peak millimeter emission is not exactly at the field center. To reduce the number of parameters to be fit, we fix the offset from the center as determined from a simple Gaussian fit to the 1.2 mm data: $`\mathrm{\Delta }`$R.A., $`\mathrm{\Delta }`$Dec. $`=(2\stackrel{}{\mathrm{.}}23\pm 0\stackrel{}{\mathrm{.}}08,0\stackrel{}{\mathrm{.}}79\pm 0\stackrel{}{\mathrm{.}}12`$), giving coordinates for the peak of $`\alpha =`$19:37:00.89, $`\delta =`$7:34:10.0. Experimenting with offsets $`0.1^{\prime \prime }`$ different from these values showed no discernable change in the fitting results. The program performs an inverse FFT (to choose the positive sign in the exponential), and adjusts phase to take account of the offset position. This procedure results in an array of the visibility function sampled on a regular grid of $`u`$ and $`v`$. The program then interpolates both the real and imaginary parts of the visibility to the exact $`(u,v)`$ points sampled by the IRAM PdBI observations for comparison with the data.

### 3.3 Model Selections

As detailed above, the model visibilities are derived from the assumed (1) mass density distribution $`\rho (r)`$, (2) dust temperature distribution $`T_d(r)`$, and (3) specific mass opacity of the dust $`\kappa _\nu `$. We consider the expected form of these quantities, and their uncertainties.

#### 3.3.1 Density

The ultimate goal is to constrain the density distribution of B335, given realistic choices for the temperature and opacity. The most basic models are simple power laws with radius, with the aim to constrain the power law index $`p`$ in the inner region. As an extension we also consider the physical models that arose in the Harvey et al. extinction study: (1) a Bonnor-Ebert sphere with $`\xi _{max}=12.5`$, and (2) an inside-out collapse model with $`R_{\mathrm{inf}}=0.03`$ pc ($`25^{\prime \prime }`$ at 250 pc). We also consider the effect of the bipolar outflow on the millimeter data, adopting the crude hollow cone description that successfully matched the near-infrared extinction data.

#### 3.3.2 Temperature

For a optically thin dust envelope heated by a central source, the dust temperature distribution is expected to follow a power law with index determined by the frequency dependence of the opacity, i.e. $`T_d(r)L^{q/2}r^q`$, where $`q=2/(4+\beta )`$ (Doty & Leung 1994). Breakdowns in this approximation will occur at small radii where the envelope becomes optically thick to the bulk of the short wavelength emission (or there is a circumstellar disk), and at large radii where heating from the Interstellar Radiation Field (ISRF) becomes important.

For B335, the inner breakdown of the power law approximation is expected at radii $`100`$ AU, and would manifest as a steeply increasing temperature gradient. For the IRAM PdBI observations, the longest baseline of $`60`$ k$`\lambda `$ corresponds to a linear resolution of $`900`$ AU. Hence, any optically thick region at small radii, or the presence of a disk, would appear as an unresolved component in these data. The nature and extent of the breakdown at large radii depends on the local strength of the ISRF. The modelling of Shirley et al. (2002) suggests that for typical envelope density distributions $`\rho r^p`$, with $`p1`$–2, the temperature follows very closely a power law description until falling to just under 10 K, where it flattens until at large radii external heating induces a rise to $`12`$–15 K. The detailed models of Shirley et al. (2002) for B335, using $`p=1.8`$ and an opacity index $`\beta =1.8`$, indicate that the temperature falls to about 10 K at a radius of $`5000`$ AU. For the model PdBI visibilities, any (small) variations in temperature beyond this radius have little effect because the mass density has already decreased to $`8\times 10^4`$ cm-3 (only 6% of its value at 1000 AU), and the interferometer begins to filter emission on angular scales much larger than the FWHM at 1.2 mm, $`20^{\prime \prime }`$ ($`5000`$ AU). For our standard model, we adopt the simple power law, with a minimum temperature of 10 K, i.e.

$$T_d(r)=\{\begin{array}{cc}10(5000\mathrm{AU}/r)^{0.4}\mathrm{K}\hfill & \mathrm{for}\mathrm{r}5000\mathrm{AU}\hfill \\ 10\mathrm{K}\hfill & \mathrm{for}\mathrm{r}>5000\mathrm{AU}\hfill \end{array}$$
(2)
In section 5.1.2 we investigate the uncertainties in the density profile that result from the temperature distribution in the model fitting process.

#### 3.3.3 Mass Opacity

The mass opacity of dust grains in the millimeter region of the spectrum in protostellar envelopes is uncertain, but generally assumed to follow a power law with frequency, $`\kappa _\nu \nu ^\beta `$. The power law index varies depending on the dust properties, but tends to be bounded by a small range, roughly 1–2 (Ossenkopf & Henning 1994). For our standard model, we follow Looney, Mundy & Welch (2000) and adopt $`\kappa _\nu =0.1(\nu /1200`$ GHz$`)`$ cm2 g-1, a power law with dust emissivity index $`\beta `$ equal to unity. Since the globule is optically thin at the wavelengths of the PdBI observations, especially at 3.0 mm, the power law index in the models can be adjusted to produce the best simultaneous match at the two observed wavelengths.

## 4 Method for Fitting Model Parameters and Evaluating Fit Quality

The basic procedure is to maximize the probability distribution:

$$P(\text{Model}|\text{data})=\underset{i}{}e^{(Z_if(x_i;p,m))^2/2\sigma _i^2}\times e^{(mm_0)^2/2\sigma _m^2}$$
(3)
where the $`Z_i`$ are the visibility data points with uncertainty $`\sigma _i`$, $`f(x_i;p,m)`$ are the model data points, $`p`$ a free parameter in the models, and $`m`$ a model parameter about which we have a constraint (namely that it is a Gaussian random variable with mean $`m_0`$ and standard deviation $`\sigma _m`$). Maximizing the probability distribution is equivalent to minimizing the logarithm of its inverse. Taking account the fact that the $`Z_i`$ are by nature complex visibilities, we want to minimize a modified $`\chi ^2`$:

$$\stackrel{~}{\chi }^2=\underset{i}{}\frac{|Z_if(u,v;p,m)|^2}{\sigma _i^2}+\frac{(mm_0)^2}{\sigma _m^2}$$
(4)
It is useful to be more explicit about the parameters $`p`$ & $`m`$. The free parameter $`p`$ is used to describe the shape of the model, e.g. the index of the power-law density distribution. The parameter $`m`$ allows us to include observational uncertainties, e.g. the 20% uncertainty in the normalization of the model derived by matching the 1.3 mm flux measured by Shirley et al. (2000). This is achieved by allowing the 1.2 mm model visibilties to be scaled by a constrained parameter $`m`$, a Gaussian random variable with mean 1.0 and standard deviation 20%. For the 3.0 mm model visibilities, we use two different approaches: (1) scale the 3.0 mm data by the same factor $`m`$ as the 1.2 mm data to force a dust opacity index of unity; (2) allow the 3.0 mm data to be freely scaled, to measure the dust opacity index and as a general check for consistency (by interpreting the derived value $`\beta `$ in the context of the realistic range described above).

We select a model description and perform a fit for the model parameters by minimizing the modified $`\chi ^2`$ distribution which includes data at both wavelengths. Since the models are non-linear in the fitting parameters, we analyze the uncertainty in the best-fit model parameters using a Monte Carlo technique known as the bootstrap method (Press et al. 1992, p. 691). To describe this method, consider the dataset $`S`$ used in the fitting. The first step is to construct a new dataset of equal size $`S^{}`$, by randomly selecting $`N`$ times from the original dataset. This new sample is then analyzed in the same way as the original dataset, and the fitted model parameters recorded. This process is then repeated $`n`$ times, where $`n`$ is sufficiently large that the resulting distribution of best-fit model parameters is insensitive to its exact value (for our models, typically $`n=200`$). If the fitting parameters may be successfully constrained by the data, then the distribution of best-fit parameters is very close to Gaussian, and the standard deviation of the distribution provides an estimate for the uncertainty of the parameters that best fit the original dataset.

The bootstrap method proves especially powerful for this analysis because small variations in the value of the modified $`\chi ^2`$ as defined above do not necessarily represent small variations in fit quality. The visibility dataset comprises a very large number of very low signal-to-noise measurements, which has the effect that the $`\chi ^2`$ wells are extremely shallow, and two models may seem almost equally good, despite the fact that there is ample signal to distinguish them. This occurs because the individual visibilities do not distinguish between the two models (represented in the $`\chi ^2`$ well being shallow) but the combination of all the individual visibilities does (represented by a lower $`\chi ^2`$ corresponding to a better fit). This issue could be circumvented by binning the visibilities e.g. radially, to increase the signal-to-noise, and then performing a $`\chi ^2`$ fit to the binned values. We avoid this solution because the binning loses information on the asymmetry of the globule that is hidden in the dataset by the vectorial nature of visibility data. For the spherically symmetric models, we bin the data for graphical comparison with the models.

The shallowness of the $`\chi ^2`$ wells results in the minimum values of the reduced $`\chi ^2`$ for various models to all have similar absolute values, in the range $`\chi _\nu ^2=1.4`$$`1.5`$. While a lower value of $`\chi _\nu ^2`$ may be interpreted as providing a better fit, the small range of variation confuses the differences in the fit quality. Therefore, we adopt graphical techniques to express the superiority of one model fit over another.

## 5 Results and Analysis

Table 1 summarizes the results from the $`\chi ^2`$ fits to several spherically symmetric models, as well as axisymmetric models. For fits marked ‘#a’, the normalization of the 3.0 mm visibilities is constrained by fixing the dust opacity index $`\beta `$ to unity. For fits marked ‘#b’, the dust opacity index is included as an additional fitting parameter.

### 5.1 Spherically Symmetric Models (Fits I–IV)

Although the B335 dense core shows clear departures from spherical symmetry when viewed both in dust extinction and molecular emission, for simplicity and insight we start by investigating spherically symmetric density models. In Fits Ia & Ib, we use a power law model for the density distribution. In Fit IIa, we use a power law density model but include an additional point source component to the intensity distribution. We also consider Bonnor-Ebert spheres (Fit IIIa; Bonnor 1956, Ebert 1955), and the Shu (1977) inside-out collapse model (Fit IVa).

#### 5.1.1 Single Power Laws

In the simplest model, the density distribution is described by a single power law. The best-fitting model has a power law index of $`p=1.65\pm 0.05`$ (Fit Ib) with a dust opacity index $`\beta =0.81\pm 0.12`$ that is consistent with unity. Fit Ia, which assumes $`\beta =1`$, gives a density power law index that agrees within uncertainties with that for Fit Ib. Figure 1 shows a plot of binned visibility amplitude vs. $`(u,v)`$ distance for both wavelengths. The binning is logarithmic and oversamples the data by a factor of 2; filled symbols and non-filled symbols are not completely independent. The amplitudes are the vectorial mean of the complex visibilities in each bin (corrected for the offset phase center), the error bars represent one standard deviation in the mean. For the dust emission at both wavelengths, the decrease in visibility amplitude with baseline length may be well described by a power law. The plots also show the visibility curves for the best fit model ($`p=1.65`$) and models that deviate by $`\pm 2\sigma `$ ($`p=1.55`$ and $`p=1.75`$). These models have not been scaled by the fitting parameters since this term contributes to the $`\chi ^2`$ and its exclusion allows the differing quality of each fit to be better seen (note that the $`\pm 2\sigma `$ models require a scaling at 1.2 mm of 1.19 and 0.75 for $`p=1.55`$ and $`p=1.75`$, respectively).

#### 5.1.2 Sensitivity to Temperature Assumptions

These results use the temperature distribution described in section 3.3.2, which is based on a physical argument. To understand the sensitivity of the density power law results to the temperature, we investigate how subtle quantitative changes affect the results by repeating the fitting procedure with temperature distributions that (1) have different (constant) temperature in the outer envelope, and (2) have different power law indices.

We find that the fitted density power law index is insensitive to (small) changes in the envelope temperature: variations of $`\pm 2`$ K lead to variations in the inferred density slope of much less than the quoted uncertainty. The main reason is that the antenna FWHM is only $`20^{\prime \prime }`$ at 1.2 mm (which has the higher signal-to-noise ratio of the two wavelengths), and the data are not very sensitive to the temperature structure at larger radii. Therefore, uncertainties in the strength of the ISRF are relatively unimportant in this study.

Changes in the power law index of the inner temperature distribution produce more significant changes in the fitting results. The globule envelope is optically thin at both wavelengths; if the R-J approximation were also to hold, then the inferred density power law index $`p`$ would be related to the temperature power law index $`q`$ by the relation $`p+q=`$ constant. However, the R-J approximation breaks down, especially at 1.2 mm, with the consequence that $`p`$ is more sensitive to changes in $`q`$. We find that varying $`q`$ by $`\pm 0.1`$ leads to variations in $`p`$ at the level of $`0.15`$.

#### 5.1.3 Sensitivity to an Unresolved Central Component

As already noted, there may exist a central region of B335 for which the envelope will be optically thick at infrared wavelengths, causing the temperature profile to depart from a simple power law. In this region, the temperature is expected to climb very sharply with radius, resulting in increased emission at long baselines beyond the extrapolation of a simple power law. The presence of an unresolved disk, effectively a warm central point source, would have a similar effect. While the PdBI data samples insufficient baseline lengths to directly measure the flux of any such compact component, it will, if present, flatten the visibility profiles of the various models and therefore affect the fitted density power law index. Fit IIa quantifies the effect using the combination of the simple power law for density, temperature and mass opacity from Fit Ia, but including an additional point source component with $`F\nu ^3`$ (i.e. the spectrum of a source where the majority of the emission is optically thin at 1.2 mm, with a dust opacity index of unity). For a point source flux $`F=12\pm 7`$ mJy at 1.2 mm, the most allowed by the measurements at long baselines, the fitted density index is reduced to $`p=1.47\pm 0.07`$. This point source flux corresponds to an implied disk mass of $`M2\times 10^3\mathrm{M}_{}(60`$ K$`/T_{\mathrm{disk}})`$ for our adopted opacity law. This is at the lower end of the 0.002–0.3 M range of disk masses that have been observed around T Tauri stars by Beckwith et al. 1990. Since accretion disks result from conservation of angular momentum during the collapse process, the small inferred mass of the disk may be related to the fact that B335 has a relatively low rotation rate (Frerking, Langer & Wilson 1987). However, accretion disks may also be connected with the driving source of bipolar outflows, and the B335 outflow is not unusual, with size, momentum and energy typical for a low mass protostar (Goldsmith et al. 1984).

When a compact source of emission is included in the analysis, the fitted density power law index is reduced by $`0.2`$. Thus, possible point source “contamination” is the dominant source of systematic uncertainty in the fitted density structure. This is also the case in single dish measurements, where a point source can add a systematic uncertainty of up to $`0.5`$ to the inferred density index (e.g. Shirley et al. 2002).

Table 2 summarizes various sources of systematic error and the level of uncertainty that the variations contribute to the fitted density power law index in the spherically symmetric models.

#### 5.1.4 Comparison With Previous Studies

The best fit density power law index of $`p=1.65`$ is slightly less steep than the value of $`p=1.8`$ found by Shirley et al. (2002) from SCUBA data. It is perhaps significant that the results from the two instruments are consistent to within the uncertainties considering the resolution of the JCMT at 850 $`\mu `$m is $`15^{\prime \prime }`$, roughly encompassing the entire range of radii probed by the the IRAM PdBI. However, even taking account all sources of uncertainty, the best-fit power law index of $`p=1.65`$ is significantly less steep than the $`p=2.0`$ power law index of the hydrostatic isothermal sphere that so successfully matches the outer envelope of B335 in near-infrared extinction (Harvey et al. 2001) and dust emission shape (Shirley et al. 2002) This result suggests that the density profile becomes less steep in the inner regions ($`r25^{\prime \prime }`$), a conclusion also suggested by the extinction data, where two models with this character (inside-out collapse and the Bonnor-Ebert sphere) provided better fits than a single power law.

#### 5.1.5 Broken Power Laws, Inside-Out Collapse, Bonner-Ebert Sphere

We have also investigated models with density distributions that flatten in the inner region: (1) a broken power law, with fixed index of $`p=2`$ in the envelope, with model parameters the break radius and inner power law index; (2) an inside-out collapse model with model parameter the infall radius (this is approximately a special case of the broken power law); and (3) a Bonnor-Ebert sphere, where the model parameter is the dimensionless outer radius of the sphere $`\xi _{\mathrm{max}}`$. None of these models fit the data nearly as well as a single power law.

For the broken power law model, it proved impossible to obtain robust results for the two fitted parameters, which are coupled and have comparable effects in the simulated datasets. Decreasing the break radius, or increasing the power law index in the inner region; each steepens the apparent density profile, and it is not possible to distinguish between the two effects.

The inside-out collapse model and the Bonnor-Ebert sphere that matched the extinction data both provide poor fits to the visibility data that sample the envelope closer to the protostar. Figure 2 shows a comparison between the binned visibility amplitude vs. $`(u,v)`$ distance and the curves produced from these models: the inside-out collapse model with $`R_{\mathrm{inf}}=0.03`$ pc, and the Bonnor-Ebert sphere with $`\xi _{\mathrm{max}}=12.5`$. There are large discrepancies between these models and the data. Because these models so poorly describe the data, the $`\chi ^2`$ statistic provides a poor discriminant of model parameters. At the $`\chi ^2`$ distribution minimum, the bootstrap calculation bounces between the end points of the grid at each iteration, no matter how wide a range of model parameters are considered. Essentially, the data prefer a model with a single power law index $`p1.65`$, and neither the inside-out collapse model nor the Bonner-Ebert sphere model can produce an equivalent slope over a sufficiently large enough range in spatial scales.

The inclusion of a point source component does not rectify the problems in fitting with these two physical models, despite effectively reducing the slope of the fitted density distribution. The Bonner-Ebert sphere has a relatively flat density distribution at small radii, and therefore the visibilities drop sharply at long baselines. Adding a point source to the Bonnor-Ebert sphere model adds intensity at long baselines but still cannot produce a good fit over the large range in spatial scale sampled by the data. For the inside-out collapse model, the fit fails to constrain the infall radius because this parameter does not affect the slope of the density distribution in the region where the data provide constraints (e.g. $`R_{\mathrm{inf}}0.024`$ pc $`=20^{\prime \prime }`$ at 250 pc).

A combination of the inside-out collapse model with a point source also produces a poorer fit to the data than does a simple power law with a point source. The $`\chi ^2`$ of this model is equivalent to power laws that are nearly $`3\sigma `$ away from the best fit (i.e. the fit quality is equivalent to the power law models with $`p=1.30`$ or $`p=1.65`$, when the best fit is for $`p=1.47\pm 0.07`$). This seems a surprising result, since the standard impression of the inside-out collapse model is that it behaves like a power law with index very close to 1.5 in the infalling region. But, while this is asymptotically true at very small radii, the local power law index of the inside-out collapse solution in fact decreases as one moves outwards in the infalling region, attaining a minimal value close to unity just within the infall radius. This means that for the same total flux within a $`40^{\prime \prime }`$ diameter, the inside-out model has greater intensity at larger radii than does a simple $`p=1.5`$ power law. The effect is that while the slopes of the two visibility profiles are similar in the spatial range where the PdBI data are sensitive, the primary beam attenuation (FWHM $`20^{\prime \prime }`$ at 1.2 mm) causes the inside-out model visibilities to be 25% fainter than the single power law model. The required increase in the scaling parameter, $`m`$, to best match the data is significant and contributes to the increased value of $`\chi ^2`$, resulting in a worse fit for the inside-out collapse model. Essentially, the inside-out collapse model cannot simultaneously match the Shirley et al. (2000) flux measurement and the normalization of our observed visibility profile.

### 5.2 Axisymmetric Models (Fit V)

In Figure 3, left column of panels show images of B335 made from the 1.2 mm and 3 mm visibilities, respectively. These images show that the B335 dense core is not spherically symmetric at this size scale, but rather is more like an oblate spheroid, flattened along the axis of the bipolar outflow (nearly East-West on the sky). In Figure 3, the second column of panels shows images of the best fit spherically symmetric models, made by processing the synthetic visibilities in the same way as the data. The aspect ratio of the emission is clearly different between the models and the data. The asymmetry may be caused by the action of the outflow, or it may reflect a flattening in the intial configuration of the core.

For the dust extinction data at scales $`3500`$ AU, the asymmetry was well reproduced by a model that assumed that the outflow had hollowed out a bipolar cone, with best fit constant semi-opening angle $`\alpha =41^{}\pm 4^{}`$. To extend the analysis of the visibility data, we have added the same simple representation of the bipolar outflow to the power law density distribution models to determine if this model can also successfully explain the asymmetry in the observed dust emission. In these models, we assume no unresolved point source of significant flux, having already quantified the effect of this additional component in section 5.1.3. We also assume that the temperature remains a simple radial power law. This is an extreme simplification as the emptied outflow cavities provide a less obscured view to the hotter inner region, and as a result the temperature of dust near to the boundary of the outflow cavity will likely have been underestimated.

Fit V (a & b) repeats the recipe of the $`\chi ^2`$ analysis, using a density model that comprises a radial power law with index $`p`$ and a hollow-cone outflow with constant semi-opening angle $`\alpha `$. As with Fit Ia, Fit Va assumes a dust opacity index of unity, whereas Fit Vb includes the dust opacity index as an additional parameter. As listed in Table 1, the two fits give nearly identical results, both with each other and with the outcome of Fit I. The best fit power law index (Fit Vb) is $`p=1.65\pm 0.05`$, with outflow semi-opening angle $`\alpha =40^{}\pm 5^{}`$, and dust opacity index $`\beta =0.74\pm 0.10`$.

In Figure 2, the additional columns of panels show images at both wavelengths made from the synthetic visibilities of the best fit model, and those with parameter values that differ by $`2\sigma `$ from the best fit, i.e. $`p=1.55`$, $`\alpha =30^{}`$ and $`p=1.75`$, $`\alpha =50^{}`$. The best fit spherically symmetric model is also shown to demonstrate the effect of adding the outflow cones. The quality of the match to the data is visible in (1) the extent of the apparent asymmetry, and (2) the number and spacing of the intensity contours. Note that the density gradient parameterized by $`p`$ and the degree of asymmetry parameterized by $`\alpha `$ are nearly independent. The outflow model with $`p=1.65`$ and $`\alpha =40^{}`$ provides the best match to the imaged data, and has lower reduced $`\chi ^2`$ despite the additional fitting parameter.

The best fit density power law index is essentially unchanged from the spherically symmetric analysis. That these results should be so similar is partly a consequence of the adopted outflow geometry. The main effect of the outflow cone is to change the normalization of the intensity profile along the direction through the cone, but not to change the shape or power law index of the profile. This would not be the case for an outflow model without a constant opening angle. For instance, for a model where the opening angle decreases with distance from the protostar (e.g. a parabola), the intensity profile along the outflow axis becomes less steep as relatively less material is removed from the core at larger radii. The fitted power law index in such a model would be higher than for the simple cone model. As an example, consider an outflow geometry where the semi-opening angle is $`50^{}`$ at 500 AU, decreasing to $`30^{}`$ at 5000 AU; in this model, the intensity profile along the outflow axis has power law index that is $`0.15`$ less steep than for the constant opening angle case, or the profile along an axis perpendicular to the flow. While we have not analyzed this class of model in detail, since the outflow covers less than two quadrants of the sky, spatial averaging of the regions affected by the outflow and the regions that are unaffected suggests an uncertainty of roughly $`\delta p0.15/20.08`$, a small effect on fitted power laws.

The good agreement in the fitted opening angles for the simple cone model in this work and the extinction study of Harvey et al. (2001) may not be especially significant. While the consistency is encouraging, it most likely reflects the fact that the crude outflow model creates a detectable degree of oblate flattening only for semi-opening angles near $`45^{}`$. For much largers opening angles, too much of the dense core is removed by the outflow and the core becomes disk like; for much smaller angles, the effect on the core is small, except in a small region at the poles.

### 5.3 Discussion

The analysis presented in the previous sections provides strong evidence that the density distribution power law index is less than 2.0 in the central $`5000`$ AU ($`20^{\prime \prime }`$) region where the PdBI observations are sensitive. The best fit single power law for the density distribution has index of $`p=1.65\pm 0.05`$. This value should be considered an upper limit since including a central point source — the dominant contributor of systematic uncertainty — can reduce the value to $`p=1.47\pm 0.07`$. The remaining sources of systematic uncertainty considered, of which the most important is the radial dependence of the temperature distribution, likely contribute a total uncertainty at the level $`\delta p0.2`$.

The Harvey et al. (2001) extinction study demonstrated a density index of very close to 2.0 ($`p=1.91\pm 0.07`$) over a $`3500`$$`25000`$ AU range in radius, indicative of an isothermal hydrostatic envelope, a result consistent with several different analyses of dust emission from single dish observations that probe similar spatial scales (Shirley et al. 2002, Motte & André 2001). The new interferometric dust emission data, which probes smaller scales, provides strong evidence that the density profile becomes less steep in the inner $`5000`$ AU of B335. This conclusion agrees with the hints from detailed modelling of dust extinction indicative of a turnover in the density profile at a radius of about $`6500`$ AU, as well as application of the Shu (1977) inside-out collapse model to molecular line observations (Zhou et al. 1993, Choi et al. 1995).

The inferred density profile of the B335 core conforms closely to the predictions in standard theories of isolated star formation (Shu 1977, Larson 1969, Penston, 1969) where an envelope with density power law index of $`p=2.0`$ surrounds an inner region with a $`p=1.5`$ distribution characteristic of material freely falling under the force of gravity. However, even with the contribution of a unresolved point source component, the specific Shu (1977) inside-out collapse model provides a worse fit to the data than does a simple power law. The local density index within the turnover (infall) radius in the inside-out collapse model is significantly less steep than $`p=1.5`$, ranging from roughly unity to $`1.5`$. Consequently, the inside-out collapse model cannot successfully match both the normalization of the visibility profile in the IRAM PdBI observations and the total flux within a large aperture measured with a single dish by Shirley et al. (2000). The constraint provided by the single dish measurement — effectively the visibility information at zero-spacing that cannot be measured directly by the interferometer — is crucial to this analysis; without it, the density distribution power law index can be constrained only to within $`\delta p\pm 0.3`$, not including contributions from systematic errors.

The images of B335 presented in Figure 3 demonstrate clear departures from spherical symmetry. The inner part of the globule appears significantly extended perpendicular to the outflow axis. This is in contrast to images of dust emission with lower angular resolution that show nearly spherical symmetry. While there are no robust predictions for the effect of the outflow on the density distribution, combining the power law density model with a very simple bipolar cone to describe the outflow improves the match to the data, but does not alter significantly the inferred shape of the radial variation in density.

Observations of B335 at longer baselines are needed to constrain the flux from any point source component that may be present at the center of the dense core. Such a component might be expected in the form of a circumstellar disk. Indeed, the maximal 1.2 mm flux of an unresolved component ($`12`$ mJy) lies on the low end of the range expected for the B335 disk based on the subarcsecond submillimeter fluxes for Class 0 sources observed with the CSO-JCMT interferometer (Brown et al. 2000), perhaps related to the relatively low rotation rate of B335 (see Zhou 1995).

In the future, high resolution observations of B335 at submillimeter wavelengths may help to mitigate the uncertainties associated with the temperature distribution. Spatial information on the optically thin emission distribution at a third wavelength that is sufficiently short to be outside the Rayleigh-Jeans regime for the majority of the dense core should allow for partly breaking the degeneracy between temperature, density and opacity. Such observations may be obtained with the Submillimeter Array (SMA), now under construction, or the Atacama Large Millimeter Array (ALMA).

### 5.4 A Physical Density Distribution for B335

Throughout the analysis presented so far, attention has been paid only to the power law index of the density distribution; the normalization has been absorbed in satisfying the constraint placed by the single dish 1.3 mm flux measurement of Shirley et al. (2000). To allow the results to be used in future modeling of B335, we now discuss the physical scale of the density distribution. The present study, in concert with the Harvey et al. (2001) extinction work, affords the unique position to propose a density distribution that satisfies both dust emission and extinction data over a range in radius from 500 AU to 25,000 AU.

For the adopted temperature distribution and mass opacity at 1.2 mm, the density within $`r20^{\prime \prime }`$ (5000 AU) is:

$$n_H(r)\{\begin{array}{cc}8.1\times 10^5\mathrm{cm}^3\left(\frac{r}{10^3\mathrm{AU}}\right)^{1.65}\hfill & \mathrm{for}\mathrm{no}\mathrm{point}\mathrm{source}\hfill \\ 7.6\times 10^5\mathrm{cm}^3\left(\frac{r}{10^3\mathrm{AU}}\right)^{1.47}\hfill & \mathrm{for}12\mathrm{mJy}\mathrm{point}\mathrm{source}\hfill \end{array}$$
(5)
where $`n`$ represents the total number density in the gas ($`n=n(H_2)+n(He)+`$$`=\rho /(\mu m_H)`$, $`\mu =2.29`$). Note that including the outflow as a hollow cone with semi-opening angle $`40^{}`$ increases by 14% the coefficient of the density distribution power law (and leaves the index unchanged).

For the envelope beyond $`r26^{\prime \prime }`$ (6500 AU), the Harvey et al. (2001) extinction work suggested an isothermal, hydrostatic, structure with $`r^2`$, consistent with the results of the dust emission studies of Shirley et al. (2002) and Motte & André (2001). This envelope density distribution can be adopted into our model without affecting the inferred structure in the inner regions, since emission from the envelope is a negligible contribution at small impact parameters ($`20^{\prime \prime }`$), and is filtered out at large impact parameters by the antenna pattern ($`20^{\prime \prime }`$ FWHM at 1.2 mm). However, the correct normalization of the density distribution in the outer envelope relative to the inner regions is difficult to assess because the uncertainties involved in the extinction work (gas-to-dust ratio, reddening law) are different from those that arise in the dust emission studies (specific mass opacity, temperature structure), and the previous dust emission studies do not resolve the inner regions.

The analysis of the interferometer dust continuum data allows us to draw some conclusions about the behaviour of the density profile in the transition region between the two regimes, and hence to infer a consistent scale for the envelope. In order for any model to simultaneously match the normalization of the observed visibility profile and the Shirley et al. (2000) measurement of the flux within $`r<20^{\prime \prime }`$, the local power law index in the transition region cannot decrease from the value measured in the inner region. The power law index must remain constant or become steeper. A simple continuation of the inner power law until reaching the envelope implies a density distribution for $`r>26^{\prime \prime }`$ (6500 AU):

$$n(r)\{\begin{array}{cc}1.6\times 10^4\mathrm{cm}^3\left(\frac{r}{10^4\mathrm{AU}}\right)^{2.0}\hfill & \mathrm{for}\mathrm{no}\mathrm{point}\mathrm{source}\hfill \\ 2.2\times 10^4\mathrm{cm}^3\left(\frac{r}{10^4\mathrm{AU}}\right)^{2.0}\hfill & \mathrm{for}12\mathrm{mJy}\mathrm{point}\mathrm{source}\hfill \end{array}$$
(6)
For the case of no point source, this composite density distribution gives an average power law index of $`p=1.90`$ over the 3500–25000 AU radius range sampled by the dust extinction data, precisely reproducing the extinction result ($`p=1.91\pm 0.07`$), and it gives an average power law index of $`p=1.80`$ over the 1000–25000 AU radius range, which reproduces the result of Shirley et al. (2002). Remarkably, the normalization of this envelope distribution is also within 6% of the singular isothermal sphere.

## 6 Summary

We present a study of the density distribution of the protostellar collapse candidate B335 using continuum observations at 1.2 mm and 3.0 mm made with the IRAM PdBI. In summary:

1. We perform a detailed analysis of the interferometer measurements directly in the visibility domain, not in the image domain. Though computationally intensive, this approach avoids the limitations of standard Fourier inversion and deconvolution algorithms.
2. The PdBI visibility data strongly constrain the density distribution at radii from $`500`$ AU to $`5000`$ AU from the protostar. Within this inner region, the density distribution is well described by a single power law. The best fit power law index $`p=1.65\pm 0.05`$ (1 $`\sigma `$) is significantly less steep than values close to $`p=2.0`$ derived for larger radii from analyses of dust emission and extinction. We consider the effects of various sources of systematic uncertainty on the derived value of $`p`$. Including a central source of point-like emission in the model, as might come from an unresolved accretion disk, reduces the power law index to $`p=1.47\pm 0.07`$. The remaining sources of systematic uncertainty, of which the most important is the radial dependence of the temperature distribution, likely contribute a total uncertainty at the level of $`\delta p0.2`$.
3. The inferred density profile for B335, with a power law index close to $`r^{1.5}`$ near the protostar within an $`r^2`$ envelope, matches the generic paradigm of gravitational collapse for isolated, low mass star formation. However, the specific inside-out collapse solution of Shu (1977) does not produce as good a fit to the dust emission data in the inner region as does a simple single power law, largely because of the shallow slope ($`p1`$) of the inside-out collapse solution just within the infall radius.
4. Images made from the visibility data show clear departures from spherical symmetry, with the inner part of the B335 core elongated perpendicular to the outflow axis. Including the outflow in the model fitting as a hollowed bipolar cone improves the match to the data and does not significantly change the derived radial density power law index from the value obtained in the spherically symmetric analysis. An outflow geometry where the opening angle decreases with distance from the protostar location would result in a slightly larger power law index, roughly $`\delta p0.05`$.
5. Observations that probe the subarcsecond millimeter continuum structure of B335 are needed to constrain the central point source contribution to the dust emission.

DWAH thanks Jonathan Jenkins for many useful discussions concerning statistics. We acknowledge the IRAM staff from the Plateau de Bure and from Grenoble for carrying out the observations and for their help during the data reduction. We are especially grateful to Roberto Neri for his assistance. Partial support for this work was provided by NASA Origins of Solar Systems Program Grant NAG5-6266.