text:
# Statistical properties of ultraluminous IRAS galaxies from an HST imaging survey

## 1 Introduction

For more than ten years, much effort has been made to understand the properties of ultraluminous infrared galaxies (hereafter ULIRGs, see Sanders & Mirabel 1996 for a review), and it is now widely accepted that this population of objects represents an important evolutionary stage triggered by strong galaxy interactions/mergers, especially major mergers between gas-rich spirals with mass ratio smaller than 3:1 (e.g., Bendo & Barnes 2000). These ULIRGs may represent an important stage in the formation of QSOs and powerful radio galaxies, as well as an essential step in the formation of elliptical galaxies (e.g., Sanders et al. 1988a; Melnick & Mirabel 1990). Although the nature of ULIRGs is believed to be understood in general terms, it is still not clear whether these objects originate from the interactions/mergers of galaxy pairs or galaxy groups.

On the theoretical side, numerical simulations have successfully reproduced global and even some fine structural properties of elliptical-like galaxies as the results of pair mergers (e.g., Toomre & Toomre 1972; White 1978; Farouki & Shapiro 1982; Negroponte & White 1983; Barnes 1992; Hernquist 1992; Heyl et al. 1994; Barnes & Hernquist 1996; Weil & Hernquist 1996; Walker et al. 1996; Mihos & Hernquist 1994; Springel & White 1998). However, the pair merger model can not account for all the observed properties of elliptical galaxies, which leads people to consider the multiple merger scenario (e.g., Mamon 1987; Barnes 1984, 1985, 1989; Schweizer 1989; Weil & Hernquist 1996). In fact, numerical simulations of multiple mergers do produce remnants different from those produced by pair mergers in both morphological and kinematic properties (Weil & Hernquist 1996). On the other hand, there is increasing observational evidence to support the multiple merger scenario. For example, the co-existence of three OH maser components in the typical ULIRG Arp 220 is of possible multiple merger origin (Taniguchi & Shioya 1998). And Lipari et al. (1999) reported detailed evidence for multiple merger in the luminous infrared galaxy NGC 3256 based on high resolution imaging by HST and ESO-NTT. Furthermore, Borne et al. (1999a, 1999b) carried out a thorough ULIRGs imaging survey by HST WFPC2/NICMOS, of which they displayed about 20 ULIRGs either with multiple nuclei or representing interacting groups (Borne et al. 2000). Based on their results, Borne et al. (2000) suggested that a possible progenitor of the ULIRG population is compact groups of galaxies.

In optical and near infrared images of luminous merging galaxies, such as ULIRGs, there always exist many kpc or sub-kpc scale intense star-forming regions (Lutz 1991; Ashman & Zepf 1992; Holtzman et al. 1992; Zepf & Ashman 1993; Surace & Sanders 1999). It has often been argued that these star clusters may form during violent collisions between gas clouds induced by tidal interactions (Hutchings 1995), and that their formation may be linked directly to intense starbursts occurring at the very centers of ULIRGs (Taniguchi et al. 1998). Compared with numerous star-forming regions in nearby interacting systems or starburst galaxies such as NGC 4038/9 (Whitmore & Schweizer 1995), NGC 3310 and NGC 2415 (Gallagher et al. 2000), star clusters in ULIRGs are generally more luminous and massive, with bolometric luminosity, $`L_{\mathrm{bol}}`$, sometimes ranging up to $`2\times 10^9L_{}`$ and their masses are typically a factor of 100 greater than Galactic globular clusters (Scoville et al. 1999). A likely explanation is that these objects are associations of star clusters seen in more nearby systems (Surace et al. 1998), since the typical intercluster separation in nearby starburst regions is too small ($`18\mathrm{pc}`$, see Meurer 1995) to be resolved at high redshift. Since both such star cluster associations and galactic nuclei may appear as local brightness enhancements on ULIRG images and the luminosities of the largest associations can be comparable with those of galactic nuclei, it is necessary to find a feasible criterion to distinguish these two types of objects of different physical nature. Hereafter for convenience, we adopt a nomenclature in which the luminous star-forming regions alone are named as knots or star cluster associations, while they are named as clumps together with galactic nuclei.

Although there exists observational evidence for multiple merging from investigations on both individual targets and large samples, most of these works are based on qualitative analysis and use somewhat subjective criteria. Therefore, it is important to investigate this issue on a quantitative basis; this becomes more requisite considering the possibility of mistaking star cluster associations as galactic nuclei in ULIRGs, since they are similar in many aspects. This comprises the main subject of the present paper; in our work, we perform a careful photometric analysis on a large HST snapshot imaging survey sample of ULIRGs (Borne et al. 1999a, 1999b), in order to find evidence for multiple merger based primarily on a quantitative criterion of I-band luminosity. On the other hand, there are still no detailed comparisons of multiple mergers with other populations of ULIRGs (e.g., double mergers and single nucleus galaxies, see Sec. 5). This provides the second motivation for our paper. The structure of this paper is as follows: Sec. 2 includes basic information about the ULIRG sample and their observations. In Sec. 3, we describe details of our data reduction. In Sec. 4, we detail our quantitative criterion for identifying putative galactic nuclei. In Sec. 5, we present our main results, i.e., evidence for multiple mergers and a comparison of different ULIRG populations. Finally, we discuss and summarize our results in Sec. 6. Throughout this paper we assume an Einstein-de Sitter cosmology ($`\mathrm{\Omega }_0=1`$) and adopt $`H_0=75\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. At the typical redshift ($`z=0.1`$) of the sample galaxies, $`0.1\mathrm{}`$ (about 1 pixel) corresponds to roughly 200 pc.

## 2 Sample and observations

For our purpose, we use archive data from an HST snapshot imaging survey of ULIRGs (Borne et al. 1999a, 1999b). The sample images were taken using the Wide Field Planetary Camera 2 (WFPC2) in the I-band pass (F814W), including 120 targets all centered in the WF3 chip, of which the resolution is 0.0996″per pixel and the FOV is $`800\times 800`$ pixels. For each target in this survey, two 400s exposures were taken. Most of the observed targets were selected from the QDOT sample (Leech et al. 1994; Lawrence et al. 1999), the rest were selected from the brightest samples of Sanders et al. (1998a, 1998b) and Melnick & Mirabel (1990). Apart from some sources which are saturated, mis-focused or contaminated by bright foreground stars, we make use of 97 independent targets for our study. In addition, all these sources are selected to possess projected nuclear separations less than 20 kpc (for galaxies with more than two putative nuclei, only minimum nuclear separations are considered). This is because generally speaking, strong interactions between galaxies are induced in systems with nuclear separations below this value, i.e., the distance adequate for two interacting galaxies to obtain a substantial disturbance down to their nuclear radii of 1-2 kpc within the lifetime of the encounter (e.g., Lin et al. 1988). Basic information for all the targets is listed in Table 1 in the order of increasing RA, including HST visit No., target name, RA, Dec, redshift, far-infrared luminosity and our classification (see Sec. 5). The redshifts of these objects range from 0.04 to 0.35.

## 3 Data reduction

The sample images were first preprocessed through the standard Space Telescope Science Institute (STScI) pipeline, using a standard WFPC2-specific calibration algorithm and the best available calibration files (Holtzman et al. 1995). The preprocessing consists of mask and analog-to-digital correction, bias and dark subtraction, as well as flat-fielding, etc.. Our post-pipeline calibrations include detection and removal of warm pixels and cosmic ray events. The additional steps of data reduction are discussed below.

### 3.1 Selection for clumpy structures

We use the standard IRAF task imexamine to carefully browse the sources in the ULIRG sample. For each target, three plots are displayed and compared with each other – the contour plot, the surface plot and the snapshot image. The image is adjusted to different greyscale levels until it renders the best visualization. The surface plot is always checked from several different angles to get a full view. A set of ceiling values are tried until the contour plots yield a clear view of the weak structures of the targets. By careful comparison, common local brightness enhancements appearing on each plot are marked as putative clumps. On surface plots, such clump candidates are identified on the basis of detectable brightness peaks; while on contour plots, they are identified on the basis of closed isophotal contours. As for snapshot images, it is hard to resolve them by greyscale adjustment to present all the clumps detectable in surface or contour plots, except for those very bright or compact clumps with high signal-to-noise ratio. However, snapshot images are especially useful in that interacting structures in ULIRGs are more perceptible in these plots, and such structures are very useful for determining whether or not several distinct separate regions near the target position belong to a single dynamically related system. This is done based on the existence of detectable tidal streamers connecting separate regions. Otherwise, we identify the object nearest to the IRAS position as the target we need, and its nearby objects are considered to be either foreground or background ones (see also Sec. 6.1). Hence our method will give a (conservative) lower limit on the fraction of multiple merger objects.

In order to distinguish true physical clumps from Poisson fluctuation peaks in the images, we make rough FWHM (full width at half maximum) measurements on the marginally detectable brightness peaks, and fluctuations are identified to be those peaks with FWHMs smaller than the typical value of observed PSF (point spread function) of HST WF3 chip with F814W filter. Furthermore, it can not be ruled out that foreground stars may coincidently fall onto the target and are thus mistaken as ULIRG clumps. A rough estimation of this probability is made using all the 7 ULIRGs in the sample at relatively low galactic latitude (i.e., $`|b|<15^{}`$). We use the standard IRAF task daofind to search for star-like objects in these 7 ULIRGs which are brighter than the minimum clump flux we have identified in the sample (about 22.5 mag). And by multiplying the star number by the ratio of the target area to the whole field area, we obtain an average contamination number of about 0.1 for one galaxy, hence the influence of star contamination on our results is negligible. Here the area occupied by each target is determined by counting the number of pixels with ADU values at least $`3\sigma `$ above the sky background. This is a conservative upper limit, since the sub-sample we treat is contaminated by foreground stars to the greatest extent due to low galactic latitude.

### 3.2 Surface photometry

We perform surface photometric measurements on the selected clumps to determine their I-band fluxes and luminosities. For well-separated clumps, we apply standard aperture photometry programs in which aperture widths are set to different values according to their actual sizes, ranging from 0.2″ to 1.0″. The size is estimated from the surface brightness profile when it reaches an approximately constant value, i.e., the local background value. In addition, estimates are made of the underlying background galaxy flux by using the mean of the pixels in a 1-pixel annulus immediately outside the photometric aperture. However, some clumps in ULIRGs are so close that their profiles overlap notably, hence standard procedures do not function well. In such cases we perform surface photometry in the following way: we apply a large aperture size to encircle all the overlapping clumps and the magnitude, $`m`$, measured this way is considered to be the combined I-band magnitude of encircled objects; at the same time we apply a small aperture size to measure the magnitude, $`m^{}`$, of the protrudent part (i.e., non-overlapping region) of each encircled clump, and the flux ratio of these parts is approximated as that of the corresponding true clumps, assuming that these overlapping clumps possess similar profiles. Using $`m`$, $`m_1^{}`$ and $`m_2^{}`$, we can estimate the true magnitude of clumps from the following expression,

$$m_1=m+2.5\mathrm{log}\left(1+10^{\frac{m_1^{}m_2^{}}{2.5}}\right),$$
(1)
$$m_2=m+2.5\mathrm{log}\left(1+10^{\frac{m_2^{}m_1^{}}{2.5}}\right).$$
(2)
In all our surface photometry, the photometric calibration is performed using published photometric solutions for WFPC2 and F814W filter (Holtzman et al. 1995).

Simulations are performed in order to estimate the uncertainty of the measured fluxes for overlapping clumps. Using the standard IRAF task mkobjects, we automatically produce two objects (artificial clumps) with different flux ratios and separations. Their FWHMs are fixed to be 5 pixels. For convenience, only de Vaucouleurs brightness profiles are considered. In addition, Poisson noise is added to the simulated objects and the effects of detector gain, readout and PSF are all modeled using the appropriate parameters for the HST WF3 chip. Table 2 lists details of each simulation run. From the table, we can see that the largest photometric uncertainty occurs when the separations between clumps are small and the corresponding magnitude difference of simulated clumps is high. Since there are no overlapping clumps in this sample with very different I-band magnitude, the typical photometric uncertainty introduced by our technique is less than 0.2 mag, and the largest uncertainty in our measurements, corresponding to run 12 in Table 2, is around 0.5 mag. Hence the introduced photometric uncertainty does not affect our results significantly.

## 4 Identification of putative nuclei

As mentioned in Sec. 1, one of the remarkable features of ULIRGs is that there always exist many kpc or sub-kpc scale luminous clumpy structures; and the key issue in our work is to distinguish putative nuclei from star cluster associations. It is found that typical masses for both galactic nuclei in ULIRGs and the giant elliptical cores are several times of $`10^9M_{}`$ (Sakomoto et al. 1999; Lauer 1985). This is typically larger than the mass found for star formation knots. In detailed studies of nearby starburst galaxies such as M82 (de Grijs et al. 1999), no star-forming knots have been observed to possess masses as high as $`10^9M_{}`$. And for the nearest ULIRG, Arp 220, its most massive knot was found to possess a mass of $`6.6\times 10^8M_{}`$ (Shioya et al. 2000). On the theoretical side, complexes with masses greater than $`10^9M_{}`$ may suffer from disintegration (Noguchi 1999). In addition, Taniguchi et al. (1998) suggested that a large mass of $`10^9M_{}`$ corresponds to a maximum star formation efficiency of 1 for knots evolving from superclouds in gravitationally unstable gas disks. Since the real star formation efficiency must be smaller (e.g., 0.1 from observations of the Galaxy), observed knots are always less massive than $`10^9M_{}`$. Given these, we conservatively adopt $`1\times 10^9M_{}`$ as the lower mass limit to pick out putative nuclei from bright clumps in the first step.

To express the above threshold in an observational quantity, the lower mass limit is converted to B-band luminosity using a constant mass-to-light ratio of 6.5 appropriate for spheroids (Fugukita et al. 1998). And the resultant value of $`1.54\times 10^8L_{}`$ corresponds to $`M_\mathrm{B}=15.0`$. Adopting the statistical results of Surace et al. (1998) that the typical $`BI`$ color index for putative nuclei in ULIRGs is 2.0 with a root mean square of 0.7, we obtain the lower $`M_\mathrm{I}`$ limit of $`17.0`$ mag as a quantitative criterion to pick out putative nuclei from bright clumps. This is consistent with the argument that bright knots may have an upper luminosity limit for intense star-formation in a star cluster (Hutchings 1995).

We caution that this quantitative criterion may not be sufficient to pick out ULIRG nuclei in some cases. This is mainly due to two reasons. First, the brightest knots and relatively faint nuclei in ULIRGs may have comparable luminosities. Second, the nuclear and circumnuclear regions of ULIRGs are heavily obscured – as a result, a single nucleus may emerge as multiple condensations due to complex obscuration patterns. In order to eliminate the influences of these two effects, we take further morphological considerations, which are discussed below.

### 4.1 Morphological considerations

Generally speaking, in order to distinguish those extremely bright knots from the putative nuclei, color information is needed (e.g., Surace et al. 1998, Scoville et al. 1999). However, we do not incorporate color index in this paper, because the ULIRGs imaging survey sample we used was performed only in the I-band. In addition, the results of $`BI`$ color show that this index covers a rather wide range from $``$0.3 to 5.2 for star-forming knots (Surace et al. 1998). Considering the relatively narrow span of 0.5 to 3.1 for putative nuclei (Surace et al. 1998), it is obvious that knots can be either bluer or redder than nuclei, as was also pointed out by Hutchings (1995). Therefore, the color index is not an ideal indicator to distinguish putative nuclei from star formation knots.

On the other hand, although the luminosities of some super star clusters are comparable to or even greater than those of putative galactic nuclei, this luminosity overlap is relatively insignificant compared with the overlap of color index. Among all the 57 knots in Surace’s sample, less than 30% are brighter than our I-band luminosity threshold of $``$17.0 mag (see Fig 1 for the magnitude distribution of all these 57 knots). Furthermore, we suspect that some of these very bright knots in Surace’s sample are possible galactic nuclei (see Sec. 6.1).

Star-forming knots with large I-band luminosities (i.e., comparable to those of putative nuclei) may also exist in the current ULIRG sample. To rule them out on a single wavelength basis, we take into account the dynamical information provided by interacting signatures of the targets. This technique led Lipari et al. (1999) to detect three galactic nuclei in NGC 3256 based on HST and ESO-NTT observations. Therefore, after applying the quantitative criterion of $`M_\mathrm{I}`$, we carefully examine the local (such as nuclear arms and disks) as well as global (such as tidal tails, plumes, rings) environments of the sample targets to confirm whether or not clumps brighter than $`M_\mathrm{I}=17.0`$ are actual descendents of progenitor bulges. Both numerical simulations and observations show that each piece of such signatures is connected to at least one nucleus residing at its starting position (e.g., Toomre & Toomre 1972; Wright 1972). Because of this, some very bright clumps (i.e., brighter than $`17.0`$ mag) are still treated as massive star-forming regions if they are located in the outer parts of the sample targets (e.g., at the tips of tail structures). Such objects are more likely to be tidal galaxies produced by merging (Fellhauer & Kroupa 2000), hence they are not included in our work. In our sample, there are 22 such clumps brighter than $`17.0`$ mag in the I-band which are removed due to these morphological considerations. Compared with the total number of putative nuclei of nearly 200, this step, although important, should not affect our statistical results.

This morphological consideration works fairly well, since in most of our cases, such interacting signatures are easy to identify and their geometrical connections to bright clumps in the same images are explicit. Thus this technique eliminates at least some of the uncertainties in our identification of putative nuclei due to the luminosity overlaps between nuclei and knots. However, it should still be notified that some levels of subjectivity cannot be avoided in the morphological considerations mainly due to two reasons. The first is that faint interacting signatures may be buried in the background noise and are not detectable, especially for the targets at relatively high redshift. The second is that the starting positions of these structures are unclear in some cases. e.g., it is not clear how to pick out the putative nuclei from a group of clumps located along a tidal ring, purely from morphological signatures. Since some features of interacting systems may be due to quirks of pre-encounter disks and do not lead to fundamental insights (Barnes 1999), it is not surprising that we sometimes have ambiguous circumstances. In such cases, we simply apply our quantitative magnitude criterion to identify putative nuclei.

### 4.2 Obscuration effects

Since the high far-infrared luminosity of ULIRGs results from thermal re-radiation from dust (Sanders & Mirabel 1996), there is no doubt that ULIRG morphologies tend to be affected by dust obscuration. Images taken at longer wavelength are needed to examine the extent to which the morphologies of the sample targets are influenced. For the near infrared images of 9 LIGs plus 15 ULIRGs taken by HST NICMOS camera (at 1.1, 1.6 and 2.2 $`\mu `$m), it was stated that although new interacting features as well as new super star clusters could be found at longer wavelengths, no new nuclei were seen in these images (Scoville et al. 1999). This fact suggests that global morphologies are not greatly affected by dust obscuration effects at optical wavelengths, and hence our identification of putative nuclei should be reasonably secure.

In fact, to further minimize the above-mentioned effect, some targets in the sample were treated separately, as in the cases of IR09427$`+`$1929 , IR06361$``$6217 and IR04024$``$8303 , etc. Although these targets show signs of separate nucleus components, the combined shapes of their brightness profiles and loss of any nearby signatures of on-going interactions suggest that they are more likely to be a single nucleus split by a foreground dust lane. Surace et al. (1998) have reported such cases in their warm ULIRG sample, e.g., IR05189$``$2524. This procedure is also consistent with Borne et al.’s criterion that only cases with clearly separated optically luminous galactic components were selected (Borne et al. 2000).

To summarize, the identification of putative galactic nuclei in our reduction includes two steps. The first step is to apply directly an I-band luminosity criterion, in which all the ULIRG clumps brighter than $``$17.0 mag were considered as nucleus candidates for further consideration. The second step is to include morphological information (primarily, the existence and characteristics of tidal features), in order to remove some super star clusters with I-band luminosities comparable to those of putative nuclei. Only nuclei which pass through these two criteria are considered as putative galactic nuclei. In the next section, we present statistical analyses of this group of objects.

## 5 Results

Based on the quantitative criterion discussed in Sec. 4 as well as morphological structures, we pick out all the putative galactic nuclei from clumpy structures in each sample galaxy. And all the 97 ULIRGs can be classified as three categories according to the number of their nuclei:

* Multiple nuclei ULIRGs (multiple mergers)
* Double nuclei ULIRGs (double/pair mergers)
* Single nucleus ULIRGs (single remnants)

Note that for the double nuclei category, two galaxies with projected nuclear separations of more than 20 kpc are excluded (IR13156$`+`$0435 and IR15168$`+`$0045). Hence these three categories compose a nearly complete sample of strong interacting systems and their remnants with ultraluminous far-infrared luminosities. Note that here we adopt a working assumption that every nucleus present in merging systems corresponds to a single progenitor nucleus (bulge), i.e., we assume that dynamical processes (such as tidal disruption) cannot produce compact regions with masses beyond $`10^9M_{}`$. We caution that this assumption, while reasonable, needs to be verified using numerical simulations.

Photometric results of all the putative nuclei in the sample are listed in Tables 3, 4 and 5 for the three categories, respectively, including positions and I-band absolute magnitudes. In all these tables, nucleus positions are expressed by RA and Dec relative to the corresponding target positions given in Table 1, which are fixed to pixel coordinate of (420, 424) in the HST WF3 chip. For multiple and double mergers (see Tables 3 and 4), projected nuclear separations as well as luminosity ratios are included, and we only give minimum separations for the multiple nuclei category. Furthermore, some basic information about the targets are also listed in these tables, including tail lengths, brief comments on morphology and available spectral types.

Our results show that among all the 97 ULIRGs of the treated sample, 17 are multiple nuclei systems, 38 are double nuclei systems, while the remaining 42 have only one identifiable galactic nucleus. This gives plausible fractions of ULIRGs with multiple, double and single nucleus/nuclei as 18%, 39% and 43%, respectively. These results evidently support the multiple merger scenario as a possible origin of ULIRGs, in addition to the widely-accepted pair merger picture. This argument was put forward by Borne et al. (2000) based on the same sample. However, our results are obtained in a more objective way due to our mainly quantitative criterion. Comparisons between Borne’s results and ours are presented in Sec. 6.2. In addition, further statistical studies on this large sample reveal some differences among the three ULIRG categories in our taxonomy, some of which give hints to a possible evolutionary sequence from multiple mergers to double mergers and then to single merger remnants (hereafter M$``$D$``$S sequence, see Sec. 5.2). These results not only support our classification scheme, but also give some insights into the dynamics of galaxy interactions/mergers (see Sec. 6.2).

Besides putative nuclei, there often exist several bright knots in most of the sample galaxies. One of the significant characteristics of these bright knots is that they tend to distribute around nuclear regions, or in regions between/among separate nuclei. This is consistent with the argument that star-forming concentrations are preferentially situated along the overlapping areas of interacting galactic disks (Scoville et al. 1999). Besides these nuclear knots, there are also cases with circum-nuclear knots, or even with knots residing at the tips of tidal tails. As pointed out by Weilbacher et al. (2000), a possible descendent of such objects at the tips of tidal tails is tidal dwarf galaxies. Another tendency is that bright knots are less frequently seen in single nucleus ULIRGs, compared with multiple/double mergers. We plan to study the properties of these bright knots and investigate their formation mechanisms and evolutionary fates in a future work.

### 5.1 Multiple nuclei ULIRGs

Table 3 gives the information on each nucleus for all the 17 multiple mergers, while the HST I-band images, surface plots and contour plots are shown in Figs. 1-17 (arranged from left to right, respectively). We put plus marks on each image to indicate the positions of putative nuclei, and the orientation of the image is indicated by the headed north arrow and the unheaded east arrow. Since pixel spacings are not always the same for the two axes, the north and east directions may sometimes be oblique. The orientation of each contour plot is selected to be the same as the snapshot image for direct comparison. For the three-dimensional surface plot, we rotate it along both horizontal and vertical viewing directions until it renders the best sight to view each nucleus clearly. In addition, for each multiple merger candidate, the snapshot image and the contour plot both contain distance scale rulers, and the scale of the surface plot is selected to be the same as that of the contour plot. For the 17 multiple merger systems, the morphologies, the relative positions and numbers of putative nuclei as well as knots are very different. We give a description of each target in turn.

This galaxy is mainly composed of two parts separated by 4.3 kpc. The northeastern part contains a single nucleus, while the southwestern part consists of two close nuclei with a projected separation of only 570 pc. Bright knots of the system preferentially distribute to the southwest of each nucleus. There is a 45 kpc tidal tail starting from the northeastern nucleus and curving to the south. In addition, the contour plot of IR00161$``$0850 presents distinct arm structures connected to both parts of the system. These structures may be remnants of progenitor spiral arms which have not been completely destroyed by the merging process.

The western part of IR02364$``$4751 contains the brightest nucleus of the system, while the eastern part is composed of two close nuclei separated by 540 pc. The projected distance between the two parts is about 2.6 kpc and several bright knots can clearly be seen between them. A 50 kpc tidal tail starts at the eastern nucleus pair and curves from the east to the west. A faint plume structure can also be seen starting from the western nucleus and stretching to the southwest. A possible evolutionary scenario is that two spirals merge first and produce the long tail we see; before they finally evolve into a single object, a third galaxy approaches from their western side, which makes up a system with three putative nuclei.

The morphology of this galaxy is like that of Mrk 273. A very long tidal tail of nearly 40 kpc starts from the central star-like object and stretches from the northeast to the southwest. Surface and contour plots of the target clearly show that the central object is in fact composed of three closely separated nuclei all brighter than $`20.5`$ mag in the I-band. According to the spectral information, the central active galactic nuclei (AGN) phenomenon is triggered in this target. The three progenitor galaxies of this system are likely to have begun merging almost simultaneously.

The morphology of this galaxy is asymmetric and complex. Three bright putative nuclei of about $`19.0`$ mag exist in the system, with a roughly linear distribution. The projected separation is about 800 pc (between the center and western nucleus) and 1.1 kpc (between the center and eastern nucleus). Two plume structures are detectable in the system, with the western one being 9 kpc long and the northeastern one 5 kpc long. Several bright knots distribute just along these two plumes. Furthermore, two much fainter plumes can also be seen from the snapshot image, stretching to the southwest and east, respectively. These fainter plumes are also much longer, the southwestern one is nearly 28 kpc long while the eastern one is about 23 kpc long.

From the snapshot image, this galaxy consists of two parts located to the northwest and southeast, respectively. The projected separation between these two parts is about 6.5 kpc. The southeastern component embeds two close nuclei separated by 2.2 kpc, and a tidal plume starts here stretching as far as 10 kpc. This gives further support to the merger origin of the southeastern component.

This galaxy is a system with peculiar morphology. The three putative nuclei we identified approximately distribute along a line. The central nucleus is 0.5 kpc and 1 kpc away from the eastern and western nucleus, respectively. In the global environment, the three closely separated nuclei reside at the southern edge of a tidal ring with 11 kpc diameter, and several bright knots are scattered along the ring perimeter. These putative nuclei are also connected with a 14 kpc tail curving from the east to the west. A possible explanation of this peculiar morphology is that two galaxies interact to produce a tidal tail, and before their bulges finally merge together, a third galaxy impact on them perpendicular to their disk planes, which results in the ring structure we see now.

This system contains two main parts separated by about 15 kpc, and a material stream can be detected connecting these two parts. The center of the western part is a compact nucleus with I-band absolute magnitude about $`19.3\mathrm{mag}`$. Two distinct spiral arms extend from it and stretch to the north and south, respectively. Along the two arms there are several bright and compact knots, two of which are almost as bright as $`18.0`$ mag. However, since these very bright clumps reside at the tip of the southern spiral arm, they are more likely to be progenitors of tidal dwarf galaxies, and hence are not identified as putative nuclei in our work (see Sec. 4). In addition, two other nuclei brighter than $`19.0`$ mag are located along the two spiral arms.

Surface and contour plots of IR14575$`+`$3256 show clearly that there are three distinct putative nuclei in the system. A close pair separated by 1 kpc is located to the east, while a fainter nucleus is located to the west. The projected separation between these two parts is about 2.4 kpc. Although its three nuclei are closely distributed, IR14575$`+`$3256 only presents moderate interacting features. A 10 kpc arm starts from the two eastern nuclei and curves to the southwest, while a 6 kpc arm starts from the western nucleus and curves to the north. Several bright knots preferentially distribute along these two arms as well as around the two eastern nuclei.

This is an object with very peculiar morphology which may be hard to explain as a pair merger product. The snapshot image shows a round center component, which is connected with two other bar-like structures to the northwest and southwest, respectively. Each of these three components encompasses a bright nucleus. In fact, there are two closely separated clumps in the southwestern bar which are both brighter than $`18.0`$ mag. However, due to their contour shapes, these two clumps appear to be a single nucleus split by a foreground dust lane (as in the case of IR05189$``$2524, see Surace et al. 1998). Hence we conservatively take them as a single object, and the $`M_\mathrm{I}`$ value listed in Table 3 is the combined magnitude of these two objects. A 11 kpc tidal plume starts from the central component and curves to the south. In addition, another tidal tail as long as 19 kpc is also detectable connected to the northwestern bar and stretches from the west to the east.

The morphology of this galaxy is rather complicated. A strong tidal ring can be seen from the HST snapshot image and its diameter is about 20 kpc. The main part of this galaxy is located at the northern edge of the ring perimeter. The galaxy is clearly composed of four compact nuclei, as shown in the surface and contour plots. These four putative nuclei are very bright, each with I-band absolute magnitude around $`20.0\mathrm{mag}`$. 6 kpc away to the north of these four nuclei is located the fifth nucleus of the system, and a 12 kpc plume structure just starts here stretching from the south to the north.

This is an example of a QSO hosted in an interacting galaxy system. Besides the central brightest nucleus with I-band absolute magnitude of nearly $`23.0`$, there exist another two bright nuclei which are separated by 1.1 kpc and located to the southeast. A distinct ring structure with a diameter of about 19 kpc locates to the north of the system. This galaxy resembles 3C 48 and IR04505$``$2958 in morphology and it provides us with direct evidence that the AGN phenomenon can be triggered by galaxy merging.

From the HST morphology, this is clearly a merging pair with two parts to the north and south, respectively. These two parts are separated by 6 kpc and some bright knots can be seen located around each part. Both the northern and southern parts consist of two distinct nuclei, separated by 1.6 kpc and 3.1 kpc, respectively. In addition, tidal plume structures can be seen connected to these two parts. From its morphological features, it is likely that the progenitor of this system is a group of four galaxies; each pair of galaxies merge first, while the whole system begins its encounter before each pair finishes its own merging process.

From the HST snapshot image, this galaxy is composed of two main parts to the south and north, respectively. The projected separation between these two parts is about 5 kpc. The southern part encompasses a single nucleus, while the northern part has four close nuclei in a roughly linear distribution, whose high I-band luminosities suggest that they are not very likely to be star-forming knots. Only weak interacting features are detectable in the system. A faint plume stretching from the south to the north is connected with the four northern nuclei. In addition, another faint tidal ring is detectable to the east of the main body.

This ULIRG has a rather peculiar morphology which consists of two parts from the snapshot image. The northeastern part has a ring-like structure with a diameter of more than 8 kpc. Four separate nuclei are concentrated at its south side, while a number of bright knots can also be clearly seen along the ring. The southwestern part of this system consists of two close nuclei separated by 600 pc and their connecting line is perpendicular to the northeastern ring plane. In addition, there exists a tidal bridge connecting these two parts.

This is a three nuclei ULIRG with very peculiar morphology. The bright object located in the center is composed of two close nuclei separated by 1.2 kpc, while another much fainter nucleus resides to their northeast. A clear tidal bridge can be seen connecting these two parts, and several bright knots are distributed just along the bridge. Two plume structures of length 6 kpc and 18 kpc start from the central nucleus pair and curve in opposite directions. In addition, connected to the northeastern nucleus, there is a 25 kpc tidal tail stretching directly from the east to the west.

This system contains three putative nuclei with very different luminosities. The one located at the center has I-band luminosity as high as $`2.6\times 10^9L_{}`$. The nucleus located 1.2 kpc to the southeast of the central bright nucleus and the one 2.4 kpc to the west are much fainter. In addition, bright star-forming knots of the system are preferentially distributed around the western nucleus. A multiple merger scenario of this system is also supported by its complex morphology with a 11 kpc tidal tail and two other plume structures.

From the snapshot image, this is a galaxy with a relatively simple and regular morphology. A distinct tidal tail about 19 kpc long starts from the central bright object and stretches from the northwest to the southeast. The central object is in fact composed of three close nuclei all brighter than $`19.5`$ mag, and around them there are several bright star-forming knots of the system. The three progenitor galaxies of this system are likely to have begun to merge almost simultaneously.

As described above, all of the targets listed in Table 3 encompass at least three putative nuclei and most of them clearly show signs of strong interaction activities. Such large scale tidal features coined as tidal traumas (Barnes & Hernquist 1992) include tails, bridges, rings, plumes, etc. In some systems with very peculiar morphologies, different types of tidal features may co-exist, such as several tails or plumes (e.g., IR22206$``$2715 in Fig. 15). Borne et al. (2000) suggested that merging pairs of spirals alone can not account for all the dynamical diversities of ULIRGs, and numerical simulations also support that such complex structures tend to appear in multiple merger systems (Barnes & Hernquist 1996).

### 5.2 I-band magnitudes of putative nuclei and spectral types

The I-band magnitude distributions of putative nuclei for the whole sample and the three different categories are presented in Figs. 18a-d, respectively. The faint side cutoff in Figs. 18b and 18c (multiple and double nuclei galaxies) is set by our criterion. In Fig. 18d (single nucleus galaxies), the faint side cutoff is at $`18.0`$ mag. This might be due to the merger-induced formation of massive and bright bulges in these systems. On the other hand, there exists a distinct high luminosity tail above $`M_\mathrm{I}=21.0`$ in Fig. 18d. This corresponds to a Seyfert 1/QSO excess in single nucleus ULIRGs, since most nuclei brighter than $`21`$ mag (9 out of 13 with available spectra) belong to galaxies with Seyfert 1/QSO spectra (see Tables 3, 4 and 5). In addition, Seyfert 1/QSO nuclei can occasionally be found in multiple/double mergers (e.g., IR20037$``$1547 and IR00204$`+`$1029 ), and there is either likelihood that these bright nuclei are primordial objects before merger or in fact evolutionary descendents of previous merging steps.

The magnitude distributions of multiple/double mergers and single nucleus galaxies give a median I-band magnitude of $`19.3`$ ($`2.3\times 10^9L_{}`$) and $`20.1`$ ($`4.7\times 10^9L_{}`$), respectively. The high median luminosity for the single nucleus category is partly due to the merger-induced formation of massive, bright bulges, and partly due to the triggering of active galactic nuclei (AGN) phenomenon in a substantial fraction of merger remnants, as mentioned above. In addition, these three categories do not show much difference in the $`M_\mathrm{I}`$ distribution width. Although magnitude distributions of multiple and double mergers do not present much differences, it is obvious from the comparison of Figs. 18b and 18c that bright nuclei (i.e., from $`19.5`$ to $`20.5\mathrm{mag}`$) tend to appear more often in double nuclei systems (nearly twice as often than in multiple mergers). This gives some hints to a possible evolutionary sequence from galaxy groups to merger remnants with a single nucleus, along which the double nuclei systems remain as intermediate products (M$``$D$``$S sequence).

To check this M$``$D$``$S merger sequence, we calculate the luminosity ratio of each nucleus to the brightest one in the same target. Our results show that nearly 40% of the double merger ULIRGs contain two nuclei with mass ratio more than 3:1 (assuming a constant mass-to-light ratio for bulge). Although spiral galaxies might possess rather different bulge-to-disk mass ratios, it cannot be ruled out that a substantial fraction of putative nuclei in pair merger are likely to be merger remnants of previous galaxy encounters, since ULIRGs are postulated to arise from mergers of two comparably massive gas-rich spiral galaxies (e.g., Mihos 1999). This gives further support to a possible M$``$D$``$S evolutionary sequence.

Further investigation based on the optical spectra available in the literature (Lawrence et al. 1999; Veilleux et al. 1999) as well as the information from NED <sup>1</sup><sup>1</sup>1The NASA/IPAC Extragalactic Database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. and our observations shows that about 64% (7/11) Seyfert 1/QSO in this HST snapshot imaging survey sample belong to galaxies with a single nucleus. This gives an important hint that the Seyfert 1/QSO phenomenon tend to appear at the final stage of galaxy interactions/mergers when separate nuclei have merged together. Table 6 gives the detailed statistics of different spectral types for the three categories. It shows that the proportions of H II region-like spectra are 64% (9/14), 35% (8/21) and 20% (5/25) for multiple, double and single nucleus/nuclei galaxies, respectively. While the AGN proportions for these three categories are 36% (5/14), 65% (13/21) and 80% (20/25). This clear trend of gradually changing proportions of different spectral types suggests that the M$``$D$``$S sequence is also one with changing energetics: from starburst-dominated to central AGN-dominated. Of course, it should be mentioned that the H II-like spectra in some objects may be due to the super star clusters rather than the identified putative nuclei, which is consistent with the fact that the images of these targets are distributed by knotty structures in their circumnuclear regions (e.g., IR14202$`+`$2615 and IR18580$`+`$6527 ).

### 5.3 Nuclear separations and tail lengths

Tables 3 and 4 list the projected nuclear separations for multiple and double nuclei ULIRGs. We present the corresponding distributions in Fig. 19a. Given the HST WFPC2 resolution limit, our measurements cannot probe separations below several hundred parsecs, except for some low redshift galaxies. At the maximum redshift in the sample, $`z=0.35`$, our measurements can only probe nuclear separations larger than 1.5 kpc. Therefore the distribution of nuclear separations less than this value is incomplete. This incompleteness is particularly serious at separations below 0.5 kpc, since very few targets in this sample are at low redshift (i.e., less than 20% with redshift smaller than 0.1). This is shown by the sharp decrease below 500 pc in Fig. 19a. We also give the same distribution based on results from the literature (Surace et al 1999; Evans 1999; Rigopoulou et al. 1999; Murphy et al. 1996; Clements et al. 1996) in Fig. 19(b) for comparison. One significant feature of our results is that there is a high fraction (3/4) of nucleus pairs with projected separations less than 5 kpc, which is much higher than the result from the literature (1/3). The reason for this large difference is mostly due to the low resolution limit of ground-based observations in previous works. This result supports that systems at late merging stage preferentially host ULIRG galaxies (Mihos 1999).

Fig. 20 displays the correlation between nuclear separations and tail lengths. Most of the data points are distributed in the lower left conner of the plot. However, considering that both the minimum separations and tail lengths applied here are projected values, the distribution profile of data points is more instructive. Fig. 20 shows a weak trend of smaller nuclear separations coupled to longer tidal tails. Simple qualitative analysis suggests a picture that when galaxies first encounter at a relatively large distance, their mutual tidal field is not strong enough to produce long tail structures; and when they are drawn closer as the merging process proceeds, disk stars in each galaxy are catapulted farther out and tidal tails continue to lengthen. This picture is consistent with numerical simulations (e.g., Toomre & Toomre 1972; Wright 1972). However, there do exist several examples with both long tails and large nuclear separations (e.g., IR21547$``$5823 ). This might imply that at least one nucleus in such systems is in fact an evolutionary descendent of a previous merging process, or merely that the resolution limit of HST prevents us from probing fine structures at its very center. Considering this, the multiple merger fractions we determined at the beginning of Sec. 5 might be an under-estimate of the true fraction. In addition, relatively short tail structures can occasionally be found in single nucleus ULIRGs. This further strengthens the case that these systems are merger remnants of previous galaxy encounters. Of course, our measurements of tail length carry an obvious bias: the observable tidal tails at a given surface brightness in distant galaxies tend to be shorter than their nearby counterparts even if they have the same physical lengths, due to the $`(1+z)^4`$ surface brightness dimming. However, since the sample targets do not cover a wide range of redshifts (0.04$``$0.35), this effect should not influence our statistical results significantly.

## 6 Discussion and summary

### 6.1 Multiple merger fraction

Based on our quantitative criterion as well as morphological features, we obtain the fractions for multiple, double and single nucleus/nuclei ULIRGs as 18%, 39% and 43%, respectively. A question is why there is no report of multiple merger cases in other sample studies on ULIRGs, e.g., in the sample of Surace et al. (1998). To test the consistency between these two results, we directly apply an upper mass limit of $`10^9M_{}`$ on their results of ULIRG knots as well as putative nuclei, and find that two targets in their sample (IR12071$``$0444 and IR15206+3342) should be classified as multiple mergers, namely, there are at least three clumps too massive to be explained as star cluster associations in each of these two targets. Considering this, the fraction of multiple nuclei ULIRGs in Surace’s sample is 22%, very close to our result of 18%. The main reason why Surace et al. (1998) did not choose these objects as multiple merger systems might be that only the one or two brightest clumps in their sample were believed to plausibly represent putative nuclei. This means the pair merger origin was considered to be the only evolutionary picture for ULIRGs. This assumption may sometimes give rise to some inconsistencies or ambiguities. This is especially true in the case of IR12071$``$0444. For this system, there are three bright objects that have comparable masses (several times $`10^9M_{}`$), but only one of them was selected as a putative nucleus. Furthermore, this identified nucleus also possesses $`M_\mathrm{B}`$, $`M_\mathrm{I}`$ and $`BI`$ values just between those of the other two bright clumps (see Table 2 of Surace et al. 1998). This means that at least this putative nucleus does not show distinguished features compared with some other clumps in the same target. Thus it is possible that any of the three brightest clumps in IR12071$``$0444 represent putative nuclei. On the other hand, in our work, the putative nuclei were identified independent of the evolutionary origin for ULIRGs, in order to check the possibility of the multiple merger scenario proposed by some researchers.

Based on the same sample, the multiple merger candidates selected by Borne et al. (2000) include 9 ULIRGs classified as multiple nuclei galaxies and 8 ULIRGs classified as interacting groups (three multiple mergers in Borne’s classification are absent from our subsample, and are therefore not discussed here). The HST snapshot imaging survey sample does include a fraction of ULIRGs which are very likely to be located within galaxy groups from their images, since there are several bright or dwarf galaxies surrounding them. However, we did not consider them as multiple merger systems except when the central target contains at least three putative nuclei or is connected with some of its companions by obvious interacting signatures. This is because further redshift observations need to be performed for confirmation. These group candidates include IR00335$``$2732 , IR01031$``$2255 , IR06268$`+`$3509 , IR13342$`+`$3932 and IR22546$``$2637 . In fact, we have carried out some spectroscopic observations using the 2.16 m telescope at the Xinglong Station of Beijing Astronomical Observatory and identified several luminous/ultraluminous IRAS galaxies located in galaxy groups (Zou et al. 1995; Wu et al. 1998). However, for two multi-merger (group) candidates in Table 1 (IR02459$``$0233 and IR12202$`+`$1646 ), neither of them belong to an interacting group since the surrounding galaxies are either foreground or background objects (Zheng et al. 1999). Considering that the evolutionary sequence from galaxy groups to ULIRGs and then to ellipticals is attractive from both theoretical and observational points of view, further detailed multi-wavelength observations on these group candidates have been scheduled.

As for the remaining 12 multiple mergers in Borne’s classifications, six of them are consistent with our classifications (IR13539$`+`$2920 , IR16007$`+`$3743 , IR18580$`+`$6527 , IR20100$``$4156 , IR20253$``$3757 and IR23515$``$2421 ). We carefully examined the six discrepant targets, of which five of them contain only one or two clumps brighter than $`17.0\mathrm{mag}`$ in the I-band. Therefore they are classified as single or double nucleus/nuclei ULIRGs. The case for IR14337$``$4134 is somewhat confusing. This target seems to reside in a group composed of several bright galaxies as well as star-like objects, of which the biggest galaxy indeed contains three obvious nuclei. However, the star-like object located in the center is the nearest target to the IRAS position. In addition, the QDOT redshift survey (Lawrence et al. 1999) reveals that its optical spectral type is Seyfert 1 and the calculation of far-infrared luminosity for IR14337$``$4134 is also based on the redshift of this Seyfert galaxy. Considering this, Borne et al. (2000) may have identified the wrong target on the WFPC2 snapshot image. Putting aside IR14337$``$4134 , we suggest that the remaining five multiple mergers in Borne’s identification which are absent from our results are due to different selection criteria involved. Borne et al. (2000) identified several targets with very complex tidal features to be multiple merger cases with no regard for the exact numbers of their observable nuclei. Although it is likely that a very complex morphology does point to a possible multiple merger origin, we identify multiple mergers simply based on the exact number of clumps brighter than $`17.0`$ mag, because any tight correlation between morphologies and nucleus numbers needs to be verified by further numerical simulations of merging processes.

### 6.2 Some hints for merging dynamics

From luminosity distributions in Sec. 5.2, it is obvious that putative nuclei in single nucleus ULIRGs tend to be more luminous than those in multiple/double mergers. Statistical studies on available spectral information also reveal that there is a Seyfert 1/QSO excess in the single nucleus category. Therefore, our results strongly support the argument of Kauffmann & Haehnelt (2000) that bulge and supermassive black holes may both grow in galaxy merging. And our results also give some hints to a possible close relation between central black holes and the bulges of their host galaxies. (Gebhardt et al. 2000).

Although only a small fraction of snapshot images are given in this paper, there is no doubt from morphological comments in Tables 3, 4 and 5 that most ULIRGs with single nuclei have weak interacting signatures (in many cases, only weak plume structures are detectable). On the other hand, peculiar morphologies frequently emerge in double/multiple nuclei ULIRGs, such as long tidal tails and distinct ring structures. This indicates that strong interacting features always thin out at final merging phase, which is consistent with results from numerical simulations (e.g., Barnes & Hernquist 1996).

From the detailed descriptions in Sec. 5.1, we can see that multiple merger processes may be very complicated and the morphologies of their intermediate products may be very diverse. This indicates that the evolutionary history of interactions/mergers in galaxy groups must vary due to different initial conditions. When a group of galaxies begins to merge, the central massive spiral may swallow its satellite galaxies one by one; on the other hand, merging processes may also happen between several galaxy pairs or among sub-groups simultaneously, then these separate parts start a second merging step. Major mergers between comparably massive spirals alone cannot account for the dynamical diversity of the ULIRG population (Borne et al. 2000).

We would like to express our gratitude to the staff members at the Beijing Astronomical Center, especially W.-P. Lin and H. Wu, for generous help in data reduction. We also wish to thank Z. Zheng and the anonymous referee for useful comments and suggestions, which improved the paper. We are also grateful to Martin Smith for a careful reading of the text. This research was based on observations obtained with the NASA/ESA Hubble Space Telescope through the Space Telescope Science Institute, Astronomy, Inc., under NASA contract NAS5-26555. This project was supported by the NSFC and NKBRSF G19990754.

Figure captions


format: text

text:
# Differentiation and Replication of Spots in a Reaction Diffusion System with Many Chemicals

## Abstract

The replication and differentiation of spots in reaction diffusion equations are studied by extending the Gray-Scott model with self-replicating spots to include many degrees of freedom needed to model systems with many chemicals. By examining many possible reaction networks, the behavior of this model is categorized into three types: replication of homogeneous fixed spots, replication of oscillatory spots, and differentiation from ‘multipotent spots’. These multipotent spots either replicate or differentiate into other types of spots with different fixed-point dynamics, and as a result, an inhomogeneous pattern of spots is formed. This differentiation process of spots is analyzed in terms of the loss of chemical diversity and decrease of the local Kolmogorov-Sinai entropy. The relevance of the results to developmental cell biology and stem cells is also discussed.

A simple chemical reaction-diffusion system, the Gray-Scott (GS) model, has recently been studied extensively. This model exhibits self-replicating spots within some range of parameter values. Similar replicating spots have been found in experiment.

In each model system studied to this time, only a single type of spots appear, and there is no differentiation to create different types of spots. To search for possibility in complex pattern formation as in biological systems, it is interesting to study a reaction-diffusion system displaying not only the division but also differentiation of such cell-like spots. There it is more necessary to consider chemical reaction dynamics within each spot which are more complex than those adopted in the GS model. Also, in order for spots to continue replication and differentiation, their division process and their internal complex chemical dynamics must somehow organize a proper relationship. We consider a model with such sufficiently complex dynamics and study how this relationship is organized. In so doing, we find the connection between self-replicating spots and cell differentiation in dynamical systems models assuming the cell itself. We do this by adopting a many-degree-of-freedom version of GS model, with the goal of realizing spot differentiation.

As a variant of the GS model we study the following model, possessing the many degrees of freedom needed to model a complex catalytic reaction network:

$`{\displaystyle \frac{u_i(x,t)}{t}}=D_u^2u_i(x,t)+A(1u_i(x,t))`$
$`u_i(x,t){\displaystyle \underset{j=1}{\overset{Q}{}}}{\displaystyle \underset{k=1}{\overset{Q}{}}}W_i^{(j,k)}v_j(x,t)v_k(x,t)`$
$`{\displaystyle \frac{v_i(x,t)}{t}}=D_v^2v_i(x,t)Bv_i(x,t)`$
$`+v_i(x,t){\displaystyle \underset{j=1}{\overset{P}{}}}{\displaystyle \underset{k=1}{\overset{Q}{}}}W_j^{(i,k)}u_j(x,t)v_k(x,t).`$
Here $`u_i(x,t)`$ denotes the concentration of the $`i`$th inhibitor and $`v_i(x,t)`$ that of the activator chemical. Each inhibitor is produced at a constant rate, $`A`$, and both the activators and inhibitors decay at rates proportional to their concentrations (with different coefficients). The diffusion constants of inhibitors and activators are denoted by $`D_u`$ and $`D_v`$, respectively. We assume that these constants are equal for all chemical species for simplicity. This model is a variant of that proposed by Cronhjort and Blomberg , and is reduced to the GS model by setting $`P=1`$, $`Q=1`$, $`W_0^{(0,0)}=1`$. Here we fix $`P=3`$ and $`Q=20`$, and choose the reaction matrix $`W`$ randomly, under the constraint that the reaction for each activator’s replication is catalized by $`k`$ (0$``$$`k`$$``$$`K`$) randomly chosen activators. (Here we have explored the cases in which $`K=4,5,6`$.) The specific numbers of chemicals $`P`$ and $`Q`$ are not important, but we note that it is very difficult to realize the differentiation we study, when these numbers are small. While here we present results for the case of one spatial dimension, numerical results for the two-dimensional case indicate similar behavior.

The behavior of our model depends on the choice of the reaction network $`W`$ (which may be regarded to represent a prototype of a complex intracellular autocatalytic reaction network). We carried out numerical computations using a variety of randomly chosen reaction networks $`W`$, and with these we classified the possible types of dynamics of our system. Although we have studied a large number of networks, the observed behaviors belong to each of the four classes as will be shown.

For each set of simulations, this matrix was fixed, while the parameter values were set throughout as $`D_u=1.0`$, $`D_v=0.010`$, $`A=0.020`$, and $`B=0.060`$ or $`0.070`$. (Below we discuss bifurcation with the change of the parameters $`A`$ and $`B`$.) These parameters were chosen so that the spot structure displayed by the GS equation is displayed by the presently considered model also. We used $`\mathrm{\Delta }t=0.010`$ and $`\mathrm{\Delta }x=1.0`$ for the numerical integration, while we have confirmed that the numerical results are unchanged by using smaller values for $`\mathrm{\Delta }t`$.

In general we begin from some initial conditions from which two simple spots are formed. These spots then produce additional spots, and the process continues until eventually spots are distributed throughout the entire system. When such a state is realized, spot division ceases. Here, we first classify the (transient) dynamics of our model starting from such initial conditions. These results were obtained by considering several thousand randomly chosen reaction matrices.

To characterize the spot dynamics quantitatively, we have introduced two quantities measuring the diversity of these dynamics. One is the chemical diversity $`S_i(nT)`$, which is defined by
$`S_i(nT)=_{j=1}^QP_i^j(nT)log(P_i^j(nT))`$,
with $`P_i^j(nT)=\overline{v_i^j}(nT)/_{j=1}^Q\overline{v_i^j}(nT)`$,
and $`\overline{v_i^j}(nT)=(1/T)_{t=(n1)T}^{nT}v_i^j(t)`$, ($`n=`$1,2,..),
where $`v_i^j`$ is the $`j`$th chemical concentration at the center of the $`i`$th spot. The interval $`T`$ used for the average is chosen to be on the order of a time scale for spot division.

The second quantity is the ‘local KS-entropy’ $`h_i`$ of the $`i`$th spot. The local KS entropy here is defined as the sum of positive local Lyapunov exponents , obtained by using the tangent vectors corresponding to the chemical concentrations of the spot in question .

The spot dynamics are classified into the following four types.
(i) fixed-point case
In this case, the set of concentrations within every spot converge to the same fixed point. Although the spatial pattern itself is not completely homogeneous, each spot has identical, fixed chemical concentrations, and they are separated by equal distance.
(ii) oscillatory case
This case is further divided into non-chaotic, intermittent, and highly chaotic cases by the choice of reaction networks. In the non-chaotic oscillation case, the spatial pattern of spots is fixed in time. In some reaction networks, there exist propagating waves(See Fig.2a,2b). In spatio-temporal intermittency, the spatial pattern of spots is clearly separated into laminar and burst regions to form some characteristic patterns (See Fig.2c). Those patterns are expected to be in the same class that had been studied extensively in CML. Contrastingly, in the chaotic case, the concentrations of each spot change chaotically in time, around the heteroclinic orbit(See Fig.2d). At any given time, the chemical concentrations vary from cell to cell, but their averages over time are almost identical for each cell. In this case, spots not only replicate but also sometimes annihilate.
(iii) case-I differentiation
In this case, with time, spots differentiate into different types, namely, inner and outer types(See Fig.1a). For the inner type, the set of chemical concentrations converge to a fixed point, with less resource chemicals. Outer type spots exhibit either fixed-point, periodic, or chaotic oscillations, depending on the reaction matrix and the parameter values.
Here, the diversity is larger for the outer type. The inner type has smaller chemical diversity and null local KS entropy, as the dynamics fall onto fixed points. For the parameter values we used, the local KS entropy of the outer part can be either positive or zero, depending on the network (See Fig.3a). (The sign also changes depending on the parameter values.)
(iv) case-II differentiation
Here, spots differentiate into two types, as in case I. Spots of the initial type exhibit chaotic oscillations. The division of these spots produces either the same type of spots with chaotic oscillations or different types with fixed-point dynamics. Here, differentiation is not governed by the inner/outer distinction, but, rather spots with chaotic oscillations appear periodically, with some interval(See Fig.1b). Thus in this case there is pattern formation on two distinct spatial scales, that of the spot size and that of the average distance between two chaotic spots. In this case, the chaotic spots have a large chemical diversity, while other types have much smaller diversities. Complex dynamics are stable within spots of the first type, with positive local KS entropy, while for the other types, the local KS entropy is zero (See Fig.3b).

Let us study the differentiation process in more detail. In both cases, spot types with smaller diversities and fixed-point dynamics are differentiated from the initial type, with high diversity.

In fact, there exist initial conditions for which there results no differentiation of spots. Depending on the initial conditions and their chemical diversity, the first spot will either be a chaotic type or a fixed-point type. In the former case, there will exist differentiation, but in the latter case, there will be none, as such a spot can only replicate. Thus, the differentiation from the first type is irreversible. For initial conditions with sufficiently large chemical diversity, a chaotic spot will appear. However, for smaller initial chemical diversities, the basin structure is so complex that it is not feasible to predict which type of spot will appear initially. If a chaotic spot does appear, the differentiation process we have been discussing will take place.

The difference between the two cases I and II is regard to the ability for spontaneous differentiation. In case I, spots located in the outer region are maintained only through the flow of chemicals from the outmost part. The differentiation from the outer region into the inner region continues if the dynamics in the outer region are chaotic. However, if the dynamics of the outer spots converge to a fixed point, the differentiation is terminated, and as a result growth of the inner part ceases.

On the other hand, in case II, the number ratio of the two types of spots remains nearly fixed (with small fluctuations). In this case, differentiation from chaotic spots continues with some rate, and the spot is distributed in space at some rate. The pattern formed by the different types of spots is independent of the boundary conditions.

The differentiation here is reminiscent of that exhibited by stem cells in biological systems. In a biological stem cell system, initial cell types can either replicate or differentiate into different types. These differentiated cells can themselves differentiate, or they can replicate. Upon such repeated differentiation, multipotency comes to be lost, and eventually cell types that can only replicate are produced. In Ref., such stem cell is found to be natural characteristic of dynamical systems in which there exists cellular structures and cell division. Here it is shown that such stem cell systems naturally appear in a reaction-diffusion system, without assuming the cellular structure in advance.

In case I of our model, these ‘stem-type’ spots are located only in the outer region, while in case II, such ’stem-type’ spots exist within the inside of the pattern. In both cases in which differentiation exhibits in our model, the diversity and local KS entropy decrease as the differentiation progress . Differentiation into a spot with fixed-point type dynamics is irreversible, as mentioned above. This is consistent with the loss of multipotency observed in real cellular systems.

Finally, we examine the parameter dependence of each type of dynamics. Here we consider altering the parameters $`A`$ and $`B`$, while keeping the other parameter values fixed. In the homogeneous case, the bifurcation that results is qualitatively similar to that for the original GS equation. Typical phase diagrams for differentiation in cases I and II are shown in Fig.4. As we see, as the decay rate for $`v_i`$ is decreased, there are bifurcations from a uniform system to a system with a single type of spots to a system with turbulent spatial structure and no clearly defined spots. This set of bifurcations exists quite generally, for a variety of reaction networks $`W`$.

There has been a great deal of effort dedicated to relating reaction-diffusion systems to morphogenesis in biological systems, since the pioneering study of Turing . However these models have not yet been able to describe the complexity observed in biological pattern formation, which includes cell differentiation from stem cells, determination of fixed types, irreversibility and robustness in development. In this work, we have found that the behavior which can be interpreted as representing the irreversible differentiation from stem cells to fixed cell types is displayed quite generally for a particular set of reaction diffusion equations. With regard to the broader context of the general study of reaction-diffusion equations, it is interesting to note that a pattern with two distinct spatial scales is organized in case II. Also it is interesting to note that the behaviors from a variety of reaction networks are classified just into the four cases. It will be important to classify possible types of spatiotemporal patterns that appear when we increase the complexity of the internal reaction dynamics.

Since our model includes only reaction and diffusion, without any other type of mechanism, the present study is also relevant to the study of prebiotic evolution leading to proliferation and diversification of cells. Indeed, the autocatalytic nature of the reaction networks studied here is similar to that of the hypercycle studied by Eigen and Shuster , while the relevance of the spatial structure in resisting parasites and short cuts of network structure  has been discussed  . In contrast to previous studies, the present study explains not only the robust replication process but also the diversification to different cell types.

We would like to thank T.Yomo and C.Furusawa for stimulating discussions. This work is supported by Grants-in-Aid for Scientific Research from the Ministry of Education, Science and Culture of Japan (11CE2006, Komaba Complex Systems Life Project, and 11837004).