text:
# Exact bounds on the critical frustration in the Shastry–Sutherland Model

## I Introduction

We study a two-dimensional Heisenberg model with additional frustrating interactions of strength $`J_1`$ on every second diagonal bond (see Fig. 1) given by the Hamiltonian

$$H=J_2\underset{}{}𝐒_𝐢𝐒_𝐣+J_1\underset{Dimer}{}𝐒_𝐢𝐒_𝐣,$$
(1)
where $`𝐒_𝐢`$ denotes a spin operator for spin $`S`$ at site i. The model, commonly referred to as the Shastry–Sutherland model , displays a rich zero-temperature phase diagram as a function of $`J_1`$ and $`J_2`$, showing long range ferromagnetic, antiferromagnetic and helical order, as well as short range spin liquid behaviour. Also for large frustrating coupling $`J_1`$ the model has an exactly known ground state, built up of uncoupled dimers on the diagonal bonds. It was the existence of this exact ground state, which has first attracted attention to the model.

Recently the Hamiltonian (1) for $`S=\frac{1}{2}`$ has been widely discussed, because it represents essential features of SrCu2(BO3)2, a newly discovered spin gap substance . SrCu2(BO3)2 appears to have a ground state of dimer singlets and some of its unusual magnetic properties are well explained by the Shastry–Sutherland Hamiltonian.

In the classical limit ($`S=\mathrm{}`$) of the Hamiltonian (1) the nature of the ground state phases and their boundaries are easily analyzed (see Fig. 2). The model shows two long range helical phases for $`0<|J_2|<J_1`$, which are separated by an antiferromagnetic dimer phase along the line $`0=J_2<J_1`$. In the regime $`|J_2|>J_1`$ the ground state is ordered antiferromagnetically for $`J_2>0`$ and ferromagnetically for $`J_2<0`$. For $`J_2=0>J_1`$ a phase of independent ferromagnetic dimers occurs.

The ground state energy of the helical phase can be found by minimizing the energy of a cluster of only three spins. From this basic three spin entity the ground state of the whole system is constructed. Depending on the initial choice of spin directions one finds a spin helix with a four-fold degenerate direction and with a twist angle of $`\theta =\pi \pm \mathrm{arccos}(J_2/J_1)`$ between neighbouring spins.

Away from the classical limit quantum effects transmute the classical antiferromagnetic dimers for $`J_1>0`$ into unique singlet dimers and stabilize the resulting dimer phase. In what follows we will often use the variable $`x:=J_2/J_1`$, which is the inverse of the frustration. For $`S<\mathrm{}`$ and $`J_2>0`$ the dimer phase occupies a finite region of the phase diagram and at a certain critical value $`x_c^a(S)>0`$ there is a phase transition of first order  to a new phase, whose nature is controversially discussed at the moment. For even larger $`x`$ the system changes to the antiferromagnetic regime.

In the ferromagnetic region ($`J_2<0`$) the situation is slightly simpler, since for $`S=\frac{1}{2}`$ there is a first order phase transition directly from the dimer to the ferromagnetic phase at $`x=1`$. This transition marks one of the exactly known points in the phase diagram. For $`S>\frac{1}{2}`$, like in the antiferromagnetic regime, the system crosses at a negative $`x_c^f(S)`$ from the dimer to an intervening phase. Finally at $`x=1`$ the ferromagnetic region is reached for all $`S`$.

For $`J_1<0`$ the classical phase diagram is not changed by quantum effects. The ferromagnetic dimer phase for $`J_2=0`$ becomes a phase of independent spin-$`2S`$ dimers with macroscopic degeneracy. It may be interesting to note, that the spin-$`S`$ square lattice antiferromagnet for $`0=J_1<J_2`$ is continuously connected to the spin-$`2S`$ square lattice antiferromagnet as $`J_2>0>J_1\mathrm{}`$.

There are two major problems connected with the phase transition between the short range spin liquid phase and the adjacent phases. Firstly the nature of the phases adjacent to the singlet dimer phase \- both on the ferromagnetic and on the antiferromagnetic side - is not clear and a point of intense investigation at the moment. Secondly for arbitrary $`S`$ the boundary of the dimer phase cannot be determined by a simple argument like in the case $`S=\mathrm{}`$.

In this paper we concentrate on the second problem and give rigorous lower (Sec. II) and upper (Sec. III) bounds on the phase boundaries of the dimer phase by using various versions of a variational ansatz for finite clusters.

Various approximative methods have been used recently to determine this phase boundary. In Ref. the phase boundary was calculated by means of Schwinger boson mean-field theory. By the flow equation method in Ref. a value of $`x_c^a(S=\frac{1}{2})=0.63`$ was obtained and in Ref. an intervening plaquette phase was suggested for $`S=\frac{1}{2}`$ occupying the regime $`0.861>x>0.677`$. In Ref. by an extension of $`SU(2)`$ to the symplectic groups $`Sp(2N)`$ an intervening phase with helical and incommensurate order was found between a dimer phase and a region with collinear commensurate order. Also in Ref. it was found for a two-dimensional model with frustration, that the dimer and the antiferromagnetic phase are separated by an intervening regime, which is characterized as a weakly incommensurate spin density wave.

## II Lower bounds on $`x_c^a`$ and on $`x_c^f`$

To obtain lower bounds on $`x_c^a`$ we decompose the Hamiltonian $`H`$ of Eq. (1) into cluster terms $`H_i^N`$ with

$$H=\underset{i}{}H_i^N$$
(2)
and calculate the lowest eigenvalues $`E_N^0`$ of the Hamiltonian $`H_i^N`$ of the clusters (N is the number of spins in the cluster). The decomposition of $`H`$ is chosen in such a way, that the clusters cover the whole lattice without overlapping bonds. Taking the ground state of H as variational state for the finite clusters  it follows that $`E_N^0`$ is always smaller than or equal to the ground state energy $`E_{\mathrm{}}^0`$ of $`H`$. To obtain a bound on $`x_c^a`$ we choose finite clusters (see Fig. 3) which have the dimer state as an eigenstate and calculate their critical value $`x_{c,N}^a`$. It is obvious, that $`x_{c,N}^a`$ is always smaller than the $`x_c^a`$ of the infinite system.

As the simplest possible system we consider a plaquette with four spins (first entry in Fig. 3). In this case the energy of the second lowest state, which for growing x is competing with the dimer state, crosses the energy of the dimer ground state at

$$x_{c,4}^a(S)=\frac{1}{2S+1}.$$
(3)
From what was said above $`x_{c,4}^a(S)`$ is a strict lower bound for $`x_c^a`$. (Note that for $`S>1/2`$ the first estimate Eq. (3) is better than the bound $`J_1/J_2>2(1+S)`$ given in Ref..) Better limits on $`x_c^a`$ were obtained by calculating the ground state energies of the systems shown in Fig. 3 using the Lanczos method.

For a system with 31 sites we thus obtain a best lower bound of 0.5914 for $`x_c(S=\frac{1}{2})`$. We want to point out, that even for $`S=\frac{1}{2}`$, where systems of up to 31 spins were calculated, it is difficult to make a good finite size analysis with the results of the finite clusters, because not only the system size, but also the shape of the cluster influences the ground state energy. On these grounds the significance of the recent estimate of $`x_c(S=\frac{1}{2})=0.7\pm 0.01`$ in Ref., which was obtained by an extrapolation of three systems with different shapes, appears questionable.

To show the tendency of a possible extrapolation we plot in Fig. 4 the results for $`x_c`$ for $`S=\frac{1}{2}`$ as a function of $`1/\sqrt{N}`$ and connect by a straight line the results of systems with $`N=12,24`$ spins which have corresponding shapes. The point $`x_c=0.65`$, where this line meets the ordinate shows a reasonable value for $`x_c`$, but obviously cannot be taken as a serious extrapolation. The best bounds we obtained for $`S=1,\frac{3}{2},2`$ are given in Fig. 3.

Let us now consider the ferromagnetic regime. We find, that the dimer state of a four spin plaquette is the ground state for $`x\frac{1}{2S}`$ which implies

$$x_c^f(S)\frac{1}{2S}.$$
(4)
This bound is exact for $`S=\frac{1}{2}`$, since $`x_c^f=1`$ coincides with the boundary of the ferromagnetic phase. For $`S>\frac{1}{2}`$ we again find an improvement on the result of the four spin system by considering larger clusters. Some results for $`S=1,\frac{3}{2},2`$ are shown in Fig. 3.

## III Upper bounds on $`x_c^a`$ and on $`x_c^f`$

Upper bounds for the stability of the dimer phase can be obtained by considering a variational ansatz for finite clusters. We only briefly sketch the idea here, since it is widely used in the literature (see e.g. Ref. and references therein).

The Hamiltonian is split into clusters without common spins and external bonds connecting the clusters,

$$H=\underset{cluster}{}H_{cluster}^N+H_{bond},$$
(5)
as indicated for four site clusters in Fig. 5. We use as variational state $`_{cluster}|\mathrm{\Psi }_{cl}`$, where $`|\mathrm{\Psi }_{cl}`$ is the ground state of $`H_{cluster}^N`$. If the total spin of the clusters is zero in the ground state the expectation value of the external bonds vanishes

$$\underset{cluster}{}\mathrm{\Psi }_{cl}|H_{bond}|\underset{cluster}{}\mathrm{\Psi }_{cl}=0.$$
(6)
We have calculated the ground state energies $`E^N`$ of $`H_{cluster}^N`$ for clusters with up to N=32 spins and obtained upper bounds on $`x_c^a`$ by comparing them with the energy of the dimer state.

For a four spin plaquette and arbitrary spin $`S`$ we thus find

$$x_c^a\frac{S+1}{2S+1}.$$
(7)
For stripe configurations of the type $`2\times \frac{N}{2}`$ the resulting $`x_c^a`$ for $`S=\frac{1}{2}`$ are shown as a function of one over system size in Fig. 6. Note that results for two different types of stripes (as indicated in the figure) are given and that both types extrapolate to almost the same point 0.7126(1) for $`N\mathrm{}`$. A similar extrapolation for stripes of width four (see. Fig. 7) gives an upper bound of $`x_c^a=0.6955(4)`$ for $`S=\frac{1}{2}`$, which is better than the bound obtained for stripes of width two. Again three different shapes converge to about the same result as N approaches infinity. The lowest value in Fig. 7 was obtained for a $`N=32`$ system yielding the exact bound $`x_c^a0.7050`$. An extrapolation of stripes of width two for S=1 and systems of up to 16 spins gives a best upper bound of $`x_c^a(S=1)0.618`$.

In the limit $`S\mathrm{}`$ the upper bound (7) derived from the four spin system does not provide the correct behaviour expected for the classical Heisenberg model. We can however find another bound, which yields a better upper limit for $`S>1`$ by using a helical product state as a variational state.

The ground state energy for $`S=\mathrm{}`$ in the helical phase  is given by

$$E_{hel}=\frac{|J_1|}{2}\frac{J_2^2}{J_1}$$
(8)
where the angle between two neighbouring spins is $`\theta =\pi \pm \mathrm{arccos}(J_2/J_1)`$ and the length of the spins is normalized to 1.

We consider as variational ansatz a system of (quantum) spins polarized along the directions of the classical spins. Comparing the expectation value of this state with the energy of the dimer state we find

$$\left(\frac{|J_1|}{2}+\frac{J_2^2}{J_1}\right)S^2<J_1S(S+1)$$
(9)
from which

$$x_c^a<\frac{1}{\sqrt{2S}}$$
(10)
results as a criterion for the instability of the dimer state. Thus for $`x>1/\sqrt{2S}`$ the dimer state is no longer the lowest state, because the helical product state has a lower energy. Since the ground state energy of the helical state does not depend on the sign of $`J_2`$ Eq. (10) also gives an upper bound on $`x_c^f`$ in the ferromagnetic regime.

## IV Summary

In this paper we have discussed the phase diagram of the Shastry–Sutherland model for arbitrary spin $`S`$. In the regime $`J_1<|J_2|`$ the phase diagram for the quantum models ($`S<\mathrm{}`$) fully agrees with that of the classical model ($`S=\mathrm{}`$). Quantum effects do, however, modify the classical phase diagram in the regime $`0<|J_2|<J_1`$. As emphasized already by Shastry and Sutherland  the dimer phase is strongly stabilized by quantum effects and exists in a regime $`x_c^f(S)<J_2/J_1<x_c^a(S)`$ which grows with decreasing spin S. The nature of the intervening phases originating from the classical helical phases is not totally clear, but their regime shrinks as the dimer phase expands. For $`J_2<0`$ the intervening phase between the dimer and the ferromagnetic phase is certain to exist for $`S>1/2`$, but completely disappears for $`S=1/2`$. For $`J_2>0`$ exact statements on the existence of an intervening phase between the dimer and the Néel phase are not available, but in Ref. rather strong evidence in favor of its existence for $`S=1/2`$ is provided.

The emphasis of the present work was put on the derivation of exact upper and lower bounds on the boundaries $`x_c^f(S)`$ and $`x_c^a(S)`$ of the dimer phase. Our results on this topic are summarized in Fig. 8. The dark grey area shows the region where the dimer state is certain to be the ground state and the white area is the region where the dimer state is certainly not the ground state. From the results presented in this paper it follows, that the exact boundary must be located in the light grey area. The black lines in Fig. 8 represent curves discussed in the text. The upper bounds for $`x_c^a(S)`$ were obtained by a helical variational ansatz ($`\frac{1}{\sqrt{2S}}`$), by a finite cluster ansatz ($`\frac{S+1}{2S+1}`$) and for $`S1`$ by extrapolating results for series of finite clusters (see Figs. 6 and 7). The lower bounds on $`x_c^a(S)`$ were derived from four spin plaquettes ($`\frac{1}{2S+1}`$) and from the finite systems shown in Fig. 3 for $`S2`$. In the ferromagnetic regime we have a lower bound ($`\frac{1}{\sqrt{2S}}`$) from the helical ansatz, which is shown together with the upper bounds obtained from the largest tractable finite clusters (Fig. 3). The intervening phase on the ferromagnetic side extends down to $`x=1`$ (bottom of Fig. 8) where the ferromagnetic regime begins for all $`S`$.

## V Acknowledgments

The authors gratefully appreciate useful discussions with K. Fabricius, G. Uhrig, and J. Zittartz.


format: text

text:
# Integration of GRACE and PYTHIA 11footnote 1Presented by S. Odaka at ACAT2000, Fermilab, Oct. 16 - 20, 2000.

## Abstract

We have successfully developed a technique to integrate an automatic event-generator generation system GRACE and a general-purpose event generator framework PYTHIA. The codes generated by GRACE are embedded in PYTHIA in the created event generator program. The embedded codes give information on parton-level hard interactions directly to PYTHIA. The choice of PDF is controlled by the ordinary parameter setting in PYTHIA. This technique enables us to create easy-to-handle event generators for any processes in hadron collisions. Especially, in virtue of large capability of GRACE, we can easily deal with those processes containing many (four or more) partons in the final state, such as multiple heavy particle productions. This project is being carried out as a collaboration between the Japanese Atlas group and the Minami-Tateya group, aiming at developing event generators for Tevatron and LHC experiments.

GRACE  is a software package for automatic calculation in high energy physics, developed by the Minami-Tateya group. The core part of the package is a program to generate Feynman diagrams relevant to specified initial and final states. It then generates FORTRAN codes for calculating the corresponding cross section on the basis of the amplitude of each diagram. The GRACE package also includes an integration and event-generation program called BASES/SPRING . Hence, it provides us with a very powerful environment for developing event generators for studies of high energy physics.

Since the basic building block is the amplitude, GRACE has an advantage in those processes composed of many coherent diagrams. (The number of elements increases quadratically if the calculation is based on the matrix element.) This turns out to be an advantage in multi-body production processes, since in general the number of contributing diagrams increases as the number of final-state particles increases. The program grc4f is a good example. This is known to be the most reliable event generator at the tree level for four-fermion productions in electron-positron collisions at LEP2 energies.

In future hadron-collision experiments such as LHC, interesting processes (Higgs-boson productions, SUSY-particle productions etc.) in many cases result in multiple production of heavy particles ($`Z/W`$ and/or top quark) or cascade decays. They produce many particles in the final state, and are composed of many coherent diagrams. For example, we have to take into account 144 diagrams in total to evaluate the process $`ppH^0b\overline{b}+Xb\overline{b}b\overline{b}+X`$. Detailed studies of such processes require an exact (coherent) treatment of the whole production and decay reactions. GRACE is expected to be a powerful tool for such studies.

However, GRACE can treat perturbative hard-scattering reactions only. We need to add non-perturbative partonic structures of hadrons (PDF) and QCD evolutions in the initial and final states, in order to construct realistic event generators. The most straightforward way to accomplish it is to connect GRACE to a general purpose event generator, such as PYTHIA , ISAJET  and HERWIG . We chose PYTHIA in this work.

There would be two ways for the connection. One way is to interface them using an external data file, in which hard-scattering event data provided by GRACE are described. Such a system can be flexible and portable. The coding may be easier since only the I/O routines are necessary to be coded. However, since the generation procedure is separated to two steps, a special care or an appropriate software assistance is necessary to keep the consistency in the parameters. This method is used in GRAPE , a previous work of the Minami-Tateya group for electron-proton collisions. A similar method is also adopted in the CompHEP for hadron collisions .

The second way is to embed the GRACE output codes in PYTHIA. Although the coding may be more complicated, the event generation can be a single-step job in this method. Since necessary parameters can be set through the parameter setting in PYTHIA, any inconsistency in the parameter choice can be automatically avoided.

We chose the latter method in the present work. All the program codes describing the hard scattering are embedded in the subroutine PYUPEV, prepared by PYTHIA to install user-defined processes. This subroutine is called during the event generation if it is requested using PYUPIN in the initialization stage.

For the event generation, it is important to give an appropriate definition of ”kinematics”, the mapping of kinematical variables to a set of uniform random numbers. We have tried two methods for it. In the first method, the mapping is fully defined by user-defined analytic functions. The CPU time may be saved if the functions are appropriately defined. The integration/event generation package BASES/SPRING is fully utilized in the second method, where a variable-grid mapping is implemented. Since the grids are optimized by BASES, users need not to care about very details of the mapping. Of course, in both methods, users need to choose an appropriate set of variables to avoid non-diagonal singularities.

The subroutine PYUPEV that we have coded can be separated to two parts, an initialization stage and an event generation stage. The initialization stage is called by PYUPIN once in each job. The maximum of SIGEV (see below) is searched to set the return argument SIGMAX in the functional-mapping method. In the grid-mapping method, BASES is called to estimate the total cross section and SIGMAX is set to be equal to it. The grids are optimized here.

PYTHIA calls the second stage in the event-generation loop with a frequency corresponding to SIGMAX. In this stage, first of all, one of the ”sub-processes” is chosen. In most cases in hadron collisions, every ”process” of interest is composed of several incoherent sub-processes, since hadrons are composite and parton species in the final state are hard to identify.

The following procedure is different for two different methods of ”kinematics”. In the functional-mapping method, a set of kinematical variables defining an event is determined from a set of uniform random numbers according to the defined ”kinematics”. The differential cross section of this event is calculated using the GRACE output codes. The return argument SIGEV is calculated from the differential cross section, PDF and the Jacobian of the ”kinematics”.

Instead, the event generator SPRING is called in the grid-mapping method. SPRING generates an event using the GRACE output codes, PDF and the grid information optimized by BASES. SIGEV is always equal to SIGMAX.

In both methods, the above procedure is followed by a Lorentz boost to the laboratory frame and a determination of the color flow. It should be noted that the color flow can be determined automatically using the information from BASES.

After returning to PYTHIA, an event sampling is done according to the weight SIGEV/SIGMAX. The sampling is dummy (i.e., all events are accepted) in the grid-mapping method. After that, the initial- and final-state parton radiations are added by PYTHIA to simulate the effect of QCD evolution, resulting in an underlying hadronic activity and a finite transverse momentum of the hard-scattering system.

The coding has been refined by developing actual event generators. So far, we have developed generators for the following processes:

$$pp(\mathrm{or}p\overline{p})q\gamma +X,$$
(1)
$$pp(\mathrm{or}p\overline{p})Wg+X\mu \nu g+X,$$
(2)
$$pp(\mathrm{or}p\overline{p})H^0W+Xb\overline{b}\mu \nu +X,$$
(3)
$$pp(\mathrm{or}p\overline{p})b\overline{b}(\mathrm{QCD})W+Xb\overline{b}\mu \nu +X,$$
(4)
$$pp(\mathrm{or}p\overline{p})H^0b\overline{b}+Xb\overline{b}b\overline{b}+X,$$
(5)
$$pp(\mathrm{or}p\overline{p})b\overline{b}b\overline{b}(\mathrm{QCD})+X.$$
(6)
The results were compared with those from other existing generators, PYTHIA built-in generators and CompHEP, to examine the coding. We found reasonable agreements in all cases.

As an example, the performance of a developed event generator is compared with that of a PYTHIA built-in generator in Table 1 for one of the simplest cases, Process (2). In this study, the weak boson was not made to decay and the simulations of both initial- and final-state parton radiations were turned off. The two ”kinematics” methods were tried.

The cross section is in very good agreement, showing that the integration of GRACE and PYTHIA is done successfully. The functional-mapping method gives a better performance in CPU time than the grid-mapping method, as expected. Although both GRACE + PYTHIA generators consume appreciably longer CPU time than the PYTHIA built-in generator, the difference is not very serious since the simulation of parton radiations and hadronization, which are not implemented in this study and common to all three generators, takes much longer time.

The most noticeable advantage of the GRACE + PYTHIA system is in the fact that using this technique we can easily develop event generators for those processes which are not and/or hard to be implemented in PYTHIA. Process (5) is one of such processes. The total cross section estimated by using the developed generator is shown in Table 2 for three assumed Higgs-boson masses, and compared with the result of CompHEP. We have applied the grid-mapping method only in this case. The agreement with CompHEP is quite good.

We plan to make some improvements in order to make the development easier. So far, routines interfacing PYTHIA and GRACE are written by hand. Since the functionality of these routines are almost common to all generators, we will be able to make them automatically generated by changing the libraries referred by GRACE.

Processes in hadron collisions are in most cases composed of several incoherent sub-processes. The present GRACE system cannot handle such processes automatically. Some modifications and additions by hand are necessary now. We would like to automate these tasks.

The difference between the sub-processes is, in most cases, the difference in quark species in the initial and/or final states. If we can treat quark masses and their couplings as variables, these sub-processes can share an identical code for the calculation. The Cabbibo-Kobayashi-Masukawa matrix can be implemented automatically if such a treatment is realized.

In summary, we have established a technique to embed GRACE output codes in PYTHIA. This technique allows us to develop new hadron-collision event generators easily. We have applied the technique to some two-, three- and four-body production processes. Obtained results are in good agreement with the results from existing generators. The generator for four b-quark productions, Processes (5) and (6), will be released in a few months. Since GRACE is expected to be advantageous in multi-body production processes, we would like to go to five- or six-body production processes as the next step. We also have an automatic next-to-leading order (NLO) calculations for hadron-collision processes in our view.

This project is being carried out as a collaboration between the ATLAS-Japan group composed of Japanese members of the ATLAS collaboration at LHC, and the Minami-Tateya group, aiming at developing event generators for Tevatron and LHC. K.S., S.T. and S.O. are from the ATLAS-Japan group, and J.F., T.I. and Y.K. from the Minami-Tateya group. There are many other people from these two groups who have contributed to this work. Among them, the authors wish to thank here Y. Takaiwa, S. Kawabata and K. Kato for their educational contributions, and T. Abe for discussions.


text:
# Modeling the spectral energy distribution of galaxies. II. Disk opacity and star formation in 5 edge-on spirals

## 1 Introduction

Optical and Far-Infrared (FIR)/sub-millimeter (submm) data from galaxies contain complementary information about the distribution of stars and dust, from which intrinsic quantities of interest - the star-formation rate (SFR) and star-formation history - can in principle be extracted. This is especially relevant to systems having intermediate optical depths to starlight, or to inhomogeneous systems with optically thin and thick components. Many, perhaps most, of normal (non-starburst) gas-rich galaxies in the local universe may fall into these categories.

On the one hand, optical data probes the colour and spatial distribution (after correction for extinction) of the photospheric emission along sufficiently transparent lines of sight. This is particularly useful to investigate older, redder stellar populations in galaxian disks with scale heights larger than that of the dust. A radiation transfer modeling technique that can be applied to edge-on systems, where the scale height of the stars and dust extinction can be directly constrained, was introduced by Kylafis & Bahcall (kylafis (1987)) and subsequently applied on several galaxies (Xilouris et al. xilouris1 (1997); xilouris2 (1998); xilouris3 (1999), hereafter referred to as X97, X98, X99 respectively). Radiative transfer codes in combination with observations of nearby edge-on galaxies were also used by Ohta & Kodaira (ohta (1995)) and Kuchinski et al. (kuchinski (1998)).

On the other hand, grains act as test particles probing the strength and colour of ultraviolet (UV)/optical interstellar radiation fields. This constitutes an entirely complementary constraint to studies of photospheric emission. In the FIR, grains are moreover detectable over the full range of optical depths present in a galaxy. At least part of this regime is inaccessible to direct probes of starlight, especially at shorter wavelengths, even for face-on systems. This particularly applies to light from young stars located in, or close by, the dust clouds from which they formed, since a certain fraction of the light is locally absorbed. Furthermore, there is at least a possibility that most of the remaining UV and even blue light from young stars that can escape into the disk might be absorbed by diffuse dust there. Observations from IRAS, IRAM (e.g. Neininger et al. neininger (1996); Dumke et al. dumke (1997)), SCUBA (e.g. Alton et al. alton (1998); Israel et al. israel (1999); Dunne dunne (2000); Bianchi et al. 2000a ) and ISO (e.g. Haas et al. haas (1998); Stickel et al. stickel (2000)) can provide information on the quantity and spatial distribution of the dust within spiral galaxies.

A combined analysis of the whole UV-optical/FIR/submm/mm spectral energy distribution (SED) of galaxies seems to be a promising way to constrain the problem. Radiative transfer codes for an assumed “sandwich” configuration of dust and stars were applied by Xu & Buat (xu1 (1995)) and Xu & Helou (xu2 (1996)) to account for the energy budget over the whole spectral range. They considered in detail the relative contribution of the non-ionizing UV photons and the optical photons in heating the grains. However these calculations did not incorporate a model for the dust grain emission, nor for radial variations in the absorbed radiation in the disk, and therefore could not account for the exact shape of the FIR SED. Recently, there have been several works modeling the SED of galaxies from UV to submm (Silva et al. silva (1998); Devriendt et al. 1999), by applying photometric and/or spectrophotometric and chemical evolution models of galaxies. While these models are adequate in describing the general shape of the volume-integrated SED, they make use of many free parameters and of a simplified geometry.

Bianchi et al. (2000b ) attempted to model NGC 6946 from the UV to FIR using a 3D Monte Carlo radiative-transfer code for a simplified geometry of emitters (a single stellar disk). They concluded that the total FIR output is consistent with an optically thick solution. However, their model did not consider a grain size distribution for grains, stochastic emission of small grains, and the contribution of localized sources within star-forming complexes. This resulted in a poor fit of the FIR SED and a failure to reproduce the IRAS flux densities.

Another work on modeling the UV to submm emission was presented by Efstathiou et al. (efstathiou (2000)) for star burst galaxies which were treated as an ensemble of optically thick giant molecular clouds centrally illuminated by recently formed stars. This modeling technique successfully reproduced the observed SED of M82 and NGC 6090. Such a technique obviously cannot be applied to normal “quiescent” disk galaxies dominated by emission from the diffuse interstellar radiation field.

In Popescu et al. (popescu (2000), hereafter referred to as Paper I) we gave a detailed description of a new tool for a combined analysis of the optical-FIR/submm SED of galaxies. This tool included solving the radiative-transfer problem for a realistic distribution of absorbers and emitters and by considering realistic models for dust, taking into account the grain size distribution and stochastic heating of small grains as well as the contribution of HII regions. We applied this tool to the edge-on system NGC 891.

In short, we used the intrinsic spatial distribution of dust and stars (in the B, V, I, J, K bands) derived for this galaxy by X99 as a constraint on the old stellar population and part of the dust distribution. In order to fully explain the optical-FIR-submm SED, it was found necessary to add to the components of the X99’s solution both a young stellar population (to correct for a shortfall in FIR brightness in terms of re-radiated non-ionizing UV) and more dust mass (to correct for a shortfall in submm brightness). Our solution explained both the observed SED in the FIR and submm, as well as the observed radial profile at 850 $`\mu `$m. We found that the dust is predominanly heated by the young stellar population.

Although NGC 891 was chosen for study in Paper I as the prototypical and most extensively observed edge-on galaxy, it exhibits some extreme features, raising the possibility that it might not in fact be representative of a “typical” spiral. In this paper we extend our SED modeling technique to four additional edge-on systems - NGC 5907, NGC 4013, UGC 1082 and UGC 2048 - with the aim of examining whether the features of the solution we obtained for NGC 891 might be more generally applicable. Radiation transfer solutions accounting for the appearance of these four galaxies in the optical bands have been obtained by X99, X98 and X97. We use these solutions to constrain our model for the optical-FIR SEDs in the same fashion as for NGC 891.

The paper is arranged as follows: In Sect. 2 we overview the model and its application to the objects in the present work. In Sect. 3 we present predictions for the FIR SED of all five galaxies studied. These are obtained using our “standard model”, i.e. a model with radiation fields derived from the radiation transfer solution for the optical appearance supplemented by a hidden population of young stars embedded in HII regions, but with no additional dust. In Sect. 4 we describe and discuss the SED of one of our four galaxies, namely NGC 5907, in terms of the more complex “two-dust-disk model” and compare this with the corresponding solution obtained for NGC 891 in Paper I. In Sect. 5 we discuss the dependence of the derived star-formation rates on the assumptions of the model, in particular the assumed geometries for the dust and stars. We show that after normalization to disk area, NGC 891 has comparable star-formation characteristics to the other four objects. We also compare the disk-area normalized massive star-formation rate in the five edge-on galaxies, obtained from the optical-FIR SED analysis, with the same quantity derived from Hα measurements of statistical samples of face-on galaxies given in the literature. In Sect. 6 we give a summary of our work.

## 2 Model calculations for the five galaxies

A full description of our model and its application to NGC 891 can be found in Paper I. Here we give an overview of the basic idea and its specific application to the four galaxies.

To determine the propagation of light in a galaxy we take both emitters and absorbers to be smoothly arranged as a superposition of simple cylindrically symmetric geometrical distributions. The assumption of an axisymmetric model and the neglect of the spiral structure does not seem to introduce a systematic error on the estimate of the overall opacity, at least when a statistical sample of edge-on galaxies is considered (Misiriotis et al. misiriotis (2000)). The emitters are divided into a young (predominantly UV-emitting) stellar population, distributed in a disk with a small scale height, an old (predominantly optical-NIR emitting) population, distributed in a disk with a larger scale height and a bulge. The absorbers are arranged in “young” and “old dust disks”, associated with the old and young stellar populations. The “old dust disk” has a larger scale height than the “young dust disk”, though not so large as the old stellar disk. Both the stellar- and dust-disk distributions are specified as exponential density distributions in radius and height, while the stellar bulge is described by a de Vaucouleurs law. All radial distributions have a common truncation at three scale lengths of the “old dust disk”.

The first step in the calculation of the optical-FIR SED is to determine the geometric and amplitude parameters for the old stellar disk and bulge independently in each band pass ranging from the B band to the NIR for a common “old dust disk” (whose geometry and opacity is also determined). This can be achieved by a radiation transfer analysis matching as closely as possible the predicted and observed brightness distributions, following Kylafis & Bahcall (kylafis (1987)) and including inclination as a free parameter. For the five galaxies discussed in this work we adopt the solutions given in X97 and X99 on the basis of observations in B, V and I, and in the case of NGC 891 also in J and K bands.

This analysis of the optical-NIR emission also yields the extinction coefficient as a function of wavelength for the “old dust disk”. For all five galaxies studied, the measured extinction law was found to be close to the predicted law derived from the dust model used in the calculation of the IR emission. This corresponds to the graphite-silicate mix and dust size distribution determined for the Galaxy by Mathis et al. (mathis (1977)) (see Eqs.  and  of Paper I) with optical constants taken from Laor & Draine (laor (1993)). When calculating the IR emission we take grains in the “young dust disk”, if present, to have the same size distribution and optical properties as those in the “old dust disk”.

To this fixed basis for the old stellar population and associated dust disk, the model allows the addition of a young, UV-emitting stellar population and, optionally, associated dust. The additional stars and dust are taken to be distributed in a common exponential disk with a scalelength equal to that of the already determined intrinsic B-band population, and a scale height fixed at 90 pc (as for the Milky Way). Since this scale height is typically a factor of several times smaller than the scale heights derived by X97 and X99 for the “old dust disk”, this constitutes the simplest assumption for the distribution of these extra components, which hide UV-optical indicators of the young stellar population and associated additional opacity. As emphasized in Paper I, a clumpy distribution for the dust associated with the young stellar population is also possible, though very difficult to calculate, and, due to the lack of high angular resolution FIR images of spiral galaxies, at present not directly verifiable. The luminosity of the young stars is a free parameter of the model, which we express as a recent star-formation rate SFR, based on the population synthesis models of Bruzual & Charlot (bruzual (2001)) for $`Z=Z_{}`$, a Salpeter IMF with a cut-off mass of $`100\mathrm{M}_{}`$ and $`SFR=(1/\tau )\mathrm{exp}(t/\tau )`$ with $`\tau =5`$Gyr.

The model also incorporates a parameter $`F`$ representing the fraction of the non-ionizing UV luminosity from the young stellar population, which is locally absorbed in star-formation regions by optically thick fragments of their parent molecular clouds. This fraction of the luminosity is considered to be re-radiated in the FIR according to a spectral template matching the observed spectrum of galactic HII regions. This template generally has a warmer 60/100 $`\mu `$m colour than that calculated for the diffuse disk.

Given the intrinsic distributions of stellar emissivities in the UV, optical and NIR, and the distribution of absorbers in the “old” and (if present) the “young dust disks”, we then proceed with a radiative transfer calculation to determine the UV-NIR energy density of the radiation field throughout the galaxy. In the absence of solutions in the J and K bands, we extrapolate the energy density from the B, V, I bands assuming a black-body spectral distribution and a colour temperature of 4000 K. In calculating the heating of grains placed in the resulting radiation field, we include an explicit treatment of stochastic emission as detailed in Paper I. The IR-submm emission from grains is then calculated for a grid of positions in the galaxy. Subsequently we integrate over the entire galaxy to obtain the IR-submm SED of the diffuse disk emission. Prior to comparison with observed FIR-submm SEDs, the spectral template for the HII regions, scaled according to the value of $`F`$, must be added to this calculated spectral distribution of diffuse FIR emission.

Due to the precise constraints on the distribution of stellar emissivity in the optical-near infrared (NIR) and the distribution of dust in the “old dust disk” yielded by the radiation transfer analysis of the highly resolved optical-NIR images, coupled with the simple assumptions for the distribution of the young stellar population and associated dust, our model has a maximum of three free parameters from which the FIR-submm SED can be fully determined. These are the SFR for massive stars, the parameter $`F`$ representing the fraction of the non-ionizing UV luminosity from massive stars which is locally absorbed, and the mass of dust in the “young dust disk”. In general terms, the last parameter is principally constrained by the submm emission, the non-ionizing UV luminosity by the bolometric FIR-submm output and the factor $`F`$ (in the absence of high resolution images) by the FIR colour. In the analysis of the five galaxies in subsequent sections, we refer to the following variants of the model:

1. The simple model: The simplest possible model takes into account only the information from the optical observations, i.e. the dust and old stellar distribution in the old disk without massive UV-emitting stars. The sole purpose of the “simple model” is to demonstrate that the old stellar population can not adequately heat the dust to account for the observed FIR luminosity.
2. The standard model: This is the “simple model” supplemented by a disk population of UV-emitting stars without associated dust, but with localised FIR emission. The free parameters are the SFR and the factor $`F`$.
3. The two-dust-disk model: This is the “standard model” supplemented by the presence of dust in the young stellar disk. The free parameters are the SFR, the factor $`F`$, and the dust mass in the “young dust disk”.

The predictions of the model calculations are compared with available FIR, submm and mm data from the literature. In Table 1 we summarise the data used for comparison in the present work. We adopted $`20\%`$ errors in the IRAS and IRAM data. For the SCUBA data we adopted the errors quoted by Alton et al. (alton (1998)). In Table 1 we also give some basic information on the properties of the sample galaxies. One factor in the choice of these galaxies for the radiation transfer analysis by X97 and X99 was their prominent dust lanes.

## 3 Results for the “standard model”

To demonstrate the need for a population of massive stars to account for the observed FIR luminosities in all five galaxies, we first calculated SEDs for the “simple model”, in which the stellar and dust distributions are completely determined from the optical-NIR data as described in Sect. 2. These are plotted in the left column of Fig. 1 In all cases, the predicted FIR flux densities in the 60 and 100 $`\mu `$m IRAS bands fall short of the observed luminosities by factors of between 5 and 10, showing that the heating of the dust from the old stellar population cannot account for the FIR emission. The results of the “simple model” are summarised in the second column of Table 2 ($`L_{\mathrm{simple}}`$), where we give the FIR luminosity of the galaxies that can be attributed to the old stellar population.

Assuming that the optically determined dust content we used for each galaxy is exact, the “standard model” can account for the observed FIR luminosity by transforming non-ionizing UV radiation from massive stars into the FIR through a combination of local absorption and absorption in the diffuse old disk, controlled by the factor $`F`$. If all the young stellar luminosity is absorbed locally ($`F=1`$), then the efficiency of the transformation from non-ionizing UV to the FIR is unity, as for starburst galaxies. In the total absence of local absorption ($`F=0`$) the global absorption efficiency is lower, with a value depending on the relative geometries of the diffuse dust and the massive stars. In principle, the factor $`F`$ can be determined from the FIR colours.

We investigated the combinations of $`F`$ and SFR for which the “standard model” was consistent with the IRAS 60 and 100 $`\mu `$m measurements. For NGC 891 and UGC 2048, good agreement between model and data could be found for values of $`F`$ well above zero. However, for the remaining galaxies, the IRAS colours required $`F`$ values of zero or close to zero. This would imply high current star-formation rates and an amplitude in the UV stellar emissivity which would not join up smoothly with the empirically determined B band stellar emissivities. Consequently, we have plotted in Fig. 1 the model prediction corresponding to the maximum value of $`F`$ still consistent with the IRAS colours. This also gives a monotonically increasing SED for the UV-optical stellar emissivities. Corresponding FIR bolometric luminosities ($`L_{\mathrm{standard}}`$), SFRs and factors $`F`$ are given in Table 2, columns 3, 4, 5 respectively. We emphasise however, that in principle one could have a non-steady SFR, so that the tabulated SFR are at the same time lower limits in SFR and upper limits for the factor $`F`$. Therefore we can discuss the results for the “standard model” also in terms of a range of possible values for $`F`$ and SFR, which will be compatible with the range in the IRAS errors.

We remark that a comparison of the values in the second and third column of Table 2 reveals that, according to the standard model, only 10% to 30% of the dust emission is powered by the old stellar population.

In the above we have made the assumption that the “standard model” is a sufficiently realistic model of the galaxies. Some evidence that this may not be the case is apparent from Fig. 1 (right column), where all the model spectra peak at around 120 $`\mu `$m, well shortward of the typical peak near 150 $`\mu `$m seen in spiral galaxies. This may be an indication that the overall optical depth of the galaxies is underestimated by the “standard model”, indicating the presence of a further dust component in the plane of each galaxy. This component may either be smoothly distributed or be in the form of optically thick clumps or some combination of these. For the two galaxies that submm and/or mm measurements are available (NGC 891 and NGC 5907), this is explicitly confirmed by the undershoot of the “standard model” prediction in the submm and mm regime. In Paper I, we used the available submm fluxes to quantify the missing dust and introduced the “two-dust-disk model”, which successfully reproduced the spectrum and the spatial distribution of the FIR emission in NGC 891. In the next section we will apply the “two-dust-disk model” to NGC 5907 to determine the total dust content of this galaxy and examine the effect of this more realistic case on the derived SFRs.

## 4 Application of the “two-dust-disk model” to NGC 5907

In Paper I we proved that in NGC 891 the amount of dust (on the basis of the grain optical constants adopted) quoted by X99 is underestimated by a factor of two, and that the addition of a second dust disk, that follows the spatial distribution of the young stellar population, reproduces very well the observed FIR luminosity and the radial profiles observed by Alton et al. (alton (1998)). For NGC 5907, the measurement at 1200 $`\mu `$m by Dumke et al. (dumke (1997)) can be used to constrain the quantity of dust in the second disk.

Inspection of Fig. 1 (right column) shows that the “standard model” solution for $`F=0.08`$, which is in agreement with the IRAS points, crassly undershoots the 1200 $`\mu `$m point. As for NGC 891, the only way to account for both the IRAS and the 1200 $`\mu `$m observations is to include more dust in our model, as simply increasing the star-formation rate would provide too much luminosity in the IRAS range.

The second-dust-disk’s spatial distribution was constrained to follow that of the young stellar population, which is taken to be exponential with scalelength and scaleheight $`h_\mathrm{d}=5020`$ pc and $`z_\mathrm{d}=90`$ pc respectively (see Sect. 2). To find the solution, the total amount of dust in the second dust disk was varied jointly with the SFR until the FIR-submm SED of our model fitted all three measurements. In Fig. 2 we present the resulting spectrum overlaid with the data. This spectrum corresponds to a SFR$`=2.2\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$ and $`F=0.10`$. The total dust mass in the second dust disk is $`4.5\times 10^7\mathrm{M}_{\mathrm{}}`$. By comparison, the dust mass quoted for the old disk of NGC 5907 by X99 is $`1.5\times 10^7\mathrm{M}_{\mathrm{}}`$.

Thus, we have had to add 3 times more dust over the quantity implied by the observed dust lane to obtain agreement with the data on the basis of the “two-dust-disk model”. The derived value of $`F`$ is still very low, but the spectrum of the diffuse component is too warm (NGC 5907 is more optically thin than NGC 891) to allow the addition of a warm component related to HII regions. While we could have produced a colder spectrum for the diffuse emission, allowing some room for HII regions, this would have been at the cost of adding yet more dust in the second disk. A larger quantity of dust in the second disk does not seem very likely, as the total dust mass for the fitted $`F=0.10`$ case already implies a gas-to-dust ratio of 130 by mass (instead of 520 according to X99). As a comparison, the gas-to-dust ratio in the Milky Way is 133.

The total IR-submm re-radiated luminosity of NGC 5907, obtained by integrating the $`F=0.10`$ “two-dust-disk model” SED, is $`50.5\times 10^{35}`$ W out of which $`27.0\times 10^{35}`$ W is attributed to heating from the young stellar population. Thus, about $`40\%`$ of the dust emission is powered by the old stellar population. More dust in the “young dust disk”, a lower SFR and a higher factor $`F`$ would increase the disk opacity and thus the contribution of the old stellar population, but we excluded this solution due to the high dust-to-gas ratio. The $`40\%`$ contribution from the old stellar population represents an upper limit, since a possible clumpiness in the dust distribution would decrease the absorption efficiency in the optical band.

On the basis of the model’s inherent assumption, that all the dust (in both disks) of the galaxy is distributed smoothly, the central face-on optical depth in the optical band is $`\tau _\mathrm{v}=1.4`$ for NGC 5907.

The major difference between NGC 891 and NGC 5907, on the basis of the “two-dust-disk model”, is that the spectrum of the former apparently allows for the existence of a larger contribution from HII regions. Presumably NGC 5907 must have HII regions, and there is no obvious reason why the local properties of these HII regions, in particular the probability of local absorption of the UV radiation, should differ between the two systems.

One potential way out of this difficulty might be to suppose that the spectral template for the HII regions used in our model underestimates the associated submm emission. As stated in Paper I, we did not attempt to include potential cold dust emission components that might be expected from “parent” molecular clouds in juxtaposition to their “offspring” HII regions. These might contribute a significant fraction of the observed submm emission, but be relatively faint in the IRAS range. If so, the dust content and the luminosity of the diffuse disk would be lower for both galaxies, but more so for the more optically thin diffuse disk of NGC 5907. This would make room for larger values of $`F`$ in the model fits, particularly for NGC 5907.

An alternative scenario, also discussed in Paper I, would be the presence of truly quiescent cold dust clouds in the disk not associated with HII regions. In that case, the difference between the two galaxies might indicate a greater star-formation activity in the disk of NGC 891, in the sense that a larger proportion of the cold dust clouds in the disk had been triggered into massive star formation.

## 5 Discussion

Our results confirm that it is a common tendency for edge-on galaxies to hide a significant fraction of their dust from optical extinction studies. In two cases where sub- or near-mm data are available (NGC 891, as analyzed in Paper I and NGC 5907 in this work), we have modeled the FIR-submm SED assuming that the additional dust is distributed in a second diffuse disk associated with the young stellar population. We have found that the SEDs could be fitted in terms of total dust masses of respectively twice and four times the masses inferred by X99 from the optical radiation transfer analysis for NGC 891 and NGC 5907. The corresponding SFRs and local non-ionizing UV absorption factors $`F`$ in the best fits were SFR $`=3.8`$ M$`{}_{}{}^{}\mathrm{yr}_{}^{1}`$ and $`F=0.22`$ for NGC 891 and SFR $`=2.2`$ M$`{}_{}{}^{}\mathrm{yr}_{}^{1}`$ and $`F=0.10`$ for NGC 5907.

The main uncertainty in the inferred SFRs is the value of the factor $`F`$. In principle, $`F`$ can be determined from the colour of the FIR SED, especially the 60/100 $`\mu `$m colour ratio. However, as discussed in Sect. 4, the value of the factor $`F`$, which can be fitted, in practice is constrained by the assumption of the “two-dust-disk model” that all the additional dust needed to fit the submm measurements is diffuse. In reality, the geometry of the additional dust could be, as emphasised in Paper I, in the form of optically thick clumps associated with the HII regions. In this case the inferred SFRs would, depending on the factor $`F`$, differ from those determined on the assumption of a purely diffuse dust distribution without local absorption. However, at present, the available data does not allow us to quantify the clumpy component of UV-absorbing dust.

The SFRs and factors $`F`$ inferred from the “standard model” for NGC 891 and NGC 5907 in Sect. 3 (where the dust content and distribution is determined from X99’s optical analysis) hardly differs from the values determined from the “two-dust-disk model” applied to NGC 891 and NGC 5907 in Sect. 4. This is at first sight surprising, as one would expect lower SFRs due to the higher opacity of the diffuse disk in the “two-dust-disk model”. However, the SFR determined from the “standard model” was obtained by fitting the SEDs exclusively to the IRAS 60 and 100 $`\mu `$m measurements common to the sample of the 5 galaxies. The reason that the SFRs determined from the “standard model” are so similar to those derived from the “two-dust-disk model” is simply that the additional submm luminosity arises from the extra dust (from the second dust disk) rather than from additional SFR. As remarked in Sect. 3, the rather warm SEDs for all the “standard model” fits suggest that additional dust is present in the disks of all 5 galaxies.

To statistically evaluate our results for SFR, we make the assumption that the SFRs determined from the “standard model” fits to the IRAS 60 and 100 $`\mu `$m measurements for our sample of 5 galaxies give, as is the case of NGC 891 and NCG 5907, a reasonable estimate of the SFR which would have been obtained with an analysis based on the “two-dust-disk model”. We can then compare the SFR characteristics of the 5 edge-on galaxies with the larger sample of 61 galaxies with inclinations less than 75 degrees which were studied by Kennicutt (kennicutt (1998)) on the basis of Hα measurements.

Fig. 3 depicts the relation between the disk-averaged surface density in SFR ($`\mathrm{\Sigma }_{\mathrm{SFR}}`$) determined from the “standard model” as a function of the average gas surface density ($`\mathrm{\Sigma }_\mathrm{g}`$) for the 5 galaxies. The plotted data is summarized in Table 3. The $`\mathrm{\Sigma }_{\mathrm{SFR}}`$ is calculated from the SFRs determined from the “standard model” and listed in Table 2. The full range in conceivable SFRs is given by the vertical error bars. The upper and lower limits for the SFRs are calculated such that the predicted SED is still consistent with the IRAS colours (within the 20$`\%`$ IRAS error bars). Lower error bars are not given when the plotted SFRs represent lower limits for the SFR (maximum values for the factors $`F`$, see discussion in Sect. 3). In these cases lower limits would be possible only if we allowed for different sources of uncertainties, like variations in the spectral shape of the template used for the HII regions. However this is hard to quantify, and in the following we assume that the errors of the SFR are given only by the uncertainties in the IRAS data.

The $`\mathrm{\Sigma }_\mathrm{g}`$ were calculated from the gas masses taken from García-Burillo et al. (garcia (1992)) and Rupen (rupen (1991)) respectively for H2 and HI in NGC 891, from Dumke et al. (dumke (1997)) for H2 and HI in NGC 5907, from Bottema (bottema (1995)) and Gomez & Garcia (gomez (1997)) for H2 and HI in NGC 4013, from Giovanelli & Haynes (giovanelli (1993)) for HI in UGC 1082 and from Huchtmeier & Richter (huchtmeier (1989)) for HI UGC 2048. In the last two cases, where only data for HI was available, we estimated the total gas mass by doubling the HI gas mass. The horizontal error bar corresponds to the uncertainty in the gas masses for which we adopted an average 0.2 dex error. The surface area of the disk was calculated for $`R_\mathrm{o}=(3\pm 0.5)h_\mathrm{d}`$, where $`h_\mathrm{d}`$ is the intrinsic radial scalelength determined from the radiation-transfer modeling in the I band. In their analysis of surface photometry of the outer regions of spiral disks, Pohlen et al. (pohlen (2000)) show that the disk boundaries are typically in this range.

The points for the 5 galaxies in Fig. 3 lie within the area of the diagram occupied by the galaxies in the Kennicutt sample. The match is even better for those members of the Kennicutt sample with Hubble types Sb to Sc. This agreement is quite reassuring, bearing in mind the several factors which could introduce a systematic difference between the SFRs inferred for a sample of nearly face-on systems from Hα measurements compared with the present technique for edge-on systems based on an analysis of broad-band non-ionizing UV re-radiated in the FIR-submm range.

Firstly, the Hα analysis is sensitive to the most massive stars and in particular to the assumed mass cut of 100 M. Whereas the FIR-submm modeling also assumes the same mass cut in the conversion of SFR to non-ionizing luminosity (see Sect. 2 and Paper I), our model is less sensitive to this effect.

Secondly, whereas the Hα is sensitive to the star-formation history of the last 107yr, our broad-band FIR-submm SED analysis samples approximately the last 108yr. Thus, our analysis is consistent with the basic hypothesis (see Kennicutt kennicutt (1998)) for “normal” spiral galaxies of a steady star-formation activity. In principle, we could extend our analysis based on our determinations of the intrinsic populations and use the determined intrinsic colours to determine more accurately the SFR history of the galaxies, though this is beyond the scope of this paper.

The assumption of a steady-state star-formation rate is also broadly consistent with the timescales for the exhaustion of the current gas supply under the derived SFRs. The dotted, dashed and dot-dashed lines in Fig. 3 represent star-formation efficiencies corresponding to consumptions of 100, 10 and 1 percent of the gas mass in 108 yr. We note that in constructing Fig. 3, we have assumed that the HI line emission, from which the gas masses for the 5 edge-on galaxies were in part determined, have not been affected by self absorption in the 21cm line. The effect of self-absorption would be to shift the points for the 5 edge-on galaxies to the left of the positions corresponding to their actual gas surface densities in Fig. 3.

Thirdly, the SFRs derived from Hα were corrected by a single factor for extinction, despite the varying orientations. As well as possibly affecting the vertical position of the galaxies on the plot, this may induce some scatter, especially if all the dust were diffusely distributed. The systematic effect may be expressed in terms of the factor $`F`$: an overestimation of the factor $`F`$ is equivalent to an overestimation of the local extinction in the star-formation regions (statistically averaged over the population of HII regions in a disk). Thus, while moving to higher factors $`F`$ would move the points for the 5 galaxies towards lower SFRs, it would have the opposite effect for the SFRs determined from the Hα.

Lastly, we remark that NGC 891 does not appear as an exceptional system compared with the other 4 galaxies in our sample (and with Kennicutt’s  normal galaxy sample) on the basis of SFR normalized to disk area. Our work thus provides no evidence that this galaxy’s exceptional layer of extraplanar Hα-emitting diffuse ionizing gas (e.g. Hoopes et al. hoopes (1999)) and surrounding X-ray-emitting hot gas (e.g. Bregman & Houck bregman (1997)) is attributable to unusual star-formation activity.

## 6 Summary

* This paper is the second part on a series of papers dedicated to the modeling of the SEDs of disk galaxies. In Paper I we described a new tool for the analysis of the UV to the submm SED and applied this tool to the well known nearby edge-on galaxy NGC 891. In the present paper we have extended the analysis from Paper I to a sample of four additional edge-on galaxies (NGC 5907, NGC 4013, UGC 1082, UGC 2048) and proved that the solution obtained for NGC 891 is generally valid. We have also shown that NGC 891 is not an exceptional system in terms of its SFR density.
* In all the systems the dust is predominantly heated by the young stellar population. The contribution of the old stellar population can account for at most $`40\%`$.
* The SFRs derived from our “standard model” are consistent with the SFRs derived from the more sophisticated “two-dust-disk model” and are also consistent with the Kennicutt’s (kennicutt (1998)) Schmidt law for disk galaxies.

###### Acknowledgements.

One of us (A.M.) thanks the Max Planck Institut für Kernphysik for financial support and hospitality during his visit to Heidelberg.


text:
# Collective Modes in a Dilute Bose-Fermi Mixture

## Abstract

We here study the collective excitations of a dilute spin-polarized Bose-Fermi mixture at zero temperature, considering in particular the features arising from the interaction between the two species. We show that a propagating zero-sound mode is possible for the fermions even when they do not interact among themselves.

PACS numbers: 03.75.Fi, 67.60.-g, 67.55.Jd

Recent experimental progress in atomically trapped gases has led to a resurgence of interest in quantum fluids. A particular notable feature is the number of systems available, ranging from single component bose gas in the original experiments where Bose-Einstein Condensation (BEC) was first achieved  to binary bose mixture , spinor condensate in optical traps  and degenerate fermi gas . Other systems have also received much recent attention, in particular bose-fermi mixture. This last mentioned system occurs naturally if ‘sympathetic cooling’ is employed to reduce the kinetic energy of the fermions . There have already been several studies on the properties of this system. Questions addressed include stability against phase separation  and collective excitations .

Although Bose-Fermi mixtures have been studied intensively in low temperature physics in the context of 3He-4He mixtures , the atomically trapped gases offer many additional possibilities. By the choice of atoms, concentration of the various components, or the control of interaction strength among them by external fields , one can unmask phenomena previously unobservable. In this paper, we shall study one example of this by considering the density oscillations of a bose-fermi mixture at low temperatures. We shall show that a variety of novel phenomena can arise due to the coupling between the two components for suitable parameters such as the ratio of the sound velocity of the bose gas to the fermi velocity of the fermions. In particular, we shall show that it is possible to have a propagating fermionic sound mode even in the absence of interaction among the fermions themselves. Sound propagation has also been considered in ref  which however did not investigate the effects being studied here. We shall comment on this later.

We shall then consider a mixture of weakly interacting bose and fermi gases at zero temperature. Both gases are assumed to be spin-polarized such as would be the case usually in magnetic traps. For a dilute mixture, interaction among the bosons themselves and between the bosons and fermions can be characterised by the scattering lengths $`a_{bb}`$ and $`a_{bf}`$ in the s-wave channels. The fermions however do not interact among themselves since they are spin-polarized. For simplicity we shall consider a uniform system. We shall further assume that the gas is stable against phase-separation unless explicitly specified. We are interested in the density waves of this system. As we shall see in general the modes may be damped. Also since the density oscillations are likely to be studied by exciting the systems with external potentials, we shall instead consider the density responses of the system under external perturbing potentials. Collective modes of the system will show up as resonances of these responses.

The Hamiltonian density is given by

$`H`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m_b}}\psi _b^{}\psi _b\mu _b\psi _b^{}\psi _b+{\displaystyle \frac{\mathrm{}^2}{2m_f}}\psi _f^{}\psi _f\mu _f\psi _f^{}\psi _f`$ (3)
$`+{\displaystyle \frac{1}{2}}g_{bb}\psi _b^{}\psi _b^{}\psi _b\psi _b+g_{bf}\psi _b^{}\psi _f^{}\psi _f\psi _b`$
$`+\psi _b^{}\psi _bV_b^{\mathrm{ext}}+\psi _f^{}\psi _fV_f^{\mathrm{ext}}`$
where the subscripts $`b`$ and $`f`$ denote bosons and fermions respectively, $`\psi _f`$, $`\psi _b`$ are the field operators, $`m_b`$, $`m_f`$ the masses, $`\mu _b`$, $`\mu _f`$ are the chemical potentials, $`V_b^{\mathrm{ext}}`$ and $`V_f^{\mathrm{ext}}`$ the external potentials. All $`\psi `$’s and $`V^{\mathrm{ext}}`$’s are implicitly at the same physical point $`\stackrel{}{r}`$ in space. The interaction parameters $`g_{bb}`$ and $`g_{bf}`$ are related to the scattering lengths $`a_{bb}`$ and $`a_{bf}`$ by $`g_{bb}=4\pi \mathrm{}^2a_{bb}/m_{bb}`$ and $`g_{bf}=2\pi \mathrm{}^2a_{bf}/m_r`$ where $`m_r`$ is the reduced mass \[ $`m_r^1m_f^1+m_b^1`$\].

We shall treat the interaction $`g_{bb}`$ and $`g_{bf}`$ within the Bogoliubov and random phase approximation respectively . The results can be written in the physically transparent form:

$`\delta n_b(q,\omega )`$ $`=`$ $`\chi _b[g_{bf}\delta n_f+V_b^{\mathrm{ext}}]`$ (4)
$`\delta n_f(q,\omega )`$ $`=`$ $`\chi _f[g_{bf}\delta n_b+V_f^{\mathrm{ext}}]`$ (5)
expressing the response of the bosons and fermions to the potentials due to the other species and the external perturbations (the terms in the square brackets). Here $`\delta n_b(q,\omega )`$, $`\delta n_f(q,\omega )`$ are the deviations of the bosonic and fermionic densities from equilibrium at wavevector $`q`$ and frequency $`\omega `$,

$$\chi _b=\frac{1}{g_{bb}}\left[\frac{c_b^2q^2}{\omega ^2c_b^2q^2(q^2/2m_b)^2}\right]$$
(6)
and

$$\chi _f=N_f\left[1\frac{\omega }{2v_fq}\mathrm{ln}(\frac{\omega +v_fq}{\omega v_fq})\right]$$
(7)
are the ($`q`$ and $`\omega `$ dependent) responses of the pure bosons and fermions systems respectively to effective external potentials. $`N_fp_fm_f/2\pi ^2`$ is the density of states for the fermions. ($`p_f=(6\pi ^2n_f)^{1/3}`$ is the fermi momentum, $`v_f=p_f/m_f`$) For simplicity, in eq (7) I have already left out terms that are small if $`q<<p_f`$. $`\omega `$ should be interpreted as having a small and positive part.

eq(5) can be re-arranged as

$$\left(\begin{array}{cc}1& g_{bf}\chi _b\\ g_{bf}\chi _f& 1\end{array}\right)\left(\begin{array}{c}\delta n_b\\ \delta n_f\end{array}\right)=\left(\begin{array}{c}\chi _bV_b^{\mathrm{ext}}\\ \chi _fV_f^{\mathrm{ext}}\end{array}\right)$$
(8)
Then finally

$$\left(\begin{array}{c}\delta n_b\\ \delta n_f\end{array}\right)=\frac{1}{1g_{bf}^2\chi _b\chi _f}\left(\begin{array}{cc}1& g_{bf}\chi _b\\ g_{bf}\chi _f& 1\end{array}\right)\left(\begin{array}{c}\chi _bV_b^{\mathrm{ext}}\\ \chi _fV_f^{\mathrm{ext}}\end{array}\right)$$
(9)
In the case where $`g_{bf}=0`$, $`\delta n_b=\chi _bV_b^{\mathrm{ext}}`$ and $`\delta n_f=\chi _fV_f^{\mathrm{ext}}`$ and the responses thus reduce to those of the pure bose and fermi gases. The corresponding formulas for $`\chi _b`$ and $`\chi _f`$ were already given in eq (6) and eq (7) above. Before we proceed we shall recall the behavior of these responses  and thus the collective modes. For simplicity we shall restrict ourselves to small wavevectors, i.e. $`q<<m_bc_b`$ and $`p_f`$, and without loss of generality $`\omega >0`$. The bosonic response $`\mathrm{Im}\chi _b`$ consists of a delta function at the excitation frequency $`\omega =c_bq`$. This is due to the Bogoliubov mode which is purely propagating and undamped. For the fermions however, there is no collective behavior. The absorptive part, $`\mathrm{Im}\chi _f`$, is finite for a whole range of frequencies $`|\omega |<v_fq`$, known as the particle-hole continuum, arising from the many possibilities of independent particle-hole excitations. $`\mathrm{Re}\chi _b`$ is simple. It is given by $`g_{bb}^1`$ at $`\omega =0`$ and diverges to $`\pm \mathrm{}`$ as $`\omega c_bq`$ from below and above respectively. $`\mathrm{Re}\chi _f`$ is given by $`N_f`$ at $`\omega =0`$. It decreases with increasing $`\omega `$, changes sign at around $`\omega 0.83v_fq`$ and approaches $`\mathrm{}`$ as $`\omega v_fq`$ from both above and below. For $`\omega >v_fq`$, it remains negative with its magnitude gradually approaching zero as $`\omega \mathrm{}`$.

Now we return to the bose-fermi mixture. The response $`\delta n_b`$ to an external potential $`V_b^{\mathrm{ext}}`$ acting on the bosons only is given by $`\chi _b/(1g_{bf}^2\chi _b\chi _f)`$. The existence and the dispersion of the bosonic collective mode are determined by the solution to the equation $`(\chi _b)^1g_{bf}^2\chi _f=0`$, i.e.

$$[\omega ^2+c_b^2q^2+(\frac{q^2}{2m_b})^2]\left(\frac{g_{bf}^2}{g_{bb}}\right)(c_b^2q^2)\chi _f=0$$
(10)
It will be convenient to discuss the normalized response

$$\stackrel{~}{\chi }_bg_{bb}\chi _b/(1g_{bf}^2\chi _b\chi _f)$$
(11)
$`\stackrel{~}{\chi }_b=1`$ in the static limit ($`\omega =0`$, $`q0`$) when there is no boson-fermion interaction ($`g_{bf}=0`$).

Similiarly the fermionic response to an external potential acting on the fermions alone is $`\chi _f/(1g_{bf}^2\chi _b\chi _f)`$. We shall discuss the behavior of the normalized quantity

$$\stackrel{~}{\chi }_fN_f^1\chi _f/(1g_{bf}^2\chi _b\chi _f)$$
(12)
The normalization is chosen such that $`\stackrel{~}{\chi }_f=1`$ in the static limit ($`\omega =0`$, $`q0`$) when there is no boson-fermion interaction ($`g_{bf}=0`$).

Before proceeding let us first examine the responses at $`\omega =0`$. Stability requires that the density responses $`\chi _b/(1g_{bf}^2\chi _b\chi _f)`$ and $`\chi _f/(1g_{bf}^2\chi _b\chi _f)`$ be positive. Using the $`\omega =0`$ values of $`\chi _b`$ and $`\chi _f`$ above, these necessary conditions can be rewritten as $`g_{bb}>0`$ and $`WN_fg_{bf}^2/g_{bb}<1`$. Using the expression of $`N_f`$ given earlier, the last inequality gives $`n_f^{1/3}g_{bf}^2<\frac{2}{3}Ag_{bb}`$ where $`A\frac{\mathrm{}^2}{2m_f}(6\pi ^2)^{2/3}`$ as defined in . These conditions were derived earlier in  and  using slightly different considerations. For bosons and fermions with similar masses, we shall see shortly that $`W`$, a dimensionless parameter, serves as a useful measure of the coupling between the bosons and fermions. If the bosons and fermions have similar masses, $`|W|`$ is of order $`|a_{bf}^2/a_{bb}n_f^{1/3}|`$ and thus typically small for dilute gases unless $`|a_{bf}|>>|a_{bb}|`$. We shall limit ourselves only to the cases where $`|W|`$’s are small.

We shall discuss now the behavior of $`\stackrel{~}{\chi }_b`$ and $`\stackrel{~}{\chi }_f`$ in turn. The results are qualitatively different depending on whether $`c_b\genfrac{}{}{0pt}{}{>}{<}v_f`$. The velocity ratio $`uc_b/v_f`$ can be re-expressed as $`u=\frac{m_f}{m_b}\frac{(4/3)^{1/3}}{\pi ^{1/6}}\frac{(n_ba_{bb})^{1/2}}{n_f^{1/3}}`$. The value of $`u`$ can basically be arbitrary without violating any stability criterion (not only the linear stability condition above but also others derived in )

Bosonic Response:

1. $`c_b>v_f`$: In this regime a propagating bosonic mode exists. It can be easily verified (e.g. graphically) that the mode frequency $`\omega `$ satisfies $`\omega >c_bq`$ ( $`>v_fq`$). The original bosonic mode at $`\omega =c_bq`$ is pushed upwards by the particle-hole ‘modes’ lying below. Some examples are shown in Fig 1. This mode ‘repulsion’ is generally expected (c.f. coupled harmonic oscillators). It is however of interest to examine the microscopic nature of the mode. At the mode frequency both $`\chi _b`$ and $`\chi _f`$ are negative. Thus, e.g., if $`g_{bf}>0`$, $`\delta n_b`$ and $`\delta n_f`$ are of the same sign (see eqn (8)). The repulsion between the two species provides the enhanced restoring force and oscillation frequency. This frequency shift is typically small since usually $`W<<1`$.

2. $`c_b<v_f`$: In this case the original bosonic mode lies inside the particle-hole continuum of the fermions. The bosonic mode is thus Landau damped. For weak-coupling the damping, thus the width of the response, can be estimated easily using eq(10) to be $`\left[\frac{\pi N_fg_{bf}^2}{4g_{bb}}\right]\left[\frac{c_b}{v_f}\right](c_bq)`$. Examples are shown in Fig 2. There is a small shift of the mode due to $`\mathrm{Re}\chi _f`$. The shift is towards higher frequency for $`u`$ sufficiently close to $`1`$ but opposite otherwise ($`\mathrm{Re}\chi _\mathrm{f}<(>)0`$ for $`\omega /v_fq>(<)0.83`$. )

3. It is also of interest to study the bosonic mode for $`g_{bb}<0`$. This is in fact the case for the 6Li-7Li mixture investigated in ref , where the 7Li bosons have a negative scattering length of $`1.5\mathrm{nm}`$. In this case the original bosonic system is unstable, and the Bogoliubov mode has an imaginary frequency for sufficiently small wavevector ( $`q<q_c=2m_b|c_b|/\mathrm{}`$, here $`|c_b|[|g_{bb}|n_b/m_b]^{1/2}`$). Since $`N_fg_{bf}^2>0>g_{bb}`$, the system is still unstable in the presence of the fermions  (see also above). Of interest is the effect of the fermions on the unstable mode. Now for imaginary frequencies $`\omega =i\alpha `$, $`\chi _f(q,i\alpha )=N_f\{1\frac{\alpha }{v_fq}[\frac{\pi }{2}\mathrm{tan}^1\frac{\alpha }{v_fq}]\}`$ is purely real and positive. $`\chi _f`$ decreases monotonically with $`\alpha `$ from $`\chi _f=N_f`$ at $`\alpha =0`$ to $`0`$ as $`\alpha \mathrm{}`$. It can be easily verified that there is a real solution for $`\alpha `$ to the dispersion relation (c.f. eq (10) )

$$[\alpha ^2|c_b|^2q^2]\left(\frac{g_{bf}^2}{|g_{bb}|}\right)(|c_b|^2q^2)\chi _f(q,i\alpha )=0$$
(13)
for sufficiently small $`q`$ (which includes the physically most relevant region where $`\alpha `$ attains its maximum, i.e., the fastest growing instability). Thus the instability is not damped by the particle-hole degree of freedom. In fact it can be verified easily that, for given $`q`$, $`\alpha `$ is increased in the presence of the fermions. The system has become even more unstable. This mode has $`\delta n_b`$ and $`\delta n_f`$ of opposite signs and corresponds to phase-separation as expected.

Fermionic Response:

1. $`c_b>v_f`$: In this case the fermionic response for $`0<\omega <qv_f`$ is only slightly modified. A new feature appears near $`\omega c_bq>v_fq`$ due to the coupling to the bosonic mode. An example is as shown in Fig 3.

2. $`c_b<v_f`$: In this regime there are two important features of the fermionic response. If $`u=\frac{c_b}{v_f}`$ is sufficiently close to $`1`$, the imaginary part contains a sharp resonance at $`\omega `$ above the particle-hole continuum (Fig 4). There are two ways of understanding this mode. It can be regarded as the continuation of the situation from $`c_b>v_f`$, i.e., it is due to the bosonic mode which is itself slightly pushed up in frequency (c.f., Fig 3, note in particular the result for $`u=1`$). Alternatively, this mode can be considered as a zero-sound mode induced by the bosons. The form for $`\stackrel{~}{\chi }_f`$ in eq (12) is precisely that of an interacting fermi gas with s-wave interaction $`g_{ff}`$ (and therefore necessary with more than one spin species, where the response is given by $`\chi _f/(1+g_{ff}\chi _f)`$) though with an effective frequency dependent interaction $`g_{ff}g_{bf}^2\chi _b`$, i.e. an effective s-wave Landau parameter given by $`F_0W/[\left(\frac{\omega }{c_bq}\right)^21]`$. The bosonic mode $`\omega c_bq`$ for $`c_b`$ sufficiently close to but below $`v_f`$ will thus induce a zero-sound mode for the fermions just like an interaction among the fermions will.  Note however there cannot be a real s-wave interaction among the fermions as they are of equal spins. Thus this mode cannot be obtained by considering the effective interaction among the fermions as in ref .

The frequency of this propagating mode can be estimated by using the well-known dispersion relation of the zeroth sound $`\omega /v_fq1+2e^{2(1+\frac{1}{F_0})}`$ with the effective $`F_0W/[(\frac{1}{u})^21]`$ as suggested above. In order for the velocity of the mode to be say $`1\%`$ above $`v_f`$, then $`c_b`$ has to be within around $`7\%`$ of $`v_f`$ if $`W=0.1`$. This estimate agrees with the numerical results of Fig 4.

The second interesting feature is that near the original bosonic mode frequency $`\omega c_bq`$, there is a reduction in the absorptive part $`\mathrm{Im}\stackrel{~}{\chi }_f`$ (see Fig 6). In fact $`\mathrm{Im}\stackrel{~}{\chi }_f0`$ as $`\omega c_bq`$. This, as well as the corresponding behavior of $`\mathrm{Re}\stackrel{~}{\chi }_f`$, can be seen easily mathematically from eq(8) due to the resonance nature of $`\chi _b`$ at this frequency. Physically this can be regarded as due to mode-mode repulsion – the bosonic mode has pushed away the particle-hole ‘modes’ near $`\omega c_bq`$. This feature is present even for small coupling $`W`$. A larger $`W`$ mainly increases the width of this ‘transparent’ region. Thus in fact the frequency dependence of $`\mathrm{Re}\stackrel{~}{\chi }_f`$ is actually stronger for smaller $`W`$’s.

The energy absorption by the bose-fermi mixture from an external perturbation acting on the fermions is thus substantially reduced for frequencies within this ‘transparent region’. The width of this region can be estimated by using the observation that the fermionic response is roughly reduced by the factor $`1+W\frac{(c_bq)^2}{\omega ^2(c_bq)^2}`$ for these frequencies. For the fermionic response to be reduced to say $`1/2`$ of its bare value, then $`|\omega c_bq|/c_bq<W/2`$. This estimate agrees very roughly with the numerical results in Fig. 6.

In conclusion, I have investigated the collective modes of a Bose-Fermi mixture, and have shown that there is important mode-mode coupling effects, especially if $`v_fc_b`$.

I thank David Edwards for helpful correspondences. This research was supported by the National Science Council of Taiwan under grant number 89-2112-M-001-105.