The Large Amplitude X-ray Variability in NGC 7589: Possible Evidence for Accretion Mode Transition

We report the discovery of large amplitude X-ray variability in the low luminosity AGN (LLAGN) MGC 7589, and present possible observational evidence for accretion mode transition in this source. Long-term X-ray flux variations by a factor of more than 50 are found using X-ray data obtained by Swift/XRT and XMM-Newton over 17 years. Results of long-term monitoring data in the UV, optical and infrared bands over ~20 years are also presented. The Eddington ratio increased from $10^{-3}$ to $\sim0.13$, suggesting a transition of the accretion flow from an ADAF to a standard thin accretion disc. Further evidence supporting the thin disc in the high luminosity state is found by the detection of a significant soft X-ray component in the X-ray spectrum. The temperature of this component ($\sim19^{+15}_{-7}$eV, fitted with a blackbody model) is in agreement with the predicted temperature of the inner region for a thin disc around a black hole (BH) with mass of $\sim10^{7}\,M_{\mathrm{Sun}}$. These results may indicate that NGC 7589 had experienced accretion mode transition over a timescale of a few years, suggesting the idea that similar accretion processes are at work for massive black hole and black hole X-ray binaries.


INTRODUCTION
It is believed that Active Galactic Nuclei (AGNs) are powered by accreting material onto massive black holes (BHs) at the centre of galaxies. Observational evidence of an optical thick, geometrical thin standard accretion disc (hereafter thin disc, Shakura & Sunyaev 1973) has been found in the multi-band data of luminous AGNs (e.g. Seyferts and QSOs). For instance, the big blue bump observed in the spectral energy distribution (SEDs) of AGNs is explained as the thermal radiation from the thin disc (e.g. Elvis et al. 1986;Czerny & Elvis 1987;Shang et al. 2005). In addition, the broad Fe Kα line, which is generally thought to originate from the inner region of an accretion disc and broadened due to relativistic effects of the central BH (Fabian et al. 1989(Fabian et al. , 2000, has been detected in the X-ray spectra of some AGNs (e.g. Tanaka et al. 1995;de La Calle Pérez et al. 2010), supporting the thin disc model with an inner radius extended to a few gravitational radius (R g = GM/c 2 ) to the BH. The high energy tail of the thermal radiation from the thin disc, which can be modelled by a blackbody with temperature of several tens of eV (depends on the BH mass), has also been discovered in the X-ray spectra of some AGNs 1 (Yuan et al. 2010;Sun et al. 2013;Shu et al. 2017). On the other hand, an optical thin radiation Contact e-mail: liuzhu@nao.cas.cn 1 However, the spectra in the soft X-ray band of the vast majority of AGN show a component with a higher temperature (∼ 0.1 − 0.2 keV) than the prediction of the thin disc -the so-called soft X-ray excess whose origin is still under debate (e.g. Arnaud 1996; Gierliński & Done 2004) inefficient accretion flow (RIAF) or advection dominated accretion flow (ADAF, Narayan & Yi 1994;Abramowicz et al. 1995;Yuan & Narayan 2014) has successfully explained the observed property of low luminosity AGNs (e.g. LLAGNs and LINERs, Ho 2008) that normally has very low accretion rates.
Theoretical calculation suggests that transition between different accretion mode will take place once the accretion rate reaches certain critical values (Meyer et al. 2000a). Evidence for accretion mode transition has long been found in Galactic X-ray binaries (XRBs), e.g. Zhang et al. (1997). In the high/soft state, a thermal component with temperature of ∼ 1 keV seen in the soft X-ray band is widely accepted as the thermal radiation from the thin disc. While in the low/hard state, the hard X-ray spectrum are consistent with the prediction from ADAF (e.g. Esin et al. 1997). Tentative evidence for accretion mode transition has been reported in a few AGNs (e.g. Quataert et al. 1999;Yuan & Narayan 2004;Xu & Cao 2009;Xie et al. 2016). However, strong observational evidence for accretion mode transition in individual AGN has not yet been found, though it has been invoked to explain the large amplitude X-ray variability in a few AGNs (e.g. NGC 7589 in , NGC 7213 in Xie et al. 2016) and changing-look AGN (Noda & Done 2018).
Observational evidence for accretion mode transition in AGNs is important not only because it proves that the accretion theory can be universally applied to stellar and supermassive BH-accretion system, but also because it may yield important insights into many yet unclear problems in AGNs. For instance, the expected SED from various accretion models are dramatically different, thus accretion mode transition can potentially explain the changing-look phenomenon found in some AGNs (e.g. MacLeod et al. 2016;Ruan et al. 2016;Yang et al. 2018). Noda & Done (2018) found that the SED of the changing-look AGN Mrk 1018 is consistent with typical Seyfert 1 AGNs in the high state, while it can be modelled with ADAF in the low state. They thus proposed that accretion state transition can explain the changing-look phenomenon in AGNs. Moreover, such transition system can provide crucial clues, such as the critical accretion rate and the transition time-scale, to improve our understanding of the physical process that causes the transition.
NGC 7589, at a redshift of 0.0298, is optically classified as Seyfert 1.9/LINERs and a LLAGN based on the weak broad Hα line, and the BH mass is estimated to be ∼ 10 7 M ). Using the archival Einstein, ROSAT and XMM-Newton data,  found that this source showed large amplitude X-ray variability during the period from 1980 to 2001, i.e. the 0.5-2.4 keV flux changed by a factor of 10 on a time-sale of months to years. The lowest Eddington ratio 2 (the 1995 ROSAT pointing observation) for this source was conservatively estimated to be less than 10 −4 , suggesting that the source was accreting via the RIAF/ADAF mode at its low luminosity state. While the accretion mode was not clear in the relatively high luminosity state observed in 2001 by XMM-Newton, the estimated accretion rate was about a few percent, suggesting that the accretion may possibly be in a transition state between RIAF and the thin disc.  thus proposed that a transition of accretion state may have taken place in NGC 7589 if the X-ray flux had ever reached a peak much higher than the flux state observed in 2001. However, as noted by the authors, a partial covering model which may also explain the observed variability cannot be ruled out.
In this work, using all the new serendipitous X-ray observations of NGC 7589 from XMM-Newton and the X-ray telescope (Burrows et al. 2005) onboard the Neil Gehrels Swift Observatory (Swift/XRT), we show that indeed NGC 7589 can reach an even higher flux than the 2001 'high' luminosity state. A low luminosity state comparable to the 1995 ROSAT observation was also detected in the 2018 Swift/XRT observations. Results from the analysis of the X-ray spectra and multi-band variability suggest that the large amplitude variability is due to the change of intrinsic X-ray flux, rather than (partial covering) obscuration, and may indicate that a transition of accretion mode had taken place in NGC 7589. Throughout this paper, we adopted a flat ΛCDM cosmological model with H 0 = 69.3 km s −1 , Ω m = 0.29 and Ω Λ = 0.71. All the quoted uncertainties correspond to the 90 per cent confidence level for one interesting parameter, unless specified otherwise.

MULTI-BAND DATA REDUCTION
NGC 7589 had been observed in X-ray band by Einstein in 1985, ROSAT in 1992, XMM-Newton in 2001June and November, and 2006 June. It was also serendipitously observed by Swift/XRT extensively from 2006 April till 2018 July. In this work, only the X-ray data from XMM-Newton and Swift/XRT observations were analysed (please refer to  for the data analysis of X-ray observations taken before 2000). NGC 7589 is located in the sky region covered by the SDSS Stripe 82 Survey. Optical and UV observations were also frequently carried out by SDSS and GALEX in the past decades, respectively. In this section, we described the X-ray, optical, UV, and MIR multi-band data analysis. 2 λ Edd = L bol /L Edd , where L Edd = 1.3 × 10 38 (M BH /M ) and L bol is the bolometric luminosity.

XMM-Newton
The Observation Data Files (ODF) were downloaded from the XMM-Newton Science Archive. The ODFs were then reduced using the XMM-Newton Science Analysis System (SAS) software (version 16.1, Gabriel et al. 2004). The SAS tasks and were used to generate the event lists for the European Photon Imaging Camera (EPIC) MOS (Turner et al. 2001) and pn (Strüder et al. 2001) detectors, respectively. High background flaring periods were then identified and filtered from the event lists. A circular region with radius of 32 and 35 arcsec was selected as the source region for the two 2001 observations and the 2006 observations (the source was much brighter and off-axis), respectively. For the background region, an annulus (concentered with the source) with inner radius of 40 arcsec and outer radius of 120 arcsec were chosen for the 2001 EPIC MOS observations, while an circular region with radius of 70 arcsec was selected for the EPIC pn data (only available for the second 2001 XMM-Newton observation). We selected an annular region with inner radius of 60 arcsec and outer radius of 140 arcsec for the 2006 EPIC MOS1 observation. In the case of MOS2, an annulus is impossible, thus a circle with radius of 100 arcsec was chosen as the background region. The and tasks were used to generate the response files.
The OM data with the filter UVW1 were available for two observations (obsID:0066950301, 0066950401). There are four exposures (in total 2400 s) in the first observation and 10 exposures in the second observations (total exposure time of 8000 s). The SAS task was used to generate the OM images. The photometric measurements were performed for each exposure in the two observations with the task . A circular region with radius of 6 arcsec was selected for the source aperture, while an annulus with inner radius of 9 and outer radius of 14 was chosen to estimate the background. No significant variability was found within each observation, thus the mean magnitude of all the exposures in each observation was then used. The standard deviation was considered as the 1 σ uncertainty. The details of all the XMM-Newton observations can be found in Table 1.

Swift observations
NGC 7589 was serendipitously observed 56 times by the Swift/XRT from 2006 to 2018. The light curve as well as the X-ray spectra were generated using the XRT online data analysis tools 3 (Evans et al. 2009). The 3 σ flux upper limits (assuming an absorbed powerlaw model, see Section 3.2.1) were given for observations that the source was not detected. NGC 7589 was in a low X-ray luminosity state in 2018, resulted in low data quality of the X-ray spectra. To increase the S/N, we generated one combined X-ray spectrum using all the 2018 Swift/XRT observations. The details of the Swift/XRT observations can be found in Table 1 (note that all the observations taken in 2018 were combined and named as 2018 combined data).
The Swift/UVOT data were available for several observations (obsIDs: 00049538003, 00049538004, 00049538005, 00049538007, 00049538009, 00049538011, 00049538013, 00049538014, 00049538016). The source counts were extracted from a circular region with radius of 5 arcsec, while a 15 arcsec circle was choose as the background region. The task was used to perform photometric measurements. The filters and exposure times in each observation can be found in Table 1.

SDSS
NGC 7589 is located in the sky area covered by the SDSS Stripe 82 Survey, and was frequently observed in u, g, r, i, z bands with SDSS. We searched in the SDSS Stripe 82 Sky Survey at the position of NGC 7589 with a searching radius of 0.5 arcsec. This resulted in a total of 77 observations. The measured photometric magnitudes and observation time of the u, g, r, i, z bands for each individual observation were obtained from the SDSS Stripe 82 Sky Survey Database 4 . We noted that the contribution from host galaxy was not subtracted.

GALEX
NGC 7589 was observed by Galaxy Evolution Explorer (GALEX) 11 times with a total exposure time of 8874 seconds in the NUV band, and 8 times with 6886 second exposure time in the FUV band. The P (Million et al. 2016a,b) Python package 5 was adopted to analyze the GALEX NUV and FUV data. Intensity maps were generated for the two bands by running the M task, and were then used to select the source and background regions. A circle with radius of 7.4 arcsec, which excluded most of the radiation from the host galaxy, was selected as the source region for both NUV and FUV bands. An annulus with inner radius of 58 arcsec and outer radius of 68 arcsec was chosen as the background region for the FUV data. For the NUV data, a proper background region without contamination from a nearby source was not possible (only a concentered annulus background region is allowed in P ). Thus background subtraction was not applied for the NUV data. We noted that this will not affect our conclusion on the NUV variability. The A task was then used to calculate the AB magnitudes for both the NUV and FUV data.

WISE
The Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) is a satellite that surveys the sky in mid-IR bands. We searched the ALLWISE (Mainzer et al. 2011) and the NEOWISE Reactivation data release catalogues (Mainzer et al. 2014) at the source position of 4 http://cas.sdss.org/stripe82/en 5 https://archive.stsci.edu/prepds/gphoton NGC 7589 with a matching radius of 1 arcsec. This results in a total of 22 and 119 single exposure observations from the ALLWISE and NEOWISE, respectively. The W1 (3.4 µm) and W2 (4.6 µm) Vega magnitudes (contribution from the host galaxy was not subtracted) and uncertainties were obtained from the NASA/IPAC Infrared Science Archive (IRSA) 6 .

X-ray variability
The long-term intrinsic 0.5-10 keV flux ( f 0.5−10 keV ) variability of NGC 7589 is shown in the upper panel of Fig. 1. The f 0.5−10 keV were estimated by fitting the X-ray spectra with an absorbed powerlaw model (see section 3.2) for observations of which NGC 7589 was detected. While the 3σ upper limits 7 , calculated using the best-fitting parameters obtained by fitting the X-ray spectra of the nearest detection, were given for observations that the source was not detected.
From Fig. 1, it is clear that NGC 7589 showed large amplitude X-ray variability on a time-scale of ∼ 10 years. For instance, the highest observed X-ray luminosity (2006 April by Swift/XRT) is more than 50 times higher than the lowest X-ray luminosity observed in 1995 by ROSAT and 2018 by Swift/XRT. Flare-like variations, over a time-scale of around half a year, were also found, e.g. one at around 2001 June and the other one at 2006 April. During those flaring phases, the X-ray luminosity changed by a factor of 3. We note that flaring-like X-ray variability with time-scale of a few month may have also occurred in NGC 7589, as suggested from the multi-band optical/UV data (see Section 3.3).
The fractional rms variability amplitude F var (e.g. Edelson et al. 2002) is used to estimate the short-term (minutes to hours) variability. It is given by:  wherex is the mean counts rate, σ 2 err is the mean of the square of the 1σ errors, We calculated the fractional rms variability amplitude for the three XMM-Newton observations. No significant short-term variability is found in all the three observations, i.e. with F var less than a few per cent.

X-ray spectral analysis
The X-ray spectra were analysed using X (version 12.10, Arnaud 1996) with the Cash statistic (Cash 1979, wstat in X ). The energy band above 8 keV is dominated by background, we thus performed X-ray spectral modelling in the 0.3 − 8.0 keV energy range for all the data. For each of the XMM-Newton observations, we jointly fitted all the spectra available from the EPIC cameras. A normalisation factor was thus added to account for the calibration differences between the instruments. Galactic and host galaxy absorptions were included in all our analysis using the TB and TB model (Wilms et al. 2000, abundances are set to wilm in X ). The column density of the Galactic absorption was fixed at 3.84 × 10 20 cm −2 (Kalberla et al. 2005). We defined an absorbed power-law model (TB * TB * in X ) as our baseline model. In the following spectral analysis, we define three different luminosity states (see the upper panel of Fig. 1): the high luminosity state (L 0.5−10 keV > 10 43 erg s −1 ), the intermediate luminosity state (10 41.5 erg s −1 < L 0.5−10 keV < 10 43 erg s −1 ), and the low luminosity state (L 0.5−10 keV < 10 41.5 erg s −1 ). The two dot-dashed grey lines in the upper panel of Fig. 1 mark the boundaries of the three luminosity regimes. The Eddington ratios are calculated from the 2-10 keV X-ray luminosity with a bolometric correction of 28 for the high and intermediate luminosity state, and 16 for the low luminosity state (Ho 2008).

High luminosity state: the 2006 April Swift/XRT observation
Swift/XRT observed NGC 7589 in its highest X-ray flux state observed so far on 2006 April 30 (one snapshot in observation 00035365002). The signal-to-noise ratio (S/N) of the X-ray spectrum is low (total background subtracted photon counts: 29) due to the short exposure time (412 second). The best-fitting photon index is 1.31 +0.48 −0.47 when fitted with the baseline model. An upper limit of 5 × 10 20 cm −2 is estimated for the column density of the host galaxy. A significant excess in the soft X-ray band is shown in the ratio of data to the best-fitting baseline model (see the upper and bottom panels of Fig. 2). We thus fitted the spectrum with the baseline model plus a blackbody component. The likelihood ratio test ( ) method was adopted to estimate the significance of the blackbody component. To perform the , we simulated 10000 spectra using the best-fitting baseline model ( in X ). The change in C stat was then calculated by fitting the simulated spectra using both the baseline model and the model consisting of the baseline model and a blackbody component. The result suggested that the blackbody component is significantly detected at more than 3σ level (99.9 per cent). The best-fitting temperature of the blackbody component is 19 +15 −7 eV, while the best-fitting photon index is 1.03 +0.83 −0.77 . An upper limit of 6.5 × 10 21 cm −2 can be given for the host galaxy absorption. The predicted temperature of the inner region of a thin disc, e.g. 3 − 5 R S , close to the innermost stable circular orbit (ISCO), is ∼ 18 eV for a non-spin BH with mass of 10 7 M and accretion rate of 0.1, consistent with the best-fitting temperature of the blackbody component. The spectrum of a thin disc can generally be described as a multicolour blackbody, rather than a simple blackbody. We thus tried to fit the XRT spectrum with the model. The best-fitting temperature (18 +7 −9 eV) is consistent with that found in the simple blackbody model, though the normalization of the model can not be well constrained due to the low S/N and the narrow energy band. Simultaneously multiband observation 8 in the future could give a better constraint on the parameters in the model. For source with relatively low photon counts, the Swift/XRT spectrum could be affected by residual hot pixels or bright earth, both of which affect the low energy end of the spectrum (Andy Beardmore, private communication). No hot pixels are reported by the task, suggesting that the spectrum is not affected by residual hot pixels. Evidence for bright earth is found in the 2006 April Swift/XRT observation (Obs. ID: 00035365002). However, we note that NGC 7589 was observed in only one snapshot of this observation during which no sign of significant contamination from the bright earth was shown in the Swift/XRT image.
To further assess the potential contribution from bright earth to the Swift/XRT spectrum, we extracted a Swift/XRT spectrum with the grade 0 events only 9 . Except for the lowest energy bin  , and MIR (brown triangle for W1 and grey triangle for W2) light curves for NGC 7589.
(< 0.32 keV), which is slightly weaker (but still consistent within uncertainty), no obviously difference was found comparing to the grade 0-12 spectrum. A soft X-ray component which can be modelled by a blackbody component with temperature consistent with the value obtained from the grade 0-12 spectrum, is still detected at 3σ confidence level. Those results indicate that the soft excess shown in the 2006 April Swift/XRT spectrum is unlikely caused by the bright earth.

Intermediate luminosity state: the 2001 and 2006 XMM-Newton observations
The X-ray spectra of the first and second XMM-Newton observations can be well fitted with the baseline model. The best-fitting results are Γ = 1.68 +0.13 −0.08 and N H,host < 3 × 10 20 cm −2 for the first observation, Γ = 1.73 +0.17 −0.10 and N H,host < 5 × 10 20 cm −2 for the second observation, consistent with that reported in . The unabsorbed 0.5-10 keV X-ray flux estimated using the best-fitting baseline model for the first and second observations are 9.5 × 10 −13 and 3.5 × 10 −13 erg cm −2 s −1 , respectively. Evidence for an emission line profile at around 6.4 keV, which was not reported in  due to the use of binned data, was found in the first observation. We then fitted the X-ray spectra with a model consisting of the baseline model and a Gaussian component, i.e. TB * TB *( + ) in X . The line width σ of the Gaussian component was fixed at 1 eV, i.e. a narrow emission line. The best-fitting line energy of the emission line profile is 6.46 +0.10 −0.10 keV, consistent with the narrow neutral Fe Kα line that has been ubiquitously found in the X-ray spectra of AGN (e.g. Nandra et al. 2007;de La Calle Pérez et al. 2010). The equivalent width of this line is EW = 523 +619 −288 eV. The best-fitting results can be found in Table 2. the host galaxy absorption. The line energy is in agreement with the neutral Fe Kα line, i.e. E = 6.43 +0.05 −0.05 keV, and the EW of the line is 523 +362 −219 eV which is consistent with EW found in the the 2001 observation. The significance of this line is higher than 3σ as inferred from the method. Leaving the line width as a free parameter does not improve the fitting, and the best-fitting line width is σ < 115 eV, suggesting a narrow neutral Fe Kα line. The EWs of the narrow Fe Kα line, though with large uncertainties, are higher than the typical value (∼ 100 eV, e.g. Shu et al. 2010) found in broad line Seyfert galaxies , which may indicate a high Fe abundance or a large covering factor of the torus.
A soft X-ray excess component is often detected in the X-ray spectrum of AGNs (e.g. Arnaud et al. 1985;Gierliński & Done 2004). This component can normally be fitted with a blackbody with temperature of ∼ 150 eV (e.g. Page et al. 2004;Gierliński & Done 2004). The method is used to test whether a soft X-ray component is shown in the 2001 and 2006 XMM-Newton observations. We found that the soft X-ray component is not significantly (less than 2σ) detected in the 2006 XMM-Newton observation. Evidence for the soft X-ray excess was also not found in the 2001 XMM-Newton observations.

Low luminosity state: the 2018 Swift/XRT observations
As mentioned in section 3.1, the lowest X-ray flux was detected in the 1995 ROSAT and the combined 2018 Swift/XRT observations. Hereafter, we refer those observations as the low luminosity state. The combined 2018 Swift/XRT spectrum, with a total exposure time of 25 ks, can be fitted with the baseline model (see Table 2), though with large uncertainties due to low S/N (26 net source counts in the 0.3-8 keV band). The column density of the host galaxy cannot be well constrained, with a 90 per cent upper limit of 2.7 × 10 21 cm −2 . The best-fitting photon index is rather flat but with large uncertainty, i.e. Γ = 1.07 +0.67 −0.53 , which is still consistent with the photon index estimated at the intermediate and high X-ray luminosity states (see Section 3.2.2 and 3.2.1).  found that the X-ray spectra of the two 2001 XMM-Newton observations can also be well fitted with a partial covering model (see their table 2). We also tried to fit the XMM-Newton and Swift/XRT data with the partial covering model, i.e. TB *TB * in (or TB * * for ionized partial covering model).

Alternative model: partial covering
We first independently fitted the three XMM-Newton and the 2006 April XRT spectra. Both the column density and the covering factor parameters cannot be constrained for the 2018 XRT spectrum, thus this observation was not fitted with the partial covering model. We found that the (ionized) partial covering model can fit the Xray spectra. However, it does not improve the fitting significantly comparing to the simple baseline model and some of the parameters cannot be well constrained (only upper/lower limits can be given, see Table 3).
If indeed the observed X-ray variability was due to partial covering absorption, then we may expect that the intrinsic primary emission did not change significantly throughout the years covering all the X-ray observations. We thus fitted the three XMM-Newton and the combined 2018 Swift/XRT observations 10 simultaneously with the partial covering model. The parameters in the power-law model, i.e. the photon index Γ and the normalization, were linked for the X-ray spectra of the four observations, while the column density and the covering factor were fitted independently for each individual observation. The best-fitting results of the joint spectral analysis can be found in Table 3. The (ionized) partial covering model can simultaneously fit the X-ray data well. However, the inferred covering factors are relatively large, especially for the 2018 Swift/XRT observation of which a covering factor close to 1 is required.

Variability in optical, UV and mid-infrared
In the bottom panel of Fig. 1, we also show the long-term light curves in optical (light blue: SDSS u band, magenta: SDSS g bands), UV (GALEX data, lime and green hexagon for FUV and NUV, respectively), and mid-infrared (WISE, brown circle: W1, gray circle: W2). Large variability with ∆m > 0.5 mag is clearly seen in the mid-infrared to optical/UV bands.
The long-term SDSS u and g bands light curves are shown in 10 The parameters cannot be constrained for the 2006 XRT data when jointly fitted with all the other X-ray spectra. We thus did not include this observation in our joint X-ray spectral analysis. the lower panel of Fig. 1. The light curve in the r, i, and z bands are similar to the g band but with smaller amplitude of variability. From The characteristic of the NUV and FUV data are consistent with the optical light curve. As in the optical, the NUV and FUV magnitudes of NGC 7589 were relatively low in the 2003 observation, while they increased to an intermediate level during [2006][2007][2008]. In addition, a rapid rise (probably a flare) within 1.3 months is clearly shown in the 2009 observations. The source went into a low flux state in NUV in 2009. Weak variability (∆m < 0.2 mag) are also shown in the UV light curves measured from the Swift/UVOT and XMM/OM. We caution here that such a small magnitude of variability may caused by systematic uncertainties in the aperture photometric measurements, especially for the Swift/UVOT measurements of which the host galaxy may dominate the UV radiation as the source was in a very low luminosity state during the observations.
The WISE W1 and W2 mid-infrared magnitude measured at around the end of 2010 is about 0.2 mag higher than that observed half a year earlier. No significant variability is seen in the data taken after 2014, which are ∼ 0.5 mag dimmer than previous observations. Those data suggest that, though the sampling is not good due to the lack of data from 2011 to 2014, a 'flaring' like profile is probably also shown in the WISE W1 and W2 data. The W1 and W2 magnitudes did not vary significantly after 2014, which may be due to the fact that the radiation in the two bands were dominated by the host galaxy.
The 'averaged' luminosity is 1.5 × 10 41 erg s −1 in 0.5-10 keV band, as estimated from the measured narrow Hα luminosity using the Hα-X-ray luminosity relation ). The long-term multi-band light curves indicate that probably NGC 7589 normally behaves like a LLAGN, and underwent several outbursts during the last three decades.

DISCUSSION
Using only the 1995 ROSAT and 2001 XMM-Newton observations,  suggested that the large X-ray amplitude, with a change in X-ray flux by a factor of > 10, can be explained by intrinsic variability, i.e. the change of accretion rate. However, a partial covering mode can also explain the observed variability, and cannot be ruled out. In this work, using more X-ray data including new observations from XMM-Newton and Swift/XRT, we found that the NGC 7589 can reach an even higher X-ray luminosity, resulting in an increase in X-ray luminosity by a factor of more than 50 comparing to the 1995 ROSAT and the 2018 Swift/XRT observations. In addition, flare-like features are also shown in the optical, UV, and MIR long-term light curves, suggesting that NGC 7589 underwent several 'outburst' in the past decades.

Intrinsic variability vs partial covering absorption
Those new X-ray observations and multi-band data allowed us to further test the (ionized) partial covering scenario. In this scenario, we assumed that the flux of the primary emission was the same during different X-ray observations. We then jointly fitted the 2018 Swift/XRT observations (low luminosity state) with the 2001 and 2006 XMM-Newton observations (intermediate luminosity state). The (ionized) partial covering model can fit the X-ray data well (see Table 3). However, the inferred covering factors are relatively large, especially for the 2018 Swift/XRT observation of which a covering factor close to 1 is required. The result indicates that the primary emission was nearly fully obscured by an approximately Comptonthick absorber during the 2018 Swift/XRT observation. However, this is at odds with the best-fitting results obtained by fitting the X-ray spectrum of the combined 2018 XRT observations with the baseline model of which no evidence for heavily absorption is found, indicating that the primary emission is unlikely to be constant.
Further evidence against the partial covering model can be found from the long-term MIR variability. As mentioned before, NGC 7589 showed flare-like features in MIR bands (see Fig. 1 and Section 3.3). The MIR emission in AGN is believed to originate from the reprocessing of the primary emission (mainly the UV/optical radiation 11 ) by the dusty torus. In the partial covering scenario, the covering factor along the line-of-sight (LOS) can change dramatically (thus cause the X-ray variability). However, the global covering factor of the absorber (as well as the dusty torus) is not expected to vary significantly. So the fraction of the primary emission that will be reprocessed and re-radiated in the MIR band by the dusty torus should not show large amplitude variability in a relatively short time-scale. The significant variability found in the MIR data of NGC 7589 is then not expected in the partial covering scenario. Thus the large amplitude X-ray variability in NGC 7589 is unlikely caused by partial absorption.
As mentioned in Section 3.3, 'flare-like' features are also seen in the long-term UV/optical light curve. The sudden increase of the MIR flux in around 2011 was more likely to be attributed to the echo of a burst in the UV/optical band. Evidence for MIR echos had been found in the changing-look AGNs (e.g. Sheng et al. 2017), and in the TDEs (e.g. Dou et al. 2016Dou et al. , 2017Jiang et al. 2016Jiang et al. , 2017. We conclude that the observed multi-band variabilities were caused by the change of the intrinsic radiation rather than partial covering absorption.

Accretion mode transition in NGC 7589?
Theoretical calculation suggests that a transition in accretion mode will occur once the accretion rate reaches the critical values (Meyer et al. 2000a). Such a transition has been observed in Galactic XRBs (Tanaka et al. 1995;Zhang et al. 1997). The predicted time-scale for the accretion mode transition in AGN is much longer than that in XRBs, making it difficult to find such a transition system among AGNs. Conclusive evidence for accretion mode transition in AGN has not been found yet.
NGC 7589 is spectroscopically classified as Seyfert 1.9 or LINERS from the 2000 SDSS spectrum (see ). In the above subsection, we argue that the large amplitude X-ray variability was caused by the change of intrinsic radiation rather than absorption. The unabsorbed intrinsic 2-10 keV X-ray luminosity in the low luminosity state is then well below ∼ 10 42 erg s −1 , i.e. in the LLAGN regime.  proposed that the accretion proceeded via RIAF/ADAF during the low X-ray luminosity state observed by ROSAT in 1995, and that the large X-ray amplitude, with an increase by a factor of > 10 from the 1995 ROSAT to the 2001 XMM-Newton observations, can be explained by a transition in accretion mode, i.e. from a RIAF/ADAF to a standard thin disc or vice verse. During the 2001 XMM-Newton observations, the accretion flow in NGC 7589 was possibly in a transition state between ADAF and thin disc.  also predicted that a transition to thin disc might have taken place if the flux could reach a much higher peak value than the 2001 XMM-Newton observations. In this 11 A positive correlation between the UV and X-ray flux is found in some AGNs during the flare phase, e.g. NGC 1556 (Oknyansky et al. 2019) and Mrk 335 (Gallo et al. 2018). Thus a flare in the X-ray band may also imply a burst in the UV/optical band. work, using latest X-ray data from XMM-Newton and Swift/XRT, we found that NGC 7589 can indeed reach an even higher X-ray flux, resulting in an increase in X-ray flux by a factor of more than 50.
NGC 7589 was in a high X-ray luminosity state during the 2006 April Swift/XRT observation, with an estimated λ Edd of ∼ 0.13 (assuming a X-ray bolometric correction factor of 28, Ho 2008). The estimated λ Edd , which is > 100 times the λ Edd found in the 1995 low luminosity state, is well beyond the critical values predicted by theoretical calculation, indicating a thin disc mode in NGC 7589 during the 2006 April observation. Evidence for the emergence of a thin disc is further supported from the X-ray spectral analysis. A soft X-ray component which can be fitted with a blackbody model is significantly detected in the 2006 April Swift/XRT spectrum. Unlike the soft X-ray excess (normally has a temperature of ∼ 150 eV when modelled with a blackbody) that commonly found in Seyfert 1 galaxies, the best-fitting temperature of this component, which is ∼ 19 +15 −7 eV, is in agreement with the expected temperature of the inner region of a thin disc with BH mass of ∼ 10 7 M . Such a component has been reported for only a few sources so far, e.g. RXJ1643+49 (Yuan et al. 2010, see also Shu et al. 2017;Sun et al. 2013 for the other sources). We conclude that a transition from ADAF to the thin disc had taken place in the 2006 April observation.
NGC 7589 was observed to exhibit a low X-ray luminosity state, i.e. comparable to the 1995 ROSAT observation, during the 2018 Swift/XRT observations. The λ Edd is conservatively estimated to be ∼ 10 −3 (assuming a X-ray bolometric correction factor of 16, Ho 2008) using the 2-10 keV luminosity measured from the combined 2018 Swift/XRT observation. This may indicate an accretion flow via ADAF, suggesting a transition from the thin disc to ADAF. More evidences to support the accretion mode transition, such as the broad band SED in the low luminosity state, the reflection spectrum, and the broad Fe Kα line, might be found in future X-ray and multi-band observations with better data quality.
Our results may imply that accretion mode transition can indeed take place in AGN, supporting the idea proposed by . It also proves the idea that the accretion theory can apply to BH-accretion system across different BH mass scales. Future simultaneously multi-band observations of NGC 7589, especially optical spectroscopic data, will potentially help us understanding some important but yet unclear questions in AGN, such as the long-term variability of AGN (e.g Cao & Wang 2014) and the changing-look AGN.

Transition time-scale
The physical mechanism of the state transition from the low state to the high state (ADAF to thin disc) or vice versa (thin disc to ADAF) could be understood in the framework of the disc evaporation model (e.g. Liu et al. 1999;Meyer et al. 2000a,b;Liu & Taam 2009;Taam et al. 2012;Qiao et al. 2013). The disc evaporation model predicts an evaporation-curve, in which the evaporation rate increases with decreasing the radius until a maximum evaporation rate M max is reached and then the evaporation rate decreases with decreasing the radius. So in the disc evaporation model, if the mass accretion rate from the outer region of the disc is less than M max (∼ a few percent of the Eddington accretion rate M Edd = L Edd /c 2 , where L Edd = 1.3 × 10 38 M BH /M erg s −1 ), the disc will be truncated at a radius where the mass accretion rate equals the evaporation rate. From the truncation radius inwards, the accretion flow will be existed in the form of the ADAF. Then the geometry of the accretion flow will be an inner ADAF plus an outer truncated disc. Such a geometry of the accretion flow is often used to explain the spectral features of the low state of black hole X-ray binaries (e.g. Esin et al. 1997) and low-luminosity AGNs (Ho 2008;Nemmen et al. 2011). Whereas, if the mass accretion rate from the outer region of the disc is greater than M max , the disc cannot be completely evaporated at any radius, and the disc will extend down the ISCO of the black hole. Such a geometry of the accretion flow is widely used to explain the spectral features of the high state of black hole X-ray binaries and luminous AGNs.
The time-scale of the state transitions might be helpful in understanding the physical process (e.g. the disc evaporation model) for accretion mode transition. Assuming that the flare-like feature showed in the multi-band data were due to the change of intrinsic radiation, and a transition of accretion had taken place, we can then roughly estimate the transition time-scales: Low state to high state: the time-scale for the transition from low state to the high state, or from ADAF to thin disc (hereafter, t ADAF→TD ), estimated from the X-ray data alone is likely to be several months to less than 4 years, as suggested from the 2001 and 2006 X-ray observations. More stringent constraints on this transition time-scale could be obtained from the optical and UV data. As can be seen from Fig. 1, the variability time-scale is likely to be longer than ∼ 2 months and shorter than ∼ 1 year estimated from the optical (the raising phase in 2005) and NUV/FUV (the raising phase of the UV flare in 2008) data. However, the relatively small amplitude variability (∼ 0.5 mag) of those flares, comparing to the X-ray variability and the decline phase of the 2008 NUV flare (∼ 1 mag), may indicate that the observations did not cover the whole transition process, and/or the transition was not complete. The transition time-scale t ADAF→TD is then roughly estimated to be a few years (< 4 years) by combing the X-ray and optical/UV data. High state to low state: a time-scale of ∼ 12 years for the transition from high state to low state, or from thin disc to ADAF (hereafter, t TD→ADAF ), can be inferred from the long-term X-ray variability observed by Swift/XRT from 2006-2018. However, a much shorter time-scale can be inferred from the NUV data. From the bottom panel of Fig. 1, it is clear that the NUV flux peaked at around late 2008, then it gradually declined to a low flux state (comparable to the lowest NUV flux detected in late 2003) till late 2009, probable implying a transition to the ADAF accretion flow. If this is the case, then the transition time-scale t TD→ADAF could be as short as ∼ 1 year.

Transition from low state to high state
In the present work, initially when the source is in the low state, the disc is suggested to be truncated at a radius of a few hundreds of Schwarzschild radii as predicted by the disc evaporation model (Taam et al. 2012). The disc component detected in the 2006 April Swift/XRT observation suggests that the inner region of the accretion disc extends probably to the ISCO. The time-scale of the transition from the low state to the high state may correspond to the viscous time-scale of the truncated disc extending down to the ISCO of the black hole as the mass accretion rate in the disc exceeding M max . Such a viscous time-scale τ vis can generally be calculated as τ vis ∼ R/v R = R 2 /(αc s H) (Frank et al. 2002;Kato et al. 2008), where α is the poorly known viscosity parameter; H is the scale height of the accretion disc; v R and c s are the radial velocity and the sound speed at radius R, respectively. Following Cao & Wang (2014), we estimate the viscosity time-scale at a given truncation radius R tr using the following equation: where r tr = R tr /R g is the inner truncation radius (in units of gravitational radius) of the accretion disc, while H/R is the disc height to radius ratio which is normally less than 0.1 for a thin disc. By modeling the SED of a sample of LLAGNs using an accretion-jet model, Nemmen et al. (2011) suggested that the thin disk in LLAGNs is truncated at 60 − 450 R g . Assuming the disc is truncated at a few hundreds R g (e.g. 150 R g ) in NGC 7589, a H/R ∼ 0.05 is required so that the predicted viscosity time-scale (4 years for a disc with α = 0.5) matches with the estimated t ADAF→TD (a few years) for NGC 7589. Although such a large H/R is not expected in the thin disc model (H/R should be ∼ 10 −3 ), it can be obtained in the magnetic pressure supported accretion disc model proposed by (Dexter & Begelman 2019). In this model, the disc is geometrically thick at all luminosities which gives a viscous propagation time-scale as short as a few years, consistent with the observed time-scale. Alternatively, the thermal-viscous accretion disc instability model (DIM) (Meyer & Meyer-Hofmeister 1981), which is caused by the ionization of the hydrogen at the outer disc, has been successfully explained the transient outbursts observed in X-ray binaries and cataclysmic variables (CVs). DIM is proposed in (Noda & Done 2018) as one possible explanation for the luminosity change in Mrk 1018. In the DIM scenario, the jump of the temperature of the outer disc, triggered by ionization disc instability when the outer disc temperature is high enough to ionize hydrogen, will eventually results in a heating front propagating inward, leading to an increase in flux/luminosity. The time-scale is determined by the travelling speed of this heating front, which is shorter than the viscous timescale by a factor of H/R. However, as noted by (Noda & Done 2018) (see also the simulations in Hameury et al. 2009), the time-scale will still be too long to match the observed time-scale in NGC 7589 and Mrk 1018. A large H/R is still required to explain the observed time-scale.

Transition from high state to low state
The transitional time-scale from the high state to the low state can be estimated in the framework of the disc evaporation model as the mass accretion rate decreasing just below M max . When the mass accretion rate is just below M max , the disc will first be truncated at a critical radius R crit corresponding to M max , leaving an inner disc and an outer disc at the same time. The inner disc will be completely swallowed into the black hole within the viscous timescale. Such a viscous time-scale (hereafter τ dic ) is related with the critical truncation radius of the disc R crit , which can be from a few tens to a few hundreds of Schwarzschild radii depending on the viscosity parameter α, and the strength of the magnetic field (Qian et al. 2007;Taam et al. 2012). The estimated τ disc could be as short as ∼ 2 months for R crit = 20 R g , assuming α = 0.3 and H/R = 0.05 (again, such a large H/R can be obtained in the magnetic pressure supported disc model). Observationally, the disc component shown in the 2006 April Swift/XRT observation was not significantly detected (less than 2σ) in the XMM-Newton observation carried out one and half month later. This time-scale is consistent with the estimated τ disc for R crit = 20 R g , indicating a small critical truncation radius at M max for NGC 7589. Here we note that the α and H/R parameters can also have a large impact on the viscosity time-scale. A more precisely measurement of the transition time-scale in the future may be helpful to roughly constrain these parameters.
On the other hand, the inner radius of the outer disc (R in,outer ) may move towards to a much larger radii from R crit as the mass accretion rate decreases, providing that the evaporation rate is higher than the mass accretion rate at R in,outer of the outer disc. Eventually the outer disc will truncate at a radii (R tr,outer ) where the evaporation rate equals to the mass accretion rate. If this is the case, then the transition time-scale depends on not only the evaporation rate but also the decline curve of the mass accretion rate, both of which are subject to large uncertainties. Thus it is difficult to give a quantitate estimation of the high state to low state transition time-scale from the disc evaporation theory based on the current data of NGC 7589. Nevertheless, we may be able to test the disc evaporation model and constrain some of the key parameters in the model with better sampled multi-band light curves and more precisely measured transition time-scale in the future observations.
In the DIM scenario, the decrease in luminosity/flux can be explained as the propagation of the cool front, which is triggered when the temperature of the outer disc drops below 10 4 K, i.e. the ionized disc becomes neutral. Again, similar to the time-scale for the propagation of the heating front, a large H/R is required to explain the observed decline time-scale.

The soft X-ray excess
The soft X-ray excess has been commonly found in the X-ray spectra of AGNs. However, its origin is still unclear. The temperature of this component, when modelled with thermal Comptonization model, is ∼ 0.12 keV (Gierliński & Done 2004) which is much higher than the predicted temperature of the accretion disc with BH mass > 10 6 M but is consistent with thermal Comptonization of a warm corona (Magdziarz et al. 1998). Evidence for a soft X-ray component originated from a thin accretion disc, however, has been found in a few AGNs (Yuan et al. 2010;Sun et al. 2013;Shu et al. 2017). The temperature of the soft X-ray excess in NGC 7589 is ∼ 19 eV (modelled with a blackbody) in the high luminosity state (with λ Edd higher than 0.1), suggesting that this component was from the accretion disc during the 2006 April Swift/XRT observation. However, no significantly soft X-ray excess is detected in the 2001 (June and November) and 2006 (June) XMM-Newton observations during which the source was in an intermediate luminosity state, implying that the soft X-ray excess component could disappear within two months in NGC 7589.
The property of the soft X-ray excess for NGC 7589 is very similar to the results reported in Noda & Done (2018) for the changinglook AGN Mrk 1018. The soft X-ray excess component in Mrk 1018 also emerged in the bright state (λ Edd > 0.02), and it disappeared in the faint state (λ Edd > 0.02). Unlike NGC 7589, the temperature of the soft X-ray component in the bright state in Mrk 1018 is consistent with the Comptonization of a warm corona. While the emergence and disappearing of the soft X-ray component in Mrk 1018 is explained as the increase and decrease of the size of the warm corona region (Noda & Done 2018), respectively, the non-detection of the soft X-ray component in NGC 7589 in the intermediate state can be explained by the truncation of the inner accretion disc under the framework of the disc evaporation as discussed in Section 4.3. Whether a warm corona has ever formed in NGC 7589 is unknown due to the lack of observational data.
It is worth noting that the property of the soft X-ray excess also changed within two months for Mrk 1018 as reported in (Noda & Done 2018). These results indicate that the soft X-ray excess in AGNs can be a transient phenomenon. The most impressive evidence for the transition/variability of the soft X-ray excess component came from the recently discovered quasi periodic eruptions in GSN 069 (Miniutti et al. 2019, see also RX J1301.9+2747, Sun et al. 2013;Miniutti et al. 2019). During the late 2018 to early 2019 X-ray observations, a cold soft X-ray component with temperature of ∼ 50 eV is shown in the quiescent state of GSN 069, while a warm soft X-ray component with temperature of ∼ 120 eV is found at the peak of the eruptions which last for several kilo-seconds with an recurrence time of ∼ 30 ks. Moreover, the cold component can increase its temperature to ∼ 80 eV within one month. We note that the physical mechanism that triggers the emergence of the warm component and the transition between the cold and warm components may be very different in GSN 069, as the time-scale is much shorter than that observed in Mrk 1018 and NGC 7589. Nevertheless, these results indicate that there may be a diversity in the origin of the soft X-ray excess. A detail temporal analysis of the soft X-ray excess, on both short and long time-scales, may provide important insights on the origin of the soft X-ray excess and the mechanism triggering the accretion mode transition.

SUMMARY
In this paper, we report the discovery of large amplitude X-ray variability, i.e. by a factor of more than 50, of the LLAGN NGC 7589 using the archival X-ray data span several decades. The long-term X-ray, optical, and UV variability (see Fig. 1) indicate that NGC 7589 may have undergone several outbursts in the past decades. A detail analysis of the X-ray spectra in different X-ray luminosity state, combining with the MIR variability, suggests that the large amplitude X-ray variability is likely caused by the change of the intrinsic X-ray flux, rather than due to the (partial) absorption. Possible evidence for accretion mode transition in NGC 7589 is found by modelling the X-ray spectra at different luminosity state. At low X-ray luminosity state (during the 2015 ROSAT and 2018 Swift/XRT observations), the λ Edd is estimated to be less than a few 10 −3 , implying accretion via ADAF. While λ Edd increases by a factor of 100, to ∼ 0.13, at the high luminosity state, suggesting a thin accretion disc. Further evidence for a thin disc is supported by the significant detection of a soft component in the X-ray spectrum in the high luminosity state. The temperature (∼ 19 +15 −7 eV) of this component is consistent with the predicted temperature of the inner region of a thin disc around a ∼ 10 7 M . The time-scales of the accretion mode transitions are estimated to be several months to ∼ 1 years from the multi-band data. Our results may indicate that the accretion mode transition can take place in AGN over a timescale of months to a few years. Future intensive multi-band monitor of NGC 7589 can possibly give a much better estimation of the timescales at different transition stage which may help us understanding the physical process driven the accretion mode transition.