text:
# Spiral Galaxies with HST/NICMOS. I. Nuclear Morphologies, Color Maps and Distinct Nuclei1footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by Association of Universities for Research in Astronomy, Inc. (AURA), under NASA contract NAS5-26555

## 1 Introduction

The nuclear regions of spiral galaxies may hold important information concerning the growth of bulge-like structures within the disks, and more generally the secular evolution of the disks themselves. Little was known about the $`10`$-100 pc scale regions of spiral galaxies before we started our HST program, a cycle 6 WFPC2 F606W ($`V`$) survey of 75 mostly early- and intermediate-type disk galaxies selected randomly from a complete sample of 134 galaxies (Carollo et al. 1997; Carollo, Stiavelli & Mack 1998; Carollo & Stiavelli 1998; Carollo 1999). The random selection is due to the snapshot nature of the program. The unprecedented angular resolution of the HST optical images allowed probing the bulges down to absolute $`V`$ magnitudes $`M_V14`$, with half-light radii $`R_e`$ as small as $`100`$pc, and unveiled an astonishing richness of structure in the nuclear regions of $`60\%`$ of the sample: bars, spiral-like dust lanes, star-forming rings or spiral arms. Even more surprisingly, a similar fraction of objects contains photometrically-distinct nuclear compact sources (hereafter ‘nuclei’). While a few of these nuclei appear point-like in the WFPC2 images, many are resolved with $`R_e0.1^{\prime \prime }0.2^{\prime \prime }`$, corresponding to linear scales of up to $`20`$ pc. The nuclei were identified primarily by visual inspection of the images, since they appear as regular, smooth structures significantly brighter than their surroundings. For the subsample of objects with a light distribution smooth enough to allow the derivation of an isophotal fit, the morphological identification of the nuclei was confirmed by inspecting the surface brightness profiles. The nuclei embedded in those hosts whose dusty/star-forming morphology did not allow us to obtain a meaningful isophotal fit are typically much brighter than the nuclei embedded in the smoother hosts, so that their identification as photometrically-distinct systems remains unambiguous. Carollo et al. (1997,1998) presented a collection of $`V`$-band images and light profiles which clearly showed the photometric decoupling between the nuclei and their underlying host galaxies. Faint nuclei embedded in bright hosts with steep nuclear stellar cusps would remain undetected (this is discussed quantitatively in the references above).

The relation between the nuclei and their galactic surroundings is unclear. In principle, the nuclei could be triggers for —and/or left-overs of— the mechanisms which produce central structure within the disks, and could thus provide additional constraints for identifying such mechanisms. Interestingly, the WFPC2 survey has shown that many of the nuclei at the faint-end of their luminosity distribution ($`9>M_V>12`$) are often hosted by pseudo-bulges with an exponential —rather than $`R^{1/4}`$-law— fall-off of the light distribution (hereafter exponential bulges; Courteau, de Jong & Broeils 1996 and references therein; see also the review by Wyse, Gilmore & Franx 1997). In fact, every exponential bulge in our WFPC2 sample hosts a distinct nucleus in its center or slightly offset from it. Distinct central structures significantly different from the massive, $`R^{1/4}`$-law elliptical-like spheroids are rather common in the intermediate-type Sb-Sc disks. In the latter, in addition to the exponential bulges, pseudo-bulges have been found with cold kinematics (Kormendy 1993) or a peanut-shape morphology (Kuijken & Merrifield 1995; Bureau & Freeman 1999). The presence of the nuclei in the centers of the exponential bulges does not necessarily imply a causal connection between nuclei and bulges. At least at the heuristic level, however, such a connection may be indirectly supported on a theoretical basis, since evolutionary scenarios for bulge formation advocate for example the dissolution of a progenitor bar by means of the accumulation of a $`1\%`$ central mass (Pfenniger & Norman 1990; Hasan, Pfenniger & Norman 1993; Pfenniger 1993; Norman, Sellwood & Hasan 1996, and references therein).

A first step toward understanding the nucleus-host relationship is to establish the properties of the nuclei as a function of the Hubble type and other physical properties of their parent galaxies. The findings based on the optical WFPC2 data could potentially be severely affected by large amounts of dust and of recent star formation that may be present in the central regions of spirals. To secure results less influenced by these sources of ‘pollution’, we have followed up the WFPC2 survey with a NICMOS Camera-2 F160W ($`H`$; in some cases, also F110W, i.e., “broad” $`J`$) snap-shot survey for the same sample of galaxies.

The analysis of our NICMOS survey is reported in this and in a companion paper (Seigar et al. 2001, paper II). In this paper we: (i) Introduce the NICMOS sample and describe our technique of data reduction; (ii) Describe the NIR images and the optical-NIR color maps of spirals at HST resolution; (iii) Detail the methodology used to measure the sizes and luminosities of the nuclei from the NICMOS images, and describe the tests performed to assess the systematic errors on these measurements; (iv) Discuss the statistical NIR and optical-NIR properties of the nuclear regions of spirals, including the distinct nuclei, as a function of Hubble type. In paper II we present the isophotal fits to the galaxy nuclear surface brightness distribution, and the analytical fits to these profiles. Paper II also reports the derived nuclear stellar cusp slopes and the optical/NIR color profiles (for those galaxies for which both WFPC2 and NICMOS images are available). The data and measurements presented here and in paper II are used in an additional paper (Carollo et al. 2001) to discuss in more detail the nucleus-bulge connection. There we present the optical/NIR color-distribution, color-magnitude and color-color diagrams for the $`R^{1/4}`$ and exponential bulges and for their nuclei, and discuss the possible implications for the formation and evolution of the central regions of early/intermediate-type disk galaxies.

In §2 we describe the NICMOS sample, the observations, the data reduction procedure and the approach adopted to derive the radii and luminosities of the nuclei. In the same section we present the tests that we performed to assess the accuracy of the measurements for the nuclei, and show the optical-NIR or NIR-NIR color maps for the target galaxies. We discuss in §3 the sizes and magnitudes of the nuclei in the NIR, their colors, as well as the structure and color of their surrounding environments. We sumarize in Section 4. Consistent with our earlier work, we adopt $`H_0=65`$ km/s/Mpc throughout this paper.

## 2 Sample and Data Analysis

The NICMOS sample presented here constitutes a random selection from the sample of 134 galaxies discussed in Carollo et al. (1997, 1998). The objects were extracted from the UGC catalog for the northern hemisphere (Nilson 1973) and the ESOLV catalog for the southern hemisphere (Lauberts & Valentijn 1989) for having angular diameter larger than 1 arcmin, redshift $`<`$ 2500 km s-1, inclination angle smaller than 75 degrees, and (non-barred) morphological type from Sa to Sbc. A few objects that did not satisfy these selection criteria (e.g., NGC 1365 and NGC 4565) were however included in the final sample in substitution of objects which we were not allowed to observe due to the HST rules on duplication of targets. A posteriori, the WFPC2 images showed that the RC3 catalog (de Vaucouleurs et al. 1991) gives generally a better description of the galaxy types. Therefore, consistent with our previous WFPC2 work, we adopt in this analysis the RC3 Hubble types (with the numerical conversion given in Table 2 of RC3: 0=S0a, 1=Sa, 2=Sab, $`\mathrm{}`$, 9=Sm).

### 2.1 The Dataset and the Basic Data Reduction

All NICMOS observations were performed with fine-lock guiding. 78 of the 134 systems were observed in snap-shot mode with NICMOS Camera-2 in F160W during the period June 1997 to February 1998.

Nine objects, NGC1365, NGC2903, NGC3031, NGC4536, NGC4565, NGC4569, NGC6217, NGC7479 and NGC7742 were observed with WFPC2 in the context of other programs; consequently, we were not allowed to ‘duplicate’ the optical observations for these 9 systems. The archival images for these galaxies, however, were not taken in the F606W filter which we had chosen for our own WFPC2 survey. Given the arbitrary assumptions that would underlie the conversion from the “archival”-F160W into the F606W-F160W color (a conversion necessary for studying these systems together with the other objects of our sample), we have not included these 9 galaxies in the statistical study presented in this paper, to guarantee the homogeneity of the results. The $`H`$ images of these systems at HST resolution are presented in Appendix A. The remaining 69 objects are considered in this study. At the time of the NICMOS Phase II submission, no WFPC2 data were available either in the archive or as part of our snapshot program for 46 of the targets. For these objects a F110W image was also included in our NICMOS observing plan; 41 of these 46 systems were actually observed in our NICMOS survey with both the F160W and the F110W filter. These F110W and F160W exposures consist of two MULTIACCUM sequences with STEP256 and NSAMP=10, for a total exposure time of 256 seconds per filter. For 30 of these 46 galaxies WFPC2 F606W images were obtained, as part of our WFPC2 snapshot program, after the NICMOS phase II deadline. Thus, images in the three filters F606W, F110W and F160W, are available for these 30 systems. For the remaining targets F606W WFPC2 images were already in hand at the time of the NICMOS phase II deadline. For these objects we acquired three F160W MULTIACCUM sequences with STEP256 and NSAMP=10, which results in a total exposure time of 384 seconds. Images in both F160W and F606W are available for a total of 58 galaxies. Table 1 lists the 69 galaxies considered in this study together with the list of filters in which they have been observed, the galaxy distance in Mpc, morphological type, apparent total $`B`$ magnitude, $`B`$ extinction, as well as the information on whether an $`R^{1/4}`$ or an exponential bulge is present, and the type of central spectrum observed at ground-based resolution.

for these objects, including whether they were also observed within our program in F110W and/or in F606W.

The NICMOS images were reduced with the standard pipeline software, using the on-orbit flats and darks identified in the HST archive as the best-reference frames. The various NICMOS image-anomalies were corrected on a case-by-case basis. All images were affected by the so-called ‘pedestal’ anomaly, a time-variable bias of unclear origin. In most of the images the pedestal was removed using the NICMOS Pedestal Estimation Software (developed by R.P. van der Marel). In about 10$`\%`$ of the images the pedestal effects were very severe, and they had to be removed iteratively. In some galaxies (ESO240G12, ESO404G3, ESO443G80, ESO548G29, ESO549G18, ESO572G22, IC1555, NGC1345, NGC1800, NGC3455, NGC4806, NGC4980, NGC7188 and NGC7241), even after removal of the pedestal, a constant gradient across the chip remained visible. This was removed by subtracting a linear fit to the edges of the chip, using the ‘imsurfit’ routine in IRAF. In about 30$`\%`$ of the images additional anomalies such as shading and bias jumps/bands were present. These were also removed from the images before performing the measurements.

The combined images were flux-calibrated adopting the pipeline flux calibration zeropoints, ZPF110W = 23.248 for the F110W filter and ZPF160W = 23.110 for the F160W filter (AB system, Oke & Gunn 1983). Magnitudes for our WFPC2 data were originally given in “Vegamag”; we converted them to the AB system by using V606,AB = V606,Vega+0.117. The conversion between AB and Vegamags was computed by us using the Synphot synthetic photometry package (in IRAF) with the spectrum of Vega as input.

$`VH`$ color maps were derived for the NICMOS galaxies for which WFPC2 F606W images were also available, including 8 objects for which the optical images were obtained after the publication of the WFPC2 results (asterisks close to the galaxy name in Table 1). To account for the difference in Point Spread Function (PSF) between the $`H`$ and $`V`$ images, the images in a given band were convolved with the PSF of the other band. All PSFs were obtained with TINYTIM (Krist 1997). The adopted extent of the PSF images generated by TinyTim was 3′′. The PSF FWHM was 0.127′′ for the F160W filter, and 0.091′′ for the F110W filter, respectively; a $`5\times 5`$ sub-pixel scale was used. The $`V`$ and $`H`$ images were aligned after sub-pixel rebinning before computing the color maps. For the 11 galaxies for which the only additional images available were those in the F110W filter, we computed the $`JH`$ color maps following the same approach just described for the optical-NIR maps.

The identification of the nuclei in the NIR images was based primarily on visual inspection of the images, and was checked on the radial light profiles for those systems for which an isophotal fit could be derived (paper II). Similar to what was done for the visual WFPC2 data, the nuclei were selected based on the criterion that their compact, smooth, regular-morphology excess of light was well-isolated from the remaining circum-nuclear structure. In Figure 1 a we show the $`H`$ images and $`VH`$ color maps for 38 of the 39 galaxies that host a distinct nucleus and that were imaged both in F160W and F606W (the image of the 39$`th`$ object, ESO549G18, had a S/N sufficient for measuring the parameters of the nucleus, but not for studying the underlying galaxy); in Figure 1 b we show instead the $`H`$ and $`JH`$ images for the 5 galaxies which host a distinct nucleus for which the only other filter available in addition to F160W is F110W. Figures 2 a and b show the similar images for those galaxies which either do not host a distinct nucleus, or in which the presence of a nucleus is uncertain. In Figures 3 a and b we show the $`H`$-band light profiles of two galaxies representative of the two subsets presented in Figures 1 and 2, respectively: the first galaxy hosts a distinct nucleus embedded in an exponential pseudo-bulge (panel a; NGC 2082); the other galaxy hosts an $`R^{1/4}`$-law bulge (panel b; NGC 2460). In the NGC 2082 the photometrically-distinct nucleus is clearly identifiable in the $`H`$ light profile inside the innermost few tenths of arcsecond, similarly to what was found in the visual band. In contrast, in the NIR as in the $`V`$ band, the $`R^{1/4}`$-law bulge shows the steep luminosity profile typical of this class of objects (Carollo et al. 1998; Carollo & Stiavelli 1998; paper II).

### 2.2 Radii and Luminosities of the Nuclei in the NIR

The NIR measurements of the magnitudes and radii of the nuclei are affected by the uncertainty on the light contribution from the underlying galaxy. To quantify the luminosities of the nuclei in the NIR we adopted the three approaches previously optimized to extract similar information from the WFPC2 visual images, namely: (1) we fitted the contribution from the nucleus with a Gaussian, assuming that the underlying galactic light continuum is well represented by the asymptotic value of the Gaussian wings; (2) we modeled the underlying galactic light with a fourth-order polynomial using the task “fit/flat” in ESO/MIDAS, subtracted this model from the image, and attributed the excess counts to the nucleus; (3) we subtracted from the original frame an iteratively smoothed and nucleus-masked version of it, and performed the measurement on the resulting image containing only the nucleus. These three methods are complementary: they sample the galactic background in different regions, and hence their combined use allows us to estimate the uncertainty in the derived values contributed by the poorly constrained underlying galactic light. As done for the optical measurements, we adopted as final the averages of, and as errorbars the mean of the differences between, these three estimates for the luminosity of the central sources. The sizes of the compact sources were taken equal to the FWHM of the best fitting Gaussians of method (1), corrected for the instrumental width of NICMOS assuming that the intrinsic profile of the nucleus is described by a Plummer law (see Eq. 2). Consistent with the WFPC2 study, in the following we adopt for the half-light radii of the nuclei the ratio between half-light radius and FWHM valid for a Plummer law (see below). This gives a better representation than, e.g., a Gaussian shape, given that all compact sources show wings in excess of a Gaussian. An assessment of the systematic uncertainties associated with the adopted radii and luminosities of the nuclei is given in Appendix B.

In our WFPC2 study of the nuclei (Carollo et al. 1997; Carollo et al. 1998), for those galaxies with an isophotal/analytical fit available, we also estimated the magnitudes and luminosities of the nuclei by using two additional approaches: (a) We adopted the best fitting (exponential) models of the bulges to describe the underlying galaxian light, subtracted these models from the images, and attributed the residual nuclear counts to the nuclei; (b) We fitted the entire radial light profiles, including the innermost points, with multi-component analytical laws representing the bulge with an exponential law and the nucleus with either an $`R^{1/4}`$-law, or a Plummer law, or a modified Hubble profile. These alternative approaches provided results in agreement with —and well within the range spanned by— those obtained using the methods (1), (2), and (3) discussed above. To guarantee a homogeneous and self-consistent set of measurements, however, in our WFPC2 analysis we only used the results obtained from the latter three methods, since the additional two approaches could only be applied to a subset of the sample. Due to the limited radial extension of the NICMOS data, no bulge/disk decomposition was performed on the NIR data, and therefore no exponential fits are available for the bulges in $`H`$ or $`J`$. Thus, these additional tests could not be performed on the NIR profiles. Nonetheless, the results of the tests conducted on the WFPC2 data give us confidence that also for the NICMOS measurements the error bars derived from the “primary” methods (1), (2) and (3) do reliably represent the uncertainty associated with disentangling the unknown contribution from the underlying galaxy.

Finally, we computed the location of the nuclei by measuring their centroids, and the centers of their host galaxies by measuring the centers of the isophotes computed at $`1^{\prime \prime }`$ (the precise radius was slightly increased/decreased in some galaxies to avoid the effects of prominent structures). The elongation of the nuclei was established by visual inspection.

## 3 Results and Discussion

### 3.1 The Central Regions of Spirals in the NIR

A large structural complexity is observed in the $`H`$ band in the central regions of many spiral galaxies; on the other hand, it is clear that the NICMOS images provide a much cleaner picture of the underlying galaxy structure as compared to the heavily dust-obscured WFPC2 visual images. In the extreme cases of NGC 3067, NGC 4750 and NGC 6239, the NICMOS images unveil the location of the true galaxy center, which we had been unable to identify in the WFPC2 images due to strong dust features and star formation. The NICMOS images have revealed in two of these galaxies (NGC 4750 and NGC 6239) the presence of a distinct compact nucleus, which subsequently, with the acquired knowledge of its precise location, we have been able to identify and parametrize also in the WFPC2 images. The optical measurements for these nuclei were not reported in our previous papers; these measurements are reported here in Table 2, identified by the “$``$” symbol. The possibility of better “penetrating” the dust in the NIR with respect to the optical data has also led to the identification, based on the visual inspection of the $`H`$ images, of nuclear bar-like features in ESO499G37, NGC 4384, NGC 4750, NGC 5448, NGC 5806, NGC 7280 (identified by ‘S’ —standing for small— in the column ‘Bar’ in Table 3), and possibly also in IC4390, IC 5271, IC 5273, NGC 289, NGC 2339, NGC 3455, NGC 3898, NGC 3949, NGC 4102, NGC 4219, NGC 4806, NGC 5678, NGC 6810, NGC 7177, NGC 7188 and NGC 7690 (these ambiguous cases are identified by ’S?’ in Table 3). In IC 4390, IC 5271, NGC 3455, NGC 3949, NGC 4219, and NGC 5448 the nuclear bar-like features become evident in the $`VH`$ or $`JH`$ maps by their distinct colors (either redder or bluer than the surroundings). Consistent with the studies which have looked for nuclear bars in the context of the feeding of an active non-thermal nucleus, we find that even in the less dust-sensitive NIR images the nuclear bar-like features are not ubiquitous. Interestingly, spiral structure down to the center/nucleus is observable in the $`H`$-band images of several systems. In some galaxies, ESO205G7, NGC3928, NGC5678, a small-scale ring is observed around the centre/nucleus; these rings are likely not an artefact of PSF mismatch (since changing the sizes of the PSFs by 10% does not affect the result), and are likely due to circum-nuclear dust or star formation.

In summary, the NIR images, coupled with the optical images, show that there is no galaxy in our sample that is not “polluted” either by dust or by recent star formation or both in its nuclear/circum-nuclear region: all galaxies contain some nuclear/circum-nuclear feature overlying the smooth/old stellar populations. Excluding a few objects for which the signal-to-noise ratio of the NIR images is too low to accurately measure the underlying light outside the bright nucleus, the galaxies in our sample can be broadly grouped into a four classes, based on their nuclear NIR- and optical/NIR-color morphology:

Class 1. Galaxies with Concentrated Nuclear Star Formation Mixed with Dust. Judging from the $`H`$ images and the optical-NIR color maps, many galaxies (ESO 498 G5, ESO 499 G37, IC 4390, IC 5273, NGC 406, NGC 986, NGC 1483, NGC 1688, NGC 1800, NGC 2082, NGC 2566, NGC 2964, NGC 3067, NGC 3928, NGC 4102, NGC 4219, NGC 4384, NGC 5678, NGC 5806, NGC 6000, NGC 6329, NGC 6951, NGC 7217, NGC 7259) host in their innermost few-hundred pc2 concentrated star formation in the form of compact blue knots, intermixed with strong red dust features. With the exception of NGC 2566, NGC 3928, NGC 7217, and NGC 7259, where the star forming and dust features have a regular spiral morphology, and of NGC 6951, where they form a circum-nuclear ring, in all the remaining objects the nuclear star formation and dust have a rather irregular morphology. Only some of these systems have been found to host a nuclear bar or a central nucleus, so that the presence of these features is not a unique characteristic of this family of objects. All but four (ESO 499 G37, IC 5273, NGC 1688, NGC 1800) of these galaxies are early-to-intermediate type spirals in the Hubble range between Sab and Sbc. From these data it seems that star formation in the central regions of early-type spiral galaxies does not occur in the rare, transient, and intense “starburst mode”: low-level star formation appears to be common in these systems. A question to answer is whether this low-level star formation represents the low-luminosity end of the starburst phenomenon, or rather the fading evolutionary descendant of a major starburst itself (Heckman 1997).

Class 2. Galaxies with Diffuse Blue Nuclear Regions. Several systems, namely ESO 240 G12, NGC 972, NGC 1345, NGC 1398, NGC 2339, NGC 2748, NGC 3949, NGC 4527, NGC 5188, NGC 5377, NGC 5448, NGC 6810, NGC 7177, and NGC 7690 also host significant amounts of circum-nuclear dust mixed with rather blue surroundings; the latter however have a more diffuse appearance than the knot-like features typical of the previous group. For example, in NGC 5377 the blue region is embedded in a completely-obscuring dust torus and has the morphology of regular disk surrounding the distinct central nucleus. Among these systems there are a few (NGC 5377, NGC 5448, NGC 6810) which are known to host an AGN; for them, it is possible that the blue regions are related to non-thermal activity from the central engine.

Class 3. Galaxies with Regular Nuclear/Circum-Nuclear Dust. There are some galaxies that do not show evidence for star formation (with the possible exclusion of the nucleus), but host in their central regions strong dust lanes with a regular, well-defined spiral or ring morphology. These systems are ESO 205 G7, NGC 488, NGC 772, NGC 2460, NGC 3259, NGC 3277, NGC 4750, NGC 6384, NGC 7013. Some of these objects show some ionized-gas emission on ground-based scales (Table 1), but this must be distributed smoothly in the central region not to show up clearly in the images in the form of identifiable structures. Apart from being all early-to-intermediate type systems (S0a to Sbc), these galaxies show no other characterizing feature: some belong to the group of Figure 1 and thus host a distinct nucleus, others belong to the group of Figure 2; some have an exponential bulge, others host instead an $`R^{1/4}`$ bulge. The nuclei which are embedded in some of the galaxies in this class tend to have relatively red $`VH`$ colors, ranging from 1.3 to 2.1 magnitudes.

Class 4. Galaxies with Irregular Nuclear/Circum-Nuclear Dust. Several more objects in the sample also show a quiescent, non-starforming morphology with strong dust features, but the latter are distributed in a much more irregular and patchy way (although in some cases, in addition to the irregular patches, spiral-like dust is also present). These objects are IC 5271, NGC 289, NGC 1892, NGC 2196, NGC 2344, NGC 2758, NGC 3455, NGC 3898, NGC 3900, NGC 4806, NGC 4980, NGC 5985, NGC 6340, NGC 7162, NGC 7188, NGC 7280, NGC 7421, NGC 7513. Similar to the case above, some of these systems host a central distinct nucleus, others do not; with the exception of NGC 1892, they are early-to-intermediate type systems, with Hubble types between S0 and Sc. The $`VH`$ colors of the four nuclei detected in this class of systems range from 1.4 to 3.1 magnitudes; the latter high $`VH`$ value clearly indicates the presence of a significant amount of dust even on the nucleus itself. In several of these systems the presence of a nucleus is uncertain, due to the dust which obscures the central region.

Additional information is concisely reported in Table 3 where we summarize, object by object, the observed properties of the circumnuclear $`VH`$ or $`JH`$ light distribution.

### 3.2 The Nuclei

Table 2 lists the F160W sizes and luminosities for the nuclei detected in the NICMOS sample, measured as described in §2.2. When available, the sizes and luminosities of the nuclei in the F110W passband are also listed, together with, for an easy comparison, the F606W measurements published in Carollo et al. (1997, 1998) for the NICMOS galaxies observed within our program also with WFPC2. The F606W properties of the nuclei for the 8 NICMOS targets which were not included in our WFPC2 papers due to the late acquisition of the WFPC2 images (§2.2) are also given in the Table. The comparison of the sizes derived from the WFPC2 and F110W images with those derived from the F160W images is shown in Figure 4. The radii measured with NICMOS and WFCP2 are in good agreement within the error bars.

Compact star-cluster-type nuclei in bulgeless, weakly- or non-active late-type disks have been the subject of a few detailed studies (Kormendy & McClure 1993; Phillips et al. 1996; Matthews & Gallagher 1997; Lauer et al. 1998). In the visual band, these late-type nuclei have been found to cover a very narrow range in luminosity ($`M_I11\pm 1`$) despite the large range in luminosity of their host galaxies (from $`M_V=15.1`$ to $`M_V=18.2`$; Matthews et al. 1999). Our WFPC2 and NICMOS surveys show that photometrically-distinct nuclei are very frequent also in intermediate-type spirals. In contrast with the four nuclei in our sample embedded in the late-type galaxies, which lie typically ‘naked’ in the centers of the disks (as reported also by Matthews et al. 1999), in the earlier-types the nucleus is embedded either in a smooth pseudo-bulge with an exponential radial light profile, or, in the case of the brightest nuclei ($`12>M_V>16`$), in a complex structure of strong dust lanes and star-forming knots in irregular, spiral or ring-like patters typically still observable also in the NIR (see also Carollo et al. 1997; Carollo et al. 1998; Devereux et al. 1997). Furthermore, the new data show that:

(i) Within our sample, in all cases for which a nucleus was detected in the WFPC2 images and for which the NICMOS data are available, the nucleus is also clearly seen in the NIR: this suggests that all nuclei, independent of Hubble type, are relatively massive, given that the NIR light around 1.6$`\mu `$ is a rather good tracer of the stellar mass. In the NIR as much as in the visual passband, however, the nuclei embedded in the bright early- to intermediate-type hosts span a much larger magnitude range than the nuclei embedded in the fainter late-type disks, as shown in Figure 5 a, where the absolute $`H_{AB}`$ magnitude of the nuclei is plotted against their half-light radius derived in the same passband: The four nuclei embedded in the Scd and later-type disks have nuclear properties quite similar to those reported in the quoted literature for other late-type nuclei, and are therefore likely representative of that category of objects.

(ii) In our sample we find nuclei that are either embedded (e.g., NGC 6000) or not embedded (e.g., NGC 2964) in large-scale ( $`>`$$``$ 1 kpc) bars. We also find that some of the nuclei sit within smaller bars, with scales of up to a few 100 pc (e.g., ESO499G37). These small-scale bars could be particularly important within the context of bar-driven secular evolution scenarios for bulge formation, since they might provide a mechanism for bringing the disk material from the outer regions down to the pc scales in the shallow gravitational potential of the progenitor disks (Shlosman 1994; Shlosman & Robinson 1995; Shlosman, Begelman & Frank 1990; Friedli & Martinet 1993). The nucleus/nuclear-bar connection is however not a straightforward one since, as stressed above, nuclear bars are not a universal feature of nucleated spiral galaxies.

(iii) The brightest nuclei embedded in the early/intermediate-type galaxies tend to have larger radii than the fainter nuclei. An increase of radius with increasing luminosity of the compact nucleus in the early/intermediate types is detected in the $`H`$ band (linear correlation coefficient of -0.76, corresponding to a probability $`>99.9\%`$). This was already noticed in the F606W passband (see Figure 5 b, where the absolute magnitude versus half-light radius relation for the nuclei is given for the $`V`$ band, for a comparison with Figure 5 a, where the same relation is given the $`H`$ band). The same finding in the NIR, however, strengthens the reliability of this result, given the much simpler circumnuclear background structure in the $`H`$ band as compared to the complex one in which (especially the brightest) nuclei were found embedded in the visual images. The measurements of the half-light radii of the nuclei have a rather large uncertainty, mostly because at the distance of our sources the extended light component of the nucleus must be separated from the wings of the PSF (see also Matthews et al. 1999). The galaxies of our sample are however restricted to $`cz<2500`$ km/s, which ensures a rather high degree of homogeneity in the derived set of measurements. In contrast with the simulations, which showed a modest brightening of the nucleus corresponding to a modest increase in radius, the nuclei observed in the centers of early/intermediate-type spirals span an entire decade in radius, over which they show an increase in luminosity of about eight magnitudes. It is therefore likely that a substantial part of this correlation is driven by a real physical effect.

(iv) The resolution of the data is not appropriate to detect color variations within individual nuclei, but it is adequate to compare the nuclei to their surroundings. The early/intermediate-type spirals’ nuclei are found to cover a large range of colors, from $`VH`$ -0.5 to 3. Previous HST studies on the colors of the central regions of spirals report very red (dusty) centers for the early-type bulges (Peletier at al. 1999), and rather blue colors for the distinct nuclei of late-type disks (Matthews et al. 1999). The few late-type nuclei in our sample show a slight tendency for having bluer colors than the nuclei of the earlier-type hosts, as is shown in Figure 6, where the $`VH`$ color of the nuclei is plotted against the RC3 Hubble type of the host galaxy. Figure 6 also shows however that the spread in $`VH`$ color is rather large within the early/intermediate types themselves: The early/intermediate nuclei have an average $`VH1.44`$ with a sigma of about 0.82 magnitudes. In the Figure, the nuclei are identified with different symbols according to the galaxy-class to which their hosts belong (among the four described above; symbols are explained in the figure caption). No significant trend is observed with galaxy-class, with only a possible hint for the nuclei embedded in galaxies with regular nuclear/circum-nuclear dust to be slightly redder than the average. A comparison of each nucleus with its own galactic surroundings shows that the nuclei are in some cases redder and in others bluer than their hosts (see Table 3, where we also provide a summary of the properties of the nuclei in the $`H`$ band, and of their optical-NIR or NIR-NIR color properties). Statistically, the distribution of $`VH`$ colors is broader for the nuclei than for the surrounding galactic structure, as seen in Figure 7 which compares the $`VH`$ color distribution for the nuclei with the $`VH`$ color of the host galaxies. This Figure contains only those galaxy-nucleus pairs for which an isophotal fit to the underlying light of the host galaxies in both the $`V`$ and the $`H`$ band could be performed, so that their $`VH`$ color could be derived by integrating the $`V`$ and $`H`$ radial light profiles within 3′′ (see Carollo & Stiavelli 1998 and paper II for the $`V`$ and $`H`$ fits, respectively). The so-obtained colors are much less affected by the effects of dust or recent star formation with respect to aperture measurements. No isophotal fit could be performed for the hosts of some of the nuclei of Figure 6, and in particular to the hosts of the reddest and bluest nuclei; this explains the smaller range in $`VH`$ of Figure 7 with respect to Figure 6. Even so, the Kolmogorov-Smirnov, Kendall and Spearman tests give a probability of 0.996, 0.999 and 0.9997 that the two distributions are statistically different. This likely indicates ubiquitous presence of either ‘activity’ or dust reddening (or both) in the nuclei population. Figure 8 a shows the $`V`$ absolute magnitude versus the $`VH`$ color for the nuclei; Figure 8 b shows their $`JH`$ versus $`VH`$ color-color diagram. Different symbols in these Figures (are as in Figure 6 and) are used to identify nuclei embedded in hosts belonging to a different galaxy-class among the four discussed above. A more detailed discussion of these color-magnitude and color-color diagrams, and in particular of the constraints that they put on the stellar population properties of the nuclei, is given in Carollo et al. (2001); what we remark here is that there is no obvious trends for any ‘segregation’ on these diagrams among nuclei embedded in hosts with different circum-nuclear properties.

(v) Within the caveats imposed by the angular resolution of our data and the distance of our targets, no significant difference is found with Hubble type concerning the exact location of the nuclei with respect to the underlying galactic structure, and their projected shape, round versus elongated. In Table 3 we report, for each galaxy, whether the compact nucleus is centered on, or offset from, the galaxy isophotal center, and whether the nucleus and more generally the nuclear region is round or elongated. There are several cases where our NICMOS data suggest that the nuclei are not located precisely at the isophotal centers of the host galaxies, as found e.g. for M33 (de Vaucouleurs & Freeman 1970; Minniti et al. 1993). There are also a few cases where the nucleus shows some degree of elongation, as also found for, e.g., the nucleus of M33 (Kormendy & McClure 1993; Lauer et al. 1998; Matthews et al. 1999). With the current data, the displacement and elongation, however, appear not to be correlated with either the luminosity and color of the nucleus, or the galaxy type. It remains an open issue whether some nuclei are disk-like features, and whether the nuclei of galaxies with an AGN are intrinsically different from those of non-AGN galaxies. On the other hand, the fact that many of the best-resolved nuclei are typically round structures in projection supports the working hypothesis that (at least in non-AGN galaxies and) for all Hubble types these nuclei are typically cluster-like systems.

## 4 Summary

We have presented the analysis of HST NICMOS data in the F160W ($`H`$) filter for the central regions of a sample of 78 mostly early-to-intermediate type spiral galaxies. For 58 of these systems we have used WFPC2 F606W ($`V`$) images from our previous study to investigate their $`VH`$ nuclear colors. For 11 systems we have used additional NICMOS images in the F110W filter to investigate their $`JH`$ nuclear colors. We have focused our attention on the structural and color properties of the photometrically-distinct nuclei that (similarly to what we had previously found in the visual band) dominate the nuclear light in a large fraction of the sample in the NIR, and on their galactic circum-nuclear regions.

Although the nuclear complexity that was unveiled in many of these galaxies by the WFPC2 images is often still present in the $`H`$ band, the latter provides a cleaner view of the nuclear galactic structure, and allows to see through areas of recent star formation and through dust patches typically present in the central regions of these system. This strengthens the findings which we had previously suggested on the basis of our HST optical images, and adds further clues for understanding the nature and origin of galactic nuclei. In particular, several nuclear bar-like features have been detected in the $`H`$-band, including some which have a $`VH`$ or $`JH`$ color distinct from that of their surroundings. Nuclear bars are however not ubiquitous in this kind of systems. Any relation between the presence of a nucleus and of a barred component, if present, is not a simple one, as nuclei are found in both barred and non-barred hosts. In the $`H`$ band, consistent with what was found in the $`V`$ band, the nuclei embedded in the bright early- to intermediate-type hosts span a much larger magnitude range (about eight magnitudes) peak-to-peak than the nuclei embedded in the fainter late-type disks. In the $`H`$ band, the brighter the nucleus the larger its radius, a correlation that was found in the $`V`$ band and is confirmed by the NIR data. Even excluding the few bluer-than-average late-type nuclei in our sample, the nuclei of early- to intermediate-type spirals span a large range in $`VH`$ colors, some being significantly redder (and very likely dust-reddened), others being significantly bluer (‘polluted’ either by a non-thermal component or by recent star formation), than their surrounding circumnuclear regions. Several nuclei embedded in the early- to intermediate-type hosts are offset (typically by a few tens of pc) from the isophotal centers of the host galaxies; several are flattened structures; however, most of the nuclei, at HST resolution, are located at the galaxy centers, and appear to be round, cluster-like structures.

A quantitative analysis of the circum-nuclear galaxy structure and strength of the nuclear stellar cusp slopes in the NIR is reported in paper II.

We thank Roeland van der Marel for providing the pedestal-removal software, and the anonimous referee for his/her careful review that has helped us improve the presentation of our results. This research has been partially funded by Grants GO-06359.01-95A and GO-07331.02-96A awarded by STScI, and has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, Caltech, under contract with NASA.
Appendix A. NICMOS $`H`$ Band Imaging for Nine Galaxies  In Figure 9 we present the NICMOS F160W images of the nine galaxies in our sample for which neither WFPC2 F606W nor F110W images are available. The data reduction was done as described in §2. NGC1365 (a well-known AGN), NGC2903, NGC3031, NGC4569, and NGC6217 host photometrically-distinct nuclei; NGC 4569 also hosts a nuclear disk. The presence of a nucleus is uncertain in NGC4536 and NGC4565; no nucleus is present in NGC7479 and NGC7742. Table 4 summarizes some properties of the nine galaxies, including the $`H`$ measurements of sizes and magnitudes for the nuclei.
Appendix B. Assessment of Systematic Effects  Measuring the half-light radius and luminosity of a distinct nucleus embedded in a galaxy core entails assumptions on how to partition the observed light between the underlying galaxy and the nucleus. The different assumptions that we have used in §2.2, while all yielding acceptable decompositions, provide however different results, since it is not easy to establish where the nucleus ends and the galaxy begins. The situation is therefore one where the systematic error arising from separating the galaxy light from the nucleus light outweighs the error on the measurement of the radius and the luminosity of the nucleus itself. To assess the systematic errors on these quantities, we performed two sets of tests. The results for the F160W filter are described below. A similar analysis was performed for the F110W and the WFPC2 F606W filters, reaching similar conclusion as for the F160W filter. First, we tested the sensitivity of our nucleus-extraction procedure in the region of parameter space covered by those nuclei embedded in hosts for which an isophotal fit could be performed, and whose nuclear surface brightness profiles could be modeled with a Nuker-law (Lauer et al. 1995; see paper II for the description of our isophotal and Nuker fits). The Nuker law reduces to a power law when the characteristic radius $`r_b`$ is much larger than the nucleus’ radius. We built model galaxy images starting from the best-fit analytical descriptions of the light underlying the nuclei in these galaxies, and added to each of these model galaxies a contribution from a nucleus of size and luminosity equal to those measured with our procedure in that specific system. Two different analytical forms were adopted to describe the nuclei, a “generalized” modified-Hubble law, which reads:
$$I(R)=\frac{L}{2\pi b^2}\left(1+\frac{R^2}{b^2}\right)^{5/2}$$
(1)
and a Plummer law, which reads:

$$I(R)=\frac{L}{b^2}\left(1+\frac{R^2}{b^2}\right)^2.$$
(2)
In both Eq. (1) and (2) $`b`$ is the scalelength and $`L`$ is the total luminosity. We then convolved the model galaxy-plus-nucleus images with the NICMOS PSF in the appropriate filter (obtained with TINYTIM, as done for the data), and applied our nucleus-extraction procedure to these PSF-convolved model-images. The radius was corrected for PSF broadening in the same way as was done for the actual data. This test provides an estimate of how accurately we can disentangle the nuclei from their hosts for galaxy-plus-nucleus pairs of the kind that we have detected (similar sizes/luminosities of the clusters; similar contrasts with the underlying/surrounding galaxy).

Figures 10 a and b show, for these galaxy-nucleus pairs, the distribution of the ratio between the recovered flux and the input flux, and the distribution of the ratio between the recovered radius and the input radius, respectively. Figure 10 c shows $`\mathrm{\Delta }mag=(m_{recover}m_{input})`$, the difference between the recovered and the input magnitude, as a function of the logarithm of the ratio between recovered and input radius (after correction for systematic shifts). The figures show that 90% of the measurements lie within $``$0.3 magnitudes from the expected values, and that our extraction procedure applied to the intrinsic nuclear light profiles systematically underestimates the flux on average by $``$0.4 magnitudes, with a scatter of $``$0.2 magnitudes. By contrast, the mean radius appears to be overestimated in average by $``$45%, with a scatter of $``$25%; the 50%-error includes $``$80% of the points. A weak correlation between the derived radius and magnitude of the nuclei is detected, but the maximum absolute variations are rather modest, about a factor 3 increase in radius for a $``$1 magnitude brightening of the nucleus.

Secondly, we performed another set of tests which covered a larger area of parameter space to assess the degree of degeneracy in the parameters provided by our galaxy-nucleus decomposition. These additional tests allow us to estimate how much and often intrinsically different galaxy-nucleus pairs can result in similar measurements when analysed with our procedure. We built a grid of model galaxies described by either a power-law cusp of slope $`\gamma `$ ranging from 0 to 0.7, or an exponential profile, and we added to these model galaxies the contribution from a central nucleus, described by a Plummer law. The nucleus was assumed to have $`H_{AB}=19.6`$; the surface brightness of the underlying galaxy was taken similar to the typical one observed for real objects. In particular, a mean $`H_{AB}`$ surface brightness within $`20^{\prime \prime }`$ equal to 19.7 mag/arcsec2 was considered for the power-law galaxies, corresponding to a central surface brightness after convolution with the PSF of 19.2, 18.4, 17.4 or 16.4 mag/arcsec2 for $`\gamma =0.1`$, 0.3, 0.5 or 0.7, respectively. For the two highest $`\gamma `$ values of 0.5 and 0.7, also a mean surface brightness within $`20^{\prime \prime }`$ equal to 20.5 mag/arcsec2 was considered, corresponding to a central surface brightness after convolution with the PSF of 18.2 or 17.2 mag/arcsec2, respectively. A central surface brightness equal to 18.5 or 19.3 mag/arcsec2 was considered instead for the exponential galaxies. We convolved this set of galaxy-nucleus models with the appropriate PSF, and applied our procedure to recover the input parameters.

Figure 11 a shows the $`\mathrm{\Delta }mag`$ (the difference between the recovered and input $`H_{AB}`$ magnitude of the nucleus) versus the nuclear cusp slope $`\gamma `$ for the power-law galaxies, and versus the exponential scale-length for the exponential galaxies. Figure 11 b shows the recovered- versus the input-FWHM of the central nucleus. Notice how difficult it is to recover the radius of the nucleus when a large nucleus lies on top of an underlying steep cusp. Figure 11 c shows the $`\mathrm{\Delta }mag`$, the difference between the recovered and the input magnitude, versus the recovered radius of the nucleus. The results of this set of tests demonstrate that when all the galaxy-nucleus pair models are taken into account, the flux is underestimated by $``$0.27 magnitudes, with a scatter of 0.60 magnitudes; 80% of the points lie within the 50% flux error. The radius appears to be overestimated by $``$30% with an rms of $``$40% and with the 50%-error including 90% of the points. Also these general models - especially the power-law ones - show a weak trend between derived radius and magnitude of the nuclei. The observed maximum variation is again about a factor 3 of increase in radius, which in this case corresponds to $`1.5`$ magnitude brightening of the nucleus.

Summarizing, both the tests performed on specific galaxy-nucleus pairs and the more general ones indicate that: (1) the reported fluxes are likely underestimated by $``$0.3 magnitudes, and are accurate to better than 0.6 magnitudes; (2) the reported radii are likely overestimated by about 40% and are accurate to within 40%; (3) the galaxy-nucleus degeneracy and possibly the adopted extraction technique may produce a mild trend for brighter nuclei to be slightly more extended than fainter nuclei, although it is unlikely that such an artifact could produce a correlation between luminosity and size of the nuclei involving several magnitudes of difference in luminosity.

We finally checked whether there was any major systematic effect affecting our measurements depending on the luminosity of the nucleus embedded in a given galaxy. We computed an additional set of about 200 galaxy-models with a central surface brightness of 18.4 mag arcsec-2 and nuclei with a range of luminosities; we used both power-laws and exponentials models with cusps and exponential lengths covering the same intervals as in the models of Figure 10. While the systematic error in the measurement increases, there is no systematic trend on the recovered luminosity of the nucleus. This is illustrated in Figures 12 a and b where we plot, as a function of the input magnitude of the nucleus, the difference between the recovered and the input magnitudes for both power-law and exponential galaxy-models, respectively. For each luminosity bin the plot shows the mean shift and the error; the latter is dominated by systematic effects.

Although these experiments do provide an estimate of the size of our systematic errors, we have chosen not to apply these “corrections” to the radii and luminosities of the nuclei in our analysis, since: (i) the absolute amount of the systematic shifts is not such to affect in any significant way our conclusions, and (ii) the estimated correction would be strictly valid if all nuclei were “homologous”, i.e., had similar intrinsic profiles, and could instead introduce additional scatter if the intrinsic light profile changed from nucleus to nucleus.

The combined results for the WFPC2 and NICMOS tests were used to assess the uncertainties associated with the colors of the nuclei. On the assumption that nuclei and underlying galaxies have similar intrinsic light profiles in the visual and NIR and do not have large unresolved color gradients, the same analytical models would properly describe both the $`V`$ and the $`H`$ band data. In this case one expects the error on the $`VH`$ color of the nuclei to be much smaller than the total formal error obtained by quadratically adding the errors on each band. We checked this assumption on the galaxy-nucleus pairs of Figures 10 a and 10 b. In Figure 13 we show the histogram of the error on the $`VH`$ color obtained adopting for the nuclei a modified Hubble profile (solid line) or a Plummer law (dotted line): The figure shows that although the systematic effect is not entirely cancelled, the resulting error is reduced to 0.4 magnitudes instead of the formal 0.8 magnitudes.


format: text

text:
# A novel approach for extracting time-delays from lightcurves of lensed quasar images

## 1 Introduction

The time delay between light rays from gravitationally lensed quasar images is a measurable parameter directly related to the gravitational potential and to the Hubble constant $`H_0`$ (Refsdal Refsdal (1964)). Accurate time delays obtained from multiply lensed QSOs can hence be used i) to determine $`H_0`$ providing that the lens mass distribution is known, or ii) to constrain the mass distribution in a given lens, once $`H_0`$ has been determined from other methods (or other lensed QSO systems). Much effort has therefore been devoted to observations of lensed QSOs during the last years, and in particular to the photometric monitoring of lensed QSO images.

Measuring the time delay between the images of lensed quasars is, for several reasons, not a trivial task. First it requires regular monitoring of a target over a long period (substantially longer than the time delay). Second, the sampling is crucial and has to be determined on basis of the intrinsic variations in the quasar and the estimated time delay. Third, most objects are not observable during the whole year, i.e., they will be below the horizon for certain periods. Fourth, there are nights with bad weather conditions when no data are obtained. Finally, variations in lensed quasars may not only be due to intrinsic fluctuations of the quasar, but also to microlensing by compact objects along the line of sight. Such a microlensing signal, depending on its timescale and amplitude, can be used to constrain the size of the continuum and the line emitting regions in the quasar, and the distribution of compact matter in the lens galaxy (Paczyński Pacz (1986), Kayser et al. Kayser (1986)). However, as long as these external variations are not clearly distinguished from the intrinsic variations, microlensing remains a nuisance for time delay measurements. For the above reasons, advanced statistical methods have to be used to measure time delays from quasar light curves.

We have developed a method based on $`\chi ^2`$ minimization of the data and a numerical model light curve. The aim was to develop a method based on simple principles, and to be able to properly detect and correct for possible microlensing effects. It has already been successfully applied to several time delay measurements (Burud et al. Burud (2000); Hjorth et al. Hjorth (2001)). In this paper we will present the principles of the algorithm (Sect. 2) apply it to various simulated light curves (Sect. 3) and to a public dataset of QSO 0957+561 from Serra-Ricart et al. (Serra-Ricart (1999)) (Sect. 4).

## 2 Method

Let us assume that we have light curves for two images, A and B, of a lensed quasar. There are N data points in each of the time dependent light curves $`a(t)`$ and $`b(t)`$ with measurement errors $`\sigma _a`$ and $`\sigma _b`$. The two curves are identical except for a shift in time, $`\mathrm{\Delta }t`$, and in magnitude, $`\mathrm{\Delta }m`$. It is therefore possible to model the two curves with one model curve $`g(t)`$, and the parameters representing the time delay and the magnitude shift. In some cases the light from one or both of the images may be microlensed by individual stars in the lensing galaxy. Except for the cases of high amplification events, which are of short duration, microlensing effects are often slow variations that can be modeled to first order as a linear variation with slope $`\alpha `$.

As model curve we choose an arbitrary curve with a fixed number of equally spaced sampling points, $`M`$. This model may be $`\chi ^2`$ minimized to the two observed light curves. The minimization is done only for the N observed data points. In this way only the model curve is interpolated and not the data. We thus minimize the following function:

$``$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left[{\displaystyle \frac{a(t_i)g(t_i)}{\sigma _{a_i}}}\right]^2`$ (1)
$`+{\displaystyle \underset{i=1}{\overset{N}{}}}\left[{\displaystyle \frac{(b(t_i\mathrm{\Delta }t)(\mathrm{\Delta }m+\alpha t_i))g(t_i)}{\sigma _{b_i}}}\right]^2`$
The determination of the time delay from sampled light curves always relies on the implicit assumption that intrinsic variations are continuous in time and slow enough to be measured, given the adopted frequency of the observations. Consequently with this assumption, if the typical sampling interval is $`\tau _1`$, we can thus smooth the model curve $`g(t)`$ on the same time scale $`\tau _1`$ (e.g., typically 7 days when observations are obtained once per week). This is done by introducing a smoothing term which minimizes the square of the difference between the original model and the same model convolved with a Gaussian $`r`$ whose Full Width at Half Maximum (FWHM) is $`\tau _1`$. This smoothing is performed for all the points $`M`$ on the model curve $`g(t)`$ so that the dates without data points are smoothed out on the model curve and gaps in the observed light curves will not contribute significantly to the result. The smoothing term is multiplied by a Lagrange parameter $`\lambda `$. This parameter can be chosen so that the model curve matches the data correctly in a statistical sense (i.e., we get the correct $`\chi ^2`$). In practice however, excessive smoothing will prevent the minimization from converging. Therefore the smoothest possible solution which allows the minimization to converge is usually adopted.

In a $`\chi ^2`$ minimization it is assumed that all the data points are independent. However, a light curve often consists in groups of nearby points and more isolated points. In regions where the time sampling is much shorter than the intrinsic quasar variations, the different measurements do not really bring new information on the shape of the light curve. Rather, they increase the precision with which the quasar magnitude is known at that moment. They might even be affected by common systematic errors, i.e. coming from a microlensing variation with a time scale comparable to the interval covered by these neighbouring points.

In such cases of strongly uneven time coverage, we find that the best results are often obtained when a weight $`W_i`$ is given to each data point in the following manner. For a point at a time $`t_i`$ we calculate the relative distance $`t_it_j`$ to all the other points in the curve. A Gaussian function with a FWHM$`=2\sqrt{2\mathrm{ln}(2)}\tau _2`$ chosen by the user is centered at $`t_i`$. The inverse of the sum of the ordinates of the Gaussian for each of the points $`t_j`$ will give the weight for the point $`t_i`$. The expression of $`W_i`$ can be written:

$$W_i=\frac{1}{_{j=1}^Nexp(\frac{t_it_j}{\tau _2})^2}$$
(2)
The weight is normalized so that the maximum allowed weight for a data point is $`1`$, and this will occur only when one point is present within the time interval defined by $`\tau _2`$. The ideal choice of $`\tau _2`$ is the approximate time scale of variations in the quasar. For some systems this time scale will be long, e.g., several hundred days, for other systems it can be shorter e.g., 50 days. We can now write the final function to be minimized:

$``$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}W_i\left[{\displaystyle \frac{a(t_i)g(t_i)}{\sigma _{a_i}}}\right]^2`$ (3)
$`+{\displaystyle \underset{i=1}{\overset{N}{}}}W_i\left[{\displaystyle \frac{(b(t_i\mathrm{\Delta }t)(\mathrm{\Delta }m+\alpha t))g(t_i)}{\sigma _{b_i}}}\right]^2`$
$`+\lambda {\displaystyle \underset{i=1}{\overset{M}{}}}\left[g(t_i)\left[rg\right](t_i)\right]^2`$
where $`rg`$ signifies a convolution between $`r`$ and $`g`$.

### 2.1 Microlensing effects and iterative version of the algorithm

Slow microlensing variations in one of the light curves are modeled as a linear term with the parameter $`\alpha `$. If higher order variations are present we can use the method in an iterative way. This is done by splitting the light curves into several parts and analyzing them separately. The $`\chi ^2`$ method is applied to determine the amount of linear external variations ($`\alpha `$) for each separate part and for a range of input time delay values. First we run the programme with the first time delay $`dt_1`$ in the chosen range and we estimate $`\alpha `$ in each of the separate parts. We then correct each part with the corresponding $`\alpha `$ and add the parts back together in order to obtain a set of microlensing-corrected light curves for the given input time delay $`dt_1`$. We can now run the $`\chi ^2`$ method to measure the time delay for these “corrected” curves. This procedure is repeated for all the input time delay values $`dt_2`$, $`dt_3`$, … $`dt_n`$ in the range to be studied. We finally obtain a time delay measurement for each of the $`n`$ input values, and all these measurements will generally converge towards the real time delay value.

## 3 Time delay measurement on simulated light curves

### 3.1 Data simulations

We apply the method to simulated light curves in order to understand and estimate the errors on the measured values. Various real quasar light curves were used as models for typical time variations and sampling of data points. The measurement errors in the curves represent simulated photon noise errors, so that the faint component will have larger errors than the bright one. We present the results from two sets of simulated light curves. The first set (1a, b and c in Table 1) consists of two curves, A and B, each of 60 points, dispersed over a time interval of 2 years (see Fig. 1 (top) and 2 (top)). The B curve is shifted by 145 days in time and has a magnitude offset of 1.95 mag. In set 1b we have added a slow linear time dependent variation with a slope $`\alpha `$ (Eq. 3) in the B curve in order to simulate long-time-scale microlensing effects (Fig. 1 middle). A higher order microlensing variation ($`\alpha t+\beta t^2`$) is added to the B curve in set 1c (Fig. 1 bottom). In the second set (2a, b, and c) we have simulated another kind of intrinsic variation and sampling. The curves span an interval of two years and they have 36 data points each. The time delay is 110 days and the magnitude offset is 0.67 mag (see Fig. 3 top). In sets 2b and c, we added linear (b) and higher order (c) variations in the A curve (see Fig. 3 middle and bottom).

### 3.2 Time delay measurements

The $`\chi ^2`$ algorithm was applied to all the simulated sets of light curves. The fit was performed several times in order to check for the influence of various parameters such as: (i) the number of data points in the model curve ($`M`$ in Eq. 3), (ii) the FWHM ($`\tau `$) of the Gaussian used to smooth the model curve ($`r(t)`$ in Eq. 3), (iii) the FWHM ($`2\sqrt{2ln(2)}\tau _2`$) of the Gaussian that defines the time scale of the variations in the curve (used to determine the weight $`W_i`$ in Eq. 3), and (iv) the Lagrange parameter ($`\lambda `$ in Eq. 3). For the first set the best results were obtained when including a weight ($`W_i`$ in Eq. 3) whereas for the second set the best results were obtained with no weight. This is because the data points in the second data set are more regularly distributed over time. The smoothing parameters were chosen so that the model curve is smooth over time, but still allowing variations large enough to fit the data. The results are shown in Table 1.

We performed Monte Carlo simulations to estimate the errors in the results. Sets of 1000 light curves were simulated with error bars determined randomly from a Gaussian distribution with a $`\sigma `$ equal to the measurement errors ($`\sigma _a`$ and $`\sigma _b`$). One set was made with the same sampling as the original curves (the simulated examples) and another set of 1000 curves was created with the same number of sampling points as the original, but randomly distributed.

When applying the iterative algorithm it is natural to think that the choice of the number of individual parts $`n`$ is important for the results. Various tests on our simulated curves indicate that if $`n`$ is too small, there is little difference from the total curve and higher order microlensing effects can not be properly corrected for. If on the other hand $`n`$ is too high, each part will contain very few data points and for some parts (depending on the sampling) the overlap between the two curves will be too small to determine the external variations. The choice of $`n`$ will thus depend on the number of data points and on how they are distributed (e.g., presence of gaps in the curves). The optimal values for the curves in the simulated curves, i.e., the values resulting in the best convergence of $`\mathrm{\Delta }t`$, were found to be 3 and 2 in sets 1 and 2 respectively. The results shown in Table 1 correspond to the mean and standard deviation of the values obtained.

### 3.3 Results

All the time delay measurements obtained in sets 1a, b and 2a, b are very well determined within the estimated errors (see Table 1). As expected, the results depend on the sampling of the light curves and the error estimates from the Monte Carlo simulations with random sampling are larger than the ones obtained with a fixed sampling. In sets 1b and 2b where a slow linear variation is added in one of the curves, the time delay is still well determined whereas in sets 1c and 2c where higher order variations are added there are small systematic errors in the time delay estimated with the direct method. In these cases, the iterative method gives more accurate results with smaller uncertainties. Fig. 2 displays how the microlensing effects in data set 1c are corrected for with the iterative method. In the sets with no microlensing (1a and 2a) we note that the iterative method yields larger errors than the direct method. When no external variations (simulated microlensing effects) are present the best constrained time delay value is obtained by setting $`\alpha =0`$, hence removing one of the free parameters. The iterative method should therefore preferably only be applied where external variations are clearly present.

The magnitude shifts between the two curves are very well determined in sets 2a, b and c whereas in example 1a, b and c there are small systematic errors of $`0.010.05`$ mag.

## 4 Application to light curves of QSO 0957+561

Our method was also applied to a public data set of QSO 0957+561 published by Serra-Ricart et al. (Serra-Ricart (1999)) (see Fig. 4). The R-band lightcurves contain 214 data points and span the time interval from February 1996 to July 1998. We removed eight points that deviate by more than 0.05 mag from neighbouring points in both components. We applied our $`\chi ^2`$ algorithm to the data and, as was done for the simulated data, we estimated the error bars using Monte Carlo simulations. A time delay estimate $`\mathrm{\Delta }t=423\pm 9`$ days and a magnitude shift $`\mathrm{\Delta }m=0.063\pm 0.007`$ were obtained. Considering the two gaps in the data set, JD 2450242 to 2450347 and JD 2450637 to 2450729, Serra-Ricart et al. have divided the data set into two parts, DS I and DS II. Since the second part (DS II) show the largest activity they have determined a time delay value on DS II only, corresponding to the autumn 1997/spring 1998 season. They measure a time delay value of 425$`\pm `$4 days and a magnitude shift of 0.06 mag on these curves. Measurements with two other methods, the dispersion spectra and discrete cross-correlation techniques give values of $`426\pm 12`$ and $`428\pm 8`$ days respectively on the same data set (Serra-Ricart et al. 1999). In order to compare the results we therefore also apply our method to DS II. On the curves in DS II we measure $`\mathrm{\Delta }t=424\pm 22`$. Both time delay values, the one measured on the entire data set and the one measured on DS II, are compatible with the results from Serra-Ricart et al. However our errors are larger than the ones found by Serra-Ricart et al. In particular we note that our errors are considerably larger when using only DS II than when using the entire data set. Serra-Ricart et al. do not give any time delay measurement from the total curve but their method seems to be well suited for time delay measurements on a continuous curve (DS II), even when the time span of observation is only of the order of the time delay value. With our method the time delay is better constrained with longer time series, even if gaps are present in the curves.

An additional reason for the differences in the results might be that whereas we have discarded 8 points form the total light curve, Serra-Ricart et al. have discarded 23 points. Given our limited knowledge of the data points we found no reason to exclude more points.

## 5 Conclusion

We have developed a method to determine time delays between light rays from lensed quasar images from their respective light curves. As other $`\chi ^2`$ based methods (e.g., Press et al. Press (1992) and Barkanna Barkanna (1997)) we assume that the statistics of the measurement errors are Gaussian, and we require smoothness for nearby points in the light curves. In our method we have included an optional weight to each of the data points, depending on the relative distance to neighbouring points. In this way one point isolated in time will receive more weight than each individual point in a cluster of nearby points.

Since microlensing effects are often encountered in real light curves of lensed quasars we have adapted our method in order to correct for such variations. A linear term in one of the curves is included to correct for slow long-term variations. Realistic simulations show that the time delay and the magnitude shift between two light curves are determined to within a few percent for cases with no or slow microlensing effects. For faster microlensing variations we find that running the algorithm in an iterative way yields better time delays. This confirms what was proposed in the time delay measurement of B1600+434 (Burud et al. Burud (2000)) for which four methods were applied: the minimum dispersion method (Pelt Pelt (1996)), the SOLA method (Pijpers Pijpers94 (1994), Pijpers (1997) ), the method described in the present paper and its iterative version. The iterative method yielded the best constrained time delay in this data set due to external variations in one of the components.

We also applied our method to a public data set of QSO 0957+561 and have shown that our results are in agreement with the published time delay.

###### Acknowledgements.

IB and SS are supported by contract Pôle d’Attraction Interuniversitaire, P4/05 (SSTC, Belgium). JH is supported by the Danish Natural Science Research Council (SNF).


text:
# 𝐵 meson 𝐵-parameters and the decay constant in two-flavor dynamical QCD Talk presented by N. Yamada.

## 1 Introduction

In recent lattice calculations, there have been several indications that the presence of sea quark increases the $`B`$ meson decay constant $`f_B`$ significantly, i.e. by about 10–15% . Then, the study of the sea quark effect on the $`B`$ parameters has a great phenomenological importance, because the relevant physical quantity to the $`B\overline{B}`$ mixing is the combination $`f_B^2B_B`$.

In this work we perform unquenched lattice calculations of both $`f_B`$ and $`B_B`$ using a consistent formulation. We use the nonperturbatively $`O(a)`$-improved Wilson quark action for both sea and valence light quarks, and the NRQCD action for heavy quark. Using the same dynamical QCD configurations, a search for sea quark effects in the light hadron spectrum is also being performed . Quenched calculations of $`f_B`$ and $`B_B`$ using the NRQCD action have been presented previously .

In full QCD, the effect of virtual pion loops leads to chiral logarithms, a specific dependence of physical quantities on the sea quark mass. Hence, it provides a non-trivial test of the sea quark effect in lattice simulations. For $`f_B`$ and $`B_B`$ chiral perturbation theory (ChPT) is extensively studied in  including its extension to the quenched and the partially quenched cases, where sea and valence quark masses can be different. We test the chiral logarithm using our numerical data.

## 2 Simulation details

We generate a set of unquenched gauge configurations with the plaquette gauge action and two-flavors of nonperturbatively $`O(a)`$-improved Wilson fermions . Simulations are performed on a $`20^3\times 48`$ lattice at $`\beta `$=5.2, which roughly corresponds to the lattice spacing $`a`$ 0.1 fm. We take five different sea quark masses covering $`m_{PS}/m_V`$=0.6–0.8. For comparison we also perform a quenched simulation at $`\beta `$=6.0 on the same lattice size. Further details of our simulations are given in .

We use the NRQCD action for heavy quark including all $`1/m_Q`$ corrections consistently in the action and operators. In the calculations of $`f_B`$ and $`B_B`$ we follow the methods adopted in the previous works . In particular, the $`B`$ parameters are computed with four different prescriptions, all of which have equivalent accuracy up to higher order contributions we neglect. We may estimate the systematic uncertainty taking the differences among them .

## 3 Decay constant

Figure 1 shows the chiral behavior of the heavy-light decay constant $`\mathrm{\Phi }_{f_B}f_B\sqrt{m_B}`$ with $`m_Qm_b`$. Both horizontal and vertical axes are normalized with a proper power of $`r_0`$ and converted to the physical unit using the lattice result for $`r_0m_\rho `$ combined with physical input for $`m_\rho `$. If we compare the unquenched result (filled triangles) with the quenched data (open diamonds), we clearly find the increase of about 15% with the unquenching as found previously

We consider the following three forms for the chiral extrapolation

$`\begin{array}{cc}\text{“lin”}:\hfill & a_0(1+a_1X),\hfill \\ \text{“qua”}:\hfill & b_0(1+b_1X+b_2X^2),\hfill \\ \text{“log”}:\hfill & c_0(1+c_1X+c_2X\mathrm{ln}X),\hfill \end{array}`$ (4)
with $`X=r_0^2m_\pi ^2`$. The “log” form is motivated by the ChPT result for heavy-light mesons . For two-flavor QCD, it predicts

$$r_0^2c_2=\frac{3}{4}\frac{1+3g^2}{(4\pi f)^2},$$
(5)
with $`f`$ the pion decay constant. The coupling $`g`$ describes $`B^{}B\pi `$ interaction, which may be evaluated from $`D^{}`$ decays using heavy quark symmetry. If we take a range $`g`$=0.2–0.7 rather conservatively covering $`D^{}D\pi `$ decay  and the recent $`D^{}`$ width experiment , then we obtain $`r_0^2c_2`$=$``$0.5(2) GeV-2.

We find that both “qua” and “log” fits describe our unquenched data well, while the linear fit gives somewhat poorer $`\chi ^2`$ than the other two. For the “log” fit we obtain $`r_0^2c_2`$=$``$0.4(3) GeV-2, which is consistent with the above phenomenological estimate.

On the other hand, in quenched QCD, no deviation from linearity is observed in our data, and the three fit forms lead to almost identical result in the chiral limit. This behavior contradicts the prediction of quenched ChPT (QChPT) , which gives even stronger logarithmic dependence. There are, however, many unknown parameters in QChPT, so it is still possible to find a parameter region which mimics the linear behavior in our light quark mass regime.

Since we are not able to decide the fit form with our current unquenched data, the uncertainty in the chiral limit is quite large. Figure 2 shows the $`1/M_P`$ dependence of $`\mathrm{\Phi }_{f_{B_d}}`$ at the physical light quark mass. At the physical $`B`$ meson mass, we observe that the two-flavor result (filled symbols) is 14–25% larger than the quenched result (open circle), depending on the fit form.

## 4 $`B`$ parameters

We also examine the chiral behavior of $`B_B`$ by applying the fit forms in (4), but, because the light quark mass dependence itself is quite small, we do not observe any significant dependence on the fit functions. A weak dependence on sea quark mass is also suggested by ChPT, which gives

$$r_0^2c_2=\frac{3}{2}\frac{13g^2}{(4\pi f)^2}=0.1(4)\text{GeV}^2$$
(6)
for $`B_B`$.

Figure 3 shows the $`1/M_P`$ dependence of $`B_B(m_b)`$ for both quenched (open symbols) and unquenched QCD (filled symbols). Four different symbols correspond to different prescriptions defined in . Within the uncertainty indicated by the different prescriptions, we find no significant effect of unquenching. $`B_{S_s}(m_b)`$ is also analyzed in a similar manner and we arrive at the same conclusion.

## 5 Conclusions

Our preliminary results in two-flavor QCD are summarized as follows,

$`f_{B_d}`$ $`=`$ $`190(14)(07)(19)\text{MeV},`$
$`f_{B_s}/f_{B_d}`$ $`=`$ $`1.184(26)(20)(15),`$
$`B_{B_d}(m_b)`$ $`=`$ $`0.872(39)(04)(73),`$
$`B_{B_s}/B_{B_d}`$ $`=`$ $`0.999(12)(04),`$
$`B_{S_s}(m_b)`$ $`=`$ $`0.858(33)(07)(72),`$
where errors are from statistics, chiral extrapolation, and systematics. We discard the linear fit from the final results because of its poor $`\chi ^2`$ value. The central value is obtained by taking an average over “qua” and “log” results, and their difference is taken as a systematic error from the chiral extrapolation. Some of the systematic errors in the $`B`$ parameters are estimated by the scatter with different prescriptions as in . Other uncertainties are estimated by a naive dimensional counting, assuming $`\mathrm{\Lambda }_{\mathrm{QCD}}`$=500 MeV and $`\alpha _s`$=$`\alpha _V(1/a)`$=0.267.

Our data for $`f_B`$ and $`B_B`$ are consistent with the prediction of ChPT. This is different from the light hadron sector where we find that ChPT fails to reproduce the lattice data for pion mass and decay constant . To investigate these issues we are currently accumulating further statistics, which also helps to reduce the statistical error in our results.

This work is supported by the Supercomputer Project No.66 (FY2001) of High Energy Accelerator Research Organization (KEK), and also in part by the Grant-in-Aid of the Ministry of Education (Nos. 10640246, 11640294, 12014202, 12640253, 12640279, 12740133, 13640260 and 13740169). K-I.I and N.Y are supported by the JSPS Research Fellowship.


text:
# Spin-Glass and Chiral-Glass Transitions in a ±𝐽 Heisenberg Spin-Glass Model in Three Dimensions

## Abstract

The three-dimensional $`\pm J`$ Heisenberg spin-glass model is investigated by the non-equilibrium relaxation method from the paramagnetic state. Finite-size effects in the non-equilibrium relaxation are analyzed, and the relaxation functions of the spin-glass susceptibility and the chiral-glass susceptibility in the infinite-size system are obtained. The finite-time scaling analysis gives the spin-glass transition at $`T_{\mathrm{sg}}/J=0.21_{0.02}^{+0.01}`$ and the chiral-glass transition at $`T_{\mathrm{cg}}/J=0.22_{0.03}^{+0.01}`$. The results suggest that both transitions occur simultaneously. The critical exponent of the spin-glass susceptibility is estimated as $`\gamma _{\mathrm{sg}}=1.7\pm 0.3`$, which makes an agreement with the experiments of the insulating and the canonical spin-glass materials.

It has been an unsettled issue whether a short-ranged spin-glass model can explain the spin-glass transition observed in real materials. The existence of the spin-glass transition in the Ising model has been widely recognized: various types of investigations  gave consistent values of the transition temperature. The recent large-scale simulation study  evaluated the critical exponents which quantitatively agree with the measurements of the Ising-type spin-glass material, Fe0.5Mn0.5TiO3. However, the spins and the magnetic interactions of many spin-glass materials are mostly isotropic, and the anisotropy effect can be regarded as a secondary effect. Therefore, the Heisenberg model with random interactions is considered to be the simplest theoretical model. The problem is whether this model exhibits the spin-glass transition or not. Previous numerical investigations suggest that the finite-temperature transition does not occur.

In order to solve this contradiction, Kawamura  introduced the chirality mechanism. The chirality is a variable whose absolute value corresponds to a volume of a parallelepiped spanned by three spins, and whose sign represents the handedness of the spin structure. In the chirality mechanism, the chirality is considered to freeze randomly at a finite temperature without the conventional spin-glass order. The experimentally-observed spin-glass anomaly is considered as an outcome of a mixture of the chirality with the spin by the weak magnetic anisotropies. The existence of the chiral-glass transition has been confirmed by several numerical analyses. However, it was pointed out recently by one of the present authors (S. E.) that the finite-temperature spin-glass transition may exist by examining the stiffness of the system when all the spins on a surface are rotated and reversed, respectively. The same evidence was also reported by a dynamical Monte Carlo simulation which removes the global spin rotation.

In this Letter, we focus on the existence or the non-existence of the spin-glass and the chiral-glass transitions. By employing the non-equilibrium relaxation method, we have revealed that both transitions may occur simultaneously at a finite temperature $`T/J0.21`$. This result possibly settles the discrepancies between the experiments and the theories in the Heisenberg spin glass in three dimensions.

In the simulational studies of the phase transitions, one obtains equilibrium quantities in finite-size systems, and the finite-size scaling is utilized to extract the thermodynamic properties from them. However, the dynamics of the simulations become very slow due to the randomness and frustration in the spin-glass models, which makes it hard to realize the equilibrium state in large systems. The system size accessible by the equilibrium simulations is restricted to linear sizes $`L30`$ up to now. In this situation, reliability of the finite-size scaling might become doubtful, which causes controversial conclusions of the existence  or the nonexistence  of the spin-glass transition depending on the sizes and the temperatures of the data used. In this Letter, we follow a completely alternative approach to the thermodynamic limit, i.e., we observe the relaxation of the infinite-size system to the equilibrium state. In order to extract the equilibrium properties, the finite-time scaling analysis is utilized instead of the conventional finite-size scaling. This approach is known as the non-equilibrium relaxation (NER) method.  Actually, we prepare a very large lattice with the paramagnetic state ( perfectly random spin configuration), and observe the relaxation of the susceptibility. The simulations are stopped before the finite-size effect appears. Time dependences of the susceptibility diverge algebraically as $`\chi (t)t^{\gamma /z\nu }`$  in the critical region, while it converges to a finite value in the paramagnetic region. We detect the phase transition by this difference. Another advantage of this strategy is an absence of a trivial exponential decay due to the global spin rotation. Since the spin autocorrelation function is essentially a one-body function about the spin, one may encounter this exponential decay even in the ordered phase after a time scale of the global spin rotation.

We consider the $`\pm J`$ Heisenberg spin-glass model on a simple cubic lattice of $`N=L\times L\times (L+1)`$ under the skew periodic boundary conditions, $`=_{i,j}J_{ij}𝑺_i𝑺_j.`$ The exchange interactions $`J_{ij}`$ take a value of $`+J`$ or $`J`$ with the same probability, and $`𝑺_i`$ is a three-component vector spin with $`|𝑺_i|=1`$. The summation runs over all the nearest-neighbor spin pairs. We evaluate the spin-glass susceptibility, $`\chi _{\mathrm{sg}}`$, and the chiral-glass susceptibility, $`\chi _{\mathrm{cg}}`$, through the following relations.

$`\chi _{\mathrm{sg}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i,j}{}}\left[𝑺_i𝑺_j^2\right]_\mathrm{c}=N{\displaystyle \underset{\mu ,\nu }{}}\left[q_{\mu ,\nu }^2\right]_\mathrm{c},`$ (1)
$`\chi _{\mathrm{cg}}`$ $`=`$ $`{\displaystyle \frac{1}{3N}}[({\displaystyle \underset{i,\mu }{}}C_{i,\mu }^{(a)}C_{i,\mu }^{(b)})^2]_\mathrm{c},`$ (2)
where $`q_{\mu ,\nu }=(1/N)_iS_{i,\mu }^{(a)}S_{i,\nu }^{(b)}`$ is a replica overlap between the $`\mu `$ component of a spin $`i`$ on a replica $`(a)`$, $`S_{i,\mu }^{(a)}`$, and the $`\nu `$ component of the spin on a replica $`(b)`$, $`S_{i,\nu }^{(b)}`$. The bracket $`[\mathrm{}]_\mathrm{c}`$ denotes the configurational average, while $`\mathrm{}`$ denotes the thermal average. The indices, $`\mu `$ and $`\nu `$, stand for three components of spins: $`x,y,`$ and $`z`$. The local chirality $`C_{i,\mu }^{(a)}`$ is defined by three neighboring spins as $`C_{i,\mu }^{(a)}=𝑺_{i+\widehat{𝒆}_\mu }^{(a)}(𝑺_i^{(a)}\times 𝑺_{i\widehat{𝒆}_\mu }^{(a)}),`$ and $`\widehat{𝒆}_\mu `$ denotes a unit lattice vector along the $`\mu `$ axis.

We use the single-spin-flip algorithm with the heat-bath probability.  For a given random bond configuration, eight replicas are prepared with different random initial spin configurations. Each replica is updated in parallel by a different random number sequence, and twenty eight replica-overlaps are calculated and averaged. Note that these data are correlated with each other, and only seven among them are independent. Then, we take averages of $`\chi _{\mathrm{sg}}(t)`$ and $`\chi _{\mathrm{cg}}(t)`$ over the bond distributions. Typical numbers of the bond configurations are tens to thousands depending on the sizes, the temperatures, and the time ranges.

From the non-equilibrium relaxation functions of $`\chi _{\mathrm{sg}}`$ and $`\chi _{\mathrm{cg}}`$, we extract the equilibrium properties and estimate the transition temperature by the finite-time scaling analysis,  which is a direct interpretation of the conventional finite-size scaling through the relation, $`\tau \xi ^z`$. We determine the correlation time $`\tau (T)`$ at each temperature and a combined exponent $`\lambda (\gamma /z\nu )`$ so that the scaled functions, $`\chi (t)t^\lambda `$, for various temperatures fall onto a single curve when plotted against $`t/\tau (T)`$. Then, we extrapolate $`\tau (T)`$ with an assumption that they diverge algebraically when the temperature approaches the transition temperature $`T_\mathrm{c}`$ as $`\tau (T)(TT_\mathrm{c})^{z\nu }.`$ We obtain the most probable estimate for $`T_\mathrm{c}`$ by the least-square fitting of this form.

In the finite-time scaling, it is crucial to use the data free from the finite-size effect. We estimate a characteristic time that the finite-size effect appears by changing the system size from $`L=9`$ to $`L=59`$ at each temperature. The typical example is shown in Fig. 1. The temperature is within our estimate of the spin-glass transition temperature obtained afterwards. The finite-size effect of $`L=19`$ data already appears at $`t=1000`$ steps, after which a long relaxation to the finite-size equilibrium state is waiting. This “finite-size” characteristic time is roughly scaled as $`\tau L^z`$ with the dynamic exponent $`z5`$. By this relation, we confirm that the relaxation of $`L=59`$ before $`t=100000`$ steps is not influenced by the finite-size effect, and thus can be considered as that of the infinite system. At each temperature, we check this time scale and use only the data with the approval.

In this figure, we notice that the appearance of the finite-size effect is in opposite ways between the spin-glass susceptibility and the chiral-glass susceptibility. The former is always underestimated, while the latter is overestimated showing upward bending. This tendency continues to the equilibrium values; the chiral-glass susceptibility always decreases as the system size is enlarged. This size dependence is quite strange as the susceptibility, but it can be understood by considering a definition of the chirality. Since the chirality corresponds to a volume of a parallelepiped spanned by three neighboring spins, it becomes larger when the relative angles between the neighboring spins are enhanced by the finiteness of the lattices. Thus, we must pay much attention to the finite-size effects in the relaxation of the chirality. For example, at paramagnetic temperatures, the chiral-glass susceptibility converges to a finite value in the large size limit, while it shows an upward bending in the small size limit. There exists a crossover size in which the non-equilibrium relaxation looks algebraically divergent as if it exhibits the criticality. This crossover size is considered to correspond to the chiral correlation length. If one tries to detect $`T_{\mathrm{cg}}`$ solely by finding the algebraic relaxation behavior, he may be misled to a wrong temperature where $`\xi =L`$ instead of $`\xi =\mathrm{}`$. Therefore, the chiral-glass transition temperature tends to be overestimated, and the spin-glass transition temperature is to be underestimated, if the system size is insufficient.

As shown in Fig. 1, the initial relaxation of the spin-glass susceptibility continues until $`t10^3`$ steps, then it exhibits the critical relaxation. In the finite-time scaling analysis, we must discard this initial relaxation data. It is also noticed that the finite-size effects of the small sizes ($`L19`$) appear before the initial relaxation ends. The critical phenomena may not be extracted from the equilibrium states of such small systems.

Results of the finite-time scaling of the spin-glass susceptibility and the chiral-glass susceptibility are shown in Fig. 2. The combined exponent $`\lambda `$ and the correlation time $`\tau (T)`$ are estimated by using the data of $`L=59`$ until $`t=70000`$ steps at seven temperatures ranging $`0.24T/J0.32`$. The data in the initial relaxation are discarded. The scaling function converges to a $`t^\lambda `$ line in the equilibrium region ($`t/\tau (T)\mathrm{}`$), where the susceptibility converges to an equilibrium value in the infinite system. On the other hand, it becomes a flat line in the non-equilibrium limit, $`t/\tau (T)0`$, since the susceptibility behaves as $`\chi (t)t^\lambda `$ in the critical region. Therefore, this curve interpolates the non-equilibrium and the equilibrium limits.

The correlation time is extrapolated by $`(TT_\mathrm{c})^{z\nu }`$, and we obtain the transition temperature that gives the best least-square fitting. The exponent is estimated to be $`z\nu =4.7\pm 0.8`$. Combining this value with $`\lambda =\gamma /z\nu =0.36\pm 0.02`$, the exponent for the spin-glass susceptibility is estimated.

$`T_{\mathrm{sg}}/J=0.21_{0.02}^{+0.01},\gamma _{\mathrm{sg}}=1.7\pm 0.3.`$

This transition temperature agrees well with that obtained by the defect free-energy method, and that obtained by the dynamical simulations. The exponent is consistent to that of the non-linear susceptibility in real spin-glass compounds: $`\gamma =2.2\pm 0.1`$ (CuMn & AgMn) and $`\gamma =2.3\pm 0.4`$ (Cd Cr2×0.85In2×0.85S4). This may suggest that the present theoretical model is relevant to the spin-glass materials. Estimates and the comparison of other exponents are left for a future investigation. It requires more detailed and accurate calculations of the Binder ratio, which diverges in the critical region as $`t^{d/z}`$. Preliminary calculations of the Binder ratio, however, suggest that the size effects and the sample dependences are much severer than those of the susceptibility. The algebraic divergence is only observed in lattices larger than $`L=19`$, below which the finite-size effects appear in the initial relaxation. The details will be reported elsewhere.

Analyses of the chiral-glass transition are done in the same procedure. Note that the amplitude of $`\chi _{\mathrm{cg}}`$ is smaller than that of $`\chi _{\mathrm{sg}}`$ about two digits, which causes rather large numerical errors. A possible value of $`\lambda `$ ranges as $`\lambda =0.23\pm 0.04`$. The transition temperature and its exponent are obtained as

$`T_{\mathrm{cg}}/J=0.22_{0.04}^{+0.01},\gamma _{\mathrm{cg}}=1.0\pm 0.5.`$

Note that this value is consistent with the estimate obtained by the equilibrium simulations which gives $`T_{\mathrm{cg}}0.21`$.  The value of $`\gamma _{\mathrm{cg}}`$ does not contradict with that of the model with the Gaussian bond distributions. The exponents $`\lambda `$ for $`\chi _{\mathrm{sg}}`$ and $`\chi _{\mathrm{cg}}`$ estimated in the scaling also coincide with the slope of the relaxation functions at the transition temperature in Fig. 1. This evidence supports the validity of the present analysis. By the results presented in this Letter, it is conjectured that the spin-glass transition and the chiral-glass transition occur at the same temperature. Since the chirality must freeze if the spin freezes, the transition is possibly driven by the spin degrees of freedom.

We have observed the simultaneous finite-temperature spin-glass transition and the chiral-glass transition in the $`\pm J`$ Heisenberg model in three dimensions. This was made possible by observing the non-equilibrium relaxation of the infinite-size system started from the paramagnetic state. We consider the inconsistency of the present results with those of the previous investigations  can be explained by the differences of the lattice size and the temperature range investigated. Generally, one expects to observe the critical behavior as long as the lattice size is comparable with or larger than the correlation length, $`L\xi `$. Contrary, when the correlation length exceeds the size, $`L\xi `$, a very large but a finite cluster swallows the whole system, and only a trivial size effect is observed. As the temperature is lowered from the paramagnetic phase, the correlation length exceeds a lattice size at a temperature above the transition temperature. If one mixes the data above and below this crossover temperature, the finite-size scaling analysis may mislead to a wrong conclusion. Generally, frustration causes a longer correlation length, and the size of the system is sometimes essential to detect the phase transition particularly in frustrated systems. In the present non-equilibrium relaxation analysis, we have used only data of the temperatures and the time ranges where they can be considered as the infinite system. The time range is limited within $`10^5`$ steps, however, it is enough to equilibrate the system at higher temperatures, $`T/J0.3`$, where we start the scaling analysis. We gradually lower the temperature so that the scaled relaxation functions overlap with each other to ride on a single curve. At lower temperatures, $`T/J0.24`$, $`10^5`$ steps correspond to a time scale that the finiteness of the correlation length begins to appear, i.e., the relaxation function begins to bend. In order to see a longer relaxation, we need to enlarge the system size to eliminate the size effect. This must become much more time-consuming and we leave it for the future study. However, it should be noticed that the finite-size effect appears quite early in the relaxation process as shown in Fig. 1, and what is slow is the relaxation to the equilibrium state after this “finite-size” characteristic time scale. The relaxation of the infinite system might be simple: after an initial relaxation of $`10^3`$ steps, the critical relaxation continues to infinity. As for the relation between the model and real spin-glass materials, the critical exponent obtained by the scaling, $`\gamma _{\mathrm{sg}}1.7`$, agrees with the experiments.

The authors would like to thank Professor Fumitaka Matsubara for guiding them to this field. The author T.N. would like to thank Dr. Y. Ozeki, Dr. K. Hukushima, and Professor H. Takayama for fruitful discussions and comments. Use of a random number generator RNDTIK programmed by Professor N. Ito and Professor Y. Kanada is gratefully acknowledged. The computation has been done partly at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo.