text:
# Chandra and HST Confirmation of the Luminous and Variable X-ray Source IC 10 X-1 as a Possible Wolf-Rayet, Black-Hole Binary

## 1. Introduction

IC 10 is a metal-poor ($`Z0.15Z_{}`$; Lequeux et al., 1979), barred dwarf irregular in the Local Group. Although hampered by uncertain reddening corrections, reliable optical and infrared distance estimates place IC 10 at 0.6–0.8 Mpc e.g., Saha et al., 1996; Sakai et al., 1999; Borissova et al., 2000; we adopt 0.7 Mpc, indicating that it is the nearest starburst galaxy to the Milky Way. IC 10 is notable for its vigorous star formation $``$0.03 $`M_{}`$ yr-1 kpc-2; e.g., Hunter & Gallagher, 1986; Thronson et al., 1990; Hunter et al., 1993 and unusually large massive (OB) star population, including a galaxy-wide surface density of Wolf-Rayet (WR) stars $``$$`4`$ times that observed in the most active regions of star formation in any other Local Group galaxy (Massey & Johnson, 1998; Royer et al., 2001; Crowther et al., 2003). This ongoing star formation in IC 10 is expected to result in copious X-ray emission, both from compact high-mass X-ray sources and supernova-heated gas. The galaxy was observed twice using the ROSAT HRI in 1996 (70.3 ks combined exposure), although the X-ray flux in the HRI band is diminished by a factor of $``$5 due to the Galactic column of $`4.8\times 10^{21}`$ cm-2 (Stark et al., 1992) and the expected absorption internal to IC 10 of $``$2$`\times 10^{21}`$ cm-2 (e.g., Yang & Skillman, 1993). One highly significant source was previously detected (hereafter X-1) within the optical extent of IC 10 (Brandt et al., 1997). This source was found to vary by a factor of $``$3 on $``$1 day timescales, with an average absorbed HRI flux of $`4\times 10^{13}`$ erg cm-2 s-1 (or an unabsorbed $`L_{0.12.5\mathrm{keV}}7\times 10^{37}`$ erg s-1). The HRI-derived position (5″ rms) placed X-1 in a region with intense star formation as well as the most massive H I cloud in IC 10. Notably, X-1 was found to lie $``$2″ from \[MAC92\] 17 ($`V=21.76`$), an emission-line source identified as a WR star by Massey et al. (1992) and Crowther et al. (2003), and within $``$8″ of the centroid of a non-thermal radio superbubble ($``$45″/150 pc diameter; Yang & Skillman, 1993). The combination of the high X-ray luminosity and strong variability argued that X-1 is a powerful X-ray binary (containing a neutron star or black hole).

In this Letter we report on Chandra and HST follow-up observations of IC 10 that further characterize the nature of X-1. In particular, these observations more accurately determine the X-ray luminosity and spectrum of X-1 (both of which were uncertain since the ROSAT HRI had essentially no spectral capability and was affected by the large absorption) and confirm the likely association of X-1 and the WR star \[MAC92\] 17 (there are several potential optical counterparts to X-1 within the 5″ radius HRI error circle). In $`\mathrm{\S }`$2 we outline the relevant X-ray and optical observations and analyses, while in $`\mathrm{\S }`$3 we summarize our findings and discuss X-1 in the broader context of black-hole binaries (BHBs).

## 2. Observations and Analysis

### 2.1. Chandra Observations

IC 10 was observed on 2003 March 12 with ACIS-S (the spectroscopic array of the Advanced CCD Imaging Spectrometer; Garmire et al., 2003). To limit telemetry saturation in the case of high background, only chips S2 and S3 in the ACIS-S array were operated, while a 500-pixel subarray was used to lower the frame time to 1.7 s and thus limit potential pileup of X-1. Events were telemetered in Very Faint (VF) mode, and the CCD temperature was $`120^{}`$ C. The optical extent of IC 10 ($``$6$`\stackrel{}{\mathrm{.}}`$3$`\times `$5$`\stackrel{}{\mathrm{.}}`$1) is only marginally cropped by the 2-chip, 500-pixel subarray configuration (16$`\stackrel{}{\mathrm{.}}`$8$`\times `$4$`\stackrel{}{\mathrm{.}}`$1).

Analysis was performed using the CIAO software (v2.3) provided by the Chandra X-ray Center (CXC), but also with FTOOLS (v5.2) and custom IDL software. The data were corrected for the radiation damage sustained by the CCDs during the first few months of Chandra operations using the Charge Transfer Inefficiency (CTI) correction procedure of Townsley et al. (2002).111For details see http://www.astro.psu.edu/users/townsley/cti/. Following CTI-correction, the data were reprocessed with the CIAO tool ACIS\_PROCESS\_EVENTS to remove the standard 0$`\stackrel{}{\mathrm{.}}`$5 pixel randomization and to flag potential ACIS background events using VF mode screening. We performed standard ASCA grade selection in the 0.3–8.0 keV band, excluding all bad columns, bad pixels, VF-mode background events, and cosmic-ray afterglows. We used only data taken during times within the CXC-generated good-time intervals. The background varied by $`<`$30% during the observation, with the average background rate being $`2.2\times 10^7`$ counts s-1 pixel-1 and $`7.0\times 10^7`$ counts s-1 pixel-1 on chips S2 and S3, respectively (in agreement with ACIS quiescent-background calibration measurements).222See http://cxc.harvard.edu/contrib/maxim/bg/index.html. The total net exposure time is 29,191 s.

X-ray sources were identified using the CIAO wavelet algorithm WAVDETECT (Freeman et al., 2002) on the 0.3–8.0 keV (full), 0.3–2.0 keV (soft), and 2.0–8.0 keV (hard) images with a probability threshold of $`10^6`$. In total, 43 sources were detected. The positions, counts, and spectra for these sources were derived using ACIS\_EXTRACT (AE; Broos et al., 2002) with 95% encircled-energy radii \[based on point spread functions (PSFs) from MKPSF and the PSF library\]; the adopted X-ray source positions were found using the matched-filter method within AE.

### 2.2. HST Observations

To compare the X-ray and optical emission, we observed IC 10 with the HST Advanced Camera for Surveys (ACS; Ford et al., 1998) on 2002 October 12. IC 10 was imaged along its major axis ($`50^{}`$) using two ACS pointings separated by 130$`\mathrm{}`$. Each pointing consisted of $`0.5`$-orbit exposures in the $`BVI`$ bands (F435W, F555W, F814W). The standard pipeline data products were used, although we mosaiced the two pointings ourselves due to a problem with the pipeline MULTIDRIZZLE procedure. Optical source positions and magnitudes were measured using SEXTRACTOR (v2.2.1; Bertin & Arnouts, 1996).

### 2.3. Astrometry of Optical and X-ray Images

Due to the small number of X-ray sources that overlap the ACS field of view, we aligned both datasets to the set of high-quality ground-based images available from the NOAO archive for IC 10 (PI: Massey). Briefly, the NOAO images were matched to 55 GSC2.2 sources within 10$`\mathrm{}`$ of IC 10 ($`\mathrm{\Delta }\alpha =0`$$`\stackrel{}{\mathrm{.}}`$20, $`\mathrm{\Delta }\delta =0`$$`\stackrel{}{\mathrm{.}}`$26), giving an rms scatter of 0$`\stackrel{}{\mathrm{.}}`$19; this alignment was confirmed using 11 Tycho 2 stars (0$`\stackrel{}{\mathrm{.}}`$26 rms). Given the optical source density in the IC 10 region and apparent small, systematic field distortions in the NOAO images, we only matched 332 HST sources brighter than $`B=23`$ within a 1$`\stackrel{}{\mathrm{.}}`$5 region of X-1 ($`\mathrm{\Delta }\alpha =5`$$`\stackrel{}{\mathrm{.}}`$32, $`\mathrm{\Delta }\delta =0`$$`\stackrel{}{\mathrm{.}}`$27). This gives an astrometric accuracy of 0$`\stackrel{}{\mathrm{.}}`$10 rms for the HST positions within the vicinity of X-1. Finally, 14 of the 43 significant X-ray sources (including X-1) have plausible optical counterparts, whereas only two false matches are expected. Aligning the X-ray and optical images ($`\mathrm{\Delta }\alpha =0`$$`\stackrel{}{\mathrm{.}}`$44, $`\mathrm{\Delta }\delta =0`$$`\stackrel{}{\mathrm{.}}`$53) yields an astrometric accuracy of 0$`\stackrel{}{\mathrm{.}}`$28 rms for the X-ray positions.

### 2.4. X-ray Analysis

IC 10 X-1 has 4403 counts in the full band and, based on our astrometric solutions, it has a J2000 position of $`\alpha =`$00h20m29$`\stackrel{\mathrm{s}}{\mathrm{.}}`$09, $`\delta =+`$59°16′51$`\stackrel{}{\mathrm{.}}`$95. X-ray contours of X-1 are shown in Fig. 1 overlaid on the HST ACS F814W image. The Fig. 1 inset indicates the Chandra positional uncertainty of X-1, which is an order of magnitude smaller than found previously with the ROSAT HRI. X-1 lies 0$`\stackrel{}{\mathrm{.}}`$23 from J002029.1+591652.1 an $`I_{AB}=21.8`$ WR star identified as \[MAC92\] 17 by Massey et al., 1992 and \[MAC92\] 17A by Crowther et al., 2003, 0$`\stackrel{}{\mathrm{.}}`$33 from J002029.0+591651.6 ($`I_{AB}=24.4`$), 0$`\stackrel{}{\mathrm{.}}`$37 from J002029.1+591651.7 ($`I_{AB}=24.0`$), and 0$`\stackrel{}{\mathrm{.}}`$42 from J002029.0+591651.8 an $`I_{AB}=21.8`$ OB supergiant identified as \[MAC92\] 17B by Clark et al., 2003. Given the astrometric uncertainty in both the HST and Chandra coordinate frames (i.e., 0$`\stackrel{}{\mathrm{.}}`$30, much larger than the individual centroiding errors), X-1 is most likely associated with \[MAC92\] 17A. Although we cannot exclude entirely an association with the latter sources (or fainter undetected sources) based on the astrometric errors, the probability of X-1 lying within 0$`\stackrel{}{\mathrm{.}}`$23 of a WR star by chance is $``$0.05–0.2%.333Narrow-band photometric observations of IC 10 have identified $`100`$ WR stars, of which 26 have been spectroscopically confirmed (Massey & Armandroff, 1995; Royer et al., 2001; Crowther et al., 2003).

In addition to an X-ray point source in Fig. 1, we found faint extended X-ray emission roughly centered on X-1 and co-spatial with the Yang & Skillman (1993) non-thermal radio superbubble (see Fig. 1). This X-ray emission extends out to $``$20$`\mathrm{}`$ in the adaptively smoothed image. Given that the X-ray emission appears azimuthally symmetric to $``$15$`\mathrm{}`$, we extracted soft and hard-band counts around X-1 in 1-pixel annular bins. Due to the large foreground absorption column of IC 10, we used the 0.8–2.0 keV band to reduce background where we detect no signal. For comparison, we generated the on-axis ACIS-S PSF using CHART,444See http://asc.harvard.edu/chart/. extracting counts in an identical manner to the data. The radial surface-brightness distributions and 1$`\sigma `$ errors for the data are shown in Fig. 2 and demonstrate that X-1 lies systematically above the background out to radii of $``$35 pixels ($`17\mathrm{}`$) and is clearly extended compared to the PSF model. The hard-band profile is consistent with a point source and background only.

We find 175 (129) counts in the full (soft) band between radii of 7–35 pixels (note that contamination from X-1 should be $`<1`$% at such radii). The fact that the extended emission is present almost exclusively in the soft band further supports the reality of the extent. An emission model consisting of both a point source and a uniform disk with a $`17\mathrm{}`$ (60 pc) radius convolved with the PSF represents the 0.8–2.0 keV surface-brightness distribution well, although there is clearly still some residual scatter due to possible clumping. The flux ratio of the extended disk component to the point source is $`2`$%. The relative softness of the extended emission is consistent with a physical origin as supernova-heated hot gas (see also the X-ray spectral analysis below).

Brandt et al. (1997) found that X-1 varied by a factor of $``$3 during their ROSAT HRI observation but did not have the statistics or the time sampling to constrain the variability further. With $``$4400 counts over a contiguous 29.2 ks, our light curve probes the short-term temporal properties better. A Kolmogorov-Smirnov test indicates that X-1 is variable at the $`>99.9`$% confidence level. To derive the amplitude and characteristic timescale of variability, we used the Bayesian Block method provided in the contributed CIAO software package SITAR.555See http://space.mit.edu/CXC/analysis/SITAR/index.html Fig. 3 shows a binned light curve of X-1, overplotted with a series of contiguous ”blocks” within which the count rate is modeled as a constant at the 99.9% confidence level. X-1 has peak, trough, and average count rates of 0.182 ct s-1, 0.088 ct s-1, and 0.151 ct s-1, respectively, and exhibits rapid variability by as much as a factor of $``$2 in $``$3000 s. The hardness ratio, however, shows no significant spectral variability. Since a significant fraction of X-ray binaries are X-ray pulsars with pulse periods ranging from 0.01–1,000 s (e.g., White et al., 1995), we generated power spectra to search for pulsations using the robust $`z_1^2`$ method (e.g., Buccheri et al., 1983). These power spectra also allowed us to search for other periodic phenomena (e.g., eclipses or regular X-ray bursts). No periodic behavior was found.

The X-ray spectra were analyzed using XSPEC (v11.2; Arnaud, 1996). Unless stated otherwise, spectral parameter errors are for the 90% confidence level, assuming one parameter of interest. The X-ray fluxes and absorption-corrected luminosities for X-1 were calculated from the spectral fitting. For the detected X-1 count rates above, pile-up is estimated to be $``$15%. Thus, spectral analysis of this source must proceed with care. To model the spectrum in the presence of such pile-up, we used the forward-modeling tool LYNX developed at PSU (Chartas et al., 2000). The ACIS spectrum of X-1 was initially fit with an absorbed power-law model, giving best-fit values of $`N_\mathrm{H}=(5.1\pm 0.4)\times 10^{21}`$ cm-2 and $`\mathrm{\Gamma }=1.77\pm 0.08`$ ($`\chi ^2=97.6`$ for 82 degrees of freedom). The $`\chi ^2`$ value and apparent systematic residuals (bottom panel in Fig. 4) indicate the model provides a poor fit to the data. Given that we detect a low-level extended component around X-1, we added an absorbed thermal plasma component (mekal; e.g., Mewe et al., 1985) to our spectral model. The resulting best-fit values are $`N_\mathrm{H}=(6.0_{0.08}^{+0.20})\times 10^{21}`$ cm-2 (in agreement with the expected Galactic$`+`$intrinsic absorption), $`\mathrm{\Gamma }=1.83_{0.11}^{+0.07}`$, and $`kT=1.49_{0.12}^{+0.15}`$ keV ($`\chi ^2=61.6`$ for 80 degrees of freedom). The abundances were not well constrained and were fixed at $`Z=0.15Z_{}`$ (see $`\mathrm{\S }`$1). While this model provides a significant improvement to both the $`\chi ^2`$ value ($`>`$99.99% based on the $`F`$-test) and the apparent systematic residuals below $``$2 keV, there are still notable residuals around 2–4 keV (see middle panel in Fig. 4). These residuals, particularly around 2.25 keV and 3.65 keV, are perhaps due to marginally-resolved emission lines from an X-ray-photoionized wind (e.g., as seen in Cyg X-3; see $`\mathrm{\S }`$3). The 0.5–8.0 keV absorbed flux and unabsorbed luminosity of X-1, corrected for pileup, are $`1.57\times 10^{12}`$ erg cm-2 s-1 and $`1.50\times 10^{38}`$ erg s-1, respectively. The thermal component contributes $``$20% of the total absorbed flux and $``$27% of the total unabsorbed luminosity, which is substantially more than was estimated from fitting the radial profile. This suggests that either much of this component arises from X-1 itself or that the extended component is strongly centrally peaked rather than uniformly distributed. We also analyzed the spectrum of the extended X-ray component (3$`\stackrel{}{\mathrm{.}}`$5–17″annulus around X-1). Given the small number of counts in this region, we modeled its spectrum using the Cash statistic (Cash, 1979) instead of $`\chi ^2`$. An absorbed mekal model with fixed $`N_\mathrm{H}=5.9\times 10^{21}`$ cm-2 (as found for X-1) gives abest-fit temperature of $`kT=0.87_{0.21}^{+0.17}`$ keV. The 0.5–8.0 keV absorbed flux and unabsorbed luminosity of this extended component are $`1.73\times 10^{14}`$ erg cm-2 s-1 and $`3.24\times 10^{36}`$ erg s-1, respectively. This flux is consistent with our radial-profile fitting and supports an origin as hot gas associated with the radio superbubble.

## 3. Discussion

This Chandra observation of IC 10 has resulted in improvements by factors of $``$15 in terms of spatial resolution and astrometric accuracy and $``$10 in terms of statistics for X-1, compared to previous ROSAT HRI exposures, allowing substantially better constraints to be placed on the nature of this object. The combination of the high luminosity and strong variability clearly demonstrates that X-1 is a powerful X-ray binary, containing a black hole or neutron star. Its likely physical association with \[MAC92\] 17A implies that the progenitor of X-1 must have evolved even more rapidly and been more massive than \[MAC92\] 17A (WNE stars typically have $`M40`$–50$`M_{}`$ and $`t_{\mathrm{life}}5`$ Myr; e.g., Maeder & Meynet, 1994). Thus, we expect X-1 to be a BHB (although we note that there are clear counter examples; see, e.g., Kaper et al., 1995). Optical spectroscopic monitoring of X-1 to search for variability would be one of the best methods to prove the WR-BHB interpretation; Clark et al. (2003) find hints of variability in the $`\lambda `$4686 feature of \[MAC92\] 17A, but longer dedicated observations are needed to secure this finding.

The X-ray spectrum of X-1 is acceptably fit with a $`\mathrm{\Gamma }=1.83`$ power law and $`kT=1.49`$ keV thermal plasma (potentially indicative a BHB in a hard state; e.g., McClintock & Remillard, 2003). The residuals of X-1 in the 2–4 keV range hint at further spectral complexity, perhaps from emission lines or absorption features. Longer Chandra or XMM-Newton observations of X-1 would also be useful to quantify the origin of the spectral residuals, as well as to place additional constraints on the nature of the X-ray emission (e.g., to rule out further the presence of pulsations or bursts).

The properties of X-1 bear several similarities to the Galactic X-ray binary Cyg X-3 in terms X-ray luminosity, spectrum, and variability (e.g., Kitamoto et al., 1994; Liedahl & Paerels, 1996; Predehl et al., 2000) as well as donor type (WNE vs. WN7; van Kerkwijk et al., 1996; Crowther et al., 2003). Moreover, Cyg X-3 has a complex X-ray spectrum including numerous emission lines and a 9 keV edge that are attributed to an X-ray-photoionized wind; the residuals seen from X-1 may have a similar origin. Thus IC 10 X-1 is the only other confirmed example of the short-lived WR X-ray binary.

The discovery of soft extended emission co-spatial with the large non-thermal radio superbubble surrounding X-1 is important for determining the nature of the bubble, although the poor statistics obtained here do not allow strong constraints. Obvious comparisons can be made with 30 Doradus (Dennerl et al., 2001), as well as with SS433/W50 (Safi-Harb & Oegelman, 1997) and IC 342 X-1 (Roberts et al., 2003), all of which have similar physical extents and radio/X-ray luminosities as the bubble in IC 10. While the radio and X-ray emission from the IC 10 bubble is consistent with multiple supernovae, the superbubble does lie off the $`\mathrm{\Sigma }`$$`d`$ relation (Yang & Skillman, 1993), suggesting that something else (perhaps X-1) may be powering the expansion. For instance, in the case of W50, SS 433 is believed to contribute to the expansion. High-resolution radio imaging of IC 10 could be used to search for any jets associated with X-1 and explore whether X-1 has had any influence on the nonthermal superbubble.

We thank L. Townsley and M. van Kerkwijk for several helpful discussions, G. Chartas for help with the forward-modeling of the X-ray spectrum, M. Muno for help with variability constraints, S. Clark and P. Crowther both for interesting discussions regarding X-1 and communicating their results prior to publication. We gratefully acknowledge the financial support of CXC grant GO3-4112X (FEB, WNB), STScI grant HST-GO-09683.01A (FEB, WNB) and NASA LTSA grant NAG5-13035 (WNB).


format: text

text:
# SPECTRA FROM FORMING REGION OF THE FIRST GALAXIES : THE EFFECT OF ASPHERICAL DECELERATION

## 1 INTRODUCTION

The epoch of cosmic reionization is one of the most important epochs to understand in the investigation of the history of the universe, because it contains the imprint of important information on cosmic structure formation and evolution. This is because the reionization of the intergalactic medium (IGM) occurs promptly (e.g. Gnedin & Ostriker 1997) and the cosmic ionizing background radiation can penetrate the dense HI condensations of virialized objects (Miralda-Escudé, Haehnelt, & Rees 2000; Barkana & Loeb 1999).

The epoch of the reionization observationally relates to the Gunn-Peterson trough (Gunn & Peterson 1965). The significantly ionized IGM in the observed spectra of QSOs at $`z<5.6`$ (Dey et al. 1998; Hu, Cowie, & McMahon 1998; Spinrad et al. 1998; Weymann et al. 1998) implies that most of the intergalactic hydrogen in the Universe has been reionized before $`z5.6`$. Some standard cosmological models predict that reionization takes place around a redshift of $`z6`$. However, according to the results of the WMAP project, Kogut et al. (2003) may have discovered the high redshift reionization of the Universe ($`z20`$). To resolve this discrepancy, the direct detection of the emission from the first objects is very important. This is because that emission reflects the physical and chemical properties of the first star clusters (Shchekinov 1991; Kamaya & Silk 2002) and/or protogalaxies. Once we know the origin of the first emission, the ionization rate at that time is investigated. Therefore, new detailed studies and confirmation of the previous studies of the first emission in another way are important for the subsequent progress of the research fields.

It is interesting to detect just the first Ly$`\alpha `$ photons since the star formation rate and/or the ionization structure of the interstellar medium are examined from that intensity. Using standard cosmological models, cosmic recombination takes place at $`z10^3`$. After that epoch, the Universe becomes predominantly neutral. When the first galaxies and/or star clusters form, they are surrounded by the neutral IGM. One might naively expect the damping wing of the Ly$`\alpha `$ trough to eliminate any trace of a Ly$`\alpha `$ emission line in the observed source spectra. However, the Ly$`\alpha `$ photons are scattered and diffused redward in frequency owing to the Hubble expansion of the surrounding HI IGM (Loeb & Rybicki 1999 ; hereafter LR99). When the cumulative frequency-shift becomes significantly large, they can escape and travel freely toward the observer.

Although we investigate a similar problem to that in LR99, the radiative transfer model is improved in the two ways : (1) Formulation of the geometrical scattering (i.e. test photons scatter spatially three dimension); and (2) Formulation of the Ly$`\alpha `$ opacity. Thanks to these two improvements and the confirmation of the previous studies, we can study more complex physical conditions in the expanding HI gas. Since the large scale structure (LSS) forms filaments, this suggests the initial stage of the gravitational deceleration for the LSS favors an asymmetrical expansion of the HI halo. Indeed, the inhomogeneity of the self-gravitating system can grow in a free-fall time (e.g. Lin, Mestel, & Shu 1965). As long as the expansion time-scale is comparable to the free-fall time scale, the asymmetrical deceleration of the IGM occurs in the time-scale of the Hubble expansion, although the initial asymmetry is obviously very weak. As the first and the simplest hypothesis to investigate the effect of the inhomogeneity via the gravitational deceleration for the LSS formation, we adopt the aspherical expansion model of HI IGM. The physical implication drawn from the current analysis will suggest to us further research projects for the first Ly$`\alpha `$. Since the lumpiness and/or small-scale irregularity may be important for a realistic treatment of the first Ly$`\alpha `$, we will investigate these effects in subsequent papers.

This paper is organized as follows: In §2 we describe the formalism of the Ly$`\alpha `$ line transfer problem commenting on the differences between our approach and LR99. In §3 we outline our Monte Carlo simulation method. In §4 we discuss our numerical results including the effects of a more complete representation of the Ly$`\alpha `$ opacity, the effect of asphericity in the expanding IGM, and the sensitivity of our results to redshift and cosmological parameters. In §5 we summarize our study.

## 2 FORMULATION

We formulate the problem according to LR99, while the HI IGM is hypothesized to expand aspherically. In general, the exactly spherical expansion is not expected as long as the IGM decelerates to form LSS. This is recognized because the free-fall time, $`t_{ff}`$, is comparable to the Hubble expansion time, $`t_h`$, as explained in §1. Then, to find the effect of the asymmetry, we shall set the most simple assumption to the recession velocity, spheroidal expansion: $`v=H_s\left[\epsilon ^2\left(x^2+y^2\right)+z^2\right]^{1/2}`$, where $`v`$ is recession velocity, $`H_s`$ is the Hubble expansion rate at the source redshift of $`z_s`$, and $`\epsilon `$ is ellipticity or decreasing rate of the recession velocity to the pure Hubble expansion. We shall estimate the allowed order of $`\epsilon `$. The ratio of the two time-scales is denoted as

$$\frac{t_{ff}}{t_h}\frac{\pi }{2}\mathrm{\Omega }_M^{1/2}\delta ^{1/2}.$$
(1)
Here, $`\delta `$ is the ratio of the density of the halo to the mean density of the Universe at $`z`$. Of course, this ratio must not be different from unity very much as long as IGM expands. Since the density of the aspherically decelerating HI halo is larger than that of the spherically expanding halo, we find $`\delta =1/\epsilon `$ with mass conservation. Fixing the ratio of Eq.(1) to be unity and postulating $`\delta 6`$, for example, we estimate $`\epsilon 0.2`$ for $`\mathrm{\Omega }_M=0.4`$. Thus, we recognize that order of $`\epsilon `$ is about 0.1. The case of $`\epsilon =0.1`$ is considered as an extreme case. Obviously, $`\epsilon =1`$ corresponds to LR99.

Although $`\epsilon `$ should approach to unity differentially from small scale to large scale, we assign $`\epsilon `$ a single value for simplicity. Fortunately, the size of the concerned region, $`r_{}`$, which is a characteristic distance from the Ly$`\alpha `$ photon source as defined later, is $``$ Mpc as estimated in LR99. That is, $`r_{}`$ corresponds to the size of groups of galaxies and the region at $`r_{}`$ can be regarded as decelerating region although the allowed deceleration should be small (but not zero). Expansion law of the spheroid approaches to that of the spherical model at larger scale than $`r_{}`$. Properties of test photons are also determined at $`r_{}`$. Hence, the single value hypothesis for $`\epsilon `$ can present the physical insight into the effect of the deceleration for the LR halo. In realistic condition, by the way, $`\epsilon `$ inside $`r_{}`$ should decrease toward the center differentially. Then, our finding characteristic effect of the deceleration may be enhanced because the number of scattering can increase owing to the deceleration as shown in this paper.

Let $`I=I(\nu ,𝛀,𝐫)`$ be the specific intensity (in photons $`\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{Hz}^1`$), where $`\nu `$ is comoving frequency, $`𝐫=(x,y,z)`$ is the radial vector from the source galaxy, and $`𝛀`$ is the scattering direction. The spatial components of the direction cosines of $`𝛀`$ are $`(u,v,w)`$, respectively. With the assumption of isotropic coherent scattering, the comoving transfer equation for a resonance line is given by

$$𝛀I+\alpha \psi (\epsilon ,w)\frac{I}{\nu }=\chi _\nu \left(JI\right)+S.$$
(2)
Here, $`\nu `$ is the shift of the frequency $`\nu _\alpha \nu _{\mathrm{photon}}`$ where $`\nu _\alpha `$ is the resonant frequency and $`\nu _{\mathrm{photon}}`$ is the photon frequency; $`\chi _\nu `$ is the scattering opacity at $`\nu `$ (eq. (4) in LR99); $`J=(1/2)_1^1I𝑑\mu `$ is the mean intensity; $`S`$ is the emission function (in photons $`\mathrm{cm}^3\mathrm{s}^1\mathrm{sr}^1\mathrm{Hz}^1`$, eq. (3) in LR99); $`\alpha H_s\nu _\alpha /c`$ is the increasing rate of the frequency shift per unit distance in the spherically expanding IGM (eq. (2) in LR99); and

$$\psi (\epsilon ,w)=\left[\epsilon ^2+\left(1\epsilon ^2\right)w^2\right]^{1/2}$$
(3)
represents the effect of aspherical expansion.

In LR99 the scattering opacity $`\chi _\nu `$ was approximated as $`\chi _\nu \beta /\nu ^2`$. In the current paper, we formulate the opacity in the following more complete form:

$`\chi _\nu `$ $`=`$ $`{\displaystyle \frac{\beta /\epsilon }{\mathrm{\Lambda }_\alpha ^2/16\pi ^2}}\mathrm{for}\nu ^2\mathrm{\Lambda }_\alpha ^2/16\pi ^2,`$ (4)
$`=`$ $`{\displaystyle \frac{\beta /\epsilon }{\nu ^2+\mathrm{\Lambda }_\alpha ^2/16\pi ^2}}\mathrm{for}\nu ^2\mathrm{\Lambda }_\alpha ^2/16\pi ^2,`$ (5)
$`=`$ $`{\displaystyle \frac{\beta /\epsilon }{\nu ^2}}\mathrm{for}\mathrm{\Lambda }_\alpha ^2/16\pi ^2\nu ^2\&\nu \nu _\alpha ,`$ (6)
$`=`$ $`\left(1{\displaystyle \frac{\nu }{\nu _\alpha }}\right)^4{\displaystyle \frac{\beta /\epsilon }{\nu ^2}}\mathrm{for}\nu \nu _\alpha ,`$ (7)
$`=`$ $`0\mathrm{for}\nu >\nu _\alpha ,`$ (8)
where $`\nu _\alpha =2.47\times 10^{15}\mathrm{Hz}`$ is the Ly$`\alpha `$ frequency, $`\mathrm{\Lambda }_\alpha =(8\pi ^2e^2f_\alpha /3m_ec\lambda _\alpha ^2)=6.26\times 10^8\mathrm{s}^1`$ is the rate of spontaneous radiative decay from $`2p`$ to $`1s`$ energy levels of hydrogen, $`f_\alpha =0.4162`$ is the oscillator strength, $`\lambda _\alpha =1216\mathrm{\AA }`$ is the wavelength of the Ly$`\alpha `$ line, and $`\beta `$ is $`1.5\mathrm{\Omega }_\mathrm{b}h_0^2(1+z_\mathrm{s})^3`$ cm-1 Hz2 (eq. (5) in LR99). In the initial scattering phase, Eq.(4) is important since the frequency shift can be smaller than $`\mathrm{\Lambda }_\alpha /4\pi `$. When the photon diffuses and the frequency shift is comparable to $`\nu _\alpha `$, Eq.(7) must be adopted as $`\chi _\nu `$.

We normalize the frequency shift and radius by convenient quantities like LR99. An appropriate frequency scale, $`\nu _{}`$, is estimated from $`\tau _\nu 1`$. Here, $`\tau _\nu `$ is the Ly$`\alpha `$ optical depth from the source to the observer, and is defined as

$$\tau _\nu =_0^{\mathrm{}}\chi _{\nu +\alpha \psi r}𝑑r.$$
(9)
The characteristic frequency shift of $`\nu _{}`$ is measured along the direction of $`z`$ (i.e. $`w=1`$). Eq. (9) with Eq.(7) can be described as

$$_0^{l_0}\left(1\frac{\nu _{}+\alpha r}{\nu _\alpha }\right)^4\frac{\beta /\epsilon }{\left(\nu _{}+\alpha r\right)^2}𝑑r=1,$$
(10)
where $`l_0(\nu _\alpha \nu _{})/\alpha `$ is a characteristic distance for the photon frequency to reach zero. We note that $`\nu _{}`$ always depends on $`\beta `$ and $`\epsilon `$. The appropriate radius, where the frequency-shift due to the Hubble expansion produces the frequency shift of $`\nu _{}`$, is given by $`r_{}=\nu _{}/\alpha `$. For example, with $`\epsilon =1`$, we find $`\nu _{}=9.64\times 10^{12}\mathrm{Hz}`$ and $`r_{}=0.77`$ Mpc at $`z_s=10`$ with a set of cosmological parameters of $`\mathrm{\Omega }_b=0.05,\mathrm{\Omega }_M=0.4,\mathrm{\Omega }_\mathrm{\Lambda }=0.6,`$ and $`h_0=0.65`$. 111We adopt these cosmological parameters because of the comparison to LR99. Although the recent results of WMAP present slightly different cosmological parameters (Spergel et al. 2003), fortunately our results are not altered very much as shown in §4.4. The normalized frequency shift $`\stackrel{~}{\nu }`$ and distance $`\stackrel{~}{r}`$ are then defined to be $`\stackrel{~}{\nu }\nu /\nu _{}`$ and $`\stackrel{~}{r}r/r_{}`$.

For the reader’s convenience, we present the Lyman$`\alpha `$ opacity $`\chi _\nu `$ $`\mathrm{cm}^1`$ in Fig.1 as a function of $`\stackrel{~}{\nu }`$. The solid line denotes our $`\chi _\nu `$, and the dashed line denotes the previously adopted one from LR99.

According to the definition of $`\stackrel{~}{\nu }`$, the photon frequency $`\nu _{\mathrm{photon}}`$ becomes small as $`\stackrel{~}{\nu }`$ increases and reaches zero at $`\stackrel{~}{\nu }=\nu _\alpha /\nu _{}2.56\times 10^2`$. At $`\stackrel{~}{\nu }=0`$, the frequency of the photon is the Ly$`\alpha `$ frequency $`\nu _\alpha `$. The differences depicted in the figure 1 are reflected in our results, while they are usually small.

## 3 Algorithm of the Calculation

LR99 solved a dimensionless transfer equation in the diffusion regime, which is valid in the region $`\stackrel{~}{\nu }\stackrel{~}{r}`$. They also studied the numerical solution of the transfer equation using a Monte Carlo technique and found that it agreed with the diffusion equation in the low frequency shift region. We also solve the transfer equation numerically by using a Monte Carlo method very similar to LR99. The number of test photons is $`10^8`$, the same as in LR99, and permits us to compare our results directly to theirs. During the calculation two photon decay is always neglected since it is difficult to destroy the Lyman $`\alpha `$ photons transferred in the surrounding HI. This is because the event rate of two photon decay is determined by the process of a transition between $`2s`$ \- $`2p`$ owing to the radiative collision with an electron (Breit & Teller 1940; Spitzer & Greenstein 1951). As long as the relic electron number density is very small, the event rate of two photon decay should be negligible.

LR99 adopted the initial shift method in the first step of each simulation. That is, the initial shift of $`\stackrel{~}{\nu }`$ is 0.01 and $`\stackrel{~}{𝐫}`$ is determined stochastically by means of the diffusion solution of LR99. We have confirmed independently that the initial shift method is very reasonable. This is because the initial optical depth is always much larger than unity. Hence we also adopt that method in all first time steps of the simulations.

1. 3D scattering
LR99 assumed isotropic, coherent scattering through a spherically symmetric, neutral medium such that the scattering could be described by a single angular parameter $`\mu =\mathrm{cos}\theta `$ relative to the radius $`r`$. We generalize to three dimensional scattering. The scattering processes in 3D coordinates is described by two parameters of $`\mu =2R_11`$ and $`\varphi =\pi \left(2R_21\right)`$. Here, $`\varphi `$ is the change in azimuth, and $`R_1`$ and $`R_2`$ are random numbers on the interval of (0, 1). Defining $`\mu `$ as $`\mathrm{cos}\alpha `$, we measure $`\alpha `$ in the plane determined by the old and new photon paths.

The direction cosines of $`u_{i+1},v_{i+1},w_{i+1}`$ after the scattering are determined by using the direction cosines of the incident photon, $`u_i,v_i,w_i`$:

$`u_{i+1}`$ $`=`$ $`\mathrm{sin}\alpha \left(\mathrm{cos}\varphi w_iu_i\mathrm{sin}\varphi v_i\right)/\left(1w_i^2\right)^{1/2}+\mathrm{cos}\alpha u_i,`$
$`v_{i+1}`$ $`=`$ $`\mathrm{sin}\alpha \left(\mathrm{cos}\varphi w_iv_i+\mathrm{sin}\varphi u_i\right)/\left(1w_i^2\right)^{1/2}+\mathrm{cos}\alpha v_i,`$ (11)
$`w_{i+1}`$ $`=`$ $`\mathrm{sin}\alpha \mathrm{cos}\varphi \left(1w_i^2\right)^{1/2}+\mathrm{cos}\alpha w_i.`$
As $`|w_i|`$ approaches unity, according to Witt (1977), the set of equations (11) breaks down and is replaced by

$$u_{i+1}=\mathrm{sin}\alpha \mathrm{cos}\varphi ,v_{i+1}=\mathrm{sin}\alpha \mathrm{sin}\varphi ,w_{i+1}=\mathrm{cos}\alpha w_i\mathrm{for}1|w_i|10^{10}.$$
(12)
2. Opacity selection
LR99 used a single approximate form for the opacity as a function of frequency (i.e. our equation (6) for $`\epsilon =1`$). We use the more complete form for the opacity by Equations (4-8). Selecting the appropriate opacity in Eqs.(4-8) according to $`\stackrel{~}{\nu }`$, we evaluate the path length along the ray, $`\stackrel{~}{l}`$, and the maximum optical depth: $`\tau _{\mathrm{max}}_0^{\mathrm{}}\chi _{\nu +\alpha \psi l}𝑑l`$. Once the optical depth $`\tau `$ for a test photon becomes larger than the maximum optical depth $`\tau _{\mathrm{max}}`$, we say that it escapes from the LR halo. Here, $`\tau `$ is determined as $`\mathrm{ln}R_{\mathrm{sca}}`$ where $`R_{\mathrm{sca}}`$ is a random number between zero and one.
3. Aspherical expansion
LR99 assumed that the neutral IGM surrounding the Ly$`\alpha `$ point source undergoes uniform Hubble expansion. Locally we do not expect this to be the case. The cumulative frequency-shift and $`\tau _{\mathrm{max}}`$ differ depending on the scattering direction. According to Eq.(3), the increment of the frequency shift of the scattering photon is represented by

$$\stackrel{~}{\nu }_{i+1}=\stackrel{~}{\nu }_i+\psi (\epsilon ,w_i)\stackrel{~}{l}.$$
(13)

## 4 Results

We first compare the results of our 3D simulation technique and improved approximation to the Ly$`\alpha `$ opacity to those of LR99. We then expand upon that previous work by considering the aspherical expansion of the IGM. Recently Spergel et al. (2003) have used WMAP data to determine more precise values for the cosmological model parameters. WMAP data also suggest that the epoch of reionization may be at a redshift larger than 10 (Kogut et al. 2003). Thus we conclude with a discussion of the effect of redshift and the cosmological model parameters on the observable Ly$`\alpha `$ line profile.

All figures presenting the results are made by counting the number of escape photons in the same interval of normalized frequency shift of $`\stackrel{~}{\nu }`$ in log scale. We choose the interval of $`10^3`$ from $`\mathrm{log}\stackrel{~}{\nu }=1.0`$ to $`2.0`$. For a clear presentation of the effect of the aspherical expansion, we count the photons all over the complete solid angle of 4$`\pi `$.

### 4.1 Geometrical effect

In Fig.2 the geometrical effect of using our 3D numerical method for the scattering is compared to the 2D treatment of LR99. The same scattering coefficient used in LR99 is adopted with $`\epsilon =1`$. The horizontal axis is the normalized frequency, $`\stackrel{~}{\nu }`$, and the vertical axis is the escape fraction of test photons. The 2D treatment of scattering by LR99 is the solid line; while our 3D treatment of the scattering is the dashed line. We observe no difference between the lines. This agreement comes from a simple approximation that the neutral IGM around the Ly$`\alpha `$ source has spherical symmetry. Furthermore, it shows our 3D scheme has sufficient resolution to compare directly with the results of LR99. Hereafter, we compare our results with the spherically symmetric case in 3D coordinates, which is called 3D LR99.

### 4.2 Opacity evaluation effect

How about the effect of opacity? Fig.3 shows the results with the fully formulated opacity and $`\epsilon =1`$ (solid line) and that of 3D LR99 (dashed line), using the same definition of both coordinates as in figure 2. It is found that there are almost no difference between the line profiles of ours and LR99’s. This result confirms that the Ly$`\alpha `$ opacity is approximated well in LR99.

Transforming $`\stackrel{~}{\nu }`$ to the observational wavelength of $`\lambda _{\mathrm{obs}}`$ by

$$\lambda _{\mathrm{obs}}=\left(1+z_s\right)\lambda _\alpha \left[1\left(\frac{\nu _{}}{\nu _\alpha }\right)\stackrel{~}{\nu }\right]^1,$$
(14)
we compare our model with 3D LR99. In Eq.(14), the term $`\left[1\left(\nu _{}/\nu _\alpha \right)\stackrel{~}{\nu }\right]^1`$ represents the effect of the cumulative frequency-shift. The bin of the wavelength for photon counts corresponds to that of the frequency. We show the normalized luminosity of the Ly$`\alpha `$ line per unit wavelength, $`\stackrel{~}{L}_\alpha (\lambda )`$, in Fig.4. This can be converted into the observational $`L_\alpha (\lambda )`$ by multiplying the steady emission rate of Ly$`\alpha `$ photons at the source, $`\dot{N_\alpha }`$ (in photons $`s^1`$), by the Ly$`\alpha `$ photon energy, $`h\nu _\alpha `$.

The observed Ly$`\alpha `$ wavelength emitted at $`z_s=10`$ is $`\lambda _{\mathrm{obs}}=1.337\mu \mathrm{m}`$. As shown in the figure, the peaks of the line profiles shift to $`\lambda _{\mathrm{peak}}1.339\mu \mathrm{m}`$. This shift from the line center is $`2\mathrm{nm}`$. We need a wavelength resolution of $`R7.8\times 10^2`$ if it is to be detected. The full widths at half maximum (FWHM) of our line profile and 3D LR99 are $`5.27\times 10^3\mu \mathrm{m}`$ and $`4.94\times 10^3\mu \mathrm{m}`$, respectively.

### 4.3 Aspherical expansion effect

In this section the effect of the aspherical expansion of the HI IGM is summarized. Fig.5 presents the normalized luminosities of Ly$`\alpha `$ line profiles for various values of the parameter $`\epsilon `$ that characterizes the asphericity of the IGM. Here, the solid line is $`\epsilon =1.0`$, long dashed line is $`\epsilon =0.9`$, short dashed line is $`\epsilon =0.7`$, dotted line is $`\epsilon =0.5`$, long dash-dotted line is $`\epsilon =0.3`$, and short dash-dotted line is $`\epsilon =0.1`$.

What is the distinct property of the line profile with $`\epsilon 1.0`$? The characteristic quantities are summarized in Table 1, where some results of other $`\epsilon `$’s are also presented. The last column in Table 1 is the necessary resolution to distinguish $`\lambda _{\mathrm{peak}}`$ from the Ly$`\alpha `$ line center of $`(1+z_s)\lambda _\alpha `$.

According to table 1, the distinct effects of the asphericity are the shift of the peak and the broadening of the FWHM of the line profiles. These are clearly seen in the line profile with $`\epsilon =0.1`$. If there is strong asphericity of the expanding HI IGM, we will detect broadening of FWHM and the shift of the peak wavelength of the line profile. This could be observable with the JWST and signal the presence of the LR halo.

For increasing asphericity, i.e. a decreasing $`\epsilon `$, the optical depth for photons to scatter becomes larger than in the spherical model. Thus, although the typical shift per scattering due to the expansion of the universe becomes small, the number of scattering events increases, thus increasing the probability that a photon will be shifted far from the line center.

### 4.4 The Effect of Redshift and Cosmological Parameters

WMAP data suggest that cosmic reionization may take place at $`z_s20`$ (Kogut et al. 2003). In Fig.6 we show the results of Ly$`\alpha `$ line transfer on the Ly$`\alpha `$ line profiles from a source at $`z_s=20`$ (solid line) compared to a source at $`z_s=10`$ (dashed line). The coordinates are the same as in Fig.3.

We find the line profile of $`z_s=20`$ coincides with that of $`z_s=10`$ very well. This comes from the fact that most of the Ly$`\alpha `$ photons can escape from the IGM when their cumulative frequency shift becomes $`\nu _{}`$, and it takes place at the distance of $`r_{}`$ from the source.

While $`\nu _{}`$ and $`r_{}`$ depend upon redshift, the shape of the two line profiles as functions of the normalized $`\stackrel{~}{\nu }`$ are very similar. The effect of redshift is not significant except in the normalizations. This result presents the useful implication that when we examine line profiles at a variety of redshifts, it is sufficient for us to know the characteristic frequency shift and radius at each redshift. Hence, we summarize the redshift dependence of $`\nu _{}`$ in Fig.7. Here, the solid line is $`\nu _{}`$ of Eq.(10), while the dashed line is $`\nu _{}^{\mathrm{app}}`$ which is adopted in LR99 and determined by Eq.(6) and Eq.(9). For simplicity, we present the cases with $`\epsilon =1`$. We hope these results are useful for the heterorelation of future observations.

We also simulate the line transfer with the cosmological parameters found by Spergel et al. (2003) in their analysis of WMAP data, who suggest $`h_0=0.71_{0.03}^{+0.04}`$, $`\mathrm{\Omega }_\mathrm{M}h_0^2=0.135_{0.009}^{+0.008}`$, $`\mathrm{\Omega }_\mathrm{b}h_0^2=0.0224_{0.009}^{+0.009}`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=1.0_{0.02}^{+0.02}\mathrm{\Omega }_\mathrm{M}`$. Fig.8 presents the results at $`z_s=10`$. The solid line adopts the cosmological parameters of LR99, while the dashed line uses those of Spergel et al. (2003). We find the effect of the cosmological parameters on the line profile is small.

## 5 SUMMARY

In this work we extended LR99 by improving the numerical treatments of the Ly$`\alpha `$ opacity and scattering processes. As part of these improvements, the effect of asphericity in the expanding IGM on the Ly$`\alpha `$ profiles was examined. The results confirm LR99 very well when the aspherical deceleration is small. The effect of aspherical expansion on the line profile is reflected in a broadening of the FWHM and a shift of the peak intensity to longer wavelengths over that in the spherical model. The effect of the redshift of reionization is reflected in the normalization of the frequency shift and escape distance from the source rather than in the shape of the normalized line ratio. Our results were insensitive to small changes in the cosmological parameters, with little difference in the results found using the cosmological model of LR99 and that derived recently from WMAP data. To detect both properties will be important to find the initial stage of the formation of LSS. Fortunately, it requires a spectroscopic resolution of $`R1001000`$, which is within the design characteristics for future observational projects like JWST.

We are grateful to the referee for his/her advisable comment and excellent refereeing, which improve our paper very much in both content and presentation. We thank Profs. S. Inagaki, S. Mineshige, T. Totani for their continuous encouragement.


text:
# Adaptability and Diversity in Simulated Turn-taking Behaviour

## 1 Introduction

Dynamical systems can theoretically simulate behaviour produced over time with interactions between various entities. This approach, based on embodied cognition , has a different perspective from the traditional AI approaches. That is, representations are not given as symbols in advance but are only realized, by the dynamics, over time . Cognitive structure is characterized by geometrical and flow patterns in an adequate phase space. As well as being characterized by attractor types (e.g., fixed point, limit-cycle, and strange attractors) they are also characterized by chaotic itinerancy and other novel concepts, such as open-ended evolution/dynamics, that describe their inherent behaviour.

Richness and the potential of the dynamical systems approach encourage us to go beyond merely adaptive behaviour. The higher functions, such as intention, motivation, emotion and consciousness, are within the scope of this study. Grey Walter has started the discussion of emotional, or play-like, behaviour by synthesizing artificial creatures . A wheeled vehicle containing a simple electric circuit can show unexpected and complex behaviour, comparable to that of living creatures. Without making real robots, Braitenberg made conceptual robots to discuss the higher functions . In his thought experiments, he designed vehicles using simple hard-wired electrical connections from sensory inputs to motor outputs. His vehicles gradually showed more complex cognitive behaviours by providing more complex internal structures. For example, the most primitive behaviour is a sense of “aggression”, which is simply given by attraction to a light source with a crossed sensory–motor connection. However, to simulate more complex behaviour, such as association and concept formation, he has to implement new wires, such as mnemotorix and ergotorix wires, with some Darwinian-type selections. Grey Walter and Braitenberg have one thing in common, in claiming that any apparently complex cognitive behaviour can be built up from simple sensory–motor coordination. That is, agents can be cognitive by having physical constraints. We basically agree that any meaningful cognition should be embodied, but focus on different aspects.

In this paper, we focus on the cognitive behaviours of turn-taking and imitation, caused by interactions between two or more humans, in which it is thought that the sharing of mental states and intentions with others is important. There are many ways to understand psychological phenomena by computer simulations and robot experiments rather than by studying human behaviour directly . We conducted computer simulations of two agents with internal dynamics, which were implemented by an artificial recurrent neural network, as a model of turn-taking behaviour. In our previous works, cognitive behaviours were explained from the dynamical systems perspective by coupling between agents with rich internal dynamics . Here, we generalize from turn-taking behaviour to autonomous role-changing, such as games of tag among children, and investigate the generic underlying mechanisms using the dynamical systems method. Therefore, this study focuses on different perspectives from those of fixed role-playing games (e.g., a pursuit-evasion game ). Here we take turn-taking as the simplest example that shows the diversity of dynamics. For turn-taking behaviour, it is necessary for roles to be exchanged autonomously, within a context constructed by the entities’ behaviours, e.g., chaser–evader and speaker–listener. When taking turns in a two-person conversation, people usually avoid overlapping or interrupting each other’s speech without setting some explicit cue to switch speakers. Some cues for this include eye contact and the detection of intonation changes. It is considered that turn-taking is established by coordination between predictions and the internal neural networks that compute the output from the inputs. Therefore, coupling between agents means a coupling of anticipatory systems with intrinsic dynamics.

By introducing neural architecture, evolutionary algorithm and a turn-taking game in §2 and 3, we explore four topics in the simulation. The first topic is dynamic repertoire. We describe how turn-taking is established with different forms of motion. In particular, we argue in §4.1 that regular motion behaviour evolves into chaotic behaviour. The second topic is predictability. Each agent has to predict the other’s future behaviour one step ahead. Interestingly, prediction precision decreases when the turn-taking role switches from one to the other. This will be discussed in §4.2. The third topic is ongoingness of interactions. Agents become robust against sensor noise; however, the turn-taking performance is established only when agents synchronize their dynamics precisely. This point is discussed in §4.3. The last topic is adaptability. As discussed in the section on dynamic repertoire, the turn-taking pattern appears to be different for different evolutionary generations. In section §4.4, we investigate the emergence of new spatio-temporal patterns by coupling agents from different generations. In §5, we discuss the potential linkage between these simulation results and the psychological experiments conducted by C. Trevarthen . A concept of intersubjectivity is also discussed.

## 2 The Model

We modelled the playing of a tag game in which the role of chaser, or evader, is not given to players in advance. There are some game models in which the roles are not predefined. Reynolds also showed that the abilities of chasing and evading evolve simultaneously by genetic programming in a game of tag, which is a symmetrical pursuit-evasion game . The variety in the behaviour of agents adapting to their environments is worth noting. In Reynolds’ game, switching between evader and chaser is predefined as happening when both agents come into physical contact. The difference between Reynolds’ model and ours is the spontaneous emergence of behaviour. Whether an agent plays the role of a chaser or an evader is dynamically determined in our model. On the other hand, Di Paolo modelled and studied social coordination with agents interacting acoustically . To avoid misperceiving the acoustical signals, their emission timings were entrained in an anti-phase state; the resulting behaviour resembles a turn-taking process.

There is a difference between Di Paolo’s turn-taking and ours. Both turn-taking behaviours are established by the coordination of agents through the history of their interactions. Di Paolo modelled turn-taking as the result of anti-phase signals to avoid signal interference; however, we modelled turn-taking behaviour as a result of coupling between richer internal dynamics. Therefore, in this paper, we pay more attention to the diversity of behaviour patterns.

### 2.1 Game and Environment

Here each agent has a circular body of radius R, with two diametrically opposed motors (Fig. 1). The motors can move the agent backwards and forwards in a two-dimensional unstructured and unlimited arena. The motion is described by the following equation of motion in terms of an agent’s heading angle ($`\theta `$) and its speed ($`v`$) in that direction.

$`M\dot{v}+D_1v+f_1+f_2=0,`$ (1)
$`I\ddot{\theta }+D_2\dot{\theta }+\tau (f_1,f_2)=0,`$ (2)
where $`f_1`$ and $`f_2`$ are the forward driving force, and $`\tau `$ denotes the torque. $`D_1`$ and $`D_2`$ express the resistance coefficients, and the agents have mass ($`M`$) and inertia ($`I`$). We solve the equations iteratively using the Runge–Kutta method. At each time step, the agents compute the forces from the inputs using the internal neural nets described below.

We assume there is no collision between agents because we focus on the internal states of the agents that generate turn-taking. Two agents try to coordinate their turn-taking behaviour, each trying to get behind the other. Because they cannot get behind each other simultaneously, the turn-taking cannot be achieved if both agents play chaser. Naturally, mutual turn-taking cannot be achieved if both agents play evader either. Therefore, it is necessary to have spontaneous symmetry break down so that one plays the role of chaser and the other plays the role of evader. However, mere symmetry breakdown is insufficient: temporal role changing is also required. By using recurrent neural networks, we focus on how the turn-taking dynamics are self-organized.

### 2.2 Agent Design

We designed the agents to have recurrent neural networks (Fig. 2). Inputs to an agent are the other agent’s position, distance and heading angle, relative to the agent. Agents move freely in the arena using two motors, the outputs of which are computed at every time step of the game. Each agent predicts the other’s next relative position, which is assigned to three output neurons. The dynamics of the recurrent neural network are expressed by the following equations at each time step $`t`$,

$`h_j(t)`$ $`=`$ $`g({\displaystyle \underset{i}{}}w_{ij}y_i(t)+{\displaystyle \underset{l}{}}w_{lj}^{}c_l(t1)+b_{j1}),`$ (3)
$`z_k(t)`$ $`=`$ $`g({\displaystyle \underset{j}{}}u_{jk}h_j(t)+b_{j2}),`$ (4)
$`c_l(t)`$ $`=`$ $`g({\displaystyle \underset{l}{}}u_{jl}^{}h_j(t)+b_{j3}),`$ (5)
$`g(x)`$ $`=`$ $`1/(1+\mathrm{exp}(x)),`$ (6)
where $`y_i,z_k,h_j`$ and $`c_l`$ represent input, output, hidden and context nodes, respectively. The respective number of nodes in these layers is set to $`(I,K,J,L)=(3,5,10,3)`$ throughout this paper. The symbols $`w_{ij},u_{jk},w_{lj}^{}`$ and $`u_{jl}^{}`$ denote the weights from input to hidden, hidden to output, context to hidden, and hidden to context neurons, respectively, and the parameter $`b`$ is a bias node. In this paper, we do not consider the results of predictions, which are discussed in . This network architecture evolves using a genetic algorithm as explained in the following section.

## 3 Genetic Algorithm and Noisy Environment

### 3.1 Genetic Algorithm

We update the weights according to the turn-taking performance. In practice, the weight set of the neural networks has a vector representation of the real weight values, which evolve using a genetic algorithm (GA).

We use a GA to evolve two separate populations, to avoid agents of a single genotype from dominating, in which case turn-taking is played between genetically similar agents. As a result, a player has to play against itself, which we wish to avoid. Each population contains $`P`$ individuals. The performance of all $`P^2`$ paired agents from the separated populations are evaluated at each generation. Agents that can exchange turns equally are evaluated as having greater fitness. At first, individuals in each population are initialized with random weight values. Then we calculate the fitness of each individual, based on its performance.

The highest value is given when both agents take their turn alternately and the agents can predict each other’s behaviour. A one-sided (i.e., role-fixed) behaviour is associated with lower fitness values. Practically, the fitness of an agent $`𝒂`$ from a population ($`A`$) against an agent $`𝒃`$ from the other population ($`B`$) is calculated as follows. Below, we define a total fitness $`F`$ as the sum of two fitnesses associated with prediction and turn-taking, respectively. When one agent gets behind the other, by definition the other agent has its turn and the rear scope is specified as $`RS`$, which is parameterized by two parameters $`r`$ and $`\varphi `$ (see Fig. 1). The agent in this state is said to be having its turn and is rewarded. A spatial position of agent $`𝒃`$ at time step t is represented by $`Pos_b(t)`$. This is compared with agent $`𝒂`$’s prediction value $`Pos_{ab}`$. Therefore the squared difference (Eq.(11)) is the measure of the precision of agent $`𝒂`$’s prediction.

$`F_a`$ $`=`$ $`{\displaystyle \frac{1}{P}}{\displaystyle \stackrel{P}{}}\left(s_1\times F_a^{turn}+s_2\times F_a^{predict}\right),`$ (7)
$`F_a^{turn}`$ $`=`$ $`{\displaystyle \underset{t}{\overset{T}{}}}g_a\left(t\right)\times {\displaystyle \underset{t}{\overset{T}{}}}g_b\left(t\right),`$ (8)
$`g_a(t)`$ $`=`$ $`\left\{\begin{array}{cc}1\hfill & Pos_a(t)RS_b(t)\hfill \\ 0\hfill & Pos_a(t)RS_b(t)\hfill \end{array}\right\},`$ (11)
$`F_a^{predict}`$ $`=`$ $`{\displaystyle \underset{t}{\overset{T}{}}}P_a\left(t\right)\times {\displaystyle \underset{t}{\overset{T}{}}}P_b\left(t\right),`$ (12)
$`P_a(t)`$ $`=`$ $`(Pos_b(t)Pos_{ab}(t))^2.`$ (13)
The performance of turn-taking is evaluated for different lengths of time ($`T=500,1,000`$ and $`1,500`$), so that agents cannot tell when the evaluation time is over. Evaluating the turn-taking performance at each GA generation, we leave the best $`E`$ individuals in each population and let them reproduce with specified mutation rates. The GA proceeds by repeating this procedure, and the recurrent neural networks evolve.
In addition, the following points should be noted.

### 3.2 Two Time Scales

Two time scales exist: the vehicle navigation time scale ($`\mathrm{\Delta }T_1`$), and the neural computation time scale ($`\mathrm{\Delta }T_2`$). The time evolution of the vehicle navigation is computed using the 4th order Runge-Kutta method, where $`\mathrm{\Delta }T_1`$ is set to 0.01. The basic process is that the neural net receives the sensor inputs and computes the motor outputs. By assuming that the vehicle navigation motion is faster than the internal neural time scale, we take $`100\mathrm{\Delta }T_1=\mathrm{\Delta }T_2`$. For simplicity, the neural net produces the outputs every 100 Runge–Kutta time steps. When the network structure evolves by GA, the time scale ratio is implicitly reflected in the net structure. Therefore, we believe that the same behaviour structure can be obtained, at least qualitatively, for a different scale ratio.

### 3.3 Noisy Environment

Living systems are involved in a fundamentally noisy environment. We know that our perception has to deal with noisy inputs. However, it is not possible to discriminate noise from other signals. We, as living systems, behave adaptively, cooperatively or selfishly while handling the problem. Therefore, we simulated the agents’ interacting with each other in a noisy environment. Noises are added to the input neurons at every game step during each run in the GA. The strength of noise is provided by uniform random numbers between zero and almost the maximum distance the agent can move during one game step. In the next sections, spatial patterns of turn-taking are studied as simulation results. If there is no excuse, those patterns are generated under a noise-free environment to clarify the intrinsic dynamics of the agents.

## 4 Simulation Results

Simulation was performed with a GA using 15 individuals ($`P=15,E=4`$). After several thousand GA generations, turn-taking is established between the two agents. The basic dynamics of the turn-taking was observed as follows. Two agents adjust their speeds and make turns automatically to switch from the role of evader to chaser and vice versa. In the following subsections, we investigate the turn-taking pattern realized from the dynamic repertoire, predictability, adaptability and evolvability concepts.

### 4.1 Diversity of Dynamic repertoire

First, the evolutionary algorithm effectively functions to improve the turn-taking performance. The development of the performance as a function of GA generations is depicted in Fig. 3. The resulting turn-taking patterns are sensitive to some of the settings. In particular, they are sensitive to the division of the agent population into two. In previous work, we encoded the pair of agents’ structures on the same gene . Then we encoded them separately but used a single population. That algorithm can also develop turn-taking behaviour but with much less diversity than the present algorithm. When the agents are on the same gene, it is difficult to show diversity as their net structures are too correlated. With a single population, development of an agent that can take turns with itself (its relatives) is enhanced. Therefore, there is a strong probability that the dynamics of the turn-taking may be tuned for self-turn-taking. To avoid this situation, we used the two-population structure.

Figure 4 shows examples of the spatial trails of an agent from different GA generations with different initial population structures. For the sake of clarity, a single agent’s trail is depicted. A paired agent tends to show the same trail with different phases.

We can classify these trail patterns approximately into regular, chaotic and others based on their appearance in space and time. When spatial trails consist of regular curves, and turns are exchanged almost periodically (which corresponds to a turning point on the figures), we call them regular turn-taking. On the other hand, if spatial trails have irregular curves with non-periodic turn-taking, we call them chaotic turn-taking. The remaining unclassified patterns are discussed below.

In the earlier GA generations, agents with regular turn-taking evolve to yield higher performance (Fig. 4(a) and (b)). The behaviour structure is as follows. One agent follows the other and passes it; then it slows as does the other agent; then both agents simultaneously turn around quickly. This returns the agents to the first phase. A series of behaviour patterns repeats almost periodically and the envelope curve of these trails constitutes a circle by fixing the centre location. In the later GA generations, more chaotic patterns emerge (Fig. 4 (c) to (h)). In contrast to the regular patterns, the turns are exchanged in different places with irregular time intervals. Therefore, the spatio-temporal pattern becomes chaotic and agents move around the entire space.

The evolution of turn-taking type from regular to chaotic is explained as follows. The evolutionary pressure of GA at first allows the agents to move stably in the noisy environment. A structured turn-taking behaviour can only be built up on stable motion dynamics that are insensitive to random noise. As argued briefly in the introduction, noise and intentional action is difficult to distinguish when the agents’ motions become chaotic. However, when their actions appear regular, we can interpret that the agents can more easily distinguish noise from the other agent’s intentional motion as they show different performance with and without partners’ adaptive motions (see §4.3). Therefore, the regular type emerges earlier than the chaotic motion. As shown in Fig. 5, regular turn-taking occurs at almost the same spatial location with different noise series. However, the chaotic type is sensitive to the noise series. The total performance of turn-taking remains high in both cases.

That is, regular turn-taking pattern suppresses a variety of dynamic repertoires. By doing so, it becomes robust against sensory noise. On the other hand, chaotic turn-taking pattern has the potential to develop dynamic repertoire, and therefore it becomes more adaptive, which is studied in §4.4.

Intuitively, agents who can take turns in the presence of noise can take turns perfectly without noise. However, this does not hold for some agents found in later GA generations. As shown in Fig. 6, agents can only take turns when there is sensory noise. We call this phenomenon Noise-induced turn-taking. As shown in the figure, there is a strong attractor to a circular motion without exchanging turns. The two agents have different neural structures, and the resulting turn-taking behaviour is generally asymmetrical. Without noise, one agent is never able to take its turn. In addition, it forms an attractor in the sense that adding a small noise cannot break this one-sided behaviour. True turn-taking only emerges above a certain noise level (Fig. 7). In another case, there exist three attractors when there is no sensory noise. One is that agent A chases the rear of agent B closely. Another is the opposite, and the last is that in which both agents chase each other. Each of the three attractors consist of circular orbits. The transition between attractors is caused by noise. Without noise, agents are trapped by one of the attractors.

Compared with these noise-induced behaviours, chaotic turn-takers can spontaneously establish turn-taking behaviour without noise. Even if noise is introduced into the system, chaotic turn-takers can establish turn-taking behaviours independent of the low noise level. That is, they do not utilize noise but suppress the effect of noise to perform turn-taking. Conversely, noise-induced turn-takers need noise to perform turn-taking.

### 4.2 Prediction Capability and Role Switching

These observations were analysed in terms of prediction capability of agents. The agents, after thousands of GA generations, are able to predict their partner’s future movements while turn-taking. Three outputs of the recurrent network simulate the other agent’s future location and heading from the current input. Fig. 8, shows the precision of predictions and the associated turn-taking patterns. In earlier GA generations, one agent’s prediction is far better than the other’s. In later generations, both predictions are improved. However, through entire GA generations, the predictions almost periodically break down when their turns (roles) are exchanged. As indicated in the figure, the prediction is also perturbed by noisy inputs. However, the effect is much smaller than that of the other agent’s action.

It should be noted that these prediction outputs are not designed explicitly to do anything in generating action sequences. However, because they depend on the common context neurons that also control the motion patterns, simulating each other’s behaviour and generating the motor outputs have indirect correlations. The correlation between prediction breakdown and the turn-taking performance will be reported elsewhere.

### 4.3 Ongoingness of Interactions

The inherent adaptability of each turn-taking pattern can be studied using its stability in the presence of noise. In other words, we study an agent’s ability to discriminate between noise and the adaptive behaviour of the other agent. In this section, we compare the behaviour of “live interaction” with “recorded interaction”. The “live interaction” is normal interaction between evolved agents, and the “recorded interaction” is that between an agent and a virtual agent, defined below.

First, we selected the two best agents, A and B, from each population. Turn-taking between these agents was studied without introduced noise. This is what we term “live interaction”. The trails of the agents were recorded during the run. Then, turn-taking between agent A and the recorded trail of agent B (i.e., a virtual agent) was conducted. This is what we term “recorded interaction”. We perturb the recorded trail and simulate the changes in the turn-taking dynamics.

Figure 9 (a) shows the growth of a discrepancy between A-virtual B and A-perturbed virtual B (chaotic turn-takers). During the initial few hundred steps, no discrepancy was observed. The behaviours are similar as shown in the figure. However, a small noise was amplified and the orbit drastically changed from the original orbit at approximately 800 time steps. In terms of the turn-taking behaviours, the adaptive agent can no longer recover harmonization with the perturbed virtual agent. The agent approaches the trail and tries dynamically to resume the original turn-taking behaviour.

Another example (the agents at 3,000 generations) is shown in Fig. 9 (b). These agents established regular turn-taking. In this case, the agents could cope with the perturbed virtual agent. Note that agents that have constructed regular turn-taking behaviour do not always, but frequently do, have a tendency to cope with a perturbed virtual agent, although this varies with the timing and strength of the perturbation. Sometimes turn-taking behaviour breaks down when more noise is added to the recorded trail. However, there are some examples in which turn-taking recovers after a period of discrepancy.

### 4.4 Evolution of Adaptability

Another novel feature of adaptability was examined. We show here that adaptability can generate novel dynamics by constructing new couplings. We examine the behaviours of new couplings between two agents from different GA generations as follows. After the turn-taking performance had attained a satisfactory plateau, we selected two individuals from different generations to play. This was to examine how they performed turn-taking without having the common experience of co-evolution. Taking agents from generations 10,000 and 3,000 as examples, we evaluated the performances of the new pairs for each generation (Fig. 10). In fact, the novel pairs often failed to sustain the same performance as the original pairs. However, the synthesized dynamics often showed novel structures. The examples can be found in Fig.11. Agents that perform chaotic turn-taking after 10,000, 8,000, and 7,000 generations (Fig. 11 (a),(c) and (e)) are coupled with agents from each different generation. As is seen in the figure, the newly coupled agents also show chaotic turn-taking but with a different kind of motion (d). Coupling of generations 1,000–7,000 and 8,000–7,000 shows a similar pattern to that by the agents from generation 7,000, which is shown in (b) and (f).

In summary, (i) Novel structures sometimes inherit the original pattern of one of the agents but not always, (ii) Agents that readily exhibit chaotic turn-taking pattern lose the original pattern and adapt to the other agent’s pattern, and (iii) conversely, regular turn-takers simply retain their original pattern and show little adaptability to a new partner.

The last point is clearly shown in Fig. 12. The regular turn-takers can only achieve higher performance with agents from near generations (Fig. 10) Our hypothesis is that chaotic turn-takers are more adaptive than regular ones. The observation here confirms the hypothesis, but we should note that performance sometimes differs significantly between populations A and B from the same GA generation. Figure 13 illustrates how turn-taking performance varies from generation to generation. We deduce from this figure that they are basically symmetrical for populations A and B. Sometimes there are notable exceptions—e.g., population A from generation 8,000–10,000 compared with population B from generation 10,000–12,000. It should also be noted that genetically closer agents can collaborate better than more distantly related agents. However, qualitatively, beyond generation 6,000, agents become more adaptive than those of earlier generations.

The turn-taking pattern resulted from the collaboration of two agents. Therefore, a neural structure in the body of a single agent alone cannot explain the turn-taking dynamics. This is an interesting part of the present study, but at the same time, a gap between the two agents may develop. That is, when one population becomes very adaptive against many others, it is not necessary for the other population to become very adaptive; it may simply become a test data set for the former population to become “universal” turn-takers. As far as we know, such universal turn-takers are yet to evolve. Here we notice that chaotic turn-taker is better at eliciting coordinated behaior from the partner.

It is also worth noting that the “experience” of two agents interacting with each other is a prerequisite for better turn-taking. The history, or the experience, of how agents have collaborated to perform turn-taking determines with whom an agent can take turns. In new pairs, responses of one agent to the other often occur at the wrong time, whereas the original pairs show complete synchronization of turn-taking. That is, we insist that it is not the neural structure but the collaboration of timing and patterning that is responsible for the better turn-taking behaviour. This is true not only for this special type of interaction—i.e., turn-taking—but may be true for cognitive interaction in general. We will argue this point in the final section of this paper.

## 5 Discussion

It was found in the virtual agents experiment (§4.3), that chaotic turn-takers are much more sensitive to the difference between live and recorded inputs. Their turn-taking patterns are driven by the ongoing interaction. On the other hand, regular turn-takers are relatively insensitive to the difference. As evolution continues, chaotic turn-taking replaces regular turn-taking in the GA simulations. This may be due to regular turn-takers’ being less adaptive than chaotic turn-takers in the sense that they can only cope with fewer agents. This is clearly seen in the new coupling experiment §4.4. The coupling with regular turn-takers only generates circular patterns but chaotic turn-takers show various patterns. In summary, we claim that chaotic turn-taking is less robust in the presence of noise but has more adaptability, compared with regular turn-taking.

This complementary relationship between adaptability and robustness has some implications in some empirical experiments. Let us introduce Trevarthen’s double-monitor experiments between a baby-infant and its mother , and Nadel’s mutual imitation experiments . In Trevarthen’s experiment, mother and baby-infant only communicate through videos that display their faces to each other. For the baby-infant to engage with the mother, correct style and timing are required. If the recorded video of the mother is displayed to the baby-infant, the baby-infant becomes withdrawn and depressed. Nadel studied how the mutual imitation game progresses between children and discussed a non-affordant means of using objects to trigger the interaction. Children regularly switch between the roles of imitating and being imitated, by having new imitation patterns.

Trevarthen’s experiments show that it is not necessarily important for the baby-infant that the mother be displayed on the monitor. It can be assumed that the most important clue during interactions is the ongoing anticipation of a partner. The baby-infant performs some actions and anticipates the mother’s reactions reflecting the baby-infant’s actions, and this is also true with respect to the mother’s anticipation of the baby-infant. Interactions in social behaviour, including turn-taking, can be established when these anticipations are mutually formed dynamically. Furthermore, it is shown by Nadel’s experiment that an affordant way of using objects can maintain interaction—i.e. some form of novelty/unpredictability is required. In our simulations, when an agent calculates outputs, this calculation simultaneously affects the internal dynamics. That is, the actions performed form its internal dynamics as much as actions form anticipations in the statement above. The agent receives a partner’s actions as inputs that reflect the agent’s own actions. We maintain that turn-taking is established when these structures are mutually organized. Turn-taking is therefore broken in the simulation with virtual agents. However, our simulations also show that unpredictability is found when turn-taking occurs. We therefore claim that mutually adaptive coupling of actions and internal dynamics between agents is essential for the establishment of cognitive interaction, which may be related to intersubjectivity.

Acknowledgements: We thank Ryoko Uno and Gentaro Morimoto for their fruitful discussions. This work is partially supported by Grant-in aid (No. 09640454 and No. 13-10950) and also by a grant-in-aid from The 21st Century COE (Center of Excellence) program (Research Center for Integrated Science) of the Ministry of Education, Culture, Sports, Science, and Technology, Japan.