text:
# The warm absorber of the Seyfert 1 galaxy H1419+480

## 1 Introduction

X-ray spectral studies of radio-quiet Active Galactic Nuclei (AGNs), reveal very frequently the presence of ionised material. This was first noted by Halpern (1984), based on Einstein IPC data, where it was also noted that these absorbers vary. Evidence for ionized absorption along the line of sight was provided by Nandra & Pounds (1992), based on the detection of an absorption edge in $`ROSAT`$ observations of the Seyfert 1 galaxy MCG-6-30-15. Fabian et al (1994) were able to disentangle the OVII and OVIII K-edges of this warm absorber using the much better spectral resolution of $`ASCA`$.

In the years that followed, it became evident that warm absorbers are ubiquitous among Seyfert 1 galaxies. Reynolds (1997) and George et al (1998) showed that more than 50% of Seyfert 1s contain a warm absorber. These are produced by partially ionised gas with total column density in the range $`N_H10^{21}10^{23}\mathrm{cm}^2`$, probably located at or outside the Broad Line Region (Reynolds & Fabian 1995).

With the advent of higher resolution grating spectrographs on board Chandra and XMM-Newton much more detailed studies of ionized absorbers have become possible. X-ray spectra obtained with the Reflection Grating Spectrograph (den Herder et al 2001) on board XMM-Newton have shown that the absorbing gas has often various components with different outflowing velocities and clearly distinct ionisation states (see Sako et al 2001 for IRAS 13349+2438 and Blustin et al 2002 for NGC 3783). One of the most remarkable features discovered in these spectra is the presence of Unresolved Transition Arrays (UTAs) of Fe M lines, which could mimic the shape of absorption edges when observed at the lower spectral resolution typical of CCDs.

Detailed studies of particular AGN which exhibit X-ray ionized absorption features also demonstrate that there is frequently a narrow-line absorber related to the AGN (i.e. an associated absorber) in the ultraviolet band, normally showing up as a CIV $`\lambda \lambda 1548,1550`$ absorption feature (Mathur 1994, Mathur et al 1994, Mathur, Elvis & Wilkes 1995). The same photoionized outflowing gas can provide the UV absorption lines and the ionized absorption features seen in X-ray spectra. More recently, Crenshaw et al (1999) found in a systematic study of UV absorption lines towards nearby radio-quiet AGN, that all objects showing UV associated absorption systems have a corresponding warm absorber in the X-rays.

Variability has also been found in some of these warm absorbers. Fabian et al (1994) reported different absorption optical depths in the best-studied warm absorber towards MCG-6-30-15. Later, Otani et al (1996) showed, using a 4-day long ASCA observation of this particular Seyfert 1, that the OVIII K edge (at a rest frame energy of 0.871 keV) varied on scales of $`10\mathrm{ks}`$, while the OVII K-edge (at a rest frame energy of 0.739 keV) remained constant. Otani et al (1996) interpreted this in terms of two different ionized absorbers, a high ionisation (OVIII) absorber located within the BLR and a lower ionisation absorber (OVII) probably associated to the Narrow Line Region (NLR).

H1419+480 was discovered in the Modulation Collimator - Large Area Sky Survey (MC-LASS) conducted with the HEAO-1 observatory (Wood et al 1984). The inferred 2-10 keV flux was $`2\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, assuming a power law spectrum with a canonical $`\mathrm{\Gamma }=1.7`$ (Ceballos & Barcons 1996). Although the modulation collimator helped to pin down the position of the X-ray source, the flux could be severely affected by the contribution of unresolved sources within the collimator field of view. The $`ROSAT`$ All-Sky Survey (RASS) also detected this source (see Schwope et al 2000), with a 0.5-2 keV flux (corrected for Galactic absorption) of $`7\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ and an additional PSPC pointed observation conducted in 1992 found a flux of $`3\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. In both cases the PSPC hardness ratio was $`0`$ which excluded a significant photoelectric absorption by cold gas.

Remillard et al (1993) identified H1419+480 with a broad-line type 1 AGN at $`z=0.072`$ ($`RA=14^h21^m29.60^s`$, $`DEC=+47^{}47^{}27^{\prime \prime }`$). Appendix A presents our own optical spectroscopic observations of H1419+480, which confirm that this source is a broad-line AGN with $`z=0.07229`$, as derived from the \[OIII\] emission doublet.

Ceballos & Barcons (1996) analyzed a sample of sources detected both in the MC-LASS sample and by $`ROSAT`$ (which included H1419+480). The large decrement of the flux from hard to soft X-rays is suggestive of heavy photoelectric absorption, which is at odds with the value reported for the PSPC hardness ratio. This was interpreted by Ceballos & Barcons (1996) as evidence for complex absorption, that could be modelled either in terms of an ionized absorber or in terms of a partial covering cold absorber.

In this paper we present XMM-Newton EPIC observations of H1419+480. We find that the source has changed the soft X-ray flux again and that its 2-10 keV flux is almost a factor 3 lower than the HEAO-1 one. A warm absorber is unambiguously seen in the XMM-Newton data, through mostly OVII absorption, but we argue that this warm absorber was weaker or absent in the previous pointed $`ROSAT`$ observation. An associated CIV absorption line is seen in an International Ultraviolet Explorer (IUE) spectrum of this source taken more than 20 years ago, but it is not seen in a similar observation taken one year later. We find that both X-ray and UV absorption can be explained in terms of the same amount of gas but varying ionising conditions.

## 2 X-ray properties

### 2.1 XMM-Newton observations

H1419+480 was observed by XMM-Newton (Jansen et al 2001) for 27 ks on the 27th of May of 2002, during revolution 451 (obsid=0094740201), as part of the Guaranteed Time programme of the Survey Science Centre. In this paper we analyze only the data obtained by the EPIC MOS (Turner et al 2001) and EPIC pn (Strüder et al 2001) cameras. The source is clearly detected by the RGS spectrographs, but it is far too faint to deliver any scientifically interesting results. The OM (Mason et al 2001) was used to obtain images in the filters U, B and UVW1.

All 3 EPIC cameras had the ’Thin 1’ filter on, and they were operated in partial window mode (MOS1 and MOS2) and small window mode (pn). Exposure times (excluding overheads) were 13 ks for the MOS cameras and 20 ks for the pn camera. However, most of the exposure time was lost due to high background flaring. After cleaning out these intervals, good-time intervals of $`8`$ ks were left for the MOS1 and MOS2 detectors, and 3 ks were left for the pn detector. Fortunately H1419+480 is bright enough (5 cts/s in the EPIC pn camera) to provide a large enough number of counts for spectral analysis.

Event files for the 3 EPIC detectors were taken from the distributed pipeline products, which were obtained by processing the observation data file with SAS v5.3.3. We filtered out high-background intervals, keeping only the most reliable single and double events and, in the case of EPIC pn, by removing those with low spectral quality. Calibration matrices were generated by using SAS v5.4.1, which we found to solve some problems at energies 0.5-1 keV with respect to the SAS v5.3.3 calibration.

X-ray spectra for H1419+480 and background, were extracted from the 3 EPIC cameras. The background spectrum was extracted from regions in the same chip as the source but not illuminated by H1419+480. This was also done in the MOS data, as opposed to source-free regions in the outer chips, to prevent vignetting affecting the background subtraction. The source plus background pn, MOS1 and MOS2 spectra were binned individually into bins containing at least 20 counts each. To fit spectral models to these data, xspec version 11.2.0 was used (Arnaud 1996). Bins outside the 0.2-12 keV bandpass were ignored. At the time of performing this analysis, a calibration problem (probably produced for an incorrect redistribution matrix) of the MOS detectors below 0.5 has been reported. In this paper we have taken the conservative approach of ignoring all MOS data below 0.5 keV.

### 2.2 The 2-12 keV X-ray spectrum

We first fit the 2-12 keV spectrum to determine the underlying continuum. A single power law fit gave $`\chi ^2=305.96`$ for 311 degrees of freedom (two parameters fitted), but with obvious residuals around 6 keV. Adding a redshifted gaussian emission line improved the fit substantially to $`\chi ^2=290.17`$ for 308 degrees of freedom. The F-test assigned a significance to the detection of this line of 99.90%. Parameters resulting from the fit are listed in table 1, while Fig. 1 shows the confidence contours in the line flux versus line dispersion parameter space. The line equivalent width derived from these parameters is $`240_{110}^{+125}\mathrm{eV}`$. All parameters are entirely consistent with those from other Seyfert 1 galaxies (Nandra & Pounds 1994). The EPIC pn spectrum shows a hint of a residual towards low energies, reminiscent of relativistic effects in a disk line, but since the MOS1 and MOS2 spectra do not show such resudials we do not explore this point any further. This emphasizes again the importance of calibration in the detection of weak features (such as relativistic line profiles) in moderate signal-to-noise data. In what follows the energy of the Fe line and its width have been fixed to their best-fit values (see Tab. 1) in order to prevent runaway solutions with unphyiscal parameters (e.g., extremely broad lines).

### 2.3 The 0.2-2 keV X-ray spectrum

We then extend the fit to 0.2 keV, including photoelectric absorption from the Galaxy ($`N_H=1.65\times 10^{20}\mathrm{cm}^2`$). The fit is not good, with $`\chi ^2=1140.45`$ for 795 degrees of freedom. There are obvious negative residuals as shown in fig. 2, where the 0.5-1 keV band has been excluded from the fitting region. We have explored a number of model fits in trying to understand this feature.

#### 2.3.1 Absorption plus soft excess

We first try to fit the spectrum with an absorbed power law plus a soft excess in trying to mimic the negative feature seen in the residuals. The soft excess is modelled as a steep power law (a black body model did not converge to a better value of the $`\chi ^2`$). The best fit yields a minimal reduced $`\chi _\nu ^2=936.03/792`$, corresponding to a $`\mathrm{\Gamma }_{\mathrm{soft}}=2.40\pm 0.09`$ and to an intrinsic absorption (at the redshift of the Seyfert 1 $`z=0.07229`$) of $`N_H=(5.7_{1.0}^{+1.2})\times 10^{21}\mathrm{cm}^2`$. Although the improvement in the $`\chi ^2`$ is very significant, strongly correlated negative residuals in the 0.5-1 keV band persist. One of the models suggested by Ceballos & Barcons (1996) to explain the unusual X-ray colours of this source on the basis of a partial covering cold absorber, yields an even worst $`\chi ^2=983.78`$.

#### 2.3.2 Single ion absorption

The absorption trough seen in the residuals is clearly reminiscent of an absorption edge or a UTA. The “equivalent width” of the absorption feature seen in the data is $`0.7`$ Å, which is twice larger than the values typically found in the Fe M UTAs. Besides, the absorption feature in H1419+480 spreads over $`5`$ Å, while the Fe M UTAs spread over $`1.5`$ Å in the RGS. The spectral resolution of EPIC at the energies of interest is $`0.61`$ Å, and therefore would be unable to spread a UTA up to 5Å. A further argument against the bulk of this feature being due to a Fe M UTA comes from the fit to a multi-ion photoinised absorber (see subsection 2.3.3) where all Fe M lines are included and yet the absorber’s parameters are very similar to those fitted by assigning the feature to (mostly) the OVII absorption edge. Therefore we conclude that, although UTAs can contribute to the detected absorption feature, they cannot account for all of it. We therefore start by attempting to interpret the negative spectral feature in terms of an absorption edge.

Indeed, adding a single absorption redshifted edge to the data results in a minimal $`\chi _\nu ^2=854.11/793`$ (i.e., substantially better than the soft excess model). The edge energy is $`E_{edge}=715_{19}^{+10}\mathrm{eV}`$ with a maximum depth of $`\tau _{edge}=0.58_{0.05}^{+0.07}`$. If interpreted as an intrinsic OVII absorption K-edge, at a rest frame energy of 739 eV, the absorber should be infalling towards H1419+480 at a velocity $`v=9100_{3800}^{+7200}\mathrm{km}\mathrm{s}^1`$. This velocity shift is significantly larger than the velocity outflow found in the best studied cases with the X-ray gratings.

One fact that might contribute, at least partly, to this shift is the presence of resonance absorption lines from the same ion, which are unresolved by EPIC, but still contributing to the shape of the absorption feature. This fact has been noted, among others, by Lee et al (2001) in the analysis of the Chandra data of the Seyfert 1 galaxy MCG-6-30-15. We therefore prefer to model the single ion absorption by the siabs model released in the PHOTOION package (Kinkhabwala et al 2003), under the assumption that the feature corresponds to OVII. The fit results in a $`\chi _\nu ^2=851.41/793`$ corresponding to a column density $`N_{OVII}=(2.1\pm 0.3)\times 10^{18}\mathrm{cm}^2`$ and an infall velocity of $`v=6000_{7350}^{+5400}`$ towards the source. Fig. 3 shows the confidence contours for these two parameters. We note that the fit is even marginally better than the one provided by the K-edge solely and that the velocity shift is now much more moderate (and consistent with zero), as in better-studied warm absorbers.

#### 2.3.3 Multiple-ion absorption

Although the single ion absorption fit does not leave obvious residuals in the X-ray data, we have attempted to add an OVIII absorber with the same recession velocity than the OVII absorber fixed at the best fit value of $`6000\mathrm{km}\mathrm{s}^1`$. The $`\chi _\nu ^2`$ further decreases to $`846.68/792`$ by just adding one parameter (the OVIII column density). The F-test statistic yields a significance of only 96.4% for this new component. Fig 4 shows the confidence contours in the OVII and OVIII column density parameter space. Fig. 5 shows the EPIC spectrum together with the 2 single ion absorption spectral fit, along with the featureless residuals.

We have also explored the possibility of describing the absorption feature in terms of a solar abundance photoinized gas, by using the model xiabs from the PHOTOION package (Kinkhabwala et al 2003). The best fit, by fixing the infall velocity at $`6000\mathrm{km}\mathrm{s}^1`$, is reached for $`\mathrm{log}\xi =1.30\pm 0.1`$ and a total H column density of $`N_H(5.1\pm 0.7)\times 10^{21}\mathrm{cm}^2`$ with a $`\chi _\nu ^2=871.32/793`$. Contour plots are shown in Fig. 6. However the fit is not as good as that for the OVII and OVIII absorption only, presumably due to unmatched features from other elements.

The xiabs model used does not include thermal (collisional) ionisation in the absorbing gas. We illustrate the effects of collisional ionisation by computing the values of $`N_{OVII}`$ and $`N_{OVIII}`$ for a photoionised gas at a temperature $`T=3\times 10^5\mathrm{K}`$ (typical of warm absorber models) and assuming a gas density of $`10^{10}\mathrm{cm}^3`$ (e.g., as in the Mathur et al 1994 models), although this last parameter is only marginally relevant as the gas is optically thin. We have used XSTAR (version 2.1) along with solar abundances and the results are shown in Fig. 4 for a few values of the ionisation parameter. In this case, some of the ionisation takes place by collisions and therefore the required ionisation parameter is marginally smaller ($`\mathrm{log}\xi 1.15`$). The required value for the column density is also around $`N_H5\times 10^{21}\mathrm{cm}^2`$. With this we conclude that both the ionisation parameter and the amount of gas along the line of sight are fairly well established.

### 2.4 Variability

Although the XMM-Newton exposure is relatively short, we have also explored whether the source varied within the observation. Counts have been accumulated in the 0.5-2 keV and 2-12 keV band. The reason for this resides in the fact that the 0.5-2 keV band is dominated mostly by the absorbed component and the 2-12 keV by the underlying continuum (plus Fe line).

Counts have been grouped into 100 s bins, which produced an EPIC pn count rate of 1.7 ct/s in the 0.5-2 keV band and 0.6 ct/s in the 2-12 keV band. Although it is a minor correction, the lightcurve of the background has been scaled and subtracted from the source’s lightcurve. Fig. 7 shows the resulting lightcurves of H1419+480. In order to look for the signficance of any variability, we fit a constant value to every lightcurve and compute the minimum $`\chi ^2`$. For the 0.5-2 keV band, $`\chi ^2=48.58`$ for 41 points, which yields a probability of the source being variable of only 83%. In the 2-12 keV band $`\chi ^2=60.9`$, with a probability of the source being variable of 98%. We have further checked whether any marginal variations are different in both bands, by fitting a constant to the difference between the 2-12 keV and 0.5-2 keV count rates. The $`\chi ^2`$ fit to a constant yields $`\chi ^2=32.0`$ for 40 points, which is entirely consistent with no X-ray colour variations. We therefore conclude that there is no significant variability of the source within the XMM-Newton observation and that in any case variations preserve the X-ray colours of the source. Therefore there are no constraints on the physical models for this source derived from variability.

### 2.5 Comparison to previous X-ray observations

Using the best-fit model with two single ion absorption (OVII and OVIII), we compute the flux of H1419+480, corrected for Galactic absorption, to be $`4.5\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the 0.5-2 keV band and $`7.3\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the 2-10 keV band. Adopting currently fashionable cosmological parameters ($`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$), the luminosity of H1419+480 is $`5.7\times 10^{43}\mathrm{erg}\mathrm{s}^1`$ in the 0.5-2 keV band and $`9.2\times 10^{43}\mathrm{erg}\mathrm{s}^1`$ in the 2-10 keV band.

The XMM-Newton 2-10 keV flux is almost a factor of 3 smaller than the HEAO-1 flux, a fact that can be at least in part explained by source confusion in the HEAO-1 collimators (see Barcons et al 2003 for a discussion on a similar situation for another source).

There are two $`ROSAT`$ observations of H1419+480 recorded from this source: one from the $`RASS`$ (reported in Schwope et al 2000) and a second one performed on January 1992 (observation sequence ROR700038), used in Barcons & Ceballos (1996). Table 2 lists some data of these observations, along with expectations from the XMM-Newton observation. We have folded our best-fit model to the XMM-Newton data through the $`ROSAT`$ PSPC-B response to compute the expected spectral shape in the $`ROSAT`$ observations and to convert count rates to fluxes.

The first obvious conclusion is that the 0.5-2 keV flux has changed significantly between the 3 different observations. The expected spectral shape in the $`ROSAT`$ observations appears similar (but not identical) to the measured one. The PSPC Hardness and Softness Ratios (see caption of table 2 for definitions) expected in the $`ROSAT`$ data from our best-fit model to the XMM-Newton spectrum have been computed with and without the OVII and OVIII absorption components. In particular the value measured for $`SR_1`$ in the 1992 observation appears consistent with the absence (or weakening) of the absorption features, but largely inconsistent with the presence of the absorption features. We conclude that the absorption features were much weaker or absent during the $`ROSAT`$ observation of H1419+480 in 1992.

## 3 An associated CIV absorber

The compelling evidence that X-ray warm absorbers share a common origin with associated narrow-line absorbers in the UV (Crenshaw et al 1999) urged us to search in available archives for UV spectroscopic observations of H1419+480. Unfortunately the Hubble Space Telescope (HST) has not observed this object with any of the UV spectrographs, but the IUE did observe it several times. Table 3 lists the observations conducted with the SWP short wavelength camera (operated at low dispersion $`6`$Å), which encompass both the Ly$`\alpha `$ $`\lambda 1216`$ and the CIV $`\lambda 1549`$ emission lines. We analyzed the INES (IUE Newly Extracted Spectra) data from both of these observations.

Fig. 8 shows the spectra obtained by IUE during both observations, around the CIV emission line. Only channels which have not been flagged for any reason are plotted. The overall flux shift between both spectra could be either true variability or a calibration problem. The most obvious discrepancy between both spectra occurs at around $`\lambda 1651`$Å, where the earlier observation (SWP17265) shows a dip that could be an absorption line.

To further investigate this, and keeping in mind that detecting absorption lines on top of emission lines is a difficult task, we have taken the 8 channels ranging from 1643 to 1656 Å in the SWP17265 spectrum. A linear function was fitted to that range, resulting in a $`\chi ^2=17.65`$ for 6 degrees of freedom. Multiplying by a gaussian absorption line and freezing its width to the resolution of the spectrograph (i.e., searching for an unresolved line) the $`\chi ^2`$ improved to $`\chi ^2=0.99`$ for 4 degrees of freedom. The F-test significance of that feature is 99.7% ($`3\sigma `$). The central wavelength is $`\lambda _{abs}=1651.1\pm 0.5`$Å, corresponding to an outflowing velocity of $`1800\pm 90\mathrm{km}\mathrm{s}^1`$. The equivalent width is also very uncertain, due to the poor sampling, but formally the best fit is $`3_{0.9}^{+0.7}`$Å.

We have also explored the wings of the Ly$`\alpha `$ $`\lambda 1216`$ and MgII $`\lambda 2800`$ (see fig. 9), but no other obvious absorption lines are evident. We conservatively estimate that a 3$`\sigma `$ upper limit for any Ly$`\alpha `$ absorption line present in this spectrum is about $`3`$Å.

At the 6Å resolution of this setup, the CIV doublet ($`\lambda \lambda 1548,1550`$) is unresolved and, obviously, we cannot measure directly the absorbing column density. Fig 10 shows the curve of growth (i.e., equivalent width versus column density) for the whole CIV doublet and a range of velocity dispersion parameters. As usual a Voigt profile has been assumed for each line in the doublet, with velocity width parameter $`b=\sqrt{2}\sigma `$ ($`\sigma `$ is the velocity dispersion of the gas, assumed Maxwellian), damping constant $`2.64\times 10^8\mathrm{s}^1`$ and a total oscillator strength of 0.28 (which is the sum of the oscillator strengths for the two lines). The approximation introduced by Whiting (1968) to the Voigt profile (accurate to better than 5%) has been used.

The amount of CIV absorbing gas is very uncertain, as the measured equivalent width falls in the saturated part of the curve of growth, but values around $`N_{CIV}10^{1516}\mathrm{cm}^2`$ would be appropriate for Doppler velocity parameters around $`100\mathrm{km}\mathrm{s}^1`$. In Fig. 11 we also plot the curve of growth for HI Ly$`\alpha `$, using an oscillator strength of 0.416 and a damping constant of $`4.7\times 10^8\mathrm{s}^1`$, along with the upper limit of $`3`$Å. We can see that relatively large HI column densities (up to $`10^{1718}\mathrm{cm}^2`$) could go undetected for similar values of the Doppler velocity parameter.

## 4 Discussion

Establishing a link between the X-ray and the UV absorbing gas is difficult in this case, as the observations are not simultaneous. We can, however, test whether or not it is plausible that the same or a similar amount of gas (perhaps in a different ionisation state) produced both absorption features.

The X-ray absorbing gas has a fairly well determined column of $`N_H5\times 10^{21}\mathrm{cm}^2`$ and it is highly ionized ($`\mathrm{log}\xi 1.11.3`$). Its velocity with respect to the emission redshift is only poorly defined by the EPIC data but consistent with zero. The absorber was significantly weaker or absent in the 1992 $`ROSAT`$ observation.

The UV gas detected as a CIV associated absorber in the 1982 IUE observation of H1419+480, has a fairly well determined outflowing velocity of $`1800\pm 90\mathrm{km}\mathrm{s}^1`$, but its CIV column density is only poorly constrained $`10^{15}N_{15}\mathrm{cm}^2`$ with $`N_{15}110`$. No HI Ly$`\alpha `$ of MgII absorption are found, although the IUE spectra at the Ly$`\alpha `$ region are noisier than around the CIV and MgII lines. The total column density implied by the CIV absorption is

$$N_H=2.4\times 10^{18}N_{15}f_{CIV}^1\mathrm{cm}^2$$
(1)
where $`f_{CIV}`$ is the fraction of C which is in $`CIV`$ ionisation state. Indeed this is a completely unknown parameter, as we cannot determine the ionisation parameter of this gas which applies to the 1982 observation from a single UV unresolved absorption feature. We must restrict ourselves to check whether a value of $`N_H5\times 10^{21}\mathrm{cm}^2`$, as in the X-ray absorber, can be obtained for reasonable ionisation conditions.

Using XSTAR (version 2.1) we have computed ionization fractions for CIV ($`f_{CIV}`$), HI ($`f_{HI}`$) and MgII ($`f_{MgII}`$) for a range of values of $`\mathrm{log}\xi 11.5`$, assuming the gas temperature fixed at $`T=3\times 10^5\mathrm{K}`$ and with a spectral shape given by a $`\mathrm{\Gamma }=1.84`$ power law. As pointed out by, e.g., Mathur et al (1994) this might result in a very inaccurate shape of the UV ionizing continuum and instead a more accurate shape should be used. However, for the purposes of this illustrative exercise, where the ionising parameter is essentially unknown, this rough modelling of the ionising spectrum is good enough.

For this range of parameters, $`f_{CIV}`$ ($`f_{HI}`$) decreases from $`7.1\times 10^3`$ ($`3.4\times 10^5`$) for $`\mathrm{log}\xi =1.0`$ to $`4.4\times 10^5`$ ($`4.0\times 10^6`$) for $`\mathrm{log}\xi =1.5`$. $`f_{MgII}`$ is zero for this range of ionisation. If we adopt the amount of gas seen in the X-ray absorber ($`N_H5\times 10^{21}\mathrm{cm}^2`$), HI Ly$`\alpha `$ would be difficult to detect in all the explored range of ionisations, as the HI column density would be below $`10^{17}\mathrm{cm}^2`$. However, $`N_{CIV}`$ would range from $`1.6\times 10^{16}\mathrm{cm}^2`$ for the lower ionisation states explored down to $`10^{14}\mathrm{cm}^2`$ in the opposite extreme. In conclusion, there is a range of ionisation parameters, in fact with values around those measured from the X-ray absorber ($`\mathrm{log}\xi 1.3`$), where it is plausible to observe the CIV absorption line and no HI Ly$`\alpha `$ or MgII absorptions. A slightly higher ionisation brings also the CIV absorption below detection limits.

## 5 Conclusions

The bright Seyfert 1 galaxy H1419+480 has been shown to be variable in X-rays by comparing our XMM-Newton observations with earlier $`ROSAT`$ data. A warm absorber is clearly seen mostly through OVII absorption, but also with $`2\sigma `$ evidence for OVIII absorption, in the XMM-Newton data at the emission redshift. The absorber was weaker or absent in a $`ROSAT`$ pointed observation performed in 1992. The absorber contains a gas column of $`N_H5\times 10^{21}\mathrm{cm}^2`$, and the ionisation parameter is around $`\mathrm{log}\xi 1.151.3`$ depending on whether or not ionisation by collisions is included.

An IUE observation of H1419+480 conducted in 1982 shows $`3\sigma `$ detection for a CIV associated absorber, outflowing with a velocity $`1800\pm 90\mathrm{km}\mathrm{s}^1`$. This confirms the close link between X-ray ionized absorbers and UV associated absorption reported, among others, by Crenshaw et al (1999). This absorption line was not present in an IUE observation obtained roughly 1 year later.

We find that a plausible common origin for all these variable absorption phenomena can be understood in terms of the same (or similar) amount of photoionised gas, but with varying ionisation state. Indeed, for values of the ionisation parameter similar to those inferred from the X-ray observations, we find that the X-ray absorbing gas could have produced the CIV absorption feature without any detectable HI Ly$`\alpha `$ or MgII lines. A slightly higher ionisation would explain the lack of the CIV absorption feature in the IUE observation in 1983.

We therefore conclude that H1419+480 has a variable ionised absorber. Both X-ray and associated UV absorbers are detected at different epochs in this object, along the lines of the sample studied by Crenshaw et al. (1999). Furthermore, we have shown that both absorbers are plausibly related.

## Appendix A The optical spectrum of H1419+480

H1419+480 was observed in the 4.2m William Herschel Telescope at the Observatorio del Roque de Los Muchachos in the island of La Palma (Canary Islands, Spain), on February 26, 1998. We used the ISIS double spectrograph with 600 line/mm gratings on both the blue and red arms, with the wavelengths centered at 5200 and 7000 Å respectively, in order to observe the H$`\beta `$+\[OIII\] region in the blue and the H$`\alpha `$+\[NII\]+\[SII\] in the red. Weather conditions were good and probably photometric, but the seeing was around 1.5 arcsec. Two 300 sec observations (co-added in the reduction process) were carried out with the slit aligned to paralactic angle. For details on the reduction process, see Barcons et al (2003) as the same setup and procedure was used.

Fig. 12 shows the optical blue and red spectrum of H1419+480. The \[OIII\] lines have some structure, with two peaks separated $`60\mathrm{km}\mathrm{s}^1`$. This is the main limiting factor in the precission of the redshift, that we measure by fitting both lines to a common redshift and find $`z=0.072296\pm 0.000004`$.

## Acknowledgements

The work reported herein is based partly on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA member states and the USA (NASA). It is also based on INES data from the IUE observatory. We thank Enrique Solano for assistance with the INES data and an anonymous referee for important suggestions. The $`WHT`$ telescope is operated by the Isaac Newton Group on the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. We acknowledge financial support by the Ministerio de Ciencia y Tecnología (Spain), under grants AYA2000-1690 and ESP2001-4537-PE.


format: text

text:
# XMM-Newton observation of SN1993J in M81

## 1 Introduction

M81 at the moderate distance of 3.6 Mpc Freedman94 (Freedman et al. 1994) is one of the best observed galaxies. When end of March 1993 the supernova SN1993J exploded in that galaxy, there was good hope to learn about the progenitor and follow, as with SN1987A in the LMC, the evolution of the early phases for years with high quality spectra in many wave bands. 5 days after the detection in the optical bands (Ripero93 (Ripero & Garcia 1993)), radio emission from SN 1993J was observed at 22.4 GHz Weiler93 (Weiler et al. 1993) by the VLA, demonstrating the interaction with circumstellar material (CMS), an important prerequisite for X-ray emission. 1 day later, at day 6, the SN was observed with the ROSAT satellite at soft X-ray energies around 1 keV (12 Angstrom). The immediate detection of a strong X-ray signal (Zimmermann et al.1993a, b) at these early times was unprecedented up to that date, even though in the following years a few of the 18 presently known X-ray supernovae111http://www.astro.psu.edu/$``$immler/supernovae\_list.html have been detected at even earlier epochs. The X-ray satellite ASCA with an extended spectral coverage Tanaka93 (Tanaka et al. 1993) and the hard X-ray instrument OSSE onboard GRO Leising94 (Leising et al. 1994) extended a few days later the spectral bands in which emission of SN1993J was observed.

X-ray observations of the supernova were continued with the ROSAT observatory up to 1998 and with the ASCA satellite up to 1996. After a break of 2 years the CHANDRA observatory observed SN1993J in May 2000 Swartz03 (Swartz et al. 2003) followed by the XMM-Newton observation in April 2001, reported here.

## 2 The Supernova Scenario

The fact that 93J began as a type II SN, characterized by prominent H lines in the optical spectra and few months later changed to the Ib type, characterized by no or few hydrogen but rich in He, is interpreted that the progenitor had lost all but a small amount of its hydrogen layers, thus that the shock heated effective photosphere could quickly sink through the thin H layer into the deeper He layers during the initial expansion and cooling phase.

In the case of SN1993J, optical plate material taken before the explosion, the optical light curve and theoretical considerations indicated, that the progenitor was probably a K0 I red supergiant with an initial main sequence mass of about $`15M_{}`$ , that had kept less than $`1M_{}`$ of H at the time of the explosion Nomoto93 (Nomoto et al. 1993)Podsiadlowsky93 (Podsiadlowsky et al. 1993, Ray et al. 1993)Bartunov94 (Bartunov et al. 1994)Utrobin94 (Utrobin 1994), Woosley94 (Woosley et al. 1994), Aldering94 (Aldering et al. 1994).

Loss of the outer hydrogen layers of the star may have been caused by either a strong stellar wind or by mass (Roche lobe) overflow to a companion star in a binary system. It is known that red supergiants blow slow winds, in the order of 10 km/s and carrying between $`10^6`$ to $`10^4M_{}`$ per year, into space. It is unclear whether a wind alone can account for the extreme loss of hydrogen in a star of $`15M_{}`$, therefore there is some tendency to prefer a binary scenario for SN1993J with mass overflow to a companion, while additional material escaped as a wind. In a binary scenario one might expect disk-like asymmetries, which indeed have been found in the profiles of some optical lines Matheson00 (Matheson et al. 2000). On the other hand radio images of the SN1993J explosion nebula Bietenholz01 (Bietenholz et al. 2001) show nearly perfect circularity with no indications for asymetries ($`<6\%`$).

In the collaps of the inner core of a supergiant star an extreme shockwave traverses the layers of the supergiant, driving the upper layers of the star into a circumstellar material producing strong shocks in that medium. Observations of the shock interaction in the radio, optical and X-ray bands provide inputs for understanding the explosion dynamics as well as element composition and mixing. Its physical basis is usually described in the standard or 2-shock model first worked out by Chevalier82 (Chevalier 1982).

In this model the expanding supernova drives a forward shock with a typical velocity of $`10^4`$ km/s into the circumstellar matter giving rise to temperatures in the order of $`10^9`$ K, $`100`$ keV. The circumstellar matter is assumed to originate from the stellar wind of the supergiant progenitor star. The resultant density profile is then a function of the mass loss rate by the wind, the wind velocity and has a radial dependance of $`r^s`$, where s=2 if the wind is isotropic and constant. In the model also a reverse shock is forming by the interaction of the shocked material with the SN ejecta, from which temperatures in the order of 1 keV are expected. Initially the reverse shock region may be radiative and can form cool material that may absorb the X-radiation of that region. The density profile of the SN ejecta is also described by an $`r^n`$ dependance. One can then derive a self similarity solution for that scenario where n and s, describing the density profiles of the ejecta and the circumstellar medium, play an important role.

The reason that early X-rays are so rarely seen, up to now only from 18 events, is mainly caused by the circumstellar environment. If the density is very high, then the initial radiation may not be able to ionize all the circumstellar material and X-rays may be fully absorbed. If, on the other hand, the circumstellar matter density near to the SN is very small then the intensity of the X-rays from the shocks may be too weak to be detected with present day instruments, as was the case with SN 1987A.

## 3 XMM-Newton observations

XMM-Newton was pointed towards the center of M81 - SN1993J is at 3 arcmin distance - on April 22/23, 2001, for a total of 132 ksec, from that data of 90 ks with the PN camera in small window mode and 83 ks with the MOS2 camera in image mode were taken. The medium filter was applied both in the PN and the MOS2 instrument setup during the observation. The MOS1 camera was set to timing mode and not used in this analysis.

Due to high and variable particle background in parts of the observing period, a temporal screening had to be applied. While for the MOS2 camera only the last few 1000 s were screened off, leaving 77 ksec of good data, in the PN observation about half of the observation was affected by bursts in the particle background. In order to maximize the signal to noise ratio in the PN spectrum, we determined the rates in intervals of about 1000 s separately for the source plus background field and the background field and eventually excluded all time intervals that deteriorated the signal to noise ratio in that band, leaving 69 ks of good data. The optimization relative to the cut radius of the source plus background field and the deselection of background affected time intervals was performed in different energy bands. The final selection of accepted intervals was done in the band above 1 keV to optimize the signal in the part of the spectrum where the lines clearly stand out in the spectrum. The cut radius chosen was 300 pixels, corresponding to 15 arcsec. The spectrum was binned such that each bin contained at least a fixed number of 5 (background subtracted) source counts.

## 4 The XMM-Newton spectrum

The PN spectrum of SN1993J shows emission lines of highly ionized Mg, Si, S, Ar, Ca and very cleary the complex of the Fe lines at high energies. Spectral fits were first performed using only the PN camera data offering the best statistics. Fitting was tested with different model components. Good fits were achieved over the whole energy band between 0.3 up to $``$ 11 keV, using a 2-component thermal model with variable element abundances (vmekal in the nomenclature of XSPEC). We also tried different shock models and found that a sedov model with variable abundances (vsedov in XSPEC) produces also acceptable $`\chi ^2`$ values, although all tested shock models tend to reproduce not so well the different line complexes visible in the spectra above 1 keV.

An overview on the best fit parameters with both the PN and the MOS2 camera and for different models is given in Tab. 1. The PN camera spectrum and the best fit 2-component thermal model (2vmekal) is shown in Fig. 1.

The spectral differences found between the PN and MOS2 fits can be understood as uncertainties in the spectral cross calibration of the two instruments (as stated in the document ’XMM-EPIC status of calibration and data analysis’, XMM-SOC-CAL-TN-0018 2003, valid for the software version SAS 5.4 applied in this analysis). In the following we therefore take the best fit results from the PN camera with its better statistics (7515 counts compared to 2593 from MOS2), especially at higher energies more representative for the overall behaviour, as the basis for the spectral discussion.

The Fe complex near to 6.7 keV is represented in the PN camera with about 33 counts (in the MOS camera about 5 counts). The width of the Fe complex in the PN spectrum (extending from 6.5 to 7.2 keV) is generally better fit with a higher ’high temperature’ component, while other line complexes, or better the regions between them tend to be better fit by somewhat lower temperatures. But the statistics in the spectrum do not reasonably suggest to introduce more components with additional temperatures.

Let us first look at the results from the 2-component thermal model fit (2vmekal). The low temperature component of $`kT_10.34\pm 0.04`$ keV dominates the spectrum up to energies of about 2 keV, but a higher temperature component is definitely needed to explain the spectrum at higher energies.

Given the fact that the high temperature component takes into account just the band between 3 keV and $``$8 keV, where the spectral fit anchors on the Ca IXX, XX and the prominent Fe XXIII, XXIV, XXV line complex, the temperature is likely to be not very well defined. Indeed the temperature of the high temperature component of $`kT_26.4\pm 4`$ keV is poorly restricted, because every change in the temperature component is compensated by changes in the abundances of Fe, Ca and Ar. Due to this dependence all the element abundances in the fits have errors typically much larger than the parameter values themselves.Therefore, the abundance values should be considered with some care.

One major difficulty with fitting the vsedov model is based in its enormous demands on computational power. It is therefore difficult to properly scan the parameter space. In the vsedov model the best fit temperatures are somewhat higher than those in the vmekal models. The fit requires additional intrinsic absorption. The lines at higher energies are not very well reproduced by the model.

## 5 Discussion

How do the X-ray observations fit into the scenario for SN1993J ?

### 5.1 Earlier X-ray spectra

We have reevaluated the ROSAT data in a homogeneous way Zimmermann97 (Zimmermann et al. 1997) and fitted the spectra with a thermal 1-component model (vmekal in XSPEC) where the element abundances were fixed to either solar or to the values obtained from a similar fit to the XMM-Newton PN spectrum. Due to the high temperatures in the early ROSAT observations there is almost no difference between fitting either solar or XMM-Newton based element abundances. In the November 1993 measurements the fit with the thermal model and solar abundances turned out to be poor but improves significantly (see Tab. 3) with the XMM-Newton based element abundances.

The first ROSAT PSPC spectra do only allow to set a lower limit of $`>17`$ keV for the temperature of a thermal bremsstrahlung, Raymond & Smith or mekal model. From the ASCA measurements a lower limit of $`>\mathrm{\hspace{0.17em}10}`$ keV was reported. Around the same time the hard X-ray instrument OSSE on the GRO satellite could determine a temperature of $`82_{29}^{+61}`$ keV at day 10 to 15 and, statistically already very weak, of $`74_{40}^{+120}`$ keV about 1 month after the event Leising94 (Leising et al. 1994). ROSAT measurements half a year later (Zimmermann et al. 1993d, 1994a) revealed a strong decrease from the initial high temperatures to temperatures around 1 keV. An ASCA measurement Uno02 (Uno et al. 2002) confirmed that strong temperature drop.

In terms of the standard 2-shock model it is suggested that the observed emission initially originated from the fast forward shock, while emission from the reverse shock region was blocked at that time by absorbing material due to fast cooling processes in the denser environment. It has been assumed that the initial absorption disappeared on a time scale of the order of 100 days so that half a year later the measured flux was then dominated by the emission from the reverse shock region Zimmermann94 (Zimmermann et al. 1994, Zimmermann et al. 1996, Fransson et al. 1996).

Unfortunately the X-ray lightcurve of SN1993J is not covered between days 50 and 200, so that the transition from the forward shock dominance to the reverse shock remains an open issue.

The Chandra spectrum of May 2000 was fitted by Swartz03 (Swartz et al. 2003) with a 3-component 2vmekal+mekal model. Their best fit model (kT’s of 0.35, 1.01, and 6.0 keV) does not fit the XMM-Newton PN data very well, especially in the energy regime below 0.7 keV, the low and high temperatures are almost identical.

### 5.2 The X-ray light curve

The lightcurve in Fig. 2 shows the development of the X-ray luminosity in the energy range 0.3 - 2.4 keV as determined from the 19 ROSAT observations, the Chandra observation from May 2000 and the new XMM-Newton data. The tendency over the first half year is characterized by a decline with $`t^{0.30}`$ (Zimmermann et al. 1993c, d), where t is the time since the outburst. Thereafter the lightcurve shows a bump and 5 years after the outburst the luminosity appears to return to the initial decline profile. The CHANDRA observation of 2000 and the XMM-Newton observation in 2001 are very close to the general $`t^{0.30}`$ decline.

The mass loss rate of the wind of the progenitor can be expressed as $`\dot{M}=4\pi \rho v_wr^s`$, where $`v_w`$ is the wind velocity. For a constant and homogeneous wind $`s`$ has a value of 2. Under the assumptions that the density $`\rho `$ of the circumstellar material at radius $`r`$ is dominated by that wind and the shock velocity stays constant, the observed luminosity is roughly proportional to the square of the density integrated over the emitting volume. In case of a constant and homogeneous wind the luminosity should show a time dependence of $`L_xt^1`$. The flatter time profile $`t^{0.30}`$ of the ROSAT luminosity observed during the first half year, corresponds for a spherically symmetric scenario to a density profile with s=1.65. In a more detailed consideration in the framework of the formalism of the standard model the development of the X-ray luminosity produced in the forward shock is expressed by: $`L_xT_e^{0.16}t^{(32s)(n3)/(ns)}`$ Fransson96 (Fransson et al. 1996), where $`r^n`$ describes the radial dependence of the density of the SN ejecta. For $`n>>10`$ it follows $`L_xt^{(2s3)}`$ and taking the observed slow decline with $`t^{0.30}`$ the value for $`s1.65`$ is confirmed, in good agreement with results published earlier by Immler01 (Immler et al. 2001).

In the picture of the standard model the shock front, while running outward, is probing increasingly older periods of the progenitor wind with a velocity roughly 1000 times the wind velocity. Thus in 5 years ROSAT scanned about 5000 years of wind history. In our simple formalism the measured time dependence of the luminosity requires that the ratio ’mass loss rate/wind velocity’ decreased during the last 5000 years in the life of the progenitor by more than a factor of 7. This scenario has been investigated in more detail by vandyk94 (Van Dyk 1994) using radio data and by Immler01 (Immler et al. 2001) using the ROSAT data to calculate absolute mass loss rates. The latter authors also speculate that the change of the wind properties might indicate a transition from a red to a blue supergiant phase where mass loss rates are at least an order of magnitude lower and wind velocities may reach 1000 km/s.

The bump in the X-ray light curve suggests a local increase in density above the general power-law profile of the circumstellar matter. It could be attributed to a change in the wind parameters of the progenitor or to an asymmetry caused by a possible binary scenario.

Interesting is a correlation with the expansion velocity of the SN shell that has been derived from high resolution radio images of SN1993J. In Fig. 6 of Bartel02 (Bartel et al. 2002) at about 350 days after the outburst the expansion rate of the radio image size shows a break and continues thereafter at a slower rate. Between day 369 and 1655, the leading edge and the plateau of the bump, the expansion velocity decreases in the order of 28%. The total mass which gave rise to the observed X-ray luminosity bracketing the bump is estimated to $`1M_{}`$. The total mass of the expanding shell covered up to day 350 was of the order of $`3to4M_{}`$. Both features, the radio deceleration break and the bump in the ROSAT light curve, independently indicate a flattening of the density profile in the circumstellar medium. Because the break in the expansion curve does not visibly affect the circularity of the radio images one may argue that the cause more likely can be attributed to a density change in the wind instead of an asymmetry by a binary scenario that would possibly affect the circularity of the images.

The closer inspection of the X-ray light curve shows a rapid decline between day 225 and day 370, followed by a rapid increase until day 935. A similar behaviour could have occured between day 1126 and day 1464, but there is just one data point in that period, the upper limit of which is consistent with no change. The rapid decline after day 1655 is very evident and recovery indicated by the Chandra data point at around day 2500 has occured. The existence of the latter two dips is statistically not too convincing, but they may be real as we show below. We suggest that each of the three dips indicates a density depression which is followed by a rise of the density. If the shock wave runs into a sufficiently low density regime no more matter will be shocked and heated, the matter heated so far will simply expand. If the expansion is isothermal the light curve will go down with t-2; if the expansion is adiabatic the decline will go like t-8/3, as long as the temperature is higher than about 2 keV. This temperature is about what we would expect at around day 350 (cf. Fig. 3). The falling dashed lines in figure 3 represent such expansion. After some time the shock wave encounters a jump in density to significantly higher values. The density increase raises the emission measure and the X-ray luminosity. Initially the luminosity increases linearly with t when it hits a step like density jump. The rising lines in figure 3 represent such density jumps. The intersection between the rising and falling lines denote the position or epoch of the density jump. We have compared these epochs with those where the highest, discontinuous deceleration of the shock wave has occured, which have been derived from the expansion rate of the radio images (cf. Fig. 6 of Bartel02 (Bartel et al. 2002)). There are also just three such events in the radio images. The two sets of epochs agree surprisingly well with each other to within less than 30 days. We conclude that these changes in the X-ray light curve reflect changes of the density in the ambient matter profile, which appears to show some repetitive pattern. Each of these density increases is preceded by some density decrease like in a wave. Whether this characterizes the activity of the progenitor star as far as the mass loss is concerned remains to be seen.

### 5.3 Spectral development

We can now combine the XMM-Newton temperatures with the ROSAT, ASCA and GRO results. Let us assume for a moment that the initial high temperature, as measured by GRO, originated from the forward shock, while half a year later the flux was dominated by emission from the reverse shock because of the higher electron densities in that region. One can use the similarity solution of the standard model to express in general terms the ratio between the temperatures in the forward and the reverse shocks that depends only on the parameters s and n, $`T_r/T_f=(3s)^2/(n3)^2`$ Fransson96 (Fransson et al. 1996).

With the wind density profile of s=1.65 derived from the light curve, and assuming that the 2 temperatures determined from the thermal 2-component fit to the XMM-Newton spectrum represent the temperatures in the forward and reverse shock regions, we obtain n=8.9 for the density profile of the ejecta. Taking into account the uncertainties of the temperatures as quoted for the PN 2vmekal model, then n results between 6.4 and 10.9.

Entering s and n into the formula that describes the temporal behaviour of the shocks, $`T_rT_ft^{[2(3s)/(ns)]}`$ Fransson96 (Fransson et al. 1996), we can extrapolate the XMM-Newton temperatures backward in time and compare them with the ROSAT, ASCA and GRO results (see Fig. 3). The temperature development of the forward shock component is roughly consistent with the temperatures deduced from the GRO data, and also the low temperature development, matching the ROSAT data at day 210, is consistent with its interpretation as coming from the reverse shock region.

But it has to be repeated that, while the parameter s appears reasonably well determined from the light curve (at least for certain periods), the value of n is derived under the assumption of thermal equilibrium between electrons and ions in the emission region. Estimates (see Fransson96 (Fransson et al. 1996)), however, show that even near to the time of the explosion equipartition is not expected in the wind shock.

One can not exclude that we see just emission from the forward shock. The hot material heated in the early days has cooled as it expands, both radiatively and adiabatically with the continuum temperature dropping faster than the higher ionization levels indicate, i.e. we have a delay of the recombination because of the long recombination time scales involved. Thermal equilibrium by just Coulomb collisions and ionization equilibrium is no more reached at times later than about a few days after outburst. Stellar wind matter heated at later stages as the forward shock proceeds is not in thermal equilibrium if not processes other than Coulomb collisions are at work. Therefore the electrons would stay relatively cool and also ionization equilibrium is far from being reached leaving the plasma underionized, which would explain the observation of the low temperature component.

With the sedov model we have tested such a scenario with emission from just one shocked region. But the relatively poor fitting results for the line complexes at higher energies show also that this model is not quite adequate.

## 6 Conclusions

The XMM-Newton spectrum is well fitted by a thermal 2-component model with ionization equilibrium. The predictions of the standard model describe reasonably well the observed development of the X-ray spectral temperatures, involving both forward and reverse shocks. The equally acceptable fit with a vsedov model shows that at this stage one cannot exclude that only one shock component is at work.

The strong interdependence between the temperature and the values for the abundances of the dominating elements, as far as the high temperature component is concerned, makes it very difficult to interpret the abundance results. But there is clearly no need to involve high values for the metal abundances, which would point to ejecta heated by the reverse shock.

The new XMM-Newton data point in the X-ray lightcurve gives arguments to the assumption of a general decline rate of $`L_xt^{0.30}`$ and it suggests that the bump in the light curve is probably due to an intermediate density fluctuation. The detailed inspection of the X-ray light curve shows that a density increase did not occur just once, but in some repetitive fashion, There are three dips in the X-ray light curve each one followed by an increase in luminosity. Each of the dips occurred very close to the times when the forward shock velocity underwent a dramatic deceleration as indicated by the expansion rate of the radio images.

The observations show that spatial changes of the powerlaw density profile have no significant impact on the morphology of the radio images, that remain very much circular. This indicates that the bumps in the X-ray lightcurve are more likely caused by density changes in a spherically symmetric wind rather than by asymmetries introduced by a postulated binary companion. For most of the time we need a wind that in the past left more material per volume unit in the circumstellar environment than near to the time of explosion.


text:
# Creation of twistless circles in a model of stellar pulsations

## 1 Introduction

During the past twenty years there has been increasing interest in the study of the coexisting regular and chaotic motions in Hamiltonian systems with few degrees of freedom. A basic problem of study was the transition to chaos in two degrees of freedom Hamiltonian systems which are slight perturbations of an integrable system, whose Hamiltonian is nondegenerate. In order to ease this task, most of the results on the dynamics of such Hamiltonian systems were obtained exploring the behaviour of a Poincaré map associated to a cross section in a fixed energy surface. A Poincaré map is an area preserving map (APM) defined on an annulus, that is, depending on $`(r,\theta )`$, where $`r`$ is a momentum–like coordinate and $`\theta `$ is an angular variable. Considering the lift of such a map, namely computing the angle without the restriction modulo $`2\pi `$, one associates to an orbit $`(r_n,\theta _n)`$ the rotation number which is the limit:

$$\rho =\underset{n\mathrm{}}{lim}\frac{\theta _n\theta _0}{2\pi n},$$
(1)
if it exists. If $`\rho `$ is an irrational number, then the orbit densely fills an invariant circle (KAM circle), while an orbit of rational rotation number $`p/q`$ — with $`p,q`$ coprime integers, $`q>0`$ — is a $`q`$–periodic orbit.

Typical questions of mathematical and physical interest include the persistence of the KAM circles after perturbation and determination of the threshold at which a circle of a given rotation number breaks up. The main results which answer these questions, such as the KAM theorem, the Poincaré–Birkhoff theorem, and the Moser twist theorem, are based on the validity of the twist property of the area preserving map transforming $`(r,\theta )(r^{},\theta ^{})`$, namely:

$$\frac{\theta ^{}}{r}0,r$$
(2)
The twist condition implies the monotonic change of the rotation number $`\rho `$ with $`r`$. The APM satisfying this condition is called monotonic twist map or simply a twist map. During the last twenty years the monotonic twist APMs were studied extensively . Consequently, their dynamics has gained almost complete understanding. In particular, this statement is a consequence of the late understanding that Hamiltonian systems have universal behaviour which can be revealed by studying explicit symplectic maps of the plane.

Recently, due to the increasing number of physical phenomena that are modeled by nontwist maps, this class of dynamical systems has finally captured the attention of the scientific community. The few studies on nontwist maps existent in the literature have revealed an unexpected variety of local and global bifurcations: dimerized island chains, periodic orbits collision, separatrix reconnection, vortex islands .

In the present paper we study the Poincaré map associated to a periodic perturbation, both in space and in time, of a linear Hamiltonian system. The dynamical system originates from a model of stellar pulsations previously introduced in  for low mass stars ($`M/M_{}<8`$, where $`M_{}`$ is the solar mass) and extended to intermediate mass stars ($`M/M_{}[8,11]`$) in . The latter work included some preliminary results concerning the nontwist character of the map. Here, we show that this map is an APM with an oscillating rotation number function. In spite of the system being a perturbation of a linear Hamiltonian, we argue that due to both the acquired nontwist character and the periodic character of the perturbation, the map presents the typical features of a generic class of nontwist APMs, such as reconnection and meandering, with the nongeneric scenario of creation of vortices. The nonmonotonic property due to the triplication bifurcation is shown to be different from that exhibited by the cubic Hénon map, which can be considered as the prototype of APMs, and which undergoes a triplication followed by the twistless bifurcation. Our study exploits the reversibility property of the initial system, which induces the time-reversal symmetry of the Poincaré map. The paper is organized as follows. In §2 we introduce our model of stellar pulsations. Before getting into the detailed analysis of it, we consider worth recalling the most important features of nontwist maps. This is done in §3. In this section we will also emphasize the consequences of both the triplication bifurcation and of the symmetry properties. Later, in §4 we study the dynamics of our simple model of stellar pulsations. Finally in section 5 we summarize our major findings and we draw our conclusions.

## 2 Description of the model

Pulsating stars are perfect tests of the theory of stellar evolution since the comparison between the theoretical pulsational characteristics (periods, amplitudes, and growth time scales) and the observed ones can help substantially in the fine tuning of the stellar models. Low- and intermediate-mass stars are prone to experience recurrent thermal instabilities and substantial mass loss during certain intervals of their evolution. Moreover, a substantial amount of pulsating stars in the Galaxy shows irregular behaviour and this makes them interesting targets in the field of dynamics systems. The present model embodies the essential of stellar oscillations by assuming an oscillatory driving originating in the stellar interior in the form of sinusoidal pressure waves. This model belongs to the family of the so-called one-zone models which treat the star as a rigid core surrounded by a dynamic, homogeneous gas shell. The pressure waves propagate through a transition zone characterized by a certain transmission coefficient until they hit the stellar envelope and dissipate. No back reaction of the outer layers on the inner ones is considered. The dissipation of the pressure waves induces a fluctuation in the radius and velocity of the outermost layer and this variability makes the object of our study.

The variation of the interior radius, $`R_\mathrm{c}`$, around an equilibrium value, $`R_0`$, is given by $`R_\mathrm{c}=R_0+\alpha R_0\mathrm{sin}\omega _\mathrm{c}\tau `$, where $`\tau `$ is the time in years and $`\alpha `$ and $`\omega _\mathrm{c}`$ are the fractional amplitude and the frequency of the driving, respectively. In absence of any driving force, the equation of motion is

$$\frac{d^2R}{d\tau ^2}=\frac{4\pi R^2}{m}P\frac{GM}{R^2},$$
(3)
where $`M`$ is the stellar mass and $`R`$ and $`P`$ are the radius and pressure in the shell of mass $`m`$, respectively. For convenience, we introduce the nondimensional variables $`xR/R_{}1`$ and $`t\omega _\mathrm{m}\tau `$, where $`R_{}`$ is the equilibrium stellar radius and

$$\omega _\mathrm{m}\sqrt{\frac{GM}{R_{}^3}}$$

is the characteristic frequency of the star.

Considering that the additional perturbative acceleration is proportional to the driving acceleration with a proportionality coefficient $`Q`$ we obtain

$$A=Q\frac{d^2R_c}{d\tau ^2}|_{\mathrm{ret}}$$
(4)
and after some algebraic manipulation  the final equation of motion reads

$$\begin{array}{ccc}\dot{x}\hfill & =\hfill & y\hfill \\ \dot{y}\hfill & =\hfill & xϵ\mathrm{sin}(\omega x\omega t\alpha \omega ^{1/3}\mathrm{sin}\omega t),\hfill \end{array}$$
(5)
where the characteristic frequency of the system, $`\omega `$ results from the use of the dimensionless time unit and from the assumption that $`R_0`$ encompasses almost the entire stellar mass. It is defined as

$$\omega \frac{\omega _\mathrm{c}}{\omega _\mathrm{m}}=\left(\frac{R_0}{R_{}}\right)^{3/2}.$$
(6)
In the action-angle coordinates $`(r,\theta )`$ resulting from the change $`x=\sqrt{2r}\mathrm{cos}\theta `$ and $`y=\sqrt{2r}\mathrm{sin}\theta `$, the Hamiltonian of the system is

$$H(r,\theta )=r\frac{ϵ}{\omega }\mathrm{cos}(\omega \sqrt{2r}\mathrm{cos}\theta \omega t\alpha \omega ^{1/3}\mathrm{sin}\omega t).$$
(7)
The parameters of the system are $`ϵ`$, the total driving amplitude, $`\alpha `$, the fractional driving amplitude and $`\omega `$, the characteristic frequency of the system. From the point of view of the parametric study, we take $`\omega `$ to be the most important parameter fixing the type of star to be studied and its evolutionary status. It is a measure of the ratio between the internal radius where the driving is produced and the equilibrium stellar radius. From the astrophysical perspective, we are interested in the dynamics of the system from Eq.(5) corresponding to $`\omega 3`$, more exactly $`\omega =3.0146`$. This low value of $`\omega `$ is associated with evolved stages of intermediate mass stars, as it results from stellar evolutionary models . Since we are interested in modeling the evolution of this type of stars, we fix $`\omega `$ during our simulations. We have performed a thorough parametric study in the space ($`ϵ`$,$`\alpha `$) as it assures a fine tuning of the strength of the perturbation. As in , we are interested in $`ϵ<100\%`$ and $`\alpha <40\%`$, which allow the perturbative approach.

## 3 Nontwist maps

The prototype of nontwist area preserving maps is the quadratic standard map :

$$\begin{array}{ccc}r_{n+1}\hfill & =\hfill & r_nk\mathrm{sin}\theta \hfill \\ \theta _{n+1}\hfill & =\hfill & \theta _n+2\pi \omega r_{n+1}^2(\text{mod}\mathrm{\hspace{0.17em}\hspace{0.17em}2}\pi ).\hfill \end{array}$$
(8)
The twist condition is violated along the curve $`r=k\mathrm{sin}\theta `$. For $`k=0`$ we get an integrable map, whose orbits lie on the circles of constant $`r`$. The rotation number function is $`\rho (r)=\omega r^2/(2\pi )`$. The circle $`r=0`$ has the maximum rotation number $`\rho _{\mathrm{max}}=\omega `$. It is called the twistless circle or shearless circle. A slight perturbation leads to the persistence of the twistless circle and nearby circles having a diophantine rotation number — that is, an irrational number badly approximated by rationals . Hence, the nontwist standard map, defined above, has an invariant circle of maximum rotation number among the rotation numbers of the nearby orbits. For a fixed perturbation parameter $`k`$ and for $`\omega `$ chosen in such a way that the map has a pair of Poincaré–Birkhoff chains containing $`p/q`$-type periodic orbits, as $`\omega `$ varies the two chains approach each other and at a threshold $`\omega =\omega _r(k)`$ the hyperbolic points of the two chains enter the twistless circle. This is called the reconnection threshold . At the reconnection threshold the two distinct hyperbolic orbits, of rotation number $`p/q`$, are connected by heteroclinic arcs.

There are many ways in which the monotonic twist condition can be violated. The vertical lines $`\theta =\theta _0`$ can be mapped into curves having a single extremum — quadratic twist  — or multiple extrema — cubic , quartic  or sinusoidal twist . Moreover, it has been proved  that whenever the elliptic fixed point of an APM of the plane passes through a triplication, a twistless bifurcation occurs or, equivalently, the rotation number as a function of the distance from the elliptic fixed point becomes a nonmonotonic function. The triplication of the elliptic fixed point of an APM $`f_\mu `$ of the plane, occurs at the value $`\mu \mu _3`$ at which the multipliers of the elliptic point cross the values $`\lambda =\mathrm{exp}(\pm 2\pi i/3)`$. At the triplication threshold $`\mu =\mu _3`$ an unstable period–3 orbit emerges from the elliptic fixed point in both directions, that is for $`\mu <\mu _3`$ and for $`\mu >\mu _3`$. As $`\mu `$ increases the elliptic point crosses the $`1/3`$ resonance and the triangular shape figure changes side after shrinking to elliptic fixed point (Figure 1).

The symmetry properties are very useful for the explanation of the behaviour exhibited by the system under study. We next recall some notions and results concerning time–reversal symmetry of a dynamical system . A system of differential equations

$$\dot{x}=F(x),x^n,$$
(9)
or equivalently, the vector field $`F`$ is called reversible if there is a smooth involution $`R:`$ ($`RR=\mathrm{Id}`$) such that $`DRF=FR`$, where $`DR`$ is the differential of the map $`R`$. $`R`$ is called the reversor of the vector field $`F`$. If $`\mathrm{\Phi }_t`$ is the flow of the vector field $`F`$, then the reversibility means that $`\mathrm{\Phi }_tR=R\mathrm{\Phi }_t`$. This is tantamount to say that the reflection of the trajectory $`x(t)=\mathrm{\Phi }_tx_0`$ is also an orbit of the system. A diffeomorphism $`f:^2^2`$ is called reversible with respect to a smooth involution $`R`$ of $`^2`$ if $`f^1=RfR`$. If $`f`$ is an $`R`$-reversible diffeomorphism then $`I=fR`$ is also an involution $`(II=\mathrm{Id})`$, and $`f=IR`$. This factorization of the map $`f`$ is very useful for the study of its dynamical properties. It can be shown that $`f^n`$ is a $`R`$-reversible map, too, $`n`$. Denote by $`I_n=f^nR`$, $`n`$, and by $`\mathrm{\Gamma }_n`$ the fixed point set of the involution $`I_n`$. $`\mathrm{\Gamma }_n`$ are called symmetry lines of the map $`f`$. A point in the intersection $`\mathrm{\Gamma }_j\mathrm{\Gamma }_k`$ is a periodic point of the map $`f`$, whose period divides $`|jk|`$. The symmetry lines $`\mathrm{\Gamma }_k`$ are transformed by $`f^n`$ into other symmetry lines in the following way:

$$f^n\mathrm{\Gamma }_k=\mathrm{\Gamma }_{2n+k}$$
(10)
In order to interpret the nontwist–type dynamics which will be encountered in our system when studying the triplication of the elliptic fixed point, we consider necessary to illustrate the twistless bifurcation in the case of the cubic Hénon map

$$\begin{array}{ccc}x_{n+1}\hfill & =\hfill & y_n+x_n^3+\mu x_n+0.7\hfill \\ y_{n+1}\hfill & =\hfill & x_n,\hfill \end{array}$$
(11)
which has an elliptic fixed point $`(x_0,y_0=x_0)`$ that can be easily found. At $`\mu _3=1.1505`$ this point undergoes a triplication, while at $`\mu _\mathrm{s}=1.12189`$ the saddle center collision occurs. In Figure 2, we show that as $`\mu `$ increases beyond $`\mu _3`$, the elliptic and hyperbolic period-3 orbits collide and disappear in a saddle–center bifurcation at $`\mu =\mu _\mathrm{s}`$.

In order to better identify the twistless bifurcation — the birth of a twistless circle — in Figure 3 we show the rotation number $`\rho (x,\mu )`$ as a function of the points $`x`$ lying on the symmetry line $`\mathrm{\Gamma }_0:=\{(x,y)|y=x\}`$ of the cubic Hénon map, which is a reversible map with respect to the reversor $`R(x,y)=(y,x)`$. More precisely, the rotation number is computed at the points $`x\mathrm{\Gamma }_0`$ and satisfying $`x<x_0`$. Figure 3 shows that between the triplication threshold and saddle–center collision threshold of the two period–3 orbits, the twistless circle is the circle of rotation number $`\rho _{\mathrm{max}}=1/3`$ (the maximum points from right to left correspond to $`\mu =1.135`$, $`\mu =1.13`$, and $`\mu =\mu _\mathrm{s}=1.12189`$). Observe that as $`\mu `$ increases up to $`\mu _\mathrm{s}`$, the twistless circle moves away from the elliptic fixed point.

Much more important (and difficult) is to detect the location of the twistless circle after the collision–annihilation of the two period–3 orbits. In order to get some insight into the manifestation of the twist propertiy beyond the threshold $`\mu _\mathrm{s}`$, we have computed numerically the rotation numbers for points starting near the elliptic fixed point of the cubic Hénon map and on the symmetry line $`\mathrm{\Gamma }_0`$. The rotation number functions $`\rho (x,\mu )`$ are plotted for $`\mu =\mu _\mathrm{s}`$, and some parameter values obtained varying $`\mu `$ in the direction given by the vector $`v=\mu _\mathrm{s}\mu _3`$ (Figure 4). Note that the cubic Hénon map has a twistless circle of maximum rotation number among the rotation numbers of nearby orbits. Varying the parameter $`\mu `$ as specified above, $`\rho _{\mathrm{max}}`$ decreases as a result of the collision–annihilation of the orbits , and at some value, $`\mu _\mathrm{t}`$, the rotation number becomes decreasing. At this point the map is no longer nontwist. This is the reason why $`\mu =\mu _\mathrm{t}`$ was called the threshold of the twistless bifurcation in  where the above presented scenario was first presented although in the reverse order.

## 4 Dynamical properties of the Poincaré map

Due to the the complexity of the perturbative term of the Hamiltonian of Eq.(7), the associated system in Eq.(5) is analytically intractable. In order to get some insight into its dynamical behaviour, we study the Poincaré map (in fact the stroboscopic map) associated to this time-periodic Hamiltonian system. The unperturbed Hamiltonian $`H_0(r,\theta )=r`$ is globally degenerate, which means that $`^2H_0/r^2=0,r`$. The dynamical consequence of this degeneracy is that all the orbits are periodic having the same period. This is in contrast to KAM theory, where the unperturbed system is nondegenerate. We will show that after a periodic perturbation (periodic in time and in space) the corresponding Poincaré map exhibits local and global bifurcations which in some cases are typical of the class of nontwist APMs which are perturbations of locally degenerate integrable APMs of an annulus, while in other cases are not.

### 4.1 The nontwist property of the map

Since the system given in Eq.(5) is time-periodic of period $`T=2\pi /\omega `$, we define the function $`s:𝕊^1`$, $`s(t)=t(\text{mod}\mathrm{\hspace{0.17em}\hspace{0.17em}2}\pi /\omega )`$, where $`𝕊^1`$ is the circle identified with the interval $`[0,2\pi /\omega )`$. Hence, we have the autonomous system of differential equations:

$$\begin{array}{ccc}\dot{x}\hfill & =\hfill & y\hfill \\ \dot{y}\hfill & =\hfill & xϵ\mathrm{sin}(\omega x\omega s\alpha \omega ^{1/3}\mathrm{sin}\omega s)\hfill \\ \dot{s}\hfill & =\hfill & 1.\hfill \end{array}$$
(12)
We denote by $`\mathrm{\Phi }_t^ϵ`$ its flow. $`\mathrm{\Phi }_t^ϵ`$ associates to each triplet $`(x_0,y_0,s_0)\times 𝕊^1`$ the position at the time moment $`t`$ of the orbit starting at $`t=0`$ from $`(x_0,y_0,s_0)`$. The plane

$$\mathrm{\Sigma }=\{(x,y,s)^2\times 𝕊^1|s=0\}^2$$
(13)
is transversal to the flow and the map $`P_ϵ:\mathrm{\Sigma }\mathrm{\Sigma }`$ defined by $`P_ϵ(x,y)=\mathrm{\Phi }_{2\pi /\omega }^ϵ(x,y,0)`$ is the associated Poincaré map.

It is a classical result in dynamical systems theory that the Poincaré map associated to a one degree of freedom periodic time–dependent Hamiltonian system is an APM. In order to interpret the dynamics displayed by the numerically computed Poincaré map, let us deduce some geometrical properties of this map. The vector field $`F`$, associated to the system from Eq.(12),

$$F(x,y,s)=(y,xϵ\mathrm{sin}(\omega x\omega sb\mathrm{sin}(\omega s)),1),$$
(14)
with ($`b=\alpha \omega ^{1/3}`$), is reversible with respect to the involution $`:^2\times 𝕊^1^2\times 𝕊^1`$, $`(x,y,s)=(x,y,2\pi /\omega s)`$, that is:

$$F=DF.$$
(15)
The reversibility property (15) is equivalent to $`\mathrm{\Phi }_t^ϵ=\mathrm{\Phi }_t^ϵ`$, $`t`$, which means that $`t`$ the diffeomorphism $`\mathrm{\Phi }_t^ϵ`$ is $``$–reversible, because $`\mathrm{\Phi }_t^ϵ=(\mathrm{\Phi }_t^ϵ)^1`$. The Poincaré map associated to the vector field $`F`$ is $`R`$–reversible, with respect to the involution $`R:\mathrm{\Sigma }\mathrm{\Sigma }`$, $`R(x,y)=(x,y)`$. Actually, notice that the reversor of the vector field $`F`$ is related to the involution $`R`$ by:

$$(x,y,s)=(R(x,y),2\pi /\omega s).$$
(16)
It is clear that the definition of reversible maps is fulfilled by $`P_ϵ`$:

$$\begin{array}{c}P_ϵ(R(x,y))=\mathrm{\Phi }_{2\pi /\omega }(x,y,0)=\mathrm{\Phi }_{2\pi /\omega }((x,y,2\pi /\omega ))=\hfill \\ \mathrm{\Phi }_{2\pi /\omega }(x,y,2\pi /\omega )=RP_ϵ^1(x,y).\hfill \end{array}$$
(17)
The symmetry line useful in the analysis of the dynamical behaviour of the Poincaré map is $`\mathrm{\Gamma }_0=\text{Fix}(R)`$, having the equation $`x=0`$.

The Poincaré map $`P_0`$ has the elliptic fixed point $`(0,0)`$. After a slight perturbation the fixed point persists as a symmetric elliptic fixed point of the reversible map $`P_ϵ`$ — it is a point $`(0,y_ϵ)\mathrm{\Gamma }_0\mathrm{\Gamma }_1`$. In order to reveal the nontwist property of the map $`P_ϵ`$ we compute and visualize the effect of the map on different semi-lines emanating from the elliptic fixed point (Figure 5). As it results from Eq.(1), the rotation number for the map $`P_0`$ is negative, more exactly $`\rho =1/\omega `$. Note that unlike the standard–like nontwist maps, here a line of constant $`\theta `$ is not mapped onto a parabola, that is a curve with a single extremum, but onto an oscillating curve, having many minima and maxima. This fact will lead to a different behaviour of our nontwist map, in comparison to the dynamical properties of nontwist standard–like mappings.

Besides this primary way to illustrate the nontwist property of the map $`P_ϵ`$, we computed the self–rotation number of the orbits of the map $`P_ϵ`$, starting on the symmetry line $`\mathrm{\Gamma }_0`$ and using the turning angle method . The self–rotation number for the map corresponding to the parameter $`ϵ`$ in the region of study is shown in Figure 6. Hence our map $`P_ϵ`$ is a nontwist APM, whose rotation number function $`\rho (y)`$ is oscillating, having more than one extremum, that is, more than one invariant circle whose rotation number is a local extremum. Note that for the fixed parameter $`\omega =3.0146`$ of the Hamiltonian system of Eq.(5), the associated Poincaré map $`P_ϵ`$ has an elliptic fixed point which undergoes triplication. In Figure 7 we represent the triplication curve $`\alpha ^{}(ϵ)`$.

In order to emphasize the significant differences in the phase portraits before and after the triplication point, we add that in Figure 1 we illustrated the triplication bifurcation of the elliptic fixed point of the Poincaré map for $`\omega =3.0146`$, $`ϵ=0.07`$ and $`\alpha =0.3`$ — panel (a)$`\alpha =0.05`$ — panel (b) — and $`\alpha =0.02`$ — panel (c). From the physical point of view it is only of interest the case in which $`\alpha `$ is constant and equal to $`0.3`$, and $`ϵ`$ varies. Thus, the considered pairs $`(ϵ,\alpha )`$ are located above the triplication curve. According to  for this case a twistless circle can exist. Comparing the rotation number associated to the cubic Hénon map (Figures 3 and 4) and our Poincaré map (Figure 6), we conclude that the Poincaré map has an oscillating rotation number with a larger amplitude near the elliptic fixed point. This could be interpreted as a superposition of the nontwist property induced by triplication on the nontwist behaviour due to the periodic perturbation of the initial Hamiltonian system. Near the elliptic fixed point, the Poincaré map has a negative twist. The apparent discontinuities of the rotation number are due to the distance chosen between the points on $`\mathrm{\Gamma }_0`$ whose rotation number is computed. While on the branch $`\mathrm{\Gamma }_0^+=\{(x,y)|x=0,y>0\}`$, a period–3 orbit has no point, it has one on the other branch, $`\mathrm{\Gamma }_0^{}=\{(x,y)|x=0,y<0\}`$. The orbit of a point on $`\mathrm{\Gamma }_0^+`$ encircling the homoclinic loops of the hyperbolic period–3 orbits can lead to an apparent discontinuity of the rotation number, although this function is continuous. The successive changes in the monotonicity of the map lead to a distinct nontwist behaviour in comparison with the known results reported for standard–like nontwist maps .

### 4.2 Birth of the dimerized island chains. Local and global bifurcations

In our parametric study, we first considered the cases in which $`ϵ`$ is small, in order to study in detail the departure of the system from the simple harmonic oscillator as this parameter is increased. We restricted the study of the Poincaré map to a rectangle limited by initial conditions close to physically reasonable values of the radius and velocity of the outer stellar layers. We focused on the creation of a typical phase-space structure of nontwist maps: a pair of dimerized island chains. This characteristic scenario of creation of new orbits in nontwist maps is illustrated in Figure 8. Panel a exhibits the phase portrait of the near–to–integrable map, where the elliptic fixed point is surrounded by invariant circles. At $`ϵ0.03`$ the Poincaré map has an invariant curve with three cusps (Figure 8b). Each cusp is a point of a new periodic orbit created through a saddle-center bifurcation for $`ϵ0.04`$: a stable and an unstable fixed point emanate from each cusp (Figure 8c). The homoclinic loop of every hyperbolic point surrounds the corresponding elliptic one. Between two successive hyperbolic points (in cyclic order) there is a heteroclinic connection. This separatrix structure called dimerized island chain appears in red in Figure 8c. In Figure 6d we have shown the self-rotation number $`\rho (0,y)`$ for the case whose Poincaré map appears in Figure 8c. In the former figure, the plateau $`\rho (0,y)=1/3`$ corresponding to $`y(1,2)`$ represents the constant rotation number due to the newly created elliptic fixed point on $`\mathrm{\Gamma }_0^+`$. The dimerized island chain is also responsable for the minimum of $`\rho (0,y)`$ around $`y2`$.

Generically, dimerized island chains of the same period exist in pairs and are born in stages. The second period–3 dimerized island chain is created through the same bifurcation process at a higher value of $`ϵ`$ (Figure 8d). In Figure 6f, the small plateau of $`\rho (0,y)`$ around $`y1`$ explains the birth of the second dimerized island chain. Between the two island chains the invariant circles are meanders, that is, the radius along such a circle is not a univoque function of the angle.

Around the last born period–3 dimerized island chain, a sequence of local and global bifurcations occurs. This is illustrated in Figure 9, where it can be seen that two independent orbits of the same period–3 are created by saddle-center bifurcation. They evolve in such a way that, finally, they interact with the orbits which belong to the first dimerized island chain. As $`ϵ`$ increases the newly born elliptic points approach the hyperbolic points of the dimerized chain. When $`ϵ`$ reaches a value of 0.11725 a global bifurcation occurs: the newly created orbits interfere, and the hyperbolic points of the dimerized island chain become hyperbolic points with homoclinic eight-like orbits encircling the newly created elliptic points.

In order to show the successive births of period–3 orbits as $`ϵ`$ increases with $`\alpha =0.3`$ fixed, we compute and visualize $`P_ϵ^k(\mathrm{\Gamma }_0^+)`$ and $`P_ϵ^k(\mathrm{\Gamma }_0^{})`$, $`k=1,2,3`$. This process was first illustrated in Figure 5a where $`\mathrm{\Gamma }_0^+`$ and $`P_ϵ^1(\mathrm{\Gamma }_0^+)=\mathrm{\Gamma }_2^+`$ were represented. We show in Figure 10 the connection between the symmetry properties of the map and the creation of periodic orbits. A saddle-center bifurcation occurs at the value of the parameter $`ϵ`$ at which $`\mathrm{\Gamma }_6^\pm `$ has a tangential contact with $`\mathrm{\Gamma }_0^\pm `$. The corresponding Poincaré map has for such an $`ϵ`$ an invariant curve with cusps which represent the points of tangency. Increasing the perturbation, the two symmetry lines intersect at two points, one being elliptic and the other, hyperbolic.

### 4.3 Creation of chains of vortices

A global bifurcation which is generic of nontwist maps is the so-called reconnection process, which is illustrated in Figure 11 for the nontwist multi-harmonic standard map :

$`x_{n+1}`$ $`=`$ $`x_n+2\pi \omega +y_{n+1}^2(\text{mod}\mathrm{\hspace{0.17em}\hspace{0.17em}2}\pi )`$ (18)
$`y_{n+1}`$ $`=`$ $`y_n+k\mathrm{sin}[x_n+\mathrm{arcsin}(e\mathrm{sin}x_n)],`$ (19)
with $`k=0.2`$, $`e=0.38`$ and different values of $`\omega `$. On the two sides of a twistless circle, two Poincaré–Birkhoff chains (necklaces of consecutive elliptic and hyperbolic points), having the same rotation number, approach each other (Figure 11a). At a given threshold, which is called the reconnection threshold, their hyperbolic points are connected by heteroclinic arcs (Figure 11b). Varying the parameter of the system further, two dimerized island chains emerge from the configuration created by the reconnection (Figure 11c). These chains are separated by meanders.

In our case, we witness the creation of the nongeneric vortices or dipoles . Figure 12 displays the birth of the first chain of vortices. As $`ϵ`$ increases, either $`\mathrm{\Gamma }_6^+`$ or $`\mathrm{\Gamma }_6^{}`$ intersects tangentially at a new point located at a radius larger than in the previous cases. As a consequence, new period–3 dimerized island chains are created after $`ϵ`$ crosses the value of tangential contact. However, each second dimerized island chain has a different creation mechanism. Unlike the previous case, the hyperbolic point of the first dimerized island chain bifurcates into two hyperbolic points in the transversal direction. In this process, a new elliptic point is born on the symmetry line. These hyperbolic points are connected by three heteroclinic arcs: one surrounding the previously existing elliptic point, another surrounding the new elliptic point, and the third one separating the two elliptic points. Thus, a pair of vortices is created.

As $`ϵ`$ increases, the process of formation of chains of vortices continues for larger radii (Figure 13). Note that if the elliptic orbits of one chain of vortices intersect $`\mathrm{\Gamma }_0^+`$, then the next one, which is created external to it, has a pair of elliptic points on the symmetry line $`\mathrm{\Gamma }_0^{}`$. The invariant rotational circle interpolating the hyperbolic points of the chain of vortices is the twistless circle. In Figure 13 it can be seen that such circles pass almost through the extremum of the symmetry lines $`\mathrm{\Gamma }_6^\pm `$. As already said, the formation of vortices is not generic and therefore the twistless circle that appears between them is not either.

### 4.4 Formation of the stochastic sea

Our system is a clear example of weak chaos, where the perturbation itself creates a separatrix network at a certain $`ϵ_0`$ (the dimerized island chains and the vortices chains) and then destroys it as $`ϵ`$ increases by producing regions of chaotic dynamics. The case of strong chaos implies that the unperturbed Hamiltonian intrinsically has separatrix structures and the perturbation just clothes them in thin stochastic layers. In both cases, the merging of the stochastic layers in the phase space may give rise to a single chaotic network called stochastic web. The larger $`ϵ`$ is, the wider is the stochastic web. Inside the cells of the web, there exists a set of islands of regular motion called invariant web-tori  whose dimensions are inversely proportional to the perturbation strength: for strong perturbations, they are engulfed by the stochastic sea. Our system presents this kind of behaviour, which is presented in Figure 14.

It has long been known (and not yet well understood) that the boundaries of the islands of regular motion can be “sticky”, meaning that chaotic orbits can spend a long time close to the islands and then escape in the chaotic sea. This phenomenon was observed already in  and was later called stickiness . In general, it can be explained by the transformation of the web-tori into web-cantori as the perturbation is increased. Contrary to the absolute barriers constituted by the invariant tori, the cantori are only temporary barriers. That is, the orbits close to the main islands can reach the chaotic sea when they encounter the gaps in the cantori, which happens after an indefinite time interval. Recent works  also suggest that the stickiness is due to the formation of higher-order resonant islands in the transition region between the main island and the chaotic sea (Figure 15). The presence of these so-called dynamical traps appears to be generic for Hamiltonian systems and therefore their existence critically determines the large-timescale behaviour of such systems. In this context, a behaviour extremely similar to sticky orbits has been shown to result from simulations of stellar variability using full hydrodynamical codes carried on significantly long time-intervals. In the absence of a better explanation, it was attributed to long-term (secular) nonlinear effects in the stellar envelopes. In a previous work , we have shown that our simple model recovers such a behaviour as it is generic of Hamiltonian systems. We also presented clear examples of sticky orbits together with an extensive discussion of their implications in the framework of the classification of variable stars. Therefore, we mention here only the fact that if a Hamiltonian approach is intrinsic to the phenomenon of stellar variability, we can expect that different categories of variable stars (classified according to their amplitude, frequency, irregularity) may be viewed under a more unifying perspective than it is presently considered. Moreover, considering that the sporadic excursions of the sticky orbits to the chaotic sea translate into high velocities of the outer layers, mass loss is very likely to occur, again in accordance with the observations.

Another important and distinctive feature of the phase portrait of nontwist maps is the existence of meanders, as it was mentioned before when discussing Figure 8d. Meanders are created between two successively born dimerized island chains or between two chains of vortices. In nontwist standard–like maps they become usual invariant curves after the reconnection of the two chains. For the sake of clarity, we illustrate next the relationship between the meanders and the reconnection process for our map. In Figure 16 we show the creation of meanders from the reconnection of the Poincaré-Birkhoff chains of period 34. In the case of the standard nontwist map , in the reconnection of Poincaré-Birkhoff chains of even periodic orbits, periodic points of the chains having the same stability type are aligned in phase. In other words, to an elliptic (hyperbolic) point of one chain corresponds below or above also an elliptic (hyperbolic) point. Due to this alignment, the reconnection is nongeneric and corresponds to a hyperbolic-hyperbolic collision leading to the formation of vortices. In our case, the two Poincare-Birkhoff chains of even period (period 34) have periodic points of opposite stability type which approach each other in their way to reconnection (Figure 16a). Therefore, the reconnection is generic and in the subsequent dimerized islands chains, hyperbolic-elliptic collision occurs. Moreover, after the reconnection, the twistless circle turns from a graph of a function of the angular variable to a meander (Figure 16b). This meander is slightly visible in Figure 14b.

Meanders appear to be robust invariant circles acting also as separatrices between distinct regions of chaotic dynamics. For stronger perturbations, their destruction allows the chaotic orbits to reach other stochastic regions previously separated by the meander. This behaviour was observed in nontwist standard-like maps , but until now there is no explanation for this robustness.

## 5 Conclusions

In the present work we have extended the preliminary results of Ref.  concerning the dynamics of a forced oscillator as a model of irregular stellar pulsations. The driving is characterized by two parameters, the fractional amplitude of the internal perturbation $`\alpha `$ and the total amplitude of the driving $`ϵ`$. Our aim was to undertake a more detailed analysis of the transition to the nontwist property in the associated Poincaré map. We have demonstrated that for a given $`ϵ>0`$ there exists an $`\alpha ^{}(ϵ)`$ that corresponds to a triplication of the central elliptic fixed point of the map. As we were interested in values of the parameters $`\alpha `$ and $`ϵ`$ above the triplication curve, the nontwist character of the map was undoubtly present .

Special attention was devoted to the process of formation of periodic orbits and to periodic-orbits collision. The reversibility of the Poincaré map and consequently, the symmetry properties allowed a clear identification of the structure of the periodic points. More precisely, the evolution and the intersections of the symmetry lines of the system provided good tracers of the steps of creation of dimerized island chains and vortex chains. We strongly argued in favor of these particularities of the map being entirely due to the superposition of two factors: the nonmonotonicity of the rotation number function induced by the triplication of the elliptic fixed point and the nonmonotonicity of the same function due to the oscillating perturbation. The latter feature admits alternative occurrence of maxima and minima, where the twist condition is violated. This translates into the minima and maxima of the symmetry lines which transform the invariant curves passing through them into twistless circles. Using these elements, we have followed the local and global bifurcations until the formation of stochastic layers around separatrices together with the associated sticky orbits. We have proved that the reconnection of even periodic orbits is generic by following the reconnection of two Poincaré-Birkhoff chains of period 34. From the astrophysical point of view, we have dwelled on the implications of the sticky regions and of the interplay between stochastic and regular regions in the framework of the analysis of stellar variability. In order to provide quantitative characterization of sticky orbits in the framework of stellar variability, we shall concentrate in the future on the causes and probabilistic estimates of the stickiness mechanism.

Acknowledgements. This work has been supported by the DGES grant PB98–1183–C03–02, by the MCYT grant AYA2000–1785, and by the CIRIT grant 1999SGR-00257. We also would like to acknowledge many helpful discussions with G.M. Zaslavsky.