text:
# Infrared Spectral Energy Distribution Model for Extremely Young Galaxies

## 1 Introduction

Galaxy formation is undoubtedly one of the desiderata in modern astrophysics, but remains a long-standing unsolved problem. The quest for ‘primeval galaxies’ is still a strong motivation of deep surveys by large telescopes and future space missions. By now, much effort has gone into the search for galaxies at high redshift ($`z`$), some of which we believe are experiencing their first episode of star formation: Lyman break objects (e.g., Madau et al., 1996; Steidel et al., 1996, 1999); Lyman $`\alpha `$ emitters (e.g., , Thompson et al.1995; Hu et al., 1998; Steidel et al., 2000); submillimetre galaxies (e.g., Hughes et al., 1998; Eales et al., 1999, 2000; Blain et al., 1999; , Barger et al.2000; Scott et al., 2002); hyper-extremely red objects (HEROs) which are also thought to be strongly enshrouded by dust (e.g., Totani et al., 2001).

However, these objects have already produced significant amount of stars, and consequently, are chemically enriched by heavy elements. Galaxies are thought to be extremely metal-poor just after their birth, hence they may have had only a small amount of dust. Even in Lyman break galaxies however, there is clear evidence that they contain non-negligible amount of dust (e.g., Adelberger & Steidel, 2000; Sawicki, 2001; Shapley et al., 2001). It remains difficult to explore the physics of metal-deficient first objects in the early Universe because of our observational limitations.

Among related issues, the extent to which very early star formation has been hidden by dust is one of the most interesting themes. Little doubt remains about the importance of extinction properties to evaluate high-$`z`$ star formation (e.g., Steidel et al., 1999). Dust grains absorb stellar light and re-emit it in the far infrared (FIR). Indeed, the FIR spectral range represents a unique opportunity to study the dust properties and distribution, and several dedicated infrared (IR) space or stratospheric missions are planned (ASTRO-F111http://www.ir.isas.ac.jp/ASTRO-F/index-e.html, SIRTF222http://sirtf.caltech.edu/, SOFIA333http://sofia.arc.nasa.gov/, Herschel Space Observatory444http://astro.estec.esa.nl/First/, etc.).

Furthermore, how dust forms and evolves in primeval galaxies needs to be considered in order to understand the chemical and thermodynamical evolution of a metal-poor interstellar medium (ISM). Dust grains drastically accelerate the formation rate of molecular hydrogen (H2), expected to be the most abundant molecule in the ISM (e.g., Glover, 2003). Hydrogen molecules emit vibrational-rotational lines. It is the only coolant in a cool (temperature $`T8000`$ K) dense core of metal-deficient system, since the atomic cooling does not work effectively. Hence, the cooling rate of a primeval galaxy consequently depends strongly on the dust content and H2 abundance.

Recently (Hirashita et al.2002, hereafter H02) have developed a model for the evolution of dust content in very young (or primeval) galaxies. Dust is expected to be produced at the final stages of stellar evolution, and Type II supernovae (SNe II) are the dominant source for the production of dust grains in young star-forming galaxies (e.g., Dwek & Scalo, 1980). Dust formation in stellar winds from evolved low-mass stars like red giant branch stars (RGBs) or asymptotic giant branch stars (AGBs) can also contribute considerably, but the cosmic time is not long enough for such stars to evolve at high-$`z`$. On the other hand, dust is also destroyed by SN shocks (e.g., , Jones et al.1996); however, H02 have shown that the dust destruction is negligible compared to the formation rate in young ($`10^6\text{}10^8`$ yr) galaxies. Then they have modeled the evolution of FIR luminosity and dust temperature in such a young starburst on the basis of SNe II grain formation model of Todini & Ferrara (2001, TF01).

While high-$`z`$ primeval galaxies are still beyond the ability of present experiments, starbursts in the Local Universe are often assumed to be good representatives of the star forming activity at high redshifts. As such, they are cosmologically interesting as key laboratories for studying the ISM, the transport of supernova processed metals, and the chemical evolution of galaxies and of the intergalactic medium (Madau, Ferrara, & Rees, 2001). H02 have paid especial attention to the nearby, extremely low-metallicity galaxy SBS 0335$``$052, since this object may be truly experiencing its first burst of star formation ($`\text{age}10^7`$ yr; Vanzi et al. 2000). The models by H02 successfully explain its observed properties such as dust mass and temperature, and far infrared (FIR) luminosity.

In this paper, we further investigate the IR properties of very young galaxies, and construct a simple model of mid- to far-infrared (MIR and FIR) SEDs of galaxies starting from the H02 models. As already mentioned above, SNe II are expected to dominate the dust formation in very young galaxies. The size of dust grains formed in SNe II cannot be as large as $`0.1\text{}1\mu `$m (TF01), and host galaxies are too young for grains to grow in the interstellar space. TF01 have shown that the size distribution of grains produced by SN II is well approximated to be single-sized. A narrow size distribution favoring small grain sizes makes the appearance of IR spectral energy distribution (SED) of young galaxies drastically different from that of nearby aged galaxies. We, for the first time, properly consider the dust size distribution peculiar to the very early stage of galaxy evolution, and construct a model of the IR SED of very young galaxies. Then, we compare our models to the observations in SBS 0335$``$052, and continue the discussion by H02 on evolution of dust content.

In order to discuss the effects of dust in such galaxies, we must first discuss the effect of dust sublimation. Dust grains are expected to sublime in H ii regions because of the strong radiation from the central source (in this case, OB stars). Indeed, the small-sized carbon grains (or PAHs) are deficient in a strong ultraviolet (UV) radiation environment (e.g., Boselli et al., 1998; Dale et al., 2001a). However, the radius of the dust cavity caused by the sublimation is much smaller than that of the star forming region for classical size grains (Inoue, 2002, and references therein). Further, as will be seen later, although our model dust grains are smaller than the classical grains, the sublimation process is governed by the energy of each injecting photon, because of the small grain cross section. For small grains, the peak temperature cannot reach the sublimation temperature. Hence we can safely neglect the effects of sublimation in our discussion.

The paper is organized as follows: We first model the dust content of a young galaxy whose age is less than $`10^8`$ yr in §2. We show the evolution of the SEDs of young galaxies and their IR extinction curve in §3. Then, we discuss some topics related to the infrared properties of young galaxies in §4. We finally compare the model predictions with the observed quantities of the extremely low-metallicity BCD, SBS 0335$``$052, in §4.1. The lowest metallicity star-forming galaxy ever known, I Zw 18 ($`1/50Z_{}`$), is another representative of the star forming dwarfs. We present a prediction for the SED of I Zw 18 in §4.2. Based on these investigations, we estimate how many low-metallicity dwarf galaxies and blue compact dwarfs will be detected in future observational projects in §4.3. The last section (§ 5) is devoted to our conclusions.

## 2 SED Model for Extremely Young Galaxies

### 2.1 Star formation and dust production

First we briefly describe the star formation and dust production model proposed by H02. This model treats the evolution of dust content in a young galaxy whose age is less than $`10^8`$ yr. Hence dust formation in stellar winds from evolved low-mass stars like RGBs or AGBs (e.g., Schmid-Burgk & Scholz, 1980; Deguchi, 1980; Woodrow & Auman, 1982; , Gail et al.1984; Gail & Sedlmayr, 1985, 1987) is negligible because the time interval is too short for stars to evolve to such stages. We can therefore safely assume that only SN II contributes to the dust formation. For similar reasons, the contribution of Type Ia supernovae to the total supernova rate can be also neglected. Further, since dust destruction by SNe may be important on a longer timescale of $`>10^8`$ yr (, Jones et al.1996; Hirashita & Ferrara, 2002), we can also safely neglect this process. Therefore, the only contributor to the total dust mass ($`M_{\mathrm{dust}}`$) in a young ($`<10^8`$ yr) galaxy is the supply from SNe II.

The rate of SNe II is given by

$`\gamma (t)={\displaystyle _{8M_{}}^{\mathrm{}}}\psi (t\tau _m)\varphi (m)𝑑m,`$ (1)
where $`\psi (t)`$ is the star formation rate (SFR) at age $`t`$ (we define $`t=0`$ at the beginning of the star formation), $`\varphi (m)`$ is the initial mass function (IMF), $`\tau _m`$ is the lifetime of a star whose mass is $`m`$, and it is assumed that only stars with $`m>8M_{}`$ produce SNe II. H02 assume a constant SFR, $`\psi =\psi _0`$, and a Salpeter IMF ($`\psi (m)m^{2.35}`$; the mass range of stars is 0.160 $`M_{}`$) to obtain a first estimate555 We assume a relatively smaller mass cut (60 $`M_{}`$) for the stellar IMF, because TF01 only give dust production rate for stars $`35M_{}`$. Though some assume that the upper mass bound $`M_{\mathrm{up}}=120M_{}`$, our framework is insensitive to the upper stellar cut-off, because the number of such high-mass ($`60M_{}`$) stars is extremely small. The contribution of the high-mass stars to UV luminosity is also much smaller than that of the intermediate ($`3\text{}40M_{}`$) mass stars if the stellar luminosity is weighted with the Salpeter IMF. In summary, the upper mass limit does not significantly affect our conclusions.. Then the rate of increase of $`M_{\mathrm{dust}}`$ is written as

$`\dot{M}_{\mathrm{dust}}=m_{\mathrm{dust}}\gamma ,`$ (2)
where $`m_{\mathrm{dust}}`$ is the typical dust mass produced in a SN II. We adopt $`m_{\mathrm{dust}}=0.4M_{}`$ (HF02). We start the calculation from $`M_{\mathrm{dust}}=0`$ at $`t=0`$. Because dust production rate is proportional to $`\gamma `$ which is proportional to the SFR $`\psi _0`$, $`M_{\mathrm{dust}}`$ calculated by our model directly scales as $`\psi _0`$. The dust mass evolution with $`\psi _0=1M_{}\mathrm{yr}^1`$ is shown in Figure 1 of H02.

### 2.2 SED construction

In this subsection, we present the construction of the SEDs for very young galaxies. As already discussed above, in their early phase, SN II explosions dominate the supply of dust. The galaxies considered here are too young to have evolved low-mass stars like RGBs or AGBs. These evolved low-mass stars are considered to be the source of large-size grains with a power-law like size distribution (e.g., Biermann & Harwitt, 1980; , Dominik et al.1989). It should also be noted that they are too young to make small grains grow larger in their diffuse interstellar medium (see, e.g., Whittet, 1992, pp.223–224).

We recall that dust grains are schematically divided into silicate and carbonaceous ones; only these two species are considered in this work. TF01 have shown that the sizes of silicate and carbonaceous grains formed in SN II ejecta to be about 10 Å and 300 Å, respectively. Then the overall size distribution of dust is well represented by a superposition of two delta functions, which correspond to the sizes of silicate and carbonaceous grains. TF01 estimated the uncertainty of the grain sizes to be a factor of $`2`$ within a reasonable range of the relevant thermodynamical parameters. Other sources of uncertainty and their estimates are discussed in TF01. Here we should keep in mind the lack of large-size grains of $`0.1\text{}1\mu `$m in very young galaxies, while a substantial contribution of such large ‘classical’ grains is inferred from the observation of the Galactic diffuse IR emission or the IR SEDs of nearby galaxies. The discrete and small grain sizes make the appearance of the SED of young galaxies drastically different from that of aged normal galaxies, especially in the IR. We next focus on the most important radiation process of small dust grains.

#### 2.2.1 Stochastic heating of very small grains

The dust grains originating from SNe II cannot be as large as those formed in the atmosphere of evolved AGB stars or grown in the diffuse interstellar space (Kozasa & Hasegawa 1987; Kozasa et al.1989; TF01). TF01 clearly showed that the sizes of the silicate grains and of the carbonaceous grains are much smaller than that of classical large grains.

Here we note that TF01 used the classical nucleation theory for the formulation of mass accretion onto grains (Feder et al., 1966). These grain sizes might be even smaller if we use a more recent formation theory of dust grains from metal-rich gas medium (e.g., Ford, 1997; Tanaka, 2002, and references therein). This revision may introduce many subtle and perplexing effects in the dust formation process, and further investigations are awaited (cf. , Tanaka et al.2002). The effect of dust size will be addressed later in the context of the origin of near-infrared (NIR) continuum radiation.

It is well accepted that very small grains are ‘stochastically heated’, that is they cannot establish thermal equilibrium with the ambient radiation field (Purcell 1976; Aannestad & Kenyon 1979; Draine & Anderson 1985; Guhathakurta & Draine 1989; Siebenmorgen et al.1992; Li & Draine 2001; Draine & Li 2001 (DL01)): the heat capacity of very small grains is too small to maintain its temperature until the next photon impinges. Hence, the equilibrium temperature of a grain cannot be well defined because its ‘temperature’ varies violently with time, i.e., it becomes very hot each time a UV photon is absorbed by the grain, and subsequently cools rapidly to a very low temperature. Consequently, the continuum shape from very small grains becomes much wider than that of a modified blackbody, and the continuum has a strong extension to the NIR.

We adopt this picture and use it to calculate the dust emission in the NIR–MIR. Now we formulate the stochastic heating of very small dust grains. The heat capacity of a grain has been extensively studied and modeled by previous papers (DL01; Li & Draine 2001, and reference therein). We apply the multidimensional Debye model for the heat capacity of a grain, $`C(T)`$, as introduced by DL01.666As the authors have already pointed out in their subsequent paper, there is a typographic error in their Equation (10) in DL01. The specific heat of silicate grains is represented by

$`C_{\mathrm{sil}}(T)=(N_{\mathrm{atom}}2)k\left[2f_2^{}\left({\displaystyle \frac{T}{500\text{K}}}\right)+f_3^{}\left({\displaystyle \frac{T}{1500\text{K}}}\right)\right],`$ (3)
and that of carbonaceous grains is

$`C_\mathrm{C}(T)=(N_\mathrm{C}2)k\left[f_2^{}\left({\displaystyle \frac{T}{863\text{K}}}\right)+2f_2^{}\left({\displaystyle \frac{T}{2504\text{K}}}\right)\right],`$ (4)
where

$`f_n(x)`$ $``$ $`n{\displaystyle _0^1}{\displaystyle \frac{y^ndy}{\mathrm{exp}(y/x)1}},`$ (5)
$`f_n^{}(x)`$ $``$ $`{\displaystyle \frac{d}{dx}}f_n(x),`$ (6)
and $`k`$ is the Boltzmann’s constant (the subscripts ‘sil’ and ‘C’ represent the silicate and carbonaceous grains, respectively). Namely, as for silicate grains, two-thirds of the vibrational modes are distributed according to a Debye model with $`n=2`$ and Debye temperature $`\mathrm{\Theta }=500`$ K, and one-third of the modes are described by a Debye model with $`n=3`$ and $`\mathrm{\Theta }=1500`$ K. For carbonaceous grains $`N_\mathrm{C}2`$ out-of-plane modes and $`2(N_\mathrm{C}2)`$ in-plane modes are set. This model excellently reproduces $`C(T)`$ of each grain species (see Figure 1 of DL01). We adopt $`\rho _\mathrm{C}=2.26[\text{g\hspace{0.17em}cm}^3]`$ (Draine & Lee, 1984) and $`\rho _{\mathrm{sil}}=3.5[\text{g\hspace{0.17em}cm}^3]`$ (Li & Draine, 2001), i.e., $`N_\mathrm{C}=1.14\times 10^{23}[\text{cm}^3]`$ and $`N_{\mathrm{atom}}=8.5\times 10^{22}[\text{cm}^3]`$, respectively ($`N_{\mathrm{atom}}`$ is the number density of atoms). Here we identify bulk graphite grain with our carbonaceous grain.

The (effective) equilibrium temperature, $`T_{\mathrm{eq}}`$ is computed by equating grain heating by photon absorption and cooling by photon emission (see Equations (6.1) and (6.2) of Draine & Lee, 1984):

$`{\displaystyle _0^{\mathrm{}}}u_\lambda cQ(a,\lambda )𝑑\lambda =`$ $`60\sigma \left({\displaystyle \frac{hc}{\pi k}}\right)^4\times `$ (7)
$`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{Q(a,\lambda )\lambda ^5}{\mathrm{exp}\left[hc/\left(\lambda kT_{\mathrm{eq}}\right)\right]1}}𝑑\lambda ,`$
where $`u_\lambda `$ is the interstellar radiation field energy density, $`a`$ is a grain size, and $`h`$, $`\sigma `$ and $`c`$ are Planck’s constant, Stefan–Boltzmann constant, and the speed of light, respectively. We adopted the values presented in Draine & Lee (1984) for the absorption efficiencies, $`Q_{\mathrm{abs}}(\lambda )`$. In practice, the above equation can be approximated by setting $`Q`$ at UV (on the lhs of Equation 7) to be unity and $`Q`$ at FIR (on the rhs) to behave like $`Q_{\mathrm{FIR}}\lambda ^2`$. Then we obtain

$`T_{\mathrm{eq}}\left({\displaystyle \frac{hc}{\pi k}}\right)\left\{{\displaystyle \frac{945u}{960\pi (2\pi Aa)hc}}\right\}^{1/6}`$ (8)
and

$`u{\displaystyle _0^{\mathrm{}}}u_\lambda 𝑑\lambda ,Q_{\mathrm{FIR}}{\displaystyle \frac{2\pi Aa}{\lambda ^2}}.`$ (9)
Here we adopt $`A_{\mathrm{sil}}=1.34\times 10^3[\text{cm}]`$ for silicate grains (Drapatz & Michel, 1977) and $`A_\mathrm{C}=3.20\times 10^3[\text{cm}]`$ for carbonaceous grains (scaled to fit $`Q_{\mathrm{IR}}(\lambda )`$ of Draine & Lee 1984).

The actual dust temperature $`T`$ has a very wide distribution around $`T_{\mathrm{eq}}`$. As we will see later, while carbonaceous grains can be in equilibrium with the ambient radiation field, silicate grains cannot establish the equilibrium, even in the intense radiation field (high photon density) of a young galaxy (§3 and §4).

#### 2.2.2 Emission

Following Draine & Anderson (1985), with this $`C(T)`$, we calculated the temperature distribution of dust as a function of size $`a`$ with Monte Carlo simulations for incident UV photons.

We calculated $`u_\lambda `$ from the history of OB star luminosity of H02 and the size of the star forming region $`r_{\mathrm{SF}}`$. We use the OB star luminosity calculated in H02 as the total luminosity of stars with $`m>3M_{}`$. We start from $`r_{\mathrm{SF}}=30`$ pc as the size of a star-forming region in young galaxies (H02), but other sizes are also examined later.

We approximated the spectrum of the UV field by the formula of Draine & Bertoldi (1996) for simplicity and scale it with luminosity by factor $`\chi (t)`$ as

$`u_\lambda =`$ $`\left[{\displaystyle \frac{25}{6}}\left({\displaystyle \frac{10^5\times \lambda }{1\text{cm}}}\right)^2+{\displaystyle \frac{25}{2}}\left({\displaystyle \frac{10^5\times \lambda }{1\text{cm}}}\right){\displaystyle \frac{13}{3}}\right]`$ (10)
$`\times 10^9\chi (t)[\text{erg\hspace{0.17em}cm}^3\text{cm}^1].`$
The UV radiation field strength of the solar neighborhood, $`u^{\mathrm{NBH}}`$, is $`u^{\mathrm{NBH}}=4.0\times 10^{14}[\text{erg\hspace{0.17em}cm}^3]`$ (Draine & Bertoldi, 1996). The intensity factor $`\chi (t)`$ is determined as a function of the galaxy age from the OB luminosity history and $`r_{\mathrm{SF}}`$ (Equation 11) 777We note that, for a high-$`z`$ galaxy, we should include CMB in the radiation field. Small grains spend most of their time with their lowest temperature, which is determined by CMB. This slightly affect the submillimetre SED of high-$`z`$ galaxies. An extensive discussion on the heating of dust by CMB is given in Ferrara et al. (1999a)..

For the geometry, we assume that the OB stars are concentrated in the centre of the system, and dust is assumed to be distributed as a shell, so as to express the fact that dust absorbs the UV photons very effectively in star forming regions (, Hunt et al.2001). Moreover, as seen later, such a “simplest” geometry is successful enough to explain the SED of SBS 0335$``$052. We regard the radius of the sphere, $`r`$ as the size of the spherical star forming region, $`r_{\mathrm{SF}}`$. Then, the flux is calculated by $`J(t)=L_{\mathrm{OB}}(t)/(4\pi r_{\mathrm{SF}}^2)`$. Then the UV energy density incident on dust grains, $`u`$, is

$`u={\displaystyle \frac{J(t)}{c}}=\chi (t)u^{\mathrm{NBH}}.`$ (11)
The spectral UV energy density $`u_\lambda `$ is calculated from the assumed spectrum of the UV radiation field (Equation 10).

Such a shell-like geometry of dust distribution may result from the drift induced by the radiation pressure. Gail & Sedlmayr (1979) show that grains in H ii regions can be accelerated nearly to $`v_\mathrm{d}10\text{km\hspace{0.17em}s}^1`$. The crossing time $`t_{\mathrm{cross}}`$ can be estimated as

$`t_{\mathrm{cross}}3\times 10^6\left({\displaystyle \frac{r_{\mathrm{SF}}}{30\mathrm{pc}}}\right)\left({\displaystyle \frac{v_\mathrm{d}}{10\mathrm{km}\mathrm{s}^1}}\right)^1[\text{yr}].`$ (12)
Therefore, the dust grains can be swept into a shell surrounding the star-forming region within a timescale of $`10^6\text{yr}`$. Stellar winds can also sweep dust grains (see Inoue, 2002, and references therein for the creation of central dust cavity). We will go back to this issue in line with the application of our model to the local low-metallicity star forming dwarf galaxies.

The rate at which a grain absorbs a photon with energy $`EE+dE`$ is expressed as

$`{\displaystyle \frac{d^2p}{dEdt}}=Q_{\mathrm{abs}}(a,\lambda )\pi a^2u_\lambda {\displaystyle \frac{\lambda }{h^2c}}.`$ (13)
The heating is represented as follows:

$`{\displaystyle \frac{hc}{\lambda }}={\displaystyle \frac{4\pi }{3}}{\displaystyle _{T_0}^T}C(T^{})𝑑T^{},`$ (14)
where $`T`$ is the peak temperature achieved by a grain hit by a photon with energy $`hc/\lambda `$, and $`T_0`$ is the grain temperature just before absorption. If more photons are absorbed by a grain simultaneously, the lhs will become the total energy of the photons. On the other hand, the grain cools through radiation as

$`{\displaystyle \frac{dT}{dt}}={\displaystyle \frac{3\pi }{aC(T)}}{\displaystyle _0^{\mathrm{}}}Q_{\mathrm{abs}}(a,\lambda )B_\lambda (T)𝑑\lambda .`$ (15)
From Equations (13), (14), and (15), we obtain the sample path of a grain temperature as a function of time. A sufficiently large set of Monte Carlo realizations yields the dust temperature distribution (Draine & Anderson, 1985). Though there could be faster algorithms (Guhathakurta & Draine, 1989), we adopt this direct method in order to examine the cooling behavior and timescale of very small dust grains in an extremely intense radiation field of a young galaxy.

The total mass of each grain component is given by TF01. The mass ratio we adopt here is $`M_{\mathrm{sil}}:M_\mathrm{C}=0.56:0.44`$ (cf. H02). With this value and material density of each species ($`\rho _{\mathrm{sil}}=3.50\text{g}\text{cm}^3`$ and $`\rho _\mathrm{C}=2.26\text{g}\text{cm}^3`$), we obtain total numbers of grains, $`N_i`$

$`N_i={\displaystyle \frac{3M_{\mathrm{dust}}f_i}{4\pi a_i^3\rho _i}},`$ (16)
where subscript $`i`$ denotes the species of dust, sil or C, and $`f_i`$ is the mass fraction of dust of species $`i`$.

Monochromatic luminosity of a grain with radius $`a`$ and species $`i`$, $`L_{\mathrm{IR}}^{\mathrm{grain},i}(\lambda )`$, is, then, obtained by

$`L_{\mathrm{IR}}^{\mathrm{grain},i}(\lambda )=4\pi a_i^2\pi {\displaystyle Q_{\mathrm{abs}}^i(\lambda )B_\lambda (T)\frac{dP_i(\chi ,a_i)}{dT}𝑑T},`$ (17)
where $`dP_i(\chi ,a_i)/dT`$ is the temperature probability density function. In the present case, $`a_i`$ is fixed for each dust component to a specific value in the following range: $`a_{\mathrm{sil}}=6\text{}10\text{Å}=6\text{}10\times 10^8\text{cm}`$ and $`a_\mathrm{C}=200\text{}300\text{Å}=2\text{}3\times 10^6\text{cm}`$, i.e., the dust size distribution function is well approximated by a Dirac’s delta function888 The results will not be affected by the detailed choice of the functional form: e.g., Gaussian with small standard deviation, instead of delta function, changes the results little. The discrete nature of the size distribution is fundamental here, i.e., slowly varying functions such as a power-law yield significantly different results in contrast.  for each species, $`i`$,

$`{\displaystyle \frac{dn_i}{da}}(a)=n_i\delta (aa_i).`$ (18)
Total IR luminosity from species $`i`$, $`L_{\mathrm{IR}}^i`$, is

$`L_{\mathrm{IR}}^i`$ $`=`$ $`N_i{\displaystyle L_{\mathrm{IR}}^{\mathrm{grain},i}(\lambda )𝑑\lambda }`$ (19)
$`=`$ $`4\pi a_i^2N_i\pi {\displaystyle Q_{\mathrm{abs}}^i(\lambda )B_\lambda (T)\frac{dP_i(\chi ,a_i)}{dT}𝑑\lambda 𝑑T},`$
and we have imposed the constraint $`L_{\mathrm{IR}}=_iL_{\mathrm{IR}}^iL_{\mathrm{OB}}(t)`$; if $`L_{\mathrm{IR}}`$ exceed $`L_{\mathrm{OB}}`$ we set $`L_{\mathrm{IR}}=L_{\mathrm{OB}}`$. Final continuum IR emission from a young galaxy is calculated by superposing each of the continuum from silicate and carbonaceous grains as a function of starburst age, $`t`$,

$`L_{\mathrm{IR},\lambda }(t)`$ $`=`$ $`{\displaystyle \underset{i}{}}L_{\mathrm{IR}}^i(\lambda )`$ (20)
$`=`$ $`{\displaystyle \underset{i}{}}4\pi a_i^2N_i\pi {\displaystyle Q_{\mathrm{abs}}^i(\lambda )B_\lambda (T)\frac{dP_i(\chi ,a_i)}{dT}𝑑T},`$
where in the rhs starburst age is incorporated with $`N_i`$ and $`dP_i(\chi ,a_i)/dT`$ through $`M_{\mathrm{dust}}(t)`$ and $`\chi (t)`$.

#### 2.2.3 Extinction

If the dust opacity is very large, then self-absorption occurs, and even MIR radiation from dust is absorbed by the dust itself. In this work we did not try to solve the detailed radiative transfer but rather treat the extinction as given by a one-zone screen model. The resulting extinction curve is different from that of Draine & Lee (1984) because of the different dust size distribution. In this case dust opacity and extinction are expressed as

$`\tau _{\mathrm{dust}}(\lambda )`$ $``$ $`{\displaystyle \underset{i}{}}\pi a_i^2Q_{\mathrm{abs},i}(\lambda )n_ir_{\mathrm{SF}},`$ (21)
$`I(\lambda )`$ $`=`$ $`I_0(\lambda )e^{\tau _{\mathrm{dust}}(\lambda )}.`$ (22)
The absorbed light is re-emitted at longer wavelengths, mainly in the submillimetre, and consequently, the SED is deformed by the self-absorption. The final SED is obtained via this self-absorption–re-emission process with an obvious constraint of the energy emitted in the IR, $`L_{\mathrm{IR}}`$, which must be conserved.

The screen geometry might seem an oversimplification, but many observations support screen-like dust configuration in intense starbursts (e.g., Meurer et al., 1995). Moreover, as discussed in §2.2.2, dust can be swept into a shell surrounding the star-forming region. The possible clumpiness of dust may also increase the dust opacity $`\tau _\lambda `$ of the galaxy (for a discussion see Ferrara et al., 1999b), as well as it may allow the UV photons to leak out of the galaxy. Although such an anisotropy is difficult to model, it is remarkable that the simple geometry we adopted can reproduce a number of the observed properties of low metallicity star forming dwarfs as we will see later. This aspect will be further considered in §4.1.2.

## 3 Results

### 3.1 Infrared SED evolution

We present the results of the evolution of IR SEDs for a very young galaxy with the age of $`10^6\text{}10^8`$ yr in Figures 1 and 2. To illustrate the results, and facilitate comparison with H02, we concentrate on the value $`\text{SFR}=1M_{}\text{yr}^1`$. We set the silicate dust grain size to $`a_{\mathrm{sil}}=10`$ Å and $`6`$Å (solid and dotted lines, respectively), because grain size can be changed by a factor of about two depending on the conditions of grain formation. The size of carbonaceous grains can also be altered for the same reason, but it affects the SED very little because the grains are in a thermal equilibrium. Thus, we only show the results for $`a_\mathrm{C}=300`$Å in Figures 1 and 2. In the present model, the radiation field is typically $`10^{4\text{}5}`$ times stronger than in the solar neighborhood.

In the case of smaller silicate dust size, continuum radiation produced by stochastically heated grain extends toward NIR wavelengths. This effect is quite sensitive to the grain size of $`a10\text{Å}`$, in our case that of silicate grains. In contrast, the so-called ‘very small grains (VSGs)’ in the local normal (viz. solar metal abundance) star forming galaxies are generally thought to consist of carbonaceous grains or their fragments. In the early, optically thin phase of a young galaxy, these small silicate grains produce a prominent band feature around $`9.7\mu `$m in emission. It is worth mentioning that even in the extremely intense radiation field of $`\chi =\text{a few}\times 10^4\text{}10^5`$, such a small grain cannot establish thermal equilibrium with the ambient radiation, because the cooling time of a grain (typically $`t_{\mathrm{cool}}10^2`$ s at $`T500`$ K) is much shorter than the mean interval of photon injection onto the grain ($`\mathrm{\Delta }t10\text{}10^3`$s).

Plante & Sauvage (2002) suggested that the absence of unidentified infrared band (UIB) features in the intense UV environment of SBS 0335$``$052 implies that stochastically heated small grains are unimportant in such galaxies. However, as we have explained, small silicate grains considered here may be uncorrelated with UIBs. Small silicate grain is naturally expected in SNe (TF01), so we should properly take into account the effect of stochastic heating of small grains in starburst systems. Indeed, until the extinction becomes significant, their contribution is dominant in N–MIR wavelengths in our results.

We also stress that the dust temperature of carbonaceous grains is very high, $`T>100`$ K in the early phase. Because of the lack of classical size ($`0.1\text{}1\mu `$m) large grains, submillimetre continuum radiation is weak, and the peaks are located around the MIR during the first a few tens of Myrs. Such a high dust temperature is also predicted for giant forming ellipticals by (Totani & Takeuchi, 2002).

### 3.2 Opacity evolution

After the early phase characterized by the radiation from hot dust ($`\text{age}10^7`$ yr), opacity of dust becomes large, and consequently self absorption of dust gets quickly severe in the N–MIR regime. Finally, the observed dust temperature decreases to 50 K or less. In particular, the $`9.7\mu `$m silicate band causes a prominent absorption feature even for a young age $`t10^7`$ yr. The absorbed light is then re-emitted in the FIR-submm wavelengths, which results in a cooler effective temperature. Consequently, the temperature history is consistent with that of H02 which was derived simply from energy considerations.

We show the evolution of the IR extinction curves for the cases $`r_{\mathrm{SF}}=30`$ pc and 100 pc in Figure 3. We should note that, since dust is supplied only from SNe II (i.e., the dust size distribution does not vary) within the timescale considered here, shape of the extinction curve does not change in time. Hence, strictly speaking, Figure 3 represents the evolution of dust opacity, rather than the evolution of extinction curve.

By our definition, $`A_\lambda `$ is independent of the grain size $`a_i`$, since, from Equations (16) and (21), $`A_\lambda `$ is obtained as

$`A_\lambda =1.086\tau _{\mathrm{dust}}(\lambda )\left({\displaystyle \frac{M_{\mathrm{dust}}f_i}{\rho _ir_{\mathrm{SF}}^2}}\right){\displaystyle \frac{Q_{\mathrm{abs},i}}{a_i}},`$ (23)
and $`Q_{\mathrm{abs},i}`$ scales as $`Q_{\mathrm{abs},i}a_i`$. Hence, results are the same for $`a_i=6`$ Å and 10 Å.

Here we should keep in mind that the large opacity makes radiation pressure effective (Ferrara et al., 1991), and may cause the expansion of the star-forming complex and eventually, the onset of galactic wind at some stage in the galaxy evolution. Then the opacity will become smaller than shown in Figure 3. In order to take it into account, we should construct a consistent framework which can treat dynamical and chemical evolution as well as radiation. We regard dynamical modeling as the next stage and do not go further in this work. Some comprehensive discussions on galactic winds in the similar condition can be found in, e.g., Plante & Sauvage (2002). More general properties of the wind in forming dwarf galaxies are discussed in, e.g., Ferrara & Tolstoy (2000).

## 4 Discussion

### 4.1 A starbursting dwarf SBS 0335$``$052

#### 4.1.1 Observed SED

Here we adopt our model to the local BCD, SBS 0335$``$052 with its unusual IR SED and strong flux at MIR in spite of its low metallicity ($`1/41Z_{}`$). Since (Hunt et al.2001) estimated the SFR of this galaxy to be about $`1.7\text{M}_{}\text{yr}^1`$ from $`\text{Br}\alpha `$ line intensity, we scale the OB star luminosity of H02 with this value. The result is shown in Figure 4. The photometric data presented in filled and open squares in Figure 4 are mainly taken from Table 1 of Dale et al. (2001b). The points at $`10.8`$ and $`21\mu `$m obtained by Gemini/OSCIR and the data (and more strict upper limits) at 60, 65, and $`100\mu `$m obtained by ISOPHOT are from Table 1 of Plante & Sauvage (2002). The ISOCAM spectrum (blue curve) is taken from (Thuan et al.1999). In the NIR, we should correct the contribution of stellar light and nebular continuum. By assuming the stellar light fractions to be $`37`$ % ($`2\mu `$m) and $`6`$ % ($`4\mu `$m) and applying the nebular emission coefficients of Joy & Lester (1988), (Hunt et al.2001) obtained the net dust contribution of $`13`$ % and $`67`$ % at $`2\mu `$m and $`4\mu `$m, respectively. We show these values by open triangles in Figure 4. In the left panel, the data are converted to the monochromatic luminosity $`\nu L_\nu [L_{}]`$, while in the right panel we show in the unit of the flux density in Figure 4.

From the NIR spectrum, the age of the stellar population is estimated to be $`5`$ Myr (Vanzi et al., 2000). In order to reproduce an extreme NIR–FIR SED of SBS 0335$``$052 within the strict age limit, stellar energy radiated in the optical wavelengths must be effectively converted to FIR radiation. It suggests a very large dust opacity, requiring a very compact SF region ($`r_{\mathrm{SF}}=20`$ pc).

The small size of the SF region is also consistent with the observed high dust temperature $`80`$ K. SBS 0335$``$052 has a strong N–MIR continuum without UIB or PAH features (, Thuan et al.1999). In our framework, this is consistent with a small silicate grain size $`a_{\mathrm{sil}}=6`$ Å. Further, the MIR peak and lack of a prominent submillimetre component strongly support the absence of classical-size large grains. The size of the carbonaceous grains is suggested to be $`a_\mathrm{C}=200\text{Å}`$ in order to reproduce the ISOPHOT observation of Plante & Sauvage (2002). Plante & Sauvage (2002) also reached to a similar conclusion for the lack of classical size grains from a different radiative transfer model.

A significant absorption feature at $`\lambda 10\mu `$m suggests a large dust opacity. It is a natural consequence of the large IR luminosity of this galaxy, and is also consistent with visual extinction of $`A_V=12\text{}20`$ mag for the embedded star forming region (, Thuan et al.1999, Hunt et al.2001). This aspect will be discussed in §4.1.2.

As we mentioned in §2.2.2, the configuration of dust may be shell-like within a timescale of $`10^6`$ yr. The estimated age of the starburst of SBS 0335$``$052 is comparable to the crossing timescale of the system.

With these physical values and under our assumptions, we have succeeded in reproducing the overall SED of SBS 0335$``$052 with the starburst age $`t6.5`$ Myr. Thus, all the peculiar properties of SBS 0335$``$052 are consistent with those of a genuine young galaxy in the Local Universe. This success strongly supports the validity of our model, despite its simplicity.

#### 4.1.2 Extinction curve

(Hunt et al.2001) have shown that the extinction properties of SBS 0335$``$052 can be well described by the extinction curve of Lutz (1999), advocated for the extinction of the Galactic Centre (GC), assuming a simplifying assumption of screen-like dust geometry. In their evaluation of extinction, Hunt et al. assumed a superposition of two modified blackbodies for the continuum radiation, only the cooler of which was extinguished by the screen. (Thuan et al.1999) and Dale et al. (2001b) also assumed modified blackbody for the unextincted continuum to explore the properties of the IR radiation of this galaxy.

However, as we have already seen in the previous section, we cannot expect equilibrium between the UV radiation field and very small dust grains, which are thought to be the main source of the continuum radiation in the N–MIR. This is also pointed out by Dale et al. (2001b). Bearing this in mind, here we reconsider the IR extinction properties.

The IR extinction curve predicted by our model is presented in Figure 5. This curve is the one of a galaxy with the same parameters used in Figure 4.

First, it is roughly consistent with the Lutz (1999) extinction curve, scaled by a visual extinction $`A_V16.3`$ mag (, Hunt et al.2001). Indeed, the GC curve is quantitatively well reproduced by our model with $`\text{SFR}=1.7\text{M}_{}\text{yr}^1`$, the star forming region size $`r_{\mathrm{SF}}20`$ pc, and the above starburst age. Since our model provides the extinction curve only at wavelengths $`\lambda >1\mu `$m, if we extrapolate the model extinction of $`K`$-band to $`V`$-band by the ratio tabulated in Rieke & Lebofsky (1985), we obtain the visual extinction strength $`A_V18`$ mag. It is also equal to the value calculated by the relation $`A_K=\mathrm{\hspace{0.17em}0.11}A_V`$ (, Cardelli et al.1989). This value is in good agreement with the estimate of Hunt et al., although lower than the value $`A_V=30`$mag by Plante & Sauvage (2002). This may be caused by the difference of the radiation models used for calculation. Here we should keep in mind that the $`A_V`$ measured at the optically visible region of SBS 0335$``$052 is small ($`A_V0.55`$ mag; e.g., Izotov et al., 1997), typical of dwarf irregulars. So, the peculiar characteristics of SBS 0335$``$052 may be governed by the heavily embedded star-forming region.

Though the model extinction curve is roughly consistent with the GC one in the NIR regime, our model predicts a flatter slope than the GC one in the extrapolation to the optical. Our extinction curve in the NIR regime has a slope similar to (Cardelli et al.1989) ($`\lambda ^{1.6}`$). (Hunt et al.2001) found that the extrapolation of Lutz’s GC curve to the optical by the slope of $`\lambda ^{1.75}`$, close to that of Seaton (1979) (see the extension of the dashed line to the optical in Figure 5). The flatter slope of our extinction curve is a natural consequence of the fact that the continuum radiation is not produced from the dust in equilibrium with ambient UV field. The shape of the non-equilibrium continuum is much broader than that of modified blackbody radiation. The intensity of our continuum radiation in the N–MIR is stronger than that of the modified blackbody if we fix the intensity at about $`20\mu `$m, hence, we need a stronger extinction (flatter slope) at N–MIR wavelength range. The stronger extinction at the $`9.7\mu `$m broad silicate feature is also explained in a similar way. Such behaviour is required to reproduce the absorption feature in the observed SED of SBS 0335$``$052.

In order to examine these issues in more detail, we estimated the extinction with our ‘bare’ model SED, calculated without taking into account the self-absorption, rather than using the modified blackbody continuum. The deduced extinction is represented by open squares and triangles in Figure 5. We call the newly derived extinction the ‘semi-empirical extinction curve’. Unfortunately, since the NIR fluxes from dust emission contain a large uncertainty depending on the nebular and stellar contribution, we cannot strongly constrain the shape of the curve, but the NIR part of our model extinction curve shows rough agreement with the semi-empirical estimates. Around the $`9.7\mu `$m silicate feature, the model extinction curve shows very good agreement with the semi-empirical one. This result vindicates the performance of our model, as well as the simplifying assumption of a screen-like dust geometry.

However, around the $`18\mu `$m silicate feature, our model significantly underestimates the extinction. Correspondingly, our model SED overproduces the flux around $`20\mu `$m, compared with the observation of Plante & Sauvage (Figure 4). Many IR galaxies show the strong absorption feature at $`18\mu `$m (e.g., Dudley & Wynn-Williams, 1997). The discrepancy stems from the shape of the extinction curve around $`20\mu `$m. We cannot reproduce such a feature within our present framework based on the dust emissivity of Draine & Lee (1984). This riddle will be investigated in future work.

Except the above discrepancy, our model extinction curve naturally explains most important properties of the GC curve of Lutz (1999). The success of our model in reproducing the extinction curve strongly supports the composition of dust supplied by SNe II as predicted by TF01.

### 4.2 Prediction for I Zw 18

I Zw 18 has long been a record holder as the lowest metallicity star-forming galaxy in the Local Universe ($`Z=1/50Z_{}`$), as well as a textbook example of a dwarf starburst. Hence, in the context of our model, it is an intriguing object to discuss as another representative of extremely low-metallicity galaxies. Östlin (2000) estimated its distance to be 12.6 Mpc, and we use this value for our discussion. Whether I Zw 18 is really a newly formed system or not has been a matter of debate. The discovery of AGB stars (Östlin, 2000) may suggest an underlying old population, but its contribution is almost certainly not dominant (, Hunt et al.2003). For this reason, we therefore regard I Zw 18 as a forming galaxy for the first approximation.

Recently, from HST WFPC2 narrowband imaging, Cannon et al. (2002) derived various important physical properties of this galaxy. They showed that the O8 ionizing equivalent number of this galaxy is 2886 for the total galaxy. It corresponds to a total OB star luminosity, $`L_{\mathrm{OB}}4.9\times 10^8L_{}`$. From the H$`\alpha `$ image, we estimate the total size of the SF regions of $`r_{\mathrm{SF}}50\text{}100`$ pc. The SFR of I Zw 18 is $`0.04M_{}\text{yr}^1`$, and for some patches, $`A_V0.5`$ mag (Cannon et al., 2002). We use these values to constrain the model prediction.

From the measured extinction, I Zw 18 may also contain a significant amount of dust. Cannon et al. (2002) estimated the total dust mass $`M_{\mathrm{dust}}2\text{}5\times 10^3M_{}`$. This value is consistent with the dust-to-gas ratio expected from the metallicity (Lisenfeld & Ferrara, 1998; Cannon et al., 2002). The dust mass of our model reaches this value at the starburst age $`t10^{7\text{}7.5}`$ yr, an age consistent with the observationally suggested major SF age in I Zw 18 (e.g., Martin, 1996; , Hunt et al.2003). Here, we note that a smaller size of the SF region makes $`A_V`$ exceed 0.5 mag and is not permitted.

We show the resulting SED in Figure 6. The photometric data in the NIR are taken from (Hunt et al.2003). For the case of I Zw 18, interstellar conditions are quite different from those in SBS 0335$``$052, in spite of their similar metallicities. Since the visual extinction $`A_V`$ is very low even in its star forming patches, the radiation reprocessing from UV to FIR is inefficient. It naturally leads to a low IR luminosity $`L_{\mathrm{IR}}10^{7.5}L_{}`$. For the same reason, the silicate band features are observed in emission in the SED of I Zw 18.

It is worthwhile to focus on the differences between I Zw 18 and SBS 0335$``$052. Indeed, these two representative low-metallicity star-forming dwarfs are significantly different in spite of their almost same metallicities. First, the SFR of I Zw 18 is much smaller than that of SBS 0335$``$052. This means that the dust production rate of I Zw 18 is much smaller than that of SBS 0335$``$052. In spite of their optical similarity, I Zw 18 has rather normal SFR while SBS 0335$``$052 has extremely vigorous star formation. Such an active/passive star formation variety may play an important role in understanding their activity (, Hunt et al.2002). Second, the apparent SF region size of I Zw 18 is much larger than that of SBS 0335$``$052. This difference could be produced by the dust sweeping by the radiation pressure from the central stellar cluster. When the age of the system becomes older, the dust layer may be swept far away from the star-forming complex, and the opacity may get smaller. Since I Zw 18 has an older age than SBS 0335$``$052, the dust layer of the I Zw 18 can be transported further away. The small dust opacity of I Zw 18 might be understood by the combination of the above two effects.

Inspection of Figure 6 tells us that the peak of the SED of I Zw 18 is located around $`100\mu `$m both in the most optimistic and pessimistic predictions. However, the detection in these wavelengths is limited by Galactic cirrus. A shorter wavelength ($`60\mu `$m) observation will be dramatically easier. The MIR observation from space seems exciting, and crucial to examine the validity of our prediction. In the NIR, dust contribution may be negligible in the NIR. The fraction of dust emission appears to be smaller than that of SBS 0335$``$052. This suggests that either of stellar and nebular contribution is dominant in the NIR.

We discuss the possibility and prospects for future observation in the next subsection.

### 4.3 Prospects for future observations

As we have seen in the above, even in the intense UV radiation field of very young galaxies, small silicate grains suffer from stochastic heating and consequently have a broad temperature distribution. Hence, silicate grains mainly contribute to the MIR continuum emission. Carbonaceous grains have larger sizes and can be in thermal equilibrium, but the strong UV field makes their equilibrium temperature very high ($`50\text{}100`$ K), and they also contribute to the MIR. These facts indicate that both (restframe) MIR and FIR surveys are necessary to examine the dust amount of galaxies over wide range of metallicity.

Our models show that the observed flux of I Zw 18 is expected to be $`10\text{}150`$mJy at $`100\mu `$m. This agrees very well with a simple argument by Kamaya & Hirashita (2001), who showed that the expected flux for I Zw 18 is 7–100 mJy at $`100\mu `$m. Typical detection limits for the forthcoming FIR surveys are $`1\text{}100`$mJy at $`80\mu `$m and $`10\text{}100`$mJy at $`150\mu `$m. We expect $`\text{several}\times 10^2`$ detections of such local low metallicity galaxies by some space observations in the near future such as the SIRTF SWIRE Legacy program and ASTRO-F FIR all-sky survey, if we assume the local luminosity functions for star-forming dwarf galaxies. However we should be aware of the cirrus confusion limit, which makes it difficult to detect IR extragalactic sources at wavelengths $`\lambda >100\mu `$m. In order to detect these galaxies at $`\lambda =100\text{}200\mu `$m, a large aperture ($`\varphi 2\text{}4`$m) FIR space telescope like Herschel is required to obtain a small diffraction limit to overcome the cirrus confusion limits. Herschel Space Observatory will be one of the optimal instruments for this purpose (on the estimation of the FIR confusion limit for various forthcoming instruments, see e.g., Ishii et al.2002). In any case, such a sample will be useful for the systematic studies on ‘primeval galaxies’, through the investigation of other possibly ‘young’ local dwarfs (e.g., , Hopkins et al.2002; Corbin & Vacca, 2002).

Furthermore, our result for hot dust also offers a potentially important cosmological insight into high-$`z`$ galaxies: if dust in high-$`z`$ object is as cool as usually assumed, then we cannot expect a high comoving SFR density because in a realistic cosmology the cosmic IR background (CIRB) spectrum strongly constrains the high-$`z`$ IR emission (Takeuchi et al., 2001a). However, very hot dust in such high-$`z`$ object reconciles the CIRB constraint with high SFR, and allows vigorous star formation hidden by dust in the early Universe (Totani & Takeuchi, 2002). Some studies attempt to predict the number of submillimetre galaxies at high-$`z`$ ($`z=2\text{}3`$) from both empirical and theoretical viewpoints (e.g., Guiderdoni et al., 1998; Tan et al., 1999; Magliocchetti et al., 2001; Takeuchi et al., 2001b). Most models assume galaxy SEDs similar to those of local dusty giant starbursts like M82 or Arp 220. However, we advocate our theoretical framework for very young galaxies to make it possible to treat dwarf starbursts, because present galaxy formation models suggest a vast number of dwarf galaxies in the early Universe. Hirashita & Ferrara (2002) provided a comprehensive framework of very high-$`z`$ (or, say, primeval: $`z5`$) dwarfs in the submillimetre and millimetre, and showed the ALMA detectability. They also suggest the existence of dust warmer than the typical Galactic environment because of the intense radiation field. This emphasizes the importance of our SED modeling which is applicable to young galaxies with intense interstellar radiation fields. With our SED modeling, thus, the advent of ALMA will enable us to compare our models with observations of such low mass, bona fide primeval galaxies.

## 5 Conclusions

We have constructed a simple model of mid- to far-infrared (MIR and FIR) SEDs of galaxies to investigate the IR properties of very young galaxies, based on the galaxy evolution model of (Hirashita et al.2002, H02).

Small grains cannot establish thermal equilibrium with the ambient radiation field because of their very small heat capacity. We find that, even in the intense UV radiation field of very young galaxies, small silicate grains still cannot be in equilibrium, and consequently have a broad temperature distribution. Hence, silicate grains contribute mainly to the MIR continuum emission. Carbonaceous grains are larger and can be in thermal equilibrium, but the strong UV field makes their equilibrium temperature very high ($`50\text{}100`$ K); hence they also contribute to the MIR.

We calculated the SEDs of very young, low-metallicity galaxies as a function of starburst age. In the first several Myrs, such a galaxy shows very hot dust emission with a peak at $`\lambda 30\text{}50\mu `$m, because the self absorption is not very pronounced in this phase. Later, the SED changes its appearance to a more common shape, dominated by cooler dust emission, caused by the strong self-absorption of MIR light and re-emission at submillimetre wavelengths. The evolution of the dust opacity and the shape of the IR extinction curve were also given.

We succeeded in reproducing the MIR–FIR SED of SBS 0335$``$052, a low-metallicity (1/41 $`Z_{}`$) galaxy with an unusual IR SED and a strong MIR flux. SBS 0335$``$052 is consistent with a very young starburst of age $`6.5\times 10^6`$ yr and compact size of $`20`$ pc. The model results of young age and compact size are also consistent with the optical properties of the galaxy. We also successfully reproduced the extinction curve of SBS 0335$``$052, which agrees well with the semi-empirical extinction curve deduced from the observational data and our unextincted ‘bare’ continuum SED. We obtain the visual extinction strength $`A_V18`$mag, consistently with Hunt et al.’s estimate of $`A_V=16.3`$mag.

For future observations, we then made a prediction for the SED of another extremely low-metallicity galaxy, I Zw 18. As in the case of the optically visible part of SBS 0335$``$052, I Zw 18 has a low visual extinction, which suggests a small dust opacity. Hence it follows that the FIR luminosity of I Zw 18 is probably low as $`L_{\mathrm{FIR}}10^{7\text{}7.5}L_{}`$, depending on the observational constraint of dust mass. We expect the peak of the IR SED of I Zw 18 to be around $`100\mu `$m, but Galactic cirrus may hamper its detection at these wavelengths. Shorter wavelength observations will be much more effective.

We conclude that both of MIR and FIR surveys are necessary to infer the dust amount of galaxies with a wide range of metallicities. About several hundreds of local low metallicity galaxies may be detected in the such as the SIRTF SWIRE Legacy program and ASTRO-F FIR all-sky survey, and even more detailed observations will be possible by Herschel and other 2–3 m class FIR space telescopes. Such samples will be useful for the systematic studies on the local populations of very young galaxies, as well as for their counterparts at high-$`z`$.

In addition, redshifted peaks of very hot dwarf galaxies will be found in the submillimetre wavelengths ($`50(1+z)\mu `$m, and $`z5`$), and we conclude that the advent of ALMA will enable us to observe such low mass, bona fide primeval galaxies.

## Acknowledgements

First we thank the anonymous referee, whose careful comments improved the quality of this article very much. We are greatly indebted to Veronique Buat, Denis Burgarella, Alessandro Boselli, Marc Sauvage, Takashi Onaka, Akio K. Inoue, Hideyuki Kamaya, Koji S. Kawabata, Misato Fukagawa, Yoko Okada, Masato Onodera, Motohiro Enoki, and Kohji Yoshikawa for fruitful discussions and comments. We also acknowledge Trinh X. Thuan and Marc Sauvage for kindly providing their ISOCAM spectroscopic data of SBS 0335$``$052. TTT and HH have been supported by the Japan Society of the Promotion of Science.


format: text

text:
# 1 Introduction

## 1 Introduction

It is generally agreed that a quarkonium state is not only made up of a quark and an antiquark ($`Q\overline{Q}`$), but contains also a component with glue ($`QA\overline{Q}`$). There are, in addition, other contributions like admixtures of several gluons or light quark-antiquark pairs which are not considered here. We emphasize that we are not talking about hybrids, where the glue component is an essential part of the wave function, but about the admixture of glue in ordinary quarkonia. We are mainly interested in equal time wave functions here, not, say, in light cone wave functions.

The presence of a component with glue in the $`P`$-states has become of great importance due to the work of Bodwin, Braaten, and Lepage . Although this component is suppressed, it allows the $`Q\overline{Q}`$ pair in the $`QA\overline{Q}`$ component to be in a color octet $`S`$-state. The $`S`$-state annihilation of this pair is of the same order as the $`Q\overline{Q}`$ annihilation in the leading $`Q\overline{Q}`$ component of the $`P`$-state, because the latter is also suppressed by the vanishing of the wave function at the origin. Logarithmic divergences cancel in both contributions. Similar considerations hold for $`P`$-wave production.

In spite of the principal as well as phenomenological importance of the $`QA\overline{Q}`$ admixture it is hard to find quantitative statements about this part of the wave function in the literature. We will approach the problem in a systematic way here. Even if one is interested in the equal time wave functions only, it is necessary to start with the general wave functions involving different times. The reason is that quark and gluon propagators also depend on time differences, therefore a treatment involving only equal time wave functions appears not possible within our approach. We leave open the question whether our results could also be obtained within an equal time Hamilton approach. In any case some of the subtleties mentioned below can only be understood if one knows the general wave function with different time arguments.

In sect. 2 we use as starting point two Green functions, one involving $`Q\overline{Q}`$ and the other one $`QA\overline{Q}`$. In the spirit of the usual derivation of the Bethe-Salpeter equation, but with more necessary technicalities, we define several irreducible kernels. These allow a decomposition of the Green functions and a derivation of two coupled Bethe-Salpeter equations for the $`Q\overline{Q}`$ wave function and the $`QA\overline{Q}`$ wave function.

In sect. 3 we proceed by using the framework of non relativistic QCD (NRQCD) , which is appropriate for the case that the velocities of the quarks are small (order $`\alpha _s1`$). We furthermore choose the radiation gauge $`𝐀=0`$ which is known to be the gauge best suited for treating bound state systems. In leading order of perturbation theory the $`Q\overline{Q}`$ kernel becomes static and the $`QA\overline{Q}`$ wave function does not enter into the equation for the $`Q\overline{Q}`$ wave function. The latter can thus be solved in the usual way. After this has been done one can insert the solution for the $`Q\overline{Q}`$ wave function into the equation for the $`QA\overline{Q}`$ wave function and calculate the latter in lowest order.

The result shows some delicate features. The wave function with the Coulomb gluon is discontinuous at $`t=0`$, where $`t`$ is the relative time of the gluon with respect to the time of the $`Q\overline{Q}`$ system. It vanishes for positive $`t`$, while for negative $`t`$ it can be simply understood as the Coulomb field moving along with the quarks. The tree graph contribution to the transverse gluon differs by a factor 1/2 from what one would expect from classical electrodynamics. The origin of this discrepancy can be understood, it is due to the use of time ordered versus retarded propagators in the different approaches. For the transverse gluon there is a non abelian loop contribution, in which the transverse gluon couples to two Coulomb gluons while the latter couple to the quark and antiquark. This loop graph contributes to the same order as the direct coupling of the transverse gluon to quark or antiquark, because the direct coupling is suppressed by a factor of the order of $`|𝐩|/m\alpha _s`$.

In sect. 4 we finally calculate the amplitude for $`Q\overline{Q}`$ at zero distance in the $`QA\overline{Q}`$ component of the $`P`$-state. As stated above, this quantity is of importance for annihilation decays and production of $`P`$-states. The result cannot simply be written in terms of, say, the derivative of the $`Q\overline{Q}`$ radial wave function at the origin, it depends on the gluon momentum as well as on the vector index of the gluon and the magnetic quantum number of the $`P`$-state.

In the conclusions of sect. 5 we briefly discuss how an extension of our perturbative results to a wider range of applicability can be carried out. Applications to $`P`$-state annihilation decays need some further effort and will be presented elsewhere.

## 2 Coupled Bethe-Salpeter equations

We consider a mesonic bound state $`|K>`$ with four momentum $`K^\mu `$. We are interested both in the quark-antiquark amplitude

$$\psi _{Q\overline{Q}}(x_1,x_2;K)=<0|T(\psi (x_1)\overline{\psi }(x_2))|K>=\frac{e^{iKX}}{(2\pi )^4}e^{ipr}\stackrel{~}{\psi }_{Q\overline{Q}}(p;K)d^4p,$$
(2.1)
as well as in the quark-glue-antiquark component

$`\psi _{QA^\mu \overline{Q}}(x_1,x_3,x_2;K)`$ $`=`$ $`<0|T(\psi (x_1)A^\mu (x_3)\overline{\psi }(x_2))|K>`$ (2.2)
$`=`$ $`{\displaystyle \frac{e^{iKX}}{(2\pi )^8}}{\displaystyle e^{i(pr+qx)}\stackrel{~}{\psi }_{QA^\mu \overline{Q}}(p,q;K)d^4pd^4q}.`$
Spin- and color indices are suppressed, explicitly the color singlet states are understood as $`\psi \overline{\psi }(1/\sqrt{N})\overline{\psi }_\alpha \psi _\alpha `$, and $`\psi A^\mu \overline{\psi }\sqrt{2/(N^21)}\overline{\psi }_{\alpha _2}A^{\mu a}(\lambda _{\alpha _2\alpha _1}^a/2)\psi _{\alpha _1}`$, with $`N=3`$ the number of colors. There is no need of specifying the gauge at this stage.

In the two particle $`Q\overline{Q}`$ system, with masses $`m_1,m_2`$, we used the usual CM and relative variables in position and momentum space:

$`X=\eta _1x_1+\eta _2x_2,r=x_1x_2`$ ; $`x_1=X+\eta _2r,x_2=X\eta _1r,`$
$`K=p_1+p_2,p=\eta _2p_1\eta _1p_2`$ ; $`p_1=\eta _1K+p,p_2=\eta _2Kp.`$ (2.3)
Here $`\eta _j=m_j/(m_1+m_2)`$, the total and reduced mass are denoted by $`M=m_1+m_2`$ and $`\mu =m_1m_2/(m_1+m_2)`$. For the three particle system $`QA\overline{Q}`$ it is convenient to use the same CM and relative coordinates $`X,r`$ for $`Q`$ and $`\overline{Q}`$ as above. We denote by $`x`$ the distance of the gluon to the CM of the $`Q\overline{Q}`$, the momentum of the gluon is denoted by $`q`$. The coordinates then read

$`x=x_3\eta _1x_1\eta _2x_2`$ ; $`x_3=X+x,`$
$`K=p_1+p_2+q`$ ; $`p_1=\eta _1(Kq)+p,p_2=\eta _2(Kq)p.`$ (2.4)
The derivation of the bound state equations for $`\stackrel{~}{\psi }_{Q\overline{Q}}(p;K)`$ and $`\stackrel{~}{\psi }_{QA^\mu \overline{Q}}(p,q;K)`$ is a bit technical. One may skip the rest of this section, the resulting lowest order graphs in Fig. 6 and Fig. 7 of the next section are easily understood anyhow. For a systematic and general approach we consider the Green functions

$`G_{Q\overline{Q},Q\overline{Q}}(x_1,x_2,x_1^{},x_2^{})=<0|T([\psi (x_1)\overline{\psi }(x_2)][\psi (x_1^{})\overline{\psi }(x_2^{})])|0>`$
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^{12}}}{\displaystyle e^{i(pr+p^{}r^{}+K(X+X^{}))}\stackrel{~}{G}_{Q\overline{Q},Q\overline{Q}}(p,p^{};K)d^4pd^4p^{}d^4K},`$
and

$`G_{QA^\mu \overline{Q},Q\overline{Q}}(x_1,x_3,x_2,x_1^{},x_2^{})=<0|T([\psi (x_1)A^\mu (x_3)\overline{\psi }(x_2)][\psi (x_1^{})\overline{\psi }(x_2^{})])|0>`$
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^{16}}}{\displaystyle e^{i(pr+p^{}r^{}+qx+K(X+X^{}))}\stackrel{~}{G}_{QA^\mu \overline{Q},Q\overline{Q}}(p,p^{},q;K)d^4pd^4p^{}d^4qd^4K}.`$
Further Green functions, which will be needed below, are defined in an analogous way and immediately recognized from their symbols. Wave functions and Green functions in momentum space are denoted by open circles with the appropriate legs attached to them. The quark gluon vertex function is defined in the usual way as the sum of all one-particle irreducible graphs contributing to the $`QAQ`$-vertex, i.e. the graphs which cannot be separated into two different pieces by cutting only one internal line. The external legs are not included.

With the help of the vertex function the Green function $`\stackrel{~}{G}_{QA^\mu \overline{Q},Q\overline{Q}}`$ can then be decomposed as in Fig. 1, where internal lines, here as well as in the other graphs of this section, denote full propagators.

We define the totally connected Green function $`\stackrel{~}{G}_{QA\overline{Q},Q\overline{Q}}^{conn}`$, denoted with the filled circle in Fig. 1, as the sum of the graphs which cannot be divided into two separate pieces by cutting one internal $`Q`$ (or $`\overline{Q}`$) line plus the left external $`\overline{Q}`$ (or $`Q`$) line.

In order to separate the different components in the bound state we have to introduce four irreducible kernels $`\stackrel{~}{I}_{Q\overline{Q},Q\overline{Q}},\stackrel{~}{I}_{Q\overline{Q},QA^\mu \overline{Q}},\stackrel{~}{I}_{QA^\mu \overline{Q},Q\overline{Q}},\stackrel{~}{I}_{QA^\mu \overline{Q},QA^\nu \overline{Q}}`$. They are defined as the sums of all connected graphs contributing to the corresponding Green functions $`\stackrel{~}{G}_{Q\overline{Q},Q\overline{Q}},\mathrm{},`$ respectively, without external legs, and will be denoted by rectangles with the appropriate lines attached. The irreducible kernels have the following properties:

The graphs cannot be separated into two pieces by cutting two internal lines of the type $`Q\overline{Q}`$, or by cutting three internal lines of the type $`QA\overline{Q}`$.

The kernels $`\stackrel{~}{I}_{Q\overline{Q},QA^\mu \overline{Q}}`$ and $`\stackrel{~}{I}_{QA^\mu \overline{Q},QA^\nu \overline{Q}}`$ have a $`QA\overline{Q}`$ state on the right hand side, and there are three further rules for them:

The graph cannot be separated into two pieces by cutting an internal $`Q`$ (or $`\overline{Q}`$) plus the right external $`\overline{Q}`$ (or $`Q`$) line.

The graph cannot be separated into two pieces by cutting two internal $`QA`$ (or $`A\overline{Q}`$) lines plus the right external $`\overline{Q}`$ (or $`Q`$) line.

The graph cannot be separated into two pieces by cutting two internal $`Q\overline{Q}`$ lines plus the right external $`A`$ line.

The prescriptions b), c), d) only refer to lines on the right hand side. If the corresponding situation appears on the left hand side, the graph is included in the definition of $`\stackrel{~}{I}_{QA^\mu \overline{Q},Q\overline{Q}}`$ or $`\stackrel{~}{I}_{QA^\mu \overline{Q},QA^\nu \overline{Q}}`$.

The above definitions have been chosen in such a way that one is able to extract the desired wave functions while avoiding double counting in the decompositions below. Note that we do not make use of the usual two quark irreducible kernel $`\stackrel{~}{I}`$ consisting of all connected graphs for $`\stackrel{~}{G}_{Q\overline{Q},Q\overline{Q}}`$ (without external legs) which cannot be separated into two pieces by a $`Q\overline{Q}`$ cut. We need a more detailed classification of kernels for our purpose.

Consider now the sum of all graphs which contribute to $`\stackrel{~}{G}_{Q\overline{Q},Q\overline{Q}}`$. They may be collected in different classes using the above definitions. This is shown in Fig. 2.

The corresponding decomposition of $`\stackrel{~}{G}_{QA^\mu \overline{Q},Q\overline{Q}}`$ is shown in Fig. 3.

The correctness of the decompositions in Fig. 2 and Fig. 3 can be demonstrated in any order of perturbation theory by showing that there is a one to one correspondence between graphs. Clearly each graph appearing on the rhs also belongs to the lhs. The proof that each graph on the lhs appears exactly once on the rhs is given in the appendix. There we also show some examples for graphs which are forbidden by our above definitions.

The intermediate state $`|K><K|`$ leads to a pole in the Green functions at $`K^0=\pm (\sqrt{M_K^2+𝐊^2}iϵ)`$ on both sides of Fig. 2 and Fig. 3. Extracting the residue and dropping the common factor $`\stackrel{~}{\psi }_{Q\overline{Q}}^{}`$ we arrive at two coupled Bethe-Salpeter equations for the wave functions $`\stackrel{~}{\psi }_{Q\overline{Q}}`$ and $`\stackrel{~}{\psi }_{QA^\mu \overline{Q}}`$ which are displayed in Fig. 4 and Fig. 5.

## 3 Perturbative solution

Up to now the treatment was general and there was no need to specify the Lagrangian or the gauge in detail. In this section we will specialize to non relativistic QCD (NRQCD up to order $`1/m`$) with the Lagrangian

$$_Q=\psi ^+(iD^0+\frac{\stackrel{}{D}^2}{2m_1}+\frac{g𝐬_1𝐁}{m_1})\psi +\chi ^+(iD^0\frac{\stackrel{}{D}^2}{2m_2}\frac{g𝐬_2𝐁}{m_2})\chi \frac{1}{2}TrF^{\mu \nu }F_{\mu \nu }.$$
(3.1)
Here $`\psi `$ is the Pauli spinor for the quark field (annihilation operator), and $`\chi `$ for the antiquark field (creation operator). Our conventions are $`D^\mu =^\mu +igA^\mu ,F^{\mu \nu }=^\mu A^\nu ^\nu A^\mu +ig[A^\mu ,A^\nu ],E^n=F^{0n},B^n=ϵ^{njk}F^{jk}/2`$. Since we are working in lowest order only, there is no need of introducing matching coefficients in the Lagrangian (3.1). We go to the rest frame of the bound state where $`K^\mu =(|E|,\mathrm{𝟎})`$. Note that the rest masses are dropped in NRQCD, we will also suppress the $`K`$ in our notation from now on. The perturbation theory is performed in radiation gauge $`𝐀=0`$, the gauge best suited for treating bound states. From now on we distinguish between Coulomb gluons $`A^0`$ denoted by dashed lines, and transverse gluons $`𝐀`$ denoted by curly lines.

We first solve the bound state equation in Fig. 4. The leading contribution to the kernel $`\stackrel{~}{I}_{Q\overline{Q},Q\overline{Q}}`$ which appears in the first term on the rhs is given by the exchange of one Coulomb gluon between $`Q`$ and $`\overline{Q}`$ which leads to a static kernel:

$$\stackrel{~}{I}_{Q\overline{Q},Q\overline{Q}}=\frac{i(2\pi )^4g^2(\lambda _{\alpha \beta }^a/2)(\lambda _{\gamma \delta }^a/2)}{|𝐩𝐩^{}|^2}.$$
(3.2)
The corresponding contribution with a transverse gluon contains additional factors of $`𝐩/m`$ at the vertices. It is therefore suppressed if we consider only soft external momenta. Internal momenta are also essentially soft due to the non relativistic wave function given below. The last three contributions in Fig. 4 are also of higher order. In lowest order Fig. 4 therefore gives the familiar non relativistic Bethe Salpeter equation in ladder approximation and radiation gauge. Because the kernel is static, the only dependence on the relative energy variable $`p^0`$ on the rhs is in the propagators. One can therefore integrate over $`p^0`$ on both sides and derive the Schrödinger equation for the momentum space wave function $`\stackrel{~}{\psi }_{Q\overline{Q}}(𝐩)(2\pi )^{5/2}\stackrel{~}{\psi }_{Q\overline{Q}}(𝐩,p^0)𝑑p^0`$ in the usual way. After this has been solved, the dependence on the relative energy can be recovered by going back to the original equation. The well known solution then reads

$$(\stackrel{~}{\psi }_{Q\overline{Q}})_{nlm}(𝐩,p^0)=\frac{i(2\pi )^{3/2}(|E_n|+𝐩_1^2/2m_1+𝐩_2^2/2m_2)\stackrel{~}{\psi }_{nlm}(𝐩)}{(p^0\eta _1|E_n|𝐩_1^2/2m_1+iϵ)(p^0+\eta _2|E_n|+𝐩_2^2/2m_2iϵ)}.$$
(3.3)
The spin wave functions factorize. $`E_n=(C_F\alpha _s)^2\mu /(2n^2)`$, with $`C_F=(N^21)/(2N)=4/3`$ the Casimir operator in the fundamental representation, is the energy eigenvalue, and $`\stackrel{~}{\psi }_{nlm}(𝐩)=\stackrel{~}{\psi }_{nl}(|𝐩|)Y_{lm}(\widehat{𝐩})`$ the normalized Schrödinger wave function in momentum space. The normalization has been chosen such that $`\stackrel{~}{\psi }(𝐩,p^0)𝑑p^0=(2\pi )^{5/2}\stackrel{~}{\psi }(𝐩)`$, which means that $`\psi _{Q\overline{Q}}`$ in (2.1) becomes the correctly normalized Schrödinger wave function if $`x_1^0=x_2^0`$. The Bohr radius is $`a=1/\sqrt{2\mu |E_1|}`$. For the lowest states one has (it is convenient to write them in terms of the energy $`E_n`$ instead of the Bohr radius)

$`\stackrel{~}{\psi }_{10}(|𝐩|)`$ $`=`$ $`{\displaystyle \frac{4\sqrt{2}(2\mu |E_1|)^{5/4}}{\sqrt{\pi }(𝐩^2+2\mu |E_1|)^2}},`$
$`\stackrel{~}{\psi }_{20}(|𝐩|)`$ $`=`$ $`{\displaystyle \frac{8\sqrt{2}(2\mu |E_2|)^{5/4}(|𝐩|^22\mu |E_2|)}{\sqrt{\pi }(|𝐩|^2+2\mu |E_2|)^3}},`$ (3.4)
$`\stackrel{~}{\psi }_{21}(|𝐩|)`$ $`=`$ $`{\displaystyle \frac{i16\sqrt{2}(2\mu |E_2|)^{7/4}|𝐩|}{\sqrt{3\pi }(|𝐩|^2+2\mu |E_2|)^3}}.`$
Having solved the bound state equation in Fig. 4 in lowest order, one can insert the solution, together with the leading approximations for the propagators, the vertex functions, and the kernels, into the rhs of Fig. 5, dropping the higher order contributions of the second line. The power counting in the coupling constant $`g`$ requires some care, because there are not only the explicit factors $`g`$ in the vertices, but also factors $`p/mg^2`$ in the coupling of transverse gluons to quarks, as well as momentum factors in the three gluon vertex. Taking this into account, the leading contribution to $`\stackrel{~}{\psi }_{QA^0\overline{Q}}`$ is then given by Fig. 6, the leading contribution to $`\stackrel{~}{\psi }_{Q𝐀\overline{Q}}`$ by Fig. 7.

In Fig. 7 c) the kernel $`\stackrel{~}{I}_{Q𝐀\overline{Q},Q\overline{Q}}`$ has been replaced by its lowest order contribution, where the transverse gluon couples to two Coulomb gluons, while the latter couple to the quark and antiquark, respectively. For the reasons mentioned above, this graph, although of higher order in the coupling constant, contributes to the same order as the tree graphs in Fig. 7 b), because the direct coupling of the transverse gluon to the quark or antiquark in 7 b) is suppressed by a factor $`p/mg^2`$. For the Coulomb gluon the corresponding graph does not exist, because there is no coupling of three Coulomb gluons. There is no non abelian contribution in leading order in Fig. 6.

We can now insert the $`Q\overline{Q}`$ wave function (3.3) into the rhs of Fig. 6 and Fig. 7, and calculate the wave functions on the lhs. We start with the tree graphs Fig. 6 a) and Fig. 7 b). They can be simply written down because there is no integration over internal momenta. We are interested in the equal time wave function, therefore we first integrate over $`p^0`$ (or, more simply, over $`p_2^0`$ in $`a_1),b_1)`$ and over $`p_1^0`$ in $`a_2),b_2))`$ with the residue theorem. There are three denominators, one arising from the external $`Q`$ or $`\overline{Q}`$ propagator and two from the $`Q\overline{Q}`$ wave function. For $`a_1),b_1)`$ e.g. one gets

$`{\displaystyle \frac{dp_2^0}{(p_2^0+q^0+|E_n|+𝐩_1^2/2m_1iϵ)(p_2^0+|E_n|+𝐩_{}^{}{}_{1}{}^{2}/2m_1iϵ)(p_2^0𝐩_2^2/2m_2+iϵ)}}`$ (3.5)
$`=`$ $`{\displaystyle \frac{2\pi i}{(q^0+|E_n|+𝐩_1^2/2m_1+𝐩_2^2/2m_2iϵ)(|E_n|+𝐩_{}^{}{}_{1}{}^{2}/2m_1+𝐩_2^2/2m_2)}}.`$
Next we have to perform the integration over $`q^0`$. For the Coulomb gluon the propagator is independent of $`q^0`$, one thus cannot simply perform the $`q^0`$-integration, because the denominator in (3.5) contains only one power of $`q^0`$. We therefore choose an infinitesimal relative time $`t`$ of the gluon with respect to the quark and antiquark. The Fourier transform then becomes

$$\frac{e^{iq^0t}dq^0}{(q^0+|E_n|+𝐩_1^2/2m_1+𝐩_2^2/2m_2iϵ)}=2\pi i\mathrm{\Theta }(t)\text{ for }t\text{ infinitesimal}.$$
(3.6)
We will comment this result below. Finally we add the contribution of Fig. 6 $`a_2)`$ and divide by $`(2\pi )^5`$ which gives the correctly normalized equal time wave function $`\stackrel{~}{\psi }(𝐩,𝐪)`$ in momentum space. This results in

$$\stackrel{~}{\psi }_{QA^0\overline{Q}}(𝐩,𝐪)=\frac{g\sqrt{C_F}\mathrm{\Theta }(t)}{(2\pi )^{3/2}|𝐪|^2}\left(\stackrel{~}{\psi }(𝐩_2)\stackrel{~}{\psi }(𝐩_1)\right).$$
(3.7)
For a transverse gluon the propagator provides enough powers of $`q^0`$ in the denominator to make the equal time wave function well defined. The result from the graphs Fig. 7 $`b_1),b_2)`$ becomes

$`(\stackrel{~}{\psi }_{Q𝐀\overline{Q}})^{(b)}(𝐩,𝐪)`$ $`=`$ $`{\displaystyle \frac{g\sqrt{C_F}}{(2\pi )^{3/2}2|𝐪|\left(|𝐪|+|E_n|+𝐩_1^2/2m_1+𝐩_2^2/2m_2\right)}}\times `$ (3.8)
$`\left({\displaystyle \frac{[(𝐩_1)_T+i𝐪\times 𝐬_1]\stackrel{~}{\psi }(𝐩_2)}{m_1}}{\displaystyle \frac{[(𝐩_2)_T+i𝐪\times 𝐬_2]\stackrel{~}{\psi }(𝐩_1)}{m_2}}\right).`$
Here $`𝐩_T𝐩(\mathrm{𝐩𝐪})𝐪/𝐪^2`$ denotes the transverse component of $`𝐩`$.

We finally come to the non abelian graph Fig. 7 c). This involves an integration over the internal four momentum $`p^{}`$. Performing first the integrations over $`p^0,p^0`$, and $`q^0`$ gives the equal time wave function

$`(\stackrel{~}{\psi }_{Q𝐀\overline{Q}})^{(c)}(𝐩,𝐪)`$ $`=`$ $`{\displaystyle \frac{g^3N\sqrt{C_F}}{(2\pi )^{9/2}2|𝐪|\left(|𝐪|+|E_n|+𝐩_1^2/2m_1+𝐩_2^2/2m_2\right)}}\times `$ (3.9)
$`{\displaystyle \frac{(𝐩𝐩^{})_T\stackrel{~}{\psi }(𝐩^{})d^3p^{}}{|𝐩𝐩^{}\eta _1𝐪|^2|𝐩𝐩^{}+\eta _2𝐪|^2}}.`$
We come to the physical interpretation of the above results. First we have to recall that there are three different scales in NRQCD: Hard ($`m)`$ soft $`(\alpha _sm)`$ ultrasoft $`(\alpha _s^2m)`$. Quark masses are hard, quark momenta soft, and quark kinetic energies as well as binding energies are ultrasoft. Since we used NRQCD, the momenta $`p_1,p_2,q`$ in the $`QA\overline{Q}`$ wave functions (3.7) - (3.9) must also be soft. In fact these expressions are dominated by soft momenta.

The graphs in Fig. 6 describe the Coulomb field moving along with the quark and antiquark. It is easy to perform the Fourier transform of (3.7) with respect to $`𝐩`$ and $`𝐪`$. It becomes

$$\psi _{QA^0\overline{Q}}=\frac{g\sqrt{C_F}\mathrm{\Theta }(t)}{4\pi }\psi (𝐫)\left(\frac{1}{|𝐱_3𝐱_1|}\frac{1}{|𝐱_3𝐱_2|}\right).$$
(3.10)
One recognizes the product of the $`Q\overline{Q}`$-wave function with the Coulomb potential generated by the $`Q`$ and the $`\overline{Q}`$ source. (In QED one has, of course, to replace $`g\sqrt{C_F}`$ by $`e`$.) The appearance of the step function $`\mathrm{\Theta }(t)`$ is due to the time ordering prescription in the wave function (2.2). For $`t>0`$ the gluon field $`A^0(x_3)`$ has to stand to the left. One can express it in terms of the quark charge densities, and the latter annihilate the vacuum in lowest order.

The graphs in Fig. 7 b) describe the fields which arise because the quarks are moving. For large $`|𝐪|`$ the spin independent contributions of the quark and antiquark are those of the Coulomb contributions (3.7) for $`t<0`$, multiplied with $`(1/2)(𝐩_j)_T/m_j`$. With the exception of the factor $`1/2`$ in front, this is exactly what one would expect from classical electrodynamics. The appearance of this factor can be understood as follows. In classical radiation problems one uses the retarded Green function. In field theory, like in our definition of the wave function (2.2), on uses, however, the time ordered product. If we had used the retarded product instead, the denominator of the transverse gluon propagator would have been changed from $`q^2+iϵ`$ to $`(q^0+iϵ)^2|𝐪|^2`$, such that both poles in the complex $`q^0`$-plane would lie in the lower half plane. In this case one would simply close the contour around the remaining $`q^0`$-pole from the quark propagator in the upper half plane and would end up with the replacement

$$\frac{1}{2|𝐪|\left(|𝐪|+|E_n|+𝐩_1^2/2m_1+𝐩_2^2/2m_2\right)}\frac{1}{|𝐪|^2\left(|E_n|+𝐩_1^2/2m_1+𝐩_2^2/2m_2iϵ\right)^2},$$
(3.11)
removing the factor $`1/2`$ for large $`|𝐪|`$ in (3.8). The use of the advanced propagator would, on the other hand, give a vanishing result. This shows that one has to define the equal time wave function very carefully.

In the following we will stick to our definition with the time ordered product. If one defines the equal time wave function for $`A^0`$ as the average of the expressions for positive and negative infinitesimal $`t`$, the replacement $`\mathrm{\Theta }(t)1/2`$ leads to the correct relation between $`A^0`$ and $`𝐀`$ for large $`|𝐪|`$. For smaller $`|𝐪|`$ the result for $`𝐀`$ is modified, the denominator is less singular for $`|𝐪|0`$. This is, of course, due to the fact that the quarks are not free, there is an additional long range radiation field. Fig. 7 c) finally arises from the glue moving along with the quarks, it has no simple classical analog.

The expressions for the transverse gluon can be simplified. From the denominators in (3.8) and (3.9) it is obvious that these expressions are dominated by ultrasoft $`q`$. If we thus restrict to ultrasoft $`q`$ and remember that $`p_1`$ and $`p_2`$ are soft, we may approximate (3.8) by using $`𝐩_1=𝐩\eta _1𝐪𝐩,𝐩_2=𝐩+\eta _2𝐪𝐩`$, and $`𝐩_1^2/2m_1+𝐩_2^2/2m_2=𝐩^2/2\mu +𝐪^2/2M𝐩^2/2\mu `$, and get

$$(\stackrel{~}{\psi }_{Q𝐀\overline{Q}})^{(b)}(𝐩,𝐪)=\frac{g\sqrt{C_F}𝐩_T\stackrel{~}{\psi }(𝐩)}{(2\pi )^{3/2}2\mu |𝐪|\left(|𝐪|+|E_n|+𝐩^\mathrm{𝟐}/2\mu \right)}.$$
(3.12)
We dropped the spin dependent term because it is suppressed for ultrasoft $`q`$.

We next discuss the non abelian contribution (3.9). For ultrasoft $`q`$ one may again simplify the kinetic energies in the denominator, furthermore one can drop the $`𝐪`$ in the integrand. The integral can then be easily performed for a given $`\stackrel{~}{\psi }(𝐩^{})`$, e.g. by introducing a Feynman parameter. We give the results for the three lowest states.

$`(\stackrel{~}{\psi }_{Q𝐀\overline{Q}})_{100}^{(c)}(𝐩,𝐪)`$ $`=`$ $`{\displaystyle \frac{g^3N\sqrt{C_F}(2\mu |E_1|)^{3/4}𝐩_T}{(2\pi )^3|𝐪|\left(|𝐪|+|E_1|+𝐩^2/2\mu \right)(𝐩^2+2\mu |E_1|)^2}}Y_{00},`$ (3.13)
$`(\stackrel{~}{\psi }_{Q𝐀\overline{Q}})_{200}^{(c)}(𝐩,𝐪)`$ $`=`$ $`{\displaystyle \frac{g^3N\sqrt{C_F}(2\mu |E_2|)^{3/4}(𝐩^232\mu |E_2|)𝐩_T}{(2\pi )^3|𝐪|\left(|𝐪|+|E_2|+𝐩^2/2\mu \right)(𝐩^2+2\mu |E_2|)^3}}Y_{00}.`$ (3.14)
The result for the $`P`$-state consists of two parts:

$`(\stackrel{~}{\psi }_{Q𝐀\overline{Q}})_{21m}^{(c)}(𝐩,𝐪)`$ $`=`$ $`{\displaystyle \frac{i4g^3N\sqrt{C_F}(2\mu |E_2|)^{5/4}}{(2\pi )^3\sqrt{3}|𝐪|\left(|𝐪|+|E_2|+𝐩^2/2\mu \right)}}\times `$ (3.15)
$`\left({\displaystyle \frac{𝐩_T|𝐩|Y_{1m}(\widehat{𝐩})}{(𝐩^2+2\mu |E_2|)^3}}{\displaystyle \frac{1}{4}}{\displaystyle \frac{(𝐞_k)_TY_{1m}(𝐞_k)}{(𝐩^2+2\mu |E_2|)^2}}\right),`$
where $`𝐞_k`$ denote the three Cartesian unit vectors and $`(𝐞_k)_T`$ their transverse parts.

For the ground state, as well as for the first term of the $`P`$-state, we may express the result directly in terms of the wave function and compare it with the result of Fig. 7 b) in (3.12). The non abelian contribution (3.13) for the ground state is $`N/(2C_F)=9/8`$ times the result in (3.12), while for the $`P`$-state the first contribution in (3.15) is $`N/C_F=9/4`$ times the corresponding result in (3.12). Thus numerically the non abelian loop graph is more important than the tree graphs.

In momentum space the leading order of the transverse electrical field $`𝐄_T=^0𝐀`$ and the magnetic field $`𝐁=\times 𝐀`$ are obtained by replacing $`i𝐪`$ and $`^0iq_0i|𝐪|`$ and applying these factors in the expressions for $`𝐀`$. The last replacement is due to the fact that the $`q^0`$-integration was performed by closing the contour around the pole at $`q^0=|𝐪|iϵ`$. The longitudinal electric field $`𝐄_L=A^0`$ can be obtained from (3.7) by multiplying with $`i𝐪`$, or directly in position space from (3.10).

## 4 The $`Q𝐀\overline{Q}`$ wave function of the $`P`$-state at $`r=0`$

A quantity of particular interest is the amplitude for the special case of quarks with identical flavor at zero distance in the $`Q𝐀\overline{Q}`$ component of the $`P`$-state. This component, although suppressed, contributes to annihilation decays to the same order as the $`Q\overline{Q}`$ annihilation in the leading $`Q\overline{Q}`$ component. The reason is that this annihilation is also suppressed, because the wave function at the origin vanishes for $`P`$-states. On the other hand the $`Q\overline{Q}`$ in the $`QA\overline{Q}`$ wave function can be in an $`S`$-state. In one loop order one finds that logarithmic divergences are canceled by adding up both contributions. This has been discussed in detail by Bodwin, Braaten and Lepage . Let us therefore consider the wave function for quarks at zero distance

$$\stackrel{~}{\psi }_{Q𝐀\overline{Q}}(𝐪)_m(2\pi )^{3/2}\stackrel{~}{\psi }_{Q𝐀\overline{Q}}(𝐩,𝐪)_{21m}d^3p.$$
(4.1)
The calculation is simple if we restrict to the dominating region where $`q`$ is ultrasoft. From (3.12), (3.15) we find

$`\stackrel{~}{\psi }_{Q𝐀\overline{Q}}(𝐪)_m`$ $`=`$ $`{\displaystyle \frac{ig\sqrt{C_F}(m|E_2|)^{5/4}(𝐞_k)_TY_{1m}(𝐞_k)}{3\sqrt{6}\pi ^{3/2}|𝐪|\left(\sqrt{m(|𝐪|+|E_2|)}+\sqrt{m|E_2|}\right)^3}}\times `$ (4.2)
$`\left(3\sqrt{m(|𝐪|+|E_2|)}+\sqrt{m|E_2|}{\displaystyle \frac{2N}{C_F}}\sqrt{m|E_2|}\right).`$
The first two terms in the second line originate from the graphs in Fig. 7 b) and the resulting contribution (3.12), the third one is a sum of $`\left(N/C_F\right)\left\{3\sqrt{m(|𝐪|+|E_2|)}+\sqrt{m|E_2|}\right\}`$ from the first part of (3.15) and of $`\left(3N/C_F\right)\left\{\sqrt{m(|𝐪|+|E_2|)}+\sqrt{m|E_2|}\right\}`$ from the second part.

Our result (4.2) shows a remarkable feature. Recalling that $`2N/C_F=9/2`$, one finds that the non abelian term dominates the abelian one for small $`|𝐪|`$ and makes the second line of (4.2) negative. For $`|𝐪|=(13/36)|E_2|`$ the bracket vanishes, for larger $`|𝐪|`$ the positive abelian term dominates.

The wave function (4.2) determines the amplitude for annihilation of the two quarks in the color octet state into gluons which subsequently hadronize. As one might have expected, (4.2) depends on the gluon momentum $`𝐪`$, as well as on the relative orientation between the gluon and the angular momentum of the $`P`$-state. It is, however, independent of the quark spins. The $`Q\overline{Q}`$ wave function has entered in a complicated way. It is not possible to express the result simply by the derivative $`R^{}(0)`$ of the radial wave function at the origin.

## 5 Conclusions

In the foregoing we derived some results which may appear rather unexpected at first sight. One of them is the discontinuity of the wave function in the relative time of the Coulomb gluon with respect to the quarks. The other one is the relation between transverse and Coulomb gluon. Both results could, however, be understood, the point is the use of different Green functions (retarded versus time ordered) in classical radiation theory and in field theory.

One may wonder whether some of these features could be artifacts of the gauge or of the perturbative treatment. We are confident that this is not the case. One could, of course, choose manifestly gauge invariant wave functions by connecting $`Q\overline{Q}`$ or $`QA^\mu \overline{Q}`$, respectively, with path ordered strings. The $`Q\overline{Q}`$ wave function would now contain already some glue in it due to the presence of the string. This glue along the straight string is, however, not an additional degree of freedom. Everything is fixed by the positions of the quark and the antiquark. The essential feature of a $`QA^\mu \overline{Q}`$ wave function is not the presence of some glue, but the presence of an additional degree of freedom carried by the gluon.

The strings would not contribute in lowest order perturbation theory. A consideration of higher order perturbative contributions within the approach presented here, whether with or without strings, is certainly possible. It will not alter the essentials of our results. Our method could also be extended to obtain wave functions for multi gluon admixtures which arise in higher orders.

A systematic treatment of non perturbative effects is beyond the present investigation. We emphasize, however, that our whole approach can be immediately generalized to the case where we replace the kernel (3.2) by some other static kernel, say the Fourier transform of a phenomenological static potential involving confinement. This will, in turn, modify the $`Q\overline{Q}`$ Schrödinger wave function and the subsequent results. For phenomenological applications such an approach appears quite legitimate, it has been familiar and successful in other fields of quarkonium physics.

The application of our results to production or decays of heavy quarkonia goes beyond the present investigation. In the framework of NRQCD, hard annihilation processes are described in lowest order by four fermion terms of dimension six. They were not included in our Lagrangian because they are of order $`1/m^2`$ and thus not important for the calculation of the wave function. While the hard subprocess of $`S`$-wave annihilation is simple and well understood, the complicated dependence of our wave function on the gluon momentum and the magnetic quantum number of the $`P`$-state require some additional work on this topic. While we know that logarithmic divergences which show up in higher orders cancel exactly against corresponding divergences in the $`P`$-wave annihilation of the $`Q\overline{Q}`$ component, there is no information about the finite part up to now. This information can be obtained from our explicit wave function. Due to the compensation between the two contributions in (4.2) and the change in sign between smaller and larger $`|𝐪|`$, suppressions compared to naive expectations may be expected. The results will be presented elsewhere.

Acknowledgement: I thank D. Melikhov and N. Brambilla for reading the manuscript and for valuable suggestions.

## Appendix A Decomposition of Green functions

In this appendix we prove the correctness of the decompositions in Fig. 2 (the proof for Fig. 3 is analogous) by showing that every graph on the lhs appears exactly once on the rhs. Consider an arbitrary graph contributing to the lhs.

1) If the graph is disconnected it uniquely belongs to the first term on the rhs of the first line of Fig. 2.

Consider a connected graph next.

2) If the graph cannot be separated by a $`Q\overline{Q}`$ or a $`QA\overline{Q}`$ internal cut, it belongs to $`I_{Q\overline{Q},Q\overline{Q}}`$ plus external lines. It then belongs to the second term on the rhs, where the disconnected part of the Green function appears.

Let’s next take the case where the graph can be cut by an internal $`Q\overline{Q}`$ or $`QA\overline{Q}`$ cut. The propagators are always taken as part of the right hand graph. We first construct a unique preliminary cut $`C_p`$ as follows: Follow the quark line from the left until one comes to an internal propagator where the graph can be cut in the above way. If there are several possibilities to do so, let $`C_p`$ be the cut where the antiquark line is cut as far to the left as possible.

3) If $`C_p`$ is a $`Q\overline{Q}`$ cut we take it as the final cut $`C`$. In this case the graph belongs to the second term on the rhs of Fig. 2 with a connected contribution to the Green function.

If $`C_p`$ is a $`QA\overline{Q}`$ cut there are two alternatives:

4) It may be possible that one can transform the cut by keeping it’s position at the $`Q`$ (or $`\overline{Q}`$) line, while moving the position of the cut at the $`\overline{Q}`$ (or $`Q`$) to the right, such that one obtains a $`Q\overline{Q}`$ cut. This happens if one can shift the cut to the right by moving it over a vertex part, such that it does no longer cut the gluon line. If this procedure is possible at all, it is unique. In this case we take this $`Q\overline{Q}`$ cut as the final $`C`$. Note that the graph left to the cut is a contribution to $`\stackrel{~}{I}_{Q\overline{Q},Q\overline{Q}}`$; it is not allowed to cut the graph by the original $`C_p`$, because it would cut an external $`Q`$ or $`\overline{Q}`$ line now. Again we get a contribution to the second term on the rhs of fig. 2.

5) If the above manipulation is not possible, we take the original $`QA\overline{Q}`$ cut $`C_p`$ as the final $`C`$. In this case the graph on the lhs is a part of $`\stackrel{~}{I}_{Q\overline{Q},QA\overline{Q}}`$, while the one on the rhs is a part of the totally connected Green function $`\stackrel{~}{G}_{QA\overline{Q},Q\overline{Q}}^{conn}`$. Contributions to the first two terms on the rhs of Fig. 1 cannot appear by construction because our prescription would have moved the corresponding cut to the right.

We have thus shown that every graph on the lhs also appears on the rhs. In order to proof that it belongs exactly to one term on the rhs we will next show that it is not possible to move the cut $`C`$ to the right, producing thereby another admissible cut $`C_r`$, which would lead to a different decomposition. (By construction it is certainly not possible to move the cut to the left).

Consider first the case that $`C`$ is a $`Q\overline{Q}`$ cut. If we move it to the right, crossing a vertex both at the $`Q`$ line as well as at the $`\overline{Q}`$ line (irrespective whether we arrive at an admissible cut), the remaining graph on the left could be cut by the original $`C`$ and thus would not belong to an irreducible kernel. If we move the cut only over one vertex, such that it becomes a $`QA\overline{Q}`$ cut, the graph on the lhs would be inacceptable because of rule b) in the definition of the kernels in sect. 2.

Consider next the case that $`C`$ is a $`QA\overline{Q}`$ cut. If we move it to the right by crossing a vertex both at the $`Q`$ line as well as at the $`\overline{Q}`$ line (irrespective whether we arrive at an admissable cut), the remaining graph on the left could be cut by the original $`C`$ and thus would not belong to an irreducible kernel. Rule d) of sect. 2 applies in the case that the gluon line which was cut by $`C`$ is also cut by $`C_r`$ and thus would become an external line.

If we move the cut only along the $`Q`$ (or the $`\overline{Q}`$) line while keeping the position at the $`\overline{Q}`$ (or the $`Q`$) line we may either obtain a $`Q\overline{Q}`$ cut or a $`QA\overline{Q}`$ cut (or an inadmissible cut). In the first case the kernel to the left is forbidden by rule b), in the second case by c).

The somewhat complicated definition of the kernels $`\stackrel{~}{I}_{Q\overline{Q},QA\overline{Q}}`$, $`\stackrel{~}{I}_{QA\overline{Q},QA\overline{Q}}`$ with respect to external lines at the rhs was necessary to guarantee the uniqueness of the decomposition. In Fig. 8 we show some graphs which do not belong to $`\stackrel{~}{I}_{Q\overline{Q},Q\overline{Q}}`$ or to $`\stackrel{~}{I}_{Q\overline{Q},QA\overline{Q}}`$, together with the cuts and the rule of sect. 2 which prevents this.