© ESO, 2015

1. Introduction

Supermassive black holes reside at the centre of most galaxies. When accreting matter, these black holes emit a huge amount of energy, becoming active galactic nuclei (AGNs). The energy and radiation produced by matter accretion onto these black holes play an important role in determining their masses and spins and the properties of the host galaxy bulges. AGNs typically have a relatively short duty cycle and are in a low-luminosity state most of the time (Ho 2008). AGNs may become bright when major inflows of gas feed the compact object (Hopkins et al. 2006). Apart from intense episodes of accretion, the black hole is powered at a floor minimum level mostly by hot gas from colliding stellar winds in the central region. A star in the nuclear cluster can be scattered close to the black hole via dynamical encounters and be placed into an orbit too close to the black hole. When the black hole tidal force overcomes the stellar self-gravity (at the so-called tidal radius), the star starts to be ripped apart. About half of the stellar mass remains bound to the black hole, forming an accretion disc and powering a luminous, long-lasting (months to years) accretion flare (e.g. Rees 1988; Phinney 1989). Despite their elusiveness, a number of events have been reported in the literature, resulting in the complete disruption of the star (e.g. Renzini et al. 1995; Bloom et al. 2011; Gezari et al. 2012).

Tidal disruption events (TDEs) provide a unique probe to reveal the presence of otherwise quiescent black holes, allowing us to also study the accretion process in a different regime than that of AGNs (Komossa 2012; Gezari 2012). TDEs were first predicted and later on observed in the soft X-ray to ultraviolet bands where the peak of the accretion disc emission lies (Strubbe & Quataert 2009). Two relativistic TDEs were discovered based on high-energy triggers, thanks to the launching of relativistic jets along the line of sight (Bloom et al. 2011; Burrows et al. 2011; Cenko et al. 2012). Only a small number of TDEs have been observed in X-ray, UV, and optical surveys to date, however, mainly because of their low rate of occurrence and sparse observations (Donley et al. 2002).

TDEs are ultimate events. A star entering the tidal radius r_t ~ R_∗ (M_BH/M_∗)^1/3 (where M_BH is the black hole mass and M_∗ and R_∗ are the mass and radius of the star) is completely disrupted by tidal forces. The fate of a star orbiting its central black hole is then defined by its pericenter radius, r_p: if β = r_t/r_p ≳ 1, we have a TDE. What happens instead if β ≲ 1? If β is still close to unity, the star is thought to still feel the black hole tidal force but will survive the encounter. Hydrodynamical simulations show that if the passage is close (0.5 ≲ β ≲ 1), the star may lose matter that can then accrete onto the black hole (Guillochon & Ramirez Ruiz 2013). These events are probably more frequent than classical TDEs (roughly by a factor of ~10) and are repeated during the stellar orbital period (if the star is on a bound orbit), increasing their observability. On the other side, these events are less energetic because the involved mass captured by the black hole is lower.

Here, after a brief introduction on total and partial TDEs (Sect. 2), we consider the case of a nearby low-luminosity AGN, IC 3599, that showed strong X-ray activity in the past, possibly resulting from a TDE (Brandt et al. 1995; Grupe et al. 1995, see Sect. 3). In Sect. 4 we describe new Swift/XRT data showing a second strong increase in flux, together with the analysis of all X-ray data on the source. The X-ray spectral analysis is presented in Sect. 5. Based on spectral results, we converted count rates into a flux and then into luminosity. The fit of the overall 24 yr long X-ray light curve is described in Sect. 6. Optical data are discussed in Sect. 7. In Sect. 8 we apply partial TDE modelling to the light curve of IC 3599 in trying to constrain the involved star and its orbit. Our conclusions are presented in Sect. 9.

2. Tidal stripping essentials

A star orbiting to a massive black hole too closely will be torn apart by the tidal force of the compact object. This occurs at the tidal radius beyond which stellar self-gravity is not able to counteract the black hole tidal field and keep the star together. The tidal radius is defined as (1)where R_∗ and M_∗ are the radius and mass of the star, M_BH is the mass of the black hole, and the mass of the central black hole in IC 3599 (Grupe et al. 2015, see also below). The numerical value refers to a 1 M_⊙ and 1 R_⊙ star, and units are gravitational radii (r_g).

As the star is ripped apart by the tidal forces of the black hole, the debris is thrown into high-eccentricity orbits with a wide range of periods and with an energy range of (Lacy et al. 1982): (2)The distribution of mass as a function of energy is nearly constant (Rees 1988), as also shown by numerical simulations (Evans & Kochanek 1989; Lodato et al. 2009).

For a parabolic orbit, nearly half of the debris is unbound, leaving the system at high velocity. The other half will return to the black hole at different times depending on the initial eccentricity (i.e. energy). The return time of the first debris at pericenter is The return of material at pericenter continues at a rate driven by Kepler’s third law as (3)where the maximum mass inflow rate is (4)The peak rate is a function of the stellar structure (Lodato et al. 2009). After disruption, the stellar debrisis launched into very eccentric orbits and gradually returns to pericenter, where they circularise and form an accretion disc at r_circ ~ 2 r_p. The subsequent fall-back of this material onto the black hole is governed by viscous times. If the viscous time is shorter than t_min, then the fall-back of matter onto the central object is almost instantaneous and t_fb ~ t_min. Based on an α-viscosity disc prescription, we can evaluate the viscous timescale t_ν(5)where t_Kep is the Keplerian time at a given radius, and h is the disc half-height divided by the radius. Because the initial luminosity is close to (or even somewhat higher than) the Eddington limit, the disc is expected to be thick (Ulmer 1999), with h ~ 1. Thus, we can estimate the ratio of the viscous time to the t_min time as (6)and a thick disc will form and drain as the material circularises down to the last stable orbit.

The infalling star is usually assumed to be in a parabolic orbit and to undergo complete tidal disruption. Encounters may still occur with the pericenter slightly outside the tidal radius (β ≲ 1), resulting in some spill-over of stellar matter, and orbits can be bound and highly eccentric, resulting in periodic outbursts. This problem received less attention than the classical TDEs, but still a number of works exist (Guillochon & Ramirez Ruiz 2013; Macleod et al. 2012, 2013; Hayasaki et al. 2013; Ivanov & Novikov 2001). If the orbit is bounded the physics, is similar to the parabolic case in terms of the effects of black hole tides on the approaching star, especially if the eccentricity is high. In the long term, the stellar orbit may change, but the typical orbital binding energy is higher than the stellar binding energy, and the orbital angular momentum is larger than the rotational angular momentum of the star, indicating that the transfer of mass will not substantially alter the orbit (Macleod et al. 2013).

In particular, Guillochon & Ramirez Ruiz (2013) investigated the fate of a star undergoing partial stellar stripping. They considered the case of stellar polytropes of γ = 4/3 and γ = 5/3 and investigated the case for β within the interval 0.6–4.0 and 0.5–2.5, respectively. They found that total disruption occurs for a critical impact parameter different from unity, which is β_c = 1.85 for γ = 4/3 and β_c = 0.90 for γ = 5/3, respectively. They also provided fitting formulae for estimating the mass accretion peak rate, the peak time, and the total mass transferred for both cases, which are strong functions of β, especially for β ≲ β_c (see also the related errata corrige¹). We used these formulae in deriving stellar parameters.

3. IC 3599

One of the first putative TDEs occurred in the close active galaxy IC 3599 (92 Mpc, at a redshift of z = 0.021; Brandt et al. 1995; Grupe et al. 1995). It was discovered in the X-ray band as a bright, soft source during the ROSAT all sky survey in December 1990. Additional ROSAT observations (June 1992 to June 1993) found IC 3599 in a dimmer (~100) and somewhat harder spectral state (see Fig. 1). A Chandra observation in March 2002 found the source at a similar level (Vaughan et al. 2004). The very large flux decrease and the softness of the X-ray spectrum led several authors to suggest that this outburst was a TDE (Brandt et al. 1995; Grupe et al. 1995; Komossa & Bade 1999; Vaughan et al. 2004), even if IC 3599 is an active galaxy. Active galaxies are predicted to host more TDEs as a result of the perturbing presence of the disc, even if it is more difficult to reveal them because of the higher overall emission (Komossa 2012; Karas & Subr 2007). The optical spectrum is characterised by strong Hα, Hβ, and [OIII] lines, showing variations in response to the X-ray outburst and classifying IC 3599 as a type 1.5–1.9 Seyfert galaxy (Brandt et al. 1995; Grupe et al. 1995; Komossa & Bade 1999). The central black hole mass has been derived from a relation between the flux at 5100 Å and the width of the quiescent Hβ emission line to be (Sani et al. 2010). Grupe et al. (2015) argued that this estimate has been obtained using a broad line region scaling relation, which is only appropriate in case of an unabsorbed line of sight. Using instead a black hole mass to bulge K-band luminosity (Marconi & Hunt 2003) or the relation between the [OIII] velocity dispersion and the black hole mass (Nelson 2000), they found to be in the (2 − 12) × 10⁶M_⊙ range. Here we adopt a black hole mass of 3 × 10⁶M_⊙. For this mass, the source was well below the Eddington limit during the ROSAT all-sky survey observation.

Fig. 1 X-ray spectra fitted with the absorbed disc black body plus power-law model described in the text. Data were rebinned to reach a 5σ significance or by a factor of 5 for plotting purposes. In the upper panel black, red, and green (filled circles) data refer to ROSAT spectra (first, second, and third observations, respectively). Blue (open circles) data refer to the Chandra spectrum. Light blue, magenta, and orange (filled stars) to Swift spectra (first, second, and third observations, respectively). In the lower panel residuals in terms of χ are plotted with the same colour codings and symbols.

The (low-luminosity) AGN nature of IC 3599 may cast doubts on the tidal event interpretation because AGNs show flares that are related to disc activity or to the uncovering of a heavily absorbed X-ray source. We note here that the quiescent Chandra spectrum is well described by a soft power law with minimal intrinsic absorption, which considerably weakens the case for an absorbed source. At the same time, variations by a factor of ≳20–30 in AGNs are very rarely observed. We study AGN variability in detail, and based on a sample of highly variable AGNs that were intensely monitored by Swift, we can assess that a variability similar to the one observed in IC 3599 occurs by chance at ~4.5σ level (see Appendix A). In particular, the second flare is predicted to result by chance from known AGN variability at ≳4.0σ. Other models apart from AGN variability have been proposed but are not applicable to IC 3599. A binary system consisting of the central black hole and an orbiting star filling its Roche lobe at periastron can be ruled out because the stellar density required would be too low (Lasota et al. 2011). In addition, we confirm that the X-ray and optical transient position is consistent with the centre of IC 3599 (see Appendix B). Periodic optical outbursts were observed in the BL Lac object OJ 287 and explained by accretion instabilities onto a binary black hole (Tanaka 2013). However, this mechanism predicts no X-rays in quiescence because the innermost part of the disc is depleted. AGN instabilities in slim discs have also been put forward to explain this variability (Honma et al. 1991). However, the estimated mass accretion rate for IC 3599 falls well below the instability region, and the duration of the flares is much longer than theory predictions (Xue et al. 2011, see Appendix). All these findings lend support to the idea that IC 3599 underwent TDEs. Two different TDEs appear unlikely, however. The occurrence probability of a TDE in a galaxy is of the order of 10^-5 galaxy^-1 yr^-1 (e.g. Donley et al. 2002), which means that two different events in 25 yr results in a probability of ~6 × 10^-8. Even increasing the rate by a factor of 100 as a result of central binary black hole merger (e.g. Perets et al. 2007) would still result in a low probability ~6 × 10^-4. A binary disrupted by the central black hole would result in two different events, but the time delay between them is much too short (several days) to account for what we observe (Mendel & Levin 2015).

IC 3599 was also detected in the radio band by the Very Large Array in June 2012. This radio observation is difficult to accounted for and challenges the prediction of the standard blast wave model (Bower et al. 2013).

4. X-ray data preparation

X-ray data were collected over 24 years with different satellites. A log of all the observations is shown in Table 1.

4.1. ROSAT All Sky Survey data

RASS data are not straightforward to analyse. They were collected using the PSPC detector in scan mode. For this reason we decided to stick to the original analysis by Grupe et al. (1995).

4.2. ROSAT/PSPC

Data were extracted using the task XSELECT (v. 2.4c) from a circular region centred on the source with a 60 arcsec radius. The background spectra were extracted from an annular region centred on IC 3599, free of contaminating sources and support shadows, with inner and outer radii 125 and 250 arcsec, respectively. Data were retained in the 12–211 channel range, corresponding to an energy range of 0.1–2.4 keV. Because IC 3599 is always seen on-axis, we used the pre-canned response matrix relative to PSPC-B, gain-2 period (pspcb_gain2_256.rsp). Given the relatively low number of photons, we binned the data to a minimum of one photon per energy channel using the grppha tool and adopted Cash-statistics to fit the data. Different observations were grouped together to increase the number of photons, as indicated in Table 1.

Cash-statistics does not allow a reliable estimate of the goodness of fit. We estimated the goodness of fit by finding the best fit with Cash-statistics and then tested the data with χ² statistics, corrected for the low number of photons per bin through the Churazov weighting scheme (Churazov et al. 1996). This allowed us to determine the goodness of fit using χ² statistics.

4.3. ROSAT/HRI

Data were first combined into a single event file using XSELECT. Source photons were extracted from a circular region of 15 arsec radius. Background photons were extracted from a nearby circular region of 45 arcsec radius. No spectral data can be obtained from ROSAT/HRI data, and the count rate was converted into a ROSAT/PSPC count rate using PIMMS (see below).

4.4. Chandra/ACIS-S

Data were reprocessed using the CIAO 4.6 (and CALDB 4.5.9) repro task. Spectral data were extracted using the CIAO task specextract using a circular 3 arcsec region for the source and an annular region for the background with 12 and 17 inner and outer radii, respectively. The task provides the user with the corresponding response matrix (rmf) and ancillary response file (arf). Data were retained in the 0.3–10 keV energy range and were binned to a minimum of one photon per energy channel using the grppha tool.

4.5. XMM-Newton Slew survey

The field of IC 3599 was observed during one XMM-Newton slew on 2007 June 19 (MJD 54 270). The source was not detected with a (2σ) upper limit of 0.8 counts s^-1 in the EPIC-pn instrument (Grupe et al. 2015). This implies an upper limit on the 0.3–10 keV unabsorbed flux of 1.4 × 10^-12 erg cm^-2 s^-1 assuming the quiescent spectral model (see below). We did not consider this observation any further.

4.6. Swift/XRT

Data were reprocessed using xrtpipeline (v. 0.12.9) (CALDB 2014-07-30). Data were extracted from a circular radius of 71 arcsec and the background from a nearby circular region of 141 arcsec radius. Data were retained in the 0.3–10 keV energy range for spectroscopy and were binned to a minimum of one photon per energy channel using the grppha tool.

5. X-ray spectral analysis

Fig. 2 Long-term X-ray luminosity light curve of IC 3599. Error bars are at 1σ confidence level. ROSAT, Chandra and Swift count rates were converted into 0.01–10 keV unabsorbed luminosities by means of spectral fits assuming a source distance of 92 Mpc. The overall X-ray light curve has been fitted with a (t − t0)− 5/3 function repeating over a P0 time. Free parameters are the starting date, the power-law normalisation, the repetition time and a constant, setting the basic emission level of IC 3599. The quiescent level (AGN activity) has been evaluated as worsening the fit until a 10% null-hypothesis probability is attained. This is nicely consistent with the mean value and the observed variability in the power-law component used in the X-ray spectra fits.

We fitted together three ROSAT/PSPC, one Chandra/ACIS-S, and three Swift/XRT spectra (see Table 2) with the X-ray spectral fitting package XSPEC (v. 12.8.1g). All spectral fits were minimised using C-statistics, and the goodness of the fit was assessed using the Churazov-weighted χ² statistics. It is readily apparent that the overall X-ray spectra are soft, as testified by previous X-ray data analysis (Brandt et al. 1995; Grupe et al. 1995; Vaughan et al. 2004). An absorbed (using TBABS) power-law model with all the column densities tied together and the same photon index for all the observations returns a χ²-statistic value of 730.9 with 410 degrees of freedom. This power-law fit provides a null hypothesis probability of 10^-20. The power-law photon index is very soft Γ = 3.9 ± 0.2 (all errors were determined for ΔC = 1.0, i.e. 1σ errors for one parameter of interest). Leaving the power-law photon index free in each spectrum improves the fit. We obtained χ² = 493.3 with 404 degrees of freedom. The corresponding null hypothesis probability is still 0.2%. In addition, the photon index of the first ROSAT spectrum is extremely high with , and the column density N_H = (4.4 ± 0.1) × 10²⁰ cm^-2 is much higher than the Galactic column density of N_H = 1.2 × 10²⁰ cm^-2 (Kalberla et al. 2005). A free column density model with a fixed power-law photon index instead provides χ² = 619.4 with 404 degrees of freedom, corresponding to a null-hypothesis probability of 10^-11. A soft model such as a single black body provides similar results: χ² = 551.7 with 404 degrees of freedom, with free temperature and radius. The null-hypothesis probability is 10^-6. It is apparent from the fit that the black body model fails to account for a high energy tail (>2 keV) at low fluxes.

We also considered a two-component model. Given the many possibilities, we explored two different models, based on the physics of possible emission mechanisms. The first model comes from the proposal that this bursting activity comes from the partial disclosing of a heavily absorbed AGN. We thus modelled the X-ray spectra with a Galactic absorption plus a partial covering of a power-law component. The second model is physically motivated by a TDE. It consists of an accretion disc spectrum (accounting for the disc emission caused by the tidal contribution) and a power law accounting for the AGN activity (MacLeod et al. 2013).

Spectral fitting results are shown in Table 2. We explored a partial covering factor model (pcfabs within XSPEC) to determine whether the outburst observed in IC 3599 might be explained as the unveiling of a heavily absorbed source. This model has been envisaged to explain the wild erratic variations observed in WPVS 007 (Grupe et al. 2013). The fit was carried out with a model made by a fixed Galactic column density plus a partial covering model with the equal intrinsic column densities tied to all the spectra but variable fractions. This composite absorption component screens a power-law component with the same photon index but variable normalisation. We obtained a χ² = 627.4 for 304 degrees of freedom, resulting in a null-hypothesis probability of 10^-24 (see Table 2).

In the second model, the parameters were tied among different spectra. The absorbing column density was considered equal for all the spectra. The diskbb temperature was set to vary, while its normalisation was tied between the spectra (because it is related to the disc inner radius and probably corresponds to the innermost stable orbit). The power-law photon index was tied between spectra, but its normalisation was free to vary. The overall fit provides a χ² = 315.9 for 303 degrees of freedom. The null-hypothesis probability is 29%. Residuals are well distributed all over the entire energy band (see Fig. 1 and Table 2). The derived column density is somewhat in excess of the Galactic value. From the normalisation of the disc black body model we can estimate the inner disc radius (modulo the disc inclination, i, and a colour factor uncertain by a factor of f_c ~ 2). Considering a maximally rotating black hole, we can estimate it mass as ~8 × 10⁵ × (cosi/ 2)^1/2 × (f_c/ 2) M_⊙. This estimate, even if approximate, is barely consistent with the value derived from optical studies, indicating that the central black hole in IC 3599 is not particularly massive.

6. Light-curve fitting

We fitted together the luminosity light curve derived from all Swift X-ray data and archival ROSAT and Chandra data. To do this, we computed the conversion factors from count rates to unabsorbed fluxes in the 0.01–10 keV energy band for each spectrum based on the absorbed disc plus power law spectral model. These fluxes were then converted into luminosities adopting a distance of 92 Mpc for IC 3599 and used in Fig. 2. Spectra that were obtained as the sum of different observations were then split into single observations, and the same conversion factor was applied to all of them. The flux of the last ROSAT/PSPC observation was estimated based on the spectrum of the third ROSAT/PSPC observation (i.e. the spectrum closest in time). The flux of the RASS point was extrapolated to the 0.01–10 keV energy band based on the best-fit spectrum (Grupe et al. 1995). For the ROSAT/HRI count rate we adopted a different approach. We fitted the last ROSAT/PSPC spectrum with absorbed single-component models (power law, black body, and bremsstrahlung). Based on these models, we used PIMMS (v. 4.7b) to extrapolate the ROSAT/HRI rates to equivalent ROSAT/PSPC rates. An average rate was computed, and the standard deviation among the different models was added in quadrature. Then we converted these rates into fluxes as above.

The overall Swift light curve shows many similarities with the previous ROSAT light curve (Brandt et al. 1995; Grupe et al. 1995; Komossa & Bade 1999). Both bolometric light curves can be fitted with a t^{− 5/3} power law, pointing to (at least) two different TDEs. The similarities in the decay of the events suggest that we are observing a recurrent phenomenon and not the random occurrence of different accretion episodes. The non-disruptive passage of the same star provides a more comprehensive explanation.

Fig. 3 Upper panel: disc temperatures as derived from the spectral fits folded on a three-passages light curve. Error bars are at 1σ confidence level. Red points (ROSAT) refer to the first passage, the green point (Chandra) to the second passage, and blue points (Swift) to the third. Disc temperature evolution is fitted with a fixed T ∝ t− 5/12 power law (deriving from Ṁ ∝ t− 5/3 and T ∝ Ṁ1/4 typical of an accretion disc and therefore being a direct test for the disc cooling, Lodato & Rossi 2010). Lower panel: as above, but for two passages. Red dots (ROSAT and Chandra) refer to the first passage and green dots (Swift) to the second. The fit with a T ∝ t− 5/12 power law is shown as a dashed orange line.

We started by fitting the luminosity light curve with a different number of outbursts. Values of the χ² are reported in Table 3. It is readily apparent that the fit with three outbursts and a power-law shape is superior to all the others. The improvement of a free power-law index with respect to the fixed value −5/3 was evaluated by means of an F-test. The probability of a random occurrence is 12%.

In addition to the luminosity light curve, we folded the temperatures derived from the spectral analysis along with the suggested orbital period. The disc temperature is expected to evolve as T ∝ t^{− 5/12}, where Ṁ ∝ t^{− 5/3} and T ∝ Ṁ^1/4 (from accretion disc theory, Lodato & Rossi 2010). By fitting the temperature evolution for three outburst peaks, we obtain a reduced for six degrees of freedom and a null-hypothesis probability of 75% (upper panel in Fig. 3). If we fit instead the temperature evolution for two outbursts with the same model, we derive a reduced for six degrees of freedom and a null-hypothesis probability of 10^-38 (orange dashed line in Fig. 3). If we add a constant to the power-law model, mimicking the presence of an underlying quiescent accretion disc, we obtain for five degrees of freedom and a null-hypothesis probability of 4 × 10^-6 (continuous line in the lower panel of Fig. 3). The constant quiescent-disc temperature is 42 ± 6 eV. We note that such a constant temperature baseline might also be present in the first fit, but it is not required by the data. We note that the three peaks are only based on the Chandra data. If these data are strongly affected, timing and spectral variability unrelated to the TDE, a longer period (twice as long) is expected.

Based on light curve and disc temperature evolution, we found a strong indication for multiple (equal) tidal stripping events. In particular, the timing and spectral data point toward three (not two) outbursts taking place during 1990–2014, with a recurrence time of 9.5 yr (see Figs. 2 and 3). Assuming an accretion efficiency of η = 0.1, we can estimate from the outburst light curve the peak mass accretion rate to be Ṁ_peak ~ 0.01 M_⊙ yr^-1 and the accreted mass per episode to be ΔM ~ 2.5 × 10^-3M_⊙. These are clearly lower limits, missing the early stages of all the outbursts. In addition, with these new ephemerides, radio observations took place only 2.6 yr after an outburst episode and can be accounted for by emission internal to the jet (van Velzen et al. 2011).

7. Optical data

Fig. 4 Optical (flux) light curve of IC 3599 from the Catalina sky survey in arbitrary units (red dots), together with the Swift/XRT light curve (open squares) converted into flux using the spectral fit model, in arbitrary units. The Gaussian fit refers to the optical data. The vertical dashed line marks the time (MJD) of the optical peak, the dotted line the peak of the X-ray flux according to the multi-peak modelling. The continuous line marks the time at which the optical light curve rises by 10% from the quiescent level.

IC 3599 fell within the field of view of the Catalina sky survey and was monitored in the optical in the MJD 53 470–56 463 interval (Drake et al. 2009). The Catalina light curve shows a broad peak (see Fig. 4). The light curve can be fitedt with a Gaussian centred on MJD 55 151 ± 11 (1σ) and with a width of 200 ± 18 d. The peak of the optical emission occurs after the estimated peak of the X-ray emission by ~140 d. However, due to the large width of the optical flare, the optical emission starts ~280 d before the high-energy peak (assumed as the start of the optical flare when the optical flux rises by 10% over the constant value).

Fig. 5 Radius of the orbiting star and impact parameter β = rt/rp as a function of the stellar mass. The allowed regions were computed based on fitting formulae in Guillochon & Ramirez-Ruiz (2013 and related errata corrige). Upper panel: orbiting stars characterised by a polytropic index γ = 4/3. The red region in the lower panel shows the allowed β parameters. The light red region approximately excludes the low-mass end (≲1 M⊙) and the stars more evolved than the giant branch because they are not well described by a polytropic index γ = 4/3. Orange dots indicate stars in the main -sequence phase, yellow dots the Hertzsprung gap, black dots the giant branch, light blue dots the helium core-burning, and blue dots the asymptotic giant branch phase. Lower panel: polytropic star with index γ = 5/3. Colours, dots, and regions are as above. The light red region approximately excludes the high-mass end (≳1 M⊙), because these stars are not well described by a polytropic index γ = 5/3 (but includes evolved giant stars). No acceptable solutions remain for this case.

This enables us to study the outburst start in detail (and prepare a follow-up strategy for the next outburst). The rise time in a TDE is governed by the circularisation time (i.e. the time to form the accretion disc) and the viscous time (i.e. the time needed to transfer matter from the outer disc edge to the central compact source). Several papers have recently appeared on this subject (Bonnerot et al. 2015; Guillochon & Ramirez-Ruiz 2015; Piran et al. 2015). Irrespective of which is the dominant mechanism, we can safely assume that the viscous time of the disc, t_ν, is shorter than the observed rise time. Judging from the optical light curve, the time it takes to rise from 10% of the quiescent flux to the peak flux is ~420 d. Fitting the optical light curve with a symmetric exponential function, we derive an e-folding rise/decay time of 178 ± 14 d. Based on Guillochon & Ramirez-Ruiz (2013), the viscous time can be expressed as (7)

8. Constraining the orbit and star characteristics

Partial disruption of stars has been investigated by means of hydrodynamical simulations (Guillochon & Ramirez Ruiz 2013; Macleod et al. 2013, 2012). Depending on the stellar structure (modelled as a polytrope with index γ = 4/3 or 5/3), it has been shown that mass can be extracted from the star for an impact parameter β as low as 0.5–0.6, while complete destruction occurs for β> 0.90 for γ = 5/3 and for β> 1.85 for γ = 4/3, respectively (Guillochon & Ramirez Ruiz 2013). To constrain the encounter and star characteristics, we compared the mass accretion at peak rate Ṁ_peak and the total amount of matter 2 ΔM lost by the star during each passage with the corresponding quantities estimated through hydrodynamical simulations (Guillochon & Ramirez Ruiz 2013). These quantities strongly depend on the impact parameter β and, to a lesser extent, on the stellar mass and radius.

The peak mass accretion rate was estimated based on the RASS data, converting the peak luminosity into a peak mass accretion rate by assuming a 10% conversion efficiency. To constrain the uncertainties, the allowed parameter space was computed accepting values of the peak accretion rate in the interval 1–3 of observed value. The other parameter considered is the total amount of mass accreted during one flare episode. This was estimated by integrating the luminosity light curve over one orbital period. For this parameter we also searched solutions in the interval 1–3 of the observed total mass accreted. Finally, we required that the viscous time of the disc, t_ν, is shorter than 178 ± 42 d (3σ). A mass range of 0.2–100 M_⊙ and a (unrelated) radius range of 0.01 − 10⁴R_⊙ were blindly searched. We considered two different star models based on a polytropic index of γ = 4/3 (over the interval β = 0.5 − 4.0) and γ = 5/3 (β = 0.45 − 2.5). Results are shown in Fig. 5 for the two indexes separately.

For γ = 4/3 we find solutions in the range β = 0.58 − 0.72. The range of allowed values of β was obtained without requiring a priori a partial disruption event, searching for a solution in the β = 0.5 − 4.0 interval. This result therefore strongly (and independently) supports that the flaring events in IC 3599 are related to a partial tidal disruption (β<β_c = 1.85) of an orbiting star. The corresponding mass range is for M_∗> 4 M_⊙. We blindly explored the mass-radius plane, but stars do not fill this plane homogeneously. To test the consistency of our findings with stellar models, we investigated the mass-radius relation by means of the single stellar evolution (SSE) code (Hurley et al. 2000), assuming solar metallicity (we verified that a change in the metallicity does not change our results by much). With this additional constraint, we have that the range of allowed masses reduces to M_∗ = 15–45 M_⊙ for stars on the main sequence and to M_∗ = 4 − 15 M_⊙ for slightly evolved stars. The corresponding radii in the selected region are R_∗ ~ 5 − 9 R_⊙ for main-sequence stars and R_∗ ~ 4 − 7 R_⊙ for slightly evolved stars (see Fig. 5, upper panel). The allowed eccentricities range in e ~ 0.995 − 0.997 and the pericenter distance r_p ~ 84 − 93 r_g (where is the gravitational radius of the black hole in IC 3599, G the gravitational constant, and c the speed of light). The eccentricity is high but not unlikely, given the predicted distribution of eccentricities around the central black hole in our Galaxy (Gillessen et al. 2009). This eccentricity is higher than the critical eccentricity for bound orbits below which all stellar debris remains bound, and it feeds the black hole on a longer timescale (Hayasaki et al. 2013).

For a polytropic star with index γ = 5/3, we searched for solutions for β in the 0.45–2.5 range (see Fig. 5, lower panel). We did find solutions in a narrow range of β = 0.49 − 0.55. Again, the allowed range of β is well below the critical value for a complete tidal disruption (β_c = 0.9). The allowed mass range is M_∗ ≳ 1.5 M_⊙. However, polytropic stars with γ = 5/3 can satisfactorily describe Sun-like stars (or smaller), which means that we did not find acceptable solutions in this case.

9. Conclusions

We reported here on the first candidate of periodic partial tidal disruption events. These events were observed as large (>100) flux increases in the 24 yr long X-ray light curve of the close low-luminosity active nucleus IC 3599 (see Fig. 2). We discussed several possibilities to explain these strong flares and conclude that they cannot be ascribed to AGN variability: single flares are unlikely, at ≳4.0σ level, to come from known observed extremely variable AGNs, and the flaring instability region suggested by some accretion disc models (e.g. Honma et al. 1991) lies above the mean accretion rate observed in IC 3599 by a factor of ≳15. By modelling the light curve and the disc temperature evolution, as derived by spectral fits, we discovered that three flares are statistically preferable over two, resulting in a ~9.5 yr periodicity. This cannot be appreciated from the light curve where the second flare is largely missed due to sparse observations, but it is clear from the disc temperature evolution: if the standard disc model can apply to the observations of IC 3599, then the three-outbursts model is preferable over the two-outburst model at ~3.8σ. Based on the disc temperature evolution, the same is true at ~4.3σ confidence level. A tidal disruption model applied to IC 3599 provides a good match to the data and, independently, solutions for partial disruption events.

In addition, we note that the average quiescent luminosity of the AGN is just a fraction of the average flare luminosity (~20%), leaving open the possibility that the overall AGN activity of IC 3599 is entirely contributed by the orbiting star (MacLeod et al. 2013). Given the relatively high mass stripped at every passage, the phenomenon is short-lived, and based on our estimates, it can last ≲10⁴ yr. The next passage will probably occur in 2019, giving us the ability to plan a detailed monitoring campaign to explore the characteristics of the orbiting star. The very short orbital period enabled us to observe more than one passage. Other tidal disruption events presently known might be similar to this case, but just with a longer orbital period.

Acknowledgments

We acknowledge useful discussions with G. Tagliaferri, G. Ghisellini and R. Salvaterra. We thank R. Campana for useful discussions about power spectral densities and for his public python simulation software.

References

Appendix A: AGN variability

Given the low-luminosity AGN nature of IC 3599, we here investigate in detail whether the observed flares can come from this activity. We approach the problem either by studying standard AGN variability or by taking an unbiased sample of the most variable AGNs as observed by the Swift satellite.

Appendix A.1: Intrinsic AGN variability

We simulated the long-term light curve of an AGN (Vaughan et al. 2003). We assumed a standard broken power law for the power spectral density (PSD) of a long-term (years) monitoring light curve. We considered an index –1 at short frequencies and –2 at high frequencies, with a break at 10³ Hz scaled to the mass and luminosity of the black hole in IC 3599 (McHardy et al. 2006). We checked that even when these numbers were stongly changed, the final result did not depend on them. We ran a simulation that generated a number of light curves based on the PSD above, including no root-mean-square (rms) variability and no background. The long-term light curve changes in all cases are lower than a factor of a few with respect to the mean starting value (on timescales of days). We then turned on the rms variability, and this value is the main driver for the variability. Even assuming a 100% rms variability, we obtained a maximum increase by a factor of ~6–7 in the count rates (no background included) from the mean value on timescales longer than days. From this analysis we conclude that normal AGN variability is not able to produce the strong variability we observe.

Fig. A.1 Swift/XRT light curves of the four selected highly variable AGN.

Appendix A.2: AGN with strong intrinsic variations

We started from Swift/XRT light curves spanning a 10 yr basis. We first considered the AGN sample monitored by BAT, consisting of bright nearby AGNs with more than five Swift/XRT observations. In this sample there are just three AGNs showing a count rate change by more than a factor of 20. These are Mkn 335 (Δ = 42), NGC 4395 (Δ = 141) and CGCG 229–015 (Δ = 56, even if in this case the change is a dip and not a flare). One very efficient method to find highly variable AGN consists of comparing soft X-ray observations at different epochs. Grupe et al. (2001) showed that 4 out of 113 bright RASS objects that were also observed in ROSAT pointings showed dramatic variations. Grupe et al. (2012) added another four objects to this list by means of subsequent XMM-Newton and Chandra observations. All but two are narrow-line Seyfert 1 galaxies. The first is RX J1624.9+7554, a non-active galaxy for which a TDE has been invoked (Grupe et al. 1999), the second one is IC 3599. In this Grupe-Komossa sample of highly variable AGNs there are then eight objects, but only three of them show variations stronger than a factor of 30. These are 1H 0707–495 (Δ = 114, dip-like variation), Mkn 335 (as above), and WPVS 007 (Δ = 36). WPVS 007 is too faint for Swift, however, and the light curve is unusable because it has too few points (Grupe et al. 2013). We have then four objects with good enough sampled light curves that show flare-like features (Fig. A.1).

Fig. A.2 Hardness ratio as a function of the normalised count rates for the variable AGNs described above. 1H 0707–495 data are shown as light blue open stars, CGCG 229–015 as black open squares, NGC 4395 as blue open circles, Mkn 335 as green open triangles, and IC 3599 as filled red circles. Count rates are normalised to the minimum rate observed. The hardness ratio is computed as the ratio as (H − S)/(H + S),whereH is the number of counts in the 2–10 keV energy band and S the number of counts in the 0.5–2 keV band.

One conservative test is to randomly extract from seven points (the number of points comprising the first TDE flare) and six points (the number of points comprising the third TDE flare) separately from these curves and fit them with a tidal decay template (t − t₀)^5/3 + c to mimic one TDE at a time. We simulated 500 000 light curve realisations for each source by randomly sampling the observed light curves, extracting seven or six points, respectively. We then excluded simulated light curves with a rate variation lower than a factor of 30 (so as not to fit constants, note that even if the flux variation in the Swift/XRT data is ~30 the count rate variation is ~100 due to spectral variability, therefore our rate change provides a conservative estimate), and fit the remaining light curves with the above model. We adopted two different approaches. In the first one we counted how many simulations have a null-hypothesis probability higher than 5% and, based on this number, evaluated the probability of obtaining by chance a TDE-like light curve. In the second approach we took the best of all these light curves in terms of χ² and derived the probability of this χ² with respect to the number of degrees of freedom (five and four in our case, respectively). We did this because even the best selected χ² is, for some sources, not good, at variance with the fit of the IC 3599 flare. This χ² probability was then weighted for the number of trials (500 000). With these numbers we assessed the probability of randomly extracting a TDE-like event from these light curves with a count rate increase by a factor of >30 (the Swift rate variation is ~100). We obtained conservative limits in this way that are shown in Table A.1.

A different check can be made on the spectral properties of IC 3599. In Fig. A.2 we show the spectral evolution of the five sources in our sample. This is quantified in terms of hardness ratio (computed as HR = (H − S)/(H + S), where H and S are the counts in the 2–10 keV and 0.5–2 keV energy bands, respectively, so that HR can vary among −1 and 1) versus cont rates normalised to the minimum observed rate. In this plane IC 3599 stands alone, reaching complete softness (HR = −1) at maximum and showing a clear spectral evolution. 1H 0707–495 is similarly soft, but it does not show a marked spectral variability as the count rate changes. Mkn 355 and NGC 4395 show spectral marked changes but never appear as very soft sources.

We can safely conclude that the flares observed in IC 3599 probably do not come from known (observed) AGN variability.

Appendix A.3: Disc instability models

The standard geometrically thin, optically thick accretion disc model has thermally unstable regions if the accretion rate is higher than a critical value. A new branch of equilibrium solutions should exist for rate across the Eddington luminosity, named slim discs (Abramowicz et al. 1988). Some works have shown that a limit cycle might take place, resulting in a flare-like behaviour (Honma et al. 1991). Magneto-hydrodynamical simulations showed, on the contrary, that such discs could be thermally stable (Hirose et al. 2009). Observationally, Galactic transient X-ray binaries in outburst do not show any evidence for such an instability, except perhaps GRS 1915+105 (Xue et al. 2011).

The mean mass accretion rate of IC 3599 is ~0.005 in Eddington units (for an efficiency of 10%, where the peak accretion rate reaches ~0.1 the Eddington rate), far away from border of the unstable region, starting at a critical mass accretion rate >0.1 in Eddington units (Honma et al. 1991). Flare models do not make a clear prediction of the shape of the flare, but the make key prediction that the flare width is a very short fraction of the times between two different flares (Xue et al. 2011). The ratio of the FWHM of the flare with respect to the duration time between two consecutive flares must be ≲3.6%. We modelled the two main peaks of IC 5399 with a Gaussian, limiting the maximum luminosity to a few times the Eddington value, and found a FWHM ~ 500 d (leaving the maximum luminosity free to vary provides a better fit and an even larger FWHM ~ 1200 d). Taking the times among the two main flares, we obtain a ratio of ~7% (~17% for a free Gaussian shape). This is more than twice of what is expected based on models.

We also note that two possible very recent flare-like events in AGNs, XMM SL1 J061927.1–655311 (Saxton et al. 2014) and NGC 2617 (Shappee et al. 2014), involve weaker flux variations (≲20) and much more complicated X-ray light curves.

Appendix B: UVOT analysis

Fig. B.1 Swift/UVOT images of IC 3599, obtained with the uvm2 filter on 2010 Feb. 25 (left panel) and 2014 March 26 (central panel). The images are 5′ × 5′, north is up, east is left. The right panel shows the result of digital image subtraction between these two images. A clear residual is visible at a position consistent with the IC 3599 galaxy centre.

Appendix B.1: Image analysis

Together with XRT exposure, the Swift satellite took contemporaneous images of IC 3599 with the UVOT instrument. During all the six Swift observations (except for the fourth), UVOT observed IC 3599 with all its six filters (in the fourth only with the uvm2 filter). Sky images were considered and fluxes were obtained with the uvotsource task using the latest calibrations (Poole et al. 2008; Breeveld et al. 2011). For UV images a circular extraction region of 8 arcsec radius was considered and for optical image a 5 arcsec radius. A background region close to IC 3599 free of sources and of 20 arcsec radius was selected. Results are shown in Table B.1. A clear decrease in flux is apparent, especially at UV wavelengths with a decrease by a factor of ~2 (uvm2).

Appendix B.2: Astrometry

To precisely determine the position of the transient source observed in IC 3599, we carried out digital image subtraction between the first and the last UV frames (uvm2) obtained with the Swift/UVOT. In the first image, obtained on 2010 February 25, the transient source is at its maximum, while in the last epoch, obtained on 2014 March 26, it has faded below the host galaxy level (Table B.1). Both images were obtained with the uvm2 filter and have exposure times of 0.2 ks and 1.2 ks, respectively. Before subtraction, the two images were aligned. We accounted for the difference in the exposure times by multiplying the image obtained in February 2010 by a factor of six. As can be seen in the third panel of Fig. B.1, the result of the image subtraction is good, showing a clear residual at the position of IC 3599. The position of this residual is RA, Dec (J2000): 12:37:41.18, +26:42:27.0 (with an uncertainty of 0.3 arcsec, calibrated against the USNOB1.0 catalogue). At the distance of IC 3599 1 arcsec corresponds to 0.4 kpc. The residual is located at 0.8 ± 0.4 arcsec (1σ confidence level) with respect to the position of the IC 3599 galaxy centre and is fully consistent with the position of the candidate TDE radio counterpart (Bower et al. 2013). The consistency of the position of the UV transient with the IC 3599 galaxy centre (at the 2σ level) agrees with a TDE origin.

All Tables

All Figures

Fig. 1 X-ray spectra fitted with the absorbed disc black body plus power-law model described in the text. Data were rebinned to reach a 5σ significance or by a factor of 5 for plotting purposes. In the upper panel black, red, and green (filled circles) data refer to ROSAT spectra (first, second, and third observations, respectively). Blue (open circles) data refer to the Chandra spectrum. Light blue, magenta, and orange (filled stars) to Swift spectra (first, second, and third observations, respectively). In the lower panel residuals in terms of χ are plotted with the same colour codings and symbols.

In the text

Fig. 2 Long-term X-ray luminosity light curve of IC 3599. Error bars are at 1σ confidence level. ROSAT, Chandra and Swift count rates were converted into 0.01–10 keV unabsorbed luminosities by means of spectral fits assuming a source distance of 92 Mpc. The overall X-ray light curve has been fitted with a (t − t0)− 5/3 function repeating over a P0 time. Free parameters are the starting date, the power-law normalisation, the repetition time and a constant, setting the basic emission level of IC 3599. The quiescent level (AGN activity) has been evaluated as worsening the fit until a 10% null-hypothesis probability is attained. This is nicely consistent with the mean value and the observed variability in the power-law component used in the X-ray spectra fits.

In the text

Fig. 3 Upper panel: disc temperatures as derived from the spectral fits folded on a three-passages light curve. Error bars are at 1σ confidence level. Red points (ROSAT) refer to the first passage, the green point (Chandra) to the second passage, and blue points (Swift) to the third. Disc temperature evolution is fitted with a fixed T ∝ t− 5/12 power law (deriving from Ṁ ∝ t− 5/3 and T ∝ Ṁ1/4 typical of an accretion disc and therefore being a direct test for the disc cooling, Lodato & Rossi 2010). Lower panel: as above, but for two passages. Red dots (ROSAT and Chandra) refer to the first passage and green dots (Swift) to the second. The fit with a T ∝ t− 5/12 power law is shown as a dashed orange line.

In the text

Fig. 4 Optical (flux) light curve of IC 3599 from the Catalina sky survey in arbitrary units (red dots), together with the Swift/XRT light curve (open squares) converted into flux using the spectral fit model, in arbitrary units. The Gaussian fit refers to the optical data. The vertical dashed line marks the time (MJD) of the optical peak, the dotted line the peak of the X-ray flux according to the multi-peak modelling. The continuous line marks the time at which the optical light curve rises by 10% from the quiescent level.

In the text

Fig. 5 Radius of the orbiting star and impact parameter β = rt/rp as a function of the stellar mass. The allowed regions were computed based on fitting formulae in Guillochon & Ramirez-Ruiz (2013 and related errata corrige). Upper panel: orbiting stars characterised by a polytropic index γ = 4/3. The red region in the lower panel shows the allowed β parameters. The light red region approximately excludes the low-mass end (≲1 M⊙) and the stars more evolved than the giant branch because they are not well described by a polytropic index γ = 4/3. Orange dots indicate stars in the main -sequence phase, yellow dots the Hertzsprung gap, black dots the giant branch, light blue dots the helium core-burning, and blue dots the asymptotic giant branch phase. Lower panel: polytropic star with index γ = 5/3. Colours, dots, and regions are as above. The light red region approximately excludes the high-mass end (≳1 M⊙), because these stars are not well described by a polytropic index γ = 5/3 (but includes evolved giant stars). No acceptable solutions remain for this case.

In the text

	Fig. A.1 Swift/XRT light curves of the four selected highly variable AGN.
In the text

Fig. A.2 Hardness ratio as a function of the normalised count rates for the variable AGNs described above. 1H 0707–495 data are shown as light blue open stars, CGCG 229–015 as black open squares, NGC 4395 as blue open circles, Mkn 335 as green open triangles, and IC 3599 as filled red circles. Count rates are normalised to the minimum rate observed. The hardness ratio is computed as the ratio as (H − S)/(H + S),whereH is the number of counts in the 2–10 keV energy band and S the number of counts in the 0.5–2 keV band.

In the text

	Fig. B.1 Swift/UVOT images of IC 3599, obtained with the uvm2 filter on 2010 Feb. 25 (left panel) and 2014 March 26 (central panel). The images are 5′ × 5′, north is up, east is left. The right panel shows the result of digital image subtraction between these two images. A clear residual is visible at a position consistent with the IC 3599 galaxy centre.
In the text