text:
# Spin injection across magnetic/non-magnetic interfaces with finite magnetic layers

## I Introduction

Spin injection across interfaces is one of the crucial ingredients for the successful implementation of novel spintronic devices wolf ; als ; Schmidt . For instance, the spin-transistor proposal datta  relies on subjecting injected spin-polarized electrons to a controllable spin precession between the ferromagnetic source and drain. Hybrid ferromagnetic/semiconductor and dilute-magnetic-semicondutor/non-magnetic semiconductor junctions constitute basic interfaces in which to investigate spin-polarized transport. However, the efficiency of spin-injection through ideal FM/NMS interfaces turns out to be disappointingly small due to the large conductivity mismatch Schmidt1  between the metallic ferromagnet and the semiconductor. The use of spin-selective interfaces can significantly enhance injection efficiencies rashba . This can be accomplished by inserting spin-dependent tunnel barriers between the FM and the NMS layers motsnyi .

Particularly promising is spin injection from a dilute-magnetic semiconductor (DMS) into a non-magnetic semiconductor (NMS). Novel DMS/NMS junctions (i) minimize the conductivity mismatch problem and (ii) naturally incorporate spin dependency in the transmission process egues . As recently demonstrated, substantial spin injection can be achieved in these Mn-based heterostructures spin-pol . More recently, a novel large magnetoresistance effect has been observed in a DMS/NMS/DMS geometries mol . An interesting observation of Ref. mol  is the dependence of the magnetoresistance effect on the thickness of the DMS layers: the MR effect doubles with the DMS thickness.

Available formulas describing GMR-type effects in magnetic/non-magnetic junctions assume 1D geometries with semi-infinite magnetic layers Khaetskii . In this case, the expressions for the electrochemical potentials for the spin-up and the spin-down electrons contain only decaying exponentials in the magnetic regions formulas . The resistances in the parallel and antiparallel configurations are calculated assuming that the device extends for a spin-flip length into the magnetic contacts. In such an approach it is not clear how to properly take into account the voltage drop across the sample.

In this work we consider DMS/NMS/DMS and FM/NMS/FM 1D systems with *finite* magnetic layers, Fig. 1. As detailed below, we describe diffusive transport in these structures via the usual diffusion theory of van Son et al. vanson  with proper physical boundary conditions between the several magnetic/non-magnetic interfaces in the system. The formulas for the magnetoresistance obtained within our treatment are in general different from the ones previously obtained in literature, even in the experimentally relevant regime $`d\lambda _D`$ and $`x_0\lambda _N`$; $`\lambda _D,\lambda _N`$ denote the spin-flip lengths in the DMS and NMS layers, respectively, $`d`$ is the length of each DMS layer and $`2x_0`$ is the length of the NMS layer. However, in the case of a FM/NMS/FM system where the conductivity mismatch is large, this difference is not important since the correction we find is very small in this case. For DMS/NMS/DMS systems, on the other hand, the conductivity mismatch is small and the magnetoresistance effect can be large. Our formulas significantly deviate from the earlier ones in this case, especially in the regime where $`\lambda _D`$, which can be magnetic-field dependent, is comparable to $`d`$.

## II System and approach

Let us denote by 1,2, and 3 the regions corresponding to the DMS (or FM), NMS, DMS (or FM) layers, respectively. Note that in Fig. 1 $`x=0`$ corresponds to the center of the structure (NMS layer), $`x=\pm x_0`$ correspond to the DMS/NMS interfaces and $`x=\pm (x_0+d)`$ to the metal contact/DMS interfaces. Consider first the DMS/NMS/DMS case or, equivalently, a parallel configuration FM/NMS/FM. Our task is to find the resistance of the structure $`R=V/j`$, where $`j=j_{}+j_{}`$ is the total current through the structure and $`V`$ is an applied bias ($`+V/2`$ is the potential of the left metal contact and $`V/2`$ is the potential of the right contact).

Diffusive transport equations. Our starting point is the diffusive transport approach of Ref. vanson, , with the basic equations relating the spin-dependent electrochemical potentials and current densities

$$j_,=(\sigma _,/e)d\mu _,/dx.$$
(1)
which hold in all three regions 1, 2, and 3 (see Fig. 1) of our system; $`\sigma _,`$ denotes the conductivity of the spin up (or down) electrons in the respective layer. The difference in the spin-dependent electrochemical potential is vanson

$$\mu _{}\mu _{}=\lambda ^2\frac{d^2(\mu _{}\mu _{})}{dx^2},$$
(2)
where $`\lambda `$ is the spin-flip length ($`\lambda _D`$ or $`\lambda _N`$) in the corresponding layer. As usual, to solve the problem we first assume a general solution in each of the three regions and then match them at the many interfaces via proper boundary conditions. Below we outline this procedure, clearly stating the physical boundary conditions we use and the simplifying symmetries of the problem.

*General solution for* $`\mu _,`$*.*  Consistent with Eq. (2), we use the general expression for the electrochemical potentials in region 1 (DMS)

$$\frac{1}{e}\left(\begin{array}{c}\mu _{}\\ \mu _{}\end{array}\right)=(A_1x+B_1)\left(\begin{array}{c}1\\ 1\end{array}\right)+(C_1\lambda _D\mathrm{exp}(x/\lambda _D)+D_1\lambda _D\mathrm{exp}(x/\lambda _D))\left(\begin{array}{c}1/\sigma _{}^D\\ 1/\sigma _{}^D\end{array}\right).$$
(3)
Here $`\lambda =\lambda _D`$ is the spin flip length in the DMS layer. Note the presence of the exponentially-growing term in Eq. (3), since we consider a finite DMS layer. The corresponding expressions for the electrochemical potentials in the regions 2 and 3 are obtained from Eq. (3) by the replacements $`A_1A_{2,3}`$, $`B_1B_{2,3}`$, $`C_1C_{2,3}`$, and $`D_1D_{2,3}`$; in addition we should make $`\sigma _,^D\sigma =\sigma _N/2`$ and $`\lambda _D\lambda _N`$ in the NMS region 2. The chosen form of the solution guarantees conservation of the total current everywhere in the sample. From the symmetry of the problem in the parallel configuration considered here, it follows that the electrochemical potentials $`\mu _{}`$ and $`\mu _{}`$ should be odd functions of the coordinate $`x`$, see Fig. 1(b). In terms of the coefficients this means that $`A_3=A_1,B_3=B_1,C_3=D_1,D_3=C_1,B_2=0`$, and $`C_2=D_2`$. Hence we can consider only the boundary conditions at the left metal contact/DMS interface and left DMS/NMS interface.

*Boundary conditions.* At the metal contact/DMS interface we use the continuity of the *total* current and the equality of the electrochemical potentials: $`\mu _{}=\mu _{}=eV/2`$. The condition $`\mu _{}=\mu _{}`$ is consistent with assuming a metal contact with infinite conductivity. Thus we obtain the following set of equations:

$$C_1+D_1\mathrm{exp}(2(x_0+d)/\lambda _D)=0,$$
(4)
$$A_1(x_0+d)+B_1=V/2,$$
(5)
and

$$j=A_1(\sigma _{}^D+\sigma _{}^D).$$
(6)
At the DMS/NMS interface we use the continuity of the electrochemical potentials, $`\mu _{}(x_0^{})=\mu _{}(x_0^+)`$, $`\mu _{}(x_0^{})=\mu _{}(x_0^+)`$, the conservation of the total current $`j=j_{}+j_{}`$, and conservation of the $`j_{}`$ component across the interface. Then we have the following additional set of equations

$`A_1x_0+B_1+\left(\begin{array}{c}\lambda _D/\sigma _{}^D\\ \lambda _D/\sigma _{}^D\end{array}\right)(C_1\mathrm{exp}(x_0/\lambda _D)+D_1\mathrm{exp}(x_0/\lambda _D))=A_2x_02C_2\lambda _N\mathrm{sinh}(x_0/\lambda _N)/\sigma ,`$ (9)
$`A_1\sigma _{}^D(C_1\mathrm{exp}(x_0/\lambda _D)D_1\mathrm{exp}(x_0/\lambda _D))=A_2\sigma 2C_2\mathrm{cosh}(x_0/\lambda _N),`$ (10)
and

$$j=2\sigma A_2.$$
(11)
Hence we have seven equations, Eqs. (4)–(11), for seven unknown quantities $`A_1,A_2,B_1,C_1,C_2,D_1`$, and $`j`$. By solving these equations we can determine the resistance of the structure. This we do next.

## III Results

For a FM/NMS/FM and a DMS/NMS/DMS system (“parallel” configuration), we find for $`R_p=V/j`$

$`R_p`$ $`=`$ $`{\displaystyle \frac{2x_0}{\sigma ^N}}+{\displaystyle \frac{2d}{\sigma ^D}}+`$ (12)
$`{\displaystyle \frac{2\beta ^2}{(1\beta ^2)\frac{\sigma ^D}{\lambda _D}\mathrm{coth}(\frac{d}{\lambda _D})+\frac{\sigma ^N}{\lambda _N}\mathrm{coth}(\frac{x_0}{\lambda _N})}},`$
where $`\beta =(\sigma _{}^D\sigma _{}^D)/(\sigma _{}^D+\sigma _{}^D)`$, $`\sigma ^N=2\sigma `$ is the total conductivity in the NMS layer and $`\sigma ^D=\sigma _{}^D+\sigma _{}^D`$ is the conductivity in the DMS layer; $`d`$ is the length of each DMS layer and $`2x_0`$ is the length of the NMS layer. The term proportional to the magnetic-layer width $`d`$ in Eq. (12) is not present in earlier formulas in the literature. Note that the term proportional to $`\beta ^2`$ in Eq. (12) differs from that in Eq. (1) of Ref. mol  because it is $`d`$ dependent and also because of its distinctive $`x_0`$ dependence correction .

For the antiparallel configuration (opposite FM layers), we can similarly obtain an expression for the resistance $`R_{ap}`$. Interestingly, $`R_{ap}`$ can be obtained from $`R_p`$ in Eq.(12) by changing $`\mathrm{tanh}(x_0/\lambda _N)`$ to $`\mathrm{coth}(x_0/\lambda _N)`$

$`R_{ap}`$ $`=`$ $`{\displaystyle \frac{2x_0}{\sigma ^N}}+{\displaystyle \frac{2d}{\sigma ^D}}+`$ (13)
$`{\displaystyle \frac{2\beta ^2}{(1\beta ^2)\frac{\sigma ^D}{\lambda _D}\mathrm{coth}(\frac{d}{\lambda _D})+\frac{\sigma ^N}{\lambda _N}\mathrm{tanh}(\frac{x_0}{\lambda _N})}}.`$
In obtaining Eq. (13) we have also assumed that the spin-up and the spin-down conductivities of region 3 are $`\sigma _{}^D`$ and $`\sigma _{}^D`$, respectively. We emphasize that the term proportional to $`\beta ^2`$ \[Eqs. (12), (13)\] does not reduce to that in Ref. Schmidt1  for $`d\mathrm{}`$ since the boundary conditions used here are different.

Magnetoresistances. For metallic FM/NMS/FM systems we define the resistance change as usual, i.e., $`\mathrm{\Delta }R=R_{ap}R_p`$. The full expression is too lengthy to be shown here. However, in the physically relevant limit $`\lambda _Dd`$ and $`x_0\lambda _N`$ we find from Eqs. (12) and (13)

$$\frac{\mathrm{\Delta }R}{R_p}=\frac{\beta ^2(\sigma ^N\lambda _D/x_0\sigma ^D)^2}{(1\beta ^2)[(1\beta ^2)+\frac{\sigma ^N\lambda _D}{x_0\sigma ^D}+(\frac{\sigma ^N\lambda _D}{\lambda _N\sigma ^D})^2\frac{1}{(1\beta ^2)}]}.$$
(14)
In deriving the above equation we have approximated $`R_p`$ in the denominator by its dominant contribution $`2x_0/\sigma ^N`$ from the central NMS layer. Again, our Eq. (14) is different from the corresponding one \[Eq. (8)\] in Ref. Schmidt1 . However, as can be seen from comparison of Eq. (14) and Eq. (8) of Ref. Schmidt1 , this deviation is not really important for metallic ferromagnets. Indeed, $`\mathrm{\Delta }R/R_p`$ which follows from our new equations practically coincides with the quantity in Eq. (8) of Ref. Schmidt1  for realistic values of the parameters when $`\sigma ^N\lambda _D/\sigma ^D\lambda _N1`$, $`\sigma ^N\lambda _D/\sigma ^Dx_01`$, and for $`\beta `$’s not anomalously close to unity.

In DMS/NMS/DMS systems, on the other hand, only the parallel configuration is possible. In this case, we define the relevant resistance change in the system with respect to the zero magnetic-field resistance $`R_0`$, i.e., $`\mathrm{\Delta }R^D=R_pR_0`$. We find

$$\frac{\mathrm{\Delta }R}{R_0}^D=\frac{\beta ^2}{(1\beta ^2)\frac{\sigma ^D}{\lambda _D}\mathrm{coth}(\frac{d}{\lambda _D})+\frac{\sigma ^N}{\lambda _N}\mathrm{coth}(\frac{x_0}{\lambda _N})}\frac{1}{\left(\frac{x_0}{\sigma ^N}+\frac{d}{\sigma ^D}\right)}+\frac{d(1/\sigma ^D(H)1/\sigma ^D(0))}{\left(\frac{x_0}{\sigma ^N}+\frac{d}{\sigma ^D}\right)}.$$
(15)
Note that $`R_0`$ is simply the combined NMS/DMS layer resistances $`2x_0/\sigma ^N+2d/\sigma ^D`$, since no spin accumulation occurs for zero magnetic field. In Eq. (15) we have indicated explicitly the additional magnetoresistance arising from the DMS (or FM) regions through $`\sigma ^D(H)`$.

Note that $`\mathrm{\Delta }R^D/R_0`$ goes to zero for either $`d\mathrm{}`$ or $`x_0\mathrm{}`$. The exact expression for $`\mathrm{\Delta }R/R_p`$ obtained from Eqs. (12) and (13) \[not the approximate one in Eq. (14)\] also vanishes in these limits. Note that this vanishing of the resistance change is symmetrical in $`d`$ and $`x_0`$ as expected. Earlier formulas in the literature do not show this symmetry since they assume different boundary conditions.

In contrast to metallic FM/NMS/FM systems, DMS/NMS/DMS structures can have conductivities of the DMS and NMS regions which are comparable. Besides, $`\lambda _D`$ can be of the same order as $`d`$ in some experiments. Experimental results in this case, e.g., the d dependence of the magnetoresistance effect, should be interpreted with the use of Eq.(15). As we mentioned earlier, the magnetoresistance effect observed in Ref. mol  doubles with the DMS thickness. Let us find out whether our 1D formula can explain this observed thickness dependence. Note that the d-dependence in Eq. (15) appears through the $`\mathrm{coth}(d/\lambda ^D)`$ term (spin accumulation effect). Since $`\sigma ^D=\sigma _{}^D+\sigma _{}^D`$ can in general be magnetic-field dependent, there is also a possible d-dependence arising from the second term which describes the resistances of the DMS leads.

In order to investigate in more detail the d-dependence in $`\mathrm{\Delta }R^D/R_0`$, we need to know the magnetic-field dependences of the quantities entering $`\mathrm{\Delta }R^D/R_0`$ \[Eq. (15)\], particularly the magnetic field dependence of the spin-dependent conductivities $`\sigma _,^D`$. Unfortunately, these dependences are not known. However, we can still gain some insight into the problem by neglecting the magnetic field-dependence of the conductivities. In this case, the second term in Eq.(15) is zero and we can investigate the expected d-dependence arising from the spin-accumulation effect alone. Note that the quantity $`\beta `$ is null in the absence of a magnetic field and increases as a function of it. Hence, we plot $`\mathrm{\Delta }R^D/R_0`$ as a function of degree of spin-polarization $`\beta `$ for different values of $`d`$. Figures 3(a), 3(b) show plots of $`\mathrm{\Delta }R^D/R_0`$ vs $`\beta `$ for different parameters. We have chosen the values of $`d,x_0`$ and $`\sigma ^D/\sigma ^N=1/3`$ corresponding to the real experimental situation of the Ref. mol . From these plots we can see that even though the absolute value of the magnetoresistance is comparable to the experimental value, the d-dependence is not as pronounced as in the experiment of Ref. mol . It is indeed true that the magnetoresistance effect is larger for increasing $`d`$. However, this behavior only appears at values of $`\lambda _D`$ comparable to the thickness $`d`$ and even for $`\lambda _D=0.5`$ $`\mu `$m $`>d`$, Fig. 3(b), the effect is way too small to explain the experimental data. We conclude that the spin accumulation effect by itself does not explain – with realistic parameters – the observed d-dependence of the data in Ref. mol .

Let us mention here the possibility to explain the experiment by the magnetoresistance of the leads, see the second term in Eq. (15). Control experiments on a layer of the DMS injector material patterned into a Hall geometry showed $`1\%`$ negative MR mol . However, the transport problem we deal with here can be shown to reduce to that of a 1D system Khaetskii . Note that in the case when the magnetic field dependence of $`\sigma ^D`$ is well pronounced and for realistic values $`x_0d\lambda _D`$ and for $`\beta 1`$ the first term in Eq. (15) is estimated as $`\lambda _D/\sigma ^D`$ while the second term is of the order of $`d/\sigma ^D`$. Thus in the case when the magnetoresistance of the leads is not zero the $`d`$-dependence associated with it can be more pronounced than the one associated with the spin accumulation effect.

Some limiting cases. From the general equations we derived above, we can obtain expressions for some interesting limiting cases. If $`\sigma _{}^D0`$, i.e. $`\beta 1`$, we find from Eq. (12)

$$R_p=\frac{2x_0}{\sigma ^N}+\frac{2d}{\sigma _{}^D}+2\frac{\lambda _N}{\sigma ^N}\mathrm{tanh}(x_0/\lambda _N).$$
(16)
If, in addition, $`\lambda _N\mathrm{}`$, then

$$R_p=\frac{2x_0}{\sigma ^N}+\frac{2d}{\sigma _{}^D}+\frac{2x_0}{\sigma ^N}.$$

Hence the resistance of the central part (NMS) doubles as expected. Note that the resistance of the DMS part does not necessarily double.

Current spin polarization. We can also easily determine some other relevant quantities in our system. For example, the degree of spin polarization $`\alpha (x)`$ of the current at the DMS/NMS (or FM/NMS) interface ($`x=x_0`$) is equal to

$$\alpha (x_0)=\frac{j_{}j_{}}{j_{}+j_{}}=\frac{\beta }{1+\frac{\sigma ^D}{\sigma ^N}(1\beta ^2)\frac{\lambda _N\mathrm{tanh}(x_0/\lambda _N)}{\lambda _D\mathrm{tanh}(d/\lambda _D)}},$$
(17)
for parallel alignment. This result coincides exactly with Eq. (7) of Ref. Schmidt1  in the limit $`x_0\lambda _N,\lambda _Dd`$. For metallic ferromagnetic contacts we can have antiparallel alignment as well. Note that $`\alpha (x)`$ is an even (odd) function of $`x`$ for parallel (antiparallel) alignment.

For completeness, we give below some other expressions for $`\alpha (x)`$. At $`x=\pm (x_0+d)`$ (even function) we have

$$\alpha (x=x_0d)=\beta \frac{\beta }{\mathrm{cosh}(d/\lambda _D)[1+\frac{\sigma ^N}{\sigma ^D}\frac{1}{(1\beta ^2)}\frac{\lambda _D\mathrm{tanh}(d/\lambda _D)}{\lambda _N\mathrm{tanh}(x_0/\lambda _N)}]}.$$
(18)
Note that for $`d\lambda _D`$ $`\alpha (x=x_0d)`$ approaches the maximal possible value $`\beta `$. In addition

$$\alpha (x=0)=\alpha (x=x_0)\frac{1}{\mathrm{cosh}(x_0/\lambda _N)}.$$
(19)
Thus in the middle of the sample the spin current polarization is exponentially suppressed. In the case of antiparallel configuration the corresponding values of $`\alpha `$ at $`x=x_0`$ and $`x=x_0d`$ can be obtained from formulas (17) and (18), respectively, by replacing $`\mathrm{tanh}(x_0/\lambda _N)`$ for $`\mathrm{coth}(x_0/\lambda _N)`$.

In conclusion, we have derived the expressions for the resistances of FM/NMS/FM (or DMS/NMS/DMS) structures using physical boundary conditions at the corresponding interfaces. We have considered magnetic layers (DMS or FM) of finite length $`d`$. When the magnetoresistance of the contacts is negligible (i.e., when $`\sigma _D`$ is not magnetic field dependent), we find that the $`d`$ dependence of the spin-accumulation contribution to the magnetoresistance is not enough to explain the experimental results of Ref. mol . On the other hand, in systems where $`\sigma _D`$ is magnetic field dependent our formulas allow for a sizable $`d`$ dependence arising from the magnetoresistance of contacts. Interestingly, for realistic parameters $`x_0d\lambda _D`$ and with $`\beta 1`$ the magnetoresistance of the contacts is larger than that from the spin accumulation effect.

This work was supported by the Swiss NSF, NCCR Nanoscience Basel, DARPA, ARO, ONR and the DARPA Spins program.


format: text

text:
# 1 Introduction

## 1 Introduction

The spectroscopy of excited $`B`$ and $`D`$ mesons is important for our understanding of QCD. Moreover, as stressed recently by Rosner , $`B^{}`$ states also have important applications to CP studies of neutral $`B`$ mesons by the identification of their flavour ($`b`$ versus $`\overline{b}`$) through the decay chain $`B^{}B^0\pi ^\pm `$. Hence narrow $`B^{}`$ resonances will be valuable for this.

In the heavy quark limit, the $`\overline{Q}q`$ meson, which we refer to as a ‘$`B`$’ meson, will be the ‘hydrogen atom’ of QCD. Since the meson is made from non-identical quarks, charge conjugation is not a good quantum number. States can be labelled by L±, where the coupling of the light quark spin to the orbital angular momentum gives $`j_q=L\pm \frac{1}{2}`$. In the heavy quark limit these states will be doubly degenerate since the heavy quark spin interaction can be neglected, so the P- state will have $`J^P=0^+,1^+`$ while P+ has $`J^P=1^+,2^+`$, etc.

This spectrum has been studied comprehensively by lattice methods in the quenched approximation  with a rather coarse lattice spacing. With $`N_f=2`$ flavours of dynamical quark, SESAM  have explored the P- excited state. Lattice studies using NRQCD have also explored the heavy-light spectrum for $`b`$ quarks, mainly using quenched lattices .

In the heavy quark effective theory, the leading order is just the static limit. The next correction will be of order $`1/m_Q`$ and will be relatively small for $`b`$ quarks, but larger for $`c`$ quarks. One way to predict the spectrum for $`b`$ quarks is to interpolate between charmed states, where the experimental spectrum is known, and the static limit obtained by lattice QCD assuming a dependence as $`1/m_Q`$. Thus the splittings among $`B`$ mesons should be approximately 0.33 of those among the corresponding $`D`$ mesons.

The striking discovery that the $`c\overline{s}`$ states with $`J^P=0^+`$ and $`1^+`$ have very narrow widths  raises the question of whether the corresponding $`b\overline{s}`$ states will also be narrow. The main reason for the narrow width of $`D_s`$ mesons is that the transition to $`DK`$ is not energetically allowed (for the 2317 MeV state) or the state is close to threshold (for the 2457 MeV state). Thus the only allowed hadronic decay proceeds via isospin-violation (since $`m_dm_u`$) to $`D_s`$$`\pi `$ and will have a very small width. Likewise, if the equivalent $`b\overline{s}`$ states are close to or below the $`BK`$ threshold, then they will be very narrow. One of our main tasks will be to determine the energy of these excited states as accurately as possible to check this.

As well as exploring this issue of great interest to experiment, we determine the excited state spectrum of the heavy-light system as fully as possible. This will help the construction of phenomenological models and will shed light on questions such as whether there is an inversion of the level ordering (with L+ lighter than L-) at larger L or for radial excitations as has been predicted . We also compare with chiral models .

## 2 Lattice evaluation

We investigate the heavy-light meson spectrum from lattice QCD using static heavy quarks. Previous lattice studies have explored  the full spectrum (i.e. S, P-, P+, D-, D+, F) in quenched QCD. There has also been a recent determination of the P- excitation energy in full QCD .

Here we present a range of different lattice studies: with different spatial volumes, lattice spacings and light quark masses — see Tables 12. We follow the all-to-all methods used in the static-light lattice study of Michael and Peisa . Keeping their parameters, we first use a larger spatial size of lattice to check for finite size effects — Q1 vs. Q3. We are also able to correct the assignments of D+ and D- states in their work, see Q1 and Q2 in Table 2. Our major study involves using lattice configurations  which include $`N_f=2`$ flavours of sea-quark, with two different lattice spacings. We only use the unitary points, namely those with valence light quarks of the same mass as the sea quarks. The details are collected in Table 1.

To extract mass values, we use operators with the appropriate representations of the cubic group (as described by ) with different degrees of non-locality. We find that the D+- operator approximately gives the spin-average of the D- and D+ levels, so we can interpret the F+- operator as representing an average of the two F levels. Our choice of operators enables us to determine $`N\times N`$ matrices of correlations for each case, where $`N`$ can vary from 2 to 5. We then perform a fit to these correlations over a suitable $`t`$-range with a number of states allowed. The requirement is then that the $`\chi ^2`$ per degree of freedom is reasonable (not much greater than 1). We always use at least 2 states so that we have a reliable estimate of the ground state mass. To extract the first excited state as well, it is preferable to use at least a 3 state fit. We check that using a subset of our largest matrix of correlators, using different $`t`$-ranges, using one more or less state, etc. gives stable results.

To compare different lattice simulations, we form the dimensionless combination of $`r_0m`$, where $`m`$ is a mass or energy, and where $`r_0/a`$ is determined relatively accurately from the static quark potential. Our results are shown in Tables 12 and some comparisons for the P- and D-wave states are shown in Figs. 1-3 versus lattice spacing and versus quark mass.

In order to relate our lattice results to experiment we have to discuss three different extrapolations.

(i) Finite Size Effects. The lattice spatial volume should be large enough. There are several related criteria: the wavefunctions of the heavy light mesons should be small compared to the spatial size $`L_S`$, the exchange of the lightest particle (the pseudoscalar meson) around the periodic boundary should be small and mixing of the heavy-light mesonic states with two body states (e.g. $`B\pi `$ where the $`\pi `$ has a low momentum) should be small.

We can estimate the size of the heavy-light mesons from the Bethe-Salpeter wavefunctions measured for ground and excited states  and also from the more physical charge and matter distributions evaluated for the ground state (1S) meson . These results suggest that 2 fm is a sufficient size for quenched evaluations, which is confirmed by our results which extend the spatial volume of the previous quenched measurements but do not show any statistically significant differences.

Dynamical fermion configurations are more sensitive to finite size effects since more loop effects are included, in particular pion exchange around the boundary becomes important. The leading correction  for the ground state is a relative energy shift of order $`ce^{mL_S}`$, where $`m`$ is the pseudoscalar mass and $`c`$ a coefficient given by the $`B^{}B\pi `$ coupling. For the excited states, the possibility of the decay to (or mixing with) nearby two-body energy levels becomes relevant. The only excited state that couples to a low-lying two-body energy level is for the P- which has a mixing with $`B\pi `$ where the pion has momentum zero. Thus we expect an enhanced finite size effect may arise for P-. We investigate this by using two spatial sizes (called DF1 and DF2, with $`L_S`$ of 1.7 fm and 2.3 fm, corresponding to $`m_\pi L_S`$ = 6.7 and 9.0 respectively). We see some sign of a shift for D+ and F but it is not very significant statistically.

Our data set with the finest lattice spacing has a relatively small volume (1.6 fm with $`m_\pi L_S`$=4.5, 6.5) and, for this situation, some evidence of finite size effects for the nucleon has been presented . Some of our results from this finest lattice spacing are significantly different from the larger volume results described above. Since the order $`a`$ formalism used is different (NP clover rather than tadpole-improved clover) we cannot select between finite size effects or lattice artifacts (order $`a`$ etc.) as a cause. We would need to use a larger spatial volume at these parameters to evaluate finite size effects fully.

(ii) Quark mass dependence. Let us first discuss the dependence on the light valence quark mass. Experimental data on the heavy light mesons with $`c`$ or $`b`$ quarks suggest that there is little quark mass dependence of excitation energies (i.e. energy differences to the ground state pseudoscalar meson) when going from strange quarks to lighter quarks. For instance the mass splittings $`D^{}(1^{})D`$ and $`D_s^{}(1^{})D_s`$ are 141 and 144 MeV respectively, while $`D^{}(2^+)D`$ and $`D_s^{}(2^+)D_s`$ are 593 and 604 MeV respectively .

We can also explore this on a lattice and quark masses (characterised by $`[r_0m(0^+)]^2`$ where $`3.4`$ corresponds to strange quarks as discussed below) are varied in the range of 0.6 to 1.5 times the strange quark. This is shown for the P-wave states in Fig. 2 which confirms that there is no significant slope. This means that interpolation to the strange quark mass is not delicate in any way. Extrapolation to light valence quarks is less straightforward and one issue that must be addressed is that some of the excited heavy-light mesons are unstable to strong decay. Since an excited L± state will have a decay to $`B\pi `$ with angular momentum given by L$`\pm 1`$, only the P- state can decay in an S-wave which then gives the lowest threshold energy on a lattice because the pion can have momentum zero. In our case, because of the discrete momentum and unphysical light-quark mass values, we do not have any open decay channels in our lattice evaluation, but they will open when extrapolating in light quark mass and to large spatial volume.

The issue of the extrapolation in the sea quark mass is difficult to resolve. We cover the range from no sea quarks (i.e. quenched) to $`N_f=2`$ flavours of sea quark with mass corresponding to 0.6 times the strange quark. The evaluation with even lighter sea quarks is computationally too demanding.

(iii) The continuum limit. It is feasible to study the continuum limit in quenched studies, but for dynamical fermions we only have access to a relatively narrow range of lattice spacing ($`a`$ from 0.15 fm to 0.1 fm). To make best use of this limitation, we use an order $`a`$ improved clover formulation of the fermion action. The coarser lattice has a tadpole-based improvement coefficient while the finer lattice uses a non-perturbatively improved value. Because of this difference in formalism, it is not straightforward to extrapolate from these two data sets to the continuum limit. We take this into account in assigning errors.

Lattice spectrum.

We average the values discussed above of the various excitation energies, weighting relatively more small lattice spacing, large volume and quark masses close to strange. Thus we obtain $`r_0\mathrm{\Delta }E`$ of 1.07(7) for P-; 1.33(13) for P+. The next excited level is the 2S which is at 1.25(-13,+50), this is an average based on the larger volume studies but with the error reflecting our results at finer lattice spacing. For the D-waves there is also a large spread so we quote a range: for D- from 2.2 to 3.1, while for D+ from 2.2 to 3.5. For the F-wave, we only have an operator which excites both F+ and F- so that our result is for an average of these two states, with an excitation energy around 3.4 to 4.4.

We need a value of the scale $`r_0`$ appropriate to light quark spectroscopy, since the dynamics of the light quark is the main aspect of heavy-light mesons. Thus we do not use values of $`r_0`$ from heavy-heavy studies (which tend to give somewhat smaller values) but an average of those from light-light mesons which span the range of 0.5 to 0.55 fm, namely 0.525$`\pm `$0.25, for a discussion see Ref. . This value of $`r_0`$, combined with the estimate of the mass of the pseudoscalar meson made from strange quarks  of 687 to 695 MeV yields $`r_0m_\pi 1.84`$ which sets the scale for the strange quark.

In our application to the heavy-light mesons with $`N_f=2`$ flavours of sea-quark, we have used valence quarks identical to the sea-quarks, which is the case where the theory is fully unitary. This can be interpreted in two ways — firstly as applying to the spectrum of excited $`b\overline{n}`$ states (where $`n`$ is $`u`$ or $`d`$) with quark masses heavier than the physical values. Indeed in the Tables we give the pseudoscalar meson mass obtained by combining these light quarks. On the other hand, for our application to the $`b\overline{s}`$ system, we have also used valence quark masses identical to the sea-quark mass. In the real world, however, there is only one flavour of strange quark, so we can interpret our results as from one flavour of strange valence quark propagating in a sea with two flavours of light quarks whose mass happens to correspond to the valence quark mass. This is effectively treating the strange quark as partially quenched and further studies would be needed to treat fully all three flavours of light quark in the sea.

In principle one can calculate corrections to the heavy quark limit from the lattice, as discussed later. Here, however, we adopt a more modest strategy and make partial use of experimental data. Thus to interpolate to $`b`$ quarks we combine our results in the static limit with experimental data  for the $`c\overline{s}`$ system as shown in Fig. 4. For the P- state the experimental excitation energy for charm quarks is 349 MeV while we obtain for static quarks 404(31) MeV. Thus the interpolation to $`b`$ quarks involves only small shifts - leading to 386(31) MeV. This is close to the threshold for decay emitting a kaon (a mass gap of 404 MeV) and probably below it. So we do expect this $`b\overline{s}`$ scalar meson to be very narrow, as was found for the $`c\overline{s}`$ counterpart . The associated axial meson at 434(31) MeV above the $`B_s`$ will be close to the $`B^{}K`$ threshold (at 450 MeV) and should also be very narrow. The P+ states lie above the $`BK`$ threshold but since these states decay in a D-wave, the centrifugal barrier effects may cause them to have narrow widths.

For the 2S, D, F states, we do not have any $`c\overline{s}`$ counterpart available from experiment to allow this interpolation. Assuming, however, that the slopes versus $`m_c/m_Q`$ are similar to those for the P-wave case, then the static energy values will be a good approximation to those for $`b`$ quarks. Again the 2S pseudoscalar (and vector) states could be sufficiently light that they lie close to the $`B^{}K`$ ($`BK`$ for vector) threshold at 450 MeV (404MeV ) and so are very narrow.

The only experimental observation  of an excited $`B_s`$ state is the $`B_s`$(5850) which lies 483 MeV heavier than the $`B_s`$ and has a width of 47(22) MeV. This mass value is indeed in the region where we predict a rich spectrum of excited $`b\overline{s}`$ states.

As we see no sign of a significant light quark mass dependence in our excitation energies, we can use our results to predict the spectrum of excited $`b\overline{n}`$ states albeit with a somewhat larger systematic error from the extrapolation to light quarks which we are assuming to be a constant. The only experimental observation  of an excited $`B`$ state is the $`B_J^{}`$(5732) which lies 419 MeV heavier than the $`B`$ and has a width of 128(18) MeV. This may indeed be composed of several states. The mass value is indeed in the region where we predict a rich spectrum of excited $`b\overline{n}`$ states, even though they should not be especially narrow since the $`B\pi `$ and $`B^{}\pi `$ decay channels are open.

As well as predicting the spectrum, lattice methods can be used to evaluate decay amplitudes  and this is feasible for the heavy-light systems too .

## 3 Discussion

We first discuss the issue of the theoretical relationship between the static limit and realistic heavy quarks. Then we can use this discussion to organise our comparison with other lattice determinations of the heavy-light spectrum.

A precise description of heavy-light mesons is provided by the heavy quark effective theory. The leading ($`1/m_Q`$) corrections  to the static (i.e. heavy quark) limit arise from two sources: kinetic and magnetic terms. The magnetic contribution splits each static energy level H (with total light quark angular momentum $`j_q=L\pm \frac{1}{2}`$, called L± above) into two with total angular momentum $`j_1=j_q+\frac{1}{2}`$ and $`j_2=j_q\frac{1}{2}`$. They have masses given by

$`m_{H_1}`$ $`=`$ $`m_Q+\mathrm{\Lambda }_H+{\displaystyle \frac{\lambda _{H,K}}{2m_Q}}+(2j_1+1){\displaystyle \frac{\lambda _{H,B}}{2m_Q}}`$ (1)
$`m_{H_2}`$ $`=`$ $`m_Q+\mathrm{\Lambda }_H+{\displaystyle \frac{\lambda _{H,K}}{2m_Q}}(2j_2+1){\displaystyle \frac{\lambda _{H,B}}{2m_Q}},`$ (2)
where $`m_Q`$ is the mass of the heavy quark and $`\mathrm{\Lambda }_H`$ is the binding energy. On the lattice there is a self-energy term proportional to $`1/a`$, but as we only discuss mass differences, this will cancel.

Here $`\lambda _{H,K}`$ arises from the insertion of the heavy quark kinetic energy for state H i.e.

$$\lambda _{H,K}=H\overline{Q}D_T^2QH.$$
(3)
As the transverse kinetic energy is expected to be positive, this implies that $`\lambda _{H,K}`$ should be positive also. However, it is the difference of kinetic energies between states that we need. In a simple approach with a confining potential, the excited state would have larger kinetic energy than a ground state, so the mass differences between the P- and S-wave states would increase as $`m_Q`$ is decreased, but this is only a qualitative indication.

The coefficient $`\lambda _{H,B}`$ arises for state H from the insertion of the $`𝝈\mathbf{.}𝑩`$ term, where $`𝝈`$ is the heavy quark spin and $`𝑩`$ is the chromomagnetic field from the light quark. For the S-wave states ($`B^{}`$, $`B`$) , the $`\lambda _{S,B}`$ parameter can be estimated from the experimental $`B^{}`$ to $`B`$ mass splitting

$$\lambda _{S,B}\frac{1}{4}(M_B^{}^2M_B^2)=0.12\text{GeV}^2.$$
(4)
The NRQCD lattice formalism allows these $`1/m_Q`$ expressions to be evaluated. The results from several recent studies  show that essentially all excitation energies increase as $`m_Q`$ is decreased from the static limit. This is what would be expected from the kinetic energy correction above. A note of caution, however, is that the magnetic contribution from these studies can be compared with experimental data on the $`B^{}`$, $`B`$ splitting and underestimates it by almost a factor of two. This suggests that the NRQCD, as currently implemented, is not reproducing the magnetic contribution accurately. Thus predictions from NRQCD of hyperfine splitting may be underestimated. In the NRQCD method one does not take a continuum limit, but the approach can be systematically improved by including more terms in the effective action and by computing the coefficients of these terms (such as $`𝝈\mathbf{.}𝑩`$) non-perturbatively, although this is yet to be carried out to a level such that systematic errors on the hyperfine splittings can be established.

Another way to estimate the $`1/m_Q`$ corrections from lattice studies is to compare static results, such as ours, with results from relativistic propagating quarks, where a continuum limit may be taken. Here recent results for charm quarks  do give a spectrum of $`c\overline{s}`$ mesons substantially in agreement with experiment and hence support the pattern of $`1/m_Q`$ corrections we show in Fig. 4 Note that Bali  used an estimate of $`1/m_Q`$ effects by taking the difference of the quenched results in the static limit (from Michael and Peisa ) and for charm (from Boyle ) and using this to correct the $`N_f=2`$ static result from SESAM . This different procedure explains why his result for the scalar P-wave meson is heavier than ours (by $`2\sigma `$) even though he obtains a similar value for the P- energy excitation in the static limit.

We illustrate some of the above discussion by presenting a compilation of relevant lattice results in Fig. 5. Some older lattice calculations of the mass spectrum of P-wave heavy-light mesons have been reviewed recently . Improved lattice calculations with reduced systematic and statistical errors are required to get definitive answers .

Having discussed the heavy quark effective theory, we now discuss the implications of our results for other models of heavy-light mesons.

A traditional way to understand such spectra would be using a quark model with an underlying potential description . This is not strictly justified for a light quark, but may be of qualitative use. For the experimentally observed excited $`D_s`$ states, it is difficult to understand why the hyperfine splitting is sufficiently big to give a $`J^P=0^+`$ meson which is so light in such an approach . Our results enable us to discuss the possible inversion of the level ordering (with L+ lighter than L-) at larger L or for radial excitations. This inversion has been predicted  from consideration of the spin-orbit force, which at larger separation would come more from the confining interaction than the short-ranged contribution from gluon exchange. We find no evidence of a sign change in the spin-orbit splitting of P- or D-waves. Thus conventional short-distance spin-orbit effects are still relevant up to radii appropriate for D-wave states.

It is possible to discuss chiral symmetry in the heavy quark limit. This allows relationships  between energy levels and also predictions for coupling strengths. A stronger assumption of the form of chiral symmetry breaking allows to obtain results away from the static limit, such as that the $`1^{}`$ to $`0^{}`$ splitting is the same as that from $`1^+`$ to $`0^+`$. Chiral symmetry in the heavy quark limit relates the S state to P- and the D- state to P+, etc. This does not seem to be a very good approximation: the spectrum is closer to being dependent on L alone. Indeed our spectrum shows an approximately linear rise in excitation energy with L.

## 4  Conclusion

We have used lattice QCD to explore the spectrum of heavy-light mesons. Our results are evaluated for static quarks but they are very relevant to $`b`$ quarks ($`B^{}`$ states) as we have argued. We have concentrated our studies on light quarks which are of the mass of a strange quark, although we find that excitation energies are consistent with being independent of the light quark mass, and hence will apply also to light quarks which are $`u`$ and $`d`$.

In our lattice studies, we have pushed towards light sea quarks, towards small lattice spacing and towards large volume, but not towards all three requirements simultaneously. This leaves some room for systematic errors in our predictions. These can be reduced by further studies with increased computational resources.

We have determined the spectrum up to F-waves and including radial excitations, see Fig. 6. This gives a rich texture for model building of heavy-light mesons. We find no evidence of a sign change in the spin-orbit splitting of D-waves. Thus conventional short-distance spin-orbit effects are still relevant up to radii appropriate for D-wave states. Rather than the pattern given by chiral symmetry (which relates S to P- and D- to P+ etc. ) we find a spectrum which is closer to being dependent on L alone. Indeed we see an approximate linear rise in excitation energy with L, up to L=3, as 0.45 L GeV, reminiscent of Regge or string models.

We have discussed corrections to the heavy quark limit appropriate to $`B^{}`$ states and have used experimental data on $`D^{}`$ states to establish this. Our results for the P-wave excitations confirm those obtained previously  that the excitation energies are relatively small. This implies that the $`b\overline{s}`$ P-wave states will be close to the lightest hadronic decay thresholds (namely $`BK`$ and $`B^{}K`$). The P- states ($`J^P=0^+,1^+`$) have an S-wave decay but are light enough that there is little or no phase space for decay. The P+ states ($`J^P=1^+,2^+`$) are heavier but have D-wave decays and so will also have narrow widths since centrifugal barrier effects will reduce them. We also see evidence that the 2S states ($`J^P=0^{},1^{}`$) are close to the lightest thresholds and so may be narrow too. Our results for $`b\overline{s}`$ states are summarised in Fig. 4 and in Table 3. This raises the prospect that there will be many (up to 6) narrow excited $`B_s`$ states to be found experimentally.

Since we see no significant dependence of the excitation energies on light quark mass, our predictions from the Table 3 can also be used as estimates for orbital and radial $`B^{}B`$ excitation energies, although these mesons will not be especially narrow since the $`B\pi `$ and $`B^{}\pi `$ thresholds are open.

Acknowledgments

The authors wish to thank the Center for Scientific Computing in Espoo, Finland and the ULGrid project of the University of Liverpool for making available ample computer resources and thank Chris Maynard for helpful advice. Two of us, AMG and JK, acknowledge support by the Academy of Finland contract 54038 and useful discussions with Petrus Pennanen. JK thanks the Finnish Cultural Foundation and the Magnus Ehrnrooth Foundation for financial support.


text:
# I Zw 1 observed with XMM-Newton

## 1 Introduction

Aside from being the prototype Narrow-Line Seyfert 1 galaxy (NLS1), I Zw 1 is an important object having influence in many areas of AGN astronomy. Its strong and narrow optical/UV emission lines minimize blending effects and allow I Zw 1 to be used as a template for optical (Boroson & Green 1992) and UV (Vestergaard & Wilkes 2001; Laor et al. 1997) FeII emission. Taking its luminosity into consideration (MB $``$ $``$23.5; Schmidt & Green 1983), I Zw 1 is also defined as a quasar, likely making it the nearest quasar known ($`z`$ = 0.0611). Although its infrared luminosity falls just short of those of Ultraluminous Infrared galaxies (ULIRGs), I Zw 1 exhibits properties of both QSOs and ULIRGs, and hence Canalizo & Stockton (2001) consider it as a possible transition object from QSO to ULIRG. Scharwächter et al. (2003; see also Eckart et al. 2000) used VLT observations to examine the host properties of I Zw 1. They observe an interacting spiral galaxy with a young stellar population in the disc as well as young, hot stars, and supergiants in the bulge, which is indicative of strong starburst activity.

At X-ray energies, I Zw 1 displays many of the characteristics typically associated with NLS1. $`ROSAT`$ observations (Boller et al. 1996; Lawrence et al. 1997) found a steep 0.1–2.4 keV spectrum, described by a power-law of photon index $``$3, and notable variability (change by a factor of 1.5 in $``$6000 s). However, Leighly (1999a) found that an absorbed power-law, without a soft excess, was sufficient to fit the 0.5–10 keV $`ASCA`$ continuum. In addition, Leighly detected a very strong ($`EW`$ $``$ 1300 eV) Fe K$`\alpha `$ line (see also Reeves & Turner 2000), which suggested a nucleus-wide overabundance of iron in I Zw 1.

In this paper we present a 20 ks XMM-Newton observation of I Zw 1. We begin by describing the observation and the data reduction in the next section. Section 3 is devoted to the results from the Optical Monitor. In Section 4 we discuss the X-ray spectra observed with the EPIC. In Section 5 we concentrate on the timing properties. Lastly, in Section 6, we will conclude with a summary of our results.

## 2 Observation and data reduction

I Zw 1 was observed with XMM-Newton (Jansen et al. 2001) on 2002 June 22 during revolution 0464 for about 20 ks. During this time all instruments were functioning normally. The EPIC pn camera (Strüder et al. 2001) was operated in large-window mode, and the two MOS cameras (MOS1 and MOS2; Turner et al. 2001) were operated in small-window mode. All of the EPIC cameras used the medium filter. The two Reflection Grating Spectrometers (RGS1 and RGS2; den Herder et al. 2001) and the Optical Monitor (OM; Mason et al. 2001) also gathered data during this time. The Observation Data Files (ODFs) were processed to produce calibrated event lists using the XMM-Newton Science Analysis System (SAS v5.4.1). Unwanted hot, dead, or flickering pixels were removed as were events due to electronic noise. Event energies were corrected for charge-transfer losses. EPIC response matrices were generated using the SAS tasks ARFGEN and RMFGEN. Light curves were extracted from these event lists to search for periods of high background flaring. A short (lasting $`<`$ 300 s), small-amplitude flare (equivalent to less than 3% of the source count rate) was detected about 9.5 ks into the observation. The data for this time period were ignored during the analysis. The total amount of good exposure time selected was 18650 s and 21006 s for the pn and MOS detectors, respectively. The source plus background photons were extracted from a circular region with a radius of 35 arcsec, and the background was selected from an off-source region and appropriately scaled to the source region. Single and double events were selected for the pn detector, and single-quadruple events were selected for the MOS. The statistics were clearly dominated by the source at energies below 10 keV. The total numbers of good counts collected by the pn, MOS1, and MOS2 instruments were 137840, 45581, and 46829, respectively. The XMM-Newton observation provides a substantial improvement in spectral quality over the 28.3 ks exposure (92.4 ks duration) $`ASCA`$ observation in which $``$12000 counts were collected (Leighly 1999b).

High-resolution spectra were obtained with the RGS. The RGS were operated in standard Spectro+Q mode for a total exposure time of 20713 s. The first-order RGS spectra were extracted using RGSPROC, and the response matrices were generated using RGSRMFGEN.

The Optical Monitor collected data through the UVW2 filter (1800–2250 Å) for about the first 9 ks of the observation, and then it was switched to grism mode for the remaining time.

## 3 UV analysis

The Optical Monitor collected data in the fast mode through the UVW2 filter (1800–2250 Å) for about the first half of the observation. The rest-frame wavelength range observed through the UVW2 filter is roughly 1700–2120 Å. Besides continuum emission, Laor et al. (1997) attributed much of the emission in this spectral region to FeII and FeIII multiplets. There is also some strong emission from AlIII, and semi-forbidden species of SiIII\], CIII\], and NIII\].

In total, five photometric images were taken, each exposure lasting 1660 s. The apparent UVW2 magnitude is 14.36 $`\pm `$ 0.04 corresponding to a flux density of 9.94 $`\times `$ 10-15 erg s-1 cm-2 Å-1 which is about a factor of three lower than during the 1994 $`HST`$ observation (Laor et al. 1997). The approximately 2.5 hour light curve is flat with the exception of one data point which exhibits a 1$`\sigma `$ deviation from the average. We do not consider this event significant.

Taking advantage of the simultaneous UV/X-ray observations, we derive the nominal spectral slope, $`\alpha _{ox}`$, between 2500 Å and 2 keV. For the UV flux, we use the peak response of the UVW2 filter (2000 Å in the rest-frame). We note from the $`HST`$ spectrum, that the flux density at 2000 Å is approximately equal to the flux density at 2500 Å (Laor et al. 1997). By assuming a negligible continuum slope between 2000–2500 Å we are able to extrapolate the UVW2 flux at 2000Å to 2500 Å. For the model dependent X-ray flux (see Section 4) we used only the EPIC data which were collected simultaneously with the UV data (i.e. the first 10 ks of the observation). We determine a value of $`\alpha _{ox}`$ = 1.20 which is notably X-ray stronger than the value of 1.41 reported by Lawrence et al. (1997), and 1.56 reported by Brandt et al. (2000).111 Brandt et al. (2000) use a reference point of 3000 Å for the UV flux. It is also worth noting that at 3000 Å there is a relatively strong FeIII multiplet emission that will make I Zw 1 appear X-ray $`weaker`$ than if we were to consider only the continuum. The spectral index between 2500 Å and 2 keV is clearly variable with time, but always appears X-ray strong compared to the predicted value for radio-quiet quasars of this luminosity (Vignali et al. 2003).

## 4 X-ray spectral analysis

The source spectra were grouped such that each bin contained at least 40 counts. Spectral fitting was performed using XSPEC v11.2.0 (Arnaud 1996). Fit parameters are reported in the rest-frame of the object, although the figures remain in the observed-frame. The quoted errors on the model parameters correspond to a 90% confidence level for one interesting parameter (i.e. a $`\mathrm{\Delta }\chi ^2`$ = 2.7 criterion), and luminosities are derived assuming isotropic emission. The Galactic column density toward I Zw 1 is $`N_H`$ = 5.07 $`\times `$ 1020 cm-2 (Elvis et al. 1989). A value of the Hubble constant of $`H_0`$=$`70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and a WMAP-cosmology with $`\mathrm{\Omega }_M`$ = 0.3 and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = 0.7 have been adopted throughout.

Initially all three EPIC spectra were fitted separately to examine any discrepancies. The 3–10 keV spectra were first fitted with a single power-law modified by Galactic absorption. The pn and MOS2 are entirely consistent; the photon indices for the pn and MOS2 fits are 2.26$`{}_{0.06}{}^{}{}_{}{}^{+0.05}`$ and 2.18 $`\pm `$ 0.10, respectively. However, the MOS1 photon index is much flatter ($`\mathrm{\Gamma }_{MOS1}`$ = 1.95 $`\pm `$ 0.10) and does not lie within the 90% confidence range of the other EPIC instruments. This inconsistency in the MOS1 data was previously realised in an observation of 3C 273 by Molendi & Sembay (2003). Molendi & Sembay noted a difference in the MOS1 photon index compared to the other EPIC photon indices of $`\mathrm{\Delta }\mathrm{\Gamma }`$ $``$ 0.1. The difference in our observation is larger (0.13 $`<`$ $`\mathrm{\Delta }\mathrm{\Gamma }`$ $`<`$ 0.31). However, the pn data dominate the fits above 3 keV given the higher statistics afforded by this instrument. Including the 3–10 keV MOS1 data has little influence on the model results. For clarity, the figures are shown without the 3–10 keV MOS1 data, but the fit statistics will include all of the data. Cross-calibration uncertainties between pn and MOS appeared quite significant in our observation at energies below $``$ 0.7 keV; therefore, the EPIC data are only fitted together at energies above 0.7 keV. Between 0.7–3 keV, the data from all three instruments appears entirely consistent.

In Figure 1 we present the 0.3–10 keV EPIC spectra fitted with an absorbed power-law. The fit is reasonable given its simplicity ($`\chi ^2`$ = 1325.4/1076 $`dof`$); however, a number of features are immediately apparent. The total photoelectric absorption is nearly twice that expected from the Galactic hydrogen column density, indicating that I Zw 1 is intrinsically absorbed. An excess above the power-law is observed between 0.4–0.6 keV. The weak soft excess is marked on the red side by absorption which extends to the energy limit of the detector, and on the blue side by an absorption feature which dips below the power-law from 0.6–0.8 keV. Positive residuals between 6–7 keV likely suggest emission due to iron.

### 4.1 Fe K emission

We begin the spectral analysis by considering the 3–10 keV band alone. Fitting a power-law in this region results in positive residuals between 6–7 keV in the observed-frame, as is shown in Figure 1 ($`\chi ^2`$ = 333.4/295 $`dof`$). In addition, the fit demonstrates that the spectrum is steep compared to broad-line Seyfert 1s, with $`\mathrm{\Gamma }`$ = 2.24 $`\pm `$ 0.05. We attempt several model fits to describe the feature between 6–7 keV. A Gaussian model, diskline model (Fabian et al. 1989), and a Laor line model (Laor 1991) are all relatively successful in fitting the data, but none could be ruled superior over the others. In all cases, the addition of a line model results in further steepening of the power-law continuum ($`\mathrm{\Gamma }`$ $``$ 2.3). The best-fitting parameters to the simplest line model (i.e. Gaussian profile) are $`E`$ = 6.69$`{}_{0.33}{}^{}{}_{}{}^{+0.14}`$ keV, $`\sigma `$ = 316$`{}_{136}{}^{}{}_{}{}^{+268}`$ eV, and $`EW`$ = 179$`\pm `$17 eV ($`\chi ^2`$ = 309.5/292 $`dof`$).

A further intriguing possibility for the reflection component is the existence of multiple emission lines. In Figure 2, we have displayed the 5.7–7.1 keV pn data (6.0–7.5 keV rest-frame energies), to take advantage of the higher signal-to-noise data afforded by the EPIC pn instrument. On careful inspection of the residuals in the pn data we can distinguish two separate humps between 6–7 keV. We introduce a single Gaussian into the fit and fix its energy at 6.40 keV to account for emission from neutral Fe K$`\alpha `$. The resulting line width is $``$ 70 eV, and both the quality of the fit and the residuals are improved ($`\mathrm{\Delta }\chi ^2`$ = 20.2 for 2 additional free parameter). A second Gaussian is added to model the remaining residuals at slightly higher energies. This second Gaussian further improves the fit and minimises the residuals in this spectral region ($`\mathrm{\Delta }\chi ^2`$ = 6.4 for 3 additional free parameters). This second line has an energy of $`E`$ $``$ 6.87 keV and it is slightly broader ($`\sigma `$ $``$ 140 eV) than the neutral Fe K$`\alpha `$ line. The energy of this broad line is consistent with emission from a blend of H- and He-like iron. Formal modelling, with all parameters left free to vary, resulted in a similar fit quality as the step-wise process just described. ($`\chi ^2`$ = 307.4/290 $`dof`$).

We further examined the possibility that the positive residuals between 6–7 keV could be the result of three narrow Fe K$`\alpha `$ lines: neutral (6.40 keV), He-like (6.70 keV), H-like (6.97 keV) lines. We included three narrow ($`\sigma `$ = 10 eV) Gaussian profiles at the above stated energies and computed the fit. The triple line model is a good fit to the EPIC data ($`\chi ^2`$ = 308.6/293 $`dof`$). We obtain values of $`EW`$ = 32 $`\pm `$ 1eV, 52 $`\pm `$ 2 eV, and 52 $`\pm `$ 2 eV for the neutral, He-like, and H-like Fe K$`\alpha `$ lines, respectively.

Modelling of the $`ASCA`$ observation of I Zw 1 required a strong, (partially) ionised Fe K$`\alpha `$ line. Leighly (1999a) fitted a 1300 eV line with a rest-energy of 7.0$`{}_{0.3}{}^{}{}_{}{}^{+0.6}`$ keV, while Reeves & Turner (2000) established the line centre at 6.77$`{}_{0.17}{}^{}{}_{}{}^{+0.11}`$ keV with $`EW`$ = 483$`{}_{211}{}^{}{}_{}{}^{+212}`$ eV. Both line energies are consistent with our simultaneous fit of the pn and MOS data, using a single, broad Gaussian profile ($`E`$ = 6.69$`{}_{0.33}{}^{}{}_{}{}^{+0.14}`$ keV), though our line appears substantially weaker ($`EW`$ = 179 $`\pm `$ 17 eV). Utilizing the higher sensitivity of the EPIC pn detector, we find that the broad feature between 6–7 keV can be well described by emission from multiple narrower lines, showing that the evidence for a truly broad line is weak. Statistically there is little difference between the two and three Gaussian fits, both are acceptable; however, the residuals at about 6.7 keV ($``$6.3 keV in the observed frame) are slightly improved by the three narrow Gaussians model. A longer exposure of I Zw 1 with XMM-Newton, in order to boost the signal-to-noise of the MOS data, would clarify the situation in this spectral region.

### 4.2 The broad-band spectrum

As is illustrated by Figure 1 a power-law model is a reasonable starting point for the broad-band (0.3–10 keV) fit. With the exception of some complexities at energies lower than $``$ 1 keV, the spectrum is relatively smooth. For simplicity, we shall fit the iron complex (6–7 keV) with a double Gaussian model: the first Gaussian profile will be fixed at 6.4 keV to describe emission from neutral iron, and the second profile at slightly higher energies will be used to approximate emission from ionised species of iron. In addition, we have introduced a second cold absorption term. The first term is fixed to the Galactic value, while the second term will be allowed to vary and will be used to estimate the neutral absorption in I Zw 1.

With these additional parameters in place, the single power-law fit was improved at high energies, but the residuals below $``$ 1 keV were not affected ($`\chi ^2`$ = 1330.7/1071 $`dof`$). We added a second continuum component to the existing model to fit the weak soft excess at lower energies. Neither a single black body ($`\mathrm{\Delta }\chi ^2`$ = 7 for the addition of 2 free parameters), nor a second power-law ($`\mathrm{\Delta }\chi ^2`$ = 10 for the addition of 2 free parameters) were entirely effective in fitting the soft excess. A significant improvement in the overall fit is seen when the single power-law continuum is replaced with a broken power-law ($`\mathrm{\Delta }\chi ^2`$ = 67.4 for the addition of 2 free parameters). In this model the hard power-law ($`\mathrm{\Gamma }`$ $``$ 2.48) dominates the spectrum at energies higher than $``$0.75 keV. At lower energies the spectrum is fitted with a steeper power-law with a photon index of ($`\mathrm{\Gamma }`$ $``$2.86. An $`F`$-test suggests that the broken power-law fit is an improvement over the single power-law at $`>99.99`$% confidence, and that a second softer component is a required addition. An alternative model, which is thermal in nature, is a double blackbody fit ($`\mathrm{\Delta }\chi ^2`$ = 98.3 for the addition of 4 free parameters). In this case, the steep power-law ($`\mathrm{\Gamma }`$ $``$ 2.32) continues to dominate much of the spectrum; however, the cooler black body ($`kT_1`$$``$78 eV) contributes about 20% of the total unabsorbed flux, while the hotter black body ($`kT_2`$$``$213 eV) contributes about 5%.

Even though both the broken power-law and double blackbody models are dramatic improvements over the single power-law fit statistically, the residuals in the low-energy band are still unsatisfactory (Figure 3).

In both cases, the residuals in the observed frame between 0.4–0.6 keV still exhibit the up-down trend seen in the single power-law fit, while between 0.6–0.75 keV the data remain slightly overestimated. Unfortunately, the short exposure and low signal-to-noise in the RGS does not allow us to examine productively the soft spectrum of I Zw 1. A by-eye inspection of the RGS spectrum does indicate a number of possible narrow absorption and emission features, but formal modelling of the spectrum does not reveal any significant results.

One successful attempt to minimise the residuals in the soft spectrum is by the addition of a Gaussian emission line. Depending on the continuum model used, the line energy is between 535–565 eV, with a width of $`\sigma `$ = 40–55 eV and an equivalent width of 25–30 eV. The addition of the soft emission feature to the broken power-law ($`\mathrm{\Delta }\chi ^2`$ = 84 for the addition of 3 free parameters) or double blackbody model ($`\mathrm{\Delta }\chi ^2`$ = 65.6 for the addition of 3 free parameters) is significant. The measured energy is consistent with emission from the He$`\alpha `$ triplet of OVII at $``$22 Å (the resonance, intercombination, and forbidden lines), which are known to be strong in some Seyfert galaxies (e.g. NGC 3783, Kaspi et al. 2002; Behar et al. 2003). At least two “emission line-like” features are suggested by the RGS spectrum at this energy, giving credence to this model. As a word of caution, we would like to make note of the calibration issues around the O-edge in the EPIC pn detector. Although the calibration is significantly improved in SAS v5.4.1 (Kirsch 2003), and we do not believe that it has any bearing in this observation, it is difficult to make strong claims about this feature without independent data.

An alternative method to minimise the residuals in the soft spectrum is to add an absorption feature (edge or line) on the blue-side of the excess. I Zw 1 does possess UV absorption (Laor et al. 1997; Crenshaw et al. 1999), so it seems perfectly reasonable that it would also possesses some complex X-ray absorption.

In the first attempt we introduce an absorption line into the fit. The addition of the absorption feature results in a small statistical improvement to the broken power-law model ($`\mathrm{\Delta }\chi ^2`$ = 18.7 for the addition of 3 free parameters). The line energy is at $`E`$ $``$ 749 eV and it is very weak ($`EW`$ $``$ $``$6 eV). However, the addition of the absorption line is better received by the double blackbody model ($`\mathrm{\Delta }\chi ^2`$ = 49.9 for the addition of 3 free parameters). In this case, the equivalent width is about $``$22 eV and the line is centred at about 737 eV. The addition of the line also minimises the residuals in a satisfactory manner. Broad absorption features at these energies are most consistent with unresolved transition array (UTA) absorption due to Fe M-shell ions ($``$729–774 eV; Behar et al. 2001), as observed in IRAS 13349+2438 (Sako et al. 2001), NGC 3783 (Blustin et al. 2002; Kaspi et al. 2002; Behar et al. 2003), and Mrk 509 (Pounds et al. 2001).

Alternatively, the low-energy absorption line can be replaced with an absorption edge. Again, for the two continuum models we are working with, the statistics are improved: $`\mathrm{\Delta }\chi ^2`$ = 68.4 and 60 for the addition of 2 free parameters for the broken power-law fit and the double blackbody fit, respectively. For both models, the optical depth of the edge is $`\tau `$ $``$ 0.23–0.26, and the edge energy is at $`E`$ $``$ 664–689 eV. The best-fit edge energies do not agree with the strongest edge feature which is expected at 731 eV (OVII). However, this inconsistency in the edge energy alone does not permit us to discard this interpretation. The edge could very well be blended with other absorption features, perhaps even with some neutral iron absorption as was seen in MCG–6–30–15 (Lee et al. 2001), in which case our model is too simple.

The three best-fitting phenomenological models of the 0.3–10 keV spectrum of I Zw 1 are shown in Table 1. The continuum can be modelled by either a broken power-law or a double blackbody and power-law. There are primarily two difficulties in distinguishing between these apparently very different models. The first is that the soft spectrum is complex. We demonstrated that an absorption edge or an emission feature is a significant improvement to both continuum models, in addition to extra neutral absorption. The second fact is that the soft excess is weak, and the hard power-law component dominates much of the spectrum down to $``$750 eV.

All of the models have the disadvantage of possibly overestimating the photoelectric absorption in I Zw 1. The models require cold absorption on the order of 1021 cm-2 in addition to the Galactic value, which is nearly a factor of 10 higher than was estimated by the $`ROSAT`$ observations (Boller et al. 1996; Lawrence et al. 1997). Given that $`ROSAT`$ was sensitive down to 0.1 keV, it stands to reason that its value for the intrinsic absorption is reliable. However, the high value measured with XMM-Newton is most likely a proxy for complicated absorption occuring below 0.4 keV (e.g. blended absorption lines and/or edges). Given the very complex nature of the low-energy absorption in Seyferts, fitting simple models such as these can only be illustrative of a much more complicated situation. A long exposure, to utilise the high-resolution RGS, is required to probe the soft absorption in a meaningful way.

A slight advantage is given to the double blackbody plus power-law continuum in its ability to fit the iron emission between 6–7 keV. The model is successful in predicting a high-energy continuum and an iron complex which is similar to that described in Section 4.1. In contrast, the broken power-law continuum requires a very broad ($`\sigma `$$``$1.5 keV) neutral line. The neutral line width is certainly overestimated by the fit to compensate for a continuum which is poorly defined between 7–10 keV. A possible shortcoming for the double blackbody model is that the temperatures estimated for the thermal components are high ($`kT_1`$$``$198 eV and $`kT_2`$$``$94 eV) for a standard accretion disc around a 107 M black hole. However, since we do not understand the physics of the soft excess very well we cannot automatically deem the model as unphysical. For example, bulk motions in the disc (Socrates et al. 2003) is just one way in which the shape of the soft excess can be altered from that expected from a standard accretion disc. If we allow a soft excess component with curvature (rather than a power-law) the systematic residuals at low energies can be made to diminish. In addition, such curvature in the soft excess allows for better fitting of the high-energy continuum component and hence, the iron emission features. The final model residuals are shown in Figure 4.

The average 0.3–10 keV flux, estimated from the double blackbody fit, and correcting for Galactic extinction is $``$2.40 $`\times `$ 10-11 erg s-1 cm-2. The corresponding luminosity is 2.03 $`\times `$ 1044 erg s-1, and the luminosity in the 2–10 keV band is 0.72 $`\times `$ 1044 erg s-1. In comparison, during the $`ASCA`$ observation the 0.5–10 keV luminosity was 1.36 $`\times `$ 1044 erg s-1, and the 2–10 keV luminosity was 0.53 $`\times `$ 1044 erg s-1 (Leighly 1999).

## 5 Timing analysis

In this section, we will discuss the variability properties of I Zw 1. The pn and MOS light curves were both analyzed and found to be entirely consistent. We will only discuss the results from the pn, given its higher signal-to-noise. The 1$`\sigma `$ errors on the light curves were calcluated using counting statistics.

### 5.1 The broad-band light curve

In Figure 5 we present the 0.3–12 keV EPIC pn light curve of I Zw 1 in 200 s bins.

The average count rate is 6.97 $`\pm `$ 0.19 counts s-1. The minimum and maximum 0.3–12 keV count rates are 5.76 $`\pm `$ 0.17 and 9.80 $`\pm `$ 0.22 counts s-1, respectively. The most notable feature in the light curve is obviously the flare starting at $``$ 10 ks, when the count rate increases by $``$40% in 2500 s. The flare is remarkably symmetric as the count rate returns from maximum to its pre-flare value in the same amount of time. Except for the flare, the rest of the light curve is surprisingly serene. Considering only the periods between 0–10 ks and 15–20 ks, the fluctuations about the mean are only about 10%, which is comparable to the observed variability on similar time scales during the $`ASCA`$ observations (Leighly 1999b); but modest compared to many NLS1 X-ray light curves. However, an observation with $`ROSAT`$ did find that the 0.1–2.4 keV flux dropped by about 25% in $``$6000 s (Boller et al. 1996). It was not clear from that observation if the flux drop was subsequent to a rapid rise or part of a gradual decrease for an extended period of time. In the remainder of this paper we will often make reference to I Zw 1 in its high- and low-flux states. In defining the high-flux state, we have selected the entire flaring event from $``$ 10–15 ks. The low-flux state includes all of the data from the remaining time intervals.

We have calculated the radiative efficiency ($`\eta `$), assuming photon diffusion through a spherical mass of accreting material (Fabian 1979; but see also Brandt et al. 1999 for a discussion of important caveats). The most rapid rise occurs from 10.1–12.6 ks when the luminosity increases by 7.41 $`\times `$ 1043 ergs s-1. The corresponding radiative efficiency is only $`\eta `$ $``$ 1.5%, which does not necessarily require radiative boosting or anisotropic emission during the flaring event.

As a first test for spectral variability we calculate the fractional variability amplitude (Fvar; Edelson et al. 2002) to examine the degree of variability in different energy bands. The fractional variability amplitude is calculated in eight energy bins between 0.3–12 keV using 200 s binning of the light curves (Figure 6).

The Fvar spectrum of I Zw 1 is remarkably different from the typically flat spectrum seen in other NLS1 such as 1H 0707–495 (Boller et al. 2002) and Ton S180 (Vaughan et al. 2002). The variations intensify as the energy increases and maximize between 2–12 keV. The Fvar spectrum gives the sense of a variable hard power-law which is physically distinct from the softer, less variable, emission component.

Much of the spectral variability observed in Figure 6 is associated with the X-ray flare seen in the broad-band light curve. We recalculated the Fvar, but this time only for the low-flux state. The result is shown as the red squares in Figure 6. The low-flux state fractional variability in each bin is less than 10%, and consistent with a constant across the entire spectrum ($`\chi ^2`$ = 2.6/6 $`dof`$). Figure 6 demonstrates that the spectral variability is stimulated by the flux outburst in Figure 5.

### 5.2 The hard X-ray flare and flux-induced spectral variability

By investigating the light curves in various energy bands we can demonstrate that the flare is concentrated at higher energies. In Figure 7 we have plotted the normalised light curves in the 0.3–0.8 keV and the 3–12 keV bands. The light curves are normalised to the average count rate during the periods 0–8.2 ks and 17–20 ks. Clearly, the higher energy band is significantly more variable during the flaring event than the soft band. During the flare the hard count rate increases by $``$100%, whereas the soft band shows only about a 20% increase. This is the first time, to the best of our knowledge, that a hard X-ray flare has been detected from a NLS1.

The hardness ratio variability curve (Figure 8) is calculated using the formula $`(HS)/(H+S)`$, where $`H`$ and $`S`$ are the count rates in the hard and soft bands, respectively.

The hardness ratio variability curve shows similar behavior as the light curve in Figure 5. The hardness ratio exhibits only small variations for most of the observation, but then the fluctuations become enhanced between 10–15 ks, showing spectral hardening precisely during the time of the hard X-ray flare. The significance of the relation between Figure 5 and Figure 8 is measured with the Spearman rank correlation coefficient to be significant at $`>`$ 5$`\sigma `$.

The fact that the flare is stronger at high energies implies that its origin is not in the cold accretion disc. It is plausible that the flare originates in the accretion-disc corona due to magnetic reconnection (e.g. Galeev et al. 1979). The hard flare would then irradiate the disc, where the high-energy photons will be reprocessed to lower energies via Compton down-scattering. In this simple picture, we may expect to see a time lag between the low and high energy light curves if the light travel time between the corona and disc is substantial. In order to investigate this we calculate cross correlation functions using 200 s bins for a number of light curves in various energy ranges. All of the cross correlations are relatively symmetric and consistent with zero time delay. The results suggest that there are no leads or lags greater than $`\pm `$100 s between the various energy bins. Assuming the simple flare geometry discussed above, this result implies a very compact system, in which the light travel time between the disc corona and the accretion disc itself is $`<`$ 100 s. All responses to the hard X-ray flare (light travel-time, reprocessing, spectral changes) appear to occur on short time scales.

Another consideration is whether the flare, which appears to induce the spectral variability, has any long-term effect on the spectrum; that is, even after the flux has returned to its pre-flare value. To examine this possibility in a model-independent manner, we calculated the ratio between the pre-flare data ($`<`$ 10 ks) and the post-flare data ($`>`$ 15 ks). The ratio is displayed in Figure 9.

There appears to be little difference in the ratio across the spectrum, and a constant fit to the ratio is acceptable ($`\chi ^2`$ = 29.5/30 $`dof`$). Even the difference in the average flux is on the order of 2%. However, on closer examination we notice that all nine data points below $``$1.2 keV lie above the ratio average. The probability that this is a random event is 1/512. Figure 9 indicates that the soft flux prior to the flare is higher than the post-flare soft flux. The shortage of post-flare data makes it difficult to determine whether this subtle spectral difference is associated with the flare, or a coincidental event.

### 5.3 The nature of the spectral variability

We have demonstrated that the flare is indeed hard, and that it induces spectral variability. The spectrum appears to become harder as the flare intensifies, and the changes brought on by the flare are immediate. Could the spectral variability be described by a component with a constant spectral slope but variable flux? Or is it more likely that the spectral variability is a result of a pivot in the X-ray spectrum at some higher energy? Taylor et al. (2003) have presented a model-independent technique to distinguish between these forms of spectral variability. However, given the limited data we have available from this short observation it is not possible to apply their technique successfully.

Fabian et al. (2002) and Fabian & Vaughan (2003) have successfully used a difference spectrum (high-flux minus low-flux spectrum) to investigate the nature of the spectral variability in MCG–6–30–15. If the difference spectrum can be fitted with a power-law then the spectral variability between the two flux states can be attributed to a difference in the normalisation of the intrinsic power-law without requiring a change in its spectral slope. In Figure 10 we display the residuals from the difference spectrum for I Zw 1 fitted with a power-law modified by cold absorption. The fit is a reasonable approximation ($`\chi _\nu ^2`$ = 1.05).

We further considered the spectral differences between the two flux states by modelling the high and low-flux state pn spectra separately. To both spectra we re-fit the double blackbody model described in Section 4.2 (Table 1 model (a)). In both the low- and high-flux states, we found that the only significant changes in the model parameters were to the normalisation of the continuum components. Within 90% confidence the values for the temperatures of the thermal components, edge energy, and neutral absorption were consistent. Even the photon indices of the low- and high-flux state power-law components were consistent ($`\mathrm{\Gamma }_{high}`$ = 2.31$`{}_{0.04}{}^{}{}_{}{}^{+0.06}`$ and $`\mathrm{\Gamma }_{low}`$ = 2.35$`{}_{0.06}{}^{}{}_{}{}^{+0.10}`$). A noticeable difference occurred only in the relative fluxes ($`f`$ = F3-10keV / F0.3-2.0keV) during the two flux states. During the low-flux state $`f_{low}`$ = 0.20 $`\pm `$ 0.01, whereas during the high-flux state $`f_{high}`$ = 0.26 $`\pm `$ 0.02.

A further point to note is the possibility of changes in the strength of the iron line(s) between the different flux states. The spectra during the two flux states show a very similar continuum shape across the 0.3–10 keV band; however, the iron lines appear stronger during the low-flux (non-flare) state while the normalisation of the lines remains consistent. The equivalent widths of the neutral and ionised iron lines during the low state are 65$`{}_{6}{}^{}{}_{}{}^{+4}`$ eV and 143$`{}_{14}{}^{}{}_{}{}^{+10}`$ eV, respectively, while during the flare the strengths are 46$`{}_{4}{}^{}{}_{}{}^{+2}`$ eV and 70$`{}_{10}{}^{}{}_{}{}^{+4}`$ eV. Although the line parameters are especially model dependent and our exposure is short, I Zw 1 appears to be another object in which the iron lines do not vary as significantly or in response to (in a simple manner) the continuum flux variations (Fabian et al. 2002; Lee et al. 2000, Reynolds 2000, Iwasawa et al. 1996).

## 6 Summary

We have presented a 20 ks XMM-Newton observation of I Zw 1, the prototypical NLS1. Our findings are summarized below.
(1) The 0.3–10 keV continuum can be described by either a double blackbody plus a power-law or a broken power-law, modified by intrinsic cold absorption. The soft excess is weak, but notable at energies between $``$0.4–0.7 keV. Superimposed on this continuum are absorption and/or emission features below $``$ 1.0 keV, and strong Fe K$`\alpha `$ emission between 6.4–6.8 keV.
(2) The Fe K$`\alpha `$ emission can be successfully fitted with several models, and the current observation does not allow us to determine the best one. The pn data suggest that the Fe K$`\alpha `$ complex may be resolved into emission from a neutral line at 6.4 keV and an ionised component at $``$6.9 keV (blend of He- and H-like Fe), weakening the case for a truly broad line.
(3) The low-energy spectrum can be modelled with either an absorption line or edge. The current observation does not allow us to distinguish between these two models. On the other hand, a low-energy emission feature ($`E`$ $``$ 0.55 keV) can be introduced to the fit with relative success.
(4) We find an X-ray flare which is isolated mostly in the 3–12 keV band. We find that the hard X-ray flare induces spectral variability in such a way that the spectrum becomes $`harder`$ during the flare. There are no discernible lags between energy bands down to $`\pm `$100 s. Comparison of the pre-flare and post-flare spectra show that there is a subtle difference in the spectrum prior to and after the flare. It is not certain whether this spectral variability is related to the hard flare.
(5) The flux of the iron line(s) do not appear to respond to the change in the continuum flux during the flare. The lines are strongest during the low-flux (non-flare) state.

### 6.1 The low-energy spectra

The quality of the data from this short exposure does not allow us to distinguish between an absorption line or an edge, or an emission line in describing the features in the soft spectrum. If the absorption is interpreted as a line then the best-fit energy may be suggestive of an UTA feature due to Fe M absorption. An absorption edge is a slightly better fit than the line. The best-fit edge energy is not consistent with that of OVII, which would indicate that the edge is severely blended with other absorption features. The emission line interpretation is consistent with emission coming from a blend of OVII lines at $``$ 22 Å. The most likely scenario is that we have a combination of OVII emission as well as absorption, but with the resolution of the EPIC instruments it is not possible to measure accurately the contribution from each component. In any regard, we have demonstrated that the soft spectrum of I Zw 1 is complex and warrants a deep observation to utilise the RGS.

### 6.2 The hard X-ray flare

The detection of a hard X-ray flare, and spectral hardening during the flare, is entirely consistent with theories of an accretion-disc corona which is being heated by magnetic reconnection (e.g. Galeev et al. 1979). The apparent instantaneous response of the lower-energy emission mechanisms to the hard X-ray flare suggests that the light travel-time between the emission regions is $`<`$ 100 s. Deeper analysis of the spectral variability suggests that spectral changes are due to an increase in the high-energy flux relative to the low-energy flux, and not due to a change in the intrinsic shape of the continuum. The weakness of the flare at low X-ray energies may be indicating that the seed photons for the Comptonisation component are not coming from the 0.3–1 keV band, but probably from lower UV energies. This observation could perhaps be useful in constraining Comptonisation models.

This short exposure with XMM-Newton has revealed a number of fascinating traits in I Zw 1. I Zw 1 has once again proven to be an important object in allowing us to probe many interesting facets of AGN, such as: X-ray flares, warm-absorption, and complex Fe K$`\alpha `$ emission, A long exposure with XMM-Newton and Chandra would be quite rewarding for studying these traits in greater detail.

###### Acknowledgements.

The authors are thankful to Frank Haberl for help regarding the EPIC calibration issues, and the anonymous referee for helpful comments regarding the soft emission feature. Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). WNB thanks NASA grants NAG5–9924 and NAG5–9933.


text:
# A Tight Upper Limit on Oscillations in the Ap star ϵ Ursae Majoris from WIRE Photometry

## 1 Introduction

The excitation of pulsations in the lower part of the classical instability strip is still far from being understood. Some normal stars in this category show $`\delta `$ Scuti pulsations while others do not, and those that do pulsate only do so in a seemingly random subset of possible modes. Among the chemically peculiar stars, the so-called Ap stars that are the subject of this Letter, there is a similar puzzling division between the rapidly oscillating Ap stars (roAp stars) and the non-oscillating variety (noAp stars). This has lead some to propose that oscillations are present in noAp stars but at an undetectable level. Here we report one of the most precise photometric observations to date on an Ap star, which provides a strong upper limit on oscillations.

The Ap stars have peculiar spectra, with lines from elements such as Silicon, Magnesium, Mercury, Chromium, Europium and Strontium. The roAp stars have pulsations with periods about 5–15 min and typical amplitudes of a few milli-mag (Kurtz, 1982, 1990; Kurtz & Martinez, 2000). So far, thirty two objects have been classified in this subclass (Kurtz & Martinez 2000; Kurtz, personal communication). For recent theoretical works on the oblique pulsator model and the excitation mechanism in roAp stars see Bigot et al. (2000), Cunha & Gough (2000), Balmforth et al. (2001) and Bigot & Dziembowski (2002).

A detailed study of the differences between roAp and noAp stars was carried out by Hubrig et al. (2000). They found that as a group, the roAp stars are 0.57$`\pm `$0.34 mag brighter than the zero-age main sequence (ZAMS) while noAp stars are 1.20$`\pm `$0.65 mag above the ZAMS, which suggests that the latter are slightly more evolved than the roAp stars. It has also been proposed that an overabundance of the rare earth elements Nd and Pr may be a signature of roAp stars (Cowley & Bord, 1998; Weiss et al., 2000; Ryabchikova et al., 2001, 2002). Kochukhov et al. (2002) suggested, on the other hand, that all Ap stars in a certain temperature range pulsate and that the observed difference between roAp and noAp stars comes from the amplitude of the pulsation (i.e., it is simply below detectability in noAp stars).

A key point in testing this last suggestion is to press down the observational upper limits on the pulsation amplitudes in noAp stars. Kurtz et al. (2003a) observed the Ap star HD 965 over 2.7 h using a 2-m telescope and put an upper limit of 0.2 milli-mag or $``$200 parts per million (ppm) on the amplitude of oscillations in the $`B`$-band. This is the lowest limit for pulsation in a noAp star in the literature. However, the rotation period of HD 965 is uncertain and believed to be longer than two years. Kurtz et al. (2003a) thus concluded that, although it is likely that HD 965 is a noAp star, their observations cannot reject the possibility that it is nevertheless a roAp star. The reason is that according to the oblique rotator model, the amplitude of the oscillation is higher near maximum of the rotational modulation and may become zero at other phases (Kurtz, 1990). Therefore, Kurtz et al. (2003a) stated that further observations of this star over a few years are required to substantiate its status as a noAp star. Here we report observations of the Ap star $`ϵ`$ UMa, whose rotation period is short enough to avoid this problem.

$`ϵ`$ UMa (HR 4905 = HD 112185) is classified as an A0p Cr star (Bohlender & Landstreet, 1990). Its effective temperature and radius were estimated as $``$9000K and $`4R_{}`$ (Pyper & Adelman, 1985; Shallis & Blackwell, 1979). At V=1.8 it is the brightest Ap star in the sky, but it is not known to pulsate.

A period of 5.0887 d was detected by Guthnick (1931) in absorption line intensities of $`ϵ`$ UMa and was confirmed by further photometric and spectroscopic studies (Struve & Hiltner, 1943; Provin, 1953; Woszczyk & Jasinski, 1980; Woszczyk & Michałowski, 1981). Provin (1953) found that the variation has a double wave structure with the two maxima being separated by half the period. The periodicity is understood as the rotation period of the star, but the strength of the magnetic field is still uncertain (Rice & Wehlau, 1990; Donati et al., 1990; Bohlender & Landstreet, 1990).

## 2 Observations and Reduction

After the failure of the main mission of the Wide Field Infrared Explorer (WIRE) satellite, launched by NASA in March 1999, its star tracker was successfully used for photometry of bright stars (Buzasi et al., 2000; Buzasi, 2002; Poretti et al., 2002; Cuypers et al., 2002; Retter et al., 2003).

The detector is a 512 $`\times `$512 SITe CCD with 27 $`\mu `$m pixels; each pixel corresponds to $``$1’ on the sky. The camera is unfiltered and the effective response is roughly $`V+R`$. The observations of $`ϵ`$ UMa were obtained in the period 2000 June 22 – July 21, with a cadence of 0.5 seconds. Data reduction was done by applying a simple aperture photometry algorithm, which involved summing the central 4 $`\times `$ 4 pixel region of the window of 8 $`\times `$ 8 pixels that can be read. The background level, due mainly to scattered light from the bright Earth, was estimated from the four corner pixels of each image, and subtracted. Following aperture photometry, any points deviating more than $`2.5\sigma `$ in mean flux, image centroid, or background level were rejected. The data points were then phased onto the satellite orbit, extreme points were rejected and the mean shape was subtracted. A final total of 5.7 million data points (out of about 6.1 million initial measurements) were obtained after applying these clipping procedures. Further details on the reduction method can be found in Retter et al. (2003).

## 3 Analysis

The upper panel of Fig. 1 presents the smoothed light curve of $`ϵ`$ UMa, where data from each satellite orbit have been binned into a single mean value, resulting in 308 points. Around JD 2451732 there are jumps of $``$10% in the brightness of $`ϵ`$ UMa. A comparison with the background level and the location of the centroid of the star on the CCD at this time (Fig. 1, three lower panels) shows that the effect is instrumental. The light curve is also affected by this effect, but to lesser extent, around JD 2451720 and 2451738.

In an attempt to remove from the light curve the effect of minor shifts in the location of the star in the CCD and small changes in the background level, several decorrelation methods were tried (Brown et al., 1991; Robinson et al., 1995). However, unlike for the case of our observations of Arcturus (Retter et al., 2003), these techniques did not satisfactorily correct the light curve. The reason for this behavior might be that, unlike Arcturus (which stayed within $`\pm `$0.03 pixel during the observing run), $`ϵ`$ UMa moved several tenths of a pixel from orbit to orbit. Consequently, the linear or low-order polynomial approximation of sub-pixel sensitivity variations was not valid over the range occupied by the star.

Fig. 2 displays the smoothed light curve of $`ϵ`$ UMa during the last two weeks of the observations, when the instrumental effects were minimal. The light curve clearly shows the presence of the 5.1-d rotation period (Section 1), which has a double-hump shape.

To search for oscillations (presumably around 10 min) in the data from the last two weeks of the run, the slow variation from rotation was first subtracted. The upper panels of Figs. 3 & 4 show the amplitude spectrum of this high-pass-filtered time series. The graph is dominated by strong peaks that correspond to the satellite orbit, its harmonics and, at lower amplitudes, their 1-d-1 aliases. The bottom panel in Fig. 4 shows the spectral window. The 1-d-1 aliases are much stronger in the data (relative to the satellite orbit and its harmonics) compared with the spectral window. This behavior indicates that the data show instrumental variations within orbits and also with a periodicity of a day, perhaps from variable reflection from Earth on the satellite. The second panels in the figures display the amplitude spectra after fitting and subtracting all these frequencies simultaneously.

All peaks in the residual amplitude spectrum (which are not harmonics or 1-d-1 aliases) have amplitudes below 75 ppm. Most of the highest remaining peaks are, however, close to the harmonics of the satellite orbit and all appear in the spectral window. We conclude that there are no oscillations with amplitudes above 75 ppm.

The possibility cannot be ruled out that the star oscillates with a single mode that falls exactly on the satellite orbit, its harmonics or their 1-d-1 aliases. We estimated the probability of this occurring by chance by dividing the typical width of the peaks (taken as the full-width at half maximum, which is set by the length of the observing run) by 16 d-1 (the satellite orbit, for higher amplitudes) and by 1 d-1 (for lower amplitudes, since the data suffer from 1-d-1 aliases). The probabilities are about 0.5% and 7%, respectively. In other words there is a 0.5% chance that it oscillates with an amplitude below $``$1800 ppm and a 7% chance for an amplitude below $``$250 ppm.

To check whether the star has stronger pulsations at a specific phase of the 5.0887-d rotation period (see Section 1), the data from the last two weeks of the run (after removing the variation of the rotation period) were divided into ten equal phase bins. The corresponding power spectra of the bins were consistent with the power spectrum of the whole data (Figs. 3 & 4) and no significant peaks were found (besides the satellite orbital frequency, its harmonics and their 1-d-1 aliases).

## 4 Discussion

The WIRE photometry clearly shows that $`ϵ`$ UMa is variable on time scales of a few days with a period of $``$5 days (Fig. 2). Our results confirm the previous detections of the rotation period (Section 1). The noise in our data is comparable to that in the highest precision photometry done so far on roAp stars obtained during several years (Kurtz et al., 1997) or by the WET collaboration (Kurtz et al., 2003b).

Of major interest is the question whether $`ϵ`$ UMa pulsates or not. We found an upper limit of 75 ppm on oscillations in this star, unless by chance it has a single mode that falls on the satellite orbit, its harmonics or their 1-d-1 aliases. This, however, is very unlikely since roAp stars usually show several pulsation periods in their light curves (Kurtz, 1990).

We note that the observations were done using an effective response of roughly $`V+R`$ (Section 2). The amplitude of pulsations in roAp stars decreases with wavelength (Kurtz & Medupe, 1996; Medupe & Kurtz, 1998; Medupe, 2000). $`ϵ`$ UMa may, therefore, pulsate with amplitudes larger than the limit given above in narrower and bluer bands. Using the detector response, a typical spectrum of an Ap star (Le Borgne et al., 2003) and Fig. 1 of Medupe & Kurtz (1998) we estimated that the limits we can put on the amplitude of pulsations in $`ϵ`$ UMa are about 50% and 120% larger in $`V`$ and $`B`$, respectively. Kurtz & Medupe (1996), Medupe & Kurtz (1998) and Medupe (2000) found that the phase of the modulation also depends on the wavelength. This effect is, however, relatively small and the phase shift between the optical filters is about 0.5 radians. Therefore, different colors are probably not anti-phased with each other and cannot cancel the effect of pulsations. Using simulations, we found that in our case the reduction of the amplitude by this effect is negligible (less than 3%).

The rotation period of $`ϵ`$ UMa is relatively short at 5.1 d (Sections 1 & 3) and the WIRE observations covered several cycles of this periodicity. Therefore, it cannot be claimed that the non-detection of oscillations presented in this work is because of observations at a wrong rotation phase, which may be the case for noAp stars with long rotation periods (Section 1). In fact, we could not find oscillations at specific phases (Section 3). Our results thus suggest that $`ϵ`$ UMa does not oscillate at all. It is interesting to note that it was proposed that an anomaly in the abundances of Nd and Pr may be a signature of roAp stars (Section 1). $`ϵ`$ UMa does not have this anomaly and the null result on oscillations in this star thus supports this suggestion.

It was suggested that all Ap stars in a certain temperature range pulsate (Kochukhov et al., 2002). The temperature of $`ϵ`$ UMa was estimated as $``$9000K (Section 1), which is slightly higher than the range of roAp stars – 7000–8300K (Kurtz, 1990), although these limits are somewhat uncertain. It is therefore possible that the non detection of oscillations in $`ϵ`$ UMa is simply because the star is too warm for the instability strip for roAp stars.

The question of whether noAp stars do pulsate with amplitudes below the detection limit is still open. Future space missions that significantly decrease the upper limit on the amplitudes of pulsations in these stars may give an answer to this question.

We are grateful to all the people at NASA headquarters, GSFC and IPAC who made the use of the WIRE satellite possible. We thank the anonymous referee for useful comments and Don Kurtz for helpful discussions. AR, TRB and LLK are supported by the Australian Research Council and HK by the Danish Natural Science Research Council, and the Danish National Research Foundation through its establishment of the Theoretical Astrophysics Center. DLB acknowledges support from NASA (NAG5-9318).