text:
# The host galaxy and optical light curve of the gamma-ray burst GRB 980703Based on observations made with the NASA/ESA Hubble Space Telescope (HST), obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.,Based on observations made with the Nordic Optical Telescope (NOT), operated on the island of La Palma jointly by Denmark, Finland, Iceland, Norway, and Sweden.,Based on observations made at the 2.2-m telescope of the German-Spanish Calar Alto Observatory

## 1 Introduction

GRB~980703 (BATSE trigger 6891) was detected on 1998 July 3.182 UT by the All-Sky Monitor on the *Rossi $`X`$-Ray Timing Explorer* satellite (Levine et al. LMM1998 (1998)). Frail et al. (FHB1998 (1998)) and Zapatero Osorio et al. (ZOC1998 (1998)) identified a fading optical afterglow (OA), and Djorgovski et al. (1998b ) derived a redshift of $`z=0.9661\pm 0.0001`$ from absorption and emission lines in the combined spectrum of the OA and host galaxy. They also derived a star-formation rate of $`10_{\mathrm{}}\mathrm{yr}^1`$ for the host galaxy from the strength of the \[O ii\] emission line. They stressed that the star formation rate depends on the adopted extinction and on the star-formation indicator used, but concluded that the host is a compact starburst galaxy. $`R`$-band optical light curves of the OA $`+`$ host galaxy were presented by Bloom et al. (BFK1998 (1998)), Castro-Tirado et al. (CTZ1999 (1999)), and Vreeswijk et al. (VGO1999 (1999)), who all found that the optical light decayed as a single power law with a slope of $`\alpha 1.4`$, and all estimated the magnitude of the host galaxy to be $`R_{\mathrm{gal}}22.6`$. Sokolov et al. (SFK2001 (2001)) used late-time photometry of the OA $`+`$ host to estimate global star-formation rates in the host galaxy of $`15`$$`66_{\mathrm{}}\mathrm{yr}^1`$, depending on the details of the initial mass function.

The host galaxy of GRB~980703 appears to have a high extinction along the line of sight to the gamma-ray burst (GRB). Castro-Tirado et al. (CTZ1999 (1999)) found $`A_V=2.2`$ from the slope of $`X`$-ray/optical/infrared fluxes 0.9 days after the burst, Vreeswijk et al. (VGO1999 (1999)) used a similar method to find $`A_V=1.50\pm 0.11`$ 1.2 days after the burst, Djorgovski et al. (1998a ) found $`A_V0.3\pm 0.3`$ from the Balmer decrement 4.4 days after the burst, and Bloom et al. (BFK1998 (1998)) estimated $`A_V=0.9\pm 0.2`$ 5.3 days after the burst. The different (and formally inconsistent) values of $`A_V`$ may be due to different methods for measuring the extinction probing different regions of the host galaxy, or may indicate that the environment where the burst occurred is changing with time.

We have obtained new ground-based and *Hubble Space Telescope* (*HST*) Space Telescope Imaging Spectrograph (STIS) images of the host galaxy of GRB~980703 taken 17, 577, 710, and 716 days after the burst. These data are discussed in Sect. 2. The nature of the host galaxy and the location of the GRB within the galaxy are discussed in Sect. 3, and the light curve of the OA is discussed in Sect. 4. Unless otherwise noted we have assumed a cosmology where $`H_0=100h_{100}`$ km s-1 Mpc-1, $`\mathrm{\Omega }_m=0.3`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. For this cosmology a redshift of 0.9661 corresponds to a luminosity distance of $`4.43h_{100}^1`$ Gpc and a distance modulus of $`43.235\mathrm{log}_{10}(h_{100})`$. One arcsecond corresponds to $`10.9h_{100}^1`$ comoving kpc or $`5.6h_{100}^1`$ proper kpc, and the lookback time is $`5.3h_{100}^1`$ Gyr.

## 2 Observations and Data Reductions

### 2.1 Ground-based observations

We obtained $`R`$\- and $`I`$-band images of the field containing GRB~980703 using the Andalucía Faint Object Spectrograph and Camera (ALFOSC) on the 2.56-meter Nordic Optical Telescope (NOT) at La Palma on 1998 July 20.172 (16.99 days after the burst). A log of all our observations is given in Table 1. The instrumental gain was $`1.0`$ e-1/ADU and the read-out noise was $`6`$ e-1/pixel. We calibrated the $`R`$-band data using six stars from Rhoads et al. (RDC1998 (1998)), and the $`I`$-band data using three stars from Vreeswijk et al. (VGO1999 (1999)). The magnitudes were measured in apertures with radii of 15 pixels ($`=2\stackrel{}{.}82`$). The derived relations between the instrumental and standard systems are $`R=r+(8.440\pm 0.011)+(0.074\pm 0.023)(ri)`$ and $`I=i+(7.934\pm 0.038)`$. The RMS of the residuals are 0.07 and 0.09 mag respectively. We found combined host$`+`$OA magnitudes of $`R=22.40\pm 0.07`$ and $`I=21.71\pm 0.09`$ in an aperture of radius $`1\stackrel{}{.}524`$ ($`=8.5h_{100}^1`$ proper kpc) centred on the OA. The full-width at half-maxima (FWHMs) of the stellar profiles on the combined NOT image were $`0\stackrel{}{.}8`$ in the $`R`$ band, and $`0\stackrel{}{.}9`$ in the $`I`$ band.

In order to determine the magnitude of the host galaxy without contamination from the OA we obtained $`V`$\- and $`R`$-band images of the host galaxy with the Calar Alto Faint Object Spectrograph (CAFOS) at the Centro Astronómico Hispano-Alemán (CAHA) Observatory’s 2.2-meter telescope on 2000 January 5 (551 days after the burst). These data have a gain of $`2.3`$ e-1/ADU, a read-out noise of $`5.06`$ e-1/pixel, and were calibrated using the same field stars and apertures as for the NOT data. The derived relations between the instrumental and standard systems are $`V=v+(5.943\pm 0.051)+(0.079\pm 0.072)(vr)`$ and $`R=r+(6.145\pm 0.021)+(0.245\pm 0.030)(vr)`$. The RMS of the residuals in the calibrations are 0.031 and 0.027 mag in $`V`$ and $`R`$ respectively. For the host galaxy we measured $`V_{\mathrm{gal}}=23.00\pm 0.10`$ and $`R_{\mathrm{gal}}=22.57\pm 0.09`$ in an aperture of radius $`1\stackrel{}{.}524`$. The CAHA combined $`V`$-band image had FWHM $`=1\stackrel{}{.}5`$ while the combined $`R`$-band image had FWHM $`=1\stackrel{}{.}7`$.

GRB~980703 is located at $`\alpha =`$ 23:59:06.72 $`\delta =+`$08:35:07.3 (J2000) which correspond to Galactic coordinates of $`\mathrm{}^{\mathrm{II}}=101\stackrel{}{.}56`$, $`b^{\mathrm{II}}=52\stackrel{}{.}10`$, hence the Galactic extinction is small. We have adopted a Galactic reddening of $`E_{BV}=0.058\pm 0.020`$ (corresponding to extinctions of $`A_V=0.19\pm 0.07`$, $`A_R=0.16\pm 0.05`$, and $`A_I=0.11\pm 0.04`$) from the DIRBE/IRAF dust maps of Schlegel et al. (SFD1998 (1998)). The dereddened magnitudes of the host are $`V_{0,\mathrm{gal}}=22.81\pm 0.12`$ and $`R_{0,\mathrm{gal}}=22.41\pm 0.10`$.

If we assume that the host galaxy has a power law spectrum of the form $`f\nu ^\beta `$ then we derive $`\beta _{\mathrm{gal}}=1.08\pm 0.59`$, which gives $`I_{0,\mathrm{gal}}=21.91\pm 0.23`$.

### 2.2 The HST data

We used the *HST*/STIS to obtain CCD images of the host galaxy of GRB~980703 as part of the Cycle 9 proposal GO-8640. The total exposure times were 5178 seconds in the 50CCD (clear; hereafter referred to as CL) aperture and 5264 seconds in the F28X50LP (long pass; hereafter referred to as LP) aperture. The CL images were taken on 2000 June 12 (710 days after the burst) and the LP images were taken on 2000 June 18 (716 days after the burst). The CCD gain was set to $`1`$ e-/ADU and the read-out noise was $`4.46`$ e-/pixel. We used a four-point STIS-CCD-BOX dithering pattern with shifts of 2.5 pixels ($`=0\stackrel{}{.}127`$) between exposures. The data was pre-processed through the standard STIS pipeline and combined using the dither (v1.2) software (Fruchter & Hook FH2001 (2001)) as implemented in IRAF111Image Reduction and Analysis Facility (IRAF), a software system distributed by the National Optical Astronomy Observatories (NOAO).(v2.11.3)/STSDAS(v2.1.1). The “pixfrac” parameter was set to $`0.6`$, and we selected a final output scale of $`0\stackrel{}{.}0254`$/pixel. These observations are part of the Survey of the Host Galaxies of GRBs (http://www.ifa.au.dk/~hst/grb\_hosts/index.html, Holland et al. 2000a ).

We measured the total AB magnitude of the host galaxy on the CL and LP images in an aperture of radius $`1\stackrel{}{.}524`$. The zero points of Gardner et al. (GBB2000 (2000)) yield AB magnitudes of $`22.87\pm 0.02`$ in the CL aperture and $`22.50\pm 0.03`$ in the LP aperture.

## 3 The Host Galaxy

### 3.1 The morphology of the host galaxy

Fig. 1a shows a contour plot of a section of the STIS CL image containing the host galaxy. The galaxy is extremely compact and has a somewhat irregular, “egg-shaped” morphology. In order to quantitatively characterize the morphology of the galaxy we fitted elliptical isophotes to the CL and LP images of the host using Jedrzejewski’s (J1987 (1987)) algorithm as implemented in the STSDAS task isophote.ellipse. These isophotes provided a resonable fit to the galaxy, but there are significant residuals as the isophotes are not elliptical. The mean observed ellipticity of the host on the CL image is $`ϵ=0.16\pm 0.01`$ for $`0\stackrel{}{.}05r0\stackrel{}{.}51`$, and the mean position angle of the semi-major axis on the CL image is $`\theta _0=12\stackrel{}{.}7\pm 0\stackrel{}{.}2`$ east of north on the sky. We find $`ϵ=0.19\pm 0.01`$ and $`\theta _0=9\stackrel{}{.}2\pm 0\stackrel{}{.}2`$ on the LP image. The variation of the ellipticity and position angle with radius is shown in Fig. 2. The host galaxy shows some evidence for twisted isophotes as the position angle of the major axis rotates from $`+45\mathrm{°}`$ to $`45\mathrm{°}`$ between $`0\stackrel{}{.}1`$ and $`0\stackrel{}{.}5`$ from the centre of the host galaxy. This change is seen in both the CL and LP images. These values have not been corrected for the effects of the point-spread function (PSF) on the the observed shape of the galaxy. Schweizer (S1979 (1979)) and Holland (H98 (1998)) have shown that the observed shape and surface brightness profile of a galaxy can be dominated by the effects of the PSF if the diameter of the galaxy is less than $`8`$ times the width of the PSF. The FWHM of the STIS PSF in the CL and LP images is $`0\stackrel{}{.}084`$, so the shape of the PSF dominates the central regions of the host galaxy, and will effect the observed shape in the outer regions. The PSF has a mean ellipticity of $`0.11`$ and a mean position angle of $`31\stackrel{}{.}5`$ in the CL image. This is $`44\mathrm{°}`$ from the observed major axis of the host galaxy, so the PSF will act to make the host appear to be less elliptical than it actually is.

In order to determine the intrinsic (i.e., before the light passed through the *HST*/STIS optics) surface-brightness distribution of the host galaxy we fit two-dimensional PSF-convolved generalized de Vaucouleurs profiles (Sersic S1968 (1968))

$$S(R)=S_e10^{b(n)/\mathrm{ln}(10)\left[(R/R_e)^{1/n}1\right]}$$
(1)
where $`R_e`$ is the radius enclosing half the light from the galaxy and $`S_e`$ is the surface brightness at $`R_e`$, to the CL and LP images of the host galaxy using the software described in Holland et al. (HCH1999 (1999)). This software allows all of the model parameters to be fit simultaneously using a chi-square minimization scheme. The parameters of the best fitting models, along with the $`\chi ^2`$ per degree of freedom (DOF) are listed in Table 2. The number of degrees of freedom in each fit is the number of pixels used in the fit minus the number of model parameters. The constant $`b(n)`$ is chosen so that $`R_e`$ encloses half of the light of the galaxy. The mean surface brightness inside $`R_e`$ is denoted by $`\mu _e`$. The formal 1-$`\sigma `$ errorbars for each parameter do not include contributions from correlations between the parameters, and may be affected by the systematic residuals caused by the presence of substructure in the host galaxy.

The standard de Vaucouleurs $`R^{1/4}`$ profile does not provide a good fit to the host galaxy, whereas a generalized de Vaucouleurs profile ($`R^{1/n}`$) with $`n1`$ provides a significantly better fit. However, both fits are formally strongly excluded by the data, which is due to the fact that there are large residuals in the central few pixels. Fig. 1b shows the CL image of the host galaxy of GRB~980703 with the best-fitting two-dimensional PSF-convolved model subtracted. The model was computed to a radius of 64 pixels ($`=1\stackrel{}{.}626=9.1h_{100}^1`$ proper kpc) from the centre of the galaxy. The model provides a good fit for $`R0\stackrel{}{.}25`$, but there are systematic residuals in the inner $`0\stackrel{}{.}25`$ of the galaxy. However, we note that the model assumes that the galaxy has a smooth surface brightness distribution and does not include small-scale structure, such as star-forming regions, within the galaxy. Therefore, we do not expect to obtain $`\chi ^2/\mathrm{DOF}1`$ for our fits. The best-fitting model provides a good fit over most of the host galaxy except for the central few pixels, and the observed residuals are consistent with the presence of substructure in the central regions of the galaxy.

A Sersic (S1968 (1968)) model with $`n=1`$ corresponds to an exponential disk. Therefore, our modelling suggests that the host galaxy of GRB~980703 has a disk-like structure, with smaller sub-structure in the central regions. The excess light seen on the west side of the galaxy after subtracting the best-fitting model (see Fig. 1b), may be a faint spiral arm. The fitted half-light radii of our Sersic (S1968 (1968)) models correspond to disk scale radii of $`0\stackrel{}{.}196\pm 0\stackrel{}{.}002`$ on the CL image and $`0\stackrel{}{.}225\pm 0\stackrel{}{.}002`$ in LP.

The derived half-light radius ($`0.72h_{100}^1`$ proper kpc) is much smaller than the half-light radii seen in local late and early type galaxies (Im et al. ICG1995 (1995)). Phillips et al. (PGG1997 (1997)) and Guzmán et al. (GGK1997 (1997)) studied galaxies in the Hubble Deep Field North that were *selected* to be compact (i.e., had measured half-light radii smaller than $`0\stackrel{}{.}5`$ including the effect of the WFPC2 PSF). Their galaxies covered a range of redshifts from $`z0.4`$ to $`z1.0`$ and appear to have properties very similar to the host galaxy of GRB~980703, namely $`VI`$ colours in the range $`0.7<VI<1.1`$, half-light radii in the range $`0.65h_{100}^1R_e2.6h_{100}^1`$ kpc (converted to our cosmology and not corrected for the effect of the WFPC2 PSF), and observed $`I`$-band magnitudes in the range $`21I24`$ for the galaxies at redshifts similar to the redshift of GRB~980703. The spectra of these galaxies are similar to the spectrum of the host galaxy—i.e., dominated by strong emission lines from O ii, and O iii, as well as Balmer lines (Guzmán et al. GGK1997 (1997); Djorgovski et al. 1998a ). The rest-frame equivalent width (EW) of the O ii line from the host galaxy was found by Djorgovski et al. (1998a ) to be 46 Å, which is in the the middle of the range $`5EW94`$ Å found by Guzmán et al. (GGK1997 (1997)) for the compact galaxies in the Hubble Deep Field North.

### 3.2 The colour profile of the host galaxy

The host galaxy of GRB~980703 has $`V_0=22.81\pm 0.12`$ and $`(VR)_0=0.40\pm 0.15`$ measured in an aperture with a radius of $`1\stackrel{}{.}524`$ on the CAHA images. This colour is consistent with the colours predicted by Bloom et al. (BFK1998 (1998)) and Vreeswijk et al. (VGO1999 (1999)) based on fitting a single power law plus a constant flux to the optical decay of the OA. We find no evidence for azimuthal variations in colour, but the center of the galaxy is $`0.2`$ mag bluer than the outer regions of the galaxy with the step in the colour occurs at approximately the effective radius of the best-fitting model (see Sect. 3.1). Matthews & Gallacher (MG1997 (1997)) note that this sort of trend in colour is unique to late-type disc galaxies and some dwarf galaxies. Fig. 3 shows the surface brightness profiles and the azimuthally averaged colour gradient that were derived from the STIS CL and LP images.

### 3.3 Star formation in the host galaxy

At a redshift of $`z=0.9661`$ a rest-frame wavelength of 2800 Å corresponds to the observer’s $`V`$ band. Therefore, we are able to use Eq. (2) of Madau et al. (MPD1998 (1998)) to estimate the star formation rate in the host galaxy from the continuum flux at 2800 Å. For a Salpeter (S1955 (1955)) initial mass function we find a star-formation rate of $`(4.1\pm 0.5)h_{100}^2`$ $`_{\mathrm{}}\mathrm{yr}^1`$ while a Scalo (S1986 (1986)) initial mass function yields a star-formation rate of $`(6.4\pm 0.7)h_{100}^2`$ $`_{\mathrm{}}\mathrm{yr}^1`$. If we take $`h_{100}=0.7`$ (de Bernardis et al. dBA2000 (2000)) then the star-formation rate is $`8`$$`13`$ $`_{\mathrm{}}`$ yr-1, depending on the assumed initial mass function. However, this calculation assumes that there is no extinction in the host galaxy, which is probably an incorrect assumption. If $`A_V=0.3`$ in the host then the star-formation rate becomes $`10`$$`20`$ $`_{\mathrm{}}\mathrm{yr}^1`$, while if if $`A_V=2.2`$ in the host then the star-formation rate is $`65`$$`100`$ $`_{\mathrm{}}\mathrm{yr}^1`$. These values are consistent with the star-formation rate obtained by Djorgovski et al. (1998a ) using the strength of the \[O ii\] emission line ($`10`$ $`_{\mathrm{}}\mathrm{yr}^1`$ with $`A_V=0.3\pm 0.3`$).

### 3.4 The luminosity of the host galaxy

Lilly et al. (LTM1995 (1995)) find a typical magnitude of $`(M_B^{})_{AB}=21.4`$ for blue galaxies at $`0.75z<1.00`$, if $`(H_0,\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda })=(50,1,0)`$. For our adopted cosmology, this corresponds to $`(M_B^{})_{\mathrm{AB}}=21.47+5\mathrm{log}_{10}(h_{100})`$. At $`z=0.9661`$ the observed Kron–Cousins $`I`$ band is approximately half-way between the rest-frame Johnson $`U`$ and $`B`$ bands. If we assume that the host galaxy has a power-law spectrum with $`\beta _{\mathrm{gal}}=1.08\pm 0.59`$ (from the CAHA data, see Sect. 2.1) then $`L_B(3.2\pm 0.7)h_{100}^2L_B^{}`$ where $`L_B^{}`$ is the rest-frame $`B`$-band luminosity of a typical blue galaxy at $`z=0.9961`$. Most of the uncertainty in the luminosity comes from the uncertainty in the shape of the host galaxy’s rest-frame spectrum at $`\lambda 4000`$ Å. For $`h_{100}=0.7`$ the rest-frame $`B`$-band luminosity is $`(1.6\pm 0.4)L_B^{}`$. Therefore, the host is approximately an $`L^{}`$ galaxy with a total luminosity similar to that of other galaxies at $`z=0.9661`$. Our comparison of the total $`B`$-band luminosity of the host to $`L_B^{}`$ is, however, somewhat uncertain because $`M^{}`$ is highly correlated with the slope of the faint end of the galaxy luminosity function, and with its normalisation, $`\varphi ^{}`$ (Lilly et al. LTM1995 (1995)).

The specific star-formation rate per unit blue luminosity of the host galaxy of GRB~980703, if we ignore extinction within the host, is $`6.5_{\mathrm{}}\mathrm{yr}^1L_{B}^{}{}_{}{}^{1}`$, which is similar to that of other GRB host galaxies. The specific star-formation rate is $`0.6`$ times that of the hosts of GRB~970508 ($`11_{\mathrm{}}\mathrm{yr}^1L_{B}^{}{}_{}{}^{1}`$) and GRB~990123 ($`11_{\mathrm{}}\mathrm{yr}^1L_{B}^{}{}_{}{}^{1}`$) and $`1.5`$ times that of the host of GRB~990712 ($`4.4_{\mathrm{}}\mathrm{yr}^1L_{B}^{}{}_{}{}^{1}`$). If the extinction in the host is non-negligible then the true star-formation rate per unit blue luminosity may be much higher.

There is a faint ($`\mathrm{CL}=27.14\pm 0.05`$, $`\mathrm{LP}=26.29\pm 0.06`$) object $`2\mathrm{}`$ ($`=11.2h_{100}^1`$ proper kpc) south of the host galaxy (see Fig. 1a), which may be a companion galaxy. This object has FWHM $`=0\stackrel{}{.}097`$, which is slightly larger than the resolution of the CL image ($`0\stackrel{}{.}084`$) and appears to have a faint extended structure around it.

### 3.5 The position of GRB 980703 in the host galaxy

The paucity of stars in the STIS images makes it difficult to accurately determine the location of GRB~980703 within the host galaxy. Only four stars are present in both the STIS CL image and our ground-based images, and two of these are very faint in our ground based images. We used these four stars to determined the location of the OA relative to the galaxy by comparing the STIS CL image to an $`R`$-band image obtained with the NOT on 1998 July 7.163 (Pedersen et al. PKJ1998 (1998)). At the epoch of the NOT observation, the OA contributed $`25`$% of the total host$`+`$OA flux, so we registered, scaled, and smoothed the STIS image to match the seeing and intensity of the NOT image, then subtracted the STIS image from the NOT image. The OA is found to be offset $`0\stackrel{}{.}23`$ west and $`0\stackrel{}{.}02`$ north of the center of the host galaxy with an estimated 1$`\sigma `$ error of approximately $`\pm 0\stackrel{}{.}10`$. This is within approximately 1$`\sigma `$ of the location determined by Bloom et al. (BKD2001 (2001)). We find the GRB to have occurred $`0\stackrel{}{.}23\pm 0\stackrel{}{.}14`$ from the centre of the host at a position angle of $`85\stackrel{}{.}0\pm 4\stackrel{}{.}6`$, while Bloom et al. (BKD2001 (2001)) find that the GRB occurred at a distance of $`0\stackrel{}{.}11\pm 0.09`$ with a position angle of $`28\stackrel{}{.}9\pm 0\stackrel{}{.}3`$. The $`1\sigma `$ error ellipses of these two position determinations overlap, but we note that Bloom et al. (BKD2001 (2001)) used deep Keck images to obtain an astrometric solution involving 23 objects while our solution involved only four stars.

## 4 The Late-Time Light-Curve of the Optical Afterglow

Figs. 4 and 5 shows the $`R`$\- and $`I`$-band light curves for the OA based on measurements taken from the literature and on our own NOT measurements. The magnitudes have been corrected for the contribution from the underlying host galaxy and for Galactic extinction. The light from the host galaxy was removed by subtracting the total light from the host galaxy as determined from our CAHA images.

We fit both a single power law ($`f_\nu (t)=f_\nu (0)t^\alpha `$), where $`f_\nu (t)`$ is the flux in $`\mu `$Jy $`t`$ days after the burst, and a broken power law

$$f_\nu (t)=\{\begin{array}{ccc}f_\nu (t_b)(t/t_b)^{\alpha _1},\hfill & \mathrm{if}\hfill & tt_b\hfill \\ f_\nu (t_b)(t/t_b)^{\alpha _2},\hfill & \mathrm{if}\hfill & t>t_b,\hfill \end{array}$$
(2)
where $`\alpha _1`$ is the slope before the break, $`\alpha _2`$ is the slope after the break, and $`t_b`$ is the time in days after the burst that the break occurred, to the $`R`$-band data as described in Holland et al. (HBH2000 (2000)) and Jensen et al. (JFG2001 (2001)). Only data $`t10`$ days after the burst were used for the fit since a Type Ic supernova (SN) like SN1998bw would not make a significant contribution to the total flux at these times. The parameters of our best fits are listed in Table 3 and the fits themselves are shown in Figs. 4 and 5. The small $`\chi ^2/\mathrm{DOF}`$ values for both the single and broken power law fits reflect the large formal uncertainties in the late time data that arise from subtracting the OA from the host at a time when both the OA and the host had similar magnitudes. The small number of data points makes it impossible to fit a broken power law in the $`I`$-band, so we scaled the broken power law that was fit to the $`R`$-band data to match (in a chi-square sense) the early $`I`$-band photometry. The best fit occurred for $`I_0(t_b)=20.77\pm 0.05`$. The broken power law is formally a slightly better fit to the $`R`$-band light curve than a single power law is, although an F-test suggests that the broken power law is only significant at about the $`1\sigma `$ level. There are several gaps in the decay curve, and large error bars in the late-time data due to the subtracting of the light from the bright host galaxy, so we cannot unambiguously conclude that there is a break in the light curve. Therefore, we prefer the single power law fit since it requires fewer free parameters to fit the light curve. We note, however, that breaks have been seen in most well-sampled GRB light curves, and that the time and magnitude of the break for GRB~980703 is consistent with the breaks seen for several other OAs.

There is growing evidence that at least some GRBs are associated with supernovae. The prototype for this association is the temporal and positional coincidence of GRB~980425 and the unusual Type Ib/c SN1998bw (Galama et al. GVP1998 (1998)). VLT and *HST*/STIS images show that this GRB/SN occurred in a star-forming region in the sub-luminous barred spiral galaxy ESO~184$-$G82 (Sollerman et al. SKF2000 (2000); Fynbo et al. FHA2000 (2000)). Evidence for a GRB/SN association has also been seen for GRB~980326 (Castro-Tirado & Gorosabel CTG1999 (1999); Bloom et al BKD1999 (1999)), GRB~970228 (Dar D1999 (1999); Reichart R1999 (1999); Galama et al. GTV2000 (2000)), GRB~970514 (Germany et al. GRS2000 (2000)), GRB~990712 (Fruchter et al. FVH2000 (2000); Björnsson et al. B2001 (2001)), and GRB~991002 (Terlevich et al. TFT1999 (1999)). In addition, Hudec et al. (HHH1999 (1999)) have identified five more candidates for SN/GRB associations based on positional and temporal coincidence. In Figs. 4 and 5 we show the predicted light curves for GRB~980703 if it was associated with a Type Ic SN like SN1998bw. Although both the $`R`$\- and $`I`$-band data are consistent with both a pure power law decline and a broken power law decline the addition of an SN Ic like SN1998bw results in a better fit to the data. However, the large error bars in the late-time data prevent us from being able to reliably distinguish between a SN and a no-SN scenario for GRB~980703.

It has recently been suggested that a dust echo could produce an excess flux in the decay of the late-time optical light curve of a GRB similar to the flux excesses seen in GRB~970228 and GRB~980326 (Esin & Blandford EB2000 (2000)). Using their method we predict an excess $`R`$-band flux due to dust 17 days after the burst of $`F_{\nu _{\mathrm{ob}}}^E0.03`$ $`\mu `$Jy, approximately $`2\sigma `$ less than the observed excess flux at this time, assuming a single power law decay. Therefore, we are unable to distinguish between the dust echo and no dust echo cases.

The $`R`$-band slope of the best-fitting single power law is $`\alpha =1.37\pm 0.14`$, which is consistent with other values in the literature. Combining the optical slope with the observed slope of the $`X`$-ray decay ($`\alpha _X<0.91`$, Vreeswijk et al. VGO1999 (1999)), and assuming that the GRB is a collimated outflow into a region with a $`\rho (r)r^\delta `$ density profile (e.g. Panaitescu et al. PMR1998 (1998), Mészáros et al. MRW1998 (1998)), implies that $`\delta =2`$ (Panaitescu & Kumar PK2001 (2001)). This suggests that the fireball is expanding into a local medium dominated by a pre-existing stellar wind, similar to the situation that Jaunsen et al. (JHB2001 (2001)) found for GRB~980519. The optical spectrum of the early OA (1.2 days after the burst) has a slope of $`\beta _{\mathrm{OA}}=2.71\pm 0.12`$ (Vreeswijk et al. VGO1999 (1999)). For an expansion into a pre-existing stellar wind the spectral index of the OA is given by $`\beta _{\mathrm{OA}}=(1+2\alpha )/3=0.58\pm 0.09`$, which is inconsistent with the observed early-time slope of OA’s optical spectrum. This discrepancy, however, can be resolved if there are several magnitudes of extinction along the line of sight to the GRB. The extinction measurements in Sect. 1 are consistent with $`A_V`$ decreasing by $`0.3`$ mag/day from $`A_V1.5`$ one day after the burst to $`A_V0.5`$ five days after the burst. An examination of the spectrum of Djorgovski et al. (1998b ) suggests that $`\beta _{\mathrm{OA}}1.6\pm 0.2`$ in the $`R`$ band 4.4 days after the burst, which suggests that the burst evolved to become bluer at later times, which is consistent with $`A_V`$ decreasing with time. This suggests that some of the discrepancy between the observed spectral index and the spectral index which is predicted from the decay of the optical and $`X`$-ray light curves may be due to the time evolution of the OA’s spectrum that result from changes in the amount of dust along the line of sight to the burst during the first few days after the burst.

## 5 Conclusions

The host galaxy of GRB~980703 is a compact blue, star-forming galaxy with $`V_{\mathrm{gal}}=23.00\pm 0.10`$ and $`(VR)_{\mathrm{gal}}=0.43\pm 0.13`$. The central regions of the host are $`0.2`$ mag bluer than the outer regions. Its rest-frame $`B`$-band luminosity is $`L_B3.2h_{100}^2L_B^{}`$, and it has a star-formation rate, based on the integrated flux at 2800 Å, of $`8`$$`13`$ $`_{\mathrm{}}\mathrm{yr}^1`$, assuming no extinction in the host. This suggests that the host galaxy is undergoing a phase of active star formation similar to what has been seen in other GRB host galaxies. We find that the galaxy is best fit by a Sersic (S1968 (1968)) $`R^{1/n}`$ profile with $`n1`$ and a half-light radius of $`0\stackrel{}{.}13`$ ($`=0.72h_{100}^1`$ proper kpc). This corresponds to an exponential disk with a scale radius of $`0\stackrel{}{.}22`$ ($`=1.21h_{100}^1`$ proper kpc). The residuals of the fit show that the morphology of the host is irregular, but we are unable to determine if the location of GRB~980703 coincided with any special features in the host. The galaxy resembles the compact star-forming galaxies that Phillips et al. (PGG1997 (1997)) and Guzmán et al. (GGK1997 (1997)) studied galaxies in the Hubble Deep Field North.

The optical decay of GRB~980703 followed the same broad pattern established by other long-duration GRBs: a power-law decay, perhaps with a break a few days after the burst, possibly followed by an increase in flux $`10(1+z)`$ days after the burst.

###### Acknowledgements.

This work was supported by the Danish Natural Science Research Council (SNF). JG acknowledges the receipt of a Marie Curie Research Grant from the European Commission. MA acknowledges the Astrophysics group of the Physics Department of University of Oulu for support of his work. GB acknowledges partial support from the Icelandic Research Council and the University of Iceland Research Fund. We would also like to thank the anonymous referee for a careful reading of the original manuscript, and for useful comments that greatly improved this paper.


format: text

text:
# A Measurement of the Branching Fraction for the Inclusive 𝐵→𝑋_𝑠𝛾 Decays with the Belle Detector

## Abstract

We have measured the branching fraction of the inclusive radiative $`B`$ meson decay $`BX_s\gamma `$ to be

$$\mathrm{𝐵𝑟}(BX_s\gamma )=(3.36\pm 0.53\text{(stat)}\pm 0.42\text{(sys)}{}_{0.54}{}^{+0.50}\text{(th)})\times 10^4.$$

The result is based on a sample of $`6.07\times 10^6`$ $`B\overline{B}`$ events collected at the $`\mathrm{{\rm Y}}(4S)`$ resonance with the Belle detector at the KEKB asymmetric $`e^+e^{}`$ storage ring.

PACS: 13.40.Hq, 14.40.Nd, 14.65.Fy

KEK preprint 2001-3  Belle preprint 2001-2

K. Abe10, K. Abe37, I. Adachi10, Byoung Sup Ahn15, H. Aihara38, M. Akatsu20, G. Alimonti9, K. Aoki10, K. Asai21, Y. Asano43, T. Aso42, V. Aulchenko2, T. Aushev13, A. M. Bakich34, E. Banas16, W. Bartel10,6, S. Behari10, P. K. Behera44, D. Beiline2, A. Bondar2, A. Bozek16, T. E. Browder9, B. C. K. Casey9, P. Chang24, Y. Chao24, B. G. Cheon33, S.-K. Choi8, Y. Choi33, Y. Doi10, J. Dragic18, A. Drutskoy13, S. Eidelman2, Y. Enari20, R. Enomoto10,11, F. Fang9, H. Fujii10, C. Fukunaga40, M. Fukushima11, A. Garmash2,10, A. Gordon18, K. Gotow45, H. Guler9, R. Guo22, J. Haba10, H. Hamasaki10, K. Hanagaki30, F. Handa37, K. Hara28, T. Hara28, T. Haruyama10, N. C. Hastings18, K. Hayashi10, H. Hayashii21, M. Hazumi28, E. M. Heenan18, Y. Higasino20, I. Higuchi37, T. Higuchi38, H. Hirano41, T. Hojo28, Y. Hoshi36, W.-S. Hou24, S.-C. Hsu24, H.-C. Huang24, S. Ichizawa39, Y. Igarashi10, T. Iijima10, H. Ikeda10, K. Ikeda21, K. Inami20, A. Ishikawa20, H. Ishino39, R. Itoh10, G. Iwai26, H. Iwasaki10, Y. Iwasaki10, D. J. Jackson28, P. Jalocha16, H. K. Jang32, M. Jones9, R. Kagan13, H. Kakuno39, J. Kaneko39, J. H. Kang46, J. S. Kang15, P. Kapusta16, N. Katayama10, H. Kawai3, N. Kawamura1, T. Kawasaki26, H. Kichimi10, D. W. Kim33, Heejong Kim46, H. J. Kim46, Hyunwoo Kim15, S. K. Kim32, K. Kinoshita5, K. Korotushenko30, P. Krokovny2, R. Kulasiri5, S. Kumar29, T. Kuniya31, E. Kurihara3, A. Kuzmin2, Y.-J. Kwon46, J. S. Lange7, M. H. Lee10, S. H. Lee32, H.B. Li12, D. Liventsev13, R.-S. Lu24, A. Manabe10, T. Matsubara38, S. Matsui20, S. Matsumoto4, T. Matsumoto20, Y. Mikami37, K. Misono20, K. Miyabayashi21, H. Miyake28, H. Miyata26, L. C. Moffitt18, G. R. Moloney18, G. F. Moorhead18, S. Mori43, A. Murakami31, T. Nagamine37, Y. Nagasaka19, Y. Nagashima28, T. Nakadaira38, E. Nakano27, M. Nakao10, J. W. Nam33, S. Narita37, Z. Natkaniec16, K. Neichi36, S. Nishida17, O. Nitoh41, S. Noguchi21, T. Nozaki10, S. Ogawa35, T. Ohshima20, Y. Ohshima39, T. Okabe20, T. Okazaki21, S. Okuno14, S. L. Olsen9, W. Ostrowicz16, H. Ozaki10, H. Palka16, C. S. Park32, C. W. Park15, H. Park15, L. S. Peak34, M. Peters9, L. E. Piilonen45, E. Prebys30, J. L. Rodriguez9, N. Root2, M. Rozanska16, K. Rybicki16, J. Ryuko28, H. Sagawa10, Y. Sakai10, H. Sakamoto17, M. Satapathy44, A. Satpathy10,5, S. Schrenk45,5, S. Semenov13, K. Senyo20, M. E. Sevior18, H. Shibuya35, B. Shwartz2, V. Sidorov2, J.B. Singh29, S. Stanič43, A. Sugi20, A. Sugiyama20, K. Sumisawa28, T. Sumiyoshi10, K. Suzuki3, S. Suzuki20, S. Y. Suzuki10, S. K. Swain9, H. Tajima38, T. Takahashi27, F. Takasaki10, M. Takita28, K. Tamai10, N. Tamura26, J. Tanaka38, M. Tanaka10, Y. Tanaka19, G. N. Taylor18, Y. Teramoto27, M. Tomoto20, T. Tomura38, S. N. Tovey18, K. Trabelsi9, T. Tsuboyama10, Y. Tsujita43, T. Tsukamoto10, S. Uehara10, K. Ueno24, N. Ujiie10, Y. Unno3, S. Uno10, Y. Ushiroda17, Y. Usov2, S. E. Vahsen30, G. Varner9, K. E. Varvell34, C. C. Wang24, C. H. Wang23, M.-Z. Wang24, T.J. Wang12,†, Y. Watanabe39, E. Won32, B. D. Yabsley10, Y. Yamada10, M. Yamaga37, A. Yamaguchi37, H. Yamamoto9, H. Yamaoka10, Y. Yamaoka10, Y. Yamashita25, M. Yamauchi10, S. Yanaka39, M. Yokoyama38, K. Yoshida20, Y. Yusa37, H. Yuta1, C.C. Zhang12, J. Zhang43, Y. Zheng9, V. Zhilich2, and D. Žontar43

In the framework of the Standard Model (SM), the inclusive decay $`BX_s\gamma `$, where $`X_s`$ is a hadronic recoil system containing an $`s`$ quark, proceeds primarily through the electroweak penguin diagrams of the radiative $`bs\gamma `$ transition. These diagrams are of interest because they have loops where non-SM particles, such as charged Higgs bosons or SUSY particles, can contribute and potentially produce deviations from SM expectations. An SM branching fraction calculation for $`BX_s\gamma `$ that includes next-to-leading order QCD corrections has an expected precision of 10% . Experimental measurements of the branching fraction at this level of precision are useful for identifying or limiting non-SM theories .

A semi-inclusive analysis is used to reconstruct the $`BX_s\gamma `$ decay from a primary photon, a kaon and multiple pions. The monochromatic photon energy from the two-body decay $`bs\gamma `$ is smeared by gluon emission, the Fermi motion of the $`b`$ quark inside the $`B`$ meson, and the $`B`$ meson momentum in the $`\mathrm{{\rm Y}}(4S)`$ center-of-mass (CM) frame.

The $`BX_s\gamma `$ decay is studied with the Belle detector at the KEKB asymmetric $`e^+e^{}`$ storage ring . The dataset consists of a sample of $`6.07\times 10^6`$ $`B\overline{B}`$ events corresponding to an integrated luminosity of $`5.8\text{ fb}^1`$ taken at the $`\mathrm{{\rm Y}}(4S)`$ resonance ($`\sqrt{s}=10.58\text{ GeV}`$), and an off-resonance background sample of $`0.6\text{ fb}^1`$ taken $`60\text{ MeV}`$ below $`\mathrm{{\rm Y}}(4S)`$. The beam energies at the $`\mathrm{{\rm Y}}(4S)`$ resonance are $`3.5\text{ GeV}`$ for positrons and $`8.0\text{ GeV}`$ for electrons. Quantities used in this analysis are defined in the CM frame, in which the beam energy $`E_{beam}`$ is defined as $`\sqrt{s}/2`$. Inclusion of charge conjugated modes are implied throughout the text.

A full description of the Belle detector is given in ; here we briefly describe the components relevant for this analysis. Charged tracks are reconstructed inside a 1.5 T magnetic field induced by a super-conducting solenoid magnet, by means of a 50 layer central drift chamber (CDC) that covers the laboratory polar angle $`17^{}<\theta _{lab}<150^{}`$. Tracks are fitted through the CDC and a three layer double sided silicon vertex detector (SVD), with a transverse momentum resolution of $`(\sigma _{p_t}/p_t)^2=(0.0019p_t)^2+(0.0034)^2`$ ($`p_t`$ in $`\text{GeV}/c`$). For particle identification, we combine information from three detectors. The high momentum range, typically from 1 to $`3.5\text{ GeV}/c`$, is covered by a silica aerogel Cherenkov counter (ACC) system, providing threshold-type kaon-to-pion separation. The momentum range below $`1.5\text{ GeV}/c`$ is covered by a time-of-flight (TOF) counter system, and the very low and high momentum ranges are accessible with $`dE/dx`$ information from the CDC. Photons are detected by an electromagnetic calorimeter (ECL) located between the particle identification devices and the solenoid coil. The entire tracking acceptance is covered with 8736 CsI(Tl) crystals that are typically $`5.5\times 5.5\text{ cm}^2`$ in cross-section at the front surface and 16.2 $`X_0`$ in depth. The photon energy resolution is $`(\sigma _E/E)^2=0.013^2+(0.0007/E)^2+(0.008/E^{1/4})^2`$ ($`E`$ in GeV), which is better than 2% for the primary photons used in this analysis. The Belle detector is modelled with the GEANT program  to simulate the detector response. Both the data and the Monte Carlo (MC) events  are reconstructed with the same program.

Hadronic events are selected with an efficiency of 99% for $`B\overline{B}`$ final states. Events with at least three tracks from the interaction region and visible energy greater than $`0.2\sqrt{s}`$ are selected. The QED and $`\tau ^+\tau ^{}`$ events that remain are removed using ECL information and event topology .

Photon candidates are selected from ECL clusters of $`5\times 5`$ crystals. Each photon candidate is required to have a laboratory energy greater than $`20\text{ MeV}`$ with no associated charged track, and shower shape consistent with an electromagnetic shower. We accept events in which the most energetic photon is in the barrel region ($`33^{}<\theta _{lab}<132^{}`$). We require that 95% of its energy is concentrated in the central $`3\times 3`$ crystals. The efficiency of these criteria for the primary photon is $`(80\pm 4)\%`$. In order to reject $`\pi ^0`$ and $`\eta `$ mesons, we combine this primary photon candidate with any other photon (for $`\eta `$ we only use photons with $`E_\gamma >200\text{ MeV}`$) and reject the event if the invariant mass of the two photons is within $`\pm 18\text{ MeV}/c^2`$ of $`M_{\pi ^0}`$ or $`\pm 32\text{ MeV}/c^2`$ of $`M_\eta `$.

We reconstruct the recoil system $`X_s`$ in 16 different final states of one charged kaon or $`K_S^0`$ plus one to four pions which may include one $`\pi ^0`$. We combine the primary photon with every $`X_s`$ combination and calculate the opening angle $`\theta _{X_s\gamma }`$ between $`\stackrel{}{p}_\gamma `$ and $`\stackrel{}{p}_{X_s}`$, the energy difference $`\mathrm{\Delta }E=E_\gamma +E_{X_s}E_{beam}`$ and the beam constrained mass $`M_{bc}=\sqrt{E_{beam}^2|\stackrel{}{p}_\gamma +\stackrel{}{p}_{X_s}|^2}`$. Then, we require $`\theta _{X_s\gamma }>167^{}`$, $`M_{bc}>5.2\text{ GeV}/c^2`$ and $`0.15<\mathrm{\Delta }E/\text{GeV}<0.1`$.

Charged kaon candidates are selected using a kaon-to-pion likelihood ratio. For every charged track, likelihoods from ACC, TOF and $`dE/dx`$ are calculated for pion and kaon hypotheses. A combined likelihood is constructed for each hypothesis, and a tight cut is applied on the likelihood ratio. Kaons are identified with a typical efficiency of $`65\%`$ with $`2\%`$ pion fake rate. All charged tracks that are not identified as kaons are considered to be pions.

The neutral kaon candidates are reconstructed in the $`K_S^0\pi ^+\pi ^{}`$ mode from two oppositely charged tracks. The $`K_S^0`$ momentum is recalculated with a vertex constrained fit. The candidates are then required to have an invariant mass within $`\pm 8\text{ MeV}/c^2`$ $`(2\sigma )`$ of $`M_{K^0}`$ and a vertex that is displaced from the interaction point in a direction consistent with the $`K_S^0`$ momentum. The $`K_S^0`$ reconstruction efficiency is $`(65\pm 4)\%`$.

The $`\pi ^0\gamma \gamma `$ candidates are reconstructed from two photons that have an invariant mass within $`\pm 16\text{ MeV}/c^2`$ $`(3\sigma )`$ of $`M_{\pi ^0}`$. Then, the $`\pi ^0`$ momentum is recalculated with a mass constrained fit. The $`\pi ^0`$ reconstruction efficiency is $`(60\pm 4)\%`$.

When multiple candidates are found in an event, the best candidate is selected as follows: If there is at least one charged kaon or pion that forms the $`X_s`$ vertex with the constraint from the run-by-run determined profile of the $`B`$ meson decay vertices , the candidate with the largest vertex confidence level is chosen. The candidate with the largest $`\theta _{X_s\gamma }`$ is chosen in the following cases: (1) there is no charged kaon or pion to form the vertex; (2) there remains an ambiguity whether to add a $`\pi ^0`$; or (3) there are several $`\pi ^0`$ and $`K_S^0`$ candidates to choose from. After selecting the best candidate, we require $`M_{X_s}<2.05\text{ GeV}/c^2`$.

The largest background source is from continuum $`q\overline{q}`$ production where the photon originates from initial state radiation or from high momentum neutral hadron decays such as $`\pi ^0\gamma \gamma `$ with one undetected photon. In order to reduce and estimate the contribution of the $`q\overline{q}`$ background, we introduce a new variable that exploits the topological difference between the jet-like $`q\overline{q}`$ events and the spherical $`B\overline{B}`$ events. We define the following variables in the $`B`$ meson rest frame,

$$R_l=\frac{\underset{i}{}|p_i||p_\gamma |P_l(\mathrm{cos}\theta _{i\gamma })}{_i|p_i||p_\gamma |},r_l=\frac{\underset{i,j}{}|p_i||p_j|P_l(\mathrm{cos}\theta _{ij})}{_{i,j}|p_i||p_j|},$$

where $`P_l`$ is a Legendre polynomial of order $`l`$, and $`i`$, $`j`$ run over all the photons and charged tracks that are not used to form the $`B`$ candidate. Then, we combine six of them into a Fisher discriminant

$$F=\underset{l=2,4}{}\alpha _lR_l+\underset{l=1,2,3,4}{}\beta _lr_l,$$

where the six coefficients $`\alpha _l`$ and $`\beta _l`$ are optimized to maximize the discrimination between signal and background as shown in Fig. 1. The terms $`R_1`$ and $`R_3`$ are excluded from the Fisher discriminant, because either of these terms is found to have a correlation with $`M_{bc}`$. We call $`F`$ the Super Fox-Wolfram (SFW) variable, since the terms are combined in such a way as to enhance the discriminating power of the original Fox-Wolfram moments . We select events with $`F>0.1`$, which is a compromise between the statistical significance of the signal and the size of the systematic error due to background subtraction. In order to estimate the $`q\overline{q}`$ background contribution, we use the SFW sideband region of $`F<1.5`$, where the signal fraction is only $`(0.7\pm 0.2)\%`$. This translates into a contribution of $`0.7\pm 0.2`$ events in the signal region.

In addition to the $`q\overline{q}`$ background, there is a contribution from $`bc`$ decays, in which $`BD^{()}\rho `$ modes are dominant. This background, which does not contribute to the SFW sideband, is estimated by MC to be $`9.1\pm 1.8`$ events and is subtracted from the data sample.

The signal yield is obtained from the $`M_{bc}`$ distribution (Fig. 2). In Fig. 3 we show that the ratio between the number of events in the signal region and in the sideband of the SFW variable is independent of the beam constrained mass $`M_{bc}`$ according to a MC calculation assuming continuum $`q\overline{q}`$ production. This observation is consistent with the background estimated using the off-resonance data. Thus, the SFW sideband data can be used to model the background shape of the $`M_{bc}`$ spectrum in the absence of a substantial sample of off-resonance data. We fit the $`M_{bc}`$ distribution with a background shape taken from the SFW sideband data and a signal shape obtained from MC simulation.

We observe $`222`$ events in the signal region $`M_{bc}>5.27\text{ GeV}/c^2`$, of which $`115.5\pm 7.7`$ is the estimated background contribution, and is consistent with the MC expectation. An excess of $`106.5\pm 16.8`$ events is clearly seen in the $`M_{bc}`$ distribution shown in Fig. 2 at the $`B`$ meson mass with a mass resolution of $`4.4\text{ MeV}/c^2`$ and a signal-to-background ratio of 0.9.

To obtain the background subtracted recoil mass spectrum shown in Fig. 4, we use the shape of the $`M_{X_s}`$ spectrum determined from MC, whose contribution is normalized to the estimated background of the $`M_{bc}`$ distribution. We adopt this method because the ratio between the number of signal events and SFW sideband events has a finite slope as a function of $`M_{X_s}`$, especially in the region $`M_{X_s}>2.05\text{ GeV}`$, and thus the SFW sideband does not properly represent the background sample for the $`M_{X_s}`$ spectrum . Nevertheless, the slope in the $`M_{X_s}<2.05\text{ GeV}`$ region is small and compatible with a flat distribution within the statistical fluctuation, and the background subtraction method using the SFW sideband sample leads to consistent results for the $`M_{X_s}`$ spectrum and the $`BX_s\gamma `$ branching fraction.

The signal MC sample is generated as a mixture of the exclusive $`BK^{}(892)\gamma `$ mode and an inclusive $`BX_s\gamma `$ contribution for $`M_{X_s}>1.15\text{ GeV}/c^2`$ in order to separate the $`K^{}(892)`$ signal from higher resonances starting with $`K_1(1270)`$. The recoil system $`X_s`$ is modelled as an equal mixture of $`\overline{s}d`$ and $`\overline{s}u`$ states . To simulate the inclusive $`M_{X_s}`$ spectrum, we adopt the model by Kagan and Neubert  with the $`b`$ quark pole mass parameter $`m_b=4.75\text{ GeV}/c^2`$.

The fraction of $`K^{}(892)\gamma `$ decays in the $`BX_s\gamma `$ transition, $`r_{K^{}\gamma }`$, is determined from the $`M_{X_s}`$ spectrum with a consideration of the migration effects on the reconstructed $`M_{X_s}`$ due to incorrect particle assignments to the $`X_s`$ final state. The $`M_{X_s}`$ spectrum is divided into two regions: the $`K^{}(892)`$ region below $`1.15\text{ GeV}/c^2`$ and the continuum region between $`1.15\text{ GeV}/c^2`$ and $`2.05\text{ GeV}/c^2`$. The MC sample is analyzed to determine the number of events migrating between the regions. The MC results are used to unfold the migration effects and to generate a set with an improved value, $`r_{K^{}\gamma }=(11.3\pm 3.5)\%`$. Using the new MC sample, the procedures of $`M_{bc}`$ fitting, $`M_{X_s}`$ background subtraction, and $`r_{K^{}\gamma }`$ unfolding are iterated. It is checked that $`r_{K^{}\gamma }`$ is stable under the iteration. The size of the migration effects and the value of $`r_{K^{}\gamma }`$ depend on the assumed $`m_b`$ value: for $`m_b=4.75\text{ GeV}/c^2`$, 77% of the events reconstructed below $`2.05\text{ GeV}/c^2`$ are genuine and 23% are from the region above $`2.05\text{ GeV}/c^2`$. We note that due to a non-negligible amount of migration between the regions below and above $`2.05\text{ GeV}/c^2`$, this boundary does not correspond to a sharp value in $`M_{X_s}`$ or $`E_\gamma `$.

According to the authors of ref. , the $`K^{}(892)`$ contribution can be approximated by integrating the $`M_{X_s}`$ spectrum up to a threshold mass above the $`K^{}(892)`$ mass, though their model does not explicitly take resonance contributions into account. This integral is sensitive to the parameter $`m_b`$. To match the measured $`K^{}(892)`$ fraction we varied $`m_b`$ and the threshold mass and observed that $`m_b4.9\text{ GeV}/c^2`$ is disfavored. We use $`m_b=(4.75\pm 0.10)\text{ GeV}/c^2`$ as the parameter region to test the model dependence of the branching fraction measurement . We assume a somewhat larger error on $`m_b`$ than that used for the semileptonic $`B`$ meson decays in the combined LEP, SLD and CDF evaluation  based on the recent theoretical work of , in which $`m_b(1\text{ GeV})=(4.58\pm 0.06)\text{ GeV}/c^2`$ is quoted. We note that the corresponding pole mass is $`m_b4.70\text{ GeV}/c^2`$. Other theoretical input parameters have little impact on the shape of the $`M_{X_s}`$ spectrum.

The event selection efficiency is determined from a MC sample that is calibrated with high statistics control data samples for all final state particles. The reconstruction efficiency for the primary photon is obtained using photons with energies between 2 and 3 GeV from the processes $`e^+e^{}e^+e^{}\gamma `$ and $`\eta \gamma \gamma `$. The charged tracking efficiency for momenta relevant to this analysis is obtained by comparing the high momentum $`\eta `$ yields for the $`\eta \pi ^+\pi ^{}\pi ^0`$ and $`\eta \gamma \gamma `$ decay modes. The kaon selection efficiency is determined from a sample of $`\varphi K^+K^{}`$ and $`D_s\varphi \pi `$ decays. The pion rejection efficiency is also determined from the $`\eta \pi ^+\pi ^{}\pi ^0`$ sample. The $`K_S^0`$ and $`\pi ^0`$ efficiencies are tested with a sample of $`K^{}(892)`$ decays into $`K^+\pi ^{}`$, $`K^+\pi ^0`$, $`K_S^0\pi ^+`$ and $`K_S^0\pi ^0`$. Combining all 16 different channels, the selection efficiency becomes $`(12.8\pm 1.3)\%`$. The $`M_{X_s}`$, SFW and $`M_{bc}`$ cuts reduce the efficiency by $`(57\pm 2)\%`$, $`(50\pm 3)\%`$ and $`(70\pm 1)\%`$, respectively. The SFW cut and the $`\pi ^0/\eta `$ veto efficiencies are measured by applying the cuts on a sample of exclusively reconstructed $`B^{}D^0\pi ^{}`$, $`D^0K^{}\pi ^+`$ decays, in which the first pion is assumed to be massless in order to mimic the primary photon for the cuts. The size of the systematic error from the efficiency determination is dominated by the statistical error of the control sample. The final event reconstruction efficiency is $`(2.58\pm 0.29)\%`$.

The total systematic error is quoted as a quadratic sum of the systematic errors on the efficiency, the best candidate selection procedure, signal MC statistics, SFW sideband uniformity and the number of $`B\overline{B}`$ events. The sensitivity of the cuts to the experimental error on the beam constrained mass $`M_{bc}`$ is tested by changing the MC resolution function, and the effect is found to be negligible. The efficiency errors on the 16 different final states are calculated individually and are then combined. The best candidate selection is tested with a $`B^{}D^0\pi ^{}`$ sample, which is analyzed using the inclusive reconstruction method of this analysis. The effect of the SFW sideband uniformity is tested by adding a slope to the background distribution within the tolerance level of Fig. 3. The contributions to the systematic error on the branching fraction are given in Table 1.

The fitted signal yield is extrapolated over the entire $`M_{X_s}`$ region using the model of ref.  with $`m_b=4.75\text{ GeV}/c^2`$ and corrected for the final states that are not reconstructed. The analysis is repeated for $`m_b=4.65\text{ GeV}/c^2`$ and $`4.85\text{ GeV}/c^2`$ to test the model uncertainty. The signal efficiency is found to be $`15\%`$ lower and $`16\%`$ higher, respectively. This error and the error of $`r_{K^{}\gamma }`$, which typically corresponds to $`5.0\%`$ in the branching fraction calculation, are combined into the theoretical model error.

Finally, we obtain the branching fraction of $`BX_s\gamma `$,

$$\mathrm{𝐵𝑟}(BX_s\gamma )=(3.36\pm 0.53\text{(stat)}\pm 0.42\text{(sys)}{}_{0.54}{}^{+0.50}\text{(th)})\times 10^4$$

where the first error is statistical, the second is systematic and the third is the theoretical model error. Using the values of the measured $`BX_s\gamma `$ branching fraction and the fraction of the $`K^{}(892)\gamma `$, we find the $`BK^{}(892)\gamma `$ branching fraction to be $`(3.8\pm 0.9)\times 10^5`$ (statistical error only), which is consistent with the $`BK^{}(892)\gamma `$ branching fraction measurements from CLEO  and a separate Belle analysis . In Fig. 5 we show the measured photon spectrum in which the background is subtracted and the efficiency loss due to the cut on $`M_{X_s}`$ is corrected by using MC input. This spectrum may be compared with model calculations.

To summarize, we have measured the inclusive branching fraction of $`BX_s\gamma `$ decay with the Belle detector. The result is consistent with the SM prediction of $`\mathrm{𝐵𝑟}(BX_s\gamma )=(3.28\pm 0.33)\times 10^4`$ . The result is also consistent with previous measurements by the CLEO  and ALEPH  experiments.

We gratefully acknowledge the efforts of the KEKB group in providing us with excellent luminosity and running conditions and the help with our computing and network systems provided by members of the KEK computing research center. We wish to thank M. Neubert and A. Ali for helpful suggestions on the modelling of the $`BX_s\gamma `$ MC. We thank the staffs of KEK and collaborating institutions for their contributions to this work, and acknowledge support from the Ministry of Education, Science, Sports and Culture of Japan and the Japan Society for the Promotion of Science; the Australian Research Council and the Australian Department of Industry, Science and Resources; the Department of Science and Technology of India; the BK21 program of the Ministry of Education of Korea and the Basic Science program of the Korea Science and Engineering Foundation; the Polish State Committee for Scientific Research under contract No.2P03B 17017; the Ministry of Science and Technology of Russian Federation; the National Science Council and the Ministry of Education of Taiwan; the Japan-Taiwan Cooperative Program of the Interchange Association; and the U.S. Department of Energy.


text:
# Explosive Decomposition in Ultrarelativistic Heavy Ion Collisions

## Abstract

Recent results from Au+Au collisions at BNL-RHIC energy hint at explosive hadron production at the QCD transition rather than soft hydrodynamic evolution. We speculate that this is due to a rapid variation of the effective potential for QCD close to $`T_c`$. Performing real-time lattice simulations of an effective theory we show that the fast evolution of the potential leads to “explosive” spinodal decomposition rather than bubble nucleation or critical slowing down.

It is hoped that heavy ion collision experiments at BNL-RHIC and later at CERN-LHC will provide evidence for the deconfined phase of QCD . Lattice gauge simulations predict a weakly first-order deconfinement transition at a critical temperature of $`T_c0.63\sqrt{\sigma }`$, where $`\sigma 1`$ GeV/fm is the string tension at $`T=0`$. For QCD with quarks, the transition could remain weakly first-order or become a smooth crossover. Even so, lattice data for two or three light flavors  indicate a steep rise in energy density and pressure with temperature and thus a rapid shift of the global minimum of the effective thermodynamic potential for QCD within a narrow temperature interval centered at $`T_c`$. To date, it has been difficult to identify unambiguous observables providing compelling evidence for a phase transition in ultrarelativistic heavy ion collisions (URHIC). However, rapid changes in the QCD effective potential near $`T_c`$ and/or supercooling followed by spinodal decomposition can have distinctive signatures. Their experimental observation would probably represent the most direct information which URHIC events can provide regarding the QCD phase transition. The purpose of this Letter is to consider the likelihood of explosive “spinodal” decomposition in RHIC experiments.

In a first-order phase transition, the thermodynamic potential for $`T`$ slightly less than $`T_c`$ exhibits a metastable minimum corresponding to the deconfined phase. At lower temperatures, this minimum becomes a point of inflection, i.e. a “spinodal instability”. In a slowly expanding system, such as the early universe, the phase transition proceeds through the nucleation of bubbles of the “true vacuum” state via thermal activation for all $`T<T_c`$ . In ultrarelativistic collisions of heavy ions or hadrons, expansion is very fast and the metastable deconfined state may reach the inflection point and undergo spinodal decomposition before significant nucleation has taken place. Calculations using an effective chiral model  suggest that spinodal decomposition is the favored mechanism in URHIC. The condensate field in each causally connected region “rolls down” to its minimum with large attendant fluctuations in the chiral order parameter and related hadronic observables.

Effects similar to those of spinodal decomposition can also be expected, independent of the order of the phase transition, for any rapidly varying thermodynamic potential. In this case, fields are unable to follow the rapid evolution of the equilibrium potential, and the system decays into any of a large number of unstable pion and sigma modes. (By contrast, the true spinodal decomposition of a slowly varying potential proceeds via the “critical” growth of long-wavelength modes of the order parameter . Here, we shall use the term “spinodal” somewhat loosely to describe any rapid decomposition.) This was first suggested in  in the context of chiral symmetry breaking. It was argued that, if expansion is sufficiently fast, the system can cool significantly while the chiral order parameter remains near its restored symmetry value. Violent decomposition would follow as a natural consequence of such a “quench”. Subsequent numerical simulations  showed that, if chiral symmetry breaking is considered independent of confinement, the chiral effective potential changes too slowly with temperature for such effects to be seen in URHIC. If, however, the effective potential changes rapidly with temperature , the quenched scenario and effective spinodal decomposition could be an interesting possibility even for moderate expansion rates.

RHIC ($`\sqrt{s}=130`$A GeV) has already produced a wealth of new results . “Dynamical” fluctuations in the average transverse momentum, $`p_t`$, are larger than at the SPS ($`\sqrt{s}=17`$A GeV). Observed Hanbury-Brown–Twiss (HBT) radii from pion interferometry are small, $`57`$ fm, and average transverse boost velocities are large, $`0.6`$c. In particular, the fact that the HBT outward and sideward radii, $`R_o`$ and $`R_s`$, are roughly equal for pion pair transverse momenta in the range $`100<p_t<500`$ MeV may indicate explosive decomposition. These phenomena all hint at a rapid out-of-equilibrium transition. By contrast, if a first-order transition proceeds at equilibrium, a mixed phase is created, and the system should spend a long time at  and slightly below  $`T_c`$. As a result, the outward HBT radius $`R_o`$ should be larger than the sideward radius $`R_s`$ by as much as $`50\%100\%`$ .

We wish to consider the possibility that the evolution of central RHIC events is determined by rapid changes in the thermodynamic potential. This would be the case, for example, if the chiral phase transition were driven by the confinement transition. The familiar rapid change of the energy density with temperature would then be reflected in equally rapid changes in the thermodynamic potential for the chiral fields. To illustrate this point, we will perform real time simulations of a weakly first-order deconfinement transition using a recently suggested effective theory, valid near $`T_c`$, in which the free energy (or pressure) is dominated by a condensate of Polyakov loops $`\mathrm{}`$ . The potential for $`\mathrm{}`$ is taken as

$$𝒱(\mathrm{})=\left(\frac{b_2}{2}|\mathrm{}|^2\frac{b_3}{6}(\mathrm{}^3+(\mathrm{}^{})^3)+\frac{1}{4}(|\mathrm{}|^2)^2\right)b_4T^4,$$
(1)
which is invariant under global $`Z(3)`$ transformations. The form of $`𝒱(\mathrm{})`$ is dictated by usual symmetry principles . The center symmetry of $`SU(3)`$, $`Z(3)`$, is broken at high $`T`$ where the global minimum of $`𝒱(\mathrm{})`$, $`\mathrm{}_0(T)1`$, and restored at low $`T`$ where $`\mathrm{}_0(T)0`$.

The assumption that quarks are not important for understanding the form of the effective potential for QCD is motivated by the results of lattice simulations. Specifically, plots of $`P/P_{\mathrm{ideal}}`$ versus $`T/T_c`$ are remarkably insensitive to the number of flavors, $`N_f=0`$, $`2`$, $`3`$ . Thus, we neglect terms linear in $`\mathrm{}`$, which must exist in QCD with quarks. Such terms would change the first-order transition (for $`N_f=0`$) into a crossover. Nevertheless, the fact that the free energy is small below $`T_c`$ means that their numerical effects must be small. We emphasize that the results (for rapid expansion) to be described below are not driven by the order of the phase transition and would not be changed in any significant way by the inclusion of such linear terms.

The coefficients $`b_2`$, $`b_3`$, and $`b_4`$ can be chosen to reproduce lattice data for the pressure and energy density of pure glue theory for $`TT_c`$ . In practice, to reproduce the pressure and energy density we may chose constant $`b_30.9`$ and $`b_415.0`$. For $`N_f=3`$, we rescale $`b_4`$ by a factor 47.5/16, corresponding to the increase in the number of degrees of freedom. In the spirit of mean field theory, the same constant values of $`b_3`$ and $`b_4`$ are to be employed for $`T<T_c`$. In this domain, $`b_2(T)0.66\sigma ^2(T)/T^4`$ can be set using the finite-$`T`$ string tension $`\sigma (T)`$ from the lattice .

The potential in eq. (1) changes extremely rapidly near $`T_c`$, see Fig. 1 of . For temperatures even $`\pm 2\%`$ away from $`T_c`$, the potential has only a single minimum, which corresponds to the deconfined (confined) state. This rapid change in the effective potential with temperature, which is directly related to the rapid increase of $`\sigma (T)`$ below $`T_c`$ seen in lattice data , makes “explosive” decomposition inevitable. (The rapid variation of $`b_2(T)`$ about $`T_c`$ can also be understood qualitatively by considering the above mean field theory for the second-order confinement transition, $`b_3=0`$, for the quenched theory with two colors. Then, $`b_2(T)=c(T/T_c1)`$. The rapid rise of the pressure relative to that for an ideal gas about $`T_c`$ , $`P(T)/P_{\mathrm{id}}(T)=b_2^2(T)`$, implies that $`c`$ is large.)

To complete the effective theory, we add a kinetic term for $`\mathrm{}`$ and allow coupling to a chiral field, $`\varphi `$, describing the usual pions and the sigma meson . Thus,

$$=_\varphi +\frac{N_c}{g^2}|_\mu \mathrm{}|^2T^2𝒱(\mathrm{})\frac{h^2}{2}\varphi ^2|\mathrm{}|^2T^2.$$
(2)
$`_\varphi `$ is the standard Lagrangian for the chiral fields, see e.g. . The coefficient of the kinetic term for $`\mathrm{}`$ is set here by the fluctuations of $`SU(3)`$ Wilson lines in space. Perturbative corrections to that coefficient are small . Thus, we assume a Lorentz invariant structure . The coupling $`h^222`$ between the chiral field and $`\mathrm{}`$ is chosen to reproduce $`m_\pi (T)`$ . This interaction term allows the Polyakov loop condensate to drive the chiral phase transition. While the inclusion of the chiral fields enables us to calculate familiar observables, it is otherwise largely passive. The behavior of $`\mathrm{}`$, i.e. confinement physics or its $`Z(3)`$ symmetry properties control the QCD transition. (In this regard, see also .)

In our numerical simulations we solve the classical Euler-Lagrange equations derived from the Lagrangian (2) on a $`64^3`$ space-like lattice. We employ a Robertson-Walker metric $`\mathrm{d}s^2=\mathrm{d}\tau ^2a^2(\tau )\mathrm{d}(\tau _0\eta )^2\mathrm{d}x_{}^2`$ with scale factor $`a(\tau )=\tau /\tau _0`$ and a Hubble constant $`H=\mathrm{d}\mathrm{log}(a)/\mathrm{d}\tau =1/\tau `$. Here, $`\tau `$ and $`\eta `$ are the longitudinal proper time and the space-time rapidity, respectively . We can vary the expansion rate by varying the initial time, $`\tau _0`$. Further, we impose periodic boundary conditions and choose a lattice spacing $`a_L=0.25`$ fm. We draw random initial field values at each lattice site from a Gaussian distribution and then coarse grain the lattice afterwards to a correlation length of $`2a_L`$. The initial temperature is chosen to be $`T(\tau _0)=1.01T_c`$. For the chiral field fluctuations, we assume  that $`\varphi =\dot{\varphi }=0`$, $`\varphi ^2=v^2/4`$, and $`\dot{\varphi }^2=v^2/\text{fm}^2`$ ($`v`$ is the vacuum expectation value of the sigma field, $`\sigma \varphi _0`$). The $`\mathrm{}`$-field is initialized in the deconfined minimum of the potential with real phase, $`\mathrm{}0.7`$. The magnitude of the initial fluctuations, $`\mathrm{}^2\mathrm{}^20.1`$ was chosen such that the initial energy density is $`4T_c^4`$ . The results are essentially independent of the initial distribution of energy between fluctuations of $`\mathrm{}(x)`$ and $`\dot{\mathrm{}}(x)`$.

The presence of fluctuations modifies the potential in eq. (1). We renormalize it by expanding $`𝒱(\mathrm{},\varphi )`$ at $`\tau =\tau _0`$ to second order in the fluctuations $`\delta \mathrm{}`$ and $`\delta \varphi `$ and then subtract the resulting $`\delta 𝒱(\mathrm{},\varphi ;\delta \mathrm{},\delta \varphi )`$. This provides us with the renormalized equations of motion.

We show the results of simulations for two different initial times; one appropriate for “slow” expansion ($`\tau _0=2000`$ fm/c) and another ($`\tau _0=100`$ fm/c) for “fast” expansion. These two cases will suffice to illustrate the radical changes which occur in the nature of the transition as a function of the expansion rate. The slow expansion is, of course, many orders of magnitude faster than that describing the QCD transition in the early universe, while the fast expansion is slow on typical URHIC scales, and thus the indicators of non-equilibrium phenomena shown below should be even stronger in real URHIC.

The potential energy of the $`\mathrm{}`$-field in the co-moving frame is assumed to red-shift according to $`T^4/T_0^4=1/a(\tau )\tau _0/\tau `$. However, the spectrum of long wavelength fluctuations on the lattice is determined by the dynamics of the phase transition; for example, in case of slow expansion with nucleation and bubble growth those fluctuations on the lattice may not cool during the phase conversion. Only hard, short wavelength fluctuations which are integrated out and determine the coefficients in the effective potential (1), but are not treated explicitly in the classical lattice simulation, cool according to the above prescribed $`T(\tau )`$ law.

Fig. 1 shows a sequence of histograms describing the distribution of $`|\mathrm{}|`$ for slow expansion. The transition is weakly first order and, accordingly, the effective potential for $`\mathrm{}`$ exhibits two local minima for $`T`$ in the immediate vicinity of $`T_c`$. Although $`\mathrm{}(x)`$ is initially localized in the deconfined minimum at $`\mathrm{}_00.7`$, the fluctuations serve as a heat bath and trigger thermal activation over the barrier. Fig. 1 shows that the confined and deconfined states coexist in different parts of the system. Well-defined bubbles form and, as time progresses, grow, coalesce, and eventually complete the transition to confinement.

Fig. 2 shows the corresponding histograms of $`|\mathrm{}|`$ for the fast expansion. Clearly, the double peak has disappeared, and the distributions merely broaden and move towards the confined minimum as the system cools. In this case, we conclude that the phase conversion happens through rapid decomposition. For a rough estimate of the relevant time scales, note that during a small time interval $`\mathrm{\Delta }\tau `$, the hard modes cool by $`|\mathrm{\Delta }T|/T_c\mathrm{\Delta }\tau /4\tau _0`$. From classical nucleation theory  a lower bound on the time scale for nucleation is $`\mathrm{\Gamma }^{1/4}\stackrel{>}{}1/T_c`$. For the model (1) with $`b_2(T)`$ set by the string tension $`\sigma (T)`$ , supercooling by $`1\%`$ suffices for the soft modes to “roll down” to the confined minimum  rather than nucleation to occur. Setting $`\mathrm{\Delta }\tau T_c1`$, we expect that for $`\tau _0\stackrel{<}{}100`$ fm/c phase coexistence is not established.

Fig. 3 describes the time evolution of $`|\mathrm{}(x)|`$, and its root mean square (RMS) fluctuations. In the slow case, we see that the mean field makes the transition from the deconfined minimum to the confined minimum very smoothly. The corresponding RMS fluctuations increase immediately before phase conversion and illustrate the competition between the two minima. (For a second-order phase transition and very slow expansion those fluctuations would go “critical” ). After the transition, the mean field shows no temporal oscillations and only small spatial fluctuations, indicating that the system is quite homogeneous after the transition.

By contrast, fast expansion leads to strong temporal oscillations of the mean field and large spatial fluctuations following the confinement transition. Thus, the scale of spatial homogeneity is smaller in this case. Also note that the rapid phase conversion and decoherence seen for fast expansion lead to the violent and rapid decay of the Polyakov loop condensate into physical particles. This, too, suggests explosive decomposition of the source. As discussed previously, strong event-by-event fluctuations of the pion mean-$`p_t`$ should emerge as a consequence . Also, the “sudden” transition could help to understand the flavor composition in the final state .

Fig. 4 summarizes the corresponding behavior of the chiral $`\sigma `$ field. The qualitative pattern here is the same, but the effects are enhanced. This indicates that substantial energy has been flushed quickly from the Polyakov loop condensate into the chiral field. This helps justify our expectation of the approximate equality of outward and sideward HBT radii: The outward radius $`R_o^2`$ is larger than the sideward radius $`R_s^2`$ by an amount roughly equal to the variance of the pion production time $`\delta \tau ^2=\tau ^2\tau ^2`$ . For a smooth quasi-adiabatic transition, entropy conservation in a comoving volume element requires $`\delta \tau ^2`$ to be large, while the non-equilibrium transition seen in our simulations is very fast.

To summarize, real-time simulations of the Polyakov loop condensate model have demonstrated that radically different phase transition dynamics can emerge in URHIC for realistically fast expansion rates when the effective potential changes rapidly in the vicinity of $`T_c`$. The model considered is consistent with lattice QCD data for the pressure, energy density, and correlation length (or string tension) near $`T_c`$. In real time, the picture that emerges is a confinement transition via spinodal decomposition which then triggers chiral symmetry breaking. The timescales involved in URHIC seem much too short to admit the traditional near equilibrium transition via bubble nucleation. As expected in any case of phase conversion via spinodal decomposition, we observe large RMS spatial fluctuations in the chiral fields after the transition, i.e. a small spatial homogeneity length. The near equivalence of outward and sideward HBT radii then follows as a natural consequence of the almost instantaneous decay of the Polyakov loop condensate into an “exploding” source of pions . Large event-by-event fluctuations of $`p_t`$ are another outcome of such a strong non-equilibrium transition.

###### Acknowledgements.

We thank J. Borg, D. Diakonov, H. Heiselberg, L. McLerran, I.N. Mishustin, R.D. Pisarski, D. Son, and R. Venugopalan for helpful criticism and discussions. A.D. gratefully acknowledges the support of the US-DOE through Contract No. DE-FG-02-93ER-40764.


text:
# Cost functions for pairwise data clustering

## Abstract

Cost functions for non-hierarchical pairwise clustering are introduced, in the probabilistic autoencoder framework, by the request of maximal average similarity between the input and the output of the autoencoder. The partition provided by these cost functions identifies clusters with dense connected regions in data space; differences and similarities with respect to a well known cost function for pairwise clustering are outlined.

Clustering methods aim at partitioning a set of data-points in classes such that points that belong to the same class are more alike than points that belong to different classes . These classes are called clusters and their number may be preassigned or can be a parameter to be determined by the algorithm. There exist applications of clustering in such diverse fields as pattern recognition , astrophysics , communications , biology , business  and many others. Two main approaches to clustering can be identified: parametric and non-parametric clustering.

Non-parametric approaches make few assumptions about about the data structure and, typically, follow some local criterion for the construction of clusters. Typical examples of the non-parametric approach are the agglomerative and divisive algorithms that produce dendrograms. In the last years non-parametric clustering algorithms have been introduced employing the statistical properties of physical systems. The Super-Paramagnetic approach by Domany and coworkers  exploits the analogy to a model granular magnet: the spin-spin correlation of a Potts model, living on the data-points lattice and with pair-couplings decreasing with the distance, is used to partition points in clusters. The synchronization properties of a system of coupled chaotic maps are used in  to produce hierarchical clustering.

Parametric methods make some assumptions about the underlying data structure. Generative mixture models  treat clustering as a problem of density estimation: data are viewed as coming from a mixture of probability distributions, each representing a different cluster, and the parameters of these distributions are adjusted to achieve a good match with the distribution of the input data. This can be obtained by maximizing the data likelihood (ML) or the posterior (MAP) if additional prior information on the parameters is available .

Many parametric clustering methods are based on a cost function: the best partition of points in clusters is assumed to be the one with minimum cost. Often cost functions incorporate the loss of information incurred by the clustering procedure when trying to reconstruct the original data from the compressed cluster representation: the most popular algorithm to optimize a cost function is $`K`$-means . Starting from a statistical ansatz and invoking maximum likelihood leads to a cost function which has been observed to work for clustering financial time series .

It is important to stress the difference between central clustering, where it is assumed that each cluster can be represented by a prototype , and pairwise clustering where data are indirectly characterized by pairwise comparison instead of explicit coordinates ; pairwise algorithms require as input only the matrix of dissimilarities. Obviously the choice of the measure of dissimilarity is not unique and it is crucial for the performance of any pairwise clustering method. It is worth remarking that it often happens that the dissimilarity matrix violates the requirements of a distance measure, i.e. the triangular inequality does not necessarily holds.

Folded Markov chains are used in the Probabilistic Autoencoder Framework to derive cost functions for clustering . Some examples of two-stage folded Markov chains, and the corresponding algorithms for clustering and topographic mapping , are thoroughly analyzed in , where it is also shown that the cost function for pairwise clustering, introduced in , may be seen as a consequence of Bayes’ theorem and the requirement of minimal average distorsion in a probabilistic autoencoder.

It is the purpose of this work to introduce a new class of cost functions for pairwise clustering which can be obtained, in the autoencoder frame, by requiring maximal similarity instead of minimal distorsion. We show that the cost functions here introduced provide a non-hierarchical clustering of points where dense connected regions of points in the data space are recognized as clusters.

Let us now discuss autoencoders described by one-stage folded Markov chains. Let us consider a point $`x`$, in a data space, sampled with probability distribution $`P_0\left(x\right)`$; a code index $`\alpha \{1,\mathrm{},q\}`$ is assigned to $`x`$ according to conditional probabilities $`P\left(\alpha |x\right)`$. A reconstructed version of the input, $`x^{}`$, is then obtained by use of the Bayesian decoder:

$$P\left(x^{}|\alpha \right)=\frac{P\left(\alpha |x^{}\right)P_0\left(x^{}\right)}{P\left(\alpha \right)}.$$
(1)
The joint distribution of $`x`$, $`x^{}`$ and $`\alpha `$, describing this encoding-decoding process, is

$$P(x,x^{},\alpha )=P_0\left(x\right)P\left(\alpha |x\right)P\left(x^{}|\alpha \right);$$
(2)
owing to (1), the joint distribution reads:

$$P(x,x^{},\alpha )=\frac{P_0\left(x\right)P_0\left(x^{}\right)P\left(\alpha |x\right)P\left(\alpha |x^{}\right)}{P\left(\alpha \right)}.$$
(3)
The conditional probabilities $`\{P\left(\alpha |x\right)\}`$ are the free parameters that must be adjusted to force the autoencoder to emulate the identity map on the data space.

Let $`d(x,x^{})`$ be a measure of the distorsion between input and output of the autoencoder. The average distorsion is then given by:

$$𝒟=\underset{\alpha =1}{\overset{q}{}}𝑑x𝑑x^{}\frac{P_0\left(x\right)P_0\left(x^{}\right)P\left(\alpha |x\right)P\left(\alpha |x^{}\right)}{P\left(\alpha \right)}d(x,x^{}).$$
(4)
Moreover, let s(x,x’) be a measure of the similarity between input and output; the average similarity is then given by

$$𝒮=\underset{\alpha =1}{\overset{q}{}}𝑑x𝑑x^{}\frac{P_0\left(x\right)P_0\left(x^{}\right)P\left(\alpha |x\right)P\left(\alpha |x^{}\right)}{P\left(\alpha \right)}s(x,x^{}).$$
(5)
It is natural to postulate a one-to-one mapping between values of distorsion and similarity, $`s=F(d)`$, with $`F`$ a strictly decreasing function. A good autoencoder is obviously characterized by a low value of $`𝒟`$ and high value of $`𝒮`$. However we remark that the two requirements $`Min(𝒟)`$ and $`Max(𝒮)`$, for reasonable choices of $`F`$, are not generically equivalent.

Now we turn back to the clustering problem. Given a data-set $`\{x_i\}`$ of cardinality $`N`$, partitioning these points in $`q`$ classes corresponds, in this frame, to design an autoencoder, with $`q`$ code indexes, acting on data space. We choose the encoder to be deterministic:

$$P\left(\alpha |x\right)=\delta _{\alpha \sigma (x)},$$
(6)
$`\sigma (x)\{1,\mathrm{},q\}`$ being the code index associated to $`x`$. The estimate for the average distorsion (4), based on the data-set at hand, is given by $`\widehat{𝒟}=NH_d[\sigma ]`$, where we introduce the hamiltonian $`H_d`$ for the Potts variables $`\{\sigma _i\}`$:

$$H_d[\sigma ]=\underset{\alpha =1}{\overset{q}{}}\frac{\underset{i,j=1}{\overset{N}{}}\delta _{\alpha \sigma _i}\delta _{\alpha \sigma _j}d_{ij}}{_{k=1}^N\delta _{\alpha \sigma _k}},$$
(7)
where $`\sigma _i=\sigma (x_i)`$, $`d_{ij}=d(x_i,x_j)`$. It turns out that $`H_d`$ is equivalent to the cost function for pairwise clustering, influential in the clustering literature, introduced in .

The estimate for the average similarity is, similarly, given by $`\widehat{𝒮}=NH_s[\sigma ]`$, where we introduce the hamiltonian $`H_s`$:

$$H_s[\sigma ]=\underset{\alpha =1}{\overset{q}{}}\frac{\underset{i,j=1}{\overset{N}{}}\delta _{\alpha \sigma _i}\delta _{\alpha \sigma _j}s_{ij}}{_{k=1}^N\delta _{\alpha \sigma _k}}.$$
(8)
If we choose the autoencoder by minimizing the average distorsion, then the best partition of the data-set in $`q`$ classes corresponds to the ground state of $`H_d`$. If we choose it by maximizing the average similarity, then the ground state of $`H_s`$ must be sought for, instead. Since both $`\{d_{ij}\}`$ and $`\{s_{ij}\}`$ may be taken positive, it follows that $`H_d`$ is characterized by antiferromagnetic couplings between the Potts variables, while $`H_s`$ is made of ferromagnetic couplings. Denominators in both $`H_d`$ and $`H_s`$ serve to enforce the coherence among the $`q`$ clusters. In particular, without the denominator the ground state of $`H_s`$ would correspond to a single big cluster.

The form of the function $`F`$, determining the relation between $`s`$ and $`d`$, has to be specified. In what follows we consider two forms of this relation. A scale-free relation

$$s_{ij}=F_\gamma (d_{ij})=\left(\frac{d_{ij}}{d}\right)^\gamma ,$$
(9)
depending on the exponent $`\gamma `$, and a scale-dependent relation

$$s_{ij}=F_a(d_{ij})=\mathrm{exp}\left(\frac{1}{2a^2}\left(\frac{d_{ij}}{d}\right)^2\right),$$
(10)
dependent on the scale $`a`$. In the formulas above, $`d`$ is the average dissimilarity over all the pairs of data-set points. The exponent $`\gamma `$ will be restricted to assume small values so as to characterize the corresponding Potts model by long-range ferromagnetic couplings; the scale parameter $`a`$ will be bounded in $`[0,1]`$.

At this point it is worth stressing that minimization of the distorsion and maximization of the similarity yield, in the autoencoder frame, different cost functions. The hamiltonian $`H_d`$ embodies the requirement that pairs of distant points (large $`d_{ij}`$) should belong to different clusters. On the other hand, the hamiltonian $`H_s`$, for reasonable choices of $`F`$, concentrates on pairs of close points (small $`d`$) and forces them to belong to the same cluster. In other words, $`H_s`$ may be seen to implement the idea that clusters should be searched for as dense connected regions in the data space.

We describe now the application of the variational criterions for clustering, described above, to some artificial and real data-sets. We consider two optimization algorithms to find the configuration of minimum cost: simulated annealing  and mean-field annealing . Both approaches associate a Gibbs probability distribution to the functional to be optimized. Simulated annealing is a Monte-carlo technique which samples the Gibbs distribution as the temperature is reduced to zero, while mean-field annealing attempts to track an approximation, to the mean of the distribution, known as mean field approximation . We remark that an efficient mean-field annealing algorithm for cost function (7), based on the EM scheme , is described in : the generalization of that algorithm to (8) is straightforward.

In many cases cost functions $`H_d`$ and $`H_s`$ have very close global minima. For example in Fig.1a we depict an artificial data-set generated by two overlapping isotropic Gaussian distributions. In this case the natural measure of dissimilarity is Euclidean metrics, and we use $`q=2`$. In Fig.1b the corresponding ground state of $`H_d`$  is depicted: it is very close to the Bayesian solution, i.e. the solution obtained drawing the symmetry plane for the centers of the two Gaussians. A similar partition is obtained minimizing, by simulated annealing, $`H_s`$. As a measure of the difference between two partitions $`\{\sigma _i\}`$ and $`\{\eta _i\}`$, we evaluate the following quantity:

$$ϵ=\frac{1}{N(N1)}\underset{i=1}{\overset{N}{}}\underset{j=1,ji}{\overset{N}{}}\left(\delta _{\sigma _i\sigma _j}\delta _{\eta _i\eta _j}\right)^2$$
(11)
which counts the number of pairs of points upon which the two partitions disagree. Using the scale-dependent $`F_a`$, we find the ground state of $`H_s`$ to differ from those of $`H_d`$ by $`ϵ<0.01`$ varying $`a`$ in $`[0.05,1]`$. Analogously, using the scale-free $`F_\gamma `$, with $`\gamma [0.1,1.5]`$, we find $`ϵ<0.02`$ when we compare the ground state of $`H_s`$ with those of $`H_d`$. Hence, on this data set, the cost functions introduced above work similarly within wide ranges of $`\gamma `$ and $`a`$ values.

We find a similar behaviour with respect to the famous IRIS data of Anderson . This data set has often been used as a standard for testing clustering algorithms: it consists of three clusters (Virginica, Versicolor and Setosa) and there are $`50`$ objects in $`𝐑^4`$ per cluster. Two clusters (Verginica, Versicolor) are very overlapping. The clustering result, with $`q=3`$ and minimizing $`H_d`$, consists of three clusters of $`61`$, $`39`$ and $`50`$ points respectively, with $`90\%`$ of correct classification percentage. We obtain exactly the same partition by minimizing $`H_s`$ using a scale-free $`F`$ (with $`\gamma [0.15,1.45]`$), and using a scale-dependent $`F`$ (with $`a[0.25,1]`$). For $`a[0.1,0.25]`$ we obtain, in the scale-dependent case, a slightly different partition with clusters’ sizes $`58`$, $`42`$, $`50`$ and correct classification percentage $`93.3\%`$. These results show that also in the IRIS case the pairwise clustering procedures by distorsion minimization and similarity maximization are almost equivalent.

A typical situation resulting in different answers from $`H_d`$ and $`H_s`$ is depicted in Fig.2a. This two-dimensional data-set is made of an elongated cluster and a Gaussian distributed circular one. It is evident that two dense connected regions are present, and that the farthest pairs of points belong to the same connected region. This is the type of data-set such that minimizing the distorsion is not equivalent to maximizing the similarity. In fig.2b the partition we obtain minimizing $`H_d`$ is depicted: it fails to recognize the structure in the data-set. Let us now consider the ground state of $`H_s`$ with the scale-dependent $`F`$. For $`a<0.7`$ the ground state, depicted in Fig.2c, recognizes with $`99\%`$ accuracy the data structure. At $`a0.7`$ a transition phenomenon occurs: the configuration depicted in Fig.2c ceases to be the global minimum, the new ground state (Fig.2d) being very close to the solution by $`H_d`$.

In Fig.3a we depict the efficiency of the classification versus the resolution parameter $`a`$, for the scale dependent $`F`$, while in Fig.3b we consider a sequence of $`a`$-values and we plot the $`ϵ`$ between partitions corresponding to adjacent values of $`a`$. The peak at $`a=0.7`$ is the indicator of the transition between global minima. Finally, in Fig. 3c the size of the two clusters, versus $`a`$, is depicted. Concerning the scale-free $`F`$, in Fig.4 the same plots as in Fig.3 are depicted, showing that the good minimum is stable for a wide range of $`\gamma `$.

The choice of the optimization algorithm deserves a comment. All the results described above are obtained by simulated annealing; we also apply the mean-field annealing scheme, described in , and we always find a configuration very close to the one from simulated annealing, while spending less computational time. This confirms that optimization algorithms rooted on mean-field theory yield quickly a good solution on these problems .

In summary, we address non-hierarchical pairwise clustering and, working in the probabilistic autoencoder frame, we introduce a class of cost functions arising from the request of maximal average similarity between the input and the output of the autoencoder. Our simulations show that the partition provided by these new cost functions corresponds to extract dense connected regions in data space, and that a relevant discrepancy with the partition provided by the cost function introduced in  is to be expected in case of non-trivial geometry of clusters. We note that the approach to clustering here described has some similarities with the method in : indeed in both cases clustering is mapped onto a ferromagnetic Potts model with couplings decreasing with the distance. In the superparamagnetic approach, however, $`q`$ is not related to the number of classes present in the data-set and one obtains hierarchical clustering as the temperature of the Potts model is varied. In the present case $`q`$ is the number of classes, which is supposed to be known (non-hierarchical clustering), and the denominators in the hamiltonian, ensuring clusters’s coherence, leads to a non-trivial ground state which reflects data structure. We consider two classes of cost function. Scale-free cost functions depend on the exponent $`\gamma `$, while scale-dependent ones depend on the scale-parameter $`a`$. Varying $`a`$, i.e. changing the resolution at which the data-set is processed, may give rise to transitions between different partitions; in the scale-free case, the clustering output is fairly stable, with respect to $`\gamma `$, in a wide range.

Further work will be devoted to test these new cost functions on other real applications and to study related issues, such as the introduction of an adaptive relation between distorsion and similarity, i.e. the function $`s=F(d)`$ might be depending on the properties of the data-set in a neighbourhood of the pair of points under consideration. It will be also important to develop cluster-validity criterions to provide a means to choose an optimal $`q`$ value in situations where the number of classes is ambiguous.