format: text

text:
# Quantum effects after decoherence in a quenched phase transition

## I Introduction

The emergence of classical behaviour in quantum systems is a topic of great interest for both conceptual and experimental reasons . It is well established by now that the interaction between a quantum system and an external environment can lead to its classicalization; decoherence and the occurrence of classical correlations being the main features of this process (for a recent overview see ).

A seemingly unrelated physical problem where the interaction between a main system and its surrounding environment is central is in determining the dynamics of a phase transition. Usually, a change in the properties of the system or the bath, forces the system to change phase via an out-equilibrium evolution. It is natural to ask what role decoherence plays in the phase transition and conversely, how the time dependent nature of the process affects the classicalization of the system.

In this article we explore two concurrent avenues. We look at what may happen with the decoherence process when we have a time dependent setting (so far this problem has been mostly studied in kicked or driven systems; see for example Refs. and ). This is a very general question, and we use to guide us a simple toy model that naturally includes time dependent features. This model also happens to mimic some properties of a non-equilibrium second-order phase transition, giving us some clues as to what may happen in a realistic case.

The paper is organised as follows. In the next section we introduce our model and review the physical role of the different terms in the relevant evolution equations. We describe how the ‘phase transition’ is implemented and discuss estimates for the different time scales involved. In Section III we present the results of a series of numerical simulations for the evolution of an initial configuration of two de-localised Gaussian wave packets. This system is subject to a sudden quench via an instantaneous change in the frequency sign. Both the cases where the temperature of the environment is kept fixed and allowed to change at the quench time are studied. We support the numerical results with a detailed analytical analysis. Section IV contains similar results this time taking as initial condition a single Gaussian state centred at the global minimum. In Section V we discuss the time dependent evolution of the linear entropy in the model, illustrating the loss of purity of the system and clarifying the physical nature of the results previously obtained. Section VI contains final remarks and the main conclusions of the paper.

## II The model

We will start by considering a quantum anharmonic oscillator coupled to an environment composed of an infinite set of harmonic oscillators. The total classical action for the system is given by:

$`S[x,q_n]`$ $`=`$ $`S[x]+S[q_n]+S_{\mathrm{int}}[x,q_n]`$ (1)
$`=`$ $`{\displaystyle _0^t}ds[{\displaystyle \frac{1}{2}}M(\dot{x}^2\mathrm{\Omega }_0^2(t)x^2{\displaystyle \frac{\lambda }{4}}x^4)`$ (2)
$`+`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{1}{2}}m_n(\dot{q}_n^2\omega _n^2q_n^2)]{\displaystyle }_nC_nxq_n,`$ (3)
where $`x`$ and $`q_n`$ are the coordinates of the particle and the oscillators respectively. The quantum anharmonic oscillator is coupled linearly to each oscillator in the bath with strength $`C_n`$. This coupling leads to a simple quantum Brownian motion (QBM) model commonly used in the study of the quantum to classical transition . Tracing over the degrees of freedom of the environment one obtains a master equation for the reduced density matrix of the system. From this one can derive the following evolution equation for the corresponding Wigner function :

$`\dot{W_\mathrm{r}}(x,p,t)=\{H_{\mathrm{syst}},W_\mathrm{r}\}_{\mathrm{PB}}{\displaystyle \frac{\lambda }{4}}x_{ppp}^3W_\mathrm{r}`$ (4)
$`+`$ $`2\gamma (t)_p(pW_\mathrm{r})+D(t)_{pp}^2W_\mathrm{r}f(t)_{px}^2W_\mathrm{r},`$ (5)
where

$`\gamma (t)`$ $`=`$ $`{\displaystyle \frac{1}{2M\mathrm{\Omega }_0}}{\displaystyle _0^t}𝑑t^{}\mathrm{sinh}(\mathrm{\Omega }_0t^{})\eta (t^{})`$ (6)
$`D(t)`$ $`=`$ $`{\displaystyle _0^t}𝑑t^{}\mathrm{cosh}(\mathrm{\Omega }_0t^{})\nu (t^{})`$ (7)
$`f(t)`$ $`=`$ $`{\displaystyle \frac{1}{M\mathrm{\Omega }_0}}{\displaystyle _0^t}𝑑t^{}\mathrm{sinh}(\mathrm{\Omega }_0t^{})\eta (t^{}),`$ (8)
$`\gamma (t)`$ is the dissipation coefficient, $`D(t)`$ and $`f(t)`$ are the diffusion coefficients. $`\eta (t)`$ and $`\nu (t)`$, the dissipation and noise kernels, are given respectively by:

$`\eta (t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑\omega I(\omega )\mathrm{sin}\omega t`$ (9)
$`\nu (t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑\omega I(\omega )\mathrm{coth}{\displaystyle \frac{\beta \omega }{2}}\mathrm{cos}\omega t,`$ (10)
where $`I(\omega )`$ is the spectral density of the environment.

The first term on the right-hand side of Eq.(5) is the Poisson bracket, corresponding to the usual classical evolution. The second term includes the quantum correction (we have set $`\mathrm{}=1`$). The last three terms describe dissipation and diffusion effects due to coupling to the environment. In order to simplify the problem, we consider a high-temperature ohmic ($`I(\omega )\omega `$) environment. In this approximation the coefficients in Eq.(5) become constants: $`\gamma (t)=\gamma _0`$, $`f1/T`$, and $`D=2\gamma _0k_\mathrm{B}T`$. The normal diffusion coefficient $`D`$ is the term responsible for decoherence effects and at high temperatures is much larger than $`\gamma _0`$ and $`f`$. Therefore in Eq.(5), we may neglect the dissipation and the anomalous diffusion terms against the normal diffusion. It is important to note that the high-temperature approximation is well defined only after a time-scale of the order of $`1/(k_\mathrm{B}T)\gamma _0/D`$ (with $`\mathrm{}=1`$). The relevant period of evolution for our systems takes place at times comfortably larger than this time-scale, safely in the validity regime of the approximation.

Time dependence will be introduced in the Hamiltonian by imposing a sudden change of sign of $`\mathrm{\Omega }_0^2`$ (typical quench). This mass term is taken to be positive initially, the original symmetry being broken by $`\mathrm{\Omega }^2`$ becoming negative. On a second stage we will also consider the case where the temperature of the environment $`T`$ changes with time. The change in the potential leads to the formation of degenerate minima mimicking the breaking of symmetry in a second order-phase transition. In a realistic model one should address this problem in the context of quantum field theory . This is an extremely difficult problem since non-perturbative and non-Gaussian effects are relevant in the dynamical evolution of the order parameter undergoing the transition and clearly numerical simulations are out of the question. We trust that any non-trivial type of behaviour that may be a feature of our simple quantum mechanical model will also be present (and likely more strongly so) in the infinite dimensional case.

## III De-localised initial states: quantum effects after decoherence

We solve Eq.(5) numerically using a fourth-order spectral algorithm (numerical checks included carrying out simulations at different spatial and temporal resolutions). We chose $`\lambda =0.1`$, $`D=0.3`$ and set $`\mathrm{\Omega }_0^2=1.0`$ initially. In order to understand the effects of the change in the mass term on the decoherence process we look first at the evolution of the quantum superposition of two Gaussian wave packets:

$$\mathrm{\Psi }(x,t=0)=\mathrm{\Psi }_1(x)+\mathrm{\Psi }_2(x),$$
(11)
where

$$\mathrm{\Psi }_{1,2}(x)=N(t)\mathrm{exp}\left[\frac{(xL_0)^2}{2\delta ^2}\right]\mathrm{exp}(\pm iP_0x).$$
(12)
The initial $`W_\mathrm{r}`$ consists of two Gaussian peaks $`W_\mathrm{r}^{1,2}`$ separated by a distance $`L_0`$ (we chose $`L_0=2.0`$ and $`P_0=0`$) and an interference term $`W_\mathrm{r}^{\mathrm{int}}`$. This quantum initial state has been widely used in the literature to illustrate decoherence phenomena (see  or  for example) and its evolution will make clear the physical nature of the effects we will observe. In the next Section we will chose a more realistic initial condition in terms of the dynamics of a phase transition.

In order to visualise deviations from classicality effectively, we define the auxiliary quantity ,

$$\mathrm{\Gamma }(t)=𝑑x𝑑p\left[|W_\mathrm{r}|W_\mathrm{r}\right].$$
(13)
When the Wigner distribution is positive and possibly identifiable with a classical probability distribution, $`\mathrm{\Gamma }`$ is zero. However if $`\mathrm{\Gamma }`$ is positive $`W_\mathrm{r}`$ must have negative values due to quantum interference terms. We can thus use positivity of the $`\mathrm{\Gamma }`$ function as sufficient condition for non-classical behaviour.

### A Mass quench at constant temperature

We start the simulation by evolving $`W_\mathrm{r}`$ for some time with the positive mass squared potential, in the presence of the bath. During this period, the initial quantum interference terms are quickly damped by the environment. Thus, for an early time $`t_{\mathrm{D}_1}`$, the system decoheres and one is able to distinguish two classical probability distributions corresponding to the two initial Gaussian peaks evolving over phase space. Suddenly, at $`t=t_\mathrm{c}`$ we change the frequency of the system from the initial positive value $`\mathrm{\Omega }_0^2`$ to a final $`\mathrm{\Omega }_0^2`$. The evolution picture changes dramatically when the frequency becomes negative and instabilities are introduced in the system. In Fig.1 we can see the behaviour of $`\mathrm{\Gamma }`$. Starting from a large initial value, $`\mathrm{\Gamma }`$ quickly tends to zero as quantum fluctuations vanish and the systems becomes classical. The potential is quenched at $`t_c`$ and shortly after the system displays once again quantum behaviour for a period of time.

In order to understand this process we go back to early times, before the quench. From $`t=0`$ up to $`t=t_\mathrm{c}`$ the diffusion coefficient $`D`$ causes the system to decohere, destroying quantum interference terms in a time that can be estimated to be of the order of $`t_{\mathrm{D}_1}1/(4L_0^2D)`$, where $`L_0`$ is the initial space separation between the peaks of the Gaussian wave packets (see ). The normal diffusion term is dominant with respect to the quantum corrections, and thereafter the evolution is given essentially by the classical Fokker-Plank flow. For our choice of initial conditions we have $`t_{\mathrm{D}_1}0.2`$. This is roughly the time quantum interference terms in the Wigner function should fall to $`1/e`$ of their initial value (we have checked that this is compatible with the decay of $`\mathrm{\Gamma }`$ in the initial period of evolution in our simulations). As soon as the frequency becomes negative, an unstable point forms in the centre of the phase space with associated stable and unstable directions. These are characterised by Lyapunov coefficients $`\mathrm{\Lambda }`$ with negative and positive real parts respectively .

The new type of dynamics gives rise to the possibility of squeezing along the stable direction. The exponential stretching of the Gaussian packets in one of the directions due to the hyperbolic point is compensated by an exponential squeezing. This will lead to a growth of gradients in the Wigner function that will make the quantum term in Eq.(5) comparable to the others. As a consequence the system will be forced to explore the quantum regime again. In a more quantitative fashion we have that the time dependence of the package width in the direction of the momenta after the quench is given by $`\sigma _p(t)=\sigma _p(t_\mathrm{c})\mathrm{exp}[\mathrm{\Lambda }(tt_\mathrm{c})]`$, where $`\sigma _p(t_\mathrm{c})`$ is the corresponding width at the time in which $`\mathrm{\Omega }_0`$ changes sign. From this we can estimate the $`p`$derivatives of the Wigner function to grow as $`_p^nW_\mathrm{r}\sigma _p^n(t_\mathrm{c})\mathrm{exp}[n\mathrm{\Lambda }(tt_\mathrm{c})]W_\mathrm{r}`$. Clearly higher order derivatives grow faster and at some point the quantum term with its third order derivative will be of comparable magnitude to the classical terms in the Poisson brackets (which are first order). This will happen (see ) when the ratio $`_p^3W_\mathrm{r}/_pW_\mathrm{r}`$ becomes of the order of $`\chi ^2=_xV(x)/_x^3V(x)\mathrm{\Omega }_0^2/\lambda `$ which characterises the scale of nonlinear terms. From this the time at which quantum effects become relevant is calculated to be

$$t_\chi t_\mathrm{c}+\mathrm{\Lambda }^1\mathrm{ln}[\chi \sigma _p(t_\mathrm{c})].$$
(14)
In the simulation used in our example we chose $`t_\mathrm{c}=2.5`$ (later than the time when the Wigner function becomes definite positive). We evaluate $`\chi 3.2`$ and numerically estimate $`\sigma _p(t_\mathrm{c})2.7`$. We are also assuming the Lyapunov coefficient is given by the value corresponding to a linear potential $`\mathrm{\Lambda }=2\mathrm{\Omega }_0^2=2.0`$. Therefore, the time in which quantum effects start being relevant is given by $`t_\chi 3.8`$. This is in good agreement with the time at which the Wigner function displays negative values once again, as can be seen in Fig.1.

From this point onwards quantum contributions increase, their growth being limited by diffusion effects which limit the squeezing of the Wigner function. The bound on the width of the packs is given by $`\sigma _\mathrm{c}=\sqrt{2D/\mathrm{\Lambda }}`$ . We use this to estimate the second decoherence time scale. We assume that quantum effects become maximal at a certain $`t_{\mathrm{max}}`$ (when in the numerical simulation $`\mathrm{\Gamma }`$ reaches its maximum) with a corresponding pack width $`\sigma _p(t_{\mathrm{max}})`$ and that decoherence is effective after the time when squeezing becomes of the order of the limiting value. This implies

$$t_{\mathrm{D}_2}=t_{\mathrm{max}}+\mathrm{\Lambda }^1\mathrm{ln}[\sigma _p(t_{\mathrm{max}})/\sigma _\mathrm{c}],$$
(15)
which defines the decoherence time after the critical time. Using $`\sigma _p(t_{\mathrm{max}})4`$ and $`\sigma _\mathrm{c}=0.5`$ we obtain $`t_{\mathrm{D}_2}6.3`$, in reasonable agreement with the simulation time for which quantum effects are exponentially suppressed (see Fig.1).

### B Mass quench with changing temperature

The pattern of classical-quantum-classical behaviour found in the above system with explicit time dependence is observed in more generic situations. As a second example we have solved Eq.(5) allowing the bath temperature to decrease simultaneously with the change in sign of the frequency term. These conditions take us somehow closer to what would happen in a true second-order phase transition caused by a temperature quench. As a consequence, the diffusion coefficient, proportional to $`T`$, goes at $`t_\mathrm{c}`$ from an initial high temperature value $`D_0`$ up to a final lower value $`D_\mathrm{f}`$ (still in the high temperature regime in order to ensure the validity of Eq.(5)). In Fig.2 we see the effect of changing the temperature with the classical potential (except for $`D`$ all simulations parameters are the same as in Fig.1). The analysis used in the previous example can be easily reproduced for this case. Both the initial decoherence time $`t_{\mathrm{D}_1}`$ and the time for the re-introduction of the quantum fluctuations $`t_\chi `$ remain unchanged as they do not depend on the temperature of the environment. The second decoherence time $`t_{\mathrm{D}_2}`$ is larger for a weaker diffusion term (we have used $`D_\mathrm{f}=0.1`$ in Fig.2-a and $`D_\mathrm{f}=0.003`$ in Fig.2-b). We have obtained respectively $`t_{\mathrm{D}_{2\mathrm{a}}}6.5.`$ and $`t_{\mathrm{D}_{2\mathrm{b}}}7.4`$. In the lowest temperature case (Fig.2-b) the analytical prediction matches the numerical result poorly.

This is due to the fact the estimation does not take into account the oscillations in the rate of decoherence coming from different orientations of the interference fringes when the Wigner function is moving around the unstable point. As the diffusion coefficient is smaller, the second decoherence time grows and the approximation of the upside-down potential in no longer valid. In any case, the analytic result can still be used as an estimated lower limit for the second decoherence time. We have included it in our analysis in order to emphasise how dramatic the quantum effects are during the quenched transition.

It is helpful to look at the Wigner function directly in order to further clarify which regions of phase space are responsible for turning $`\mathrm{\Gamma }`$ positive. In Figs.3 we show $`W(x,p,t)`$ for the quench case corresponding to Fig.2-a.

The four plots in the left column correspond to the decoherence period before the quench. The two Gaussian peaks (light spots) rotate in phase space around the minimum of the potential while the negative components (dark patches) of the Wigner function are cleared away by the environment. When the potential changes (right column) the wave packets start spreading and exploring the new non-linear regions of phase space giving rise to the dark interference patches. For longer times decoherence takes over again and the Wigner function becomes once more positively defined.

## IV Single initial Gaussian state

As a further example we take a single Gaussian state centred at the global minimum of the quartic potential as initial condition. This is a more reasonable initial condition in terms of a realistic phase-transition, mimicking a high-temperature thermal distribution. It will also allow us to see that the above results are not an artifact of the initial state. This initial Wigner function is already classical and so we ignore the initial evolution period and take $`t_\mathrm{c}=0`$. Fig.4 and Fig.5 show the $`\mathrm{\Gamma }`$ function for the same quenches as before (without and with temperature change respectively). The initial classical configuration ($`\mathrm{\Gamma }=0.`$ for the initial time) develops quantum effects as the classical potential and the temperature change. The relevant time scales are evaluated as before and once again, the estimates are in good agreement with the simulation results. In the constant temperature case $`\sigma _p(t_\mathrm{c}=0)0.7`$ which gives $`t_\chi 0.6`$ (see Fig.4). We also have $`\sigma _p(t_{\mathrm{max}})3.2.`$ and $`\sigma _\mathrm{c}0.5`$ leading to $`t_\mathrm{D}3.2`$, which agrees with the numerical result.

Fig.5 shows the cases where the change in frequency is followed by a change in the environmental temperature (same coefficients as in the example of Fig.2). For Fig.5-a $`\sigma _p(t_{\mathrm{max}})3.3`$ and $`\sigma _\mathrm{c}=0.3`$, and therefore the decoherence time is $`t_\mathrm{D}3.5`$. This scale is in good agreement with the numerical result. The estimation for Fig.5-b gives a decoherence time $`t_\mathrm{D}4.5`$ which again (as in the case of Fig.2-b) fails to fit the numerical result. We have included it in our analysis in order to emphasise how dramatic the quantum effects are during the quenched transition.

## V Linear entropy

One of the most salient features of the quantum to classical transition concerns the production of entropy as a consequence of the entangling interactions between the system and the environment. In order to clarify the nature of the post-decoherence quantum effects in the systems simulated above we have looked at the corresponding time evolution of the linear entropy. This is given in terms of the density matrix by :

$$S_l(t)=\mathrm{ln}[\text{Tr}[\rho ^2(t)]].$$
(16)
This quantity can be easily obtained from the Wigner function giving a good measurement of the ‘loss of purity’ of the system as it interacts with the bath (see ). We found that as expected the entropy increases throughout the whole evolution. The system starts as a pure state and while interacting with the heat bath it looses coherence and simultaneously starts behaving as a classical ensemble. When the potential changes it evolves for some time as a quantum mixed system but the original ‘purity’ is never recovered. In this sense the decoherence process is irreversible. In terms of the Wigner function the linear entropy is related to the area of its non-zero component in phase space. Due to the coupling to the environment the total area is not conserved, the Wigner function keeps spreading at all times leading to permanent growth of the entropy.

Fig.6 shows the time dependent linear entropy (top plot) and its production rate (bottom plot) for the two different initial conditions considered before, in a quench with fixed environment temperature.

In the case of the double-gaussian initial state (solid line) there is an initial period of evolution up to the first decoherence time $`t_{D_1}=0.2`$, where the linear entropy grows as a consequence of diffusion effects (as during the whole evolution) and also due to the disapearence of initial interference terms which are washed away by the environment. As these vanish the entropy production rate decreases as can be seen in the bottom plot. After an oscillation in the rate caused by the rotation of the Wigner function in the phase space which generates some low amplitude interference terms, the rate reaches a minimal value near the quench time at $`t=t_\mathrm{c}`$. After $`t_\mathrm{c}`$, the entropy rate starts growing again as the system gets rid of the newly induced interference terms. Finally at $`t_{D_2}`$, the entropy rate decreases to a low, slow decaying value driven by diffusion only.

The single gaussian evolution (dashed line) confirms this picture. From a low initial value (the initial state is free from negative terms) the entropy production rate grows as the quench generates interferences. Later, the environment cleans them out leading to the final decaying rate.

We should stress that a growing linear entropy function does not imply classicality (positivity of the Wigner function is an extra necessary condition in order to have a classical probability distribution). Increase of $`S_l`$ tells us that the pure initial quantum state is evolving into a mixed state. It does not of course, tell us whether this mixed state is a classical or quantum one. In particular this is the case between the quench time and the second decoherence time. During this period quantum effects are re-introduced while the linear entropy is still growing (faster even since its production rate increases).

## VI Final remarks

We have shown, using an exact numerical evaluation of the Wigner function that quantum effects can be re-introduced after decoherence in several systems with explicit time dependence. These quantum effects are originated when the changing dynamics introduce instabilities in previously stable regions of the phase space. When this happens the dynamics of the Wigner function becomes more relevant than the decoherence effects due to the environment (and the lowest the final bath temperature the more dominant these are). The system then displays quantum behaviour for a length of time until the environment manages to catch up and force classicalization once again.

Since all examples so far were based on systems described by a double well potential one could wonder whether our results could be a consequence of possible tunnelling phenomena between the two minima. Tunnelling is possible between symmetry related eigenstates with energy below the barrier. The tunnelling time-scale for each pair is well known to be inversely proportional to the energy splitting of the symmetry related pair of eigenstates. For the parameters of our system ($`\mathrm{\Omega }^2=1`$, $`\lambda =0.1`$, $`\mathrm{}=1`$) only seven pairs of states are found below the barrier. Their energy splittings range from $`\mathrm{\Delta }E_010^{12}`$ to $`\mathrm{\Delta }E_710^2`$, and thus the tunnelling would firstly be expected after $`t100`$. Therefore and considering the time scales in which our simulations take place, tunnelling should play no role. The stretching and folding of the Wigner function responsible for the observed effects happens on both ‘sides’ of the potential well independently. This is in agreement with the conclusion invariably found in the literature (see for example Ref.) that tunnelling takes place rather slowly when compared with all natural time-scales in the system.

In order to confirm directly that tunnelling phenomena are not responsible for the effects observed, we solved numerically the problem of a single Gaussian packet cantered at $`x_0=2`$ evolving in the usual quartic potential but with its motion restricted to $`x>0`$. The resulting $`\mathrm{\Gamma }`$ is shown in Fig.7. As before quantum behaviour is swiftly recovered, the corresponding time scales being in good agreement with the analytical estimates.

Our results open up several interesting possibilities. The most obvious one would be to try to ‘maximise’ the recovering of quantum effects to the extent of making them effectively permanent. An oscillatory frequency  that would continuously force instabilities into the system could prevent classicalization or at least postpone it for a great length of time.

In terms of the specific case of the dynamics of a second-order phase transition one could expect quantum effects to be present. Though the model used is a crude simplification of what happens in a realistic phase transition, the same features of time-dependent introduction of non-linearities would be present in that case, leading to similar, probably stronger quantum effects. Critical properties of infinite dimensional systems such as critical slowing down could play an interesting role in the process.

###### Acknowledgements.

We would like to thank S. Habib, F. D. Mazzitelli, and J. P. Paz for comments and useful discussions. The work of N. D. A. was supported by a PPARC Postdoctoral Fellowship and F. C. L. was supported by CONICET and Fundación Antorchas.


text:
# BeppoSAX Observations of the Seyfert-2 Galaxies NGC 7172 and ESO 103-G35

## 1 Introduction

In recent years X-ray missions such as EXOSAT Ginga and ASCA have demonstrated that Seyfert-2 galaxies have a very complex X-ray spectrum (see Mushotzky, Pounds & Done 1993 for a recent review). In particular their spectrum contains at least the following features: a) a power-law continuum with slope $`\mathrm{\Gamma }1.9`$ absorbed by a very large column ($`>10^{23}\mathrm{cm}^2`$) (eg Turner & Pounds 1988); this obscuring screen is probably associated with a molecular torus. b) an FeK line at 6.4 keV ie originating from neutral Fe c) a spectral curvature (hump) above 10 keV; this is generally believed to originate from reflection of the power-law component on cold material (eg Lightman & White 1988, George & Fabian 1991).

In contrast to the recent progress in understanding the spectral features of Seyfert-2 galaxies, our knowledge of their X-ray variability properties remains scanty. Usually, this has been based on comparison of the X-ray flux between observations from different missions spanning time intervals of several years. Recently, Turner et al. (1997) have detected some evidence for short term variability in a handful of Seyfert-2 galaxies using ASCA . The systematic monitoring of Seyfert-2 galaxies has become possible with the RXTE mission. RXTE observations of a few objects revealed strong variability in time-scales of hours as well as evidence for spectral variability (Georgantopoulos et al. 1999, Georgantopoulos & Papadakis 2000). Indeed, time variability studies can provide strong constraints in understanding the geometry of the circumnuclear matter in these objects. For example, the time lags between the power-law continuum and the normalization of the Fe line could constrain the origin of the line. In the case where the Fe line originates in an accretion disk (as is the standard scenario for Seyfert-1 galaxies), the Fe line should closely track the variations of the intrinsic power-law.

In this paper we present an analysis of four observations of the Seyfert-2 galaxies NGC 7172 and ESO 103-G35 from the BeppoSAX archive. The two observations for each galaxy are separated by a year’s interval. Our goal is to study any long-term variations in the spectral features between the two epochs. The broad energy band of $`\mathrm{𝐵𝑒𝑝𝑝𝑜𝑆𝐴𝑋}`$ as well as its good spectral resolution, provide the opportunity to place stringent constraints on the geometry of the nuclei of ESO103-G35 and NGC7172 and any surrounding gaseous media.

## 2 The sample

### 2.1 NGC 7172

NGC 7172 is an obscured, almost edge on galaxy. It belongs to the compact group HCG90 and has a redshift of z=0.0086. It is classified as a type 2 Seyfert galaxy. Strong infrared emission from the nucleus is observed, which is variable on time scales of approximately 3 months (Sharples et al. 1984). The first X-ray observation which was performed by $`\mathrm{𝐸𝑋𝑂𝑆𝐴𝑇}`$ showed a $`\mathrm{\Gamma }`$=1.84 power-law absorbed by $`N_H10^{23}\mathrm{cm}^2`$ (Turner & Pounds 1989). Later observations with $`\mathrm{𝐺𝑖𝑛𝑔𝑎}`$ confirmed the above results and also indicated that a reflection component may be important (Smith & Done 1996). Furthermore the analysis of ASCA observations (Turner et al. 1997) showed the presence of a 6.4 keV emission line with an equivalent width (EW) of $`68_{35}^{+36}`$ eV. In contrast with the previous results they found a rather harder photon index, $`\mathrm{\Gamma }1.7`$, while the amount of the absorption column remained constant. Guainazzi et al. (1998), using 2 $`\mathrm{𝐴𝑆𝐶𝐴}`$ observations separated by one year, found an even flatter photon index of $`\mathrm{\Gamma }1.5`$ absorbed by the same column density found in previous observations. They further detect short-term variability (within a day) of about 30 percent. Until 1995, before the BeppoSAX observations, the flux of NGC 7172 has remained almost constant between $``$ 3 and $`5\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. Since May 1996 the source has reached a very low flux state of $`1\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. The flux variability history of NGC7172 is summarized in Fig. 1.

### 2.2 ESO 103-G35

The HEAO A2 hard X-ray source 1H 1832-653 was identified as the S0/Sa type galaxy ESO 103-G35, with a redshift of z=0.013, by Phillips et al. (1979). They classified the AGN as a type 1.9 Seyfert galaxy because of the weakness of $`H\beta `$ emission line and the presence of a broad $`H\alpha `$ line. $`\mathrm{𝐸𝑋𝑂𝑆𝐴𝑇}`$ observations, showed a steep photon index of $`\mathrm{\Gamma }=1.92_{0.34}^{+0.34}`$ absorbed by a high column density of N$`{}_{H}{}^{}>10^{23}\mathrm{cm}^2`$, which was found to vary within factor of 2 in a period of 90 days possibly due to motion of material in the Broad Line Region across the line of sight (Warwick, Pounds & Turner 1988). Later observations obtained by $`\mathrm{𝐺𝑖𝑛𝑔𝑎}`$ strongly suggested the presence of a narrow emission line at 6.4 keV and also of a reflection component with a photon index $`\mathrm{\Gamma }=2.19_{0.08}^{+0.08}`$ (Smith & Done 1996). $`\mathrm{𝐴𝑆𝐶𝐴}`$ observations (Turner et al. 1997) confirmed the previous $`\mathrm{𝐸𝑋𝑂𝑆𝐴𝑇}`$ results and showed strong evidence for the presence of a triplet Fe line (6.4, 6.68, 6.98 keV). They also found a decrease by a factor of two in the continuum flux as compared to $`\mathrm{𝐺𝑖𝑛𝑔𝑎}`$. Finally, Forster,Leighly & Kay (1999), analyzing three $`\mathrm{𝐴𝑆𝐶𝐴}`$ observations spanning a period of $``$2 years (1994-1996), found a decrease in the EW of the Fe line while the continuum flux increased by a factor of 2. An edge at 7.4 keV was also found in the 1996 observation. The long term variability history of ESO 103-G35 is summarized in Fig. 2.

## 3 Observations

The scientific instrumentation on board the Italian-Dutch X-ray Satellite BeppoSAX includes a Medium Energy Concentrator Spectrometer, MECS, which consists of three units, (Boella et al 1997) a Low Energy Concentrator Spectrometer, LECS, (Parmar et al 1997) a High Pressure Gas Scintillation Proportional Counter, HPGSPC, (Manzo et al 1997) and a Phoswich Detector System, PDS, (Frontera et al 1997), all of which point in the same direction. The MECS instrument consists of a mirror unit plus a gas scintillation proportional counter and has imaging capabilities. It covers the energy range between 1-10 keV with a spatial resolution of about 1.4 arcmin at 6 keV and a spectral resolution of 8 per cent at 6 keV. The LECS instrument is similar to the MECS and operates in the 0.1-10 keV band with the same spectral and spatial resolution of the LECS. The PDS is a direct view detector with rocking collimators and extends the $`\mathrm{𝐵𝑒𝑝𝑝𝑜𝑆𝐴𝑋}`$ bandpass to high energies (13-300 keV). Its energy resolution is 15 per cent at 60 keV. In our analysis we use MECS and PDS data only. Our objects have not been detected by the HPGSPC detector. We did not use the LECS data as there are nearby contaminating soft X-ray sources, most probably stars in case of ESO 103-G35, while in case of NGC 7172 the contamination is due to the presence of the Hickson compact group H90. These sources are not present in the hard MECS images. Another reason for excluding the LECS data was the poor photon statistics (see tables 1 and 2).

We analyze here two observations of NGC7172 and ESO 103-G35 both obtained with $`\mathrm{𝐵𝑒𝑝𝑝𝑜𝑆𝐴𝑋}`$ within an interval of about one year. More specifically, NGC7172 was first observed in 1996 (from October 15th to 16th) and then in 1997 (from November 6th to 7th). ESO 103-G35 was first observed in 1996 (from October 3rd to 4th) and then in 1997 (from October 14th to 15th). Tables 1 and 2 list the exposures times and net count rates (background subtracted) for all four observations. The three different (or two after May 7 1997, when MECS1 failed) MECS units are merged together in order to increase the signal to noise ratio. This is feasible because the three MECS units show very similar performance and the difference in the position of the optical axis in the three units is smaller than the scale on which the vignetting of the telescopes varies significantly ($`>`$5 arcmin).

## 4 Timing Analysis

In Fig. 3 and 4 we present the light curves for the MECS instrument for each observation separately, using a binning time of 3 ks (the errors correspond to the 68 per cent confidence level). When we perform a constant fit, the resulting $`\chi ^2`$ implies short term variability in both objects. In particular, for the first observation of NGC 7172 we find that the best fit average count rate ($`0.147\pm 0.007`$) gives a $`\chi ^2`$=36 (for 24 degrees of freedom, dof), while for the second observation we find $`0.052\pm 0.005`$ with $`\chi ^2`$/dof=64/33. In the case of ESO 103-G35 the best fit constant is $`0.305\pm 0.004`$ corresponding to $`\chi ^2`$/dof=72/29 for the first observation while for the second one we find $`0.122\pm 0.004`$ with $`\chi ^2`$/dof=10/8. The above $`\chi ^2`$ suggest the presence of statistically significant ($`>3\sigma `$) short-term variability in the 1997 and 1996 observations of NGC7172 and ESO103-G35 respectively. Our results are compatible with previous ASCA findings. More specifically, Guainazzi et al (1998) detected a flux variability of about 30 per cent, on time scales of hours in a 1996 $`\mathrm{𝐴𝑆𝐶𝐴}`$ observation of NGC 7172. Furthermore, Forster et al (1999) using $`\mathrm{𝐴𝑆𝐶𝐴}`$ data found statistically significant short term variability in a 1996 observation of ESO 103-G35.

The variability amplitude was estimated in all four observations by means of the excess variance $`\sigma _{rms}^2`$ (for a definition of this quantity see Nandra et al 1997). In table 3 we present the values of $`\sigma _{rms}^2`$ and the unobscured luminosity for each observation. Unfortunately a comparison with the results found by Nandra et al (1997) for a sample of Seyfert 1 galaxies is not straightforward. This is due to both the different exposure times and binning time we have used. For example our large binning time has the effect to suppress the $`\sigma _{rms}^2`$ value. Indeed our values are slightly lower than the values quoted in Nandra et al (see their Fig 4).

## 5 Spectral Analysis

We used an extraction radius of four arcmin. The above area encircles more than 85 per cent of the photons in the 2-10 keV energy band. The spectrum of the background was estimated from source free regions of the image. In particular, we consider an annulus with inner radius 6 arcmin and outer radius 12, for the MECS instrument. We used the Rise Time threshold background subtracted PDS spectra released from the BeppoSAX archive. We use data between the energy ranges 1.65-10 keV and 15-50 keV for the MECS and PDS detectors respectively where the response matrices are well calibrated. We do not consider the PDS data above 50 keV as in this energy range the signal to noise ratio is low. The MECS spectral files were rebinned linearly to give a minimum of 20 counts per channel while the PDS spectral files were binned logarithmically in 18 channels. The spectral fitting was carried out using the XSPEC v10 software package. All errors quoted in the best-fit spectral parameters correspond to the 90 per cent confidence level. The MECS and PDS data were fitted simultaneously. A relative normalization factor was introduced between the PDS and MECS data spectrum since the BeppoSAX instruments show some mismatches in the absolute flux calibration. We assumed a PDS to MECS normalization factor between 0.77 to 0.95 (Fiore, Guainazzi & Grandi 1999).

### 5.1 NGC 7172 (1996)

We first fit a simple model consisting of a power-law absorbed by a neutral column $`N_H`$ ( Fig. 5). This provides a good fit to the data ($`\mathrm{\Delta }\chi ^2`$=145 for 146 dof). The addition of a broad Gaussian line improves significantly the fit ($`\mathrm{\Delta }\chi ^2`$=10). The best-fit line energy was $`6.0_{0.6}^{+0.6}`$ keV. Since this value is consistent with 6.4 keV, which indicates a neutral or weakly ionized Fe, (e.g. Nandra et al 1997) we fix it at 6.4 keV in the rest of the analysis. The line width is $`\sigma _{K_\alpha }=0.8_{0.4}^{+0.3}`$ keV, so it is marginally resolved given the energy resolution of the MECS. When we fix the line width at 0.1 keV we obtain an increase of $`\mathrm{\Delta }\chi ^2`$=2.8. The best-fit parameters are presented in table 4.

In fig 5 residuals can be seen above 7 keV. We therefore attempted to include an absorption edge at $``$7.2 keV. We found that the addition of the edge significantly improves the fit ($`\mathrm{\Delta }\chi ^2`$=3.4 for the addition of one parameter). The optical depth of the edge is $`0.18_{0.12}^{+0.11}`$. The rather flat spectral index $`\mathrm{\Gamma }=1.64_{0.09}^{+0.12}`$ could suggest that a reflection continuum is present. Therefore, we further tried to add a reflection component to our model. We used the PEXRAV model in XSPEC (Magdziarz & Zdziarski 1995). We assumed that both the reflection component and the power-law are absorbed by the same column density. We fix the PEXRAV normalization to that of the power-law. We also arbitrarily fix the inclination angle (of the plane of the reflecting material) to $`i=60^{}`$, since the shape of the reflection spectrum in our energy band is relatively independent of the inclination angle. The cut-off energy of the incident power-law was set at 300 keV. Finally, the Fe and light element abundances were kept fixed at the solar abundance values. We find that a reflection component is not statistically significant ($`\mathrm{\Delta }\chi ^21.5`$ for one additional parameter). However the rather uncertain value of the reflection parameter (R$`=1.2_{0.9}^{+0.7}`$) and the steeper spectral photon index of $`\mathrm{\Gamma }=1.78_{0.19}^{+0.26}`$ ,which is closer to the canonical AGN spectral index (Nandra & Pounds 1994) suggest that a reflection component cannot be ruled out. The fraction of the intrinsic reflected flux relative to the total in the 2-50 keV band is 0.25.

### 5.2 NGC 7172 (1997)

Similar to the 1996 observation of NGC 7172, an absorbed power-law model ($`\mathrm{\Gamma }1.5`$) provides a good fit ($`\mathrm{\Delta }\chi ^2`$=112 for 109 dof, see Fig. 5). The addition of a Gaussian line improves the fit significantly ($`\mathrm{\Delta }\chi ^2=6)`$. The line energy is $`6.74_{0.42}^{+0.26}`$ keV consistent with neutral Fe and therefore hereafter we fix the line energy at 6.4 keV. The range of the line width is 0.08-1 keV. The best fit parameters are shown in table 4. Inspection of the residuals shows that there is no evidence for the presence of an absorption edge. However, we attempt to include an edge at 7.25 keV based on the previous 1996 results. We found that there is no improvement in the fit ($`\mathrm{\Delta }\chi ^2<1`$). The optical depth of the edge is $`0.08_{0.08}^{+0.15}`$. The rather flat spectral index $`\mathrm{\Gamma }=1.53_{0.16}^{+0.08}`$ could again suggest the presence of a reflection component. However, when we try to include such a component, we find that, although it does not improve the fit ($`\mathrm{\Delta }\chi ^2=1.5`$ for one additional parameter), the photon index becomes steeper $`\mathrm{\Gamma }=1.81_{0.40}^{+0.60}`$ and the reflection parameter is R=$`2.5_{1.8}^{+1.3}`$. The unabsorbed reflected flux is 40 percent of the total flux in the 2-50 keV energy band. Comparison between the spectral fits in the 1996 and 1997 observations reveals that there is no change in the photon index and the column density within the errors.

Finally, we perform joint fits in order to investigate any possible variability in the normalization of the different spectral components. We fit both NGC 7172 observations with an absorbed power law model ($`\mathrm{\Gamma }=1.62_{0.09}^{+0.07}`$ and N$`{}_{H}{}^{}=9.1_{0.6}^{+0.3}\times 10^{22}\mathrm{cm}^2`$) plus a Gaussian component at 6.4 keV to account for the Fe line. The width of the line is $`0.6_{0.5}^{+0.4}`$ keV and therefore remains unconstrained. Comparison of the two observations reveals weak evidence for spectral variability. The power-law normalization has decreased by about a factor of two between the 1996 and 1997 observations. This may be followed by a similar change in the Fe line normalization. This is evident in Fig. 6 where we present the 68, 90 and 99 per cent confidence contours for the Fe line against the power-law normalization. The contours were extracted using the above best fit model. The EW are $`280_{130}^{+140}`$ and $`233_{160}^{+190}`$ eV for the 1996 and 1997 observations respectively. However we note that when we fix the width of the line at 0.1 keV (this increases $`\chi ^2`$ by 2.5) we do not detect any significant variability in the normalization of the line.

### 5.3 ESO 103-G35 (1996)

A poor fit ($`\chi ^2`$=240 for 173 dof) is obtained with a single power-law model, emphasizing the importance of an Fe emission line. This is again evident in Fig. 7 where we plot the absorbed power-law model and residuals. Indeed when we include a Gaussian line the $`\chi ^2`$ is reduced by 66. The Fe line has an energy of $`6.37_{0.10}^{+0.09}`$ keV consistent with neutral Fe and therefore hereafter we fix it at 6.4 keV. Its width is $`\sigma =0.3_{0.1}^{+0.1}`$ keV and therefore the line remains essentially unconstrained. The photon index, $`\mathrm{\Gamma }=1.87_{0.09}^{+0.06}`$, is consistent with the typical value of 1.9 for Seyfert galaxies (Nandra & Pounds 1994). The best-fit spectral parameters are presented in table 5. We further investigate the existence of excess Fe absorption near 7.5 keV following the claims of Warwick et al. (1993). Although there is some hint for an Fe edge in Fig. 7, we find no significant improvement when we include an edge component ($`\mathrm{\Delta }\chi ^2=1.5`$ for two additional parameters). The edge is found at $`7.35_{0.35}^{+0.65}`$ keV with optical depth $`0.1_{0.10}^{+0.08}`$. Finally, based on the $`\mathrm{𝐴𝑆𝐶𝐴}`$ results (Turner et al 1997) we tried to include two emission lines at 6.68 and 6.96 keV (ionized Fe) but the reduction in $`\chi ^2`$ was less than 1. A reflection component is not statistically significant as the reduction of $`\chi ^2`$ is very low ($`\mathrm{\Delta }\chi ^2`$= 1.1 for one additional parameter). Note that the power law spectral index in the case of a reflection component is $`\mathrm{\Gamma }=2.0_{0.2}^{+0.2}`$ while R=$`1.2_{0.6}^{+0.4}`$. The intrinsic flux of the reflection component corresponds to the 20 per cent of the total flux in the 2-50 keV energy band.

### 5.4 ESO 103-G35 (1997)

An acceptable fit is obtained when we use the absorbed power-law model with $`\chi ^2`$=70 for 71 dof, (see Fig. 7). An FeK emission line is highly significant yielding $`\mathrm{\Delta }\chi ^2`$=14. The Fe line energy is $`6.5_{0.15}^{+0.10}`$ keV and therefore again we fix it at 6.4 keV. The line is unresolved: its width range is 0-0.4 keV. The photon index, $`\mathrm{\Gamma }=1.81_{0.30}^{+0.09}`$, is consistent with the previous 1996 observation. The results are presented in table 5. A reflection component cannot be ruled out as the reflection parameter is R=$`1.6_{1.2}^{+2.0}`$. However,the reduction in $`\chi ^2`$ is very low ($`\mathrm{\Delta }\chi ^2`$= 1.7 for one additional parameter). When the reflection component is included the power law spectral index steepens to $`\mathrm{\Gamma }=1.97_{0.30}^{+0.30}`$. The ratio between the unabsorded reflected flux and the total flux in the 2-50 keV energy band is 0.3. Next, we investigate the existence of an Fe edge at 7.35 keV. However, the reduction in $`\chi ^2`$ was less than 1 for one additional parameter. The optical depth of the edge is $`0.11_{0.11}^{+0.10}`$. In general, we detect no spectral variability between the 1996 and 1997 observations in the sense that both the absorption column and the photon index remain constant.

We perform a joint fit in the 1996 and 1997 observations in order to detect any possible line flux variations. We used a model consisting of an absorbed power-law plus a Gaussian line at 6.4 keV. The best-fit photon index is $`1.83_{0.09}^{+0.09}`$ and the column density is $`18.5_{0.8}^{+0.2}\times 10^{22}\mathrm{cm}^2`$. The line is unconstrained, ie the range of the line width is 0.15-0.36 keV. The EW are $`220_{50}^{+120}`$ and $`330_{140}^{+170}`$ eV for the 1996 and 1997 observations respectively. The continuum flux is reduced approximately by a factor of two. However, there is no evidence that this variation is followed by the line flux, in contrast to the case of NGC7172 (see Fig. 8 where we plot the 68, 90 and 99 per cent joint contours for the power-law and Fe line normalization).

## 6 discussion

### 6.1 NGC 7172

The spectral analysis of both observations implies that NGC 7172 is absorbed by a column density of N$`{}_{H}{}^{}=9.1_{0.6}^{+0.3}\times 10^{22}\mathrm{cm}^2`$, in good agreement with previous $`\mathrm{𝐸𝑋𝑂𝑆𝐴𝑇}`$ (Turner & Pounds 1989), $`\mathrm{𝐺𝑖𝑛𝑔𝑎}`$ (Smith & Done 1996) and $`\mathrm{𝐴𝑆𝐶𝐴}`$ (Turner et al 1997) observations. The data require the presence of an Fe line at 6.4 keV. In the 1996 BeppoSAX observation the line is marginally resolved, at less than the 2$`\sigma `$ confidence level, while in the 1997 observation (where the photon statistics are poorer) the line width is unconstrained. We note that the line has been resolved in previous ASCA observations (Guainazzi et al. 1998). There is some evidence (Fig. 6) for a variation in the normalization of the line. This is in agreement with Guainnazi at al (1998) who presented similar evidence for a variation of the line based on two $`\mathrm{𝐴𝑆𝐶𝐴}`$ observations (1995 and 1996). The line variation in conjunction with the line width has important implications for the geometry of the circumnuclear matter in NGC7172. Indeed, the Fe line must originate close to the nucleus possibly in the outer regions of an accretion disk. This would simultaneously explain the fact that the continuum flux variations are tracked by changes in the line flux, and the large width of the line. Indeed, the line broadening may be due to the large velocities in the disk similar to that witnessed in MCG-6-30-15 by Tanaka et al. (1995). The observed broadening in the case of NGC 7172 implies a velocity of $``$70000 $`\mathrm{km}\mathrm{s}^1`$. We tried to investigate further the above hypothesis including a line from a relativistic accretion disk (for an accretion disk viewed at an intermediate inclination angle with solar abundances). The fit was only marginally improved ($`\mathrm{\Delta }\chi ^21`$). We note that Georgantopoulos & Papadakis (2000) find no significant variation of the normalization of the line in NGC 7172 using several RXTE observations spanning a period of around a week. This would set a lower limit of about a light week to the size of the line emitting region.

The observed EW cannot be interpreted in terms of transmission processes in the torus alone. Indeed Leahy & Creighton (1993) using Monte Carlo simulations found that a column density of N$`{}_{H}{}^{}9\times 10^{22}\mathrm{cm}^2`$ could produce an Fe line EW of only $`60`$ eV. This independently suggests that a reflection component may be present. Indeed, in this case the EW of the line should be increased by $``$150 eV due to reflection on the accretion disk (Reynolds, Fabian & Inoue 1995, George & Fabian 1991). Alternatively the observed EW may be explained on the basis of a high Fe abundance. Indeed, excess absorption from Fe is present in the first observation of NGC 7172 at $`7.2`$ keV, confirming earlier reports by Warwick et al. (1993); however, note that superior spectral resolution observations with the $`\mathrm{𝐴𝑆𝐶𝐴}`$ SIS (Turner et al 1997) do not corroborate this result. The edge energy is consistent with a neutral or low ionization material surrounding the source (Kallman & McCray 1982). The optical depth of the edge corresponds to an absorption column of $`22.5_{14.2}^{+15.3}\times 10^{22}\mathrm{cm}^2`$. This large column density may suggest an Fe abundance well above the solar value. If this is the case, the observed EW could be explained by transmission processes alone. Finally note that cases where a strong Fe line is detected although there is no strong evidence for a reflection component are not uncommon (eg Gilli et al 2000, Reynolds et al 1995).

The power-law photon index has a rather flat value of $`1.6`$ while it remains constant in the two BeppoSAX observations. The value of the photon index is in agreement with previous ASCA observations (Guainazzi et al 1998, Ryde et al 1997) who found a photon index of $`1.47\pm 0.15`$. However, it is clearly flatter than that found in $`\mathrm{𝐸𝑋𝑂𝑆𝐴𝑇}`$ (Turner et al 1989) and $`\mathrm{𝐺𝑖𝑛𝑔𝑎}`$ (Smith & Done 1996) ($`\mathrm{\Gamma }1.8`$). Therefore the above may be suggesting long term spectral variability (Ryde et al. 1997). For example, assuming that the X-ray emission comes from Comptonization of UV photons from an accretion disk on a hot electron corona (Rybicki & Lightman 1979), a steepening of the spectral index may be attributed to a change of the temperature or the optical depth of the corona. Of course it is difficult to make a detailed comparison due to differences in the energy bandpass of the various instruments as well as in the spectral analysis: for example the addition of a reflection component would result in the steepening of the spectral index. However, we note that there is no correlation between the flux level and the photon index observed in the case of the NGC7172 observations in the past 20 years. Finally we note that despite the flat photon index found above, our data do not significantly support the existence of any reflection continuum. A recent analysis of $`\mathrm{𝑅𝑋𝑇𝐸}`$ monitoring observations of NGC 7172 (Georgantopoulos & Papadakis 2000) again shows no evidence for reflection component.

### 6.2 ESO 103-G35

The spectral analysis of the two observations of ESO103-G35 implies that a simple power-law model absorbed by a column density of N$`{}_{H}{}^{}2\times 10^{23}\mathrm{cm}^2`$ plus a Gaussian line at 6.4 keV provides a good fit to our data. EXOSAT observations of ESO103-G35 (Warwick et al. 1988) reveal a 50 per cent decrease in the column density over a 90 day period. This variation was interpreted in terms of an X-ray absorbing screen composed of broad-line clouds moving around the X-ray source. Alternatively, the N<sub>H</sub> variation could be attributed to photoionization of the neutral column as the flux increases; in this case an ionized Fe edge should be detected. However, in our analysis we detect no variation of the column within the errors. An Fe line at 6.4 keV is strongly needed to fit our data. Unfortunately, the width of the line remains unconstrained. However, ASCA  SIS observations with better spectral resolution have probably resolved the line into three components at 6.4, 6.68 and 6.96 keV (Turner et al 1997). Monte Carlo simulations (Leahy & Creighton 1993) show that a column density of N$`{}_{H}{}^{}2\times 10^{23}\mathrm{cm}^2`$ could produce an Fe line EW of $``$170 eV. However, in the case of ESO 103-G35 where the line originates possibly far away from the central region, the EW varies accordingly to the variations of the continuum. Although a comparison is not straightforward, it is possible that a large fraction of the line emission could originate from transmission processes.

Some residuals around 7.5 keV (see Fig. 7) may indicate the presence of an absorption edge due to ionized Fe. This was first suggested by $`\mathrm{𝐺𝑖𝑛𝑔𝑎}`$ results (Warwick et al. 1993) and was confirmed later with $`\mathrm{𝐴𝑆𝐶𝐴}`$ observations (Forster et al 1999). The power-law photon index has a steep value of $`1.85_{0.05}^{+0.05}`$ which is consistent with the typical value 1.9 for Seyfert galaxies (Nandra & Pounds 1994). Despite our large energy bandpass we do not find strong evidence for a reflection continuum in contrast to previous $`\mathrm{𝐺𝑖𝑛𝑔𝑎}`$ (Smith & Done 1996) and $`\mathrm{𝑅𝑋𝑇𝐸}`$ (Georgantopoulos & Papadakis 2000) observations who found evidence for such a component but in a smaller band (3-20 keV). Higher signal to noise observations are necessary in order to test the presence of such a component.

### 6.3 SPECTRAL VARIABILITY IN SEYFERT 2 GALAXIES

Previous studies of spectral variability in Seyfert-2 galaxies have shown some interesting results. In particular, Iwasawa et al (1994) and Griffiths et al. (1998) analyzed non-simultaneous Ginga, ASCA, ROSAT and BBXRT observations to investigate the spectrum of the Seyfert-2 galaxy Mrk 3 . Their analysis demonstrated that the Fe line does show some variability on time-scales of a few years. Systematic monitoring observation became feasible with the RXTE  mission. In particular, Georgantopoulos et al. (1999), Georgantopoulos & Papadakis (2000), Smith, Georgantopoulos & Warwick (2000) present monitoring observations of several Seyfert-2 galaxies (Mrk3, ESO 103-G35,IC 5063, NGC 4507, NGC 7172 & Mrk 348) spanning time periods from about seven days to seven months. In most cases, they detect spectral variability in the sense that the column density decreases with increasing flux. However, they find no evidence for Fe line variability. Here instead, we find tentative evidence for line variability in NGC 7172 in a period of one year confirming previous ASCA results (Guainazzi et al 1998). To our knowledge this is one of the very few examples of Fe line variability in Seyfert-2 galaxies. Our finding requires at least one of the line emitting regions has dimension not much more than about a light year. It would seem more likely that this scale size corresponds to an accretion disk or inner cloud structure rather than the inner extent of the putative molecular torus. Indeed, the molecular torus is probably a much greater structure. At least in the case of NGC 1068 observations of molecular hydrogen (Tacconi et al 1994) showed that the obscuring torus is located between 30 and 300 pc.

## 7 conclusions

We have analyzed four BeppoSAX observations of NGC7172 and ESO103-G35 (two observations for each galaxy separated by about a year). The goal was to search for spectral variability in a large energy band and therefore to attempt to constrain the geometry of the circumnuclear matter in these two galaxies. The spectra of both galaxies are fitted by an absorbed power-law plus a neutral Fe line in agreement with previous Ginga and ASCA observations. A neutral Fe edge is probably detected in the 1996 observation of NGC7172. Our data do not require a reflection component at a statistically significant level. We detected flux variability in both objects by about a factor of about two. We find no significant evidence for spectral variability in the sense that the power-law spectral index and the column density have remained constant in both galaxies implying a homogeneous obscuring screen and constant physical conditions in the accretion disk corona despite the continuum flux variations (but see Warwick et al. 1988, Georgantopoulos & Papadakis 2000). However, in the case of NGC7172 we detect a variation of the Fe line flux by a factor of two following the power-law flux, confirming previous ASCA results. This provides strong support to a scenario where the Fe line at least in the case of NGC7172 originates in a region close to the nucleus. These findings support the idea that one region responsible for the line emission should lie within one light year away from the nucleus. This region may be associated with the outer part of an accretion disk.

The above results clearly emphasize the strength of monitoring observations in probing the geometry of the central engine in AGN. High resolution spectroscopic observations with gratings with both Chandra and XMM are expected to revolutionize our knowledge of the physical conditions of the circumnuclear matter in AGN. These combined with future monitoring observations with large effective area instruments such as EPIC on-board XMM are expected to provide for the first time a detailed mapping of the central regions in AGN.

## 8 acknowledgements

We thank the anonymous referee for many useful comments and suggestions. This project was funded by a Greek-Italian scientific collaboration under the title “Observations of active galaxies with the Italian astrophysics mission BeppoSAX”. AA is grateful to the Bologna observatory group for their warm hospitality. We are grateful to M. Cappi, C. Vignali and R. Lazio for many useful discussions.


text:
# Critical Exponents of the Two-layer Ising Model

## Introduction

The two-layer Ising model (TLIM), as a simple generalisation of the two-dimensional Ising model and a simple model for the magnetic ultra-thin films, has been studied for a long time Ball64 ; Alla70 ; Abe70 ; Suzu71 ; Bind74 ; Oitm75 . Cobalt films grown on a Cu(100) crystal, for instance, have highly anisotropic magnetisation Oepe90  and could therefore be viewed as layered Ising-spin systems. It has been found that capping PtCo in TbFeCo to form a two-layer structure has applicable features, for instance, raising the Curie temperature and reducing the switching fields for over-writable magneto-optical disks Shim92 . Apart from various possible applications to real physical materials, the model is theoretically interesting for its rich phase structure. The model has several interesting equivalents, such as a two-species gauge invariant Ising model Li94 , a spin-$`\frac{3}{2}`$ Ising model Hori96 , a model of the dilute lattice gas and a quantum spin-$`\frac{1}{2}`$ ladder at zero temperature. The TLIM is important for the investigation of crossover from the 2-dimensional Ising model to the 3-dimensional one. In particular, it has been argued that the critical point of the latter could be found from the spectrum of the TLIM Wosi94 .

In recent years, some approximation methods have been applied to this model Hori96 ; Ange92 ; Ange95 ; Ange97 ; Lipo93 ; Lipo98 ; Hori97 ; Hu99 . A critical line has been found in all these studies. As expected, the Curie temperature is very sensitive to the inter-layer interaction. Many discussions have been directed to the shift exponent at the decoupling point. Abe Abe70  and Suzuki Suzu71  have predicted $`\gamma =7/4`$ for the isotropic model many years ago. Only in recent years the computational results have converged to the prediction with still about $`2.3\%`$ deviation Hori97 ; Lipo98 .

Apart from the shift exponent, it is also interesting to study the critical behaviour along the critical line. The model has the same critical exponents at the two ends of the critical line corresponding to the solvable decoupling limit and the infinite interlayer coupling limit. But it is clear that the decoupled system has a higher symmetry than the coupling layers, hence one cannot assume a priori that the TLIM belongs to the one-layer Ising universality class. Unusual finite-size effects have been observed in the dynamic process Li00 . It has been proposed that the critical exponents (or some of them) would vary continuously along the critical line Li94 ; Hori97 . However, there are also arguments in favour of unchanged exponents. Based on their correlation length equality method, Angelini et al. argued that the exponent $`\nu `$ is a constant Ange95 . The accurate computation of critical exponents along the critical line remains a difficult unsolved problem.

It is our purpose to provide a reliable prediction for the critical exponents based on the corner transfer matrix renormalisation group (CTMRG) method Nish96a ; Nish96b ; Nish97 .

The CTMRG method follows White’s idea of density matrix renormalisation group Whit92 ; Whit93 . It reduces the phase-space dimensions by omitting the eigenstates corresponding to small eigenvalues of the density matrix. But instead of constructing the density matrix through a series of column transfer matrices, Nishino and Okunishi made use of the corner transfer matrices. The CTMRG method has shown many merits for the two-dimensional classical models which can be mapped into interaction-around-face (IRF) models. Since the system can be easily extended to very large lattices, the finite size effect can be ignored. The error mainly comes from the finite dimensions of the renormalised phase space, which only becomes severe for huge lattices at or very close to the critical point. Numerically it has been observed that the extrapolation from the finite size scaling can provide very accurate information about the criticality.

We will focus on the TLIM composed of two identical layers. In this case, the model has a symmetry of interchanging the spins of the two layers. This symmetry is vital to the critical behaviour of the model. A breakdown of universality due to a small symmetry breaking Hori97  has been observed. In order to keep this symmetry, special care should be taken in the programming. Technically, the symmetry reduces the dimensions of the transfer matrix remarkably. But the extra degeneracy will slow down the convergence of the iteration which makes the problem much more difficult compared to the one-layer Ising model.

The critical behaviour can be easily observed in various quantities, such as the rapid increase of the magnetisation, the lambda peak of the capacity, and the anomalous slowing down of convergence of the renormalisation group iterations. However, to calculate the critical exponents, one needs to locate the critical points to high precision. This is done by a careful extrapolation of the renormalisation dimension $`m`$.

In Section 2, the TLIM is defined and transformed into an IRF model which is suitable for the CTMRG study. In the same section, the CTMRG method particular for the considered model will be introduced. The finite-size scaling and the finite renormalisation dimension effect is studied in Section 3. At the same time the critical points are located, critical exponents of the scaling law are computed. Results for the critical exponents $`\eta `$ and $`\nu `$ can be found also in Section 3. Discussions and conclusions are given in the last section.

## The CTMRG method for the TLIM

The Hamiltonian of the model is defined by

$$H=J_1\underset{<i,j>}{}s_is_jJ_2\underset{<i,j>}{}\sigma _i\sigma _j\lambda \underset{i}{}s_i\sigma _i,$$
(1)
where $`s_i=\pm 1`$ and $`\sigma _i=\pm 1`$ are Ising spins on two identical square lattices. For concreteness, let us say, $`s_i`$ is at the site $`i`$ of layer 1, $`\sigma _i`$ is at the corresponding site of layer 2 which is also labelled by $`i`$. Then $`<i,j>`$ should be understood as a nearest neighbour pair of sites either on layer 1 or layer 2, depending on which of the spin variable, $`s`$ or $`\sigma `$, is concerned. In this paper, we only discuss the symmetric ferromagnetic case with nearest neighbour coupling $`J_1=J_2=J>0`$ and with interlayer coupling $`\lambda =\rho J>0`$.

To obtain the equivalent IRF model, two species of Ising spins, $`u_{ij}=\pm 1`$ and $`v_{ij}=\pm 1`$, are introduced for each link $`<i,j>`$. The Hamiltonian (1) is generalised to

$$H^{}=J^{}\left[\underset{<i,j>}{}u_{ij}(s_i+s_j)+\underset{<i,j>}{}v_{ij}(\sigma _i+\sigma _j)+\rho ^{}\underset{i}{}s_i\sigma _i\right],$$
(2)
where $`J^{}`$ and $`\rho ^{}`$ are defined as

$`J^{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}(e^{2J}+\sqrt{e^{4J}1}),`$ (3)
$`\rho ^{}`$ $`=`$ $`{\displaystyle \frac{\rho J}{J^{}}}.`$ (4)
It is not difficult to show that, by summing over all $`u`$ and $`v`$spins, the TLIM is recovered apart from a overall irrelevant constant factor to the partition function. Inversely, if the s-spins and $`\sigma `$-spins are summed, one will obtain the IRF model with four u-spin interaction and four v-spin interaction around each vertex. The vertex weight of the IRF model is given by

$$W_{abcd}=\underset{s,\sigma =\pm 1}{}\mathrm{exp}\left\{J^{}\left[s(u_a+u_b+u_c+u_d)+\sigma (v_a+v_b+v_c+v_d)+\rho ^{}s\sigma \right]\right\},$$
(5)
where $`a,b,c`$, and $`d`$ denote the four (u,v)-spin pairs in the links joining to the vertex. Special weights in the boundary, which have links less than four, can be defined in a similar way as (5) but the values of $`s`$ and $`\sigma `$ are subject to the boundary condition. For instance, for the fixed boundary condition, the weight along the boundary, $`W_{abc}`$, is defined by

$$W_{abc}=\mathrm{exp}\left[J^{}\left(u_a+u_b+u_c+v_a+v_b+v_c+\rho ^{}\right)\right].$$
(6)
The vertex at the corner has a weight

$$W_{ab}=\mathrm{exp}\left[J^{}\left(u_a+u_b+v_a+v_b+\rho ^{}\right)\right].$$
(7)
They are schematically given in Figure 1a,b,c, respectively.

The partition function is a trace of products of all vertex weights. The trace here means to sum over all spins. In the following, whenever two vertices are connected by a link, a summation over the (u,v)-spin pair in that link is implied. The vertex weights of a quarter lattice, after summation over all inner spins and boundary spins as it has been implied, is called a corner transfer matrix(CTM), denoted by $`C_{\alpha \beta }(N)`$. Here $`\alpha `$ and $`\beta `$ denote the spin configurations of two open sides of the corner respectively, and $`N`$ is the size of the corner. A schematic representation for a CTM is given in Figure 2a.

Another matrix we need is the so-called half-row transfer matrix (HRTM) denoted by $`P_{\alpha \beta ,a}(N)`$. It is defined by the recursion relation

$`P_{\alpha \beta ,a}(N)`$ $`=`$ $`{\displaystyle \underset{c}{}}W_{abcd}P_{\alpha ^{}\beta ^{},c}(N1),`$ (8)
$`P_{bc,a}(1)`$ $`=`$ $`W_{abc},`$ (9)
where $`\alpha ^{}`$ and $`\beta ^{}`$ are spin configurations of two sides of the HRTM with size of $`(N1)`$ ; $`c`$ is the (u,v)-spin pair of the top of that HRTM; $`\alpha `$ and $`\beta `$ are spin configurations of the enlarged HRTM with $`\alpha `$ representing $`(\alpha ^{},b)`$ and $`\beta `$ representing $`(\beta ^{},d)`$ respectively. For each value of $`a`$, $`P_{\alpha \beta ,a}(N)`$ can also be considered as a matrix, denoted by $`P_a(N)`$. Figure 2b is an example of the HRTM.

In the process of the CTMRG iterations, the size of CTM is enlarged by one lattice spacing each time. This is done in this way: (1) glue a HRTM to each open side of the CTM, (2) complement the bigger corner by adding a vertex weight. The recursion relation is

$$C_{\alpha \beta }(N)=\underset{\alpha ^{\prime \prime }\beta ^{\prime \prime }}{}\underset{cd}{}W_{abcd}P_{\alpha ^{}\beta ^{\prime \prime },c}(N1).C_{\alpha ^{\prime \prime }\beta ^{\prime \prime }}(N1)P_{\alpha ^{\prime \prime }\beta ^{},d}(N1).$$
(10)
It is schematically shown in Figure 3.

By use of CTM, the partition function of an even size lattice is simply the trace of the CTM to power four

$$Z_{2N}=Tr(C(N)^4).$$
(11)
Therefore, the CTM and the density matrix have common eigenvectors. Denote the eigenvalues of $`C`$ by $`\{\omega _k|k=0,1,2,\mathrm{}\}`$, assuming that $`\omega _0\omega _1\omega _2\mathrm{}`$, the eigenvalues of the density matrix are $`\{\omega _k^4|k=0,1,2,\mathrm{}\}`$. According the spirit of DMRG, only a certain number of the biggest eigenstates will be kept. For instance, the first $`m`$ eigenstates are kept, the eigenstates with smaller eigenvalues $`\{\omega _k|k>m\}`$ are ignored. $`m`$ is the so-called renormalisation dimension.

The CTMRG iterations contains the following steps: (1) construct $`C(N)`$ of a small CTM which is small enough to be exactly diagonalised; (2) construct $`P_a(N)`$; (3) diagonalise $`C(N)`$; (4) if the dimension of $`C(N)`$ is smaller than $`m`$, all eigenstates are used as basis of the renormalised phase space, otherwise only $`m`$ eigenstates with the biggest eigenvalues are used as basis; (5) by use of the eigenvectors which have been chosen for the renormalised space, form a projection operator $`U(N)`$ that projects old spin states into the renormalised space; (6) use $`U(N)`$ to project $`P_a(N)`$ to the renormalised HRTM $`P_a^r(N)`$ ; (7) form a diagonal matrix $`C^r(N)`$ by the first $`m`$ largest eigenvalues, which is the renormalised CTM; (8) use $`P_a^r(N)`$, $`C^r(N)`$, and the vertex weight $`W`$ to form a bigger CTM, $`C(N+1),`$ and a bigger HRTM, $`P(N+1)`$, and repeat the above procedure from step (3) with $`N`$ replaced by $`N+1`$.

Since the renormalised phase space has dimension $`m`$, the CTMRG in principle can infinitely go on. Therefore the lattice can be easily enlarged to very big sizes.

Particularly for the TLIM, it is efficient to choose the basis to be symmetric and antisymmetric in the exchange of two layers. In this way, the CTM is automatically block diagonal. The symmetric block and the antisymmetric block can be diagonalised separately. This is not just allowing larger $`m`$ and high precision, but is also essential for keeping the presumed symmetry between two layers.

For small interlayer coupling, there are many approximate degenerate eigenstates due to the layer symmetry. If one does not separately consider the symmetric and antisymmetric phase space, it is easy to destroy the layer symmetry in the renormalisation procedure.

Following the method of Nishino and Okunishi Nish96a , the magnetisation $`M`$ can be calculated on an odd size lattice which is constructed by inserting a HRTM matrix between two adjacent CTM matrices and adding a proper vertex at the lattice centre. The probability that the central site has a spin state $`(s,\sigma )`$ is given by

$$M^{s\sigma }=\frac{_{abcd}X_{abcd}^{s\sigma }Tr\left(P_a^rC^rP_b^rC^rP_c^rC^rP_d^rC^r\right)}{_{abcd}W_{abcd}Tr\left(P_a^rC^rP_b^rC^rP_c^rC^rP_d^rC^r\right)},$$
(12)
where $`X_{abcd}^{s\sigma }`$ is defined by

$$X_{abcd}^{s\sigma }=\mathrm{exp}\left\{J^{}\left[s(u_a+u_b+u_c+u_d)+\sigma (v_a+v_b+v_c+v_d)+\rho ^{}s\sigma \right]\right\}.$$
(13)
The magnetisation of layer 1, $`M=<s>`$, is

$$M=M^{++}+M^+M^{}M^+.$$
(14)
Due to the layer symmetry, it is also one half of the total magnetisation. The local energy density $`E`$ can be calculated in a similar way.

$$E=\frac{_{abcd}Y_{abcd}Tr\left(P_a^rC^rP_b^rC^rP_c^rC^rP_d^rC^r\right)}{_{abcd}W_{abcd}Tr\left(P_a^rC^rP_b^rC^rP_c^rC^rP_d^rC^r\right)},$$
(15)
where

$$Y_{abcd}=\underset{s,\sigma }{}X_{abcd}^{s\sigma }\mathrm{ln}(X_{abcd}^{s\sigma }).$$
(16)
## Extrapolation and finite-size scaling

At the critical point $`J_c`$, it is expected that $`M`$ obeys the finite-size scaling form $`ML^{(d2+\eta )/2}`$, where $`L`$ is the lattice size, so

$$\mathrm{ln}M(L)=const\frac{d2+\eta }{2}\mathrm{ln}L.$$
(17)
If we calculate $`\mathrm{ln}M/\mathrm{ln}L`$ and plot it with $`\mathrm{ln}L`$ as the horizontal axis, the curve should be a horizontal line for the critical point $`J_c`$ with the vertical value of the line as $`(d2+\eta )/2`$. For $`J`$ near $`J_c`$, there should be deviations to the power-law behaviour. A simple physics consideration reveals that

the curves should be in the form as shown in Figure 4. So by searching for the “best” horizontal line, one can determine the critical point $`J_c`$ as well as the critical exponent $`\eta `$.

Things are not so easy, however. Figure 5

shows the curves of $`\mathrm{ln}M/\mathrm{ln}L`$ against $`\mathrm{ln}L`$ for $`\rho =0.1,J=0.397726`$ with different renormalisation dimensions $`m`$. This $`J`$ value is quite near to $`J_c`$. For small $`L`$, there exists a region where the lattice size is too small and finite-size scaling has not come into play. We will refer to this region as the small-size region. As $`L`$ increases, the number of states of the system increases while the renormalisation dimension $`m`$ remains unchanged, and we will suffer from the “finite-$`m`$ effect” in due time. To get a more complete knowledge of the finite-$`m`$ effect, let us look at Figure 6

where curves for different $`J`$ with $`m=60`$ are shown. $`m=60`$ is relatively small in our calculation and it has a more drastic finite-$`m`$ effect. When $`L`$ is small, all the three curves shown in Figure 6 increase, which is an indication of the small-size effect. After a certain value of $`L`$, one of the curves begins to bend down and the other two become much flatter. Here the critical behaviour dominates. However, the finite-$`m`$ effect soon shows itself by raising the curves, which is also demonstrated clearly in Figure 5. As a result, only approximately the part inside the dotted-line box resembles Figure 4. It is hard to tell exactly where the finite-$`m`$ effect begin to dominate without comparison with curves of higher $`m`$.

Nevertheless, some useful information can be drawn from the analysis. If a curve bends down for a period of $`L`$, one can judge that it corresponds to a $`J`$ that is smaller than $`J_c`$. In other words, if a curve has a local maximum, it corresponds to a $`J`$ with $`J<J_c`$. Thus the point now is to find $`J_c^{(60)}`$, above which the corresponding curve has no local maximum and below which it has. $`J_c^{(60)}`$ can be understood as a lower limit for $`J_c`$. In general, one can obtain $`J_c^{(m)}`$ as an estimate of $`J_c`$. The former will approach the latter from below as $`m\mathrm{}`$. Meanwhile the vertical value of the horizontal part of the curve with $`J=J_c^{(m)}`$ gives the estimate for $`(d2+\eta )/2`$. To make the above discussion simpler, we have limited it only to the case with similar behaviour as that of $`\rho =0.1`$. Later we will see in Figure 7,

that the small-size effect may show itself differently as $`\rho `$ changes. But similar analyses can still be developed for other cases.

In most of our calculations, we chose $`m`$ to be 60, 80, 120, and 240. For each $`m`$, we calculated the magnetisation with a set of $`J`$ values and then determined the critical point using the method described above. The results for $`\rho =0.1`$ are shown in Table 1.

The last row shows the extrapolated values for $`m=\mathrm{}`$.

An alternate and finer way for dealing with the finite-$`m`$ problem is to extrapolate the data to $`m\mathrm{}`$ first and then analyse the extrapolated data. The extrapolation should lead to reasonable and smooth curves for $`m=\mathrm{}`$, and we assume the data converge monotonically as $`m\mathrm{}`$. This defines our criterion for choosing the extrapolation formula. We still take as an example the case with $`J=0.397726`$, $`\rho =0.1`$, whose curves are plotted in Figure 5. For each value of lattice size $`L`$, we fit the points $`(1/m,M^{(m)}\left(L\right))`$, where $`m=60,80,120,240`$, with

$$y=a_5x^5+a_4x^4+a_3x^3+a_0.$$
(18)
$`a_0`$ then is considered to be the value for $`m=\mathrm{}`$ (see Figure 8).

Eq.(18) does not include terms of $`x^1`$ and $`x^2`$ because including them violates the criteria above. The extrapolated curve is also shown in Figure 5 as a dashed line.

Now that we obtained the curves for $`m=\mathrm{}`$ with different $`J`$, we can search for the $`J`$ value whose $`\mathrm{ln}M/\mathrm{ln}L\mathrm{ln}L`$ curve is closest to a horizontal line. This $`J`$ value is then the estimate of $`J_c`$. In our calculation, we use interpolation technique to help locate the “best” $`J`$, as shown in Figure 9.

The solid lines are those extrapolated to infinite $`m`$. We can obtain the magnetisation for any value of $`J`$ that is between $`J=0.397726`$, for which the curve bends down, and $`J=0.397727`$, where the line bends up, by linear interpolation. We found $`J_c`$ to be $`0.3977262`$ with $`(d2+\eta )/2=0.12507`$ .

The small-size effect for different ratios is presented in Figure 7. For small $`\rho `$, this effect is big. In the case of $`\rho =0.1`$, for instance, the small-size effect is visible with $`L`$ up to 800. However, it seems that we overcome the small-size effect with $`L=2000`$.

Having $`J_c`$ in hand, we can determine another critical exponent $`\nu `$ through the finite-size scaling form of the energy density

$$E(L)E_{\mathrm{}}L^{(d\frac{1}{\nu })}.$$
(19)
The critical point $`J_c`$, critical exponents $`(d2+\eta )/2`$ , and $`d\frac{1}{\nu }`$ for $`\rho =0.1,0.4`$, and $`1.0`$ are collected in Table 2.

We also estimated the shift exponent by taking $`\rho =0.005`$ and $`0.01`$. The critical points are found to be $`J_c=0.432507`$ and $`0.428593`$ respectively, which lead to a value $`1.773`$ for the shift exponent that consists with the theoretical prediction $`1.75`$ with $`1.3\%`$ deviation.

## Discussion and conclusions

Our numerical results suggest that the TLIM model is in the Ising universality class. That is, for a finite interlayer coupling, the critical behaviour of the TLIM is the same as that of the one-layer system. Both $`<s>`$ and $`<\sigma >`$ are order parameters in the transition. Our observations strongly exclude a phase transition for the electric order, except at the decoupling point $`\rho =0`$.

The critical line of the TLIM is in fact not a fixed line since the interlayer coupling is a relevant interaction which should flow to infinity in renormalisation group transformations. Then the TLIM resembles the one-layer Ising model in the long-range critical behaviour. To confirm this picture, we also carry out the real-space renormalizational group transformation to the model. The lattice is divided into $`3\times 3`$ blocks. All nearest neighbour interactions are kept. The RG trajactories are plotted in Figure 10.

Two attractive fix points both having infinite $`\rho `$ are observed. They are the ordered and disordered fix point, respectively. The critical line seperates two attractive fix points and obviously it is not a fix line.

However, as the interlayer coupling is small, it needs a large scale transformation to reach the infinite fixed points, so the model has extraordinary large finite-size effects as we have observed in the numerical results.

The decoupling limit shows singularities due to the existence of the severe competition between the unstable decoupling fixed point and the stable infinite interlayer coupling fixed points which, however, are far away. This supported in some sense by the unusual sensitivity of the shift exponent on the couplings Hori97 ; Lipo98  and the strange layer symmetry breaking in the effective state at the decoupling limit Ange97 . Our estimate for the shift exponent is consistent with the theoretical prediction.

The unusual small-size effects will have important consequences on the local properties of the TLIM. It is manifested in the short-time critical behaviour that we recently observed in Monte Carlo simulations which will be reported elsewhere Li00 .

###### Acknowledgements.

The authors thank B. Zheng for helpful discussions.