text:
# Flat-Band Ferromagnetism in Organic Polymers Designed by a Computer Simulation

## I Introduction

Materials design to realize desired properties in condensed-matter physics is becoming increasingly realistic. One important factor is, given high levels of computer performance, the computer simulation is an essential part in designing. Namely, computer simulations with first principles calculations should be imperative in narrowing down, or even pin-point, from a wide variety of candidate or designed materials for the interesting properties and functions we search for. Such an approach is especially promising for molecules and nano-structures. While the atom manipulation with scanning tunneling microscopy (STM)Eigler ; Shen ; Hashi1 ; Taro2  provides new possibilities for materials design, the design of molecules, especially polymers as we focus on here, should be fundamental.

In this context organic molecules and polymers are of special interest, since they have versatile structures and chemical properties that can be wider than in inorganic materials. Indeed, the discovery of a conducting organic polymer by ShirakawaShirakawa74 ; Shirakawa77  kicked off intensive studies that are paving a new way to adopt organic polymers or oligomers in realizing various functions as molecular-electronics devices such as field-effect transistorsGarnier94  or light-emitting diodesFriend99 .

Organic ferromagnets have also attracted much attention as a challenging targetAllemand91 ; Rajca . In particular, organic magnets consisting entirely of non-magnetic elements is of fundamental as well as practical interests. Ordinary ferromagnets consisting of magnetic elements exploit electrons in d- or f-orbitals that are strongly interacting. Can magnetism arise in p-electron systems that are weakly interacting? One theoretical possibility is to apply the flat-band ferromagnetism proposed by Mielke and by TasakiMielke1 ; Tasaki1 . This arises as an effect of electron-electron repulsion when the (one-electron) band structure contains a dispersionless band. The mechanism is interesting in many ways, but essential features are, first, the system is totally distinct from the “narrow-band limit” in textbooks, since the magnetism occurs when the transfer between different sites is finite.

Second, this is a band ferromagnetism rather than a magnetism arising from the exchange interaction between localized spins. So, this mechanism should be a good candidate for designing organic ferromagnets. While there are also other possible mechanisms such as the intra-molecular Hund’s coupling for spins in p-orbitalsAllemand91 , the flat band ferromagnetism is advantageous in that spins are carried by itinerant electrons which can be utilized for spin injectors in spintronics.

Here a theoretical remark on the flat-band ferromagnetism is due. The flat-band ferromagnetism was first proposed by Lieb,Lieb  who considers bipartite lattices (consisting of two sublattices) that have different numbers of A sublattice sites $`N_A`$ and B sublattice sites $`N_B`$. Lieb proved that, when we switch on the electron-electron repulsion (assumed to be short-ranged, so that we take the Hubbard model), the ground state should then be ferrimagnetic with the magnetization $`N_AN_B`$. Quantum chemically, this model contains $`N_AN_B`$ non-bonding molecular orbitals, so is similar to Mataga’s model. Shima and AokiShima93  then proposed a systematic way to realize such systems as superhoneycomb structures. MielkeMielke2 ; Mielke3 ; Mielke4  and independently TasakiTasaki4  then constructed other flat-band ferromagnetism, which is distinct from Lieb’s in that the flat band is constructed from quantum mechanical interference between the nearest-neighbor and further transfers (so the system is necessarily non-bipartite) and the ground state is now ferromagnetic.

To be more precise, the lattice is required to satisfy what is called the “local connectivity condition”. Namely, the flat band is by no means a sufficient condition for ferromagnetism, and the magnetism arises for special lattices on which adjacent “Wannier” orbitals have to overlap with each other no matter how they are combined to minimize each orbit size. So, while ordinarily a flat band implies a disjointed set of orbits, they are connected in Mielke-Tasaki lattices. This is intuitively why spins align to lower the repulsion energy from Pauli’s principle. Rigorous proof for the ferromagnetic ground state has been given for the half-filled flat band.

In actual materials, it is very hard to realize perfectly flat bands, and to make them perfectly half-filled. For the former, a rigorous proof (for Mielke-Tasaki lattice)Tasaki3  and numerical simulationsSakamoto ; Arita98  showed that a slight band dispersion does not destroy the magnetism. For the latter, there are also a proofMielke1  and a numerical simulationSakamoto  which revealed that small deviations from half-filling are allowable for the ferromagnetism.

In the context of this background, what kind of organic materials are promising for realization of the flat-band ferromagnetism? The hardest part is to realize the connectivity condition. In general we have to consider complex lattices such as Kagomé, or lattices having distant-neighbor transfers. However, the condition is easier to satisfy on one-dimensional polymers. For example, Mielke-Tasaki model is realized as a chain of triangles.Penc ; Sakamoto ; Arita3  So we can concentrate on polymers.

Next comes the choice of the monomers that should be polymerized. Even-membered rings such as benzene are disadvantageous in that antiferromagnetic order tends to occur, which should dominate over the ferromagnetism. So we should opt for odd-membered rings, and the above example of the chain of triangles is indeed an example of this. Since triangular molecules are scarce, we can conclude that we should focus on polymers based on five-membered rings, such as polypyrrole, polythiophen, etc.

For the computer design, the purpose of the present paper, there are two tasks: (i) to search for the materials that realize the flat band a la Mielke-Tasaki. Since the chain of odd-membered rings is by no means a sufficient (nor necessary) condition for the flat band, this is an important part.

(ii) Next, we should do a first-principles (spin density functional) band calculations to confirm (a) the flatness of the band and the ferromagnetism in the ground state, and (b) whether the magnetism can indeed be interpreted as the Mielke-Tasaki mechanism, for which we have to evoke an electron-correlation calculation (exact diagonalization of the Hubbard model here). (c) We also address the question of how the departure from half-filling affects the magnetism in the flat band.

For the first-principles calculation, we have employed the generalized gradient approximationPBE96  (GGA). Depending on the purpose of each calculation, we have performed either the spin-unpolarized case based on the density functional theory (GGA-DFT) or the spin-polarized one based on the spin density functional theory (GGA-SDFT). We used the plane-wave-based ultrasoft pseudopotentialsVanderbilt ; Laasonen . The energy cutoff was taken as 20.25 Rydberg. The convergence criterion of the geometry optimization was that all of the forces acting on each atom were within $`1\times 10^3`$ Hartree/a.u.

The plan of this paper is as follows. In section II, we describe the search for various materials and show some examples which exhibit possibilities of ferromagnetism. In section III, we show the result for the most promising material, poly-aminotriazole, and discuss how the theory for the flat-band ferromagnetism works on this particular material. In section IV, we discuss how the ferromagnetic state is robust against the deviation of the band filling from the half-filling.

## II Examined materials

### II.1 Flat bands in the polymer of five-membered rings

Let us first examine the simple tight-binding model,

$$H=t\underset{i,j}{}c_i^{}c_j+\epsilon _0\underset{i=2}{\overset{5}{}}n_i+\epsilon _1n_1,$$

for a chain of five-membered rings, where we first assume that all the transfer integrals within the ring and those connecting rings have the same $`t`$, all the on-site energies are the same ($`\epsilon _0=0`$), except for the one at the top of the ring ($`\epsilon _1`$). We first consider the half-filled case where one electron per site on average. This is a reasonable assumption when we consider $`\pi `$-orbital networks on such polymers.

Let us start from an observation that Mielke-Tasaki’s condition is satisfied when $`\epsilon _1=t`$ ($`>0`$), a not unrealistic condition. In this case, the third band with $`E=0`$, i.e., the half-filled band, has no dispersion. This provides a heuristic example for the Mielke-Tasaki eigenfunctions: one can construct the most compact eigenfunction which consists of two rings as depicted in Fig. 1(a), where the amplitudes of the numbered sites are given as (1, 1, 0, -1, 0, -1, 0, 1, 0, -1). The fact that the amplitudes at the sites 5 and 7 are zero ensures the localized nature of this eigenfunction. Note that we are not displaying a part of a Bloch function, but the whole eigenfunction. One can also construct another, linearly-independent eigenfunction by simply shifting the wave function from one ring to the next. The eigenvalue of the new function is the same ($`E=0`$). So the set of all these eigenfunctions can be a basis for the flat band. Because the eigenfunction is not confined in a single ring (unit cell) but extends over the two rings, two neighboring bases (eigenfunctions) always overlap with each other as shown in Fig. 1(b). One cannot remove such overlaps between the bases no matter how the linear combination of these bases are taken. This overlap is the origin of the ferromagnetic coupling between the electrons in the flat band.

### II.2 Known polymers of five-membered rings

A typical polymer of five-membered rings is polypyrrole. The atomic structure and calculated band structures are shown in Figs. 2(a) and (b), respectively. The five-membered rings alternate their directions in this compound, so that the unit cell contains two rings. X point in Fig. 2(b) has $`𝒌=(\pi /a,0,0)`$, where $`a`$ is the lattice constant of this unit cell as shown in Fig. 2(a). Because of the doubled unit cell, most of the bands are folded at X. For directions ($`y,z`$) perpendicular to the chain, we have taken a cell size large enough to avoid inter-chain interactions. We can see that polypyrrole has, not surprisingly, no flat bands around the Fermi energy.

We have also looked at other typical polymers of five-membered rings, including polythiophen (where an N-H block in polypyrrole is replaced with an S atom) and polytriazole (where all the C-H blocks in polypyrrole are replaced with N atoms). However, their band structures have no flat bands, either.

### II.3 Designing the polymer

This is exactly where the designing comes in. Our strategy consists of two approaches. One is to replace H atom bound to N atom at the top of the ring with various kinds of bases. According to the tight-binding calculation described in the previous subsection, the on-site energy ($`\epsilon _1`$) of the top of the ring should be higher than that of C atoms to have a flat band. Because the on-site energy of N is expected to be lower than that of C, the first attempt should be the replacement of the N atom with other elements with higher on-site energies. However, that should be difficult if one wants to retain the existence of the $`\pi `$-orbital and single electron occupation on it. Therefore, we have opted for controlling the on-site energy by replacing the H atom with other functional groups, by leaving the N atom intact.

Alternatively, we can replace C-H blocks on the bottom edge of the pentagons (3, 4, 8, 9 in Fig. 1(a)) with N atoms, where the lowering of the on-site energies of the bottom of the ring may effectively realize the required condition. In other words, we use polypyrrole and polytriazole as the platform to modify either the top or the bottom of the ring. Since it is not obvious which approach should be better, we consider both platforms for all the substituents considered here.

Substituents we have tested are sodium (Na), potassium (K), chlorine (Cl), fluorine (F), cyanogen (CN), nitro (NO2), sulfate (SO4), carboxyl (COOH), amino (NH2), methyl (CH3), and hydroxyl (OH). For all these candidates, we have performed first-principles (GGA-DFT) optimization of atomic structures, and obtained their band structures. Most of them have no flat bands around $`E_F`$, but some of them have turned out to have flat bands. For these we have performed further (GGA-SDFT) first-principles optimizations under the doping condition to make the flat band half-filled.

First, we have examined single atoms (alkali metal or halogen) as substituents. This is motivated by the fact that small substituents may be better than larger ones, because a bulky group at the top of the ring may introduce a coupling (transfer) with the bottom of the adjacent ring, which would destroy the basic assumption that the adjacent rings are connected by one bond to satisfy Mielke-Tasaki’s condition.

In order to raise the on-site energy of the top of the ring, we have first considered an alkali-metal atom, Na, because of its low electron affinity. Figure 3(a) shows the atomic structure of poly-sodium-triazole, where hydrogen atoms in polytriazole are replaced by sodium atoms. Figure 3(b) shows the band structure, where one can see that a nearly flat band exists above the Fermi energy. Absence of any unnecessary bands between the flat band and $`E_F`$ is good sign, because we can dope carriers into the flat band without any complications.

We have also examined the Na-substituted polypyrrole, and obtained a similar band structure. However, in this case the flat band above $`E_F`$ is much closer to a dispersive band below $`E_F`$, which can introduce unwanted inter-band interactions. As for the alkali-metal we have also considered potassium, but the case of Na has turned out to be better.

We move on to the doping. First, the flat bands always appear as a pair of bands folded at X point, except for the case of antiferromagnetic ordering. Therefore, we make the flat-band half-filled by doping two carriers per unit cell into the pair of flat-bands. Namely, one carrier is doped for each ring. In poly-sodium-triazole, we have considered the doping of two electrons per unit-cell. Here we realize the doping condition by increasing the number of electrons with a uniform positive background charge for charge neutrality.

After a geometrical optimization of doped poly-sodium-triazole, we have calculated the spin-dependent (GGA-SDFT) band structure in Fig. 3(c). One can see that the number of occupied states differs between major- and minor-spins, which means that the system is ferromagnetic. It should be noted that (the pair of) the flat band is not entirely spin-split. The difference of the number of spins is 0.65 per unit cell. The spin polarization is smaller than the expected polarization (2.0). We also found that the total energy (per unit cell throughout this paper) of the ferromagnetic (F) state is lower by 12 meV than that of the paramagnetic (P) state. We could not find antiferromagnetic (AF) solution even when we started geometry optimization from an AF electronic state. So the F state is the most stable in this material, but is only slightly lower in energy than the P state.

The reason why this material gives such an insufficient result in spite of the existence of the nearly flat band should be sought in the nature of the flat-band’s wave function. We found that the half-filled flat band here is made of an orbital localized around the Na atom, so does not satisfy the local connectivity condition. This should be why the stability of the F state is very weak, while the finite spin-polarization in GGA-SDFT should be a narrow-band effect rather than a much more robust Mielke-Tasaki effect. Indeed, the difference between the weak ferromagnetism here and a robust ferromagnetism in poly-aminotriazole, shown later, is a good demonstration of the importance of the local connectivity condition, hallmark of the flat band ferromagnetism.

Now we digress a bit, and consider the contrary case of an atom having high electron affinity such as a halogen, e.g., fluorine atom, which might be heuristic. Figure 4(a) shows the atomic structure of poly-fluoropyrrole, where the fluorine atoms are substituted for hydrogen atoms bound to nitrogen atoms in polypyrrole. Figure 4(b) is the band structure obtained by spin-unpolarized (GGA-DFT) calculation. We find a flat band which lies just in between two dispersive bands below $`E_F`$. Because of this the doping greater than one carrier per unit cell is required.

So we consider making the flat-band half-filled by doping four holes per unit-cell. Here, two holes are necessary to make the dispersive band just below the Fermi level empty, and two holes are necessary to make the folded flat bands half-filled. Figure 4(c) shows spin-dependent band structure after the geometry optimization under the doping condition. The system is seen to be ferromagnetic. The difference of the number of spins per unit cell is just 2.0. One problem of this material, however, is that doping twice as large is necessary. When the doping is performed in this material, the flat band shifts upward (downward) for minority (majority) spin. Replacing the fluorine atom with a chlorine atom results in a band structure similar to the case of fluorine. In fact an AF state, obtained by starting geometry optimization from an antiferromagnetic state, has a total energy 21 meV lower than that of F state.

We move on to functional groups as substituents to search further possibilities. We have considered an example, poly-hydroxytriazole, where H atom in polytriazole is replaced by OH with a low electron affinity. Although a pair of slightly dispersive bands exist just below $`E_F`$, these bands, when doped with two holes per unit cell, become mixed with dispersive bands underneath. As a result, the difference of the number of spins in the F state is only 0.39 per unit cell, and the paramagnetic state is more stable than the F state by 3 meV.

Figure 5(a) is the atomic structure of poly-aminopyrrole, where H atom bound to N atom in polypyrrole is replaced by NH2, whose electron affinity is low. In the band structure, Fig. 5(b), one can find a pair of flat bands below $`E_F`$. The lower half of it is very flat while upper half is slightly dispersive.

We have then doped two holes per unit cell. The band structure after geometry optimization is shown in Fig. 5(c). The difference of the number of spins is 1.0 per unit cell. As one can see from Fig. 5(c), the Fermi level intersects the middle of the pair of flat bands. Although the lower part of the pair of flat bands is quite flat, the dispersion of the upper half part is too large. Here only lower half part seems to work as a flat-band. The total energy of paramagnetic state in this material is found to be 65 meV higher than that of ferromagnetic state, while an AF solution does not exist, so the F state is the most stable.

All other substituents tested here except for the one in the next section have turned out to be inappropriate. In most of those materials, the band which we expect to be flat is either too dispersive or too deep. For the remaining ones, band structures are qualitatively similar to the examples described in this section. Poly-aminotriazole, discussed in the next section, also gives qualitatively similar results to poly-aminopyrrole shown in this section, but quantitatively much better.

## III Poly-aminotriazole

### III.1 First principles calculation

The best candidate for ferromagnetic organic polymer among the materials we have considered is poly-aminotriazoleArita1 . Figure 6(a) shows the atomic structure of poly-aminotriazole, where H atom in polytriazole is replaced with NH2. The difference between poly-aminopyrrole shown in the last section and this material is whether the platform is polypyrrole or polytriazole, i.e., the C-H blocks in the bottom of the ring is replaced with N atoms. The band structure of the poly-aminotriazole, shown in Fig. 6(b), is similar to that of Fig. 5(b), but the dispersion of the upper half of the pair of flat bands is smaller and the separation from dispersive bands below is greater than that of Fig. 5(b). These features are desirable for the flat-band ferromagnetism.

After the optimization of the atomic structure under the doping condition of two holes per unit cell, the spin dependent band structure was calculated and shown in Fig. 6(c). The result hits on the ideal situation — The pair of flat bands is made half-filled, and as a result, fully occupied by majority spins while totally unoccupied by minority spins. Reflecting this situation the difference of the number of spins is the desired 2.0 per unit cell.

In this material, we found a metastable antiferromagnetic state. The band structure of the AF state is shown in Fig. 7. This is calculated with GGA-SDFT, but the resulting band structures for up- and down-spins are the same, while the wave functions for up- and down-spins are arranged alternately on the chain of rings. This means that even though we have avoided intra-ring AF instability by choosing five membered rings, inter-ring AF can exist at least as a metastable state. The upper and lower parts of the flat band folded at X point is separated with $`E_F`$ in between. The AF state is found to be 52 meV higher in energy than the F state. Paramagnetic state is even higher (by 384 meV) than the F state. Therefore F state is the most stable in this material.

Let us discuss how the flat-band ferromagnetism is achieved in this material by combining the wave functions of the flat band obtained by first principles calculation with the calculation in the tight-binding model. Figure 8 shows the majority-spin’s wave functions in the pair of flat bands just below the Fermi level in Fig. 6(c). While the phases of the wave functions can be taken arbitrarily, we choose for clarity the wave functions of the lower and upper bands at $`\mathrm{\Gamma }`$ point (Figs. 8(a) and (f)) to be real and pure imaginary, respectively, and those at intermediate $`k`$-points to be continuous. If we first look at the real part, typically for Fig. 8(a), the wave function for each monomer (aminotriazole) has two nodal lines which divide the amplitude into three parts, corresponding to NH2, C-N-C, and N-N blocks. For the imaginary part, typically for Fig. 8(f), the wave function has one nodal line, which divides the amplitude into two, consisting of two C-N blocks. The two wave functions (real and imaginary), which are orthogonal, are mixed at the intermediate $`k`$-points.

### III.2 Tight-binding model

In order to capture the essence of the flat band, here we introduce a tight-binding model which represents $`\pi `$-orbital’s network in this system. Figure 9(a) shows its band structure. The parameters of the model displayed in Fig.9(b) are fitted so as to reproduce Fig. 6(b). Obtained values are: $`t_{\mathrm{CN}}=t_{\mathrm{CC}}=t_\mathrm{f}=2.5`$ eV, $`t_{\mathrm{NN}}=3`$ eV, $`\epsilon _0=1.4`$ eV, and $`\epsilon _1=0.5`$ eV. One can see that the tight-binding band (Fig. 9(a)) excellently agrees with the $`\pi `$-bands obtained by first-principles calculation (solid lines in Fig. 6(b)).

Figure 9(c) shows the tight-binding wave function of the eigenstate corresponding to that in Fig. 1(a). In this material, one site is added to the top of the ring which represents a nitrogen atom of NH2. The open and closed circles indicate the sign of the wave function, respectively, and their sizes amplitudes. This wave function satisfies the local connectivity condition.

Here it should be noted that the substitution of NH2 for H not only raises the on-site energy of the top of the ring but also introduces an extra $`\pi `$-orbital of the N in NH2. The existence of the extra site completely changes the connectivity condition, so that we have no longer to satisfy $`\epsilon _1=t>0`$, which has turned out in the present study to be rather difficult condition to realize.

Figure 10 shows $`k`$-dependence of the tight-binding wave functions of the flat band. Here $`k`$-space is extended to the second Brillouin zone, so the left and right halves of the figure represent the wave functions of the lower and upper parts of the pair of flat bands, respectively. Phase factors are chosen again for clarity as in Fig. 8. The label “Re(2,5)” indicates that the line represents the real part of the wave function at the site 2 and 5, where the numbering of the sites is shown in the insets. On the $`\mathrm{\Gamma }`$-points at the left- and right-edge of the figure, the wave function becomes real and pure imaginary, respectively. These two wave functions are displayed in the inset of Fig. 10, where it should be noted that we display a part of the periodically extending Bloch functions, in contrast to Fig. 9(c) which shows the entire eigenfunction.

One can see that the symmetry and the character of the wave functions shown in Fig. 10 is similar to the first-principles wave function (Fig. 8). These wave functions can be constructed by a linear combination of the eigenfunctions shown in Fig. 9(c). Figure 11 illustrates how the linear combination of the localized eigenfunctions provides the Bloch function. When the two neighboring eigenfunctions, A and B, are added, the sum at the ring on which they overlap (enclosed by dotted lines) forms the shape of the imaginary-part wave function as shown in the inset of Fig.10. When subtracted, it forms the shape of the real-part wave function. The entire imaginary- and real-part wave functions are constructed by the summation of all localized eigenfunctions (A+B+C+D+$`\mathrm{}`$) and the combination of the summation and subtraction (A$``$B+C$``$D+$`\mathrm{}`$), respectively. The real-part wave function and the imaginary-part wave function are orthogonal to each other. The gradual mixing of these two orthogonal wave functions as $`k`$ is changed is the very character of the Mielke-Tasaki flat band, which should never occur on an ordinary flat band.

### III.3 Comparison with the Hubbard model

Having constructed the tight-binding modeling we proceed to the question of whether the ground state is spin-polarized in the presence of the Hubbard interaction, $`_U=U_in_in_i`$. In Fig. 12 we show the phase diagram against $`U`$ and $`\epsilon _1`$ obtained with an exact diagonalization calculation for a 12-site (2 rings) Hubbard model for $`t_0=2.5`$ eV and various values of $`t_{\mathrm{NN}}=3.04.0`$ eV. As indicated in the inset, $`\epsilon _0`$ is chosen throughout to satisfy the condition for the flat band,

$`\epsilon \epsilon _1`$ $`=`$ $`(1\epsilon )(t_{\mathrm{NN}}^2(\epsilon \epsilon _1)^2)t_{\mathrm{NN}},`$
$`(\epsilon \epsilon _1+t_{\mathrm{NN}})/(1\epsilon )`$ $`=`$ $`t_{\mathrm{NN}}/(\epsilon _0\epsilon )t_{\mathrm{NN}}(\epsilon \epsilon _0),`$
where $`\epsilon `$ is the eigenenergy of the flat band and $`t_{\mathrm{CN}}=t_{\mathrm{CC}}=t_0`$ ($`=1`$ here) is assumed for simplicity. We can see that we have indeed a ferromagnetic phase unless the repulsion is too strong (i.e., $`U<U_c`$ with $`U_c=25`$ eV). Our preliminary quantum Monte Carlo calculation confirms that two unit cell is sufficient to roughly determine the ferromagnetic phase boundary. Peculiar destruction of the ferromagnetism above $`U>U_c`$ is another hallmark of the Mielke-Tasaki magnetism.Arita98

## IV Deviation from the half-filled flat band

### IV.1 Comparison of the total energies

In order to realize ferromagnetism in the materials we proposed here, carriers have to be doped into the flat band. We can consider several possible methods to achieve this such as field effect dopingLee , photo-induced carrier injectionHiroi , as well as chemical doping. The first two are now attracting much attention in molecular electronics, but they are very recent techniques and may need further developments. We can instead concentrate on the chemical doping, where the polymers are crystallized along with the accepters, and we found it feasible. Details are presented elsewhereArita2 . Now, if the flat-band ferromagnetism can occur not only just at the half-filling but also away from that, this would make the chemical doping much easier.

We have to know how robustly the ferromagnetism survives when the filling of the flat band deviates from just half-filled. Although there is a proofMielke1  and a numerical studySakamoto  for a Mielke-Tasaki lattice which show that slight deviations do not wash out the ferromagnetism, those conclusions may not be simply applied to other models. For the present model, calculations taking full account of the electronic correlation such as the exact diagonalization of the Hubbard model would be unsuitable, because changing the number of electrons in small systems corresponds to too big a change in the filling. We have instead performed an SDFT calculation, because it can handle the fractional number of electrons (per unit cell) if a uniform background charge is introduced.

At this point we can note the following. In SDFT calculations, electronic correlations are not fully taken into account, so we can question whether the flat band ferromagnetism may be described within the SDFT framework. This problem has been addressed in our previous publicationArita1  where SDFT and exact diagonalization of a Hubbard model have given qualitatively consistent results. This is not surprising, since one of the remarkable features of the flat-band ferromagnetism is that the magnetism exists over the whole range (infinitesimal $`<U<\mathrm{}`$) of the interaction as proved rigorously, so the weak-$`U`$ case (for which the methods such as SDFT is meant to be applicable) crosses continuously over to strong-$`U`$ case. As far as the on-site repulsion is not too strong, the comparison between the gain of the exchange-correlation energy by spin polarization (F state) and the gain of the band energy by forming AF order (or equivalently, a spin density wave (SDW) with the periodicity of two five-membered rings) accompanied by a lattice relaxation should predominantly determine which state is stable. These two energies can be accurately estimated by SDFT calculations. Therefore, we expect here that SDFT can serve our purpose, to see the band-filling dependence of the total energies of F, AF, and P states.

Figure 13 (a) shows band-filling dependence of the total energies of the F and AF states of the poly-aminotriazole measured from that of the P state. The position where the number of holes, $`n_h`$, is 2.0 corresponds to the half-filling. The results for F, AF, and P are obtained by fully relaxing the atomic geometry for each electronic state. It can be seen that F state is stable in the range $`1.1<n_h2.0`$, while P state is stable for $`n_h<1.1`$. In the region where F is stable, AF is found to be metastable. In the P-region, $`n_h<1.1`$, no AF state is found. This result shows that the ferromagnetism in this material is robust, even down to the quarter filling, against the deviation from the half-filling. The result that quarter-filling is acceptable greatly relaxes the difficulty in doping.

### IV.2 Correlation with the atomic configuration

Since the polymer proposed here has functional groups, we should also examine how the three-dimensional atomic configuration is correlated with the band-filling dependence. Namely, all the atoms do not lie on a plane: in poly-aminotriazole two hydrogen atoms in NH2 in the undoped polymer lie above the plane of the rings, while all the atoms are co-planar in the half-filled poly-aminotriazole. We can more closely look at the torsion angle, $`\varphi `$, defined as the angle with which C-N-N-H atoms are connected. Figure 13(b) shows the torsion angles in the P and F states, respectively. From Figs. 13(a) and (b), one can see that the breakdown of the ferromagnetism with the decrease in the number of holes is correlated with the change in the atomic configuration.

When the system is half-filled ($`n_h=2.0`$), all atoms become co-planar, where the N atom of the NH2 base has $`sp^2`$ hybridized orbitals, forming chemical bonds with two H atoms and a N atom at the top of the ring. The network of the $`\pi `$ orbitals in poly-aminotriazole is formed by those of the rings and that of NH2’s. When the number of holes is decreased and becomes $`n_h<1.1`$, the N atom of the NH2 forms $`sp^3`$ hybridized orbitals, i.e., the $`\pi `$ orbital at NH2 disappears. Since the chemical bonds to the N atom and two H atoms are formed in the three directions of $`sp^3`$, leaving one lone pair, NH2 has a three-dimensional structure. This time, the network consists of five $`\pi `$-orbitals of the rings and a $`sp^3`$ hybridized orbital of the lone pair. This change is disadvantageous for the scheme of the flat-band formation (Fig. 9) due to the balance of the six $`\pi `$ orbitals. This is considered why F state always takes a flat geometry while P state a non-flat geometry.

To see what part of the energy affects the change in the stable state when $`n_h`$ is decreased, we have calculated the energy difference between the F and P states by fixing the atomic structure to the optimum one for the F state at $`n_h=2.0`$. The result is plotted as a dotted line in Fig. 13(a). This shows that the value is negative over the whole range plotted here, which indicates that the electronic part of the total energy always favors F state regardless of the filling of the flat band. Therefore, the change in the stable state at $`n_h=1.1`$ can be ascribed to the change in the structure, especially of NH2. If we can find a material in which there is no such change in the atomic structure but the parameter values are similar, we can expect an occurrence of the ferromagnetism even for a low-doping region less than quarter-filling.

In summary, the present work describes an example of materials design from first principles for the purpose of realizing flat-band ferromagnetism in organic polymers. We have found a promising material, polyaminotriazole, for which why it is the best and the robustness of the ferromagnetism has been discussed.

###### Acknowledgements.

We would like to thank T. Hashizume, B. Choi, M. Ichimura and J. Yamauchi for fruitful discussions. We are also indebted to H. Nishihara and Y. Yamanoi for discussions on the chemistry of polymers. This study was performed through Special Coordination Funds for Promoting Science and Technology of the Ministry of Education, Culture, Sports, Science and Technology of the Japanese Government. The first-principles calculations were performed with TAPP (the Tokyo Ab-initio Program Package).


format: text

text:
# Use of the likelihood for measuring the skill of probabilistic forecasts

## 1 Introduction

Users of forecasts need to know:

* whether the forecasts they are receiving have been adequately calibrated
* whether the forecasts they are receiving are any better than an appropriate simple model such as climatology
* which of the forecasts they are receiving is the best

To answer these questions, a single measure of forecast quality is needed. For calibration, the measure serves as a cost or benefit function that must be minimized or maximised in order to find the optimum values for the free parameters in the calibration algorithm. For comparison with climatology or other forecasts, the measure serves as a way of deriving a ranking.

There are many standard measures of forecast quality. For example, for calibrating and comparing single-valued temperature forecasts, mean square error (MSE) is common. For binary probabilistic forecasts, the Brier score (Brier, 1950) is often used. For continuous probability forecasts, the continuous rank probability score and the ignorance have been suggested.

In this paper we will argue that likelihood-based measures provide a simple and natural general framework for the evaluation of all kinds of probabilistic forecast. For example, likelihood based measures can be used for binary and continuous probability forecasts, for temperature and precipitation, and for one lead time or many lead times simultaneously.

In section 2 we define the likelihood and discuss why we think it is a useful measure of forecast skill. In section 3 we include expressions for the likelihood for the normal distribution and in section 4 we discuss relations between the likelihood and other forecast scoring methods. Finally in section 5 we summarise and describe some areas of future work.

## 2 Probabilistic forecasts and the likelihood

How should we evaluate the skill of a probabilistic forecast? We advocate the use of a particular set of measures that are taken from classical statistics, and are all based on the *likelihood*. Likelihood is defined very simply as *the probability of the observations given the forecast*. In this phrase *the observations* refers to the entire set of observations that we have available to validate a certain forecast, and *the forecast* refers to the entire set of corresponding forecasts.

Likelihood was first used by Fisher (1912) as a method for fitting parameters to parametric distributions. Fisher proposed the likelihood as the natural benefit function that one should maximise in order to define the best-fit parameters of the distribution. This suggestion was given a mathematical basis when it was shown that the parameter values that maximise the likelihood are the most accurate possible estimates for the unknown parameters for most problems (see Casella and Berger (2002)).

Fisher’s problem, of how to evaluate the goodness of fit of a distribution to a number of samples, is exactly the same as the problem of how to evaluate a probabilistic forecast. Instead of the distribution we have the probabilistic forecast and instead of the samples we have the verifying observations.

### 2.1 Advantages of the likelihood as a measure for skill

We consider that the likelihood has the following advantages as a measure of probabilistic forecast skill:

* It has a simple definition that, from a purely intuitive point of view, seems to be a reasonable basis on which to compare forecasts
* It is mathematically optimal in the sense that estimates of parameters of calibration models fitted by maximising the likelihood are usually the most accurate possible estimates (see Casella and Berger (2002)).
* It is a generalisation to probabilistic forecasts of the most commonly used skill score for single forecasts: the RMSE (see section 4 below for a discussion of this).
* It also shows how the RMSE can be generalised to the case of autocorrelated forecast errors
* The properties of the likelihood have been studied at great length over the last 90 years: it is well understood
* It is both a measure of resolution and reliability
* Likelihood can be used for both calibration and assessment: this creates consistency between these two operations
* Use of the likelihood also creates consistency with other statistical modelling activities, since most other statistical modelling uses the likelihood. This is important in cases where use of forecasts is simply a small part of a larger statistical modelling effort, as is the case for our particular business.
* Likelihood can be used for all meteorological variables
* Likelihood can be used to compare multiple leads, multiple variables and multiple locations at the same time in a sensible way (giving a single score) even when these leads, variables and locations are cross-correlated

### 2.2 Forecast scores derived from the likelihood

A number of different scores can be derived from the likelihood.

* The log-likelihood (LL) reduces the range of values of the likelihood to a more manageable scale
* Minus the LL (MLL) has the characteristic that better forecasts have lower values: in this way it is analogous to the MSE
* The square root of the MLL (RMLL) has a further compressed scale
* All these measures can be transformed into skill scores from zero to one in the usual way

Other transformations are also possible: for instance, one might consider normalising by the number of data points.

## 3 The likelihood for the normal distribution

For a normal distribution the likelihood is given by:

$$L=\frac{1}{\sqrt{2\pi \text{det}}}\text{exp}(\frac{1}{2}(T\mu )^T\mathrm{\Sigma }^1(T\mu ))$$
(1)
where $`T`$ is the vector of observations, $`\mu `$ is the vector of means from the forecast, $`\mathrm{\Sigma }`$ is the covariance matrix of the forecast errors, and det is the determinant of $`\mathrm{\Sigma }`$.

The log-likelihood is then:

$$l=\frac{1}{2}ln(2\pi \text{det})\frac{1}{2}(T\mu )^T\mathrm{\Sigma }^1(T\mu )$$
(2)
In the case where the forecast errors can be assumed to be uncorrelated in time, the likelihood becomes:

$$L=\frac{1}{\sqrt{2\pi \text{det}}}\text{exp}\left(\frac{1}{2}\underset{i=1}{\overset{i=N}{}}\frac{(T_i\mu _i)^2}{\sigma _i^2}\right)$$
(3)
and the log-likelihood is:

$$l=\frac{1}{2}\underset{i=1}{\overset{i=N}{}}ln(2\pi \sigma _i)\frac{1}{2}\underset{i=1}{\overset{i=N}{}}\frac{(T_i\mu _i)^2}{\sigma _i^2}$$
(4)
When evaluating a forecast using the likelihood, calculating the covariance matrix is straightforward because the forecast errors are known. When calibrating a forecast using the likelihood, calculating the covariance matrix is more difficult. If it is reasonable to assume that the errors are uncorrelated in time, then this simplifies the calibration considerably. However, this is generally not the case.

## 4 Relations between the likelihood and other skill scores

Likelihood is closely related to the RMSE and the ignorance, as we see below.

### 4.1 Relation between the likelihood and RMSE

We show that the RMSE and the likelihood are consistent (i.e. give the same ranking of forecasts) in the case of two normally distributed probabilistic forecasts with different means but the same constant spreads. Likelihood is used to compare the whole distribution, while RMSE is used to compare the means.

Suppose we have two forecasts, $`A`$ and $`B`$, and suppose:

$$L_A>L_B$$
(5)
Taking logs, this gives:

$$l_A>l_B$$
(6)
Substituting in the expression for the log-likelihood for a normal distribution we see that:

$$\frac{N}{2}ln(2\pi )\frac{N}{2}ln(\sigma )\frac{1}{2\sigma ^2}\underset{i=1}{\overset{i=N}{}}(xf_A)^2>\frac{N}{2}ln(2\pi )\frac{N}{2}ln(\sigma )\frac{1}{2\sigma ^2}\underset{i=1}{\overset{i=N}{}}(xf_B)^2$$
(7)
where $`N`$ is the number of observations, $`f_a`$ and $`f_b`$ are the time varying forecasts, and $`x`$ is the time-varying observations.

Cancelling terms from both sides:

$$\frac{1}{2\sigma ^2}\underset{i=1}{\overset{i=N}{}}(xf_A)^2>\frac{1}{2\sigma ^2}\underset{i=1}{\overset{i=N}{}}(xf_B)^2$$
(8)
Cancelling more terms this gives:

$$\underset{i=1}{\overset{i=N}{}}(xf_A)^2<\underset{i=1}{\overset{i=N}{}}(xf_B)^2$$
(9)
or

$$\text{MSE}_A<\text{MSE}_B$$
(10)
and so we see that comparing these forecasts using likelihood or MSE gives the same results i.e. that forecast A is better than forecast B.

### 4.2 Relationship between the likelihood and ignorance

Roulston and Smith (2002) describe a score for the assessment of probabilistic forecasts that they call the *ignorance*, and justify its usage on the basis of information theory and use in an optimal betting strategy. They define the ignorance for a single forecast-observation pair as minus the log (base 2) of the probability of the observation given the probabilistic forecast. We see that this is equivalent to minus log (base 2) of the likelihood for that single forecast-observation pair.

Comparing forecasts using the ignorance or any of the likelihood-based scores described above will give the same results if the forecasts errors are uncorrelated in time. If the errors are correlated in time, and this is taken into account in the calculation of the likelihood, then they may give differing results.

One can consider the likelihood as a generalisation of the ignorance to a) forecasts with autocorrelated forecast errors and b) forecasts for many variables, locations or leads at once. One can consider the ignorance as a special case of the likelihood when forecast errors are taken to be uncorrelated, and when looking at only a single variable, location and lead.

## 5 Summary

We have summarised the use of the likelihood for the evaluation of the skill of probabilistic forecasts. We believe that likelihood provides a useful general framework for the calibration and evaluation of all probabilistic forecasts, for all variables. We are in the process of applying the likelihood to various forecasting situations that are relevant to our business: examples are given in Jewson et al. (2003a) and Jewson et al. (2003b).

A number of question arise that merit further investigation. These include:

* When calibrating forecasts to maximise the likelihood, what numerical methods can be used to estimate the forecast error covariance matrix?
* Is it really necessary to calculate the likelihood using the correct forecast error covariance matrix, or is it satisfactory in practice to make the assumption that forecast errors are uncorrelated? One can argue that if the covariance matrix is not correctly modelled, then forecasts with autocorrelated errors are given more credit than is their due. However, it may be that in practice the ranking of forecasts is the same whether or not the covariance is estimated accurately.
* What are the relationships, if any, between the likelihood and other skill scores apart from those discussed above?

## 6 Acknowledgements

The author has had helpful discussions on the topics discussed above with a number of people, including Anders Brix, Pablo Doblas-Reyes, Renate Hagedorn and Christine Ziehmann.

## 7 Version number

This is version 2. Version 1 (August 12, 2003) contained some incorrect statements concerning the relationship between the likelihood and the Brier score. Sorry for any confusion this may have caused.


text:
# NIR surface photometry of a sample of nearby spiral galaxies

## 1 Introduction

Surface photometry of external galaxies in the optical domain has burgeoned in recent decades with advent of large format CCD arrays and the popularity of the 2D optical cameras. Data have been amassed since then for a wide range of astrophysical studies, from morphological classification (van der Kruit & Searle 1981, 1982) to stellar population characterization (de Jong 1994, 1996a, 1996b); from star formation research (Kennicutt 1983, 1989) to the investigation of the ISM in external galaxies (Valentijn 1994). The situation is much poorer in other wavelength regimes. The near infrared (NIR) domain is particularly interesting in this respect since in this range the flux is largely dominated by the direct stellar radiation while drastically reducing the general extinction. In addition, NIR colours are very sensitive to population changes in the galaxy (Peletier et al. 1994; de Jong & van der Kruit 1994; Giovanardi & Hunt 1996; Hunt et al. 1997; Peletier & Balcells 1996, 1997).

Reasons for the scarcity of NIR extragalactic databases are to be sought in the late development of 2D detector arrays with respect to their optical counterparts and in the greater difficulty in building NIR astronomical instruments, which have to be cryogenically cooled to achieve good performances in the thermal infrared, namely from the $`K`$ band towards longer wavelengths. Furthermore, observations and data reduction are far more complicated and time consuming in the NIR than in the optical. The basic reason is that the signals have to be extracted from a largely dominant sky emission, coming mostly from the thermal radiation of the earth atmosphere, which in addition is much more time-dependent, with variation timescales ranging from several tens of minutes in the $`J`$ band, to one or two minutes in $`K`$. This means that to properly subtract the sky emission from the image of an extended object, like a galaxy, one has to nod between two adjacent sky positions with a frequency adjusted to the variation timescale of the sky background. Also, the high value of the sky flux necessarily shortens the single on-chip integration time to avoid saturation and/or to keep the counts within the linear range of the detector. The combined result of all these effects is that an observing run for an extended object ends up with hundreds of individual exposures that have to be treated carefully. It is also worth mentioning that this procedure of sky substraction permits the galaxy object to completely cover the detector, since the sky will be measured in the other nodding position. Thus, the criterion of size in the sample can be extended to the full limit imposed by the instrument’s field of view.

We have recently started an observational programme aimed at producing a database of NIR infrared images of spirals galaxies, with several primary objectives in mind:

* To characterize the presence of the relevant morphological components and compare, whenever possible, their structural parameters with those obtained in the optical.
* To check their morphological classification.
* To investigate trends in structural parameters with morphological type.

In this first paper we will not cover the last two items, the sample does not contain a wide range of galaxy types. In subsequent reports we will, try to enlarge the sample and investigate the underlying stellar population in order to delineate the presence of interstellar dust by means of NIR colours.

The galaxies presented in this paper have been selected primarily according to their isophotal $`D_{25}`$ diameter, as listed in the RC3 (de Vaucouleurs 1991), and their brightness, both being considered within the suitable range for testing the feasibility of the project. We included only spirals in this first sample since they exhibit a wider range of morphological structures.

## 2 Observations

The observations were made at the Teide Observatory (OT) on the island of Tenerife (Spain) on the 1.5 m Carlos Sánchez Telescope (CST). The images were taken with the common user NIR camera, CAIN, equipped with a NICMOS3 $`256^2`$ detector array and standard broad band NIR filters $`JHK_\mathrm{s}K`$, cooled to LN2 temperature. CAIN holds two different optical set-ups, both fully cryogenic, which provide two plate scales: the narrow camera, giving 0.4 arcsec per pixel, and the wide one, with 1 arcsec per pixel, used for this work. All the galaxies have been observed in at least two filters, $`J`$ and $`K_\mathrm{s}`$ (or $`K`$), during several campaigns from 1998 onwards. Typical seeing value was, on average, a little in excess of 1 arcsec.

As mentioned in the preceding section, this first sample was selected taking into account several prescriptions:

* Declinations between $`35`$ and +60 degrees, to be observable from the OT.
* Morphological types with T in the range $`(1,9)`$. All the galaxies are spirals (barred and non-barred) and nearly face on (inclination $`<60^{}`$) for easier decomposition of their morphological structures. The majority of the galaxies in the sample have $`T`$ in the range 4–5, see table 1.
* High surface brightness, $`B<13`$, and diameters suitable for the telescope size and the instrument FOV, respectively.

The observed sample till now contains more than the ten galaxies, NGC3344, NGC3686, NGC3938, NGC3953, NGC4254, NGC4303, NGC4314, NGC5248, NGC6384 and NGC7479, which constitute the body of this paper (table 1). The remaining ones will be presented in subsequent papers.

The observing strategy has been in all cases the classical ON–OFF method, in which alternative exposures of the galaxy and the adjacent sky are taken until completing the total exposure time. The integration time of a single frame was determined by the linear part of the well depth of the NICMOS detector and the sky flux, which is by far the dominant contribution to the measured signal. Typical integration times per frame range from 6 to 12 s in $`K_\mathrm{s}`$ and from 20 to 30 s in $`J`$. Several frames are taken per pointing position until the grand total equals the timescale of the variation of the sky background, which was set at two minutes with a sufficient safety margin. This is then the time between nodding positions (table 2). The final images showing the surface photometry of each galaxy are displayed in Figs. 1, 2, 3, 4 and 5.

## 3 Data reduction and calibration

The data have been reduced using several combinations of IRAF tasks. The read-out mode has in all cases been the standard Fowler mode, in which N non-destructive reads are taken after the detector reset and then integrated for the required time and the detector read $`N`$ times again. By so doing the determination of the pedestal level of the signal is given by the average of the first $`N`$ reads, and the image signal is defined by the difference between the average of the last $`N`$ reads and that of the first ones. The S/N ratio of the image is improved by this method .

We have followed the standard process for treating the data to correct for bias and dark current, and flatfield calibration. Then we subtracted the sky from each individual frame using a coadded image of the sky measurements taken immediately before and after. Hence, we can correct for trends observed in the background. Finally, the corrected object frames are aligned and averaged to give the final product per galaxy and filter in adu. These are transformed to flux-calibrated images using observations of several NIR standard stars selected from the UKIRT list (Hunt et al. 1998) observed during the night.

In these calibrated coadded images we can estimate the sensitivity limits defined as the signal equivalent to 3 or 5 times the standard deviation. The results are given in table 3.

## 4 Components Decomposition

One of our main purposes is to investigate the presence and geometrical parameters of the principal structural components that contribute to the observed surface brightness. Radial profiles are the appropriate tool for this purpose and can be easily constructed with the technique of ellipse fitting. This is done by means of the IRAF task ELLIPSE, which computes the isophotal contours of the surface brightness images and fits them to ellipses. ELLIPSE uses an algorithm developed by Jedrzejewski (1987). As a result, one gets radial profiles of brightness, ellipticity and position angles along the major axis of the ellipse family. This is formally equivalent to first correcting for the inclination angle by deprojection of the galaxy image and then getting those radial variations averaged over the azimuthal angle. We have limited the scope of these profiles to the extent where the brightness of the galaxies equals the sensitivity figures given in table 3. The profiles, then, will provide information about the morphological structure of the galaxies.

Outer isophotes were used for determination of global galactic parameters such as ellipticity, position angle (hereafter, PA) and inclination angle, this latter derived from the ellipticity. Galaxies with a lower inclination angle have higher errors because of the lack of a privileged direction.

One can easily infer the presence of the different structural components, discussed in the next section, from the variation in these parameters with galactic radius (Varela et al. 1996; Wozniak et al. 1995).

In table 4 we have the results of the ELLIPSE task. The ellipticity is defined as $`b/a`$, where $`a`$ and $`b`$ are the major and minor axes of the ellipses, respectively. The PA is the angle between the major axis of the ellipses and north–south axis in the sky, measured from north to east. The values of each galaxy in both filters are in reasonably good agreement, their deviations being due mainly to errors in the ELLIPSE process, possibly enhanced by the slightly different sensitivities of the images in the two filters.

We can compare our results for the inclination angle and PA of each galaxy listed in table 4 with the values quoted in RC3. The agreement is generally good. NGC 4314 exhibits the largest differences certainly because of its large size in relation with the instrument FOV. We will not extend further the comparison between the NIR images presented here and the optical images in the RC3, since the component descomposition is not performed in the optical images. The reasons for this are to be found in the non-uniformity of the optical images and the expected spurious effects due to observational circunstances beyond our control. However, it is clear that the NIR images more closely represent the true stellar distribution caused by the smaller value of the extintion. It is then expected that the geometrical features of the main structural components will be somewhat different in the two spectral regimes.

Using the results of the ellipse fitting, we have measured the integrated magnitudes of the sample in both filters (see table 4) within a circular aperture of the radius at which the surface brightness from the object equals that of the sky at $`3\sigma `$ over its mean value (see table 3). These values are in fairly good agreement with those from the 2MASS data (Skrutskie et al. 1995), taken from NED (see table 4). However, our images extend in general to larger radii than the fixed scale of 80 arcsec quoted in NED.

In those galaxies where a clear indication exists of the presence of a large central bar, either from the images themselves and/or from the sudden change in the radial profiles. We used an additional method to estimate the geometrical parameters of the bar itself. The global PA and inclination angle of the object galaxy are used to deproject the galaxy image in order to measure quantities in the proper galactic plane.. Columns 6 and 7 of table 4 show the values of ellipticity and PA for the bar taken from the ELLIPSE results, after deprojection, and measured at the galactic radius spanned by the bar.

Once the radial profiles have been obtained and the general parameters calculated, we can now proceed with the decomposition of the averaged brightness profile into different structural components, each with a separate contribution to the observed flux. There are a large number of similar researches in the literature where structural decomposition of external galaxies are attempted in various ways (Simien & Michard 1990; Prieto et al. 1990, 1992; Byun & Freeman 1995), but most of them are based on visible data that do not so closely resemble the true stellar distribution because of their higher value of the extinction. Previous NIR decomposition can be found in de Jong (1996), Peletier & Balcells (1997), Moriondo et al. (1998), Seigar & James (1998) and references therein. For the different component identification we have followed the technical approach of Prieto et al. (2001) where the changes in the radial profiles of ellipticity, PA and B4 are used to identified the radial interval over which the given component extends.

Then, the analytical functions, described in the table 5 and in Aguerri (1998), are fitted to the brightness radial profile within that interval. This technique would be of little use in the case of strong bars since these structures only occur over a narrow range of azimuthal angles. Hence, the use of azimuthally averaged profiles, as is the case for the results of the ellipse fitting, could well mask the bar against the averaged background. We will come back to these cases at the end of this section.

Our procedure to fit the several morphological components to the measured brightness profile is similar to that described in Prieto et al. (2001). Once the components have been identified in the radial profiles of several parameters coming from the ellipse fitting, we first fix the disc in the outer parts of the galaxy, where contamination from other components can be neglected. Next, this disc model is subtracted from the original profile and we fit the bulge to the residuals with the best Sersic law ($`n=1`$–4). These operations are iterated until convergence, defined as the difference between two consecutive set of parameters being less than the measured noise, is reached.

After the disc–bulge pair has been set up, we proceed with the secondary components, by fitting the residuals of the subtraction of the disc–bulge model from the brightness profile to the spatial extent indicated in the radial profiles of the B4, PA and ellipticity parameters. In galaxies with prominents bars, such as NGC 7479 or NGC 4314, we have tried two different radial profiles for the fitting. In addition of using the azimuthally averaged brightness profiles resulting from the ellipse fitting, we have also built another two averaged brightness profile by computing the radial profiles by averaging in a narrower range of azimuthal angles: either profiles along the bar direction or perpendicular to the bar. A disc–bulge pair is now fitted to the profiles averaged along the minor axis of the bar where the bar contribution is low. Then, this model is subtracted from the radial profile averaged along the major axis of the bar direction and the residuals matched with a bar model. As before, this procedure is iterated until convergence, defined as in the previous case, is reached. In the following sections we will refer to the decomposition results from the ELLIPSE fitting as isophotal parameters and as fitting parameters to those of this last technique. In the case of the fitting parameters for the bar, we can also fit a bar profile along both the major and minor axes of the bar once the disc–bulge pair is fixed.

## 5 Results

The numerical results of our component decomposition are summarized in tables 7 and 8. Historically, the standard function for approximating galactic bulges has been the $`r^{1/4}`$ law, but it has been suggested that some bulges might be better represented by an $`r^{1/n}`$ law (Caon et al. 1993; Andredakis et al. 1995). When $`n`$ increases, the luminosity profiles are more centrally concentrated. For this reason, early-type galaxies show larger values of $`n`$ than late types (Andredakis et al. 1995). In this respect, observational effects can affect the results of the fit. For instance, the seeing size, if comparable with the galaxy diameter, can distort the concentration of the surface brightness (Trujillo et al. 2001). This should not be the case in our sample since all the galaxy bulges are much bigger than the seeing. In our fits, the best results for the whole sample is obtained with a value of $`n`$ equal to 1, but we cannot correlate $`n`$ with the morphological type because we have only Sc or Sbc galaxies (late-type galaxies) in this first sample of objects.

## 6 Conclusion

Most of the galaxies in the sample are of late type (bc or c), NGC 4314 being the only example of an early-type (SBa) object. Hence, we cannot infer any correlations between morphological type and the component parameters. This will be done in subsequent works with a larger sample of morphological types. The majority of the galaxies in the sample have $`T`$ in the range 4–5, see table 1. Only NGC 4314 and NGC 6384 are earlier galaxies with $`T`$ 1.0 and 3.6 respectively.

The bulge/disc ratio ($`B/D`$), defined in equation 1 (Mihos et al. 1994),

$$B/D=\frac{1}{0.28}\frac{r_e}{h}\frac{\mathrm{\Sigma }_e}{\mathrm{\Sigma }_o},$$
(1)
is rather low (see table 6), showing that the inner parts of the galaxy are dominated by components other than the central bulge, most probably the bar and ring.

It is worth noting, however, that NGC 4314 has the highest B/D ratio, with a value around of 0.6, indicating its early type. This galaxy has a spike in the ellipticity and PA profiles around 10′′. This peak points to the presence of an additional bar in the inner parts of the galaxy (Aguerri 1998). The existence of this secondary bar in NGC 4314 has been also proposed by Erwin & Sparke (2002) analysing images from the Hubble Space telescope (see their fig. 2b). The secondary bar has been fitted to the residuals after removal of the main components in the averaged brightness radial profile.

The resulting inclinations and PAs are in agreement with the data of other authors, in particular with the RC3. The major differences are between those galaxies with low PAs, in which the measurements are more difficult because of the poorer definition of this parameter. This is equally true for the inclination angles. There are minor discrepancies between the numbers in each filter, with some departures in NGC 3938 because of its shorter exposure time in $`K_\mathrm{s}`$. On the other hand, NGC 4303 and NGC 4314 show the greatest difference in inclination compared with the RC3 values in both filters. NGC 4303 has outer arms that are fainter in the NIR than in the optical bands, so we might be underestimating the outer part of the galaxy (on the other hand, we are comparing visible and NIR data). The same reason applies to NGC 4314. In general, all the galaxy images are bigger in the optical band, the images being extracted from the DSS, than in this work. But the particular orientation of the outer arms and the central bar causes these larger departures in these two galaxies. Some galaxies present a ring in RC3 (NGC 3686, NGC3953, NGC 5248 and NGC 6384), in NIR their are very faint in the images and could be mixed with the arms.

The two methods that we have used to fit the parameters of the brightness distribution for the morphological components, namely the isophotal parameters and the fitting parameters, give very similar results, as can be seen in Fig. 6. This can be attributed to the small influence of the bar in the overall brightness distribution outside the very central regions. In addition, part of the observed differences are certainly due to innaccuracies inherent in the fitting procedures. We can see that galaxies with strong bars (such as NGC 7479 and NGC 4314) have their isophotal bar parameters overestimated. This happens because of distortions in the ellipse fitting due to the presence of the spiral arms in the interarm region near the bar (NGC 6384, NGC 3953 and NGC 3686; see Figs. 1, 2 and 5).

We have obtained the best results with $`n=1`$ in the Sersic law for all the bulge fittings. Prieto et al. (2001) have argued that $`n`$ increases with the wavelength of the filters. We cannot discuss this effect since our filters are rather close each other in wavelength. Only NGC 4314 (the early-type galaxy) presents a possible $`n=23`$ class bulge, but in the ellipticity profile there is a peak around $`r=10^{\prime \prime }`$, which indicates a secondary bar (Erwin & Sparke 2002; Aguerri 1998). After the removal of this component, the profile of the bulge is then back of type $`n=1`$. The uniqueness of this result may be questioned, since we would have obtained a similarly good agreement with the observed profile using a bigger and triaxial bulge instead of a secondary bar. However Erwin & Sparke (2002), using deep images taken with the Hubble Space Telescope, have also found these inner structures. Some authors (Andredakis et al. 1995; Moriondo et al. 1998) have reported a correlation between the exponent n in the Sersic law with the morphological type of the galaxy, which we cannot check here because of the lack of coverage of this parameter.

It is noteworthy that NGC 4314 and NGC 7479 present the strongest bars in the sample, and that both have a LINER in their nuclei. This fact is discussed by Martin (1995), Athanasoula (1982) and Friedli & Berz (1993). Finally, all the bars found in this work are best fitted assuming that they are flat-type bars that follow the equation in table 5.

###### Acknowledgements.

We are thankful to Drs J.A. López Aguerri and A. Cabrera for many helpful and useful conversations during the preparation of this manuscript. To Norberto Castro-Rodríguez for his helpful in the fitting of the profiles. The CST is operated on the island of Tenerife by the Instituto de Astrofísica de Canarias at the Spanish Observatorio del Teide of the Instituto de Astrofísica de Canarias. This article makes use of data products from 2MASS, which is a joint project of the Univ. of Massachusetts and the Infrared Processing and Analysis Center, funded by the NASA and the NSF. This work has been partially funded by the Spanish Plan Nacional de Astronomía y Astrofísica, project AYA2000-2046-C02-02 and the Universidad de La Laguna-Cajacanarias fellowship (Canary Islands, Spain).


text:
# GRB 021219: the first Gamma-Ray Burst localized in real time with IBAS Based on observations with INTEGRAL, an ESA project with instruments and science data centre funded by ESA member states (especially the PI countries: Denmark, France, Germany, Italy, Switzerland, Spain), Czech Republic and Poland, and with the participation of Russia and the USA.

## 1 Introduction

Although considerable progress in the understanding of Gamma-Ray Bursts (GRBs) has been made after the discovery of their afterglows at longer wavelengths (see, e.g., van Paradijs, Kouveliotou & Wijers (2000)) the debate on their progenitors is still open. A model involving the core collapse of a massive star (e.g. Woosley (1993)) seems to be favored for long GRBs and evidence for a Supernova-GRB association has been found recently (Stanek et al. (2003); Hjorth et al. (2003)). The situation is more uncertain for what concerns short ($`<`$ 2 s) bursts, mostly because of the lack of counterparts at other wavelengths for this class of GRBs.

A multi-wavelength approach is crucial to the understanding of the complex nature of GRBs. The short duration of the prompt $`\gamma `$-ray emission and the fading character of the afterglow impose a rapid follow-up. This can be achieved only if the positions derived from the prompt emission are immediately distributed to the scientific community. The INTEGRAL satellite, although not built as a GRB-oriented mission, can contribute to this task thanks to the INTEGRAL Burst Alert System (IBAS; Mereghetti et al. 2003b ). To date IBAS has successfully detected and located 5 GRBs, with increasing accuracy and decreasing delays (4.36 error radius after 30 s for GRB 030501; IBAS Alert 596).

GRB 021219 is the first GRB detected in real time by IBAS. This occurred only two months after the launch of INTEGRAL, while the satellite was still in its performance and verification (PV) phase. During the PV phase, the external distribution of the IBAS alerts was not enabled yet. An internal alert message, produced only 10 s after the start of the burst, reached the members of the IBAS Localization Team and a GRB Coordinates Network (GCN) circular with a preliminary position could be issued 5 hours later, after a quick analysis to confirm the GRB (Mereghetti et al. (2002)). Since during the PV phase the relative alignment between IBIS and the satellite star trackers was not well measured yet, an uncertainty of 20 was attributed to the derived position. A refined position could be derived later (Götz et al. (2002)), thanks to the presence of Cyg X-1 in the field of view.

Here we report on the results on GRB 021219 obtained with the IBIS (Ubertini et al. (2003)) and SPI (Vedrenne et al. (2003)) instruments. They are both coded mask imaging telescopes with a large field of view (29$`\times `$29 IBIS, 36 diameter SPI). IBIS is based on two detectors, ISGRI (Lebrun et al. (2003)) and PICsIT (Di Cocco et al. (2003)), operating in the 15 keV - 1 MeV and 175 keV - 10 MeV energy ranges, respectively. SPI consists of 19 Ge detectors cooled at 85 K and works in the 20 keV - 8 MeV range. The GRB was located outside the field of view of the two low energy monitoring instruments JEM-X (Lund et al. (2003)) and OMC (Mas-Hesse et al. (2003)).

## 2 Data Analysis

### 2.1 IBIS

During this observation PICsIT provided images integrated on time intervals of several minutes, which cannot be used for the study of GRBs. Owing to the telemetry limitations at satellite level this is the standard operation mode of PICsIT. Useful data for GRB studies can be obtained by PICsIT when it is operated in photon-by-photon mode (as it occurred for GRB 021125, see Malaguti et al. (2003)). Our results for GRB 021219 are therefore based only on data from the ISGRI detector.

We have analyzed ISGRI single interaction events, for which the arrival time, energy deposit and interaction pixel are known.

GRB 021219 has been detected at off-axis angles Z=1.38 and Y=10.02, in the partially coded field of view of IBIS, with a signal-to-noise ratio of 15.5 in the 15-100 keV band. The top panel of Fig. 2 shows the GRB light curve binned at 0.2 s intervals. It refers to the 15-500 keV band, and, to increase the signal to noise ratio, it has been extracted using only the pixels illuminated by the GRB for at least half of their surface. The GRB started at 07:33:55 UTC and had a T90 duration of 5.5 s. The time profile shows a single peak.

To derive an accurate position we selected a 15 s time interval starting at 07:33:55 UTC, obtaining $`\alpha _{J2000}`$ = 18h 50m 25s, $`\delta _{J2000}`$ = +31 56 23′′, with an uncertainty of 2 radius. As shown in Fig. 1, this position is consistent with the one published earlier and with the annulus derived independently with the International Planetary Network (IPN, Hurley et al. (2002)).

We have extracted the GRB peak spectrum (over a time interval of 1 s) using 128 linearly spaced bins between 19 keV and 1 MeV and then rebinned it to have at least 25 counts per channel. Since IBIS is a coded mask instrument, the background spectrum can be measured simultaneously to that of the target, by using the pixels not illuminated by the source. A fully calibrated response matrix for large off-axis angles is not available yet. Therefore, we divided the counts spectrum by the closest (in detector coordinates) spectrum of the Crab Nebula. The resulting photon spectrum (Fig. 3) is well described by a power law with photon index $`\mathrm{\Gamma }`$=1.6$`\pm `$0.1 and flux in the 20-200 keV band of 3.7 photons (3.5$`\times `$10-7 erg) cm-2 s-1.

The total (time averaged) spectrum has been obtained with the same technique and can be fitted by a power law with $`\mathrm{\Gamma }`$=2.0$`\pm `$0.1. The fluence in the 20-200 keV band (corrected for the telemetry gaps) is 9$`\times `$10-7 erg cm-2.

We have also investigated the spectral evolution with time. Four spectra of the duration of $``$1 s each have been extracted: one during the rising part of the GRB, one at the peak, and two during the decaying tail. All the spectra are well represented by a power law model. The corresponding photon indices, plotted in the lowest panel of Fig. 2, provide a significant evidence for a softening with time of the burst spectrum.

### 2.2 SPI

At the time of GRB 021219 the SPI instrument was in low-telemetry mode. Therefore, no single events have been transmitted to the ground in photon-by-photon mode. Only the events which occurred in several detectors (multiple events) and those analysed by the on-board pulse discriminator (PSD events) were kept. Due to this, there was no sensitivity below $`200\mathrm{keV}`$. Still SPI transmitted to the ground the countrates of all the events measured in each of the 19 Ge detectors in bins of 1 s (see middle panel of Fig. 2). For these events the information about the photon energy is lost and thus only a broad band count rate can be given. The background was determined using the data starting from 20 minutes before the burst occurrence. The peak flux and fluence have been then derived from the broad band count rate by comparison with the count rate measured for the Crab Nebula. This results in a peak flux of $`10.7\pm 1.1\mathrm{Crab}`$ at 07:33:56 UTC, which is consistent with the value obtained with IBIS/ISGRI. The values are consistent also for what concerns the fluence.

The imaging analysis of the SPI data, using a 3 s time interval starting at 07:33:55 UTC, yields a position $`\alpha _{J2000}`$ = 18h 49.3m, $`\delta _{J2000}`$ = +31 46.1, with an uncertainty of 30 (S/N = 8.5) localizing the GRB 17.6 from the IBIS/ISGRI position.

## 3 Discussion

GRB 021219 had a rather steep average spectrum, indicating that it can be considered an X-ray rich GRB. Usually GRB spectra are well described by a phenomenological model (Band et al. (1993)) consisting of two power laws, one at low energies (with slope $`\alpha `$) and one at high energies (with slope $`\beta `$), and a smooth break typically between about 100 and 400 keV (Preece et al. (2000)). The fact that no break is seen in the spectrum of GRB 021219 up to 150 keV indicates that the break, if any, is outside our detection range. The value of the photon index $``$2, compared to the ones typically observed in the BATSE sample $`\beta 2.5`$, suggests that we are observing the part of the spectrum above the break. We cannot exclude, however, that the break is at higher energies: this would qualify GRB 021219 even more as an X-ray rich GRB.

The time resolved spectroscopy performed with IBIS/ISGRI indicates a clear hard-to-soft spectral evolution. This is a common feature in many GRBs observed with previous satellites (e.g. Norris et al. (1986), Ford et al. (1995), Frontera et al. (2000)).

No optical counterpart for this GRB has been reported, with limiting magnitudes of R=13.7 at t-t0$``$7.5 hours (Lipunov et al. (2002)), R=18 at t-t0$``$11 hours, I=19.5 at t-t0$``$18 hours (Henden et al. (2002)) and R=20.5 at t-t0$``$34 hours (Castro-Tirado et al. (2002)). Observations in the radio band at 4.86 GHz (at t-t0$``$16 hours) have not detected any new source down to a 4$`\sigma `$ limit of 220 $`\mu `$Jy (Berger et al. (2002)).

The fact no optical transient has been detected indicates a ”dark” or at least dim nature of this GRB. A similar behaviour has been noticed before in other X-ray rich events. For example, the X-ray rich GRB 021211 (Crew et al. (2003)) was fainter than R$``$22.5 12 hours after the burst (Klose et al. (2002)). The optical time history of this GRB is compatible with the upper limits reported above for GRB 021219.

## 4 Conclusions

The successful detection and localization of GRB 021219 only one month after IBIS activation, showed that IBAS is able to derive and distribute the position of the GRBs detected in the field of view of IBIS within a few tens of seconds.

The analysis of GRB021219 reported here shows that IBIS/ISGRI is indicated for detailed spectral studies. In fact we have been able to obtain spectra integrated over only 1 s with good statistics, thus deriving the spectral shape of the GRB and its variation with time.

###### Acknowledgements.

This research has been partially supported by the Italian Space Agency. LF acknowledges the hospitality of the ISDC during part of this work.


text:
# The Spiros imaging software for the Integral SPI spectrometer

## 1 Introduction

Astronomical images are often of fields consisting of (or at least dominated by) a number of sources which are essentially points. Their angular size is much smaller than the resolution of current instrumentation \- frequently by a factor more than $`10^{10}`$. At gamma-ray energies the number of such sources is usually small. In these circumstances the model with a minimum number of parameters that is consistent with the data will consist of a list of the positions of those sources with their intensities. For an instrument with spectroscopic capabilities like INTEGRAL/SPI, a description of the variation of those intensities with photon energy, and perhaps with time, must be added. In general it is when a minimum number of parameters is sought that each of these may be obtained with the highest precision111We note that in a certain sense, this is a ‘minimum information/maximum entropy’ interpretation of the data - it can be stored or transmitted with a small number of bits and hence very high informational entropy as defined by Shannon Shannon (1966), so in these circumstances point source searching and fitting is the preferred data analysis technique.

## 2 Context and environment

The technique of ‘Iterative Removal of Sources’ (IROS, Hammersley et al. hammersly (1984)) has been widely used for coded mask instruments. A simple image of the field of view is made using a mapping technique which is optimised for finding a source assuming that the data can be explained by only that source, plus background. The mapping gives the approximately location and intensity of the source, which are then improved by maximising a measure of the goodness of fit. The residuals of the fit are used as the input for a further image reconstruction and source search. The parameters of the two source model are refitted and if the fit represents a significant improvement on the original one the process is continued with more and more sources.

SPIROS is a programme which implements this algorithm for the SPI spectrometer of INTEGRAL, a coded mask imaging instrument with a detector array comprising 19 high-purity Germanium detectors giving an angular resolution of about 2.5 over a field of view of 16 (Vedrenne, et al. instrument\_paper (2003)). SPIROS operates within the ISDC software environment (Courvoisier et al. isdc (2003)).

## 3 Optimisation criteria and background modelling

The SPI instrument differs from most coded mask instrument for which IROS has been used because

1. The number of detector elements is very limited – 19 in the basic mode, though the concept of ‘pseudo-detectors’ can increases this by up to 104, with double and triple events and pulse shape discrimination (Vedrenne et al. instrument\_paper (2003)).
2. In consequence multiple telescope pointings are analysed simultaneously in order to obtain enough information for a unique solution.

The most widely used coded mask image reconstruction techniques, based on correlations of the detector plane with a representation of the form of the mask shadow, usually by FFT, are not appropriate in such circumstances. The limited number of measurements means that effects which average out in the many-pixel limit cannot be ignored. Edge effects and differences between detectors, in their response for different directions, and in their noise levels take on a dominant role. However because of the limited number of measurements, one can use matrix techniques that would be impracticable in other cases.

SPIROS operates on data which have been already binned by (pseudo-) detector, by pointing, and by energy (strictly ‘pulse height’) bin. In most modes it treats data from one energy bin at a time, although a group of input energy bins can be combined to form one. It also reads files containing attitude information, integration times, etc, and a file containing one or more background models.

Considering, then, a particular energy, for a given set of assumed source positions, the expected count in detector (or pseudo-detector) $`d`$ during pointing $`p`$ is

$$\widehat{P}_{dp}=t_{dp}\{\underset{i}{}S_iA_{idp}+B_{dp}\}.$$
(1)
Here $`A_{idp}`$ is the effective area of the detector to a source in the direction of the $`i`$th source, which is assumed to have flux $`S_i`$ within the energy bin considered. The $`A_{idp}`$ values are obtained by interpolation of data which are stored in so-called IRF (Instrument Response Function) files generated by Monte Carlo simulations. The effective exposure time, corrected for dead-time and all other effects, is $`t_{dp}`$ and the $`A_{idp}`$ roll together all the effects of the mask coding and the detector efficiency. For simplicity, only the photopeak response and a single energy are considered here.

The possibility that the background is a combination of components which vary in different ways between pointings and from detector to detector can be included in this formalism by treating those components much like sources:

$$\widehat{P}_{dp}=t_{dp}\{\underset{i}{}S_iA_{idp}+\underset{i^{}}{}F_i^{}B_{i^{}dp}\}.$$
(2)
Here provision is made for the possibility that the form of the background variation with time may be known but that it needs to be multiplied by an unknown factor $`F_i^{}`$, analogous to a source intensity (see Section 4).

The objective is then to find the model (the combination of $`S_i`$ and, if required, $`F_i^{}`$) which, best explains the observed data (the count rates $`P_{dp}`$) in the “Maximum Likelihood” (ML) sense.

In the general case where the number of counts per bin may be small, as can be the case when short exposures or narrow energy bins are being considered, then the relevant statistic is (Cash cash79 (1979)) :

$$C=2\underset{dp}{}\{\widehat{P}_{dp}P_{dp}Ln(\widehat{P}_{dp})\}.$$
(3)
Even in the case of fixed source positions, finding the $`S`$ and $`F`$ which optimises this statistic is a non-linear problem and has to be solved iteratively.

If the counts are large, Gaussian statistics can be assumed and one can minimise the $`\chi ^2`$ statistic

$$\chi ^2=\underset{dp}{}\frac{(P_{dp}\widehat{P}_{dp})^2}{\sigma _{dp}^2}.$$
(4)
If the source positions are known (or supposed) and if in addition one makes the approximation that $`\sigma _{dp}^2=P_{dp}`$, instead of $`\sigma _{dp}^2=\widehat{P}_{dp}`$ then the problem of finding the $`S_i`$ becomes a linear one.

SPIROS can be run either using the general $`ML`$ approach or using $`\chi ^2`$ and the assumption that $`\sigma _{dp}^2=P_{dp}`$. The latter is faster and more efficient in those cases where it can validly be used. Note that the user must beware of the dangers of adopting the $`\chi ^2`$ statistic where it is not valid because for low counts per bin there is a significant probability that $`P_{dp}=0`$. The infinite weights that would result are avoided by ignoring such data, but this introduces a bias, as does the incorrect weight given to other bins with low, but non-zero counts.

## 4 Background handling

Usually the background has to be treated as an unknown, but some constraining assumptions are necessary - with $`n_p`$ pointings and $`n_d`$ detectors, it is obviously not possible at the same time to obtain information about sources in the field of view *and* make independent estimates of the background in each of the $`n_dn_p`$ combinations of detector $`d`$ and pointing $`p`$.

SPIROS reads in one or more ($`n_i^{}`$, in general) background models which are sets of $`B_{dp}`$ values222strictly $`t_{dp}B_{dp}`$ is stored, generated for a specific data set by an independent programme (called ‘spiback’). Sometimes the model components may be absolute (when they are based on data preceding or following an observation, or from energy bands just above and below the region of interest, for example). Sometimes they are simply tracers of a time variation with arbitrary scaling. The simplest case would be background which is uniform and constant, corresponding to a single $`i^{}`$ component in which all the $`B`$ have the same value. Other possibilities are components which depend on time or from detector to detector in a specific way, or which follow tracers of expected background contributions (for example the veto shield count rate is a measure of the particle flux, the rate of out of range events in the germanium detectors is a measure of high energy particles…).

In its existing form SPIROS has provides for the following possibilities :

1. The backround estimates may be taken to be absolute estimates, not requiring the fitting of a factor $`F`$
2. They may be taken to reflect the time (pointing to pointing) variation of the background but with a potentially different scaling factor $`F`$ per detector, to be found by fitting within SPIROS.
3. It may be assumed that relative, detector-detector, backgrounds have been provided, as well as the correct time ($`p`$) dependence, so that SPIROS need only fit a single normalising $`F`$.
4. In a slight generalisation of (3), one factor may be fitted for all single detectors but a different factor for pseudo-detectors corresponding to double events in particular pairs of detectors (or triples, etc.)

The most general case (2) is usually adopted. Further options may be added in the future, such as assuming that the relative backgrounds in different detectors are known but with a different unknown scaling factor to be fitted per pointing or on some other timescale.

## 5 Mapping; finding and positioning sources

A basic mapping operation consists of considering each pixel in the image successively, placing a test source at that position and establishing the intensity that it would have in order to best match the observed data, along with the uncertainty in that intensity. Sources found in previous iterations, or read in from a catalogue, are either subtracted out from the data ($`\chi ^2`$) or taken into account in the analysis ($`ML`$). For source searching the intensities and uncertainties evaluated on a comparatively coarse grid (e.g. 0.5) can be used. A smoothed linear interpolation is then used to fill in estimates on a finer pitch.

The source selected for potential addition to the list of sources is that which has the highest value of intensity/uncertainty.

Before accepting a new source as real, a simultaneous optimisation of its position and reoptimisation of all the other sources that do not have good catalogued positions is performed (section 6 below). The procedure used is an iterative one with a descent along the line of maximum slope.

## 6 Relationship of IROS to other methods

Numerous deconvolution techniques have been proposed and used for image reconstruction in astronomy. Examples are Maximum Entropy, and Richardson-Lucy; see Starck & Pantin (decon-review (2002)) for a recent review. Images obtained by such methods can be considered as vectors in an N-dimensional space, where the component of the vector along each of the axes represents the intensity in one of the N pixels. In the presence of noise, there will always be a certain volume in this space, inside which the points correspond to images that are consistent with the data. Different methods apply constraints which are not the same and so do not lead to the same choice of a point in hyperspace among all of the possible ones.

The IROS algorithm searches the space for a solution consistent with the data according to the following rules:

Starting from point in this space which corresponds to the current solution, test each axis (pixel number) to find the point along that axis with the lowest $`\chi ^2`$.

Identify which pixel gave the maximum descent in the sense corresponding to a positive source flux. Step along the axis corresponding to that pixel to the lowest point; all other components of the vector are left at their original values.

Intensity fitting corresponds to optimising in a subspace which has just those dimensions corresponding to the sources already found. The position fitting is more complex. An optimisation is made in a separate space with $`3S`$ dimensions, corresponding to the $`2`$ position coordinate and the intensity of each of $`S`$ sources. The $`S`$ directions corresponding to the optimised positions of the sources are then added to the original space to make an $`N+S`$ dimensional one.

Formally, there is no guarantee that this algorithm leads to a minimal description of the data in terms of point sources. Pragmatically, in simulations and in observations of involving known sources, it is found to be effective.

The procedure is very similar to the CLEAN method used in radio astronomy (Högbom, hogbom (1974); Schwarz, schwarz (1978)), except that in CLEAN (i) position optimisation is not normally performed, (ii) the point in the map with the highest absolute value is chosen, so iterations can add negative components to the image, (iii) only a fraction of the intensity of a source is subtracted. The resulting image will usually have a relatively large number of non-zero pixels, whereas with IROS, there is just one per source. CLEAN can be considered as an exploration of the multi-dimensional space considered above, with no positivity constraint and without the fitting stage.

The IROS algorithm differs from many image reconstruction techniques in that the possible source positions are not restricted to a fixed grid of pixel positions.

## 7 Spectral extraction

In a different mode, SPIROS is used for the extraction of spectral information. This requires an input catalogue, which may simply be a list of known sources or which may be the result of a previous run of SPIROS in imaging mode. In each energy bin in turn, the combination of source intensities (and, if required, of background parameters) which best matches the data is found. The method is as described in section 3.

In this way one obtains for each source a spectrum analogous to a ‘Pulse Height’ spectrum in that off-diagonal terms in the energy response matrix are not taken into account. However other aspects of the instrument response (detector photopeak efficiency variation with energy, for example) have been corrected for. For sources with conventional continuum spectra, this measure is already a very good estimate of the input spectrum. However for definitive results, a programme such as XSPEC (Arnaud xspecref (1996)) needs to be used to take into account the off-diagonal terms in the response of the combination, SPI+SPIROS. XSPEC-compatible response matrices for this step have been derived using Monte Carlo simulations of observations of monoenergetic sources at 100 different energies.

## 8 Other modes

Optionally, the sources may be treated as having a finite extent. For example it is possible to treat each source as a Gaussian function and find the width which best matches the data.

Extracting a light curve (intensity as a function of time) for each of the sources in an input catalogue is directly analogous to the extraction of pulse height spectra.

Although SPIROS is intended for fitting of point sources, it does have a mode in which one solves simultaneously for the intensities of the fluxes in each pixel of an image, allowing map of diffuse emission to be generated. Such inverse problems are notoriously unstable if the number of pixels is high and the coding is not ideal. Thus instead of a simply multiplying by the inverse of the coding matrix, it is modified by the addition of a diagonal Wiener term or some other smoothing constraint matrix, giving a stabilising effect by allowing the diagonal terms to dominate.

## 9 Examples of results obtained with SPIROS and Conclusions

Results obtained with Spiros can be seen in other articles in this issue and elsewhere - see, for example, Bouchet et al. (bouchet (2003)).

As an example of a slightly non-standard use of Spiros, Figure 1 shows an image obtained using *only* events interacting in 2 or 3 detectors. For such events there is no information about which of the 2 or 3 detectors corresponds to the first interaction. The source is well identified and located, though with significance lower than in the corresponding single event image (Figure 2).

The form in which the IROS algorithm is implemented within SPIROS handles for the first time the situation where no simplifying assumptions can made about position independence of the recorded mask shadow. By using the instrument response characterised by a generalised matrix and an object-oriented programme structure, additional sophistication can be introduced as necessary and as the knowledge of the instrument improves.

###### Acknowledgements.

The development of SPIROS was supported by a grant from PPARC. The work has benefited from inputs from many members of the Integral SPI Data Analysis Group (ISDAG) and the Integral Science Data Centre (ISDC).