text:
# Revealing Superfluid–Mott-Insulator Transition in an Optical Lattice

## Abstract

We study (by an exact numerical scheme) the single-particle density matrix of $`10^3`$ ultracold atoms in an optical lattice with a parabolic confining potential. Our simulation is directly relevant to the interpretation and further development of the recent pioneering experiment . In particular, we show that restructuring of the spatial distribution of the superfluid component when a domain of Mott-insulator phase appears in the system, results in a fine structure of the particle momentum distribution. This feature may be used to locate the point of the superfluid–Mott-insulator transition.

The fascinating physics of the superfluid–insulator transition in a system of interacting bosons on a lattice has been attracting constant interest of theorists during recent years . Lattice bosons is one of the simplest many-body problems with strong competition between the potential and kinetic energy, and a typical example of the quantum phase transition system. One of its great advantages is the possibility to study it by powerful Monte Carlo methods which nowadays allow simulations of many thousands of particles at low temperature with unprecedented accuracy (see, e.g., ). However, until very recently the canonical Bose-Hubbard model

$$H=t\underset{<ij>}{}a_i^{}a_j+\frac{U}{2}\underset{i}{}n_i^2\underset{i}{}\mu _in_i,$$
(1)
(where $`a_i^{}`$ creates a particle on the site $`i`$, $`<ij>`$ stands for the nearest-neighbor sites, $`n_i=a_i^{}a_i`$, and $`t`$, $`U`$, and $`\mu _i`$, are the hopping amplitude, the on-site interaction, and the on-site external field, respectively) was not particularly useful in the analysis of realistic systems. The situation has changed with the exciting success of the experiment by Greiner et al.  (originally proposed by Jaksch et al. ) in which a gas of ultracold 87Rb atoms was trapped in a three-dimensional, simple-cubic optical lattice potential. The uniqueness of the new system is that it is adequately described by the Bose-Hubbard Hamiltonian , and allows virtually an unlimited control over the strength of the effective interparticle interaction $`U/t`$ and particle density.

The characteristic feature of the experimental setup of Ref.  is the presence of the overall parabolic potential $`V(r)`$ which confines the sample. This feature could be of great advantage if one was able to directly measure the spatial density distribution in the trap. We recall the structure of the $`\mu U/t`$ phase diagram for the Bose-Hubbard system  which predicts commensurate particle density distribution for the insulating phase whenever the chemical potential lies within the Mott-Hubbard gap. The slowly varying (at the length scale of the lattice period) trapping potential effectively provides a scan over $`\mu `$ of this phase diagram at a fixed value of $`U/t`$.

Unfortunately, what is measured in the experiment is not the original spatial density distribution in the trap, but the absorption image of the free evolving atomic cloud, after the trapping/optical potential is removed; i.e. the quantity which is directly related to the the single-particle density matrix in momentum space, $`\rho _{\mathrm{𝐤𝐤}}=n_𝐤`$. \[This statement implies that in the free evolving atomic cloud the interparticle interaction can be neglected; see the discussion below.\] Now, in terms of $`n_𝐤`$ the inhomogeneous trapping potential is a disadvantage since it broadens the superfluid $`\delta `$-functional contribution at $`𝐤=0`$, and the observed picture is a convolution of the original real-space density matrix $`\rho (𝐫,𝐫^{})`$. As we show below one has then to look at the fine-structure of the central peak in the experimental data to decipher the Mott-Hubbard phase diagram (the “fading” of the Bragg peaks in the experiment has very little to do with the superfluid–Mott-insulator transitions and happens when the system is already deep in the insulating phase).

In this Letter we relate quantitatively the particle distribution in momentum space observed in experiments to the corresponding spatial density distribution in the trap. Our ultimate goal is to reveal which features (if any) in the structure of $`n_𝐤`$ indicate unambiguously the presence of the Mott phase. To this end we perform quantum Monte Carlo simulations of the single-particle density matrix for the Bose-Hubbard systems with up to $`16^3`$ lattice sites using continuous-time Worm algorithm . We find that the onset of the phase transition in the trap center should result in appearance of at least one satellite peak in $`n_𝐤`$, reflecting a shell-type form of the superfluid component. This peak was not mentioned in the experiments of Ref. . We suggest a possible explanation to this fact, and argue that by collimating the expanding atomic cloud one can render this peak observable. We also discuss the role of self-repulsion in the free expanding cloud, which can affect the simple interpretation of the absorption images in terms of the initial single-particle momentum distribution.

In Figs. 1(a-f) we present our data for the density distribution as function of the lattice cite distance from the trap center, $`r/a`$, where $`a`$ is the lattice constant. For all practical purposes one may assume the zero-temperature limit here. The simulation was done at finite but very low $`T=0.2t`$; the relevant energy parameters in this model are the bandwidth $`W=12t`$ and $`(U/t)_c35`$ . In accordance with the phase diagram of Ref. , we observe a shell-type structure of the particle density with the Mott-insulator phases visible as integer plateau regions.

Next, we relate each of the above figures to the corresponding momentum distribution function $`n_𝐤`$. By definition, $`n_𝐤=\mathrm{d}^3r\mathrm{d}^3r^{}\mathrm{exp}[i𝐤(𝐫𝐫^{})]\rho (𝐫,𝐫^{})`$, where $`\rho (𝐫,𝐫^{})=\psi ^{}(𝐫)\psi (𝐫^{})`$, and $`\psi (𝐫)`$ is the bosonic field operator. In our case of a single-zone lattice, the field operator is expanded as follows:

$$\psi (𝐫)=\underset{i}{}\varphi (𝐫𝐫_i)a_i,$$
(2)
where $`\varphi `$ is the Wannier function. We thus finally have

$$n_𝐤=|\varphi (𝐤)|^2\underset{i,j}{}e^{i𝐤(𝐫_i𝐫_j)}\rho _{ij},$$
(3)
where

$$\rho _{ij}=a_i^{}a_j,$$
(4)
and $`\varphi (𝐤)`$ is the Fourier transform of $`\varphi (𝐫)`$. From Eq. (3) it is seen that up to a trivial reweighting factor $`|\varphi (𝐤)|^2`$ the distribution is a periodic function in the reciprocal lattice. Thus without loss of generality we may restrict ourselves to the first Brillouin zone. Actually, $`\varphi (𝐤)`$ has nothing to do with the Bose-Hubbard model, being a non-universal property of the lattice cite potential; in what follows we will ignore this function altogether by formally setting it to unity.

Having calculated $`\rho _{ij}`$ with the Worm algorithm , we readily obtain $`n_𝐤`$ using Eq. (3); the results are presented in Figs. 2(a-f). In Fig. 2(a) we see a typical picture for the strongly correlated superfluid phase, characterized by a single, narrow peak at small momenta. When a domain of the Mott-insulating phase appears in the center of the trap (where the on-site filling is close to unity), a pronounced fine structure develops in Fig. 2(b). We associate this structure with the shell-type form of the condensate wave-function. To prove the point we model the situation with the pure-condensate density matrix $`\rho (𝐫,𝐫^{})=\mathrm{\Psi }_0^{}(𝐫)\mathrm{\Psi }_0(𝐫)`$, where the condensate wave-function $`\mathrm{\Psi }_0(𝐫)`$ has the shell-type form with the shell radius $`l`$. The presence of the Mott insulator is taken into account through the suppression of the $`\mathrm{\Psi }_0(𝐫)`$ in the center. The Fourier transform of such $`\mathrm{\Psi }_0(𝐫)`$ is alternating in sign, with the half-period in $`k`$ related to the shell radius as $`k\pi /l`$. Thus in the pure condensate we would see exact zeros in $`n_𝐤`$ with the typical separation between them $`\pi /l`$. Surprisingly, this naive model works extremely well and adequately describes the case of the realistic strongly correlated system close to the phase transition (in Fig. 2(b) the coupling $`U/t=32`$ is close to the critical value estimated in Refs. ). We consider the appearance of the satellite peaks as a clear signature of the Mott-insulator transition in the center of the trap.

In Fig. 2(c) the coupling strength is significantly increased to $`U/t=80`$, but the shell-type structure of the superfluid phase, and the corresponding fine-structure in $`n_𝐤`$, is still present. The crucial difference between Fig. 2(c) and Fig. 2(d) is in the suppression (almost complete) of the superfluid fraction. Now the distribution $`n_𝐤`$ has only the central peak and is extended towards large momenta, as expected for the system deep in the insulating phase. Still, it is not flat (!), which tells us about large off-diagonal correlations between the nearest lattice cites even for $`U/t`$ as large as $`2.5(U/t)_c`$.

Fig. 2(e) is similar in physics to the case (a), except for the large-momentum tail due to the Mott-insulator shell. Finally, in Fig. 2(f) we again see the fine structure of satellite peaks reflecting the appearance of the Mott-insulator phase in the center of the trap and the corresponding superfluid phase shell \[in close resemblance with Fig. 2(b)\].

We note, that momentum distributions $`n_𝐤`$ presented above, may be observed experimentally if atoms are collimated out of the expanding cloud so that the distribution in a given direction is photographed. In the current setup, the absorption images of the three dimensional distribution are taken along two orthogonal axis. This procedure reveals only the integral $`N(k_x,k_y)_{\mathrm{}}^{\mathrm{}}𝑑k_zn(𝐤)`$. It is easy to see that integration effectively erases fine-structure features of $`n(𝐤)`$—although peaks do not disappear completely, they now show up as shoulders in $`N(k_x,k_y)`$.

Finally, we would like to discuss the problem of the repulsive interaction between particles during the initial period of their free expansion. Obviously, the interpretation of the final absorption pattern in terms of the initial momentum distribution is valid only if the effect of this interaction is small enough. Meanwhile, given realistic experimental parameters, this turns out to be the case only for rather moderate system sizes. The criterion for neglecting the effect of interparticle repulsion is

$$E_{\mathrm{kin}}/E_{\mathrm{pot}}1,$$
(5)
where $`E_{\mathrm{kin}}`$ and $`E_{\mathrm{pot}}`$ are, respectively, the kinetic and potential energy per particle in the most fragile low-momentum part of the distribution $`n_𝐤`$ and at the most dangerous period of free evolution at the end of the restructuring period, when the “discrete” distribution of density transforms into the “continuous” spatial distribution with the typical size of order of the original system size (plus the corresponding replicas in higher Brillouin zones). For the potential energy we have $`E_{\mathrm{pot}}nU_0`$, where $`n`$ is the continuous number density, $`U_0=4\pi \mathrm{}^2a_s/m`$, $`a_s`$ is the $`s`$-scattering length, and $`m`$ is the atom mass. Recalling that the lattice filling factor is of order unity, we can estimate $`n1/a^3`$. The lowest kinetic energy is associated with the spatial distribution of the condensate. Estimating $`E_{\mathrm{kin}}\pi ^2\mathrm{}^2/ma^2L^2`$, where integer $`L`$ stands for the typical size of the superfluid component in units of lattice constants, we arrive at a simple requirement

$$(a/a_s)L^21.$$
(6)
In the experiment of Ref. , the ratio $`a/a_s`$ is of order $`10^2`$. Hence, we are restricted to $`L<10`$, that is to typical system sizes of our present simulation. Note also, that the condition (5) is much easier to satisfy for the atomic cloud in the second Brillouin-zone peak, where the spatial density is significantly suppressed .

Summarizing, we presented a simulation of the ground state properties of ultracold atoms in an optical lattice with confining external potential, in the regime where the Mott-insulator and superfluid phases coexist. We have demonstrated that when the insulator domain in the center of the trap is surrounded by the superfluid component, the global momentum distribution of particles features satellite peaks. This picture can be employed by the experiment as an unambiguous evidence of the Mott transition. We do not see other features of the momentum distribution that could be associated with the Mott transition: unless the ratio $`U/t`$ is not much larger than the critical one, the momentum distribution in the reciprocal lattice still has a peaked form reflecting strong local off-diagonal correlations.

This work was supported by the National Science Foundation under Grant DMR-0071767. BVS acknowledges a support from the Russian Foundation for Basic Research (RFBR) under Grant 01-02-16508 and from the Netherlands Organization for Scientific Research (NWO). VAK acknowledges a support from RFBR under Grant 00-02-17803. We acknowledge helpful discussions with H. Monien and Yu. Kagan.


format: text

text:
# Abundance gradients and the role of SNe in M87

## 1 Introduction

The X-ray emitting hot intra-cluster medium (ICM) of clusters of galaxies is known to contain a large amount of metals: for rich clusters between red-shift 0.3-0.4 and the present day the observed metallicity is about 1/3 the solar value (Mushotzky & Lowenstein, 1997; Fukazawa et al., 1998; Allen & Fabian, 1998; Della Ceca et al., 2000; Ettori et al., 2001), suggesting that a significant fraction of the ICM has been processed into stars already at intermediate red-shifts.

While the origin of the metals observed in the ICM is clear (they are produced by supernovae), less clear is the transfer mechanism of these metals to the ICM. The main mechanisms that have been proposed for the metal enrichment in clusters are: enrichment of gas during the formation of the proto-cluster (Kauffman & Charlot, 1998); ram pressure stripping of metal enriched gas from cluster galaxies (Gunn & Gott, 1972; Toniazzo & Schindler, 2001); stellar winds AGN- or SN-induced in Early-type galaxies (Matteucci & Vettolani, 1988; Renzini, 1997).

Spatially resolved abundance measurement in galaxy clusters are of great importance because they can be used to measure the precise amounts of metals in the ICM and to constrain the origin of metals both spatially and in terms of the different contributions of the two different type of SNe (SNII and SNIa) as a function of the position in the cluster. The first two satellites able to perform spatially resolved spectroscopy, *ASCA* and *BeppoSAX* have revealed abundance gradients in cD clusters (Dupke & White, 2000b; De Grandi & Molendi, 2001), in particular M87/Virgo (Matsumoto et al., 1996; Guainazzi & Molendi, 2000), and variations in Si/Fe within a cluster (Finoguenov et al., 2000) and among clusters (Fukazawa et al., 1998). Since the SNIa products are iron enriched, while the SNII products are rich in $`\alpha `$ elements, such as O, Ne, Mg and Si, the variations in Si/Fe suggest that the metals in the ICM have been produced by a mix of the two types of SNe. The exact amount of these mix still remain controversial: Mushotzky et al. (1996) and Mushotzky & Lowenstein (1997) showed a dominance of SNII ejecta, while other works on *ASCA* data (Ishimaru & Arimoto, 1997; Fukazawa et al., 1998; Finoguenov et al., 2000; Dupke & White, 2000a), still indicating a predominance of SNII enrichment at large radii in clusters, do not exclude that as much as 50% of the iron in clusters come from SNIa ejecta in the inner part of clusters.

M87 is the cD galaxy of the nearest cluster and its high flux and close location allows, using the unprecedented combination of spectral and spatial resolution and high throughput of the EPIC experiment on board *XMM-Newton*, a detailed study of the ICM abundance down to the scale of the kpc. Throughout this paper, we assume $`H_0=\mathrm{\hspace{0.17em}50}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, $`q_0=0.5`$ and at the distance of M87 $`1^{}`$ corresponds to 5 kpc.

Recently, when we were finishing the writing of the paper, we have learned of the Finoguenov et al (2001) analysis of the same data. Our work is complementary with theirs in the sense that we use 10 annular bins instead of 2, fully exploiting the quality of the *XMM* data and taking into account the multi temperature appearance of spectra in the inner regions with adequate spectral modeling.

The outline of the paper is as follows. In section 2 we give information about the *XMM* observation and on the data preparation. In section 3 we suggest the use of a unique set of “standard” abundances. In section 4 we describe our spectral modeling. In section 5 we present spatially resolved measurements of metal abundances and abundance ratios and in section 6 we discuss our results. A summary of our conclusions is given in section 7.

## 2 Observation and Data Preparation

M87/Virgo was observed with *XMM-Newton* (Jansen et al., 2001) during the PV phase with the MOS detector in Full Frame Mode for an effective exposure time of about 39 ks. Details on the observation have been published in Böhringer et al. (2001) and Belsole et al. (2001). We have obtained calibrated event files for the MOS1 and MOS2 cameras with SASv5.0. Data were manually screened to remove any remaining bright pixels or hot column. Periods in which the background is increased by soft proton flares have been excluded using an intensity filter: we rejected all events accumulated when the count rates exceeds 15 cts/100s in the \[10 – 12\] keV band for the two MOS cameras.

We have accumulated spectra in 10 concentric annular regions centered on the emission peak extending our analysis out to 14 arcmin from the emission peak, thus exploiting the entire *XMM* field of view. We have removed point sources and the substructures which are clearly visible from the X-ray image (Belsole et al., 2001) except in the innermost region, where we have kept the nucleus and knot A, because on angular scales so small it is not possible to exclude completely their emission. We prefer to fit the spectrum of this region with a model which includes a power law component to fit the two point like sources. We include only one power law component due to the similarity of the two sources spectra (Böhringer et al., 2001). The bounding radii are 0-0.5, 0.5-1, 1-2, 2-3, 3-4, 4-5, 5-7, 7-9, 9-11 and 11-14. The analysis of the 4 central regions within 3 arcmin was already discussed in Molendi & Gastaldello (2001).

Spectra have been accumulated for MOS1 and MOS2 independently. The Lockman Hole observations have been used for the background. Background spectra have been accumulated from the same detector regions as the source spectra.

The vignetting correction has been applied to the spectra rather than to the effective area, as is customary in the analysis of EPIC data (Arnaud et al., 2001). Spectral fits were performed in the 0.5-4.0 keV band. Data below 0.5 keV were excluded to avoid residual calibration problems in the MOS response matrices at soft energies. Data above 4 keV were excluded because of substantial contamination of the spectra by hotter gas emitting further out in the cluster, on the same line of sight.

As discussed in Molendi (2001) there are cross-calibration uncertainties between the spectral response of the two EPIC instruments, MOS and PN. In particular for what concern the soft energy band (0.5-1.0 keV) fitting six extra-galactic spectra for which no excess absorption is expected, MOS recovered the $`\mathrm{N}_\mathrm{H}`$ galactic value, while PN gives smaller $`\mathrm{N}_\mathrm{H}`$ by $`12\times 10^{20}\mathrm{cm}^2`$. Thus we think that at the moment the MOS results are more reliable than the PN ones in this energy band, which is crucial for the O abundance measure. For this reason and for the better spectral resolutions of MOS, which is again important in deriving the O abundance, we limit our analysis to MOS data.

## 3 Solar Abundances

The elemental abundances of astrophysical objects are usually expressed by the relative values to the solar abundances. The so-called solar abundances can be either “meteoritic” or “photospheric”.

This distinction between “meteoritic” and “photospheric” solar abundances was made in the review by Anders & Grevesse (1989). Significant discrepancies exist between the two sets of abundances quoted in that paper, particularly for iron and this has caused in the past some controversy in the discussion of the results of cluster abundances (Ishimaru & Arimoto, 1997; Gibson et al., 1997). However recent photospheric models of the sun indicate that photospheric and meteoritic abundances agree perfectly and the community has converged toward a “standard solar composition” (Grevesse & Sauval, 1998), with suggestions to the astrophysical community to accept this new state of the art (Brighenti & Mathews, 1999). For the above reason, in this paper we shall adopt the Grevesse & Sauval (1998) values. Since the solar abundance table used by default in XSPEC is based on photospheric values of Anders & Grevesse (1989), we have switched to a table taken from the data by Grevesse & Sauval (1998) by means of the XSPEC command ABUND. In general a simple scaling allows to switch from one set of abundances to the other.

## 4 Spectral modeling and plasma codes

All spectral fitting has been performed using version 11.0.1 of the XSPEC package.

All models discussed below include a multiplicative component to account for the galactic absorption on the line of sight of M87. The column density is always fixed at a value of $`1.8\times 10^{20}\mathrm{cm}^2`$, which is derived from 21cm measurements (Lieu et al., 1996). Leaving $`\mathrm{N}_\mathrm{H}`$ to freely vary does not improve the fit and does not affect the measure of the oxygen abundance, which could have been the more sensitive to the presence of excess absorption. The $`\mathrm{N}_\mathrm{H}`$ value obtained is consistent within the errors with the 21cm value.

The temperature profile for M87 (Böhringer et al., 2001) shows a small gradient for radii larger than $`2`$ arcmin and a rapid decrease for smaller radii. Moreover, as pointed out in Molendi & Pizzolato (2001) all spectra at radii larger than 2 arcmin are characterized by being substantially isothermal (although the spectra of the regions between 2 and 7 arcmin are multi temperature spectra with a narrow temperature range rather than single temperature spectra), while at radii smaller than 2 arcmin we need models which can reproduce the broad temperature distribution of the inner regions.

We therefore apply to the central regions (inside 3 arcmin) three different spectral models.

A two temperature model (vmekal + vmekal in XSPEC and model II in Molendi & Gastaldello (2001) using the plasma code MEKAL (Mewe et al., 1985; Liedahl et al., 1995). This model has 15 free parameters: the temperature and the normalization of the two components and the abundance of O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Fe and Ni, all expressed in solar units. The metal abundance of each element of the second thermal component is bound to be equal to the same parameter of the first thermal component. This model is used (e.g. Makishima et al. (2001) and refs. therein) as an alternative to cooling-flow models in fitting the central regions of galaxy clusters.

A “fake multi-phase” model (vmekal + vmcflow in XSPEC and model III in Molendi & Gastaldello (2001). This model has 15 free parameters, as the two temperature model, because the maximum temperature $`T_{max}`$ is tied to the vmekal component temperature. As indicated in recent papers (Molendi & Pizzolato, 2001; Molendi & Gastaldello, 2001) this model is used to describe a scenario different from a multi-phase gas, for which it was written for: the gas is all at one temperature and the multi-phase appearance of the spectrum comes from projection of emission from many different physical radii. A more correct description will be given by a real deprojection of the spectrum (Pizzolato et al., in preparation).

The third model is the analogue of the vmekal two temperature model using the plasma code APEC (Smith et al., 2001). This model has 14 free parameters, one less than the corresponding model using vmekal because APEC misses the Na parameter.

We can’t adopt an APEC analogue of the fake multi-phase model because the cooling flow model calculating its emission using APEC is still under development. Given the substantial agreement between 2T and fake multi-phase model (Molendi & Gastaldello, 2001), we can regard the 2T APEC results as indicative also for a fake multi-phase model.

For the spectrum accumulated in the innermost region we included also a power law component to model the emission of the nucleus and of knot A.

For the outer regions (from 3 arcmin outwards) we apply single temperature models: vmekal using the MEKAL code, with 13 free parameters and vapec using the APEC code, with 12 free parameters.

As pointed out by the authors of the new code, cross-checking is very important, since each plasma emission code requires choosing from a large overlapping but incomplete set of atomic data and the results obtained by using independent models allows critical comparison and evaluation of errors in the code and in the atomic database.

## 5 Results

### 5.1 Abundance measurements and modeling concerns

The X-ray emission in cluster of galaxies originates from the hot gas permeating the cluster potential well. The continuum emission is dominated by thermal bremsstrahlung, which is proportional to the square of the gas density times the cooling function. From the shape and the normalization of the spectrum we derive the gas temperature and density. In addition the X-ray spectra of clusters of galaxies are rich in emission lines due to K-shell transitions from O, Ne, Mg, Si, S, Ar and Ca and K- and L-shell transitions from Fe and Ni, from which we can measure the relative abundance of a given element.

In Figure 1 we show the data of the 3-4 bin together with the best fitting model calculated using the MEKAL code. The model has been plotted nine times, each time all element abundances, except one, are set to zero. In this way the contribution of the various elements to the observed lines and line blends become apparent. In the energy band (0.5-4 keV) we have adopted for the spectral fitting, the abundance measurements based on K-lines for all the elements except for Fe and Ni, for which the measure is based on L-lines. The K-lines of O, Si, S, Ar and Ca are well isolated from other emission features and clearly separated from the continuum emission, which are the requirements for a robust measure of the equivalent width of the lines and consequently of the abundances of these elements. The Fe-L lines are known to be problematic, because the atomic physics involved is more complicated than K-shell transitions (Liedahl et al., 1995), but from the very good signal of *XMM* spectra and from the experience of *ASCA* data (Mushotzky et al., 1996; Hwang et al., 1997; Fukazawa et al., 1998) we can conclude that the Fe-L determination is reliable. Some of the stronger Fe-L lines due to Fe XXII and Fe XXIV are close to the K-lines of Ne and Mg, respectively and blending can lead to errors in the Ne abundance and, to a smaller extent, to the Mg abundance (Liedahl et al., 1995; Mushotzky et al., 1996). Also the Ni measure is difficult due to the possible confusion of its L-lines with the continuum and Fe-L blend.

### 5.2 Abundance profiles

In Figure 2 we report the MOS radial abundance profiles for O (top panel), Si (middle panel) and Fe (bottom panel), in Figure 3 those for Mg (top panel), Ar (middle panel) and S (bottom panel), in Figure 4 those for Ne (top panel), Ca (middle panel) and Ni (bottom panel). We note that the measurements obtained using the two different plasma codes agree for what concerns Fe, Ar and Ca; they are somewhat different for what concerns O, S and Ni and in complete disagreement for Mg and Ne.

The temperature profile obtained with the two codes is showed in Figure 5: there are some differences in the inner regions, while in the outer isothermal bins there is substantially agreement.

The models using the APEC code give a systematically worse description than the ones using MEKAL code: for the multi temperature models with APEC the $`\chi ^2`$ range from 403 to 498 for $`\nu =218`$ (216 in the central bin due to the two additional degree of freedom of the power-law component), while for MEKAL models (2T and fake multi-phase give the same results) the $`\chi ^2`$ range from 323 to 382 for $`\nu =217`$ (215 in the central bin due to the two additional degree of freedom of the power-law component); for single temperature models with APEC the $`\chi ^2`$ range from 546 to 886 for $`\nu =222`$, while for MEKAL models the $`\chi ^2`$ range from 326 to 731 for $`\nu =221`$. In Figure 6 we compare the residuals in the form of $`\mathrm{\Delta }\chi ^2`$ between a 2T model using the MEKAL code and the APEC code for the inner bin 1-2, as well as a 1T model using the MEKAL code and the APEC code for the “isothermal” bins 3-4 and 11-14. It’s evident that in the external bins the differences in the fit between the two codes are due to APEC over-prediction of the flux of Fe-L lines from high ionization states, considering the fact that the temperature obtained by the two codes are nearly coincident. For the inner regions the differences between the two codes are further complicated by the different temperature range they find for the best fit. In general, where the temperature structure is very similar, as in the innermost bin, the difference is as in the outer bins in the high energy part of the Fe-L blend, while where the temperature structure is different, as in the 1-2 shown in Figure 6, the differences between the two codes are primary due to different estimates of the flux of He-like Si-K line.

We therefore choose as our best abundance profiles those obtained with a 2T vmekal fit for the central regions and with a 1T vmekal fit for the outer regions. In Table 1 we report the abundance profiles for O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Fe and Ni obtained in this way. Abundance gradients are clearly evident for Fe, Si, S, Ar and Ca; O, Mg and Ni show evidence for an enhancement in the central regions while only Ne is substantially flat.

### 5.3 Comparison with Finoguenov et al. (2001) results

We made a direct comparison of our results with the abundances and abundance ratios derived by Finoguenov et al (2001). Their values are consistent within the errors with ours except for the oxygen abundance, which is roughly two times higher in Finoguenov’s analysis. To better understand the origin of this discrepancy we have extracted spectra from the same radial bins, 1-3 and 8-14, and fitted a single temperature model (vmekal) in the 0.5-10 keV band, excising the 0.7-1.6 keV energy range, in order to avoid the dependence from the Fe L-shell peak, as done by Finoguenov et al (2001). Our results for the O abundance, given in units of the solar values from Anders & Grevesse (1989) to make a direct comparison with the results of Table 1 of Finoguenov et al (2001), are $`0.32\pm 0.03`$ for the 1-3 bin and $`0.20\pm 0.02`$ for the 8-14 bin, to be compared with $`0.535_{0.021}^{+0.019}`$ and $`0.386_{0.021}^{+0.025}`$. Also interesting are the results for the Ni abundance: $`2.55\pm 0.33`$ for the 1-3 bin and $`2.34\pm 0.26`$ for the 8-14 bin, to be compared with $`2.573_{0.918}^{+0.924}`$ and $`0.800_{0.800}^{+1.732}`$. We also note that a 2T modeling of the inner bin 1-3 gives a statistically better fit over a single temperature model ( $`\chi ^2`$/d.o.f of 919/440 respect to 1083/442, using the 0.5-10 keV band) breaking down the assumption of near isothermality.

We also cross compare the K- and L-shell results for Ni and Fe in our analysis, performing the spectral fits for all the radial bins in the 0.5-10 keV band, but excising the 0.7-1.6 keV energy range. The derived abundances for the two elements are consistent within $`1\sigma `$, although the K-shell Ni abundance is 20$`\%`$ higher than the L-shell abundance in the bins fitted with a single temperature model. It should be borne in mind that particularly for the outer bins the K-shell Ni abundance measure is very sensitive to the background estimate.

### 5.4 Abundance ratios and SNIa Fe mass fraction

From the abundance measurements we obtain the abundance ratios between all the elements relative to Fe, normalized to the solar value. They are shown in Figure 7, Figure 8, Figure 9 and in Table 2, together with the abundance ratios obtained by models of supernovae taken by Nomoto et al. (1997) and rescaled to the solar abundances reported in Grevesse & Sauval (1998). We use those abundance ratios to estimate the relative contributions of SNIa and SNII to the metal enrichment of the intra-cluster gas. Such estimates are complicated by uncertainties both in the observations and in the theoretical yields. Our approach is to use the complete set of ratios trying to find the best fit of the function

$$\left(\frac{X/Fe}{X_{}/Fe_{}}\right)_{observed}=f\left(\frac{X/Fe}{X_{}/Fe_{}}\right)_{SNIa}+(1f)\left(\frac{X/Fe}{X_{}/Fe_{}}\right)_{SNII}$$
(1)
where $`\left(\frac{X/Fe}{X_{}/Fe_{}}\right)_{observed}`$ is the measured abundance ratio of the $`X`$ element to $`Fe`$, given in solar units, $`\left(\frac{X/Fe}{X_{}/Fe_{}}\right)_{SNIa}`$ and $`\left(\frac{X/Fe}{X_{}/Fe_{}}\right)_{SNII}`$ are the theoretical abundance ratio by the two types of supernovae, also given in solar units and $`f`$ is directly the SNIa Fe mass fraction. The result of the simultaneous fit of the eight ratios is presented as circles in Figure 10, using the SNII model by Nomoto et al. (1997) and W7, WDD1 and WDD2 models for SNIa respectively for the three panels. Due to the large uncertainties in the yields for the SNII model, the results are strongly SNII model dependent. For comparison we use the range of SNII yields calculated by Gibson et al. (1997), also listed in Table 2, which involve only ratios for O, Ne, Mg, Si and S and as before finding the best fit for eq.(1). The results are shown as squares for the lower end of the range and triangles for the upper part of the range.

The best fits are obtained in the inner bins with a combination of the WDD2 model for SNIa and the Nomoto model for SNII, with reduced $`\chi ^2`$ which ranges from 2 in the inner bin up to 10 in the 4-5 bin. In the outer bins the fit is slightly better (reduced $`\chi ^2`$ of 8-10 instead of 12-13) with a combination of the W7 model for SNIa and the Nomoto model for SNII. This is shown in Figure 11 where the fits with the three SNIa models together with the Nomoto SNII model are reported for the 0.5-1 bin and for the 11-14 bin. It’s clear from the inspection of the residuals in terms of $`\mathrm{\Delta }\chi ^2`$ that the combination of W7 model and Nomoto SNII model fails in the inner bins because it predicts a higher Ni/Fe and O/Fe ratio, compared to the delayed detonation models. The situation is the opposite for the outer bins where a higher Ni/Fe and a lower S/Fe favors the W7 model. However the preference for the W7 model in combination with the Nomoto SNII model in the outer bins is strongly dependent on the Ni/Fe ratio: if we exclude it from the fit the WDD2 model provides the better fit to the data.

## 6 Discussion

The model emerging from the *ASCA* and *BeppoSAX* data for the explanation of abundance gradients in galaxy clusters was that of a homogeneous enrichment by SNII, the main source of $`\alpha `$ elements, maybe in the form of strong galactic winds in the proto-cluster phase and the central increase in the heavy element distribution due to an enhanced contribution by SNIa, strongly related to the presence of a cD galaxy (Fukazawa et al., 1998; Dupke & White, 2000b; Finoguenov et al., 2000; De Grandi & Molendi, 2001; Makishima et al., 2001). As a textbook example we can consider the case of A496 observed by *XMM-Newton* (Tamura et al., 2001). The O-Ne-Mg abundance is radially constant over the cluster, while the excess of heavy elements as Fe, Ar, Ca and Ni in the core is consistent with the assumption that the metal excess is solely produced by SNIa in the cD galaxy. The crucial ratio for the discrimination of the enrichment by the two types of supernovae, O/Fe, is then decreasing towards the center.

The *XMM* results for M87/Virgo question this picture. They confirm and improve the accuracy of the measure of heavy elements gradients previously found by *ASCA* and *BeppoSAX*, but they also show a a statistically significant enhancement of $`\alpha `$ elements O and Mg in the core. If we consider the inner 9 arcmin the ratio O/Fe is constant with $`\chi ^2=6.9`$ for 7 d.of. and adding a linear component does not improve the fit ($`\chi ^2=5.4`$ for 6 d.o.f). These facts points toward an increase in contribution also of SNII, since O is basically produced only by this kind of supernovae. Although there is little or no evidence of current star formation in the core of M87, the O excess could be related to a recent past episode of star formation triggered by the passage of the radio jet, as we see in cD galaxies with a radio source (A1795 cD: van Breugel et al. 1984 and A2597 cD: Koekemoer et al. 1999), nearby (Cen A: Graham 1998) and distant radio galaxies (van Breugel et al., 1985; van Breugel & Dey, 1993; Bicknell et al., 2000) (for a comprehensive discussion see McNamara 1999).
To put the above idea quantitatively, waiting for a true deprojection of our data, we use previous *ROSAT* estimate of the deprojected electron density in the center of the Virgo cluster (Nulsen & Böhringer, 1995) to calculate the excess mass of oxygen. To estimate the excess abundance we fit the inner bins, where we see the stronger increase in the O abundance, with a constant obtaining an abundance of 0.32, while for the outer bins we obtain an abundance of 0.21 so the excess is 0.11. Then we estimate the oxygen mass, $`M_\mathrm{O}`$ to be $`M_\mathrm{O}=A_\mathrm{O}y_{\mathrm{O},}Z_\mathrm{O}^{excess}M_H`$, where $`M_H=0.82n_e\frac{4\pi }{3}(R_{out}R_{in})^3`$, $`R_{in}`$ and $`R_{out}`$ being the bounding radii in kpc of the bins, $`A_\mathrm{O}=16`$ and $`y_{\mathrm{O},}=6.76\times 10^5`$ from Grevesse & Sauval (1998). We obtain a rough estimate of $`10^7\mathrm{M}_{}`$ which, assuming that about $`80\mathrm{M}_{}`$ of star formation are required to generate a SNII (Thomas & Fabian, 1990) and Nomoto SNII oxygen yield, requires a cumulative star formation of $`45\times 10^7\mathrm{M}_{}`$. This star formation is in agreement with a burst mode of star formation ($`10^7`$ yr) at rates of $`1040\mathrm{M}_{}\mathrm{yr}^1`$, as it is observed in the CD galaxies of A1795 and A2597 (McNamara, 1999).

For what concern heavy element gradients and the contribution of SNIa, M87 data suggests an agreement with delayed detonation models (in particular for the inner bins), as stressed by Finoguenov et al (2001), in contrast with the preference of W7 model set by the high Ni/Fe ratios found by Dupke & White (2000a). This is particular evident if we consider the S/Fe ratio in Figure 7. If we consider the set of theoretical values for SNII and W7 SNIa models the behavior of these ratio would indicate an increasing contribution by SNIa going *outward* to the center. We recover the correct behavior if we choose the WDD1 yield and we reduce the S SNII yield of Nomoto et al. (1997) by a factor of two to three, as was already indicated by *ASCA* data (Dupke & White, 2000a). We caution however that this is a substantial contribution larger than that allowed by the SNII models choosed by Gibson et al. (1997). The use of delayed detonation model for SNIa could also explain the over abundance of S and also Si (respect to the W7 model) found by Tamura et al. (2001) for the core of A496.
With increasing radius the W7 model gives a better fit to the data respect to delayed detonation models. This fact indicates a SNIa abundance pattern change with radius and could be taken as an independent X-ray confirmation of the conclusions of Hatano et al. (2000) on the optical spectroscopic diversity of SNIa, as suggested by Finoguenov et al (2001). However we stress that the preference of W7 over WDD models in the outer bins is entirely due to the Ni/Fe ratio which could be affected by systematic uncertainties, as discussed in section 5.1 and 5.3.

In Figure 10, we show the relative importance of SNIa, measured by the Fe mass fraction provided by this kind of supernovae. This is substantially constant through the 14 arcmin analyzed. This fraction is considerable and ranges between the 50% and 80% and it depends only slightly on the SNIa model used. Instead the uncertainties involved in using different SNII models are large and a definitive answer cannot be reached until further convergence of SNII models is achieved.

## 7 Summary

We have performed a spatially resolved measurement of the element abundances in M87, the Virgo cluster cD galaxy. The main conclusion of our work are:

* the APEC code gives a systematically worse description than the MEKAL code in modeling M87 spectra;
* we confirm the increase of Fe and other heavy elements towards the core indicating that the SNIa contribution increases;
* the increase in O abundance and a constant O/Fe ratio in the inner 9 arcmin indicates an increase also in SNII ejecta possibly from star-burst in the recent past;
* Si/Fe and S/Fe profiles favor WDD models over W7, also requiring substantial reduction of the SNII yield of S;
* the indication of a change of the SNIa abundance pattern, provided by a preference of the W7 model over delayed detonation models in the outer bins, is entirely due to the Ni/Fe ratio. Since the Ni measurement is difficult and uncertain, this indication should be taken with some caution.

S. Ettori and S. Ghizzardi are thanked for useful discussions and suggestions. We thank the referee for several suggestions that improved the presentation of this work. This work is based on observations obtained with *XMM-Newton*, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA).


text:
# Adaptive Detection of Instabilities: An Experimental Feasibility Study.

## 1 Introduction

Over the last 30 years increasing efforts have been directed towards the understanding of experimentally observed complex dynamical behavior. The development of Nonlinear Dynamics as a branch of applied mathematics provides us the set of tools necessary for understanding the response of nonlinear dynamical systems as operating parameters are varied. Characterization of the nonlinear dynamical behavior of experimental systems, however, remains a difficult task. This is mainly because real systems cannot be manipulated as easily or as cleanly as a model would on the computer when its bifurcation behavior is traced.

Typically, an experimental bifurcation diagram is obtained via repeated application of a set-and-observe procedure. The operating parameter of interest is set and the system dynamics allowed to evolve asymptotically to stationary (e.g. steady or oscillatory) behavior. The operating parameter setting is then changed (e.g. incremented slightly) and the experiment is repeated; if a qualitative difference is detected between the asymptotic dynamics reached for two consecutive parameter settings, a bifurcation is inferred to have occurred somewhere in-between. A similar approach can be applied to the computational detection of bifurcations if one is only allowed to use forward-in-time integration of the model equations (i.e. set parameter values and initial conditions and run a simulation to stationarity). In both cases there are serious drawbacks. For example, in the investigation of a period doubling bifurcation of a periodic orbit, as the critical value of the search parameter is approached the dynamics will slow down considerably, and the time required to determine whether ultimately single loop or double loop oscillations are observed increases dramatically. Furthermore, if one wishes to detect a subcritical bifurcation, or a saddle-node bifurcation, once one steps over the critical parameter value, the system will move far away from the phase space region of interest. The experiment will then need to be reset to reach the vicinity of the stable solution located just before the critical parameter value that was overstepped.

In previous work  we proposed a protocol for the experimental detection of bifurcations motivated by the procedure of computational construction of bifurcation diagrams from mathematical models. Numerical bifurcation theory in a computational environment allows the bifurcation location via a simultaneous search for criticality in (phase)$`\times `$(parameter) space. An augmented set of equations (e.g. steady state equations plus a criticality condition) is typically formed, and its fixed point (steady state values at criticality and critical parameter value) is converged upon quadratically using a Newton-Raphson contraction mapping for the augmented system of equations. Obviously in this case the states do not evolve on the computer based on the system dynamics, but based on the computational contraction mapping dynamics.

The proposed protocol for the detection of bifurcations seeks an “experimental-computational” compromise: an experiment (i.e. states evolving based on system dynamics) is coupled to a parameter evolution law which is based on criticality conditions and feedback control. The augmented, closed loop system, constructed in this manner has a stable steady state at what used to be the open loop system’s marginally stable steady state. The computational building blocks that constitute the parameter evolution law are based on the on-line identification of a “local nonlinear” dynamic model, obtained via least-squares regression , used to estimate criticality on-line, and simple feedback control algorithms that use the bifurcation parameter as the actuator in order to drive the system towards the bifurcation point. As it will be illustrated below, this algorithm enhances our ability to efficiently detect instabilities, trace segments of the bifurcation diagram (“operating diagram”), and, in general, characterize the system dynamics in multiparameter space. The “experimental” system used in our illustration is an analog realization of the Rössler system. For this system, we converge experimentally on the period doubling bifurcation of simple periodic orbits in one and two-parameters. The local adaptive model used to infer the presence of the bifurcation is a 2D nonlinear map. Our illustrations are based on (a) full state vector measurement; and (b) reconstruction of the state-space based on time series of a single state measurement and the use of delay coordinates.

An alternative motivation for studying this type of problem is the typical Ziegler-Nichols (ZN) tuning procedure for PID (proportional- integral-derivative) controllers , which motivated the work of O’Neil, Lyberatos and Svoronos  on adaptive determination of bifurcation points. In the “textbook” ZN procedure the “experiment” is the closed loop system with a proportional controller; the gain of this controller is the “bifurcation parameter”. The closed-loop system loses stability when the proportional gain $`K`$ reaches the value $`K_u`$. This critical parameter value is precisely the ZN “ultimate gain”. In practice, one obtains $`K_u`$ by slowly increasing the gain of the proportional controller until one observes incipient oscillations. Knowledge of the location of the bifurcation and the frequency of the oscillations dictates then the ZN-tuning settings for a PID controller; it is interesting that this practical controller tuning approach is based on experimental bifurcation detection.

The paper is organized as follows: first, we describe our system (the electronic Rössler circuit) and the bifurcation detection protocol before presenting our experimental results. We then provide some concluding remarks.

## 2 The Rössler Circuit

The Rössler model was proposed in the early 1970’s  as one of “the simplest” nonlinear vectorfields capable of generating chaotic behavior. It consists of the following set of three ordinary differential equations:

$`{\displaystyle \frac{dx}{dt}}`$ $`=`$ $`yx`$ (1)
$`{\displaystyle \frac{dy}{dt}}`$ $`=`$ $`x+ay`$ (2)
$`{\displaystyle \frac{dz}{dt}}`$ $`=`$ $`b+z(xc)`$ (3)
This system has been extensively studied. For $`b=2`$ and $`c=4`$ and as $`a`$ is increased, it exhibits a steady state undergoing a Hopf bifurcation followed by a cascade of period doubling bifurcations of the periodic orbits culminating in an apparently chaotic attractor around $`a=0.396`$.

In the implementation of the protocol for the adaptive detection of experimental bifurcations, we use as model system an analog realization of the Rössler system. The Rössler equations are implemented via a sequence of feedback loops of electronic signals. The circuit was manufactured at the FHI labs111 A detailed diagram of the circuit is available in the web in the following address http://w3.rz-berlin.mpg.de/pc/spatdyn/spatdyn.cgi?Publications and includes an encasing that allows both manual and digital manipulation of all three system parameters. In our experiments the parameter dynamics were computer-generated in closed loop, and were communicated to the device through its digital I/O channels. The circuit is coupled (Fig. 1) with a VMEbus System (32bit backplane system) consisting of a VMEbus CPU-board Baja40 from Heurikon Division controlled via a Motorola 68040 32-bit microprocessor chip running internally with 66 MHz, externally with 33 MHz, 32 MByte RAM and several I/O interfaces and a VMEbus IP carrier board with four IndustryPack slots. One slot was used with a TIP501 optically isolated 16 channel 16 bit ADC module from TEWS Datentechnik, input voltage +/- 10 volt, 14.5 usec data acquisition and conversion time, accuracy +/- 4 LSB, linearity +/- 4 LSB, and the other slot with a TIP550 optically isolated 4 channel 12 bit DAC from TEWS Datentechnik, output voltage +/- 10 volt, 13 usec settling time, accuracy +/- 1 LSB, linearity +/- 1 LSB. The array uses an operating system VxWorks 5.3.1 (Tornado) from Wind River System. The data acquisition components allow us to specify the accuracy of the measurements from 11 to 16 bits.

The dynamic response of the motivating set of differential equations is well reproduced. Figure 2 presents a comparison of a couple of projections of the double-loop limit cycle predicted and observed at $`a=0.35`$, $`b=2`$ and $`c=4`$. The deviation of the signal given by the circuit with respect to the one obtained by integrating the differential equations is minimal, both in the actual values of the variables as well as in the period of the oscillations. Similar results are observed within the range of operation of the circuit: $`a[0,1]`$, $`b[0,10]`$, $`c[0,10]`$ and values of the state variables inside the range \[-10,10\].

## 3 The Bifurcation Detection Protocol

The building blocks that constitute the bifurcation detection protocol are schematically depicted in Figure 3. They involve four different modules: one that identifies on-line an approximate model of the system, a second one that uses the identified model to estimate the location of the bifurcation point in phase$`\times `$parameter space; this second module relies in the information provided by a test function (third module) that contains the conditions defining the type of bifurcation sought. The fourth and final module uses the estimated critical parameter value (our estimate of the open-loop critical conditions) to devise a policy that will bring and keep the system closer to open-loop criticality. The repeated application of the algorithm, composed by these four modules, drives the system to the true open-loop bifurcation point, and holds it there. In devising the algorithm, we assume that at least one of the states of the experimental system can be observed (measured), and that at least one operating parameter is available to be manipulated (the bifurcation will be located with respect to this parameter). Furthermore, we assume that the open-loop system is in the neighborhood of a local bifurcation; thus its dynamics can be captured via a low-order nonlinear model. Using such a low-order approximation to model the dynamics of the system is motivated (and can be justified) by the existence of a normal form for the bifurcation. The discrete-time model used in the system identification module of the algorithm ((1) in Fig. 3) takes the form of a low order polynomial on the observed variables and the parameter, motivated from Taylor series expansions in the neighborhood of the bifurcation point, and the theory of Center Manifolds/Normal Forms. The polynomial model is of the type:

$`x_i(k+1)`$ $``$ $`F_i[𝐱(k);p(k)]`$
$`=`$ $`a_i+{\displaystyle \underset{j=1}{\overset{n}{}}}b_{i,j}x_j(k)+c_ip(k)+{\displaystyle \underset{j=1}{\overset{n}{}}}{\displaystyle \underset{l=1}{\overset{n}{}}}d_{i,j,l}x_j(k)x_l(k)`$
$`+{\displaystyle \underset{j=1}{\overset{n}{}}}e_{i,j}x_j(k)p(k)+{\displaystyle \underset{j=1}{\overset{n}{}}}{\displaystyle \underset{l=1}{\overset{n}{}}}{\displaystyle \underset{m=1}{\overset{n}{}}}f_{i,j,l,m}x_j(k)x_l(k)x_m(k)(i=1,\mathrm{},n)`$
where $`𝐱`$ represents the $`n`$-dimensional vector of state variables, $`p`$ is the parameter, $`a,b,c,d,e`$ and $`f`$ are the polynomial coefficients, and $`k`$ is the iteration counter. Note that the polynomial is linear in the parameter and contains cubic terms on the state variables. The normal forms of the most common bifurcations are linear in the parameter, and the normal form of the Hopf bifurcation has cubic terms in the state variables . The polynomial can be simpler for certain specific types of bifurcation; the form presented here is appropriate for the detection of local codimension-1 bifurcations. For the experimental system, we rely in Poincaré sections and/or reconstructions based on delayed measurements to form the vector of states. The identification module uses a least-squares algorithm  to fit the coefficients of the polynomial from experimental observations; other algorithms, such as the projection method, can of course also be used. The model is updated at every step as new observations become available. In order to prevent ill-conditioning of the least-squares algorithm during the experimental run, a well-known problem in the on-line identification of experimental systems, the parameter signal must be persistently excited (here, by addition of random noise) .

The second module of our algorithm ((2) in Fig. 3) uses the polynomial model identified by the first module to estimate the location of the open-loop bifurcation. This is a computational step that involves the simultaneous solution of the steady state equations along with a criticality condition appropriate for the bifurcation sought ((3) in Fig. 3). This is accomplished using an algebraic equation solver (usually a contraction mapping like a Newton-Raphson). For example, for a period-doubling bifurcation, which we detected in the Rössler circuit, the bifurcation is located at the point where one Floquet multiplier (an eigenvalue of the linearization around the fixed point of the open-loop Poincaré map model) is equal to -1. Thus, a suitable test function, for a two-variable discrete-time system, will be $`det[𝐈+𝐉]=0`$, where $`𝐈`$ is the identity matrix and $`𝐉=[_x𝐅[𝐱;p]]`$ is the Jacobian of the identified model at the steady state $`(𝐱,p)`$.

In order to find the location of the critical state ($`\{\widehat{𝐱}_{cr},\widehat{p}_{cr}\}`$), one must solve for $`\{𝐱,𝐩\}`$ such that (period doubling case):

$`x_iF_i[𝐱;p]`$ $`=`$ $`0(i=1,\mathrm{},n)`$
$`det[𝐈+𝐉]`$ $`=`$ $`0`$
The final module ((4) in Fig. 3) uses the identified model and the estimated location of the bifurcation ($`\{\widehat{𝐱}_{cr},\widehat{p}_{cr}\}`$) to find a parameter variation protocol that will drive the system to the bifurcation point. This critical step can be accomplished using several alternative “feedback laws”. Previously , we had illustrated the use of one of the simplest: a policy that will bring the system to criticality in the minimal possible number of steps (dead-beat policy). In the runs presented below, however, we rely on “softer” policies in order to overcome problems associated with the local character of the model identified, the limited resolution of measurements and the need for continuous excitation of the control signal associated with the ill-conditioning of the least squares identification procedure. In addition, the control policy should be possible to be evaluated and implemented very fast (almost instantaneously) in real time. This is because we have ignored in our experiment the delays associated with I/O tasks as well as those associated with the computation of the control action - we have effectively assumed that computation and implementation of control action is instantaneous. This avoids a number of complications and possible instabilities having to do with physical delays in the control loop. These are important issues, and will be taken up in future work; they are, however, more “implementation difficulties”, and do not directly pertain to the mathematical underpinnings of our procedure.

With these restrictions in mind, we explored a couple of plausible two parameter variation protocols: the first policy attempts to minimize the distance from the critical point in the next iteration, penalizing large parameter variations; the second protocol is based on decreasing the distance to the critical state (in phase space) in the next iteration by a fixed fraction. In the notation used above, the first protocol results from:

$$\underset{p(k)}{\mathrm{min}}\underset{i=1}{\overset{n}{}}(x_i(k+1)\widehat{x}_{i,cr})^2+W(p(k)\widehat{p}_{cr})^2$$
(4)
Note that, since the identified model (used to predict $`𝐱(k+1)`$) is linear in the parameter $`p`$, the solution of this minimization problem renders a linear equation for $`p`$. The parameter $`W`$ is a penalty weight, which we usually set in the range \[0.1,10\]. Since the model form used is only valid locally, close to the open-loop bifurcation, its predictions away from the bifurcation are not accurate. In the presence of modeling errors, a dead-beat policy may drive the system farther away from the bifurcation instead of bringing it closer. The introduction of the weight $`W`$ gives us a mechanism to moderate excessive control action, and thus, to possibly overcome the effect of such modeling errors away from the bifurcation point. A more appropriate implementation of this policy would seek only to penalize actions that bring the system outside the range of validity of the model. That is, the penalty should only be applied if the prescribed control action is larger that the estimated range of validity of the identified model. We are currently exploring this issue.

The second policy is motivated by the form of transients approaching a steady-state near a Hopf bifurcation. In such a scenario the dynamical system will slowly converge to a marginally stable steady state by (possibly fast, but slowly decaying) spiraling towards it. Thus, this second policy attempts to control the speed of the approach by requiring that the distance from the critical state location in the next iteration becomes a fraction of the current distance from criticality. In mathematical terms we solve for $`p(k)`$ such that:

$$\underset{i=1}{\overset{n}{}}(x_i(k+1)\widehat{x}_{i,cr})^2=C^2\underset{i=1}{\overset{n}{}}(x_i(k)\widehat{x}_{i,cr})^2$$
(5)
Here $`C`$ is the desired fractional reduction, in the next iteration, of the distance from estimated open-loop criticality. We typically set $`C`$ in the range \[0.9,0.95\]. Since the model to predict $`x_i(k+1)`$ is linear in $`p`$, the solutions of this equation are the roots of a quadratic equation that can be solved explicitly. The root selected is the one closer to the current value of the parameter.

Both control policies must be applied repeatedly in order to drive the closed loop system to the point that, for the open-loop system, was a marginally stable solution. Although the number of iterations to reach the bifurcation may be large (e.g. more than 100 iterations), we have found that, even with these simple policies, the closed-loop system performs well, finding the bifurcation for every experiment independently of the initial condition used.

By using the identification and control methods in the manner described, we circumvent the long settling periods associated with open-loop near-critical dynamics. As previously illustrated in  through computational modeling, this methodology is also capable of driving the system back to steady state after (mildly) overshooting the critical value of the bifurcation parameter even for “hard” bifurcations such as turning points or subcritical Hopf bifurcations. In the following section, we illustrate our approach by locating period-doubling bifurcations of periodic orbits for the Rössler circuit device.

## 4 Experimental Runs.

A period-doubling cascade of periodic orbits (or fixed points of maps or Poincaré maps) is perhaps one of the best understood routes to deterministic chaos. The Rössler system exhibits this type of transitions at $`b=2`$ and $`c=4`$ for increasing values of the remaining parameter $`a`$. The first of such period doubling bifurcations is observed for $`a0.33485`$. It is precisely this bifurcation that we set to locate using the protocol described above.

Poincaré sections provide a natural framework for observing/reporting period doubling of periodic orbits. The construction of our approximate models is performed on such a (transversally intersecting the flow) Poincaré section. Initially, we measure all three state variables from the circuit $`(x,y,z)`$ and define the Poincaré section as the plane $`x=0.2`$, crossed towards increasing values of $`x`$. For this Poincaré section the critical point is located at $`y=4.1409`$ and $`z=0.44601`$. In the proposed protocol, the model to be identified involves two state variables (the coordinates $`(y,z)`$ on the Poincaré section) and one parameter ($`a`$). Thus, our model will try to predict where the trajectory intersects the plane at the next iteration ($`y_{k+1},z_{k+1}`$) based on the current location ($`y_k,z_k`$) and the value of the parameter ($`a_k=p_k`$). For the experiment, the intersection point on the surface defining the Poincaré surface was obtained via fast (between 200 to 1000 Hz) sampling of all variables and linear interpolation between the sampled points once the crossing with the surface was detected.

For the circuit, the period of the oscillations before the bifurcation is approximately 0.62 s. In comparison with the Rössler ODE model, the arbitrary units of time translate as seconds for the circuit device and the arbitrary units of the state variables and parameters as voltages measured in Volts. Using a 16-bit data acquisition card, the measurement resolution for the state variables is approximately 0.0006 volts and 0.00003 volts for the parameter $`a`$ which has a smaller range of variation (voltage). As we mentioned earlier, in order to successfully identify a varying on-line model, one should apply persistent excitation to the actuator (the parameter signal). This is done in the form of random noise with maximum amplitude of 0.003 volts. The excitation signal prevents the system from reaching stationary states for which the measurement will be effectively constant. Such measurements will render the matrix used in the identification computations severely ill-conditioned, and the approximate model will deteriorate. With persistent excitation signals the system will be always “on the move” allowing the collection of information-rich data resulting in more accurate models. Since, however the final goal is to drive the system to the open-loop bifurcation, very large excitation signals are undesirable because they work against the control actions implemented. The excitation level selected represents a compromise between these competing effects. Depending on the system under consideration, the magnitude of the excitation signal should be adjusted. In our experience with several different simulated systems, we have found that the amplitude of the excitation signal must be at least one order of magnitude larger than the resolution available in the measurement.

The use of this relatively large excitation signal combined with the soft control policies described previously constrains the selection of the appropriate convergence criterion. The approach to the critical state may be slow, possibly also “erratic” due to the effect of the excitation; on the other hand, the system may become trapped and the limited resolution may translate in biased models and incorrect criticality predictions.

In order to prevent the latter occurrence, we opted for a two-level convergence criterion: (1) we require that the norm of the deviation from criticality in the extended phase$`\times `$parameter space is sufficiently small (i.e. less than 0.008 volts) for a relatively small number of consecutive iterations (typically 5); and (2) after complying with the first criterion, we restart the system identification procedure by gathering enough data to construct a “fresh” model, using the excitation signal without control, and then restart the feedback control until the first criterion is again satisfied. Overall convergence is declared only if two consecutive restarts reach the same point within the tolerance indicated above. The first criterion does not excessively burden the controller performance, while the second seeks to prevent declaring convergence at a possibly unconverged, incorrect prediction of the location of the bifurcation. The tolerance to declare convergence (0.008 volts in our case) should clearly be larger than the set maximum amplitude of the random noise (excitation signal, 0.003 volts); a reasonable rule-of-thumb is to set this tolerance half-an-order of magnitude above the maximum amplitude of the random noise. We note, however, that the achievable accuracy in locating the bifurcation is limited, not because of the convergence criterion, but most importantly because of the excitation signal, model sensitivity and, ultimately, limited resolution in the measurements.

With the feedback protocols implemented, we observed, during our runs, that large parameter variations (“control actions”) arose frequently. Since such variations may drive the system outside the range physically accessible to the circuit device (i.e. negative parameter values), or outside the range for which sustained oscillations are observed (thus, preventing the trajectories from intersecting the plane used in the definition of the Poincaré section) we opted for establishing hard bounds (“saturation”) for the values that the parameter can take. If the control action prescribed goes beyond these bounds, the parameter is set to the value of the bound that has been surpassed (the actuator saturates). For the run at $`b=2`$ and $`c=4`$, the minimum value allowed for $`a`$ was 0.12 and the maximum $`0.5`$. Below $`a0.125`$ the only attractor of the system is a steady state, while above $`0.5`$ the chaotic attractor observed in the range \[0.396,0.5\] disappears via global bifurcations and the state variables of the circuit tend to saturate.

With the hardware available to us, the evaluation of the protocol, including identification of an updated model, location of criticality and evaluation of the control policy, involved between 30 and 60 milliseconds, depending on the number of evaluations required for some of the subroutines such as the algebraic solver to find criticality. The effect of such “large” delay (about one tenth of the period of oscillation in the worst case) was often detrimental to the performance of the protocol, preventing the accurate location of the bifurcation. In order to overcome this situation, we modified the evaluation order of the subroutines. That is, we first calculate and implement a control policy using one-iteration-old models, then we use the collected data to update the model and prediction of criticality for the next iteration. This strategy was successful, since the evaluation of the control policies proposed require only the solution of linear or quadratic equations and thus take little time (less than 3 ms), and the slow approach to the critical state makes consecutive models very similar most of the time (depending on the conditioning of the matrix involved in the identification of the model).

Figure 4 illustrates the performance of the protocol. The parameter was initially set at $`a=0.42`$. At this parameter setting, the open-loop circuit exhibits a chaotic attractor. The experiment is left to evolve in the vicinity of this parameter value for about 30 iterations, while data are gathered to construct an initial model. After approximately 70 iterations, the system has been driven to the vicinity of the bifurcation. The last 50 iterations, for which the response of the circuit is practically constant, represents the period needed to declare convergence with one restart as described above. For this run the critical state was predicted at $`y=4.1127`$, $`z=0.4537`$ and $`a=0.3389`$. This compares favorably with the location of the bifurcation for the Rössler ODE system. As we indicated before, we anticipated that the estimation of the bifurcation location will not be very accurate because of the persistent excitation signal and soft control policies. A more representative result should be reported after gathering statistics over a number of runs. On average, over 15 runs, each run with different initial condition and initial parameter value, the bifurcation is predicted to be located at $`a=0.3309`$ with a standard deviation of $`0.0082`$ volts, with the critical state located at $`y=4.0964`$ and $`z=0.4552`$ and standard deviations of $`0.0134`$ and $`0.0014`$ volts respectively. The deviation of this average is well within the range allowed for declaring convergence in our runs and the standard deviation commensurate with the excitation signal and convergence criteria. We note that all our runs did converge to the critical state albeit occasionally at a larger number of iterations than for Figure 4.

Figure 5 presents another run that illustrates the role of hard, saturation bounds limiting the values that the parameter can take during the search for criticality. After approximately 40 steps, an incorrect criticality prediction translates into a drastic control action. Since the value prescribed for the parameter would be below the preset lower bound, the parameter is set to the bound. The performance of the controller deteriorates and the system is trapped near the lower bound for about 40 iterations before restarting its movement towards the critical state.

Figure 6 shows an example of the phenomenon commonly referred to as “bursting”. The system starts its gradual approach toward criticality, but since the coordinates of the points in the trajectory do not change much, the matrix used for the identification becomes progressively more ill-conditioned, eventually rendering a “bad” model that, in turn, dictates a large parameter excursion (“burst”). The trajectory to criticality exhibits, in this example, a couple of instances of such bursting on its route to the eventual correct detection of the critical fixed point.

For those runs initialized at parameter values lower than the bifurcation value, as the one illustrated in Figure 6, the convergence to the critical point was consistently declared for slightly lower values of the parameter than the true critical parameter value. In contrast, when the initialization took place at conditions for which the open loop system exhibited apparently chaotic oscillations (i.e. Figure 4), the collection of information-rich data allowed the accurate estimation of the critical fixed point location early on, resulting in faster convergence towards the bifurcation.

One potential limitation of the application of this strategy may be the available measurement resolution. The examples presented above rely on a 16 bit A/D card. However, as illustrated in Figure 7, the bifurcation was successfully located also for lower resolution, equivalent to that provided by a standard 11 bit A/D converter (i.e. one order of magnitude less in resolution for the parameter and state variables).

Although the examples discussed until this point involve measurements of all the states of the system, it is also possible to use only a single measurement time series. Figure 8 illustrates the performance of the protocol for the location of the same bifurcation using only measurements of the variable $`x`$ and relying on delayed measurements to reconstruct the phase space. The delay ($`\tau `$) for the reconstruction was set to 40 ms. In practice one could seek to find an optimal value of the delay using, for example, mutual information measures; this was not done here. The Poincaré section was set as before using the plane $`x(t)=0.2`$. The critical state is located for the Rössler system at $`x(t\tau )=1.1701`$ and $`x(t2\tau )=2.2020`$. On average over 5 runs, the bifurcation is predicted at $`a=0.3235`$ for $`x(t\tau )=1.1143`$ and $`x(t2\tau )=2.1320`$, with standard deviations of $`0.0080`$, $`0.0042`$ and $`0.0010`$ volts respectively. The trajectory to reach the critical state, however, exhibits a larger number of instances of bursting and saturation of the parameter and thus requires a larger number of iterations as illustrated in Figure 8. This deterioration in performance of the protocol can be attributed to the increased sensitivity of the identification resulting from the use of a sub-optimal delay and the smaller effective resolution of the sampled variable.

A final test of the capabilities of the proposed protocol for the detection of bifurcations comes by attempting to draw, experimentally, a two-parameter bifurcation diagram. That is, we attempt to trace automatically a segment of our codimension-1 bifurcation boundary in two parameters. We vary the (second) parameter $`b`$ in the range . Starting at $`b=3.0`$, we locate the value of $`a`$ for which the bifurcation is observed and then decrease $`b`$ by a fixed amount (0.05 volts). Figure 9 presents the comparison of these runs for our circuit device, and the computed boundary for the Rössler ODE system. The results presented are the average over five sweeps, three calculated decreasing $`b`$ and two increasing it. The previously found bifurcation point was used as the starting point for the new search with a different value of $`b`$. Figure 10 presents the comparison of the location of the (open loop critical) fixed points for the Poincaré section. The boundaries for the Rössler ODE system were calculated using AUTO . As can be seen in Figures 9 and 10 the segment is accurately located, as are the critical states. For both plots, standard deviations have been omitted; they are comparable to the previously indicated values for the single parameter runs.

## 5 Conclusions

We have presented a successful implementation of an algorithm for efficiently driving an experimental system to (what at open loop would have been) criticality. The control protocol was applied to an experiment in the form of an electronic circuit device that implemented the dynamics of the Rössler system. The proposed algorithm is based on adaptive system identification, separation of time scales close to criticality, and on the exploitation of the identified local model to search for the criticality conditions; the estimation of open-loop criticality is used to devise a control policy that drives the closed-loop system towards the bifurcation. The protocol seeks, in a sense, to mimic the procedure used in the numerical construction of bifurcations from mathematical models: the critical state of the open-loop system becomes a stable steady state of the augmented system in (phase) $`\times `$ (operating parameter) space. We illustrated the capabilities of the algorithm in the location of period doubling bifurcations of periodic orbits. The bifurcation was successfully located even in the presence of limited resolution. We also illustrated the applicability of this approach when only one state variable can be measured by reconstructing the phase space using delayed measurements. Our final illustration involved the tracing of a segment of the two parameter bifurcation boundary.

Our selection of an electronic circuit device with Rössler dynamics, to demonstrate the experimental feasibility of the protocol for the automatic detection of bifurcations, was intended to provide a testing ground for assessing the performance of the protocol. The circuit device provides us with a highly reproducible and easily manipulated experiment with a well charted dynamical behavior. We chose the period doubling bifurcations of the oscillations of the circuit as our representative codimension-1 bifurcation. The low-dimensionality of the model used crucially depends on the time-scale gap between the long time scale associated with the critical eigenvalue(s) and the time scale of the “next least stable” eigenvalue. It is important that all operations performed in our system (I/O as well as software computations) are fast compared to this “next least slow” time scale. Since the return time for our Poincaré map was roughly 600 ms, our relatively short (approximately 40 ms) action delay could be considered “approximately instantaneous”. An alternative that would speed up the protocol evaluation is the use of recursive least squares techniques  that may reduce significantly the time involved in the identification of a local model; however, this alternative resulted in more sensitive and prone-to-burst trajectories, even though we attempted several different strategies in the assignment of forgetting factors for the identification. Other alternatives, such as the use of projection algorithms , are also feasible, but were not investigated.

Perhaps the most limiting aspect for the widespread application of the proposed protocol involves the high sensitivity of the strategy as a result of the combination of the control policies, excitation signal and on-line identification of a model with only local validity. Such sensitivity translates in erratic and slow motion towards the bifurcation, frequent bursts, trapping of trajectories near hard bounds and ill-conditioning of the identification process.

We are currently working on posing the entire problem (system, identification and control) as an augmented dynamical system, for which the existence of an overall converged state can be demonstrated, and its stability quantified. Such an approach could also allow us to explore issues related to the performance in the presence of multiple nearby bifurcations of the same type.

Extensions of the algorithm are straightforward: higher codimension bifurcations can be detected by solving the corresponding (more than one) criticality conditions for the (open-loop) marginally stable steady state and the (more than one) critical parameters. In principle one might even be able to accommodate unstable, nonbifurcating parts of the spectrum on the right-hand plane. Several important elements can be improved upon: (a) Formulation of the criticality conditions in a less naive and computationally more efficient manner, along lines of modern numerical bifurcation theory (e.g.  and references therein); (b) Usage of more sophisticated (e.g. receding horizon) feedback control laws ; and (c) Improvement of the type of local model used. While in this problem we used polynomial models based on Taylor-series expansions and normal forms, it is obvious that other types of nonlinear dynamic models (e.g. models based on artificial neural networks) can also be used in this context; we, among many others, have demonstrated that such models are capable of accurately predicting the bifurcation structure of a nonlinear system when trained with time series data .

While we are continuing theoretical and computational work along the directions described above, we are also implementing this approach in the study and characterization of the bifurcations of electrochemical systems .

It may be interesting to note that this approach can also be coupled with large scale scientific simulators, for which the relevant long term dynamics are low-dimensional, in order to chart their bifurcation scenario  (as opposed to writing a bifurcation algorithm for them from scratch). The strategy of seeking and stabilizing a fixed point with one additional condition (here, open-loop criticality), falls in the same category as extremum seeking controllers ; controller developments along these lines would find direct application in our protocols, and vice versa. We believe that an approach like the one described here can eventually become a standard real-time tool which will assist experimentalists in efficiently detecting bifurcations and tracing operating diagrams; such a tool could drastically enhance the experimental characterization of nonlinear dynamics in multiparameter space.

## 6 Acknowledgements

We are grateful to Professor Ertl for his interest in and generous support of this work. We thank Heinz Junkes and Georg Heyne from the FHI technical staff for technical support. We gratefully acknowledge the support of CONACYT and COSNET (RRM), NATO (IGK, KK), and AFOSR (IGK) for partial funding of this research. IGK and RRM gratefully acknowledge the support from the Alexander von Humboldt Foundation.