Delphine Porquet
email: porquet@astro.u-strasbg.fr

Phil Uttley

James N Reeves

Alex Markowitz

Stefano Bianchi

Nicolas Grosso

Lance D Miller

Sandrine Deluit

Ian M George

A 100ks

Keywords: galaxies, Seyfert -galaxies, active -X-rays, galaxies -accretion discs -quasars, individual, ESO 113-G010

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The Seyfert ESO113-G010 (z=0.0257) was observed for the first time in X-rays, in December 1995, by ROSAT, revealing an X-ray luminosity of about 2.5 × 10 43 erg s -1 (assuming H 0 =75 km s -1 Mpc -1 and q 0 =0.5) and some indication of X-ray variability (Pietsch et al. 1998). An optical follow-up (2.2 m ESO/MPG telescope at la Silla) performed in November 1996 done by the same authors to characterize this bright galaxy led to the conclusion that this Seyfert is of type 1.8, and measured log L FIR =44.34 erg s -1 , log L B =43.97 erg s -1 , and log L X,ROSAT =43.40 erg s -1 . This object was observed for the first time above 2 keV with XMM-Newton on May 3, 2001 during a short exposure of about 4 ks. That observation revealed a soft excess below 0.7 keV and, more interestingly, a narrow Send offprint requests to: D. Porquet emission line at 5.4 keV (rest-frame), most probably originating from a redshifted iron Fe Kα line (Porquet et al. 2004b). The line was detected at 99% confidence, calculated via Monte Carlo simulations which fully account for the range of energies where a narrow iron line is likely to occur. Such narrow spectral features in the 5-6 keV energy range were discovered with XMM-Newton and Chandra in several other AGN: e.g. NGC 3516 (Turner et al. 2002), NGC 7314 (Yaqoob et al. 2003), Mkn 766 (Turner et al. 2004), AX J0447-0627 (Della Ceca et al. 2005). The redshifted Fe Kα line discovered in ESO 113-G010 had an EW of ∼ 265 eV, making it one of the strongest (in EW) redshifted iron lines observed to date in an AGN. The energy of the redshifted line could indicate emission from relativistic (0.17-0.23 c) ejected matter moving away from the observer, as proposed for Mrk 766 by Turner et al. (2004). Alternatively, if the classification as type 1.8 object is correct, then emission from a narrow annulus at the surface of the accretion disc is unlikely due to the very small inclination angle (i.e. less than 10 • ) required to explain the narrow, redshifted line in this intermediate Seyfert galaxy. However, emission from a small, localized hot-spot on the disc, occurring within a fraction of a complete disc orbit, could also explain the redshifted line (e.g., Turner et al. 2002;Dovčiak et al. 2004).

Here, we present a 100 ks XMM-Newton follow-up observation of ESO 113-G010. In section 2, we present the long XMM-Newton analysis performed in November 2005, and the corre- sponding data reduction used for this work. Section 3 presents the timing analysis: light curves, PSD analysis, rms-spectrum, and flux-flux analysis. Section 4 presents the spectral analysis, focusing on the study of the Fe Kα line complex between 6.4-7 keV. Finally, we summarize and discuss the main results of the present work (section 5).

XMM-Newton observation and data reduction

XMM-Newton observed ESO 113-G010 starting on November 10, 2005 (orbit 1085) with exposure times of ∼ 103 ks and ∼ 102 ks for the MOS, and PN respectively. The EPIC-MOS cameras (Turner et al. 2001) operated in the Large Window mode, while the EPIC-PN camera (Strüder et al. 2001) operated in the standard Full Frame Window mode. Background flaring was negligible during most of the observation and therefore the whole observation was used, except when mentioned in Section 3. The EPIC data were re-processed and cleaned using the XMM-Newton SAS version 6.5.0 (Science Analysis Software) package. The net resulting livetimes were 102.1 ks and 90.6 ks respectively for the MOS and PN cameras. Since pile-up was negligible, X-ray events corresponding to patterns 0-12 and 0-4 events (single and double pixels) were selected for the MOS and PN, respectively. Only good X-ray events (with FLAG=0) were included. The low-energy cutoff was set to 300 and 200 eV for MOS and PN, respectively. The source data were extracted using a circular region centered on the source position of diameter of 50 ′′1 , 25 ′′ , and 35 ′′ for MOS1, MOS2 and PN respectively. The diameter was smaller for MOS2 compared to MOS1 due to a dead pixel column about 26 ′′ from the source in the MOS2 data. ESO 113-G010 was by far the brightest X-ray source in this 30 ′ EPIC field-of-view. Background data were taken from box regions on the same CCD as the source (excluding X-ray point sources). The xspec v11.3 software package was used for spectral analysis of the background-subtracted spectrum using the response matrices and ancillary files derived from the SAS tasks rmfgen and arfgen.

Timing Analysis

PN light curves (time bin of 500 s, and error bars of 1σ) clearly show fairly rapid variability with relatively large amplitude, both in the soft (S: 0.2-2 keV) and hard (H: 2-10 keV) energy ranges (Fig. 1; upper and middle panels, respectively). The corresponding PN time-averaged count rate observed during May 2001 (duration of ∼ 4 ks) is reported (Fig. 1, dashed lines) and shows that the source was in a higher flux state for virtually all of the 2005 observation when compared to the 2001 observation. The PN average count rates in the 0.2-12 keV Fig. 1. Soft (S: 0.2-2 keV) and hard (H: 2-10 keV) 500-s resolution (background-subtracted) PN light curves of ESO 113-G010 (upper and middle panels, respectively), together with the hardness ratio between the two bands (lower panel). The grey stripes correspond to the 1σ error bars. For comparison the count rate and hardness ratio level found during the 4 ks May 2001 observation are reported (dashed lines). energy range are 2.42±0.02 cts s -1 and 3.71±0.01, for the 2001 and 2005 observations, respectively. The corresponding Hardness Ratio, HR (≡(H-S)/(H+S)) is also displayed (Fig. 1, lower panel). The error bars for the hardness ratio are obtained using Gaussian error propagation (see e.g., Pietsch et al. 2005). The 2005 HR values show that the source was on average in a softer state compared to May 2001. There was no very strong variation amplitude of the HR, except within the first 10 ks (with <∆HR>∼0.07).

To quantify this variability better, we carried out timing analysis using the power spectral density function (PSD), the rmsspectrum, and flux-flux analysis. To maximize the signal-tonoise we combined events from all three instruments, re-scaling the background regions used for the MOS detectors to match the source/background area ratio of the PN. We used only data from times between 8 ks and 95 ks from the start of the observation, to remove the times of non-negligible background flaring for secure timing analysis, and merged the good time intervals of the different instruments in order to use only data received simultaneously from all three instruments. The combined events file was used to extract light curves in different bands to make the PSDs and rms-spectra. Although the instrument responses differ, in all cases (including rms-spectra) we normalize by the mean flux in each band, so the effect of the spectral response is taken out.

PSD Analysis

We first consider the PSD for the light curve in the 0.3-10 keV energy range, which is shown in Fig. 2. The PSD represents the average of the (squared) amplitude of variations as a function of temporal frequency. The data are then fitted with a simple power law plus the flat level of variability due to Poisson noise expected given the error bars. A simple power-law plus noise is not a very good fit (χ 2 = 24.1 for 11 d.o.f.), and is rejected at > 98% confidence (for power-law slope α = -2.1). Allowing the noise level to be free improves the fit to χ 2 = 18.9 for 10 d.o.f., which is still rejected at > 95% confidence. Finally, leaving the noise level free while allowing a break in the PSD improves the fit: fixing the low-frequency slope to -1 (similar to the low-frequency PSD shapes observed in other AGN, e.g. McHardy et al. 2005), we find a break frequency of (5.9 ± 1.7) × 10 -4 Hz (errors correspond to 90% confidence for a single interesting parameter), with an upper limit on high-frequency slope of α < -2.9. The corresponding χ 2 = 11.9 for 9 d.o.f., which, according to the F-test, is a significant improvement over the unbroken power-law model at the 95% level. We also find a good fit (χ 2 = 13.9 for 10 d.o.f.) for an exponentially cut-off powerlaw, with low-frequency slope fixed to -1 and cut-off frequency (3.7 +1.0 -1.7 ) × 10 -4 Hz, which we show in Fig. 3. We caution that the possible break or cut-off in the PSD is only marginally significant, but note that, in support of these models, the product of frequency and power, which is constant on the slope=-1 part of the PSD is ∼0.01, which is of the same order as that observed in other AGN. In comparison with other AGN, the measured cut-off frequency is comparable to those measured in some other rapidly-variable Seyferts, such as MCG-6-30-15 (8 × 10 -5 Hz, McHardy et al. 2005; Vaughan et al. 2003b), and NGC 4051 (8 × 10 -4 Hz, McHardy et al. 2004). The correlation between AGN X-ray variability PSD break time-scale and black hole mass is now wellestablished. It appears that AGN PSD break time-scales can be linearly scaled down by the black hole mass to match the breaks observed in the PSDs of stellar mass black holes in X-ray binary systems (e.g. Markowitz et al. 2003;Uttley & McHardy 2005). Recently, McHardy et al. (2006) showed that the residual scatter in the mass-time-scale relation could be explained if timescales also scale with luminosity or equivalently, accretion rate. Combining AGN and stellar mass black hole data for Cyg X-1, McHardy et al. framed the best-fitting relation for break timescale logT B = 2.1logM BH -0.98logL bol -2.33, with T B in units of days and M BH and bolometric luminosity L bol in units of 10 6 M ⊙ and 10 44 erg s -1 respectively. The best-fitting PSD break timescale we measure, 3.7 ×10 -4 Hz, corresponds to T B = 0.03 days. For a range in L bol from 10 44.5 -10 45 erg s -1 (i.e. assuming a minimum equal to the sum of FIR, optical and X-ray luminosities reported by Pietsch et al. 1998, and in this work), we infer from the mass-luminosity-time-scale relation of McHardy et al. (2006) that M BH ranges from 4 × 10 6 -10 7 M ⊙ and the source is accreting at or close to the Eddington rate (assuming even higher L bol implies that the AGN is a super-Eddington accretor).

We also measured PSDs using smaller subsets of the observed energy range, but due to the relatively weaker S /N of these PSDs, the energy dependence of the PSDs could not be well constrained.

The rms-spectrum

The rms-spectrum quantifies the spectral shape of the variable part of the emission (see Vaughan et al. 2003a for a review). Specifically, one measures the root-mean-square (rms) variability in light curves made from a number of narrow energy bands or 'channels', and then normalizes each rms value by the mean in that channel to obtain the fractional rms (also known as F var ). A plot of F var versus energy yields the shape of the variable part of the spectrum relative to the mean spectrum (i.e. independent of instrumental spectral response). The error bars are determined using the formula of Vaughan et al. (2003a). We show the measured rms-spectra in Fig. 4, obtained using light curves binned to 5 ks (top panel) and 500 s (bottom panel). The rms-spectra obtained with 500 s binning naturally includes faster fluctuations than for the 5 ks binned data, and the F var measurements are accordingly larger. F var appears to decrease above about 2 keV for both bin timescales. Therefore the soft X-ray excess down to 0.3 keV appears more variable compared to the hard energy range. This is different from the typical rms-spectra consistent with a broad peak in the variability amplitude between 0.8 and 2 keV, and a much lower variability amplitude below 0.8 keV (e.g., Fabian et al. 2002;Inoue & Matsumoto 2003;Ponti et al. 2006). In addition, there is a hint of a dip in the 500 s bin rms-spectrum around the Fe Kα complex energy (∼6.4-7 keV). This could be due to a reflected component remaining constant while the power-law continuum varies, as proposed for the twocomponent model ( §3.3). An origin in a distant material is unlikely since no narrow 6.4 keV Fe Kα line is observed (cf. § 4).

Flux-flux analysis

Many AGN show a characteristic steepening of their continuum X-ray spectra at higher fluxes. A useful insight into the origin of this continuum spectral variability can be obtained using fluxflux analysis. Taylor et al. (2003) showed how a simple plot of soft X-ray flux versus a harder X-ray flux can indicate whether spectral variability is due to pivoting at high energies (where the flux-flux plot takes a power-law form), or due to changes in the strength of a continuum component with a constant spectral shape plus a constant-flux, constant spectral shape harder component (where the flux-flux plot is linear with a positive intercept on the hard flux axis). Here we carry out a flux-flux analysis of the ESO 113-G010 data, using the same combined MOS and PN events data used for the rms-spectral analysis.

We first realize a flux-flux plot to compare the 0.3-1 keV and 3-10 keV ranges. We measured the fluxes in simultaneous 500 s bins, and following the approach of Taylor et al. (2003), we bin the 3-10 keV flux as a function of the 0.3-1 keV flux, obtaining the standard error in the mean 3-10 keV flux of each bin from the spread in fluxes in the bin. The resulting flux-flux plot is shown in Fig. 5. A linear plus constant model provides a good fit to the data, with χ 2 = 5.4 for 6 d.o.f., for positive constant offset on the 3-10 keV flux axis of 0.082±0.014 count s -1 (error bars are 1σ). A single power-law (with no offset) also provides a good fit to the flux-flux plot (χ 2 = 5.4 for 7 d.o.f.), for a power-law index of 0.75 ± 0.04. Thus, the existing data cannot distinguish between spectral pivoting of the continuum and a two-component origin for the spectral softening, primarily because the data span an insufficiently broad range of fluxes to see the characteristic change in gradient expected from the power-law.

Although we cannot distinguish between spectral pivoting and a two-component spectral model, we can use the flux-flux method to determine the shape of any constant component under the assumption that the two-component model is correct. Following the approach of Taylor et al. (2003), we fit linear plus constant models to the flux-flux plots measured from light curves obtained in a number of narrow bands. For the x-axis, we use the 0.3-0.5 keV band as our 'reference' band. Thus all offsets measured in other bands are relative to this band only, under the assumption that there is zero flux in the constant component in the reference band. If there is a non-negligible flux in the reference band, the implied offsets in other bands increase accordingly. Thus, if the two-component model is correct the offsets we measure from flux-flux plots here represent a lower-limit on the strength of any constant component. To remove the effects of the instrument spectral response, we normalize each measured offset by the mean flux in each band. The resulting spectrum of A power-law has been fitted to the 2-5 keV data and extrapolated to lower and higher energies. A significant soft X-ray positive residual is seen below 1 keV, as well as the presence of a Fe Kα complex near 6-7 keV. For presentation, the data have been re-binned into groups of 15 bins for PN (10 for MOS), after group of a minimum of 100 and 50 counts per bin were used for the fit for the PN and MOS respectively. Left panel: 0.2-12 keV PN spectrum. Right panel: co-added MOS 0.3-10 keV spectrum.

fractional offsets is shown in Fig. 6. Since fractional offset increases significantly toward higher energies, the offset spectrum is harder than the mean spectrum, implying that any constant component is also significantly harder than the variable component, as expected if we are witnessing a constant reflection component. Notably the reflected flux contributes about 40% of the mean at the hardest energies, similar to what is observed in MCG-6-30-15 using the same approach (Vaughan & Fabian 2004).

Spectral analysis

In this section, we present the spectral analysis of ESO 113-G010. In this article, we focus on the analysis of the 2-12 keV hard X-ray range (not affected by WA signatures), including the study of both the continuum and the Fe Kα line complex between 6.4-7 keV. We co-added the MOS1 and MOS2 data into a single spectral file to maximize the signal-to-noise ratio. The time-averaged PN and MOS spectra are binned to give a minimum of 100 and 50 counts per bin, respectively; while the subspectra used in the time-resolved spectral analysis (see below) are binned to give a minimum of 20 counts per bin. In all subsequent fits (above 2 keV), we include the Galactic column density (N Gal H =2.74 ×10 20 cm -2 ), obtained from the coldens program using the compilations of Dickey & Lockman (1990). Note that all fit parameters are given in the source's rest frame, with values of H 0 =75 km s -1 Mpc -1 , and q 0 =0.5 assumed throughout. The errors quoted correspond to 90% confidence ranges for one interesting parameter (∆χ 2 =2.71). Abundances are those of Anders & Grevesse (1989). In the following, we use the updated cross-sections for X-ray absorption by the interstellar medium (tbabs in XSPEC) from Wilms, Allen, & McCray (2000).

The main spectral components

As a first step, we fit an absorbed power-law model over the 2-5 keV energy range where the spectrum should be relatively unaffected by the presence of a broad soft excess, a warm absorber-emitter medium, an Fe Kα emission line, and a contribution above 8 keV from a high energy Compton reflection hump. In this energy range, the PN data are well fitted by a single power-law model with Γ= 2.02±0.04 (χ 2 /d.o.f.=228.3/214). For comparison the value from the 2001 observation is 1.96±0.22. This power-law index is consistent with values found for broadline radio-quiet quasars: Γ=1.90 with a standard deviation of 0.27 (Porquet et al. 2004a). The 2-10 keV time-averaged luminosity (inferred from the parameter fit reported in Table 1) is 4.0×10 42 erg s -1 . Figure 7 displays the 2-5 keV PN spectrum extrapolated over the 0.2-12 keV broad band energy. A strong positive residual is seen below 1 keV due to the presence of a soft X-ray excess; and in addition, there is a positive residual in the 6-7 keV energy range (see § 4.2). For comparison, we also report the co-added MOS 1 and MOS 2 spectrum, which confirms the presence of a complex soft excess as well as a positive deviation near the 6-7 keV. We henceforth focus only on the PN data, since the PN CCD has a better sensitivity over a broader energy range (0.2-12 keV) compared to the MOS CCD. However we have checked that all the spectral fits are consistent with those obtained with the MOS, albeit with much lower photon statistics. We would like to note that the RGS data show the presence of several narrow blue-shifted absorption lines that reveals the presence of a warm absorber (WA) in outflow. The presence of this medium has to be taken into account in order to fit properly the overall 0.2-12 keV PN spectrum. Spectral analysis combining RGS and PN data will be presented in a forthcoming paper.

The Fe Kα line complex

Time-averaged analysis

The positive residuals between 6.5-7 keV in the quasar frame can be most probably associated with a Fe Kα complex as observed in most AGN (Porquet et al. 2004a;Piconcelli et al. 2005;Jiménez-Bailón et al. 2005). The mean iron line profile as a ratio of the continuum above 2 keV is displayed in Fig. 8, and shows a double-peaked line profile. Therefore, we fit the 2-12 keV energy range with an absorbed power-law continuum and two narrow Gaussian line profiles (Table 1). We found a much better fit, with ∆χ 2 =33.2 for 4 additional parameters (F-test > 99.9%). The line energies are 6.50 +0.06 -0.05 keV and 7.00±0.04 keV. The F-test probabilities for each line are 99.7% and 99.98% for Table 1. Best-fitting spectral parameters of PN time-averaged spectrum in the 2-12 keV energy range with an absorbed (Galactic, N H =2.74×10 20 cm -2 ) power-law (PL) component plus a line profile: zgauss: Gaussian profile; and diskline: profile line emitted by a relativistic accretion disc for a non-rotating black hole (Fabian et al. 1989) assuming an emissivity law q equal to -2. (a): R out =1 000 R g . (b): R out =1.2 × R in . The line fluxes are expressed in units of 10 -6 ph cm -2 s -1 . F-test probabilities are calculated with the corresponding power-law model as reference. Fixed parameters are indicated by (f). Fig. 8. Fe K line profile (quasar frame) when PN time-averaged data are fit with a power-law continuum. For presentation, the data have been re-binned into groups of 10 bins for PN, after group of a minimum of 100 counts per bin were used for the fit. The grey vertical line corresponds to 6.4 keV. the 6.50 keV and 7.0 keV lines, respectively. Allowing the line widths to vary does not significantly increase the goodness of fit, with ∆χ 2 = 0.5 for two additional parameters. The former energy is compatible within the error bars with a line emitted by moderately ionized iron (∼Fe xix-Fe xxii). However this line is more likely a blue or red peak of a relativistic line, indeed a genuine Fe xix-Fe xxii (L-ions) emission line complex would suffer from resonant Auger destruction (Ross et al. 1996;[START_REF] Liedahl | Diagnostics of Astrophysical Plasmas: Theory, Experiment, and Observation[END_REF]). The latter line energy corresponds to highly ionized iron line (Fe xxvi) or to a highly blue-shifted iron line from lower ionization matter. However there is no evidence for narrow emission from neutral iron, the 90% confidence upper-limit to a narrow (σ=10 eV) iron line at 6.4 keV (a 6.5 keV narrow line with a width of 10 eV is also fitted to take into account its emission contribution) is only 32 eV. If the 6.5 keV line contribution is not taken into account, we would infer EW(6.4 keV) ≤ 47eV. From the branching ratio of Fe Kβ/Fe Kα, expected to be less than 0.145 (Palmeri et al. 2003, and references therein), the maximum contribution of the Fe Kβ (from cold iron) to the 7.00 keV line is only up to 5 eV.

We also tested if these two features can be explained by a double-peaked relativistic line profile emitted by a relativistic accretion disc for a non-rotating black hole (diskline; Fabian et al. 1989). The improvement of the fit is also significant (F-test > 99.9%) and the line energy is about 6.8 keV, i.e., corresponding to a highly ionized iron ion (Fe xxv-Fe xxvi). Then, we checked whether the double-peaked like profile can be explained by emission from a disc annulus (using the diskline profile with the width fixed to 20% of R in , i.e. R out =1.2 × R in ). We fixed the disc inclination to 45 • (representative of an intermediate type Seyfert). We found a very good fit to the data with a line energy of 6.79±0.05 keV and an inner annulus radius of 191 +71 -93 R g , indicating that the line is not emitted in the inner part of the disc. In the case where the disc inclination is not fixed to 45 degrees, then the line solution is degenerate, whereby the inner annulus radius increases with the disc inclination, while the line energy decreases.

Monte Carlo simulations were done to assess the significance of detecting the emission complex in this time-averaged spectrum. These simulations, described in detail in the Appendix, demonstrate that regardless of the model used, the emission complex is indeed significantly detected.

We also checked for an upper limit to a line at 5.39 keV, since such a line was found in the 2001 XMM-Newton observation of May 2001. We fixed the energy of the line to 5.39 keV (see below), and the line width to 0.1 keV (as done in Porquet et al. 2004b). We find EW≤22 eV (at 90% confidence level for one interesting parameter).

For comparison, we have re-processed the XMM-Newton May 2001 observation with the same version of the SAS (SAS 6.5.0) used for the 2005 observation. We fit the data with an absorbed (Galactic column density) power-law in the 2-10 keV energy range and four Gaussian lines in order to take into account the emission contribution of the lines near/at 5.4 keV, 6.4 keV, 6.5 keV, and 7.0 keV. The three latter line energies are fixed as well as their widths to 10 eV. Only the energy of the first Gaussian line is let free to vary with the width fixed to 0.1 keV (as done in Porquet et al. 2004b). One should notice that in Porquet et al. (2004b), only the line at 5.4 keV was fitted. We found for the four lines: EW 1 =246±147 eV (E 1 =5.39±0.09 keV), EW 2 ≤190 eV (E 2 =6.4 keV), EW 3 ≤100 eV (E 2 =6.5 keV), and EW 4 ≤208 eV (E 2 =7.0 keV). The EW for the line at 6.4 keV, 6.5 keV, and 7.0 keV are compatible with the values found during the 2005 observation. For the line at ∼5.39 keV, the EW are not compatible and would mean that if the line detected during the 2001 observation is real, it is transient.

Time-resolved analysis

To test for possible rapid shifts in the energy and/or flux of the Fe Kα emission complex, we have created X-ray intensity maps in the energy-time plane using the PN data. Photons from the source cell were accumulated in pixels in the energy-time plane. The pixel distribution was smoothed in energy by the instrumental resolution, using a Gaussian of 140 eV (appropriate for the single plus double events with the latest calibration), and smoothed in time using a top-hat function of width 20 ks. Each time-slice was background-corrected by subtracting a time-dependent background spectrum measured in an offsource region on the same detector chip as the source. The source continuum was modeled as an absorbed power-law, of variable amplitude and slope but time-invariant absorption column density. This continuum was subtracted, leaving positive and negative residuals that comprise noise plus any emission or absorption components on top of the continuum. More information about this method is given in Turner et al. (2006). The "signalto-noise" (S/N) map presented in Fig. 9 is the ratio of the fluctuation amplitude to the calculated noise.

The highest S/N observed is 4 and there are several points of significant excess emission above the continuum. However, we note here that the energy-time map highlights fluctuations above the continuum which may be real spectral features, but does not assess the significance of variations in those features, assuming they are real. The apparent fluctuations in the energy-time map could be consistent with the complex emission features detected in the time-averaged spectrum being constant, e.g. if the average significance in one independent 10 ksec segment is 2 sigma above the continuum, then we might expect statistical variations from 0-4 sigma in a 100 ksec observation.

To quantify whether the line features are variable, we split the data set into three time intervals, 0-40 ksec, 40-70 ksec and 70-100 ksec, which are chosen to maximize the differences in the Fe K line region which are highlighted by the energytime map. Note that in the case that the lines are constant, this approach essentially maximizes the effect of noise fluctuations that are picked up by the energy-time map, so that the statistical tests we carry out are in some sense the most optimistic for detecting line variability. The line profiles for each segment (data/model ratio for a power-law continuum over the 2-10 keV energy range) are shown in Fig. 10. We then fit these three subspectra with (individual or double) Gaussian line profiles; the results are reported in Table 2. Though the line profiles in Fig. 10 appear to be different, adding a double Gaussian line profiles at 6.5 keV and 7.0 keV (as found in the time-average spectrum) are significant (F-test significance=99.5-99.97%) in all segments. We also carried out Monte Carlo simulations to verify the significance of lines fitted to the individual sub-spectra; we detail these results in the Appendix (see A.3), but note here that they confirm that the lines are significantly detected in the individual sub-spectra.

Having verified the statistical significance of the lines in each of the 3 segments, we check for significant variability of the line flux. First, we fit simultaneously all segments with a power-law model plus two Gaussian lines with line energies fixed to 6.5 keV and 7.0 keV and line widths and line normalisations are free to vary within each segment. The contour plots at 90% confidence level for two interesting parameters (∆χ 2 =4.61) are reported for each line in Fig. 11, and show that assuming two Gaussian lines there are no statistically significant line variations. We also sim-Table 2. Best-fitting spectral parameters of PN time-resolved spectra in the 2-10 keV energy range with an absorbed (Galactic, N H =2.74×10 20 cm -2 ) power-law (PL) component plus a Gaussian profile (zgauss): . The line fluxes are expressed in 10 -6 ph cm -2 s -1 . F-test probabilities are calculated with the corresponding power-law model as reference. (*) The F-test probability for freeing the line width in the fit is 92%. (a) Monte-Carlo simulations have been performed (see text). Fixed parameters are indicated by (f). Fig. 11. Width and normalisation line contour plots at 90% confidence level (∆χ 2 =4.61) for segment 1 (dotted-line), segment 2 (continuous line), and segment 3 (dashed-line). Left panel: line at 6.5 keV. For segment 2, only an upper limit is found, corresponding to the continuous line. Right panel: line at 7.0 keV. For segments 1 and 3, only an upper limit is found. ulated 100 spectra using the fakeit command, and assuming an initial model consisting of a power-law plus a single Gaussian emission line at either 6.5 or 7.0 keV. We assumed that the line is intrinsically constant, with parameters described by values, and used an exposure time of 29.8 ks. We found that the observed values of line intensity and line width σ were usually consistent with the resulting spread in 100 simulated values. This suggested that we could not rule out at high significance the hypothesis that these parameters are intrinsically constant.

Discussion

We report the data analysis results of a long 100 ks XMM-Newton observation in November 2005. We have used both model-independent techniques (e.g., rms-spectra, flux-flux plots) and spectral analysis (time-averaged and time-resolved spectroscopy).

Timing analysis

Light curves reveal that this source was highly variable on short timescales (∼500 s), and was on average in a higher/softer state compared to the May 2001 observation. The Power Spectral Density (PSD) analysis indicates (at 95% confidence) a break at a frequency of 3.7( +1.0 -1.7 ) × 10 -4 Hz. This measured cut-off frequency is comparable to those measured in some other rapidlyvariable Seyferts, such as MCG-6-30-15 and NGC 4051. For a range in L bol from 10 44.5 -10 45 erg s -1 (i.e. assuming a minimum equal to the sum of FIR, optical, and X-rays luminosities, reported by Pietsch et al. 1998, and in this work), we infer from the mass-luminosity-time-scale relation of McHardy et al. (2006) that M BH ranges from 4 × 10 6 -10 7 M ⊙ and the source is accreting at or close to the Eddington rate (assuming even higher L bol implies that the AGN is a super-Eddington accretor).

The rms-spectra obtained using light curves binned to 5 ks and 500 s show that F var appears to decrease above about 2 keV for both bin timescales. Therefore the soft X-ray excess down to 0.3 keV appears more variable compared to the hard energy range. This is different from the typical rms-spectra from other Seyferts, which tend to show a broad peak in the variability amplitude between 0.8 and 2 keV, and a much lower variability amplitude below 0.8 keV (e.g., Fabian et al. 2002;Inoue & Matsumoto 2003;Ponti et al. 2006). From the flux-flux analysis, the existing data cannot distinguish between spectral pivoting of the continuum and a two-component origin for the spectral softening, primarily because the data do not span a broad enough flux range to see the characteristic change in gradient expected from the power-law. However, in the case of the two-component model being correct for this source, the fractional offsets measured in the flux-flux plots increase significantly toward higher energies (similarly to what is observed in MCG-6-30-15), as expected if we are witnessing a constant reflection component. Because there is no evidence for distant, cold material such as a torus ( §4), then if the Compton reflection component is real, the reflection is more likely associated with the accretion disc.

Spectral analysis: The FeKα line complex

Instead of a single highly redshifted line at about 5.4 keV as observed in the 4 ks May 2001 observation, we found in the November 2005 time-averaged observation two narrow Fe Kα emission lines at about 6.5 keV and 7.0 keV. The former rest-frame energy is compatible within the error bars with a line emitted by moderately ionized iron (∼Fe xix-Fe xxii). However this line is more likely a blue or red peak of a relativistic line, indeed a genuine Fe xix-Fe xxii (L-ions) emission line complex would suffer from resonant Auger destruction (Ross et al. 1996;[START_REF] Liedahl | Diagnostics of Astrophysical Plasmas: Theory, Experiment, and Observation[END_REF]). The latter line energy corresponds to highly ionized iron line (Fe xxvi) or to a highly blue-shifted iron line from lower ionization matter. The two features can be well represented by a double-peaked line profile from the disk with a rest-frame energy of 6.75-6.80 keV. In case of an annulus diskline profile an inner radius of about 200 R g is inferred. Alternatively the Fe xxvi line could be the signature of a strongly photoionized (log ξ >3), circumnuclear matter (seen in emission) with an iron overabundance of a factor 3, and/or very high column densities with N H,warm >> 10 23 cm -2 (Bianchi & Matt 2002). There is no significant evidence for line variability, however with the current data we also cannot reject the possibility that the lines do vary in response to continuum variations. Further XMM-Newton monitoring of this source can help to establish or reject line variability at different time-scales (from day(s) to months) and also allow us to determine if the line profile is indeed double-peaked and originates in the disk or if there exist independent or transient red-/blue-shifted emission lines. Shorter time-scales will only be reached by higher collecting area (and resolution at the iron Kα energy band) instruments such as Constellation-X and XEUS.

In the May 2001 observation of ESO 113-G010 we observed a line feature at 5.39 keV, which was detected at 99% confidence (based on Monte Carlo simulations) with EW= 233 ± 140 eV. We notice that a weak feature is seen at 5.4 keV near 20 ks in the energy-time map (Fig. 9), but this line has a detection probability of only 38% using Monte-Carlo simulations. The 90% upper limit on EW for a 5.39 keV line in the total observation is 22 eV (Section 4.2.1). Therefore, if the line observed in May 2001 is real, it is transient. We note however that, if we treat the total data set on ESO 113-G010 as a set of independent ∼ 4 ksec segments, the chance of seeing a 99% fluctuation in a single segment by chance becomes quite high. However, in support of the possibility the 5.39 keV feature was real, we also point out that the flux and spectral hardness in May 2001 corresponds to the extreme low and high values observed in 2005 (see Fig. 1), which might be more favourable to transient line formation (though we also note that such a line feature was not observed during the 2005 flux minimum).

In addition, as suggested by the 2001 observation, no significant constant narrow 6.4 keV Fe Kα line is observed with an upper-limit of < 32 eV, hence excluding any dominant emission from distant cold matter such as a torus. For a neutral solar abundance reflector subtending 2π sr -1 to the line of sight, illuminated by a Γ = 1.9 continuum and inclined at 45 • , the expected iron line equivalent width is ∼ 125 eV with respect to the primary continuum (George & Fabian 1991). Thus in ESO 113-G010 any cold, Compton-thick reprocessor outside of the line of sight is likely to subtend a much smaller solid angle, of the order < π/2 sr -1 .

The optical classification as a Seyfert type 1.8 of ESO 113-G010 (observed in November 1996; Pietsch et al. 1998) suggests the presence of a reddened broad line region (BLR). For a reddening of E(B-V) = 0.8-1.0, typical for a Seyfert 1.8 (Véron-Cetty & Véron 2006), the hydrogen column density would be equal to or greater than about 4-6×10 21 cm -2 (assuming a dust-to-gas ratio equal to or smaller than the Galactic value; Bohlin et al. 1978). However, during the present 2005 XMM-Newton observation (also suggested by the 2001 observation), instead of an intrinsic neutral absorption (in addition to the Galactic column density value of 2.74 × 10 20 cm -2 ) expected about below 1 keV, a soft excess is observed. A resolution of this apparent discrepancy, without the need to invoke a very high dust-to-gas ratio is either that dust that obscures the BLR clouds could lie out of the direct line of sight toward the X-ray emitter; or dust responsible for the BLR reddening could be mixed in with the warm absorber-emitter medium detected in the RGS spectrum. Alternatively, the apparent disagreement between the (non-simultaneous) optical and X-ray classifications could be explained by a time variability of the optical and/or X-ray absorbing properties (e.g., H1320+551: Barcons et al. 2003).

Fig. 2 .

 2 Fig. 2. 0.3-10 keV PSD of ESO 113-G010, including the ratio of the data to the best-fitting power-law plus constant noise component model (model shown as solid line).

Fig. 3 .

 3 Fig. 3. 0.3-10 keV PSD of ESO 113-G010, including the ratio of the data to an exponentially cut-off power-law (with lowfrequency slope α = -1), plus constant Poisson noise component (model shown as solid line).

Fig. 4 .

 4 Fig. 4. The rms-spectra of ESO 113-G010, for light curve bin times of 5 ks (top) and 500 s (bottom). See text for details.

Fig. 5 .

 5 Fig. 5. 3-10 keV versus 0.3-1 keV binned flux-flux plot of ESO 113-G010, compared with a best-fitting linear plus constant model (dashed line, representing two-component spectral variation) and a simple power-law (dotted line, representing high energy spectral pivoting). The two models are indistinguishable over the relatively narrow flux range covered by the data.

Fig. 6 .

 6 Fig.6. Fractional offset spectrum, showing fractional constant offsets (y-axis intercept divided by mean flux) for a linear plus constant model fitted to various energy bands versus the 0.3-0.5 keV reference band. In the case that the two-component model for spectral variability is correct, the offset spectrum shows lower limits on the strength of the constant component spectrum relative to the mean spectrum.

Fig. 7 .

 7 Fig.7. Spectra of ESO113-G010 (in the observer frame) observed in November 2005. A power-law has been fitted to the 2-5 keV data and extrapolated to lower and higher energies. A significant soft X-ray positive residual is seen below 1 keV, as well as the presence of a Fe Kα complex near 6-7 keV. For presentation, the data have been re-binned into groups of 15 bins for PN (10 for MOS), after group of a minimum of 100 and 50 counts per bin were used for the fit for the PN and MOS respectively. Left panel: 0.2-12 keV PN spectrum. Right panel: co-added MOS 0.3-10 keV spectrum.

Fig. 9 .

 9 Fig. 9. Energy-time map. The energy values are in the quasar frame. The contours represent 1.5-4.0 sigma above the continuum (in 0.5 increments), with darker bins representing an excess above the continuum. The segments 1-3 (delimited by red vertical lines) correspond to the time selections used for the timeresolved spectroscopic analysis. See text for explanation.

Fig. 10 .

 10 Fig. 10. The PN data/model sub-spectra of ESO 113-G010 (quasar frame) for the three time intervals. The data have been fitted assuming a power-law underlying continuum. Upper panel: 0-40 ks sub-spectrum (segment 1). Middle panel: 40-70 ks sub-spectrum (segment 2). Lower panel: 70-100 ks subspectrum (segment 3).

For the spectral analysis ( §4) the extraction region radius was taken to be the same for MOS1 and MOS2 (i.e.

′′ ), in order to co-add both spectra.

Acknowledgements. The XMM-Newton project is an ESA Science Mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). We thank the anonymous referee for fruitful comments and suggestions.

Appendix A: Monte-Carlo simulations

A.1. Method

We carried out rigorous tests of the significance of the lines in the time-averaged spectrum and in each segment using Monte Carlo simulations; see Porquet et al. (2004b) and Markowitz et al. (2006) for full details. For the null hypothesis, we assumed that the spectrum is simply an absorbed power-law continuum, with the same parameters as the absorbed power-law model fitted to the real data in each segment. We used the XSPEC FAKEIT command to create 1000 fake EPIC-pn spectra corresponding to this model, with photon statistics appropriate for each exposure (38.6, 24.0, and 26.9 ks for segments 1, 2 and 3, respectively), and grouped each spectrum to a maximum of 20 counts per bin. For each faked spectrum, we re-fit the null hypothesis model, yielding a "modified" null hypothesis model. We ran FAKEIT a second time using this re-fit model, again with the appropriate exposure time; this process accounts for the uncertainty in the null hypothesis model itself. For each new spectrum, we re-fit the null hypothesis model, recorded an initial χ 2 value, and then added a Gaussian component to the fit, as described below. The line centroid energy was constrained to the 4.5-8.0 keV range, and the line normalization was allowed to be positive or negative. We stepped the Gaussian centroid energy over this range in increments of 0.1 keV, fitting separately each time to ensure the lowest χ 2 value was found. For each spectrum, we compared the lowest χ 2 value with that from the null hypothesis fit, to obtain 1000 simulated values of the ∆χ 2 , which we used to construct a cumulative frequency distribution of the ∆χ 2 expected for a blind line search in the 4.5-8.0 keV range, assuming the null hypothesis of a simple power-law with no line is correct.

A.2. Results for the time-averaged spectrum

The observed 6.5 and 7.0 keV lines had ∆χ 2 values of 13.2 and 19.0, significant at 97.4% and 99.7% confidence, respectively, according to the simulations. However, since we have detected two lines that are close in energy, we can quantify their detection significance as a pair. The likelihood that both lines are due to photon noise is (1-0.974)(1-0.997) = 8×10 -5 , i.e., the pair is significant at just over 99.99% confidence. As an alternative quantification of the detection significance of the pair, we ran simulations which tested for the presence of a pair of narrow Gaussians with a fixed energy separation of 0.50 keV and with both line intensities left free. The observed ∆χ 2 value of 33.2 is significant at >99.9% confidence. However, as there was no a priori expectation for the value of the energy separation, using such a line profile tests only a small subset of the full range of line pair profiles possible, and so this probability should be treated as an upper limit. Assuming the annulus profile is correct, the observed ∆χ 2 value of 28.8 corresponds to a detection probability of >99.9%.

A.3. Results for the time-resolved sub-spectra

Segment 1 (0-40ks): The observed 6.51 keV line had a ∆χ 2 of 9.4, significant at 85.8% according to the Monte Carlo simulations. The observed 7.03 keV line has a ∆χ 2 of 4.2, significant at only 18.6% confidence, i.e., we cannot significantly rule out the hypothesis that the 7.03 keV feature is due to photon noise. However, the likelihood that both lines are due to photon noise, calculated in the same manner as for the time-averaged spectrum, is 0.116, i.e., the pair of lines is significant at 89.4%. Performing Monte Carlo simulations assuming two lines with a fixed energy separation yields (an upper limit of) 92.0% confidence. Segment 2 (40-70 ks): The simulations suggest the broad Gaussian 7.0 keV line, with an observed ∆χ 2 of 17.9, is significant at 99.7% confidence. Segment 3 (70-100 ks): The simulations suggest the broad Gaussian 6.7 keV line, with an observed ∆χ 2 of 13.9, for only 3 additional parameters is significant at 99.9% confidence.

List of Objects

'NGC 3516' on page 1


The ancestor of extant land plants diverged from the lineage leading to the Zygnematophyceae algae and gradually adapted to the terrestrial environment. Fossils of 450-475 million years old spore dyads and tetrads, which are structures associated only with land plant, are the first evidence of the changes in this branch of the green lineage [START_REF] Wellman | Fragments of the earliest plants[END_REF]. Wellpreserved fossil specimens tell us that less than 50 million years later the plant lineage had already diversified, with plants displaying subaerial and underground axes, and vasculature [START_REF] Harrison | The origin and early evolution of vascular plant shoots and leaves[END_REF]. Between these two milestones: a gap of 50 million years devoid of known macrofossils, leaving a gap even wider in our understanding of plant evolution. The morphology, physiology and overall biology of the first land plants remain an absolute mystery.

Although the hunt for fossils of the first land plants is ongoing, an alternative approach would be to infer, to predict, their morphology and physiology based on extant land plant species.

For such an inference, an accurate understanding of the plant phylogeny and the precise description of the traits of interest across this phylogeny are needed [START_REF] Delaux | Reconstructing trait evolution in plant evo-devo studies[END_REF].

In this issue of Current Biology Harris et al. (XXX) combined a consolidated plant phylogeny with the phylogenetic analysis of genes linked to stomata development to determine whether stomata were present in the most recent common ancestor (MRCA) of land plants.

The monophyly of Bryophytes supported once again

Resolving the deep nodes of the land-plant phylogeny, the ones close to their origin more than 450 million years ago, has been a major goal in evolutionary biology for decades. Over the last five years, the availability of newly sequenced genomes and transcriptomes combined with a wide array of molecular phylogenetic methods allowed two hypotheses to stand out. Harris et al. conducted yet another phylogenetic analysis on a carefully curated gene set. Although the debate is likely not over yet, this work, together with a number of other very recent studies [START_REF] De Sousa | Nuclear protein phylogenies support the monophyly of the three bryophyte groups (Bryophyta Schimp.)[END_REF][START_REF] Puttick | The Interrelationships of Land Plants and the Nature of the Ancestral Embryophyte[END_REF][START_REF] Li | Anthoceros genomes illuminate the origin of land plants and the unique biology of hornworts[END_REF], brings even more support to one of the two hypotheses: the monophyly of the Bryophytes that encompasses liverworts, mosses and hornworts. In other words, the bryophyte and the tracheophytes (vascular plants) have independently evolved for over 400 million years.

One essential consequence of the Bryophyte monophyly, highlighted elsewhere [START_REF] Rensing | Plant Evolution: Phylogenetic Relationships between the Earliest Land Plants[END_REF], is that any trait present only in the tracheophytes or the Bryophytes might be either ancestral (i.e. was present in the MRCA of the land plants) and lost in one of the two lineages, or derived in the lineage where it occurs. Ambiguously, Bryophytes are often described as lower plants or even primitive by contrast with higher plants referring to angiosperms. This is not appropriate as, obviously, vascular plants, that encompass angiosperms, and Bryophytes have evolved for the exact same time since they diverged from their MRCA. Therefore, both lineages should display a mix of derived and ancestral traits. In their study, Harris et al. perfectly illustrate this concept, referred to as heterobathmy, by studying the evolution of stomata.

Tracking gene losses consolidates the model of stomata evolution

Stomata are present in the sporophytes of most vascular plants in many mosses and hornworts. Most importantly, stomata are absent from all known liverworts.

With the ever increasing body of proof in favour of vascular plants and Bryophytes as sister monophyletic clades, the parsimonious scenario for the evolution of stomata is a single gain in the extant embryophyte ancestor and multiple losses in the Bryophytes. This single gain hypothesis is supported first by the homology of molecular mechanisms found between Physcomitrella patens and Arabidopsis thaliana by the analysis of mutants of orthologous genes [START_REF] Caine | An ancestral stomatal patterning module revealed in the nonvascular land plant Physcomitrella patens[END_REF]. In addition, the phylogenetic analyses of the entire genetic toolkit of stomata development and function such as in Harris et al. and in recent other studies [START_REF] Li | Anthoceros genomes illuminate the origin of land plants and the unique biology of hornworts[END_REF][START_REF] Chater | Origin and function of stomata in the moss Physcomitrella patens[END_REF] further consolidate the conservation between angiosperms, mosses and hornworts.

Apart from gene conservation, gene loss is another powerful way of inferring the evolutionary history of a trait. The principle in which a gene linked to a specific trait is lost following the loss of the trait is called co-elimination. In plants, this principle had been for instance used to validate the multiple losses of the nitrogen fixing symbiosis in the Fabales, Fagales, Cucurbitales, and Rosales over a previous theory of predisposition and multiple gains [START_REF] Van Velzen | Comparative genomics of the nonlegume Parasponia reveals insights into evolution of nitrogen-fixing rhizobium symbioses[END_REF][START_REF] Griesmann | Phylogenomics reveals multiple losses of nitrogen-fixing root nodule symbiosis[END_REF]. In the study by Harris et al., the sampling includes several species that do not display stomata including liverworts, some mosses and hornworts the lycophyte Isoëtes tegetiformans and the marine angiosperm Zostera marina. Those last two plants are vascular plants that adapted to an aquatic environment and for which the absence of stomata is unambiguously the result of secondary losses. Both Z. marina and I. tegetiformans have lost two genes from the stomata toolkit: TMM and EPF1/2, as a likely result of co-elimination. In the stomata-less bryophyte species, similar losses were observed indicating that, as in vascular plants, the selection pressure to maintain TMM and EPF1/2 was linked to the ability to form stomata. Therefore, the most parsimonious scenario of a single gain of stomata in Embryophytes followed by multiple independent losses is fully supported. Having determined the evolution of the trait, Harris et al. went on to infer the ancestral stomata toolkit based on gene gains, losses and duplications and the mechanistic knowledges acquired in angiosperms and mosses. One outstanding result of the proposed inference is that the genetic toolkit found in stomata-bearing Bryophytes is very different from the ancestral one.

Reversely, the one known in angiosperms overlaps significantly with the predicted ancestral stomatal genetic toolkit. This indicates that vascular plants, including angiosperms, retain part of the ancestral toolkit that has derivedand even be lostin extant Bryophytes.

The stomata of the first land plants

All the recent data, and the comprehensive study by Harris et al., infer that the MRCA of land plants had stomata on the sporophytic generation which were controlled by a more complex gene network than the one observed in extant Bryophytes, allowing to conclude that "… bryophyte stomata underwent reductive evolution". In other words this exemplifies the concept of heterobathmy and shows that, yes, Arabidopsis can be genetically more ancestral than Bryophytes.

In extant Bryophytes, stomata are restricted to the sporophyte where they contribute to its dehiscence and have been maintained only in part of the moss and hornwort lineages [START_REF] Chater | Origin and function of stomata in the moss Physcomitrella patens[END_REF][START_REF] Merced | Structure, function and evolution of stomata from a bryological perspective[END_REF][START_REF] Duckett | Exploding a myth: The capsule dehiscence mechanism and the function of pseudostomata in Sphagnum[END_REF]. In extant vascular plants, besides the losses in aquatic species, stomata are extremely conserved and involved in respiration and evapotranspiration, a trait linked to a functional vasculature in the sporophyte. Given that the function of the stomata, and the selection to maintain them, differs so strongly between the two sister-lineages of extant land plants, inferring their function in the MRCA of land plants becomes extremely difficult.

One can hypothesize that such a complex structure is more likely to have arisen from selective pressure for increased air intake and had been exapted for sporophyte maturation in Bryophytes than the other way around. In this scenario, the MRCA of land plants had stomata that would improve its CO 2 intake in the sporophyte. While what would later become vascular plants became gradually larger with a proper vasculature and megaphylls, the Bryophytes took the reverse route by favouring smallness which in turn released the need for elaborate air exchange and transpiration. So, as for the genetic pathway regulating their development, the stomata function in bryophyte may represent a derived state. Fossils from 400 million years ago are remarkably diverse with some displaying stomata in both the sporophyte and the gametophyte, a state not found in extant land plants. With the bryophyte monophyletic, the affinity of these fossils with the tracheophyte or protracheophyte lineages remains even more elusive than earlier. Being very speculative, we build on the data presented in Harris et al. and propose that these fossils may reflect the diversity of the land plants before the split between bryophytes and tracheophytes. This would result in a most parsimonious scenario for the evolution of stomata involving a single gain in relatively isomorphic gametophytes and sporophytes, followed by a single loss of stomata in the gametophytic generation before the divergence of the Bryophyte and tracheophyte lineages and the exaptation toward a non-essential role in sporophyte dehiscence in bryophytes (Figure 1).

Plant

  Evolution: When Arabidopsis is more ancestral than Marchantia Mélanie K. Rich MK & Pierre-Marc Delaux Laboratoire de Recherche en Sciences Végétales (LRSV), Université de Toulouse, CNRS, UPS, Castanet Tolosan, France Summary The quest for determining how looked like the plants that first lived on lands 450 million years ago is among the most exciting challenges in evolutionary biology. Recent works indicate that they displayed angiosperm-like stomata.Text: