text:
# 1 Introduction

## 1 Introduction

In the study of the sinh-Gordon model restricted to a half-line by integrable boundary conditions, an important issue is to determine how the reflection factor depends upon the two parameters introduced at the boundary. Some progress was recently made in tackling this question . Namely, the idea of  was to compute the bound-state spectrum of the model in two different ways, firstly by using a boundary bootstrap principle, and secondly by quantizing the classical boundary breather states using a WKB approach . Each method provides an independent description of the energy spectrum of the boundary states, and comparing the two provides information relating the boundary parameters and the bulk coupling to the reflection factor. Clearly, one expects the most general situation where the two boundary parameters are independent to be technically more involved. It is therefore natural - as was done in  \- to first implement the programme outlined above in the special case where the bulk symmetry $`\varphi \varphi `$ is preserved at the boundary, by requiring the two parameters to be equal. A further step was taken in , where a perturbative calculation to lowest order in the bulk coupling and to first order in the difference of the two boundary parameters was used to make an informed guess as to the general dependence of the reflection factor on those boundary parameters. The results of the present paper, which extends the analysis in  to the case where the two boundary parameters are different, will underpin that guess. Finally, a version of a weak-strong coupling duality transformation appropriate to the model with boundary conditions was also proposed in  and the conclusions drawn from our paper give further support to that conjecture.

Following a brief introduction to the sinh-Gordon model with integrable boundary conditions in section 1, we recall the classical breather solutions in section 2 and describe the properties and facts we need for the semi-classical quantization. In section 3, we review the boundary state bootstrap and use it to derive a formula for the boundary state energy spectrum. The adapted Dashen-Hasslacher-Neveu method is described in section 4 and used to provide an alternative calculation of the boundary state spectrum in terms of the boundary parameters. Our results and some additional remarks are summarized in section 5.

## 2 The sinh-Gordon model on the half-line

The sinh-Gordon model describes a single real scalar field $`\varphi `$ in 1+1 dimensions with exponential self-interaction. The field equation is

$$_t^2\varphi _x^2\varphi +\frac{\sqrt{8}m^2}{\beta }\mathrm{sinh}(\sqrt{2}\beta \varphi )=0,$$
(2.1)
where $`m`$ and $`\beta `$ are parameters and we have used normalizations customary in affine Toda field theories of which the sinh-Gordon model is the simplest . The dimensional mass parameter $`m`$ will be set to unity for convenience.

In contrast to the sine-Gordon model, with its soliton and breather solutions, the sinh-Gordon model is at first sight relatively uninteresting. There is a constant vacuum solution $`\varphi =0`$ and, in the quantum theory, the small oscillations around this vacuum correspond to the sinh-Gordon particle. In the bulk, the spectrum of the model consists of a single species of scalar particle interacting with itself. Nevertheless, however simple the model may appear, it is not easy to analyze it directly , and much of what is known about it has been deduced from features of the lightest breather in the sine-Gordon theory.

The sinh-Gordon model is integrable which implies in particular that there are infinitely many mutually commuting, independent conserved charges $`Q_{\pm s}`$, where $`s`$ is any odd integer, and the S-matrix describing the scattering of two sinh-Gordon particles with relative rapidity $`\mathrm{\Theta }`$ is conjectured to be given by ,

$$S(\mathrm{\Theta })=\frac{1}{(B)_\mathrm{\Theta }(2B)_\mathrm{\Theta }}.$$
(2.2)
In (2.2) we have used the convenient block notation

$$\left(x\right)_\mathrm{\Theta }=\frac{\mathrm{sinh}\left(\frac{\mathrm{\Theta }}{2}+\frac{i\pi x}{4}\right)}{\mathrm{sinh}\left(\frac{\mathrm{\Theta }}{2}\frac{i\pi x}{4}\right)},$$
(2.3)
and the coupling constant $`B`$ is related to the bare coupling constant $`\beta `$ by $`B=2\beta ^2/(4\pi +\beta ^2)`$. For compactness, we will generally omit the subscript $`\mathrm{\Theta }`$ from the block notation since it is generally clear from the context what is intended.

The sinh-Gordon model can be restricted to the left half-line $`\mathrm{}x0`$ without losing integrability by imposing the boundary condition

$$_x\varphi |_0=\frac{\sqrt{2}m}{\beta }\left(\epsilon _0e^{\frac{\beta }{\sqrt{2}}\varphi (0,t)}\epsilon _1e^{\frac{\beta }{\sqrt{2}}\varphi (0,t)}\right),$$
(2.4)
where $`\epsilon _0`$ and $`\epsilon _1`$ are two additional parameters . Again, $`m`$ sets the scale but is taken to be unity for convenience in what follows. This set of boundary conditions generally breaks the reflection symmetry $`\varphi \varphi `$ of the model although the symmetry is preserved when $`\epsilon _0=\epsilon _1\epsilon `$.

Assuming factorization, the description of the sinh-Gordon particles on the half line does not only require the two-particle scattering amplitude (2.2), but also the amplitude for the reflection of a single particle from the boundary. This reflection amplitude was deduced from the lowest breather reflection amplitude in the sine-Gordon model by analytic continuation in the coupling constant (i.e., the continuation $`\lambda 2/B`$ in the notation of ). Using the breather reflection amplitudes calculated by Ghoshal , the analytic continuation leads to cccIn Ghoshal’s notation $`E=B\eta /\pi ,F=iB\vartheta /\pi `$.

$$K_q(\theta ,\epsilon _0,\epsilon _1,\beta )=\frac{(1)(2B/2)(1+B/2)}{(1E(\epsilon _0,\epsilon _1,\beta ))(1+E(\epsilon _0,\epsilon _1,\beta ))(1F(\epsilon _0,\epsilon _1,\beta ))(1+F(\epsilon _0,\epsilon _1,\beta ))},$$
(2.5)
where we are again using the block notation from (2.3) but, here, $`\theta `$ represents the rapidity of a single particle. When the bulk reflection symmetry is preserved one of the two parameters $`E`$ or $`F`$ vanishes. Without loss of generality, we shall take the vanishing parameter in the symmetric case to be $`F`$. All reflection factors satisfy the crossing-unitarity relation which, in the case of scalar reflection factors, reads,

$$K_q\left(\theta +\frac{i\pi }{2}\right)K_q\left(\theta \frac{i\pi }{2}\right)S(2\theta )=1.$$
(2.6)
Actually, (2.5) is the simplest solution to the crossing-unitarity relation using the S-matrix (2.2), taking into account the independently calculated classical limit of the reflection factors. In  the classical reflection factor was found to be given by the formula

$$K_c(\theta ,\epsilon _0,\epsilon _1,\beta )=\frac{(1)^2}{(1a_0a_1)(1+a_0+a_1)(1a_0+a_1)(1+a_0a_1)},$$
(2.7)
in which it was convenient to use an alternative expression for the boundary parameters, viz.,

$$ϵ_0=\mathrm{cos}\pi a_0,ϵ_1=\mathrm{cos}\pi a_1.$$
(2.8)
The formula is well-defined provided we select the ranges $`0a_i1`$ for $`i=0,1`$, to ensure a one to one correspondence between the alternative parameters. Clearly, (2.5) has the correct limit provided

$$Ea_0+a_1,Fa_0a_1.$$
(2.9)
Recently, on the basis of a perturbative calculation, the relationship between the various parameters was conjectured to be ,

$$E=(a_0+a_1)(1B/2),F=(a_0a_1)(1B/2).$$
(2.10)
The principal purpose of this paper is to provide further evidence for these formulae.

It was noted in  that contrary to the situation on the whole line, the sinh-Gordon equation restricted to a half-line by integrable boundary conditions has non-singular, finite energy, breather solutions. These solutions were described in some detail in , particularly for the special case in which $`ϵ_0=ϵ_1`$. Here, we shall concentrate on the solutions containing two independent boundary parameters. The existence of these special solutions and their associated states makes it clear that the sinh-Gordon model has a rather more interesting structure than one would be led to believe on the basis of bulk calculations.

## 3 Boundary Breathers

The boundary breathers may be conveniently described following Hirota’s prescription . Generally, provided we set

$$\varphi =\frac{\sqrt{2}}{\beta }\mathrm{ln}\frac{\tau _+}{\tau _{}},$$
(3.1)
and choose the two $`\tau `$-functions as follows:

$$\tau _\pm =1\pm (E_1+E_2+E_3)+(A_{12}E_1E_2+A_{13}E_1E_3+A_{13}E_1E_3)\pm A_{12}A_{13}A_{23}E_1E_2E_3,$$
(3.2)
with

$$E_p=e^{a_px+b_pt+c_p},a_p=2\mathrm{cosh}\rho _p,b_p=2\mathrm{sinh}\rho _p,A_{pq}=\mathrm{tanh}^2\left(\frac{\rho _p\rho _q}{2}\right),$$
(3.3)
then we have just enough freedom to accommodate the general boundary conditions (2.4). In effect, with general boundary conditions, the solutions have the flavour of a ‘breather’ superimposed on a stationary ‘soliton’, borrowing the language of the sine-Gordon model. In detail, it is enough to choose $`E_1`$ and $`E_2`$ to be complex conjugate partners and periodic in $`t`$, and to take $`E_3`$ to be real and independent of $`t`$; thus, setting $`\rho _1=\rho _2=i\rho `$ and $`\rho _3=0`$, one has:

$$E_1=e^{2x\mathrm{cos}\rho +2it\mathrm{sin}\rho +c}=E_2^{},E_3=e^{2x+d};A_{12}=\mathrm{tan}^2\rho ,A_{13}=A_{23}=\mathrm{tan}^2(\rho /2).$$
(3.4)
The period of the breather is $`\pi /\mathrm{sin}\rho `$. To match the boundary conditions, $`c`$ and $`d`$ need to be determined. It is convenient (actually just a change of origin in $`t`$) to take $`c`$ to be real. Then, the boundary conditions (2.4) are satisfied provided

$$e^c=\frac{s}{\mathrm{tan}\rho },e^d=\frac{r}{\mathrm{tan}^2(\rho /2)}$$
(3.5)
where

$`r={\displaystyle \frac{\mathrm{sin}\frac{\pi a_0}{2}\mathrm{sin}\frac{\pi a_1}{2}}{\mathrm{sin}\frac{\pi a_0}{2}+\mathrm{sin}\frac{\pi a_1}{2}}}`$
(3.6)
$`s^2={\displaystyle \frac{1+\mathrm{cos}\rho }{1\mathrm{cos}\rho }}{\displaystyle \frac{\mathrm{cos}\frac{\pi (a_0+a_1)}{2}+\mathrm{cos}\rho }{\mathrm{cos}\frac{\pi (a_0+a_1)}{2}\mathrm{cos}\rho }}{\displaystyle \frac{\mathrm{cos}\frac{\pi (a_0a_1)}{2}\mathrm{cos}\rho }{\mathrm{cos}\frac{\pi (a_0a_1)}{2}+\mathrm{cos}\rho }}.`$
In (3) we have made use of the alternative parametrisation (2.8). There are other ways of writing (3) which are useful when it comes to evaluating certain integrals. For example, setting $`q=\mathrm{tan}(\rho /2)`$ and $`q_\pm =\mathrm{tan}(\pi a_\pm /2)`$, where $`a_\pm =(a_0\pm a_1)/2`$, we have

$$s^2=\frac{1}{q^2}\frac{1q^2q_+^2}{1q^2q_+^2}\frac{1q^2q_{}^2}{1q^2q_{}^2},r=\frac{q_{}}{q_+}.$$
(3.7)
For symmetrical boundary conditions with $`\epsilon _0=\epsilon _1=\epsilon `$, the terms containing $`E_3`$ are not required (effectively $`d\mathrm{}`$), and the other pieces of the $`\tau `$-functions collapse to:

$$\tau _\pm =1\pm 2\mathrm{cos}(2t\mathrm{sin}\rho )e^{2x\mathrm{cos}\rho }\frac{1}{\mathrm{tan}\rho }\sqrt{\frac{\epsilon +\mathrm{cos}\rho }{\epsilon \mathrm{cos}\rho }}e^{4x\mathrm{cos}\rho }\left(\frac{\epsilon +\mathrm{cos}\rho }{\epsilon \mathrm{cos}\rho }\right).$$
(3.8)
The expression (3.8) is quite easy to analyze. In particular, the solutions ought to be real and have no singularities in the region $`x<0`$. It was pointed out in  that these requirements are met provided the parameters $`\epsilon `$ and $`\rho `$ are suitably restricted:

$$1<\epsilon <0\text{ and }\mathrm{cos}\rho <\epsilon .$$
(3.9)
It is interesting that the upper limit $`\mathrm{cos}\rho =\epsilon `$ corresponds to a choice of frequency at which the amplitude of the solution has collapsed to zero yielding the vacuum configuration, $`\varphi =0`$. This is quite unlike the bulk breathers of the sine-Gordon theory which approach the vacuum as their frequencies approach zero. On the other hand, this behaviour of the boundary breathers is closer to that of a standard harmonic oscillator of a given frequency whose classical amplitude may be arbitrarily small. Once quantized, the energy spectrum of the oscillator depends only on its frequency and the zero-point energy is there to remind us that an oscillator of arbitrarily small amplitude is distinct from the vacuum. Similarly, we expect that the sinh-Gordon boundary breathers will have a non-zero ground state energy.

For general boundary conditions the properties of the solutions are less amenable to analysis. However, numerical investigation of the general boundary breathers indicates that they are non-singular in the region $`x<0`$ provided

$$\mathrm{cos}\frac{\pi (a_0+a_1)}{2}<0,0<\mathrm{cos}\rho <\mathrm{cos}\frac{\pi (a_0+a_1)}{2},\mathrm{cos}\frac{\pi (a_0a_1)}{2}>0.$$
(3.10)
Again, the breathers have frequencies bounded below because the parameters are restricted. We shall take the parameters to satisfy $`a_0a_1`$ (since one must be larger than the other unless they are equal, we take the larger to be $`a_0`$), and to lie within the region defined by:

$$0a_0,a_11,0a_{}\frac{1}{2},\frac{1}{2}a_+1,\pi (1a_+)\rho \frac{\pi }{2}.$$
(3.11)
As we have noted, the boundary breathers for boundary conditions preserving the symmetry of the sinh-Gordon equation are included as the special case $`a_0=a_1`$.

The energy functional of the sinh-Gordon model incorporating the boundary condition (2.4) is given by

$`[\varphi ]=`$ $`{\displaystyle _{\mathrm{}}^0}𝑑x\left({\displaystyle \frac{1}{2}}\dot{\varphi }^2+{\displaystyle \frac{1}{2}}\varphi _{}^{}{}_{}{}^{2}+{\displaystyle \frac{2}{\beta ^2}}\left(\mathrm{cosh}(\sqrt{2}\beta \varphi )1\right)\right)`$
$`+{\displaystyle \frac{2}{\beta ^2}}\left(\epsilon _0(e^{\frac{\beta }{\sqrt{2}}\varphi (0,t)}1)+\epsilon _1(e^{\frac{\beta }{\sqrt{2}}\varphi (0,t)}1)\right),`$ (3.12)
but it is most easily calculated in terms of the $`\tau `$ functions as a boundary term ,

$$[\varphi ]=\frac{2}{\beta ^2}\left(\epsilon _0\left(\frac{\tau _+}{\tau _{}}1\right)+\epsilon _1\left(\frac{\tau _{}}{\tau _+}1\right)\left(\frac{\tau _+^{}}{\tau _+}+\frac{\tau _{}^{}}{\tau _{}}\right)\right)|_{\mathrm{x}=0}.$$
(3.13)
Using this, the energies of the boundary breathers were calculated in  and are

$$=\frac{4}{\beta ^2}\left(22\mathrm{cos}\rho +\left(\mathrm{sin}\frac{\pi a_0}{2}+\mathrm{sin}\frac{\pi a_1}{2}\right)^2\right).$$
(3.14)
For the symmetric case with $`a_0=a_1`$ the energy simplifies to

$$_{\text{breather}}=\frac{8}{\beta ^2}(\mathrm{cos}\rho \epsilon ).$$
(3.15)
In the quantum field theory, the continuum of boundary breather solutions is expected to lead to a discrete spectrum of boundary bound states. To estimate this spectrum we shall follow the method described in  which adapts the techniques developed by Dashen, Hasslacher and Neveu (DHN) . Although it is not obvious that this method gives an exact result, it is clear that the result is non-perturbative, in the sense of being an all orders computation in perturbation theory in terms of the bulk coupling constant and the boundary parameters.

One of the ingredients to the DHN prescription is the classical action computed over a single period of the boundary breather. This quantity is relatively straightforward to calculate for the symmetric boundary condition but, as far as we can see, is not tractable analytically for the general breather. Nevertheless, we have conjectured the result and checked it numerically in two different ways. The principal difficulty lies in calculating the integral representing the kinetic energy over a single period. We maintain this should be given by a simple expression linear in $`a_+`$, that is

$$I=_0^T𝑑t_{\mathrm{}}^0𝑑x\dot{\varphi }^2=\frac{8\pi }{\beta ^2}(\rho \pi (1a_+)),$$
(3.16)
where $`T=\pi /\mathrm{sin}\rho `$ is the period of the general breather. When $`a_0=a_1`$, eq(3.16) agrees with the result given in ; using Maple it is possible to check that the integral is independent of $`a_{}`$ by differentiating the integrand with respect to $`a_{}`$ and calculating the resulting integral numerically, obtaining $`0`$; it is possible to check directly, again by numerical integration within the ranges of parameters (3.10)), that $`I`$ depends linearly on both $`\rho `$ and $`a_+`$. Some details of these computations will be given in appendix A. The conjectured expression (3.16) seems to us to be an astonishing result for which we are quite unable to find an analytical derivation, despite its beguiling simplicity.

## 4 The boundary bootstrap

For certain ranges of the parameters $`E`$ and $`F`$ the particle reflection amplitude (2.5) has simple poles at particular imaginary values of $`\theta `$ on the physical strip, $`0<\text{Im}(\theta )<\pi /2`$. These are expected to be due to the propagation of virtual excited boundary states although there are also other potential explanations via generalized Coleman-Thun mechanisms . The reason for this is the following. Once the boundary condition breaks the bulk symmetry of the model, the lowest energy field configuration is no longer $`\varphi =0`$. This means that in the bulk there will be effective odd-point couplings in addition to the standard even-point couplings  and these may be used to construct loop Feynman diagrams ‘attached’ to the boundary, some of which may generate poles in the reflection factor $`K_q`$. This possibility is difficult to analyze and we are unable to pursue it here. We shall simply assume the poles are due to virtual excited boundary bound states without concerning ourselves with their dynamical origin. The amplitudes for the reflection of the sinh-Gordon particle from these excited boundary states are obtained by the boundary bootstrap . When the reflection factor (2.5) has a pole at $`\theta =i\psi `$ with $`0<\psi <\pi /2`$ then the reflection factor corresponding to the associated excited boundary state is calculated via the relation

$$K_1(\theta )=K_0(\theta )S(\theta i\psi )S(\theta +i\psi ),$$
(4.1)
where $`S(\theta )`$ is the two-particle S-matrix (2.2) and $`K_0(\theta )`$ is the ground state reflection factor. Also, since energy is conserved, the energy of the excited boundary state relative to the ground state is given by

$$_1=_0+m(\beta )\mathrm{cos}\psi ,$$
(4.2)
where $`m(\beta )`$ is the mass of the sinh-Gordon particle.

As a parenthetical remark it is worth mentioning the consequences of this bootstrap procedure for a free, real scalar field of mass $`m`$ with a linear boundary condition at $`x=0`$:

$$_x\varphi |_0=m\lambda \varphi .$$

Its $`S`$ matrix is unity but the reflection factor is given by

$$K=\frac{i\mathrm{sinh}\theta +\lambda }{i\mathrm{sinh}\theta \lambda }=\frac{1}{(1+2a)(12a)},$$

where $`\lambda =\mathrm{cos}a\pi `$ in the second formula. This has a pole which indicates (provided $`m<\lambda <0`$ or $`1>a>1/2`$) a boundary breather of a fixed frequency $`\omega `$, given by $`\omega ^2=m^2\lambda ^2=m^2\mathrm{sin}^2a\pi `$. Explicitly, the appropriate normalizeable solution to the Klein-Gordon equation and the boundary condition is $`\varphi =Ae^{\lambda x}\mathrm{cos}\omega t`$. Repeated application of the bootstrap equations leads to a tower of boundary states, each with the same reflection factor, but with energies given by

$$_n=_0+n\omega .$$

This is just as one would expect for a harmonic oscillator attached to the boundary although an alternative dynamical argument would be needed to determine the ground state energy $`_0=\omega /2`$.

Returning to the sinh-Gordon case, let us begin by considering the classical reflection factor (2.7) in the light of the boundary breather parameter restrictions (3.10). It is immediately clear that because of the restrictions on the parameters only one of the four factors in the denominator of (2.7), namely $`(1a_0a_1)`$ has a pole in the physical strip, (remember, we have taken the principal values $`1>a_+>1/2,1/2>a_{}>0`$). We may also recall the perturbative result reported in  which implies that at least for sufficiently small coupling $`E`$ and $`F`$ do not roam far from their classical values. Bearing in mind the classical limits (2.9), this remark suggests that we may consider the poles associated with $`E`$ alone, ignoring $`F`$.

The remainder of the discussion in this section follows closely the calculations reported in  but, for the sake of completeness and small changes of notation, they will be repeated here. The regions in $`E`$ where the amplitude (2.5) has poles on the physical strip are

$$\mathrm{I}:2>E>1\text{and}\mathrm{II}:2<E<1;$$
(4.3)
since $`0B2`$, the other factors never have poles in the physical strip. In region I, $`\psi =\pi (E1)/2`$ and, using (4.1), we derive the reflection factor for the first excited state,

$$K_1=\frac{K_S}{(1E)(1+E)}\frac{(1+E+B)(1EB)}{(1E+B)(1+EB)},$$
(4.4)
where the ‘spectator’ factors have been lumped together in

$$K_S=\frac{(1)(1+B/2)(2B/2)}{(1F)(1+F)}.$$
(4.5)
The reflection factor (4.4) in turn has a new pole at $`\psi =\pi (E1B)/2`$, provided $`B<E1`$, indicating another excited state whose reflection factor is

$$K_2=\frac{K_S}{(1E+B)(1+EB)}\frac{(1+E+B)(1EB)}{(1E+2B)(1+E2B)}.$$
(4.6)
Continuing the procedure leads to a set of excited states with associated reflection factors given by,

$$K_n=\frac{K_S}{(1E+(n1)B)(1+E(n1)B)}\frac{(1+E+B)(1EB)}{(1E+nB)(1+EnB)}.$$
(4.7)
Note, the pole corresponding to the ($`n+1`$)st state will be within the correct range provided $`E`$ satisfies $`2>E>1+nB`$. Thus, for a given $`E`$ and $`B`$ there can be at most a finite number of bound states, and possibly none. Note too that the reflection factor for scattering from the $`n`$th bound state also contains a pole corresponding to the ($`n1`$)st bound state. Since the factor $`K_S`$ is merely a spectator, at no stage will poles with $`F`$-dependent positions enter the game if they did not do so at the start.

The energies $`_n`$ of the boundary states are found by repeatedly applying (4.2). They are given by,

$$_{n+1}=_n+m(\beta )\mathrm{cos}\frac{\pi }{2}(nBE+1).$$
(4.8)
This is the result that we will compare with the quantization of the classical breather spectrum in order to determine how $`E,F`$ and $`m(\beta )`$ depend upon $`a_0`$ and $`a_1`$.

The poles in region II do not represent a new set of states. In the discussion section we shall make some comments concerning the relative rôles of $`E`$ and $`F`$.

## 5 Semi-classical quantization

The first step in carrying out the semi-classical calculation is to solve the sinh-Gordon equation linearised in the presence of the boundary breathers. Setting $`\varphi =\varphi _0+\eta `$, where $`\varphi _o`$ is a classical breather, the linear wave equations which ought to be satisfied by the fluctuations are:

$$\frac{^2\eta }{t^2}\frac{^2\eta }{x^2}+4\eta \mathrm{cosh}\sqrt{2}\beta \varphi _0=0,\left(\frac{\eta }{x}+\eta \left[\epsilon _0e^{\beta \varphi _0/\sqrt{2}}+\epsilon _1e^{\beta \varphi _0/\sqrt{2}}\right]\right)_{\mathrm{x}=0}=0.$$
(5.1)
It is convenient to solve (5.1) by perturbing (3.1). In other words, we may take

$$\eta =\frac{\tau _{}\delta \tau _+\tau _+\delta \tau _{}}{\tau _+\tau _{}},$$
(5.2)
with $`\delta \tau _\pm `$ chosen in turn by adding a pair of ‘small’ exponential terms to (3.2) as follows. Take $`e_1`$ and $`e_2`$ defined by

$$e_1=\lambda _1e^{i\omega t+ikx},e_2=\lambda _2e^{i\omega tikx},\omega ^2k^2=4,$$

where $`\lambda _1,\lambda _2`$ are infinitesimally small. We use Hirota’s method again, but with five basic ingredients this time, instead of the previous three, keeping only the terms linear in $`e_1`$ and $`e_2`$. Thus, we write,

$`\delta \tau _\pm `$ $`=`$ $`{\displaystyle \underset{p=1,2}{}}e_p(\pm 1+(\mu _{p1}E_1+\mu _{p2}E_2+\mu _{p3}E_3)`$ (5.3)
$`\pm (\mu _{p1}\mu _{p2}A_{12}E_1E_2+\mu _{p1}\mu _{p3}A_{13}E_1E_3+\mu _{p2}\mu _{p3}A_{23}E_2E_3)`$
$`+\mu _{p1}\mu _{p2}\mu _{p3}A_{12}A_{13}A_{23}E_1E_2E_3).`$
The new coefficients $`\mu _{pq}`$ are obtained from the general formulae (3.3) and given by:

$`\mu _{11}`$ $`=`$ $`{\displaystyle \frac{1}{\mu _{22}}}={\displaystyle \frac{ik\mathrm{cos}\rho \omega \mathrm{sin}\rho 2}{ik\mathrm{cos}\rho \omega \mathrm{sin}\rho +2}},\mu _{12}={\displaystyle \frac{1}{\mu _{21}}}={\displaystyle \frac{ik\mathrm{cos}\rho +\omega \mathrm{sin}\rho 2}{ik\mathrm{cos}\rho +\omega \mathrm{sin}\rho +2}},`$
$`\mu _{13}={\displaystyle \frac{1}{\mu _{23}}}={\displaystyle \frac{ik2}{ik+2}}.`$
Matching the boundary condition at $`x=0`$ fixes the ratio $`\lambda _2/\lambda _1`$ to be ddd Remark: the right hand side of the corresponding result in , eq(5.5), has a misprint and the printed formula ought to be inverted to obtain the correct result for $`K_B`$

$$K_B=\frac{\lambda _2}{\lambda _1}=\left(\frac{ik2\mathrm{cos}\rho }{ik+2\mathrm{cos}\rho }\right)^2\frac{ik2}{ik+2}\frac{ik2c_+}{ik+2c_+}\frac{ik+2c_{}}{ik2c_{}},$$
(5.5)
where $`c_\pm =\mathrm{cos}\pi (a_0\pm a_1)/2`$. The result (5.5) is not intended to be obvious; its derivation relied heavily on using a symbolic algebraic manipulation language—in our case, Maple. In the limit $`x\mathrm{}`$, the fluctuation $`\eta `$ is a superposition of left and right moving plane waves,

$$\eta \lambda _1e^{i\omega t}\left(e^{ikx}+K_Be^{ikx}\right),$$
(5.6)
and the relative phase of these waves defines the classical reflection factor corresponding to the boundary breather. Taking $`\mathrm{cos}\rho =c_+`$, the breather collapses to the static ground state solution and the reflection factor collapses to

$$K_0=\frac{ik2}{ik+2}\frac{ik+2c_+}{ik2c_+}\frac{ik+2c_{}}{ik2c_{}},$$
(5.7)
which, with $`k=2\mathrm{sinh}\theta `$, is simply an alternative form (ie without using the ‘block’ notation) of the expression (2.7).

The period $`T=\pi /\mathrm{sin}\rho `$ of the boundary breather defines the ‘stability angles’ via

$$\eta (t+T,x)=e^{i\nu }\eta (t,x)e^{i\omega T}\eta (t,x)$$
(5.8)
and the field theoretical version of the WKB approximation makes use of the stability angles together with a regulator to calculate a quantum action. If there are no boundaries, it is natural to add some artificially to render the spectrum discrete and facilitate the necessary calculations. One way, the simplest and most commonly used, would be to place the field theory in an interval $`[L,L]`$ with periodic boundary conditions and to manipulate the sum over the discrete stability angles. However, since we have one boundary already prescribed, we need to do something else as suggested in . It is convenient to treat the sinh-Gordon model in the interval $`[L,0]`$ and to impose the Dirichlet condition $`\eta (t,L)=0`$. Since the limit $`L\mathrm{}`$ will be taken eventually, the stability angles for the boundary breather ($`\nu _B`$), or the vacuum solution ($`\nu _0`$) are effectively determined by the reflection factors given in (5.6) or (5.7), respectively. A potentially more interesting calculation but one which we do not yet have the machinery to carry out, would be to consider two sets of two-parameter boundary conditions. However, the breathers we are considering are inadequate for that.

Following  we need to calculate a sum over the stability angles and use it to correct the classical action. Thus,

$$\mathrm{\Delta }=\frac{1}{2}(\nu _B\nu _0)\frac{T}{2}\left(\sqrt{k_B^2+4}\sqrt{k_0^2+4}\right),$$
(5.9)
where $`k_B`$ and $`k_0`$ are the sets of (discrete) solutions to

$`e^{2ik_BL}=\left({\displaystyle \frac{ik_B+2\mathrm{cos}\rho }{ik_B2\mathrm{cos}\rho }}\right)^2{\displaystyle \frac{ik_B+2}{ik_B2}}{\displaystyle \frac{ik_B+2c_+}{ik_B2c_+}}{\displaystyle \frac{ik_B2c_{}}{ik_B+2c_{}}},`$
$`e^{2ik_0L}={\displaystyle \frac{ik_0+2}{ik_02}}{\displaystyle \frac{ik_02c_+}{ik_0+2c_+}}{\displaystyle \frac{ik_02c_{}}{ik_0+2c_{}}}.`$ (5.10)
Once $`\mathrm{\Delta }`$ is known, the quantum action is defined by

$$S_{\mathrm{qu}}=S_{\mathrm{cl}}\mathrm{\Delta },$$
(5.11)
where the classical action is readily calculated from the kinetic energy integral (3.16) and the total energy (3.14):

$$S_{\mathrm{cl}}=_0^T𝑑t_{\mathrm{}}^0𝑑x=\frac{8\pi }{\beta ^2}\left(\rho \pi (1a_+)+\frac{1}{\mathrm{sin}\rho }\left(\mathrm{cos}\rho +1\frac{1}{2}(1c_+)(1+c_{})\right)\right).$$
(5.12)
We shall proceed along the lines described in  noting that for large $`k`$ the solutions to either of (5) are close to

$$k_n=\left(n+\frac{1}{2}\right)\frac{\pi }{L},$$

and so it is reasonable to set $`(k_B)_n=(k_0)_n+\kappa ((k_0)_n)/L`$ where, for $`L`$ large, the function $`\kappa `$ is given approximately by

$$e^{2i\kappa (k)}=\left(\frac{ik+2\mathrm{cos}\rho }{ik2\mathrm{cos}\rho }\frac{ik+2c_+}{ik2c_+}\right)^2.$$
(5.13)
Interestingly, all dependence upon $`c_{}`$ has dropped out and therefore, from this point on, the calculation is identical to the corresponding part of the calculation presented in  with $`\epsilon =\mathrm{cos}\pi a`$ replaced by $`c_+`$.

In terms of $`\kappa `$ the expression (5.9) is rewritten

$$\mathrm{\Delta }\frac{T}{2L}\underset{n0}{}\frac{(k_0)_n\kappa ((k_0)_n)}{\sqrt{(k_0)_n^2+4}}+O(1/L^2),$$

and this, in turn, as $`L\mathrm{}`$ can be converted to a convenient (but actually divergent) integral,

$$\mathrm{\Delta }=\frac{T}{2\pi }_0^{\mathrm{}}𝑑k\frac{k\kappa (k)}{\sqrt{k^2+4}}$$
(5.14)
with which we shall have to deal. Note that $`\kappa `$ vanishes when $`\mathrm{cos}\rho =c_+`$. Integrating (5.14) by parts we find

$$\mathrm{\Delta }=\frac{T}{2\pi }\left(\kappa \sqrt{k^2+4}|_0^{\mathrm{}}_0^{\mathrm{}}𝑑k\frac{d\kappa }{dk}\sqrt{k^2+4}\right),$$
(5.15)
where

$$\frac{d\kappa }{dk}=\frac{4\mathrm{cos}\rho }{k^2+4\mathrm{cos}^2\rho }+\frac{4c_+}{k^2+4c_+^2},$$
(5.16)
and we note that with a suitable choice of branch

$$\kappa \frac{4\mathrm{cos}\rho }{k}\frac{4c_+}{k}\text{as}k\mathrm{}.$$
(5.17)
Using (5.17) and recalling that $`\mathrm{cos}\rho <c_+`$, we deduce that $`\kappa `$ approaches zero from above as $`k\mathrm{}`$. Also, from (5.16) it is clear that the derivative of $`\kappa `$ is positive near $`k=0`$ but negative as $`k\mathrm{}`$. Hence, the first term in (5.15) is well-defined and the appropriate branch of $`\kappa `$ has $`\kappa (0)=0`$. On the other hand, the derivative of $`\kappa `$ is not decaying sufficiently rapidly to ensure the second term in (5.15) is finite. However, this was to be expected since a perturbative analysis of the sinh-Gordon model confined to a half-line needs mass and boundary counter terms to remove logarithmic divergences (which would be removed automatically by normal-ordering the products of fields in the bulk theory). With this in mind, the integral remaining in (5.15) should be replaced by

$$_0^{\mathrm{}}𝑑k\sqrt{k^2+4}\left(\frac{4\mathrm{cos}\rho }{k^2+4\mathrm{cos}^2\rho }\frac{4\mathrm{cos}\rho }{k^2+4}+\frac{4c_+}{k^2+4c_+^2}\frac{4c_+}{k^2+4}\right),$$
(5.18)
the first counter-term removing the bulk divergence and the second being there to remove a similar divergence associated with the boundary. In effect, we are regarding the parameter $`a`$ as describing the bare coupling which appears in the boundary part of the Lagrangian once it is written in terms of normal-ordered products of fields. The counter-terms clearly respect the symmetry and the whole expression vanishes when $`\rho =\pi (1a_+)`$. The integrals in (5.18) need to be treated carefully noting that $`\mathrm{cos}\rho >0`$ but $`c_+<0`$.

Besides the towers of real solutions to (5), there is also a discrete set of solutions for which $`k_0`$ and $`k_B`$ are pure imaginary. These were not discussed in  but including them in the argument would not have altered the conclusions, as we shall show. First, notice that as $`L\mathrm{}`$ the imaginary solutions which survive are either zeroes or poles of the right hand sides of the equations in (5), according to the signs of $`ik_0`$ or $`ik_B`$. Thus for $`k_0`$ we have $`ik_0=2,2c_+`$if $`ik_0>0`$ and $`ik_0=2,2c_+`$ if $`ik_0<0`$, while for $`k_B`$ we have $`ik_B=2\mathrm{cos}\rho ,2`$ if $`ik_B>0`$ and $`ik_B=2\mathrm{cos}\rho ,2`$ if $`ik_B<0`$. However, the signs should be disregarded because for either sign each solution for a specific value of $`|ik_0|`$ or $`|ik_B|`$ represents a single ‘bound-state’ function $`\eta `$. Taking this into account, these special solutions contribute to $`\mathrm{\Delta }`$ an additional piece:

$$T(\mathrm{sin}\rho \mathrm{sin}\pi a_+)=\pi \frac{\pi \mathrm{sin}\pi a_+}{\mathrm{sin}\rho }.$$
(5.19)
Assembling the various components leads to

$$\mathrm{\Delta }=\pi \frac{2}{\mathrm{sin}\rho }\left(\mathrm{cos}\rho +\mathrm{cos}\pi a_++\rho \mathrm{sin}\rho +\pi (a_+1/2)\mathrm{sin}\pi a_+\right).$$
(5.20)
Recalling (5.12), and using (5.20), the quantum action defined in (5.11) is given by an expression of the form

$$S_{\mathrm{qu}}=\frac{4}{B}\left(\frac{\mathrm{cos}\rho }{\mathrm{sin}\rho }+\rho \frac{\pi }{2}\right)+\frac{8\pi }{\beta ^2}\left(\pi a_+\frac{\pi }{2}\right)+\frac{\mathrm{\Gamma }(a_+,a_{})}{\mathrm{sin}\rho },$$
(5.21)
where $`\mathrm{\Gamma }`$ is independent of $`\rho `$,

$$\mathrm{\Gamma }(a_+,a_{})=2\pi \left(\frac{2}{\beta ^2}(1c_+)(1+c_{})+\frac{4}{\beta ^2}+\frac{c_+}{\pi }+(a_+1/2)\mathrm{sin}\pi a_+\right).$$

Once the quantum action is determined, the quantum energy is defined by

$`_{\mathrm{qu}}`$ $`=`$ $`{\displaystyle \frac{S_{\mathrm{qu}}}{T}}={\displaystyle \frac{\mathrm{sin}^2\rho }{\pi \mathrm{cos}\rho }}{\displaystyle \frac{S_{\mathrm{qu}}}{\rho }}={\displaystyle \frac{4}{\pi B}}\mathrm{cos}\rho {\displaystyle \frac{\mathrm{\Gamma }(a_+,a_{})}{\pi }},`$ (5.22)
and the WKB quantization condition states that

$$W_{\mathrm{qu}}=S_{\mathrm{qu}}+T_{\mathrm{qu}}=\frac{4}{B}\left(\rho \frac{\pi }{2}\right)+\frac{8\pi }{\beta ^2}\left(\pi a_+\frac{\pi }{2}\right)=2n\pi .$$
(5.23)
Where, $`n`$ is a positive integer or zero.eee Including the imaginary solutions of (5) obviates the need for making the change $`nn+1/2`$, since we shall see below the zero-point energy is automatically correct. In this sense, including the imaginary solutions leads more naturally to the expected result. Hence, the energies of the quantized boundary breather states are determined by a set of special angles $`\rho _n`$,

$$\rho _n=\frac{\pi }{2}\left(1+nB\frac{2\pi B}{\beta ^2}(2a_+1)\right),$$
(5.24)
and given by

$$_n=\frac{4}{\pi B}\mathrm{cos}\rho _n\frac{\mathrm{\Gamma }}{\pi }=\frac{4}{\pi B}\mathrm{cos}\frac{\pi }{2}\left(nB+1\frac{2\pi B}{\beta ^2}(2a_+1)\right)\frac{\mathrm{\Gamma }}{\pi }.$$
(5.25)
Notice that as $`\beta 0`$, $`\rho _n\pi (1a_+)`$ independently of $`n`$. Thus, the frequencies collapse to the lowest allowed frequency, namely $`\omega _0=2\mathrm{sin}\pi a_+`$. On the other hand, in the same limit the energies are independent of $`\beta `$ and non-zero,

$$_n\left(n+\frac{1}{2}\right)\omega _0.$$
(5.26)
This is precisely the spectrum of a harmonic oscillator vibrating at the fundamental frequency $`\omega _0`$.

Using (5.25) the corresponding differences in the energy levels are given by

$$_{n+1}=_n+\frac{8}{\pi B}\mathrm{sin}\frac{\pi B}{4}\mathrm{cos}\frac{\pi }{2}\left(\frac{2\pi B}{\beta ^2}(2a_+1)\left(n+\frac{1}{2}\right)B\right).$$
(5.27)
Comparing (5.27) with the outcome of the bootstrap calculation (4.8) allows us to identify the parameter $`E`$ which appeared in the expression for the reflection factor (2.5). Thus, from the first excited level we deduce,

$$E(a_0,a_1,\beta )=\frac{2\pi B}{\beta ^2}(2a_+1)\frac{1}{2}B+12a_+\left(1\frac{B}{2}\right)=(a_0+a_1)\left(1\frac{B}{2}\right),$$
(5.28)
and then all the other levels match up without further restriction. This is in agreement with the suggestions made in .

As was pointed out previously in , the comparison with (4.8) also permits us to deduce an expression for $`m(\beta )`$, the mass of the sinh-Gordon particle:

$$m(\beta )=\frac{8}{\pi B}\mathrm{sin}\frac{\pi B}{4}.$$
(5.29)
## 6 Discussion

From the point of view of the classical reflection factor (2.7), it is clear that the parametrisation (2.8) is the most appropriate. Moreover, the expressions for $`\epsilon _0`$ and $`\epsilon _1`$ in terms of $`a_0`$ and $`a_1`$ are invariant under independent changes of signs of either $`a_0`$ or $`a_1`$ and this symmetry is incorporated in the factors appearing in the denominator of (2.7). We expect the symmetry under reversing the signs of $`a_0`$ and $`a_1`$ should persist in the quantum reflection factor, and indeed that was part of the thinking behind the conjectured forms for $`E`$ and $`F`$ presented in (2.10). In other words, having calculated $`E`$ in (5.28) we can immediately deduce the partner expression for $`F`$.

We found that the general breathers exist as non-singular real solutions provided the parameters satisfy the constraints (3.11). It is gratifying to discover at the end of the calculation that the renormalization factor $`1B/2`$ in (2.10) for real coupling $`\beta `$ is restricted to lie in the range

$$01\frac{B}{2}1,$$

and thus has the effect of scaling the parameters but preserving the constraints (3.11). Notice that with $`E`$ and $`F`$ given by (2.10) it is impossible to have $`E`$ and $`F`$ lying simultaneously in regions (4.3) for which the bound state poles of the reflection factor lie on the physical strip. Thus, the poles depending on $`E`$ or $`F`$ are alternatives and actually lead to the same tower of states simply seen differently in terms of the parameters $`a_0`$ and $`a_1`$.

Changing the sign of either $`a_0`$ or $`a_1`$ takes us outside the principal ranges of these parameters and obliges us to reformulate the breather solutions. In fact, there are no others besides the ones we have found although the expressions for them will be a little different if we adopt a different principal region.

It is worth emphasizing that the expressions (2.10) incorporate a weak-strong coupling duality  which extends that enjoyed by the S-matrix itself. If a new triple of coupling constants $`(a_0^{},a_1^{},\beta ^{})`$ is defined by

$$(a_0^{},a_1^{},\beta ^{})=\frac{4\pi }{\beta ^2}(a_0,a_1,\beta ),$$

then

$$E(a_0^{},a_1^{},\beta ^{})=\frac{4\pi }{\beta ^2}(a_0+a_1)\frac{B}{2}(a_0+a_1)(1B/2)=E(a_0,a_1,\beta ),$$
(6.1)
and similarly for $`F`$.

We must acknowledge that the semi-classical approximation, although non-perturbative, is not guaranteed to be exact. However, for the bulk sine-Gordon model, the Dashen, Hasslacher and Neveu approach does give exact information concerning the spectrum and we would expect that the same should be true here. Unfortunately, alternative exact computations of the spectrum are not yet available for comparisons to be made.

As far as other models are concerned, it will be interesting to try out these methods in other cases where the boundary parameters form a continuous set. The simplest such example, in which there are two distinct one-parameter families of boundary parameters, is the model based on the $`a_2^{(2)}`$ root data .

Acknowledgements

One of us (AT) is supported by a Leverhulme Fellowship, the other (EC) thanks the University of Durham for a Visiting Professorship and the Yukawa Institute for Theoretical Physics for hospitality during a visit funded by a Royal Society/JSPS exchange programme; both of us have been partially supported by a TMR Network grant of the European Commission contract number FMRX-CT-960012 and we thank the Universities of Montpellier and Mons for their hospitality. We are indebted to Nigel Glover for advice on using Maple and to Gustav Delius for discussion and comments.

## Appendix A Appendix

In this appendix we shall give some information concerning the kinetic energy integral (3.16). The first step in evaluating the integral is to perform the integration over time. This is relatively easy and leads to:

$$I(x)=_0^T𝑑t\dot{\varphi }^2=\frac{8\pi \mathrm{sin}\rho }{\beta ^2}\left[\frac{2C\sqrt{B^2A^2}}{ADBC}+\frac{B}{\sqrt{B^2A^2}}\frac{2A\sqrt{D^2C^2}}{ADBC}+\frac{D}{\sqrt{D^2C^2}}\right],$$
(A.1)
where the right hand side is built from the components $`A,B,C,D`$ which are defined by

$`A={\displaystyle \frac{su}{q}}(1q^2)(1rv),C={\displaystyle \frac{su}{q}}(1q^2)(1+rv),`$
$`B={\displaystyle \frac{1}{q^2}}\left((rvs^2q^2u^2)+q^2(1s^2q^2rvu^2)\right),`$
$`D={\displaystyle \frac{1}{q^2}}\left((rv+s^2q^2u^2)+q^2(1+s^2q^2rvu^2)\right),`$ (A.2)
with

$$u=e^{2x\mathrm{cos}\rho },v=e^{2x}.$$
(A.3)
The parameters $`q,r,s`$ were defined previously in (3.7). Notice, when $`a_0=a_1`$ we have $`r=0`$ allowing all of these expressions to simplify dramatically. Consequently, the integral of the right hand side of (A.1) can be done by making a suitable change of variables. Thus, for $`a_0=a_1=a`$ we find:

$$_{\mathrm{}}^0𝑑xI(x)=\frac{8\pi }{\beta ^2}(\rho \pi (1a)).$$
(A.4)
However, in the general case we have not found a change of variables which simplifies (A.1). Therefore, we have had to proceed numerically.

One check we have made is to differentiate (A.1) with respect to $`a_{}`$ and demonstrate

$$\frac{}{a_{}}_{\mathrm{}}^0𝑑xI(x)=0,$$
(A.5)
provided $`a_\pm `$ lie within the ranges (3.11). Another is to calculate the integral directly as a function of $`a_+`$ and $`\rho `$ and demonstrate it depends linearly on each of these parameters separately. The following plots generated by Maple indicate unambiguously the linear dependence on $`a_+`$ (Figure 1)

Figure 1: The integral as a function of $`a_+`$.

or on $`\rho `$ (Figure 2)

Figure 2: The integral as a function of $`\rho `$.

In Fig(1) we have overlaid several plots. The dots represent the numerical values of the integral as $`a_+`$ varies from $`1\rho /\pi `$ to $`1a_{}`$, in accordance with (3.11), for each of five different values of $`\rho `$ (viz. $`\rho /\pi =(2+k)/16,k=1,\mathrm{},5)`$; the continuous line is simply the conjectured value of the integral given in (3.16) plotted as a function of $`a_+`$ for the same five values of $`\rho `$. In all plots, the factor $`8\pi /\beta ^2`$ has been ignored. The vertical line indicates the upper bound on $`a_+`$ for $`a_{}=1/8`$. Fig(2) presents similar plots for the integral as a function of $`\rho /\pi `$ in the range $`1a_+`$ to $`1/2`$ for five different values of $`a_+`$ (viz. $`a_+=(8+k)/16,k=1,\mathrm{},5)`$. We regard the sets of parallel lines depicted in Figs(1,2) as very convincing numerical evidence for (3.16).


format: text

text:
# Modelling the spectral energy distribution of galaxies. I. Radiation fields and grain heating in the edge-on spiral NGC891

## 1 Introduction

Dust grains can be considered as test particles for the intrinsic radiation field in galaxies. Observations of their emission in the infrared (IR), combined with optical and ultraviolet (UV) observations of the light from stars, attenuated by the grains, should, in principle, strongly constrain the intrinsic distribution of stellar luminosity and dust in galaxies. This would address the fundamental question of optical thickness in galactic disks and allow evaluation of optical and IR observational data in terms of physical (e.g., star formation history) rather than empirical (e.g., colours) parameters. Ultimately, such analysis techniques will be necessary to interpret the evolution of the UV - far-infrared (FIR)/sub-millimeter (sub-mm) spectral energy distribution (SED) of galaxies and the IR background radiation over cosmological timescales.

However, even for galaxies in the local universe, technical difficulties are considerable from both theoretical and observational stand points. From the theoretical point of view, it is necessary to solve the radiation transfer problem for a finite disk, and self consistently calculate the stochastic heating of grains placed in the resulting radiation field using realistic grain models. To that must be added the problem that disks are fundamentally inhomogeneous on different scales. From the observational point of view, it is only now becoming possible to obtain resolved images covering the mid-IR to sub-mm spectral range for nearby galaxies.

Because of these difficulties, the origin of the FIR emission in local universe spiral galaxies is still a subject of debate. Based on models of dust heating in various environments (e.g., Mezger et al. 1982; Cox et al. 1986), and from comparison of IRAS data with optical properties of galaxies, several studies concluded that a substantial fraction of the FIR emission from galaxies could be due to dust heated by the diffuse interstellar radiation field (ISRF), the remaining part being due to localised sources around OB stars in star-forming regions (e.g., de Jong et al. 1984; Walterbos & Schwering 1987). Some later papers, however, suggested that the correlation between H$`\alpha `$ and FIR emission implies that OB stars by themselves can account for all the dust heating, so that there is no need for an additional contribution to heating by the ISRF (e.g., Devereux & Young 1990). By contrast, there have been a series of papers (e.g., Xu et al. 1994) which concluded that the non-ionising UV accounts for most of the energy absorbed by dust. They showed that the FIR/radio correlation can be separated into a warm-FIR/thermal-radio correlation which is due to massive ionising stars ($`>20\mathrm{M}_{}`$), and a cool-FIR/nonthermal-radio correlation, principally due to intermediate mass stars ($`520\mathrm{M}_{}`$) - the supernova progenitors. A heating component from the old stellar population was also needed to account for the non-linearity of the correlation, but this accounts only for about a third of the heating.

In this paper we describe a new method for the analysis of the radiation energy budget in edge-on spiral galaxies, as part of an effort to develop a theoretical tool for the interpretation of the SED in spiral galaxies. The method presented here is based on some knowledge of the primary radiation field (i.e., from stars) in spiral galaxies and the dust distribution in them. For edge-on spirals, it has been demonstrated (Xilouris et al. 1997; 1998; 1999) that it is possible to determine the distribution of older stars and associated dust from analysis of optical and near-infrared (NIR) images. Thus, the diffuse radiation field in the galaxy from these optical bands can be computed quite accurately. The diffuse radiation field in the UV is generally not so well constrained by the UV data alone, due to the larger optical depth of disks in the UV. This younger stellar component is parameterised in terms of a recent star-formation rate (SFR) and it can be observationally constrained in the case that a direct, optically thin indicator, such as radio free-free emission, is available. Solving the radiation transfer problem within the disk and bulge, and calculating the heating of grains we can determine the diffuse FIR/sub-mm emission. Potentially, further contributions arising from emission from localised sources need to be added to the diffuse component. Then, because the dust is optically thin in the FIR, the predicted emission can be compared with observations.

We apply the above method to the well-known edge-on spiral galaxy NGC891, for which we assume a distance of 9.5 Mpc (van der Kruit & Searle 1981). This is one of the most extensively observed edge-on galaxies in the nearby universe, which makes it ideal for a verification of our modelling technique. The analysis of a sample of 6 more edge-on galaxies is in progress (Misiriotis et al. 2000a). Modelling edge-on spiral galaxies has several advantages, mainly when investigating them in the optical band. One advantage is that, in this view of a galaxy, one can easily separate the three main components of the galaxy (i.e., the stellar disk, the dust and the bulge). Another advantage is that the dust is very prominently seen in the dust lane, and thus its scalelength and scaleheight can be better constrained. A third advantage is that many details of a galaxy that are evident when the galaxy is seen face-on (e.g., spiral arms), are smeared out to a large degree when the galaxy is seen edge-on (Misiriotis et al. 2000b). Thus, a simple model with relatively few parameters can be used for the distribution of stars and dust in the galaxy. However, the third advantage comes with disadvantages, especially when trying to model the galaxy in all wavelength ranges, including UV and FIR/sub-mm. Thus, in edge-on galaxies it is very difficult to see localised sources (i.e., HII regions), in which the radiation can be locally absorbed and thus not contribute to the diffuse radiation field. Also, if the galaxy has a thin (young) stellar/dust disk, highly obscured by the dust lane in the plane of the galaxy, then this disk cannot be inferred from observations in the optical/NIR spectral range. In passing we mention that throughout this paper we will use the terms “thin/thick disks” to describe their scaleheights, and not their optical thickness.

The goal of this paper is to find the simplest phenomenological description of NGC891 that can adequately account for the observed optical to sub-mm SED. In particular we investigate what fraction of the FIR is accounted for by the Xilouris et al. (1999) model, namely by a disk of dust, heated by a diffuse, smooth, optically emitting, old stellar disk and bulge, supplemented by a diffuse, smooth, UV emitting disk of newly-formed stars. We refer to this model as the “standard model”. We will show that this model fails to reproduce the observed SED of the galaxy at FIR/sub-mm wavelengths. We investigate then different possibilities to explain the origin of the missing FIR/sub-mm flux.

Recently there have been several works modelling the SED of galaxies from the UV to the sub-mm (Silva et al. 1998, Devriendt et al. 1999). Their method consists of using models for photometric and/or spectrophotometric and chemical evolution of galaxies in order to fit the observed SED. For example Silva et al. (1998) used a chemical evolution code to follow the SFR, the gas fraction, and the metallicity of the galaxies. In their approach the parameters describing the star formation history as well as the geometrical parameters of the intrinsic distribution of stars and dust are left as fitting parameters. While this approach was very successful in fitting the SED of different types of galaxies, it implies the use of many free parameters. Furthermore, a more detailed description of the geometry of the galaxy is indispensable for an accurate determination of the dust content in individual cases.

This paper is organised as follows: In Sect. 2 we describe our “standard model” for the diffuse dust and radiation field, including the optical properties of the dust grains, and the derivation of the diffuse radiation field powered by the old and the young stellar populations. In Sect. 3 we give the results of our calculations, in terms of integrated SED for the “standard model” and in Sect. 4 we supplement our model with the contribution of localised sources. Since the “standard model” fails to reproduce the observed SED, we discuss in Sect. 5 the origin of the FIR emission, considering four possibilities to produce the missing IR component. From these, the two dust disk model, together with the clumpy scenario constitute the simplest solutions we have been able to identify to describe the origin of the missing infrared component. In Sect. 6 we give our results in terms of radial profiles and in Sect. 7 we discuss the contribution of radiation at different wavelengths in heating the dust. In Sect. 8 we summarise our results and give our conclusions.

## 2 The model

### 2.1 The properties of the dust grains

We consider only the spectral range $`\lambda 40\mu `$m as our dust model does not include Polycyclic Aromatic Hydrocarbon (PAH) molecules, which are important for $`\lambda 15\mu `$m and may even contribute at $`\lambda =25\mu `$m, a point that is still very uncertain. For example, Désert et al. (1990) proposed a dust model in which 48% of the IR emission at 25 $`\mu `$m (for dust in the solar neighbourhood) is emitted by the PAH molecules as continuum emission. Draine and Anderson (1985) proposed that the IR emission at 25 $`\mu `$m and shorter wavelengths can be explained by enhancing the number of small grains. Our restriction to the spectral range $`\lambda 60\mu `$m allows us to avoid the uncertainties introduced by a different dust component (PAHs or modified grain size distribution).

In characterising the dust properties we thus considered spherical “astronomical grains” for two component materials - graphite and silicate - of the Mathis et al. (1977; hereafter MRN) interstellar grain model. MRN have proposed that the observed interstellar extinction of star-light may be produced by a mixture of graphite and silicate particles with a simple power-law distribution of sizes

$`dN_i=N_i(a)da=N_ia^kda`$ (1)
where $`dN_i`$ is the number of grains of type $`i`$ with radii in the interval $`[a,a+da]`$, and $`k=3.5`$. We consider the lower and the upper cutoffs to the size distribution to be respectively $`a_{\mathrm{min}}=10\mathrm{\AA }`$ and $`a_{\mathrm{max}}=0.25\mu `$m. The graphite and silicate abundances were taken from Draine & Lee (1984, hereafter DL), and are only slightly different from those found by MRN; DL proposed a mixture of 53% silicates ($`N_{\mathrm{Si}}`$) and 47% graphites ($`N_{\mathrm{Graphite}}`$), which were chosen to fit the extinction curve in our Galaxy, and which we also adopted for NGC891.

The absorption efficiencies $`Q_{\mathrm{abs}}(\lambda ,a)`$, the scattering efficiencies $`Q_{\mathrm{scat}}(\lambda ,a)`$ and the scattering phase function $`g(\lambda ,a)`$ were taken from Laor & Draine (1993). The absorption coefficient $`\kappa _{\mathrm{abs}}`$, the scattering coefficient $`\kappa _{\mathrm{scat}}`$ as well as the averaged anisotropy parameter of the Henyey-Greenstein scattering phase function $`g`$ can then be derived by integrating $`Q_{\mathrm{abs}}(\lambda ,a)`$, $`Q_{\mathrm{scat}}(\lambda ,a)`$, and $`g(\lambda ,a)`$, respectively, over the size distribution

$`\kappa _{\mathrm{abs}}(\lambda )`$ $``$ $`{\displaystyle _{a_{\mathrm{min}}}^{a_{\mathrm{max}}}}{\displaystyle \underset{i}{}}N_\mathrm{i}a^\mathrm{k}da\pi a^2Q_{\mathrm{abs}}(\lambda ,a)`$ (2)
$`\kappa _{\mathrm{scat}}(\lambda )`$ $``$ $`{\displaystyle _{a_{\mathrm{min}}}^{a_{\mathrm{max}}}}{\displaystyle \underset{i}{}}N_\mathrm{i}a^\mathrm{k}da\pi a^2Q_{\mathrm{scat}}(\lambda ,a)`$ (3)
$`g(\lambda )`$ $`=`$ $`{\displaystyle \frac{{\displaystyle _{a_{\mathrm{min}}}^{a_{\mathrm{max}}}}{\displaystyle \underset{i}{}}N_\mathrm{i}a^\mathrm{k}da\pi a^2g(\lambda ,a)Q_{\mathrm{scat}}(\lambda ,a)}{{\displaystyle _{a_{\mathrm{min}}}^{a_{\mathrm{max}}}}{\displaystyle \underset{i}{}}N_ia^kda\pi a^2Q_{\mathrm{scat}}(\lambda ,a)}}`$ (4)
with the extinction coefficient $`\kappa _{\mathrm{ext}}=\kappa _{\mathrm{abs}}+\kappa _{\mathrm{scat}}`$.

As the $`60100\mu `$m data are expected to contain a significant amount of emission from small grains not in equilibrium with the radiation field, it is necessary to model the stochastic emission from the dust. For this calculation we adopted the heat capacities of silicate grains from Guhathakurta & Draine (1989), which were derived as a fit to experimental results for SiO2 and obsidian at temperatures $`10<T<300`$ K (Leger et al. 1985), with a simple extrapolation for $`T>300`$ K. For graphite grains the heat capacities were taken from Dwek (1986).

### 2.2 Derivation of the diffuse optical and UV radiation field.

Our first step is to estimate the diffuse radiation field at any point in the galaxy. The most straightforward way to do this is to solve the radiative transfer equation (RTE) for a given distribution of emitters and absorbers of radiation. In this section we will consider the distribution of emitters and absorbers as derived from optical wavelengths by Xilouris et al. (1998,1999), while at the UV wavelengths we will fix the parameters describing the geometry of the emitters and use population synthesis models to parameterise the amplitude of the UV radiation.

#### 2.2.1 Radiation transfer

In this subsection we will describe our solution of the RTE. Our method is based on that of Kylafis & Bahcall (1987). This method will be briefly described here for the sake of completeness. Let us introduce the RTE (Mihalas 1978) for a time independent radiation field along a line of sight, at a fixed wavelength

$$\frac{dI(s,\widehat{n})}{ds}=\kappa _{ext}I(s,\widehat{n})+\eta _0(s,\widehat{n})+\eta _s(s,\widehat{n})$$
(5)
where $`I(s,\widehat{n})`$ is the specific intensity of the radiation field at a point $`s`$ on the line of sight (defined by $`\widehat{n}`$), $`\kappa _{ext}`$ is the extinction coefficient for the wavelength considered, $`ds`$ is an infinitesimal length element along the line of sight, $`\eta _0`$ is the emissivity along the line of sight due to primary sources (stars) and $`\eta _s`$ is the emissivity due to scattering into the line of sight. For coherent scattering, $`\eta _s`$ is given by

$$\eta _s(s,\widehat{n})=\omega \kappa _{ext}I(s,\widehat{n}^{})p(\widehat{n},\widehat{n}^{})\frac{d\mathrm{\Omega }^{}}{4\pi }$$
(6)
where $`p(\widehat{n},\widehat{n}^{})`$ is the Henyey-Greenstein phase function (Henyey & Greenstein 1941), parameterised by the anisotropy parameter $`g`$, and $`\omega `$ is the scattering albedo for the wavelength considered, defined as $`\omega \kappa _{scat}/\kappa _{ext}`$. Let $`I=I_0+I_1+\mathrm{}+I_n+\mathrm{}`$, where $`I_0`$ is the intensity of photons that come directly from the primary sources, $`I_1`$ the intensity of photons that have been scattered once after leaving the primary sources and $`I_n`$ the intensity of photons that have suffered $`n`$ scatterings. The intensities $`I_n`$ satisfy the following equations (Henyey 1937)

$$\frac{dI_n(s,\widehat{n})}{ds}=\kappa _{ext}I_n(s,\widehat{n})+\eta _n(s,\widehat{n})$$
(7)
where $`\eta _n(s,\widehat{n})`$ is the emissivity due to photons scattered $`n`$ times, and is given by

$$\eta _n(s,\widehat{n})=\omega \kappa _{ext}I_{n1}(s,\widehat{n}^{})p(\widehat{n},\widehat{n}^{})\frac{d\mathrm{\Omega }^{}}{4\pi }$$
(8)
Eq. (7) can be solved for as high order of $`n`$ as allowed by the given hardware and/or software. In our method, we calculate $`I_0`$ and $`I_1`$ and approximate the higher orders according to Kylafis & Bahcall (1987). Therefore, the solutions are given by

$$I_{0,1}(s,\widehat{n})=_0^s𝑑s^{}\eta _{0,1}(s^{},\widehat{n})\mathrm{exp}[\tau (s,s^{})]$$
(9)
while for $`n>1`$,

$$I_{n+1}=I_n\frac{I_1}{I_0}$$
(10)
where $`\tau (s,s^{})`$ is the optical depth between $`s`$ and $`s^{}`$. Given the intensity of the radiation field at a point, the energy density can be calculated by (Mihalas 1978)

$$u(x,y,z)=c^1I(x,y,z,\widehat{n})𝑑\mathrm{\Omega }$$
(11)
#### 2.2.2 Optical radiation field and dust distribution

The determination of the intrinsic stellar content of a galaxy is often difficult because of the uncertainty introduced by dust obscuration. However, in the case of edge-on galaxies, dust obscuration can be used to our benefit to determine the spatial distribution of dust. In an edge-on galaxy, dust becomes visible as a very prominent dust lane along the major axis of the galaxy. The effects of the dust are strong enough to allow the existence of a fairly constrained model for the distribution of old stars and associated dust (Xilouris et al. 1997; 1998;1999; Kuchinski et al. 1998; Ohta & Kodaira 1995). This does not include, however, any stellar population with a scaleheight small enough to be effectively completely hidden by the dust. Such a stellar population will be discussed in Sect. 2.2.3.

NGC891 was modeled (Xilouris et al. 1998,1999) by fitting an artificial image produced by radiative transfer on observations in an attempt to conclude on the dust morphology. In this section we will adopt the latest model for NGC891 (Xilouris et al. 1999). In this model, the emissivity is described by an exponential disk and a de Vaucouleurs bulge, while the dust resides in a pure exponential disk. The parameters describing the stars and the dust are determined in 5 wavebands (B, V, I, J, K).

For the stellar emissivity the following formula is used

$$\eta _0(R,z)=L_s\mathrm{exp}\left(\frac{R}{h_s}\frac{|z|}{z_s}\right)$$

$$+L_b\mathrm{exp}(7.67B^{1/4})B^{7/8}$$
(12)
In this expression the first part describes an exponential disk, and the second part describes the bulge, which in projection is the well-known $`R^{1/4}`$-law (Christensen 1990). Here $`R`$ and $`z`$ are the cylindrical coordinates, $`L_s`$ is the stellar emissivity per unit volume per steradian at the center of the disk and $`h_s`$ and $`z_s`$ are the scalelength and scaleheight respectively of the stars in the disk. For the bulge, $`L_b`$ is a normalisation constant, while $`B`$ is defined by

$$B=\frac{\sqrt{R^2+z^2(a/b)^2}}{R_e}$$
(13)
with $`R_e`$ being the effective radius of the bulge and $`a`$ and $`b`$ being the semi-major and semi-minor axis respectively of the bulge. In this paper we will refer to the exponential disk described above as to the old stellar disk, in order to differentiate it from the young stellar disk emitting in UV (Sect. 2.2.3).

For the dust distribution a similar expression as that adopted for the stellar distribution in the disk is used, namely

$$\kappa _{ext}(\lambda ,R,z)=\kappa _{ext}(\lambda ,0,0)\mathrm{exp}\left(\frac{R}{h_d}\frac{|z|}{z_d}\right)$$
(14)
where $`\kappa _{ext}(\lambda ,0,0)`$ is the extinction coefficient at the center of the disk for the wavelength considered, and $`h_d`$ and $`z_d`$ are the scalelength and scaleheight respectively of the dust. The central optical depth of the model galaxy seen face-on is

$`\tau ^f(\lambda )=2\kappa _{ext}(\lambda ,0,0)z_d`$ (15)
The parameters determined for NGC891 are presented in Table 1. The total amount of dust derived from this model is $`M_{dust}=5.6\times 10^7\mathrm{M}_{}`$ (Xilouris et al. (1999).

The extinction coefficients $`\kappa _{ext}(\lambda ,0,0)`$ were derived as fitting parameters (Xilouris et al. 1999), separately for each of the five wavebands used for modelling the galaxy. They were then used to calculate the absorption of the optical light by dust, via the radiation transfer procedure. Since in calculating the emission of the dust we used the absorption efficiencies ($`Q_{abs}`$) taken from Laor & Draine (1993), as presented in Sect. 2.1, we compare in Fig. 1 the extinction coefficients derived from the fitting routine with the theoretical ones derived from the extinction efficiencies (Eqs.  and ). The result of such a comparison shows that indeed the extinction derived from observations (square symbols) is very close to the theoretical one used to calculate the IR emission (solid line). This ensures that the treatment of the absorption and of the emission is consistent. In Fig. 1 we also give the wavelength dependence of the absorption and scattering coefficients, plotted with dashed and dotted lines, respectively.

#### 2.2.3 UV radiation field

The UV radiation field cannot be derived directly from observations. For edge-on galaxies it is extremely difficult to infer something about the young stellar disk, which, due to it’s small scaleheight, is presumably highly obscured by line-of sight dust. Any young stellar population whose scaleheight is smaller than the scaleheight of the diffuse dust will be totally unprobed by UV, optical and NIR data on an edge-on galaxy. Thus for an edge-on galaxy mid to far-IR emission would be the primary observational signature for any young massive star population having a scaleheight of 100 pc or less. This stellar population would (we assume luminosity is dominated by stars) produce almost all the non-ionising and ionising UV from the galactic disk.

We thus consider the UV radiation field as a parameter and we parameterise the UV luminosity in terms of a recent SFR, based on the population synthesis models of Mateu & Bruzual (2000). We consider $`Z=\mathrm{Z}_{}`$, a Salpeter IMF, $`M_{up}=100\mathrm{M}_{}`$ and $`\tau =5`$ Gyr. This gives, for example, SFR $`=8.55\times 10^{28}L_{UV}`$ \[erg/s/Hz\] for $`\lambda =912`$Å, SFR $`=3.74\times 10^{43}L_{UV}`$ \[erg/s\] for the integrated ionising UV shortwards of 912 Å, SFR $`=1.18\times 10^{28}L_{UV}`$ \[erg/s/Hz\] at $`\lambda =1500`$ Å, etc. The diffuse UV radiation derived from the calibration above was confined in a thin disk, with a scaleheight $`z=90`$ pc (close to that of the Milky Way, Mihalas & Binney 1981) and the same scalelength as the blue disk.

The UV luminosity could, in principle, also be indirectly inferred from different indicators of the SFR, like the free-free emission or H$`\alpha `$ emission. There are however large uncertainties in deriving these quantities, especially in the case of NGC891. The first indicator, the radio free-free emission, has the advantage of being optically thin. Niklas et al. (1997) found that the total radio emission of NGC891 is dominated by synchrotron emission, while the free-free emission constitutes a thermal fraction of just 0.05 at an observing frequency of 10 GHz. However, the radio morphology of NGC891 is rather complicated, with more morphological components having different spectral indices. Allen et al. (1978) suggested the existence of two components for the radio continuum emission: a highly flattened thin disk component coinciding with the equatorial plane of the galaxy and a thick disk or halo component (see also Hummel et al. 1991). If this is the case, then the free-free emission estimated from a spectral decomposition of the integrated fluxes of the galaxy (as done by Niklas et al. 1997) could have been seriously underestimated. We conclude that the free-free emission does not constitute an accurate constraint of our model.

The second indicator of the SFR is the H$`\alpha `$ emission. It is known that correcting H$`\alpha `$ data for internal extinction is a very uncertain step, since the extinction contains a large-scale component plus extinction due to localised sources, which are difficult to account for in edge-on disks. Hoopes et al. (1999) found that in NGC891, $`8386\%`$ of the observed H$`\alpha `$ emission is in the diffuse component associated with the halo of the galaxy. This would imply that most of the emission from the disk is highly obscured, making unreliable any derivation of the SFR in the disk.

### 2.3 Calculation of the IR spectrum from grains embedded in the diffuse optical and UV radiation field

The dust heating due to photons from the diffuse interstellar radiation field was calculated following the method of Guhathakurta & Draine (1989). This method derives the temperature distribution $`P(a,T)`$ of various grain radii $`a`$ as a function of dust temperature $`T`$. Due to stochastic heating, small grains undergo significant fluctuations from the equilibrium temperature, while larger grains have narrower probability functions, eventually approaching delta functions. We calculated the temperature distribution at any point in the galaxy, based on the energy density of the radiation field derived in the previous subsections. The spectral shape of the UV-NIR energy density and its magnitude will obviously determine the shape of the temperature distributions. The grain sizes used in the calculations were chosen between $`a_{min}=10`$Å and $`a_{max}=0.25\mu `$m, for a step $`\mathrm{\Delta }\mathrm{log}a=0.05`$. In Fig. 2 we give the temperature distribution as $`\mathrm{d}p/\mathrm{dlog}T`$ for some silicate (Fig. 2a,c) and graphite (Fig. 2b,d) grain sizes embedded in the diffuse radiation field, for a solution with the UV radiation field in the diffuse component corresponding to a SFR $`=3\mathrm{M}_{}`$/yr. Fig. 2 also illustrates the dependence of the grain temperature on the position in the disk. We give two examples, for two extreme positions in the plane of the disk: in the center of the disk (Fig. 2a,b) and at the edge of the disk (Fig. 2c,d). At the edge of the disk the probability distribution converges towards lower equilibrium temperatures, as expected due to the decrease in the intensity of the radiation field. We also note the relatively large difference in the equilibrium temperature between $`30`$ K for the central region and $`15`$ K for the outer disk. This is interesting in the context of the results of Alton et al. (1998) who fitted the SED of NGC891 with 15 and 30 K grey-body components. We discuss this further in Sect. 7, in terms of the geometry and opacity of the disk. We also note that in passing from the centre to the outer regions, progressively bigger grains exhibit stochastic heating, as expected due to the lower energy density of the radiation field. For the same location in the disk, graphite grains undergo larger fluctuations than the silicate grains.

Given the temperature probability function $`P(a,T)`$, the IR flux is given by

$`F_\nu `$ $`=`$ $`{\displaystyle \frac{1}{d^2}}{\displaystyle _{a_{\mathrm{min}}}^{a_{\mathrm{max}}}}N(a)da\pi a^2Q_{\mathrm{abs}}(\nu ,a)`$ (16)
$`\times {\displaystyle _0^{\mathrm{}}}B_\nu (T)P(a,T)\mathrm{d}T`$
where $`B_\nu `$ is the Planck function, $`N(a)`$ is the grain size distribution given by Eq. (1), and $`d`$ is the distance to the galaxy.

The IR emission was calculated in each point of the galaxy and then the total emission was obtained by integrating the emission over the galaxy. In Fig. 3 we give two examples of IR spectra for two extreme locations in the plane of the galactic disk, namely in the center of the galaxy ($`R=0`$ kpc, $`z=0`$ kpc) and at the edge of the disk ($`R=15`$ kpc, $`z=0`$ kpc). At the edge of the disk the infrared emission is colder, as expected due to the decrease in the dust temperature, though the increased importance of stochastic heating in the outer regions tends to flatten the SED on the Wien side.

## 3 Results for the “standard model”

The FIR emission for our “standard model” was calculated for a grid of SFR, and the different solutions were compared with the observations. To extract the FIR luminosity of the galactic dust we integrated the SED of NGC891 between 40 and 1000$`\mu `$m, with the data taken from Alton et al. (1998). The SED was parameterised by Alton et al. (1998) by a fitting two-temperature greybody curve (with a wavelength dependence of grain emissivity in the FIR $`\beta =2`$) to the observed fluxes from IRAS and SCUBA at 60, 100, 450 and 850 $`\mu `$m. This yielded best-fit “temperatures” of 30 and 15 K to the observed spectrum, and an integrated luminosity between 40 and 1000$`\mu `$m of $`5.81\times 10^{36}`$ W. This output luminosity can be accommodated in our “standard model” with a solution for a SFR $`=7.5\mathrm{M}_{}`$/yr. The predicted FIR SED from the model is given in Fig. 4a, together with the UV-optical-NIR intrinsic emitted stellar radiation (as would have been observed in the absence of dust). The IRAS and SCUBA data (diamonds) from Alton et al. (1998) are also given for comparison. The inferred SFR is an upper limit, since we assumed that all the non-ionising UV escapes into the diffuse component. Indeed, in the absence of local absorption, all the lines of sight from each source would contribute primary photons to the diffuse interstellar radiation field. In practice, some fraction of the lines of sight - call this fraction $`F`$ \- from each source will be opaque to local absorption, giving rise to discrete IR sources embedded in a diffuse IR emitting disk. Our solution was thus calculated under the extreme assumption that $`F=0`$ for the non-ionising UV and $`F=1`$ for the ionising UV. While requiring a solution with a large SFR, the predicted SED is slightly warmer than the observed SED at 100 $`\mu `$m. However, the main difficulty of the model is that it underestimates by a factor of 2 to 4 the observed sub-mm fluxes.

## 4 Results for the “standard model” supplemented by localised sources

A more realistic solution is to consider that some fraction $`F`$ of the non-ionising UV is locally absorbed in star-forming complexes (e.g., HII regions). Such a solution requires two free parameters ($`F`$ and SFR) and can be accommodated with a lower SFR, due to the higher probability of absorption of the non-ionising UV photons. To include the contribution of HII regions we used a template SED of such forming complexes. We utilised the galactic UC HII region G45.12+0.13, with FIR observations at 60 and 100 $`\mu `$m from the IRAS Point Source Catalog and 1300 $`\mu `$m measurements from Chini et al. (1986). We fitted the observed data with a greybody curve (with a wavelength dependence of grain emissivity in the FIR $`\beta =2`$) and used it as a template SED for the compact source component for our modelling. We do not attempt here to include in this template potential cold emission components that might be expected from “parent” molecular clouds in juxstaposition to their “offspring” HII regions.

The best solution was obtained for a SFR $`=3.5\mathrm{M}_{}`$/yr and $`F=0.28`$. The predicted FIR SED is given in Fig. 4b, again together with the UV-optical-NIR intrinsic stellar radiation. In Fig. 4b we also give the contribution of the diffuse component (dashed line) and of the HII component (dotted line) to the the total predicted FIR SED (solid line). The integrated emitted luminosity predicted by the model is $`L=4.79\times 10^{36}`$ W, which accounts for $`82\%`$ of the total observed emitted FIR luminosity. Of this, the diffuse component contributes $`46\%`$ of the total observed dust luminosity and the localised sources $`36\%`$. However, the main difficulty of the model is again that the solution fails to reproduce the observations longwards of 100 $`\mu `$m.

## 5 The missing FIR/sub-mm component of our “standard model”

We have shown that our “standard model” fails to reproduce the observed SED of NGC891 in the sub-mm spectral range. The more realistic case, where the contribution from HII regions is added to the diffuse component, is also not able to predict the observed FIR/sub-mm SED. However, this may be, to some extent, because we did not include a cold dust emission component that might be expected to arise from associated parent molecular clouds.

Below we consider four more ways that might account for this discrepancy. Firstly, we discuss the effect of altering the geometry of the large-scale distribution of stars and dust. Alternatively, we consider extraplanar emission. Thirdly, the predicted apparent sub-mm/optical ratio could be increased by including more grains in the disk, but distributed in clumps, so as not to effect the optical extinction properties of the diffuse disk. Lastly, a more radical solution is presented, namely a model with two dust-disk components, which is a moderately optically thick solution.

### 5.1 The effect of altering the geometry of the large-scale distribution of stars and dust

The diffuse IR emission predicted by our model was derived for an exponential disk of stars and dust, adopted from Xilouris et al. (1999). Galactic disks are however known to be quite complex systems, where the large-scale distribution of stars and dust presents inhomogeneities in the form of spiral structures. Adding logarithmic spiral arms as a perturbation on the exponential disks for stars and dust would alter the solutions of the radiative transfer, and thus of the calculated diffuse IR emission. However, Misiriotis et al. (2000b) showed that plain exponential disk models give a very accurate description for galactic disks seen edge-on, with only small deviations in parameter values from the real ones (typically a few percent). Thus, we conclude that the geometry adopted here cannot be responsible to account for the discrepancy between our model predictions and the observations.

### 5.2 Extraplanar dust?

Observations of external galaxies have revealed extensive thickened layers of ionised gas traced by its H$`\alpha `$ emission in several edge-on spiral galaxies (e.g., Rand et al. 1990, 1992; Dettmar 1990; Pildis et al. 1994; Rand 1996, Ferguson et al. 1996; Hoopes et al. 1999). Amongst the most spectacular examples of extraplanar diffuse ionising gas (DIG) is that seen in NGC891. The recent observations of Hoopes et al. (1999) showed that the DIG layer in NGC891 extends out to at least 5 kpc from the plane, and possibly as far as 7 kpc, and has the brightest and largest DIG layer known.

In addition to detecting the high z DIG studied in earlier work, Howk & Savage (2000) also detected individual dust-bearing clouds observable to heights $`z2`$ kpc from the midplane. Moreover, they detected the presence of HII regions at large distances from the midplane (0.6-2 kpc), which suggests that on-going star formation may be present in some of the dense, high-z clouds.

In principle both the dust-bearing clouds ($`z2`$kpc) seen in absorption as well as the dust that may be associated with the DIG ($`z5`$kpc) could contribute to the FIR/sub-mm emission of the galaxy. The possibility of the extraplanar dust clouds emitting in the sub-mm was recently investigated by Alton et al. (2000b). From the upper limits in the total amount of extraplanar dust these authors found that less than $`5\%`$ of galactic dust exists outside the galactic disk, if the dust grains are not colder than 17 K, or $`9\%`$, if the dust temperature is 10 K. This small percentage cannot account for our discrepancy between the model and the observations.

We also consider it unlikely that any diffuse dust component associated with the DIG ($`z5`$kpc) could provide a substantial fraction of the sub-mm emission. In that case, if the extraplanar emission were powered by the absorption of disk photons, the optical depth of the DIG would have to be comparable to that of the disk, which would be a very extreme scenario, probably also in conflict with the constraints on submillimeter extraplanar emission.

### 5.3 The case of hidden dust in clumps

One possible explanation to account for the missing dust component could be that Xilouris et al. (1999) underestimated the dust content of the galaxy because of the adopted diffuse distribution for the dust. Since the galaxy is optically thin, any additional dust component would have to be distributed in very dense clumps instead of being diffuse. Only then would a significant fraction of the lines of sight avoid the additional dust, so as not to affect the derived optical-NIR optical depth.

Such dust clumps could either be a component of star-formation regions, or could have no associated sources - we refer to the latter hypothesis as the “quiescent clumps”. The quiescent clumps must be optically thick to the diffuse UV/optical radiation field in the disk to have an impact on the predicted submillimeter emission. They would radiate at predominantly longer wavelengths than the diffuse disk emission, in the FIR/sub-mm spectral range. One can speculate that such optically thick “quiescent clouds” could be physically identified with partially or wholly collapsed clouds that, for lack of a trigger, have not (yet) begun to form stars. However, due to the lack of intrinsic sources, a very substantial mass of dust would have to be associated with the quiescent clouds to account for the submillimeter emission missing from our model prediction.

In reality there must be also dark clouds associated with star-forming complexes. In the Milky Way HII regions around newly born massive stars are commonly seen in juxtaposition to parent molecular clouds (e.g., M17). This is thought to be a consequence of the fragmentation of the clouds due to mechanical energy input from the winds of the massive stars. Thus, warm dust emission from cloud surfaces directly illuminated by massive stars can be seen along a fraction of the lines of sight, together with cold sub-mm dust emission from the interior of the associated optically thick cloud fragments. These dense clumps (with or without associated sources) may contribute to the sub-mm emission, and thus supplement the contribution of compact HII regions.

However, the question remains whether it is possible to place such a substantial amount of additional dust in clumps without affecting the opacity of the galaxy. Kuhinski et al. (1998) modeled a sample of highly inclined galaxies using both smooth and clumpy dust and concluded that clumpiness does not effect dramatically the opacity of highly inclined galaxies. On the other hand, Bianchi et al. (1999) adopt different parameterisation for clumps and they report that clumpiness does significantly effect the opacity of edge-on galaxies. It is beyond the scope of this paper to calculate a quantitative solution for a clumpy distribution, but we qualitatively consider the clumpy scenario, as quiescent clumps or associated with star-formation regions - to be a possibility to account for the missing FIR-sub-mm component.

### 5.4 A two-dust-disk model

We have seen that the observed sub-mm SED requires more dust in the galaxy than predicted by our “standard model”. We have already discussed the possibility of including this extra dust in clumps. An alternative solution is obtained if it is postulated that this additional dust is confined in a diffuse thin disk, in which the young stellar population is embedded. (A thin young stellar disk, but with no associated dust, was considered in Xilouris et al. 1998.) If this second dust-disk component has a scaleheight comparable with the 90 pc scaleheight of molecular gas clouds in the Milky Way, it would not have been seen in the optical/NIR analysis of Xilouris et al. (1999), being totally obscured by the optically visible disk of dust. Indeed, since the galaxy is almost edge-on, adding a second thin disk of dust will probably not change the parameters for the intrinsic distributions of stars and dust derived by Xilouris et al. (1999). However, adding a second disk of dust will change the optical depth of the disk, transforming our solution into a moderately optically thick solution ($`\tau _V^f=3.1`$).

We have calculated the energy density of the radiation field from UV to NIR for such a two-dust-disk model, and subsequently the infrared emission. The best solution was obtained for a SFR $`=3.8\mathrm{M}_{}`$/yr, $`F=0.22`$ and $`M_{dust}=7\times 10^7\mathrm{M}_{}`$ in the second disk of dust. The corresponding non-ionising UV luminosity for this SFR is $`8.2\times 10^{36}`$ W. This could be provided, for example, by a population of B5 stars ($`T_{\mathrm{eff}}=15500`$ K) with a space density in the centre of the galaxy of $`7.5\times 10^4\mathrm{pc}^3`$. The predicted combined SED from the two diffuse dust-disk components is shown by the dashed line in Fig. 5. By adding the further contribution from localised sources (dotted line), as discussed in Sect. 4, we were able to fit the SED of NGC891, from the FIR to the sub-mm. The luminosity of the diffuse component is $`4.07\times 10^{36}`$ W, which accounts for $`69\%`$ of the observed FIR luminosity, and the luminosity of the HII component is $`1.82\times 10^{36}`$ W, making up the remaining $`31\%`$ of the FIR luminosity. This model can thus successfully fit the shape of the SED. In the next section we will show that the predictions of the two-dust-disk model are in agreement with the observed radial profile at 850 $`\mu `$m. Nevertheless, further modelling of face-on galaxies will be required to test the validity of the two-dust-disk model. Ultimately, the two-dust-disk hypothesis will be directly tested observationally using the new generation of sub-mm interferometers, which will resolve edge-on disks.

## 6 The radial profiles

Since we have calculated the 3 dimensional FIR radiation field in NGC891, for both our “standard model” and the two-dust-disk model, it is straightforward to integrate along the line of sight and thus produce a 2-dimensional map of the galaxy at different wavelengths, and then radial profiles. Here we consider only the diffuse component of the IR emission. Because observed radial profiles were derived by Alton et al. (2000a) using the SCUBA observations at $`850\mu `$m, we first attempt to calculate the radial profiles at this wavelength and compare it with the observations. The calculations were done at high resolution and afterwards smoothed to 16′′ resolution. Since the observed radial profile at $`850\mu `$m from Alton et al. (2000a) were derived as averages over a bin width of 36′′, and with a sampling of 3′′, we have also calculated an averaged radial profile in the same way as Alton et al. (2000a) did, and directly compare it with the observations. We have found that in the case of the two dust disk model there is a very good agreement between the model predictions and the observations and this comparison is illustrated in Fig. 6. The predicted radial profile can be traced out to 300 arcsec radius (15 kpc), as also detected by the SCUBA.

In Figs. 7a and 7c we give our prediction for the radial profile of NGC891 at 850 $`\mu `$m, this time as it would be observed within a beam of 16′′, and without any averaging. This is shown for our “standard model” and the two-dust-disk model, respectively. In Fig. 7b.d we give the radial colour profiles at the wavelengths used by the IRAS, ISO and SCUBA, normalised to the radial profile at 850 $`\mu `$m. The strong colour gradients, seen in particular between the 60 or 100 $`\mu `$m and the sub-mm regime, are a consequence of the strong gradients in the radiation field calculated in our finite disk model. It is interesting to note that the stronger contribution of stochastic emission from small grains (relative to big grains) in the outer disk at 60 $`\mu `$m, does not strongly influence the gradients. Potentially, therefore, large-scale colour gradients in galaxy disks, could act as a diagnostic for the relative contribution of the diffuse emission and compact sources, as the latter would not be expected to have IR spectra strongly dependent on position in the disk. As expected, the 450/850 colour shows only a shallow gradient, as both these wavelengths are displaced well longwards of the spectral peak, and thus depend only approximately linearly on grain temperature.

## 7 Discussion: the contribution of different stellar populations in heating the dust in NGC891

The relative contribution of optical and UV photons in heating the dust has been a long standing question in the literature. Since we have a detailed calculation of the absorbed energy over the whole spectral range and at each position in the galaxy, we can directly calculate which part of the emitted FIR luminosity from each volume element of the galaxy is due to the optical and NIR photons, and which part is due to the UV photons. The IR emission from any given volume element arises from a common local grain population, so that the IR colours from that volume element will be the same for the emitted IR photons regardless of whether the absorbed energy had an optical or UV origin. However the relative fraction of energy absorbed from the UV and optical light varies with position in the galaxy, as does the dust temperature, so that the volume-integrated IR spectrum representing re-radiated optical light will differ from that of the re-radiated UV light. In this way, some information on the stellar populations heating the grains is preserved in the IR domain, in the form of the observed IR colours. Volume-integrated IR spectral components arising from re-radiated optical and UV light are presented in Fig. 8a,b for the “standard model” and in Fig. 8c,d for the two-dust-disk model, respectively. The overall contribution of the diffuse optical and NIR radiation ($`400022000`$ Å) to the total FIR luminosity is $`1.07\times 10^{36}`$ W, or $`22\%`$ for the “standard model”. Since the predicted luminosity of this model is lower than the observed one, it is also meaningful to give the contribution of the diffuse optical radiation to the observed FIR luminosity, which is $`18\%`$. This component of the FIR SED is plotted with dashed line in Fig. 8a,b. The corresponding optically generated IR luminosity for the two-dust-disk model is $`1.85\times 10^{36}`$ W, or $`31\%`$, higher than in the “standard model”, largely because the galaxy is becoming optically thick in the optical range. This component of the FIR SED is again plotted with dashed-line in Fig. 8c,d. The contribution of different stellar components to the total FIR luminosity for the two dust models is summarised in Table 2. Table 3 gives the contribution of different stellar components to the FIR wavelength commonly used by IRAS, ISO and SCUBA.

We note that, regardless of the optical thickness of the adopted model, the diffuse optical radiation field makes only a relatively small contribution to the total emitted dust luminosity. This is in qualitative agreement with various statistical inferences linking FIR emission with young stellar populations, in particular the FIR-radio correlation which, also after normalisation to indicators of galactic mass, such as K-band flux, remains rather tight (Xu et al. 1994). Perhaps more surprising, at first glance, is the predicted predominance of UV-powered grain emission in the sub-mm range in both our dust models (Fig. 8b,d.). However, our analysis has demonstrated that the FIR colours have to be fundamentally interpreted also in geometrical terms, rather than simply as separate temperature components. Thus, the predicted increasing predominance of re-radiated diffuse UV photons over re-radiated diffuse optical photons going from the FIR into the sub-mm regime is due to an increasing fraction of the emission arising from more optically thin regions, where the ratio of optical to UV absorption is largely controlled by the relative optical and UV grain emissivities. This particularly applies in the the outer disk, which, in addition to being (relative to the inner disk) optically thin, has systematically lower grain heating and temperature.

## 8 Summary and Conclusions

In this paper we have modeled the edge-on spiral galaxy NGC891 from the UV to the FIR and sub-mm with the purpose of understanding what is the origin of the FIR/sub-mm emission. The strategy of the paper was to use a model for the intrinsic distribution of older stars and associated dust as derived from fitting the optical and NIR images of NGC891 by Xilouris et al. (1999), which, supplemented by a distribution of newly-form stars, should predict the diffuse FIR emission. By comparing the prediction of the “standard model” with the observations we tried to identify whether this model can reproduce the observed dust emission.

To calculate the FIR emission in the galaxy we needed to have a dust model and also the energy density of the radiation field at each position in the galaxy. In describing the properties of the grains we used the extinction efficiencies from Laor & Draine (1993) for spherical “astronomical grains“, for a mixture of graphites and silicates of the MRN grain size distribution. For small grains not in equilibrium with the radiation it was necessary to model the stochastic emission. For this calculation we adopted the heat capacities from Guhathakurta & Draine (1989) for silicate grains and from Dwek (1986) for graphite grains. Our theoretical dust model was found to be consistent with the extinction coefficients in the optical and NIR spectral range, derived independently by Xilouris et al. (1999) for NGC891. The similarity between the extinction coefficients derived from our theoretical model and those derived as fitting parameters assure us that we can use the model for the distribution of old stars and associated dust of Xilouris et al. (1999).

The calculation of the energy density in the galaxy was done using the solution of the radiation transfer equation, based on the method of Kylafis & Bahcall (1987). This calculation was done separately for the optical/NIR and UV radiation field. For the distribution of old and intermediate age stars we used a smooth exponential disk and a de Vaucouleurs bulge, with the parameters derived from Xilouris et al. (1999). For the distribution of newly-formed stars we used a smooth exponential disk, with fixed parameters, e.g., with the same scalelengths as that of the stars radiating in the B band, and with a scaleheight of 90 pc. The amplitude of the radiation emitted by the young stars (of the ionising and non-ionising UV radiation) was considered as a free parameter. Using the population synthesis models of Mateu & Bruzual (2000) we parameterised the UV radiation in terms of the recent SFR. Another free parameter was the fraction $`F`$ between the UV photons which are locally absorbed within star-forming complexes and those which participate in the diffuse ISRF.

The predicted SED from our “standard model” was compared with the observed one. In the extreme case that $`F=0`$, the SFR can be taken so as to exactly reproduce the total bolometric output in infrared. This is a solution with SFR $`=7.5\mathrm{M}_{}`$/yr. While requiring a high SFR, this solution fails to reproduce the sub-mm part of the spectrum. A solution with $`F>0`$, which includes also the contribution of star-forming complexes, underestimates even more the emission in the sub-mm spectral range. The failure of our “standard model” to fit the observed SED of NGC891 indicates, most probably, that there is missing dust, which could not have been detected by Xilouris et al. (1999). Two scenarios are proposed to explain this missing dust component, namely a scenario in which the missing dust is in the form of very small optically thick clumps, such that they do not affect the optical extinction derived by Xilouris et al. (1999), or in the form of a second diffuse thin disk, where the young stellar population is embedded. For the second case we give a quantitative description of the model, and we make calculations for the FIR output. We have found that the two-dust-disk model reproduces very well the observed SED, as well as the radial profile at 850 $`\mu `$m. This solution requires a SFR $`=3.5\mathrm{M}_{}`$/yr and $`F=0.28`$. This value of the SFR fits better with the preconception that NGC891 is a relatively quiescent normal galaxy.

Potentially, a further possibility exists to explain the shortfall in sub-mm emission predicted by our “standard model”. A higher level of sub-mm emission could be attained from a given grain mass if the sub-mm grain emissivity were higher than assumed in our dust model. This approach was followed by Alton et al. (2000a). However, for a given grain heating the sub-mm/FIR colour of the grain emission will become colder, in order to preserve the energy balance between the absorbed UV/optical and the radiated emission. Thus, the peak of the SED would be shifted to longer wavelengths, so in this scenario there would be a need to increase the SFR to account for the observed fluxes in the 60-100 $`\mu `$m range.

We emphasise that, on the basis of currently available observational evidence, it is difficult to distinguish between the two dust-disk-model and the clumpy scenario. As present day sub-mm and FIR telescopes generally have insufficient angular resolution to resolve the postulated geometrical emission components, the best way of distinguishing between them might be to compare the statistics of quantities derived from the UV/optical/FIR/sub-mm SEDs. Such statistics could be, for example, the derived star-formation rates from a sample of edge-on and face-on systems with comparable morphological types. With the advent of the future generations of FIR telescopes with arcsec resolution we will be able to directly distinguish between diffuse disk and cloud scenarios. Then such techniques will be needed to understand young galaxies at cosmological distances.

###### Acknowledgements.

We would like to thank the anonymous referee for his useful comments and suggestions. We kindly acknowledge Dr. Juan Mateu for providing us with unpublished results on the population synthesis models. We would also like to thank Dr. Emmanuel Xilouris for motivating discussions. This work was supported in part by Projects 50OR99140 and 50QI92014 of the Deutsches Zentrum für Luft- und Raumfahrt.