text:
# Enhancement of the Curie temperature in GaMnAs/InGaMnAs superlattices

## Abstract

We report on an enhancement of the Curie temperature in GaMnAs/InGaMnAs superlattices grown by low-temperature molecular beam epitaxy, which is due to thin InGaMnAs or InGaAs films embedded into the GaMnAs layers. The pronounced increase of the Curie temperature is strongly correlated to the In concentration in the embedded layers. Curie temperatures up to 110 K are observed in such structures compared to 60 K in GaMnAs single layers grown under the same conditions. A further increase in TC up to 130 K can be achieved using post-growth annealing at temperatures near the growth temperature. Pronounced thickness fringes in the high resolution X-ray diffraction spectra indicate good crystalline quality and sharp interfaces in the structures.

Developing ferromagnetic semiconductors with Curie temperatures near room temperature or beyond is of major interest for future spintronic devices. One material under intense investigation in recent years is the diluted magnetic semiconductor Ga1-xMnxAs which can be made ferromagnetic at present with Curie temperatures up to 110 KMatsukuratransport  in as-grown samples and up to 150 KKuannealing  in annealed samples. GaMnAs is fabricated by low-temperature (LT) molecular beam epitaxy (MBE) with growth temperatures of around 250C, and therefore, a lot of defects are incorporated into the lattice. Different defect species like As antisites, AsGa, or As interstitials, AsI, as well as Mn interstitials, MnI, are discussed in the literature,Potashnik -Blinowski\_interstitials  having one aspect in common: They are all acting as donors and are therefore compensating the Mn acceptors in GaMnAs. Hence, the defect structure influences strongly the electrical and, consequently, the magnetic properties of this material system, since the ferromagnetic coupling in GaMnAs is suggested to be due to an indirect exchange interaction between the magnetic moments of Mn++-ions mediated by holes.DietlScience287zenermodell  So, one possible way to increase the Curie temperature TC in GaMnAs is to increase the hole concentration, e.g., by reducing the defect density in the material by post-growth annealing at temperatures near or even below the growth temperature.Kuannealing ,Foxon  Another possibility to change the magnetic properties is to control the defect structure of the material even during growth.Tanaka\_epitaxial -Campion\_epitaxial

Here, we report on the increase of TC of GaMnAs due to strained (InyGa1-y)1-xMnxAs or InyGa1-yAs layers embedded into the magnetically active GaMnAs layers.

All samples studied were grown on epiready semi-insulating GaAs(001) substrates in a Riber 32 MBE system. To provide As4, an As cracker cell was used in the non-cracking mode. To supply Ga and In fluxes, conventional Knudsen cells were used, whereas a Hot-Lip cell was used for Mn. The growth procedure of the GaMnAs/InGaMnAs superlattices was as follows: First, a 100 nm thick GaAs buffer layer was grown at a substrate temperature of TS = 585C (standard GaAs growth conditions). Then the sample was cooled down during a growth break to the growth temperature of TS = 230C. The growth was continued by a 20 nm thick Ga1-xMnxAs layer followed by a 3-6 nm thick (InyGa1-y)1-xMnxAs layer. This LT procedure was repeated 15 times. Finally, the structure was capped by an additional 20 nm thick Ga1-xMnxAs layer. The In content in the InGaMnAs films was varied from sample to sample between zero and 50%. Mn-flux was kept constant during the whole LT-growth and corresponded to a Mn content of $`x`$ = 0.058 in the Ga1-xMnxAs layers. In the (InyGa1-y)1-xMnxAs layers it was somewhat less, depending on the In content in the sample. The same growth procedure was used to grow the GaMnAs/InyGa1-yAs superlattices except the missing Mn in the embedded layers. The growth was monitored in situ by reflection high energy electron diffraction (RHEED). The RHEED pattern showed (1 $`\times `$ 2) reconstruction during the growth of the whole LT part of the structure and no evidence of a second phase on the surface (i.e. MnAs-clusters) was observed. Details of the growth of GaMnAs single layers is also described elsewhere.Koeder\_gradient  Annealing experiments were performed in air using a LINKAM THMS 600 heating chamber. The annealing temperature was chosen to be 250C, and the annealing time was 30 min. Details of the annealing experiments are described elsewhere.wolfgang\_annealing

For all superlattices we performed high resolution X-ray diffraction (HRXRD) measurements. As seen for example in Fig. 1, the $`\theta `$-$`2\theta `$-scans indicate good crystalline quality and sharp interfaces between the GaMnAs and the InGaMnAs or InGaAs layers. Up to ten superlattice peaks can clearly be resolved, and even between the peaks pronounced thickness fringes are observed (see inset Fig. 1). From the superlattice fringes, the period of the superlattices can be caluclated. In the particular example, the 23.6 nm derived from HRXRD are in good agreement to what we expect from the growth rates. Finally, the sharp symmetrical peaks indicate a fully strained structure.

A series of samples containing GaMnAs/InGaMnAs superlattices were grown with varying In content from 0% In (single layer) up to 50% In. To extract the Curie temperature of the samples superconducting quantum interference device (SQUID) magnetization measurements were performed as a function of temperature (accuracy $`\pm `$ 5 K) in a small in-plane applied magnetic field (50 Oe $``$ (110)). As can be seen in Fig 2, TC is continously increasing with increasing In content (filled squares) from 60 K in the GaMnAs single layer up to 110 K in the as-grown structure containing 50% In.

In order to further analyze this effect, several reference samples were grown under the same growth conditions: First, two LT-GaAs/InGaMnAs superlattices containing 20 and 30% In, in order to clarify whether the higher Curie temperature stems from the GaMnAs or the InGaMnAs layers. These two samples, in Fig. 2 marked by filled triangles, show much lower Curie temperatures of around 30K. Therefore, the InGaMnAs films alone cannot be responsible for the high Curie temperatures in the superlattices. Second, two GaMnAs/InGaAs superlattices containing 40 and 50% In (open circles in Fig. 2). These two samples show enhanced Curie temperatures of 70 and 80 K, respectively, even though there is no Mn in the InGaAs. The reason for the differences in TC, measured for the superlattices with and without Mn in the InGaMnAs layer, is not yet known.

Obviously, the high Curie temperatures of the superlattices are due to an enhancement of TC in the GaMnAs layers by incorporating strained In-containing layers.

A further increase in TC was observed upon annealing the samples 30 min at 250C, which is near the growth temperature. Finally, a Curie temperature of 100 K was observed for the annealed GaMnAs single layer and 130 K for the annealed superlattice containing 40% In (open squares in Fig. 2).

The hole concentrations of heavily doped GaMnAs layers can be estimated by a careful inspection of the coupled plasmon-LO-phonon modes, observed in Raman spectroscopy. Details of this method are described elsewhere.LimmerRaman  In Fig. 3, the Raman spectra recorded from the GaMnAs single layer (sample B313) and from the GaMnAs/InGaMnAs superlattice with 40% In (sample B315) before and after annealing are depicted for comparison. All four spectra exhibit a broad Raman line, located near the frequency of the GaAs TO phonon, which arises from the coupled mode. With increasing hole concentration this line is known to shift from the frequency of the LO phonon to that of the TO phonon. Therefore, the Raman spectra reveal two features: First, an increase of the hole concentration upon annealing in both samples, and second, a significant enhancement of the hole concentration in the GaMnAs region of the superlattices compared to the GaMnAs single layer.

This agrees very well with the observed increase in TC. Figure 4 shows the SQUID magnetization curves of the superlattice in Fig. 3 (sample B315) before and after annealing. The increase in Curie temperature and also in saturation magnetization due to annealing can clearly be seen.

At present, the underlying mechanism of the drastic enhancement of TC in the GaMnAs/InGaMnAs heterostrucures is still unclear. However, we suggest that the incorporation of strained InGaMnAs or InGaAs films into GaMnAs leads to a reduction of the density of compensating defects, resulting in a higher hole density, and thus, in a higher Curie temperature than in GaMnAs single layers. This assumption is in agreement with the Raman data which revealed an increase of the hole density in the GaMnAs barrier layers of the superlattices compared to the GaMnAs single layer. From recent publicationsYuinterstitials ,Blinowski\_interstitials  it is known, that the defects, which are most likely responsible for carrier compensation in GaMnAs, are the Mn atoms on interstitial sites, and LT-annealing can reduce their density. Edmonds et al.Edmonds\_annealing\_interstitials  demonstrated that out-diffusion of MnI towards the surface takes place during annealing of GaMnAs single layers. In GaMnAs/InGaMnAs superlattices, the In-containing layers and/or interfaces can act as sinks for point defects like MnI and AsI during the growth. Most likely diffusion of MnI towards the heterointerfaces is responsible for the increase in carrier density and, consequently, TC. Since in our superlattices a further enhancement of TC can be achieved by post-growth annealing, too, a considerable amount of MnI ions seems to be still present in the as-grown superlattices. In contrast to GaAs/GaMnAs interfaces, which have been shown to prevent removal of MnI during annealing,Kuannealing  the embedded InGaMnAs layers obviously do not have any adverse effect on the annealing process.

In summary we have shown that the Curie temperature of epitaxially grown GaMnAs layers can be increased by embedding strained InGaMnAs or InGaAs layers into the magnetically active GaMnAs. We suggest that this is due to the influence of strain on the diffusion of defects like Mn and As interstitials into the In containing embedded layers, or to the reduction of their concentration during growth as compared to GaMnAs single films. The Curie temperature can be controlled by varying the In content in the embedded layers.

The authors acknowledge the financial support by the Deutsche Forschungsgemeinschaft, DFG Wa 860/4-1.


format: text

text:
# A new adiabatic quantum search algorithm

## 1 Introduction

Many computational search problems are surprisingly difficult to solve. Quantum computing is a promising candidate to tackle such difficult tasks. An important landmark was achieved by Grover: while an oracle search for a marked item out of $`N`$ classical unstructured items requires, on average, $`O(N)`$ steps, Grover  found a quantum algorithm that performs the oracle search in $`\sqrt{N}`$ steps. This quadratic speedup leads to the question of whether similar speedups can be achieved for other search problems.

Structured searches are natural extensions of the oracle search. They are used when the databases possess some structure. Exploiting the structure of the database will increase the performance of a search (classical or quantum). Cerf et al  were the first to study structured quantum search algorithms. Their work is based on quantum circuits. Recently Roland and Cerf  gave a quasi-adiabatic quantum version of the structured search.

To understand better the algorithm presented in  we recall that Grover’s algorithm uses an iterative improvement strategy. It starts with an equal superposition of all possible states representing items of the database, and by iterative application of unitary operations (gates) the initial state is rotated towards a state that encodes the solution of the search problem. However if the database has some structure of its own other strategies can be used. For example, one can use the divide-and-conquer strategy. In this strategy one divides the problem into subproblems of manageable size, then solves the subproblems. The solutions to the subproblems must then be patched back or nested together. For this method to produce good global solutions, the subproblems must naturally disjoint, and the division made must be an appropriate one (optimal). The optimization results in making the division such that the errors made in patching does not offset the gains obtained in applying more powerful methods to the subproblems. Roughly speaking, the algorithm of  uses a combination of the two outlined strategies. The problem is divided first into two oracle searches, then a third global search is used to patch back the two previous searches.

The present work gives an alternative scheme for dividing the original problem into subproblems. As will be seen later, our method allows one to relax some of the assumptions made in .

The paper is organized as follows: in section 2 we present the necessary tools for our algorithm, whereas the algorithm itself is presented in section 3. Using an average complexity analysis for our algorithm, we give in section 4 the running time for a variety of hard search problems. In the appendix we present a generalization of our algorithm to higher nesting levels.

## 2 Definitions

Search problems belong to the family of Constraint Satisfiability Problems (CSP). A search problem is defined as the problem of finding a satisfying assignment for a set (formula) of constraints acting on $`n`$ (qu)bits. In general, each constraint acts on a number of bits less than or equal to $`n`$. The search will terminate whenever the “program” finds an assignment, or assignments, that satisfy simultaneously all the constraints of the CSP formula. For example, the oracle search over an unstructured database corresponds to the problem of finding an assignment that satisfies all the constraints of the predefined formula, whose constraints are all $`n`$ bit constraints. Even though the number of qubits is fixed, by varying the length or the type of constraints the CSP can cover different types of problems.

For a structured database (problem) we can divide the initial search problem of $`n`$ qubits into subproblems with smaller number of variables. One can choose to test first if an assignment to $`n_A`$ variables out of the $`n`$ variables satisfies all the constraints acting uniquely on this subset of variables, labeled $`A`$. The result of this test is kept in a register. A similar test is then applied to the remaining bits, that we label as subset $`B`$. And finally one must nest both tests to find a global solution for the initial search problem. This is what is meant by a structured or nested search.

More specifically, assume that the initial search problem admits a number $`M_{AB}`$ of solutions. In a quantum search problem, a solution or a satisfying assignment is a state in the Hilbert space of the $`n`$ qubits. The set of solutions we thereafter denote by

$$_{AB}=\left\{|m_i/i=1,2,\mathrm{},M_{AB}\right\}.$$
(1)
All elements of this set satisfy the predefined formula of constraints denoted generically by $`\left\{C_{AB}\right\}`$. For our patching strategy, the constraints $`\left\{C_{AB}\right\}`$ are classified as

* $`\left\{C_A\right\}`$: the set of constraints acting only on the subset $`A`$. In general, there exist different solutions satisfying all the constrains of this set. The set of possible solutions is

$$_A=\left\{|m_i^A/i=1,2,\mathrm{},M_A\right\}.$$
(2)
* $`\left\{C_B\right\}`$: the set of constraints acting only on the subset $`B`$. Again, different solutions can be found to satisfy all the constraints of this set. The set of possible solutions is

$$_B=\left\{|m_i^B/i=1,2,\mathrm{},M_B\right\}.$$
(3)
* $`\left\{C_{A/B}\right\}=\left\{C_{AB}\right\}\left\{C_A\right\}\left\{C_B\right\}`$ is the set of constraints acting simultaneously on $`A`$ and $`B`$111This procedure is analogous to the resolution method used in Davis-Putnam procedure . Hence $`\left\{C_{A/B}\right\}`$ can then be called the reduction of $`\left\{C_{AB}\right\}`$ by $`A`$ and $`B`$..

Some elements of $`_A`$ may not satisfy one or more constraints in $`\left\{C_{A/B}\right\}`$, those elements can not give global solutions, so they are Not Solutions. The set of such elements is denoted by $`_A^{NS}`$. On the other hand, the set of elements of $`_A`$ satisfying all the constraints in $`\left\{C_{A/B}\right\}`$, and hence giving global solutions, will be denoted by $`_A^S`$. Therefore

$$_A=_A^S_A^{NS}.$$
(4)
Similarly we have

$$_B=_B^S_B^{NS}.$$
(5)
So, an element $`|m_i_{AB}`$ can be written as

$$|m_i=|m_A^S|m_B^S,$$
(6)
where $`|m_A^S`$ and $`|m_B^S`$ are elements of $`_A^S`$ and $`_B^S`$ respectively.

It is evident that for a given state $`|m_A^S`$ there could be different states $`|m_B^S`$ such that their product is an element of $`_{AB}`$. At some stages, the authors of  and  made the assumption that for a given $`|m_A^S`$ there is only one $`|m_B^S`$ such that their product is a global solution. Although one may argue that such an assumption is justified for hard search problems, where there is just one global solution, it is unlikely that it holds for a generic search problem. This assumption is not needed for our algorithm, and hence our algorithm applies to a broader class of search problems.

## 3 Nested algorithm

The present algorithm consists of two stages. In the first stage (stage $`I`$) we evolve adiabatically the quantum system from a state which is the ground state of a Hamiltonian that can be easily obtained to a state which is a product of states of $`_A`$ and states of $`_B`$. As previously mentioned, this product is not yet a global solution. In stage $`II`$ a global search, similar to what is labeled as stage C in , is applied to the output of stage $`I`$ to rotate it towards an element of $`_{AB}`$. The output of stage $`II`$ is a global solution of the problem.

### 3.1 Stage $`I`$

In this stage we use the procedure defined in . The Hamiltonian is split into two parts, each acting on one of the Hilbert spaces $`_A`$, of dimension $`N_A=2^{n_A}`$, and $`_B`$, of dimension $`N_B=2^{n_B}`$, spanned by subsets $`A`$ and $`B`$ respectively.
The most convenient initial state is the equal superposition of all possible pure states of the system

$`|\mathrm{\Psi }(t=0)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N_A}}}{\displaystyle \underset{i_A_A}{}}|i_A{\displaystyle \frac{1}{\sqrt{N_B}}}{\displaystyle \underset{i_B_B}{}}|i_B`$ (7)
$``$ $`|\mathrm{\Psi }_{0A}|\mathrm{\Psi }_{0B},`$
where $`N=2^n=N_AN_B`$. $`|\mathrm{\Psi }_0`$ is the ground state of the Hamiltonian

$$H_0=(I_A|\mathrm{\Psi }_{0A}\mathrm{\Psi }_{0A}|)I_B+I_A(I_B|\mathrm{\Psi }_{0B}\mathrm{\Psi }_{0B}|).$$
(8)
The initial Hamiltonian should evolve adiabatically into a final Hamiltonian which has the ground state

$$\frac{1}{\sqrt{M_A}}\underset{m_A_A}{}|m_A\frac{1}{\sqrt{M_B}}\underset{m_B_B}{}|m_B|\mathrm{\Psi }_{fA}|\mathrm{\Psi }_{fB}.$$
(9)
A possible final Hamiltonian is

$$H_f^I=(I_A|\mathrm{\Psi }_{fA}\mathrm{\Psi }_{fA}|)I_B+I_A(I_B|\mathrm{\Psi }_{fB}\mathrm{\Psi }_{fB}|).$$
(10)
We have chosen $`H_0`$ and $`H_f^I`$ so that the evolution is decoupled in $`_A`$ and $`_B`$. The time dependent Hamiltonian evolving $`H_0`$ to $`H_f^I`$ in a time $`T_I`$ is the linear combination

$$H(t)=(1s(t))H_0+s(t)H_f^If(s)H_0+g(s)H_f^I,$$
(11)
where $`s(t)`$ is chosen such that $`s(0)=0`$ and $`s(T_I)=1`$. When the adiabaticity condition holds, the initial ground state $`|\mathrm{\Psi }_0`$ will evolve slowly in time $`T_I`$ to the ground state of the final Hamiltonian. The meaning of “slow” evolution is quantified by the adiabatic theorem , which states that the accuracy with which the system remains in its instantaneous ground state is equal to the sum of the ratio of the transition matrix element from the ground state, with energy $`E_0`$, to any other state, with energy $`E_i`$, over the fourth power of the “radiated” energy $`\omega _{0i}=E_0E_i`$, i.e.

$$\underset{i}{}\frac{\left|E_i|\frac{dH}{dt}|E_0\right|^2}{\omega _{0i}^4}ϵ^2.$$
(12)
This gives a lower bound on the evolution time $`T_I`$.

The eigenvalues and eigenfunctions of the Hamiltonian $`H(t)`$ can be calculated analytically. We find that the only nonvanishing transition probabilities from the ground state are to the two lowest energy levels:

* $`E_A=E_1^A+E_0^B`$, i.e. the state where the subsystem $`A`$ is in its first excited state, and the subsystem $`B`$ is in its ground state
* $`E_B=E_0^A+E_1^B`$, i.e. the state where the subsystem $`A`$ is in its ground state, and $`B`$ is in its first excited state.

The corresponding transition matrix elements are given by

$$\left|E_i|\frac{dH}{dt}|E_0\right|^2=\left(\frac{ds}{dt}\right)^2\frac{M_i^2/N_i^2}{\omega _i^2}\left(\frac{N_i}{M_i}1\right)\frac{\xi _i^2}{\omega _i^2}\left(\frac{ds}{dt}\right)^2,i=A,B,$$
(13)
where

$$\omega _i=\sqrt{(fg)^2+\frac{4M_i}{N_i}fg}=\sqrt{(12s)^2+4\frac{M_i}{N_i}s(1s)}$$
(14)
is the gap (energy difference) between the first excited state and the ground state of the $`i`$-th subsystem.
By integrating eq. (12) over $`s`$ using eq. (13) the lower bound, or the minimal time, will be

$$T_I=\frac{1}{ϵ}_0^1𝑑s\sqrt{\underset{i=A,B}{}\frac{\xi _i^2}{\omega _i^6}}.$$
(15)
Note that the parallel evolution, in the two Hilbert spaces $`_A`$ and $`_B`$, of stage $`I`$ is the major feature that distinguishes our algorithm from that of . Stage $`I`$ replaces the two-stages evolution (sequential evolution) used in . A priori the sequential evolution seems more efficient, as one can use the result of the first stage to eliminate some of the no-good assignments in the second stage. However the algorithm used in  is not a tree-like search procedure, i.e. it does not pick an order in which it instantiates the variables. Hence the algorithm of  does not eliminate “bad” trees. Therefore replacing the sequential by the parallel evolution will not affect dramatically the time cost of the algorithm. On the other hand, as mentioned before , the symmetric parallel evolution allows one to analyse a broader set of structured search algorithms.

### 3.2 Stage $`II`$

The output state of stage $`I`$ can be written as the sum of a state encoding the solution to the search problem and a residual part, which is not a solution:

$$|\mathrm{\Psi }_{fA}|\mathrm{\Psi }_{fB}=\sqrt{\frac{M_AM_BM_{AB}}{M_AM_B}}|\mathrm{\Psi }^{NS}+\sqrt{\frac{M_{AB}}{M_AM_B}}|\mathrm{\Psi }^S,$$
(16)
where the “solution state” is

$$|\mathrm{\Psi }^S=\frac{1}{\sqrt{M_A^SM_B^S}}\underset{m_A^S_A^S}{}|m_A^S\underset{m_B^S_B^S}{}|m_B^S=\frac{1}{\sqrt{M_{AB}}}\underset{m_i_{AB}}{}|m_i$$
(17)
and we have used the fact that $`M_A^SM_B^S=M_{AB}`$.

At this stage one is inclined to start again an adiabatic evolution from a Hamiltonian whose ground state is $`|\mathrm{\Psi }_{fA}|\mathrm{\Psi }_{fB}`$, e.g.

$$H_i=1|\mathrm{\Psi }_{fA}\mathrm{\Psi }_{fA}||\mathrm{\Psi }_{fB}\mathrm{\Psi }_{fB}|$$
(18)
to a final Hamiltonian whose ground state is $`|\mathrm{\Psi }^S`$, e.g.

$$H_f=1|\mathrm{\Psi }^S\mathrm{\Psi }^S|,$$
(19)
in time $`O(\sqrt{M_AM_B/M_{AB}})`$. Unfortunately, this is not possible since the initial Hamiltonian $`H_i`$ is not accessible. Following  a global search can be applied to nest between the two subsets $`A`$ and $`B`$. The global search is achieved through the following procedure

* Stage $`I`$ is approximated by an evolution operator $`U`$ such that $`|\mathrm{\Psi }_{fA}|\mathrm{\Psi }_{fB}U|\mathrm{\Psi }_{0A}|\mathrm{\Psi }_{0B}`$.
* Hence $`H_iUH_0U^{}`$, note that $`H_i`$ is “replaced” by $`H_0`$ which is easily accessible.
* The adiabatic evolution from $`H_i`$ to $`H_f`$ is implemented on quantum circuits. The continuous evolution is replaced by a “step evolution” over intervals of time. In each step, at a given time $`t`$, we use the approximation

$$e^{iH_it}Ue^{iH_0t}U^{}$$
(20)
i.e. a backword evolution in time $`T_I`$, an application of $`\mathrm{exp}(iH_0t)`$, then a forward evolution in time $`T_I`$. This step requires a time $`O(T_I)`$.

The number of steps is chosen to minimize the error involved in the discretization needed during the global evolution. The number of iteration is at least (see ) of the order of $`\sqrt{M_AM_B/M_{AB}}`$. So that the total running time needed to get a global solution, is proportional to

$$T=T_I\sqrt{\frac{M_AM_B}{M_{AB}}}=T(M_A,M_B,M_{AB},n,n_A/n).$$
(21)
This expression is not yet the final answer, since there is still the problem of determining the values of $`M_{A,B,AB}`$. This will be discussed in the next section.

## 4 Complexity analysis

The total running time $`T`$ is a function of $`M_A,M_B,M_{AB},n`$ and $`x=n_A/n`$, among which $`n`$ is the only input parameter, and $`x`$ will be chosen to minimize $`T`$. Hence we are left with the three parameters $`M_A,M_B`$ and $`M_{AB}`$. A priori, these parameters are problem dependent and are difficult to determine. However, since we are interested in getting general results which depend as little as possible on the details of the problem, we can approximate these quantities using a complexity analysis.

Computational complexity theory studies the quantitative laws which govern computing. It seeks a comprehensive classification of problems by their intrinsic difficulty and an understanding of what makes these problems hard to compute. The simplified average complexity analysis derived in  can be adapted for our generic algorithm. The major simplification made in  is the approximation of independent no-good assignments222For most CSP there are classes of no-good assignments that can be deduced from each other. For example, one can use proposional reasoning to generate some no-goods from others, see for example chapter two of .. This allows for a relatively easy combinatoric analysis. Following  the unknown parameters are estimated by

$$M_j2^{n_jn\alpha (n_j/n)^k}\mathrm{for}j=A,B,AB\mathrm{with}n_{AB}=n,$$
(22)
where $`k`$ is the number of variables per constraint, assuming a constant length constraint. $`\alpha `$ represents the average difficulty of the problem, it characterizes the average number of no-good ground instances per variable. Hardest problems are found near a critical333This critical behavior is similar to a phase transition in condensed matter physics. Similar phase transition features were observed for random K-SAT problems . value $`\alpha =\alpha _c=1`$. The critical $`\alpha `$ is obtained from the number of solutions

$$M_{AB}N^{1\alpha }$$

for $`\alpha <1`$, $`M_{AB}`$ is large, i.e. the problem is under-constrained and easy to solve, while for $`\alpha >1`$, $`M_{AB}<1`$, i.e. the problem is over-constrained and probably has no solution. Finally for $`\alpha =1`$ we have one solution, this is the definition of the critical value.

Therefore, the optimization procedure reduces to a minimization of the running time as a function of the ratio $`x`$ for different values of $`\alpha `$ and $`k`$.

### 4.1 Approximate and exact scalings

The running time $`T`$ can not be calculated analytically. However by the following simple arguments it is possible to give an estimate for $`T`$:

* We first recognize that our algorithm is inefficient if $`M_i/N_i1`$. This can be easily seen by remembering that $`N_i`$ is the number of possible assignments in space $`A`$ or $`B`$ and $`M_i`$ is the number of solution in that space. For $`M_iN_i`$ any assignment will become a solution, then, one step (picking up any state) is sufficient to find a solution in the subspace $`i`$. This is reflected through the limit

$$\xi _i=\frac{M_i}{N_i}\sqrt{\frac{N_i}{M_i}1}0\mathrm{when}M_iN_i.$$

Therefore, the first stage does not eliminate any no-good assignment and consequently the second stage will take $`\sqrt{N}`$ steps to give a solution, i.e. as if the search is an oracle search. But what does this means?
Using the average complexity analysis of the previous section we write:

$$M_A/N_A=2^{n\alpha x^k}.$$

Taking $`M_A/N_A1`$ implies that $`n\alpha x^k1`$ or $`k\mathrm{ln}(n\alpha )/\mathrm{ln}(1/x)`$, on the other hand we expect444If $`x=0`$ or 1, the first stage is again useless as there is no division of the initial space. that optimization will give a value of $`x`$ not far from $`1/2`$, hence for large values of $`n`$ and near the critical $`\alpha =1`$ our algorithm breaks down for $`k\mathrm{ln}(n)`$. But the larger the $`k`$ is, the closer we are to the unstructured search where each constraint is an $`n`$ qubit constraint.
* From the previous remark we deduce that our algorithm is effective as long as $`M_i/N_i1`$.
* The running time of stage $`I`$, given in eq. (15), exhibits a near singular behavior for $`M_i/N_i1`$. This near singularity occurs when the integration variable $`s`$ is close to the value for which $`f(s)=g(s)`$, i.e. $`s=1/2`$.

These remarks imply that, for a given $`x`$, $`T_I`$ can be approximated using the near singular behavior to obtain

$$T_I(x)\sqrt{\mathrm{Max}(\frac{N_A}{M_A},\frac{N_B}{M_B})}.$$
(23)
In other words, the time $`T_I`$ scales with the square root of the largest value of the dimension of each Hilbert space divided by the number of possible solutions in that space. Multiplying by the number of iterations of stage $`II`$, the total running time reads

$$T(x)\sqrt{\mathrm{Max}(\frac{N_AM_B}{M_{AB}},\frac{N_BM_A}{M_{AB}})}.$$
(24)
Using the average complexity analysis we get

$$T(x)\sqrt{\mathrm{Max}(N^{\alpha \alpha (1x)^k},N^{\alpha \alpha x^k})}.$$
(25)
In figure 1 we compare this approximate formula with the time calculated numerically. The two results match extremely well for all values of $`x`$.

The time $`T`$ is then optimal for $`N^{\alpha \alpha (1x)^k}=N^{\alpha \alpha x^k}`$, i.e. $`x=1/2`$, which is compatible with the symmetric nature of the algorithm. This result agrees with the numerical results presented in figure 2, where $`x=1/2`$ is found to give the optimal time for different values of $`k`$. The plot also illustrates that the larger $`k`$ is the closer we are to an oracle type search, which is what we predicted before. Hence the optimal time is given by

$$T_{min}=T(x=1/2)N^{\alpha /2\alpha /2^{k+1}}2^{n\alpha (1/21/2^{k+1})}.$$
(26)
We conclude that the parameter of interest in these considerations is $`n\alpha `$. In figure 3 we plot the running time as a function of $`x`$ for different values of $`\alpha `$ and $`n`$ keeping their product constant. The fact that the different plots coincide strongly supports the scaling obtained from the analytical estimate.

Another important aspect is the (sub)-exponential growth of $`T`$ as a function of $`\alpha `$. This can be seen either from eq. (26) or from the “critical” behavior depicted in figure 4.

Finally we consider the special case $`(k=2,\alpha =1)`$. In this case the running time is $`O(N^{3/8})`$, which is better than the classical running time $`O(N^{0.5})`$ but less efficient than the result obtained using the algorithm of Cerf et al. which gives $`O(N^{0.31})`$. The $`O(N^{0.375})`$ scaling is compared with the numerical result in figure 5.

## 5 Conclusions

We have presented a new adiabatic quantum algorithm for searches over structured databases. Our algorithm is constructed with minimum assumptions about the nature of the database and the specificity of the search problem. Moreover, it results in a significant potential speedup over its classical counterpart.

The new algorithm is in fact “quasi”-adiabatic, it requires the use of quantum circuits at some stages. An interesting topic for future work would be to construct a “pure” adiabatic quantum algorithm. Another important issue is to explore the possibility of using smart classical strategies such as back tracking and constraint propagation, and to see the effect of such strategies on quantum interference and coherence. These and other related issues will be addressed elsewhere.

## Appendix A Multi-nesting verses multi-partition

The nesting procedure can be applied to subsets $`A`$ and $`B`$, resulting in higher level nesting or multi-nesting . In multi-nesting stages $`I`$ and $`II`$ are applied to subsets $`A`$ and $`B`$ separately. Then stage $`II`$ is applied to nest $`A`$ and $`B`$. Then the time $`T_I`$ is replaced by a shorter time. More nesting can be used to enhance the effectiveness of the structured search.

Another alternative is to modify stage $`I`$ by splitting the $`n`$ variables into more than two subset then

$$T_I=\frac{1}{ϵ}_0^1𝑑s\sqrt{\underset{i=A,B,C,D,\mathrm{}}{}\frac{\xi _i^2}{\omega _i^6}}$$
(27)
and hence applying the procedure of stage $`II`$ we get

$$T=T_I\sqrt{\frac{\underset{i=A,B,C,D,\mathrm{}}{}M_i}{M_{AB}}}=T(M_A,M_B,M_C,\mathrm{},M_{AB},n,n_A/n).$$
(28)
An average complexity analysis can then be applied to obtain an estimate for $`T`$.


text:
# The Changing AGN Population

## 1. Introduction

It is becoming increasingly clear that low-redshift accretion may be the dominant phase of supermassive black hole (SMBH) growth (Barger et al. 2001; Cowie et al. 2003). Although theoretical models of SMBH formation have relied on comparisons with the optical quasar luminosity function, optically accessible broad-line active galactic nuclei (AGNs) form only about 30% (Barger et al. 2003b) of the X-ray background (XRB). Explanations for the lack of broad lines in most AGN optical spectra include absorption along the line-of-sight or dilution of the emission-line signatures by stellar light (Moran, Filippenko, & Chornock 2002). Regardless of the reason, optically and soft X-ray selected samples substantially undercount the AGN population relative to hard X-ray selected samples, at least at intermediate luminosities (Cowie et al. 2003; Barger et al. 2003a).

The hardest band for which the hard X-ray luminosity function (HXLF) can presently be determined is rest-frame $`28`$ keV. The $`2`$ Ms exposure of the Chandra Deep Field-North (CDF-N) samples a large, distant cosmological volume down to a very faint hard X-ray flux limit ($`f_{28\mathrm{keV}}1.4\times 10^{16}`$ ergs cm-2 s-1; Alexander et al. 2003). However, this is substantially fainter than the X-ray flux of the sources that contribute the most to the XRB (Cowie et al. 2002), and the low-redshift volume sampled is small. To sample a large cosmological volume at low redshifts, we have mapped with Chandra an $`0.4`$ deg2 area in the Lockman Hole-Northwest (LH-NW) field to a flux level of $`f_{28\mathrm{keV}}3\times 10^{15}`$ ergs cm-2 s-1 (Yang et al. 2003), where most of the hard XRB is resolved. In this paper we use our LH-NW data, together with deep Chandra data on the CDF-N, A370, SSA13, and SSA22 fields and bright ASCA data from Akiyama et al. (2000), to construct rest-frame $`28`$ keV luminosity functions at low redshifts. We then investigate the evolution with redshift of AGN populations separated by optical spectral type. We assume $`\mathrm{\Omega }_M=1/3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=2/3`$, and $`H_o=65`$ km s-1 Mpc-1. We use $`L_x`$ to denote the rest-frame $`28`$ keV luminosity.

## 2. X-ray Sample Selection

The LH-NW data are taken from the X-ray and optical catalogs of Y. Yang et al. (in preparation) and A. T. Steffen et al. (in preparation). Here we consider only sources with fluxes above $`5\times 10^{16}`$ ergs cm-2 s-1 ($`0.52`$ keV) and $`5\times 10^{15}`$ ($`28`$ keV) to provide good, high-significance samples. We also only consider those sources that have been targeted spectroscopically. The only source property we considered when making our masks for spectroscopy was X-ray flux (sources with higher X-ray fluxes had higher priority when we had to choose one X-ray source over another due to mask conflicts). We spectroscopically observed 150 sources out of our total $`28`$ keV sample of 228. We were able to confidently identify 100 of these 150 sources (67%) based on multiple emission and/or absorption line features. Two of the identified sources are stars. A practical optical limit for obtaining redshifts for X-ray sources using the Deep Extragalactic Imaging Multi-Object Spectrograph (DEIMOS; Faber et al. 2002) on Keck II is about $`R=24`$ (see Barger et al. 2003b). Other spectra may not have been identifiable due to slits falling too close to an edge of the CCD or bright neighbor objects compromising the spectra.

The faintest X-ray data are taken from the $`2`$ Ms CDF-N catalogs of Alexander et al. (2003; X-ray) and Barger et al. (2003b; optical). Excluding stars and the small number of extended X-ray sources, the $`28`$ keV sample contains 326 sources, 68% of which have either spectroscopic (165) or reliable photometric (56) redshifts. Our other Chandra and ASCA samples are summarized in Cowie et al. (2003). Hereafter, we refer to the A370, SSA13, and SSA22 Chandra samples collectively as “other”.

In our analysis we place every X-ray source that was not spectroscopically identified for any reason into each and every redshift interval in order to obtain a strict upper bound on the number of AGNs that could lie in that redshift interval, modulo the existence of a substantial population of Compton-thick sources not detected by Chandra. We note in passing that Barger et al. (2002) found from the 1 Ms CDF-N data that a substantial population of Compton-thick sources is not needed to explain the greater than $`8`$ keV XRB, provided that the ASCA and BeppoSAX XRB measurements, instead of the lower HEAO-1 A2 measurement, are assumed at energies less than $`8`$ keV. For the LH-NW field, we only ever consider the sources that were spectroscopically observed (identified or not). This is acceptable because the spectroscopy are purely from targeted observations of the X-ray samples. For the other fields, the spectroscopy are not purely from targeted observations of the X-ray samples, so we need to include all of the sources.

We determined the solid angle covered by the combined samples at a given flux by comparing the observed numbers of sources versus flux with the averaged number counts in the $`28`$ keV band from Cowie et al. (2002). This method allows a simple treatment of the incompleteness that was modeled in computing the counts. However, the counts in Cowie et al. also include the low CDF-S counts, which may affect the normalization at the 10% level relative to the LH-NW average (Yang et al. 2003). We consider this to be a reasonable estimate of the systematic errors in the present analysis. The solid angle covered by the combined $`28`$ keV samples ranges from 0.002 deg2 at the faintest fluxes to 5.8 deg2 at the highest fluxes. At $`2.3\times 10^{14}`$ ergs cm-2 s-1 ($`28`$ keV) the solid angle is 0.52 deg2.

In Figure 1 we show redshift versus $`28`$ keV flux for the ASCA (circles), LH-NW (diamonds), CDF-N (squares), and other (triangles) samples. The LH-NW sources nicely fill in the almost an order of magnitude flux gap between the CDF-N and ASCA samples. The solid curves correspond to loci of constant $`L_x`$. Any source more luminous than $`L_x=10^{42}`$ ergs s-1 is very likely to be an AGN on energetic grounds (Zezas, Georgantopoulos, & Ward 1998; Moran et al. 1999), though many of the intermediate luminosity sources do not show obvious AGN signatures in their optical spectra.

## 3. The AGN X-ray Luminosity Function

Cowie et al. (2003) used Chandra, ASCA, and ROSAT data to determine rest-frame $`28`$ keV AGN luminosity functions in two redshift intervals, $`z=0.11`$ and $`z=24`$. They did not have very many sources in their sample with $`28`$ keV fluxes $`2\times 10^{14}`$ to $`10^{13}`$ ergs cm-2 s-1 since that intermediate flux range is not well covered by either small-area, deep Chandra observations or bright ASCA observations. With our wide-area LH-NW data, we now have good coverage of the intermediate flux range and can therefore improve the determination of the $`z=0.11`$ HXLF.

In Figure 3a we show our measured HXLF for the $`z=0.11`$ interval (solid diamonds), computed following Cowie et al. (2003), who used the traditional $`1/V_a`$ method of Felten (1977). The uncertainties are Poissonian, based on the number of galaxies in each logarithmic luminosity bin. However, since incompleteness is a potentially larger source of error, we also recomputed the HXLF (open diamonds) by assigning all of the unidentified sources to the volume-weighted centers of each and every redshift interval. Because the spectroscopic identifications are much more complete at the high X-ray fluxes, the associated systematic uncertainties are larger at low $`L_x`$. The HXLF of only the broad-line AGNs is shown by the large open squares in Figure 3a. In the high luminosity range where the HXLF is steep, we can see that the broad-line AGNs dominate the total HXLF.

The uncorrected $`z=0.11`$ HXLF can be well represented by a conventional broken power law fit of the form $`6.6\times 10^6L_{44}^{1.37}`$ per Mpc3 above $`2\times 10^{43}`$ ergs s-1 and $`1.5\times 10^5L_{44}^{0.87}`$ below, where $`L_{44}`$ is the luminosity in units of $`10^{44}`$ ergs s-1. The integrated $`z=0.11`$ light density is $`1.8\times 10^{39}`$ ergs s-1 Mpc-3. Half of the total uncorrected light in the $`z=0.11`$ interval arises in sources more luminous than $`2.1\times 10^{43}`$ ergs s-1, and about one-third arises in broad-line AGNs.

In Figure 3b we show our measured HXLF over the luminosity range $`L_x=10^{43}`$ to $`10^{45}`$ ergs s-1, where most of the sources have been spectroscopically identified, divided into two smaller redshift intervals, $`z=0.10.5`$ (squares) and $`z=0.51`$ (diamonds). For comparison, we show the local Seyfert galaxy HXLF (solid line) from Piccinotti et al. (1982; converted to $`H_0=65`$ km s-1 Mpc-1 and, using their assumed $`\mathrm{\Gamma }=1.7`$, to $`28`$ keV) and the high-redshift HXLF (dashed line) from Cowie et al. (2003; updated to include our LH-NW data and the 2 Ms CDF-N data, and to cover the interval $`z=1.53`$). We have computed the maximal high-redshift HXLF (dotted line) by assigning all of the unidentified sources to the volume-weighted centers of each and every redshift interval. The high-redshift HXLF is much tighter than in Cowie et al. (2003) because of our slightly lower redshift interval and our higher spectroscopic completeness. From Figure 3b we see that the number density of $`L_x>10^{44}`$ ergs s-1 sources steadily declines with decreasing redshift. However, the number density of $`L_x=10^{43}10^{44}`$ ergs s-1 sources is largest in the $`z=0.51`$ interval and smaller at both higher and lower redshifts (see also Fig. 5).

## 4. X-ray Fluxes and Luminosities

Prior to Chandra and XMM-Newton, various groups (e.g., Comastri et al. 1995; Gilli, Salvati, & Hasinger 2001) using XRB population synthesis models based on unified AGN schemes predicted that a mixture of unobscured type I AGNs and intrinsically obscured type II AGNs were needed to produce the bulk of the XRB. The obscured systems were expected to lie at high redshifts ($`z=23`$), with the high luminosity end (type II quasars) contributing substantially to the XRB. Contrary to the predictions, optical observations of the sources detected in deep Chandra and XMM-Newton hard X-ray surveys have found a redshift distribution that peaks at $`z<1`$ (Barger et al. 2002, 2003b; Hasinger 2003). Only a handful of type II quasars (here we define quasars as having $`L_x>10^{44}`$ ergs s-1) have been detected, and these contribute only a small fraction of the XRB (e.g., Barger et al. 2003b).

In Figure 4 we plot for the redshift intervals $`z=0.11`$ and $`z=1.53`$ the fractional number of X-ray sources in each luminosity bin, by optical spectral type, versus rest-frame $`28`$ keV luminosity. From Figure 4a we see that the dominant population at quasar luminosities, where our redshift identifications are very complete, is broad-line AGNs. However, at more Seyfert-like luminosities ($`L_x<10^{44}`$ ergs s-1), the broad-line AGN fraction declines rapidly with decreasing luminosity. Figure 4b shows a similar trend with luminosity. These results suggest a luminosity dependence in optical spectral type: type I AGNs dominate at quasar luminosities, while type II AGNs can only become an important component of the X-ray population at Seyfert-like luminosities.

To quantify this, we note that in the $`z=0.11`$ redshift interval there are 20 sources with quasar luminosities, of which 15 are broad-line AGNs. Only three of the unidentified sources would have quasar luminosities if placed in this redshift interval, so the systematic errors are small, and the bulk of the quasars are broad-line AGNs. In the $`z=1.53`$ interval there are 65 sources with quasar luminosities, of which 58 are broad-line AGNs. There are 56 unidentified sources that would have quasar luminosities if placed in this redshift interval, so the systematic uncertainties are larger. However, at least half (and probably considerably more) of the sources with quasar luminosities in this redshift interval are broad-line AGNs.

## 5. Discussion

Barger et al. (2003a) graphically showed how the $`L_x=10^{43}10^{44}`$ ergs s-1 (Seyfert-like luminosities) number density slowly increases with decreasing redshift from $`z=6`$ to $`z=0.1`$, while the $`L_x=10^{44}10^{45}`$ ergs s-1 (quasar luminosities) number density peaks at $`z=1.53`$. They also showed that the evolution of the $`L_x=10^{44}10^{45}`$ ergs s-1 number density closely matches in shape the evolution of optically selected quasars. Given these trends and our interesting result from Figure 4 that while the dominant population at quasar luminosities is broad-line AGNs, at lower luminosities the broad-line fraction drops dramatically, here we explore a different aspect of the number density evolution.

In Figure 5 we show number density versus redshift for the sources with $`L_x>10^{43}`$ ergs s-1. Rather than dividing the sources into luminosity classes, as was done in Barger et al. (2003a), we have instead divided the sources into broad-line (open squares) and non-broad-line (solid circles) AGNs. The horizontal bars show the maximal number densities for the non-broad-line AGNs, obtained by assigning all of the unidentified sources redshifts at the volume-weighted centers of each and every redshift bin and then retaining only those with $`L_x>10^{43}`$ ergs s-1. The bars are not consistent with one another because all of the unidentified sources are included in each and every redshift bin, as long as their luminosities at those redshifts exceed $`10^{43}`$ ergs s-1.

From Figure 5 we see that the number density of non-broad-line AGNs with $`L_x>10^{43}`$ ergs s-1 rises with decreasing redshift, until the last redshift bin. Strikingly, however, the number density of broad-line AGNs remains roughly constant with redshift while their average luminosities (dashed line; right-hand axis) decline at the lower redshifts. This is another example of cosmic downsizing. In the terminology of previous quasar studies, the broad-line AGNs are undergoing pure luminosity evolution, not density evolution. Pure luminosity evolution is consistent with the most recent optical determinations (Boyle et al. 2000). The increase in the number density of non-broad-line sources and the decrease in the broad-line AGN average luminosities results in non-broad-line sources increasingly dominating the light at lower redshifts, modulo the existence of a substantial population of Compton-thick sources.

We gratefully acknowledge support from NASA’s National Space Grant College and Fellowship Program and the Wisconsin Space Grant Consortium (A. T. S.), CXC grants GO2-3191A (A. J. B.) and GO2-3187B (L. L. C.), NSF grants AST-0084847 (A. J. B.) and AST-0084816 (L. L. C.), the University of Wisconsin Research Committee with funds granted by the Wisconsin Alumni Research Foundation (A. J. B.), the Alfred P. Sloan foundation (A. J. B.), and the IDS program of R. F. M.


text:
# Chandra Observations of the Dwarf Nova WX Hyi in Quiescence

## 1 Introduction

Cataclysmic Variables (CVs) are a class of interacting binaries in which a donor star transfers mass onto a white dwarf (WD) accretor (see Warner 1995 for a comprehensive review). The characteristics of the accretion flow depend mostly on the rate of accretion onto the WD and on the WD magnetic field strength. In non-magnetic systems ($`B10^4`$ G), accretion from the donor proceeds through an undisturbed accretion disk, which connects to the WD through a boundary layer (Patterson & Raymond 1985). When the WD is more strongly magnetized, the accreting material is forced to follow the topology of the field lines before it hits the surface (Aizu 1973).

As a dwarf nova of the SU UMa type, WX Hyi (discovered by Luyten 1932) belongs to the class of non-magnetic CVs. Its orbital period is 0.0748134 days, and because it does not show eclipses, only an estimated inclination angle of $`40\pm 10^{}`$ is available (Schoembs & Vogt 1981). Schoembs & Vogt (1981) derive a white dwarf mass of 0.9$`\pm `$0.3 M. Normal outbursts occur on average every 11.2 days (Ak, Ozkan & Mattei 2002). Distance estimates range from 120 pc based on the correlation between H$`\beta `$ equivalent width and mass transfer rate, $`\dot{M}`$ (Patterson 1984), and 265 pc based on the correlation between period and absolute magnitude during outburst (Warner 1987). We adopt 265 pc because the latter correlation is tighter. Patterson (1984) estimates $`\dot{M}=1.8\times 10^{16}\mathrm{g}\mathrm{s}^1`$ from the H$`\beta `$ equivalent width, which would give an X-ray luminosity of $`2\times 10^{33}\mathrm{erg}\mathrm{s}^1`$ if half the accretion energy emerges as X-rays. However, Schwope et al. (2002) find an X-ray luminosity of only $`4\times 10^{31}\mathrm{erg}\mathrm{s}^1`$ from the ROSAT All-Sky Survey data (assuming a 265 pc distance), which yields $`\dot{M}3\times 10^{14}\mathrm{g}\mathrm{s}^1`$. K. Long (2003, personal communication) finds accretion rates of order $`10^{15}\mathrm{g}\mathrm{s}^1`$ from fits of accretion disk models to the UV continuum. The optical emission lines are single-peaked (Schoembs & Vogt 1981), perhaps placing WX Hyi in the SW Sex class of CVs. Strong, rather broad UV emission lines are seen in IUE spectra, with the N V line being stronger than is typical of dwarf novae in quiescence (Hassall et al. 1983).

Here we are interested in the structure of the accretion flow in the quiescent state. In quiescence, the accretion rate in the disk is low ($`\dot{M}10^{16}`$ g s-1) and the disk is too cool to contribute any significant X-ray emission, while on the other hand the boundary layer is expected to be optically-thin and thus emit copious amounts of X-rays. According to the standard accretion theory (e.g. Lynden-Bell & Pringle 1974), the gravitational energy from the accreting material should reemerge in roughly equal proportions from the disk and from the boundary layer (BL). Whereas initial Einstein data on quiescent CVs appeared to be consistent with emission from an hot, optically-thin BL (Patterson & Raymond 1985), later observations with Ginga, ROSAT, ASCA, RXTE showed that the X-ray luminosity is actually lower than theoretically expected. Several suggestions have been made to solve this “mystery of missing BLs” (e.g. Ferland et al. 1982), such as (among others) reflection effects and cooling flows (Done & Osborne 1997), and disruption of the inner disk region by the WD rotation (Ponman et al. 1995), by irradiation from the WD (King 1997), or by “coronal” evaporation (Meyer & Meyer-Hofmeister 1994).

The presence and structure of BLs (or more generally hot flows) in quiescent CVs is therefore still a subject of debate, and wildly different models have been proposed (e.g. Narayan & Popham 1993, Meyer & Meyer-Hofmeister 1994; Mahasena & Osaki 1999; Medvedev & Menou 2001). While all the models can account for the total observed X-ray luminosity, each of them, however, makes very specific predictions for the physical conditions within the flow, and hence the resulting temperature and density profiles. As shown in a number of papers (among the most recent ones: Pradhan 2000; Mauche et al. 2001; Perna et al. 2000; Menou et al. 2001; Porquet et al. 2002; Dubau & Porquet 2002; Szkody et al. 2002), a powerful way to constrain the characteristics of hot X-ray emitting flows is through the relative intensities of emission lines produced within them. The high resolution of Chandra gives us the opportunity, for the first time, to resolve relative line intensities and perform such an analysis.

In this paper, we report Chandra observations of the dwarf nova WX Hyi in quiescence. Its X ray spectrum shows many lines from different ionization states of N, O, Ne, Mg, Si, S and Fe. First, we use line diagnostics to set overall constraints on the range of temperatures and densities spanned by the flow. Second, we compute X-ray spectra for different types of accretion flow structures, and test them against our data. The comparison shows that none of the models provides an adequate fit to the entire emission line spectrum, in that more plasma must be present at relatively low temperatures than predicted by the models. A cooling flow model (used by Mukai et al. 2003 to model spectra of non-magnetic CVs) is consistent with the short-wavelength spectra, but fails to reproduce the strength of the O VII line. This line is much stronger and broader than in the objects studied by Mukai et al., and we consider the intriguing possibility that it might be produced in a wind from the accretion flow.

## 2 Data Analysis

WX Hyi was observed on-axis by the Chandra X-ray Observatory’s ACIS-S detector and HETG grating on 2002 Jul 25 and Jul 28 (observation ids 3721 and 2760) and with effective exposure times after standard processing of 48366s and 49266s respectively. Optical data provided by the AAVSO show that WX Hyi was in quiescence during both observations, with an outburst 4 days after the second observation on August 2. The previous detected outburst was 20 days before our first observation. Examination of the background light curves shows no significant variability due to dropouts or flares. The destreak algorithm was applied to remove electronic noise from the noisier S4 chip, and the standard extraction parameters were used to obtain dispersed count spectra. The $`1`$ and $`+1`$ orders for both observations were coadded to obtain summed HEG and MEG spectra. There was no significant difference in the spectra between the two exposures. The secular degradation of the ACIS quantum efficiency was taken into account using the acis-abs model while making the response files.

The spectra were grouped to a signal-to-noise of 2 per bin. Initial attempts to model the lines in the merged spectrum using the CIAO fitting program Sherpa using both global model fits and fits to individual lines with a gaussian line model and polynomial fits to the local continuum gave unsatisfactory, poorly constrained results for the weaker lines. Instead, we determined line fluxes by simple integration of the counts in a band around each line, subtracting the continuum contribution using two methods: by comparison with counts in nearby wavelength bands and by subtraction of a broken powerlaw fit to the global continuum. These two methods gave consistent results, and are in reasonable agreement with the gaussian line fits in the case of the strongest lines. In addition, measurements of the stronger lines in the individual single order exposures give consistent results. Line flux measurements are presented in Table 1 together with probable identifications. The errors in Table 1 are 1-sigma statistical uncertainties computed for the line and continuum counts using the Gehrels (1986) approximation suitable for low count Poisson statistics.

The cleanest X-ray line is the O VIII resonance line at 18.967 Å. Its width is $`900\mathrm{km}\mathrm{s}^1`$ (FWHM) after correction for the instrument profile. Correcting for the $`40^{}`$ inclination, the intrinsic width is $`1400\mathrm{km}\mathrm{s}^1`$. While this is similar to the widths ($`1200\mathrm{km}\mathrm{s}^1`$) of the UV lines, the wings ($`4000\mathrm{km}\mathrm{s}^1`$) apparent in the UV line profiles (K. Long 2003, private communication) are absent. The inclination-corrected line width is $`1/3`$ the Keplerian velocity near the WD surface.

## 3 Constraints on the characteristics of the accretion flow from the X-ray spectrum

### 3.1 Line diagnostics

The X-ray spectrum, including both the MEG and HEG observations, is displayed in Figure 1. The spectrum shows a smooth continuum with strong emission lines, particularly from the H and He-like ions of O, Mg, Ne, Si, Fe, as well as from all the Fe L-shell complex (Fe XVII-Fe XXIV). The presence of this variety of ionization states indicates emission from a plasma distributed over a wide range of temperatures, from $`T10^6`$ K, necessary to have some contribution to the O VII line, up to $`T10^8`$ K, a temperature at which the ion Fe XXVI is most abundant in equilibrium (without strong photoionizing radiation).

From the X-ray spectrum in Figure 1, one can see that the O VII line profile appears to be wider that those of O VIII and the other lines, though it is clearly noisy. The profile, which is separately displayed in Figure 2, seems to consist of narrow O VII resonance and forbidden lines on top of a broad feature perhaps 5000 $`\mathrm{km}\mathrm{s}^1`$ wide. Some of the noise level results from the presence of bad columns in the $`+1`$ order of the MEG, which are broadened by dithering. However, careful examination of the counts in the $`1`$ order in the individual exposures indicates that the broad feature is probably real, though details of the structure are questionable. Note that the effective area changes rapidly in this wavelength range, making the number of counts in the blue-shifted part of the feature more significant. It is difficult to determine an unambiguous value for the uncertainty, in that different choices for spectral binning yield different estimates of the significance of the feature. There are several lines of calcium near this wavelength, but comparison with iron lines formed at the same temperatures indicates that the calcium lines are too weak to account for the observed structure. It should be noted, however, that there appears to be also a feature around 20.3 Å for which we have no identification.

If the broad O VII component is genuine, the lack of broad components in the other lines suggests that the high velocity material only appears in low ionization species. The lack of this component in N VII must be due to the small number of counts in that line. The line width is comparable to or larger than the widths ($`1200\mathrm{km}\mathrm{s}^1`$) of the UV Lines (K. Long 2003, personal communication), so the broad O VII might originate in the innermost part of the disk or in a wind.

As a first attempt to constrain the physical conditions in the accretion flow, we use the diagnostics provided by the He-like triplet line ratios, which have been widely used in analyses of the solar plasma (Gabriel & Jordan 1969; Mewe & Schrijver 1978; Doyle 1980; Pradhan & Shull 1981). In particular, the ratio $`Rf/i`$ between the forbidden line $`f`$ and the intercombination line $`i`$ is a strong function of the density at high densities, due to the suppression of the forbidden lines with increasing density. The ratio $`R`$ is typically on the order of a few at low densities, and rapidly drops at densities larger than a critical value $`n_c`$ that increases with $`Z`$, ranging from $`n_c6\times 10^8`$ cm-3 for C, to $`n_c3\times 10^{17}`$ cm-3 for Fe (e.g. Porquet et al. 2001). We find the 1$`\sigma `$ upper limits $`R<0.47`$ and $`R<0.92`$ for Ne111It should be noted that the intercombination line of Ne is not very clean, and might have some contamination from Fe XIX which is not possible to quantify. and Mg respectively. The 1s2s$`{}_{}{}^{3}S`$ \- 1s2p$`{}_{}{}^{3}P_{1}^{}`$ transition, however, is driven not only by collisions, but also by UV radiation, so that in gas near a 30,000 K black body the $`f/i`$ ratio will be in the high density limit for elements through Mg (Mauche 2002). We do not know the white dwarf temperature for WX Hyi. K. Long (private communication, 2003) has fit UV continua of WX Hyi observed by the Hopkins Ultraviolet Telescope and by the Hubble Space Telescope GHRS with combinations of white dwarf and accretion disk spectra, and typically finds WD temperatures above 20,000 K. The $`f/i`$ ratios of O VII, Ne IX and Mg XI are probably dominated by UV radiation, but the $`f/i`$ ratio for Si XIII should be dominated by collisions, even if the entire UV continuum in the Hopkins Ultraviolet Telescope wavelength range arises from the white dwarf. Unfortunately, we are unable to obtain a reliable measurement of the Si XIII ratio.

Another useful plasma diagnostic is the ratio $`G(i+f)/r`$, where $`r`$ is the resonance line. This is sensitive to the ionization state of the gas and to the electron temperature. When the resonance line is strong compared to the forbidden and the intercombination line ($`G1`$), the plasma is collision-dominated. On the other hand, a plasma in which photoionization is important will have a weaker resonance line. Values of $`G>4`$ are typically considered to be indicative of a photoionization-dominated plasma (e.g. Porquet et al. 2001). For O, Si and S, we found that the $`i`$, $`r`$ and $`f`$ lines were not well resolved to allow a reliable estimate of the $`G`$ value. In the cases of Mg and Ne, we were able to set the 1$`\sigma `$ upper limits $`G<0.73`$ and $`G<1.08`$, respectively. This indicates that photoionization is not dominant.

A useful density diagnostic is provided by the ratio between the 17.10 Å and the 17.05 Å lines of the ion Fe XVII, which is less sensitive to UV radiation than are the He-like line ratios (Mauche et al. 2001). The Fe XVII ratio $`I(17.10)/I(17.05)`$ is found to be $`0.6`$ (at the $`1\sigma `$ level) for WX Hyi. For a collisionally-dominated plasma, this value suggests a density $`n3\times 10^{13}\mathrm{cm}^3`$ in the relatively cooler region ($`T3`$$`6\times 10^6`$ K) where Fe XVII is produced. A comparable density is implied by the best value of the ratio between the 11.92 Å and the 11.77 Å lines of Fe XXII, which is found to be $`0.68\pm 0.4`$ in our data. Accounting for the error in the measurement, the $`1\sigma `$ confidence level for the density ranges between $`10^{13}10^{14}\mathrm{cm}^3`$. As discussed by Mauche, Liedahl & Fournier (2003), this density determination is very insensitive to temperature and photoexcitation.

### 3.2 Comparison with specific models for the X-ray emission

The physics, and hence the density-temperature structure, of boundary layers is still poorly understood. Only a few models have been proposed, and they have been primarily tested on their ability to account for the total observed X-ray and UV radiation. The Chandra observations of the X-ray emission, by resolving lines, allow us to set much tighter constraints on the models, due to the high sensitivity of the relative line intensities to the temperature and density profile characterizing the flow. Here we consider five types of models proposed in the literature to explain the origin of the quiescent X-ray emission of CVs. For each of them, we compute the expected X-ray spectrum using the temperature and density profiles that they predict, and compare the spectra with our observations. The cooling flow model has been implemented in the xspec package (Arnaud 1996). For the other models we compute the emissivity of the medium by using the atomic rate packages of the shock model code described by Raymond (1979), modified with updated atomic rates, as in Cox & Raymond (1985). This package has the disadvantage that some Fe L-shell multiplets that are resolved in our spectra are computed as single lines. The package also lacks recent improvements to the atomic rates of the more complex Fe ions. It computes an ionization balance which is quite similar to that of Mazzotta et al. (1998) in ionization equilibrium. The advantages of the code are that it allows us to compute time-dependent ionization and to include photoionization. It includes dielectronic recombination satellite lines and the contributions of recombination to excited levels for the spectra of H-like and He-like ions. Overall, the predictions should be reasonably accurate for most of the strong lines; those of H-like and He-like ions, Fe XVII and Fe XXIV. For the models presented here we included the effects of photoionization by the X-ray radiation produced within the flow itself, though this turned out to have little effect on the emission lines. For a proper comparison with the data, the computed X-ray spectra were then convolved with the response matrices and effective areas of the Chandra detectors. In particular, we used the HEG detector for $`\lambda 3`$ Å, and the MEG detector for $`\lambda >3`$ Å.

(a) Cooling flow

A cooling flow model, based on Mushotzky & Szymkowiak (1988), has been implemented in xspec with the routine mkcflow. Its basic assumption is of steady-state isobaric radiative cooling. The two main parameters are the maximum temperature, $`T_{\mathrm{max}}`$, and the overall normalization, which directly relates to the mass accretion rate. Mukai et al (2003) showed that cooling flow models provide a good representation of the X-ray spectral properties of non-magnetic CV systems. A fit to our data with the model mkcflow is shown in the top panel of Figure 3. We find that the fit is acceptable at short wavelengths but poorly constrained, with the maximum temperature parameter at a value of 20 keV, uncertain by a factor of two, and the column density $`N_\mathrm{H}=(2\pm 2)\times 10^{20}`$ cm-2 (basically consistent with zero). The value of the temperature is consistent with the Fe XXVI/Fe XXV ratio that gives an upper limit of 10 keV. Note that this ratio is about 1/7 in WX Hyi as opposed to 3/4 in U Gem (Szkody et al. 2002), implying a corresponding lower temperature of WX Hyi. The value of $`\dot{M}`$ implied by the cooling model is $`2\times 10^{14}\mathrm{g}\mathrm{s}^1`$ for $`T_{\mathrm{max}}=20`$ keV (assuming that half of the X-rays are absorbed by the WD). This is about a factor of 5 smaller than the value found by K. Long (2003, private communication) from fits of disk models to the UV continuum. While the system is by definition variable, it is not surprising that the X-ray luminosity is less than expected from $`\dot{M}`$ farther out in the disk. Possible interpretations are that $`\dot{M}`$ decreases inward in the disk as material builds up for the next outburst, that some material and energy is lost to a wind, or that some energy is advected to the WD surface and reemitted in the UV.

Overall, the cooling flow model appears to account reasonably well for the continuum and the line strengths in the short wavelength region of the spectrum, but it underpredicts the emission at long wavelengths. In particular, it is not able to account for the emission from the O VII line, which in WX Hyi is significantly stronger than in the objects studied by Mukai et al. (2003).

(b) Hot boundary layer

Narayan & Popham (1993; NP93) describe detailed models for thin disk boundary layers in CVs. A key element of their theory is the role of heat advection which, they show, allows boundary layers at low accretion rates to be radially thicker and significantly hotter than those at higher accretion rates. We have used the profiles of density and temperature corresponding to the model with $`\dot{M}=2\times 10^{15}\mathrm{g}\mathrm{s}^1`$ , which is the closest to the accretion rate in our system. The behavior of the low temperature region closest to the white dwarf surface cannot be read from the plots in NP93, but this is a small fraction of the total X-ray luminosity and may be neglected. We note that the comparison with this particular solution should be interpreted with caution, as temperature and density do not simply scale with accretion rate. Panel (b) of Figure 3 shows a comparison between our data and the X-ray spectrum that we computed for this boundary layer solution. The spectrum has been rescaled to match the continuum. If the density profile were to scale linearly with the accretion rate, the data would imply a value of $`\dot{M}`$ of about $`2\times 10^{14}`$ g s-1. As the figure shows, this model predicts substantially more emission than the cooling flow in the O VII and O VIII lines. This is a result of the fact the BL solution of NP93 predicts an increase of the density in the outer, colder regions of the flow where the oxygen emission peaks. However, the reasonable agreement between the observed and predicted low temperature lines is partly fortuitous in that Narayan & Popham assumed bremsstrahlung cooling, while the actual line cooling rate at the temperatures where the O VII, O VIII and Fe XVII lines are formed is nearly an order of magnitude higher. Thus if Narayan & Popham had used a cooling rate consistent with the lines we observe, their model would have predicted weaker emission by an order of magnitude.

(c) Coronal siphon flow

Meyer & Meyer-Hofmeister (1994) investigated a model based on the interaction between a cool disk and a hot corona above it via thermal conduction. They showed that, in these conditions, a coronal “siphon flow” can develop. This evaporates mass from the disk, which is partly accreted onto the white dwarf and partly lost to a wind. The solutions for this model are uniquely determined once the viscosity parameter $`\alpha `$ for the flow is specified222Under the conditions in which the radiative energy loss in the corona is small compared to conductive energy loss towards the lower boundary and the wind loss, the radial solutions for $`T`$ and $`\rho `$ are independent of $`\dot{M}`$ for a given white dwarf mass (Meyer & Meyer-Hofmeister 1994).. We varied the value of this parameter to find the best fit coronal flow solution to our data. The resulting X-ray spectrum predicted for the coronal siphon flow is shown in the middle panel of Figure 3. As it can be seen, this model falls short of producing enough emission in the long wavelength region of the spectrum. This is a consequence of the fact that this solution predicts a radial temperature profile $`Tr^1`$, and a run of density which is a very steep function of the radius, $`nr^3`$. This implies that there is relatively little mass at the lower temperatures where the longer wavelength emission is generated.

(d) X-ray emitting corona

Mahasena & Osaki (1999) constructed steady-state solutions of X-ray emitting corona, including the effect of thermal conduction. We have used the density and temperature profiles shown in their Figure 1 for an accretion rate $`\dot{M}=10^{15}`$ g s-1, and calculated the corresponding X-ray spectral emission. It was not possible to read numbers from their plot for the region below $`6\times 10^6`$ K, so we computed a matching solution that neglected the rotation and infall velocities, both of which are small in this region. The result is compared to the data in panel (d) of Figure 3. This comparison should also be taken with caution because the value of the accretion rate is not precisely tuned to match that of WX Hyi. This coronal model, which transports most of the accretion luminosity to a cool, thin, conduction dominated layer near the white dwarf surface and radiates it there, produces emission lines that are far too strong and a continuum that is much too soft. The higher accretion rate model of Mahasena & Osaki, $`\dot{M}=10^{16}`$ g s-1, did a much better job matching the Chandra spectrum of U Gem (Szkody et al. 2002)333Note that there is a typo in the Szkody et al. (2002) in reporting the higher value of the accretion rate used as $`\dot{M}=10^{15}`$ g s-1 rather than $`\dot{M}=10^{16}`$ g s-1. because it is closer to the cooling flow solution, which matches that spectrum quite well (Mukai et al. 2003). However, the cooling model predicts too little emission in the low temperature lines to match the WX Hyi spectrum. The density in the region where the Fe XXII lines are produced is about $`3\times 10^{14}\mathrm{cm}^3`$ in this model, a few times higher than the density sensitive line ratio indicates.

(e) Hot settling flow solutions

Following Medvedev & Narayan (2001), Medvedev & Menou (2002) presented solutions for hot accretion onto unmagnetized, rotating white dwarfs. Together with the accretion rate, the WD rotation rate is an important parameter in these models, because viscously-mediated losses of rotational energy by the WD constitute an additional source of energy powering the X-ray emission from the flow (see Medvedev & Menou 2002 for examples). To calculate the hot flow models we used the numerical relaxation code which solves height-integrated two-temperature hydrodynamic equations with high spatial resolution providing accurate solutions even deep inside the boundary layer. This is the same code which has previously been used by Medvedev & Menou (2002) except for the cooling part which now also incorporates (in addition to bremsstrahlung) the emission line cooling in the form $`q_{\mathrm{line}}^{}=6.6\times 10^{22}T_5^{0.73}n_en_H`$ erg cm-3 s-1, which is a good approximation above $`T=10^5`$ K, where $`T_5=T/(10^5\mathrm{K})`$ (Raymond, Cox & Smith 1976). At the high densities of the WX Hyi boundary layer, quenching of UV lines reduces the cooling rate below $`10^6`$ K, but here we are concerned only with lines formed above $`10^6`$ K. In our present solutions we have chosen the outer boundary of the flow to be at radius $`R=10^2R_{\mathrm{WD}}`$, which is realistic. We have computed hot settling flow models appropriate for WX Hyi. Absent any observational constraint on the WD rotation in this system, we have computed a variety of models with values for the spin parameter, $`s`$ (the angular velocity of the WD in units of the Keplerian angular velocity), ranging from 0.03 to 0.3. The accretion rate in the models was adjusted so that the observed X-ray luminosity is reproduced in each case (see Medvedev & Menou 2002 for details). The required value was on the order of $`10^{15}`$ g s-1.

We found that the X-ray spectrum predicted by the hot settling flow solution greatly overpredicts the X-ray line emission with respect to our observations. Figure 3 (e) shows the spectrum predicted by this solution for the case $`s=0.03`$. As $`s`$ increases, line emission becomes even more pronounced. Therefore, unless optical depth effects become sufficiently important to suppress line emission, this solution does not appear to reproduce the observations well.

## 4 Conclusions

The Chandra observations of the dwarf nova WX Hyi in quiescence show resolved spectral lines produced within the accretion flow. Line diagnostics have allowed us to set constraints on the characteristics of the hot flow, showing that gas densities are very high, $`n10^{13}`$$`10^{14}\mathrm{cm}^3`$, and span a wide range of temperatures, $`T10^6`$$`10^8`$ K.

The spectrum of WX Hyi shows some unusual features if compared to the spectra of the objects in the same class studied by Mukai et al. (2003). Lines in the longer wavelength region of the spectrum are relatively stronger. The O VII line especially cannot be accounted for by the cooling flow model used by Mukai et al. to fit the spectra of this type of objects.

The Chandra data allow the first tests of specific models proposed in the literature for the X-ray emission in quiescent dwarf novae. We have computed the spectra for a set of models ranging from hot boundary layers, to hot settling flows solutions, to X-ray emitting coronae. While most of these models can reproduce well the shape of the continuum, none of them can fully account for the relative line strengths over the entire spectral range. Most coronal models fall short of predicting enough radiation at longer wavelengths, though a thermal conduction-dominated boundary layer predicts too much long wavelength emission. Hot accretion solutions tend to overpredict the short wavelength emission. It is possible, but far from obvious, that some combination of these models, perhaps with thermal conduction present but reduced by a magnetic field, might match the observations.

Given the difficulties of the existing models, it is possible that somewhat different physical processes dominate the boundary layer flow. In particular, it has been suggested that magnetic flux tubes rise above the disk surface and produce X-rays in the same manner as solar flares (Galeev, Rosner & Vaiana 1979). While X-ray emission from the disk as a whole is excluded by the X-ray light curve of OY Car, which limits much of the X-ray emission to the immediate vicinity of the white dwarf (Ramsay et al. 2001), the boundary layer is likely to be magnetically dominated. The boundary layer itself is stable against the magnetorotational instability (Balbus & Hawley 1991), but magnetic flux generated just beyond the boundary layer diffuses into it, and the toroidal component is strongly amplified by the enormous shear (Armitage 2002; Steinacker & Papaloizou 2002). Lacking a more definite prediction, one might expect a power law distribution of X-ray flare energies, as on the Sun (UeNo et al. 1997) and perhaps a significant contribution from flares that reach modest temperatures (e.g. Raymond 1990). Detailed models require Monte Carlo simulations of flare energies, intervals between flares, and physical processes such as thermal conduction, and are beyond the scope of this paper. It is possible, however, that some combination of parameters can be found to roughly match the data.

It is therefore possible that the boundary layer is better described as a collection of impulsive magnetic reconnection events than as a smooth fluid flow. An implication of this idea is that a significant amount of mass may be expelled as a wind. Solar flares are often accompanied by Coronal Mass Ejections (CMEs), and on average the CMEs carry somewhat more kinetic energy than the flares emit as radiation. Thus mass ejections could easily account for the factor of 3 or so discrepancy between the accretion rates estimated from the UV spectrum and the X-ray spectrum. Thus far, winds from cataclysmic variables have only been detected in high $`\dot{M}`$ systems. It has recently been shown that these winds are not driven by radiation pressure alone (Mauche & Raymond 2000; Proga 2003). There is as yet no way to determine whether large scale fields, as in Proga’s models, or small scale CME-like eruptions dominate the wind. Winds from low accretion rate systems are difficult to detect because of the high ionization state expected.

Assuming the broad component of the O VII emission (section 2) is real, it could be formed in the wind hypothesized in the previous paragraph. A simple estimate of the mass loss rate from a surface corresponding to the white dwarf circumference and the boundary layer width ($`10^8`$ cm) with a density $`n3\times 10^{13}\mathrm{cm}^3`$ (as required to match the O VII luminosity) and a speed $`v`$ of 3000 $`\mathrm{km}\mathrm{s}^1`$ yields about $`2\times 10^{15}\mathrm{g}\mathrm{s}^1`$, or more than the total $`\dot{M}`$ inferred from the UV spectra. Thus in order to identify the broad O VII emission with a wind, we require a filling factor of order 0.1. With this filling factor one would in fact need a density $`\sqrt{10}`$ larger to keep the emission measure ($`n^2V`$) of the O VII line at the observed value. The accretion rate $`\dot{M}nvA`$ ($`A`$ being the area) would then be $`\sqrt{10}`$ smaller, consistent with the value inferred from the UV observations. One must be cautious about a wind interpretation both because of the low statistical quality of the O VII profile and because a wide O VII profile might arise from Keplerian rotation of the inner disk just outside the boundary layer. However, a wind of moderate ionization state does provide an appealing explanation for the difference between the accretion rates derived from UV and X-ray observations.

We thank the referee, Chris Mauche, for his careful review of our manuscript and very useful comments. In this research, we have used, and acknowledge with thanks, data from the AAVSO International Database, based on observations submitted to the AAVSO by variable star observers worldwide. This analysis was supported by the CHANDRA grant GO2-3032X to the Smithsonian Astrophysical Observatory and by the Celerity Foundation.


text:
# The Last Hurrah: PPN Formation by a Magnetic Explosion

## 1. Introduction

The formation and shaping of planetary nebulae (PNe) has been described by a generalized interacting stellar wind (GISW) model (Kwok et al., 1978; Balick, 1987), in which the central star of a PN produces a fast wind that interacts with the slower moving, previously ejected, asymptotic giant branch (AGB) wind. According to the GISW model, both the slow and fast winds are driven by radiation, and asymmetric PN shapes result from asymmetries in the slow, AGB wind, from magnetic effects in the shocked fast wind (e.g., Chevalier & Luo, 1994; García-Segura et al., 1999, see also García-Segura in these proceedings), or from some combination of the two.

However, some recent observations of proto-planetary nebulae (pPNe, precursors of PNe) call for a revolutionary change in our understanding of the driving and shaping of post-AGB outflows. These observations are as follows:

1. The material ejected during the pPN phase is highly structured, often containing multiple apparent ejection axes and occasionally appearing quite chaotic (e.g., Balick & Frank, 2002). Several pPNe (e.g., CRL 2688, OH231.8+4.2, IRAS 17150-3224, and IRAS 17441-2411) exhibit quadrupolar symmetry, consisting of outflowing disks—dense enough to obscure central starlight—and bipolar lobes or jets.
2. Spherical shells of reflected starlight often surrounding highly structured PNe and pPNe (e.g., Terzian & Hajian, 2000; Balick et al., 2001, see also Su in these proceedings). This suggests that the AGB wind is spherically symmetric, and so it cannot be responsible for the asymmetry in subsequent flows. Thus, pPN outflows must be self-shaped.
3. The momentum carried in pPN outflows is often a few orders of magnitude larger than can be explained by radiation pressure driving (Knapp, 1986; Bujarrabal et al., 2001, see also Bujarrabal in these proceedings), so a new and more powerful wind driving mechanism is needed.
4. The pPN phase is extremely short lived, typically only a few hundred years. Yet the observed flows often follow a “Hubble law” expansion ($`vR`$, e.g., Alcolea et al., 2001), suggesting that ejecta accelerate in a short time compared to the pPN lifetime, and afterward follow ballistic trajectories. Therefore, the flows are really outbursts, transient winds where the momentum and kinetic energy suddenly and temporarily increase.

Together, these key observations seem to require a new ejection mechanism that is quasi-explosive, is triggered at the end of the AGB phase, and produces intrinsically asymmetric flow structures.

Here we show that the rotational kinetic energy of the core of an AGB star can be converted to linear kinetic energy and transferred to the envelope. The conversion and transferral of energy happens via the twisting of magnetic fields anchored in the core, and naturally results in a complete and asymmetric ejection of the envelope. The acceleration of the envelope occurs in a very short time compared to the typical age of pPNe. We will describe the conditions necessary for such an ejection to take place and discuss how stars may achieve those conditions near the end of the AGB phase. Our model is general and may be relevant to other classes of objects, such as SNe, GRB’s, and $`\eta `$ Carinae.

## 2. The Magnetic Explosion

Consider an AGB star with a slowly contracting core and an expanding, convective envelope. If the star is rotating at all, angular momentum conservation requires that there will be a region of differential rotation at the interface between the core and the envelope. Blackman et al. (2001) have shown that the differential rotation in an AGB star can drive a solar-like magnetic dynamo, in which magnetic fields are amplified to an (time-averaged) equilibrium value of $`10^4`$$`10^6`$ Gauss at the interface. The magnetic field is mixed by convection throughout the envelope. Whether a binary interaction is required to maintain (or trigger) differential rotation in AGB stars remains an open question.

Large mass loss rates from the envelope, during the late stages of the AGB phase, result in a slow expansion of the envelope and a decrease in its density at all depths. As the weight of the overlying envelope decreases, the relative importance of the magnetic field increases. In particular, there may be a threshold at which the magnetic (plus thermal) energy within the envelope exceeds its gravitational potential energy. At that point, the envelope as a whole would be driven off in a single, short-lived event.

Due to the extreme difficulty of modeling, self-consistently, the complex interface dynamo, together with a proper treatment of the turbulent envelope and mass loss from the top of the convection zone, we consider here a slightly different and simplified situation. Rather than having an envelope with a constant magnetic energy and a decreasing density, we consider a constant density envelope in which the magnetic energy increases. Due to our simplifications, the magnetic energy is generated in a short time and only at the interface between the core and envelope, where strong differential rotation twists up pre-existing poloidal magnetic fields. Though the magnetic energy is produced only at the base of the envelope (and not mixed everywhere throughout), we expect the same qualitative behavior described above, where the magnetic energy exceeds the gravitational potential energy, and the envelope is driven off in an outburst.

For this initial study, we have carried out time-dependent, numerical magnetohydrodynamic (MHD) simulations. We use the simulation code of Matt (2002), which employs a finite difference scheme to solve the standard, adiabatic ($`\gamma =5/3`$), ideal MHD equations, plus source terms for gravity, on an Eulerian, cylindrical ($`r`$, $`z`$), nested grid. The code assumes axisymmetry (where $`\delta /\delta \varphi `$ of all quantities $`=0`$). We use “outflow” boundary conditions (BCs) on the outermost boundary (though material doesn’t leave the grid by the end of the simulations), and assume reflection symmetry about the equator ($`z=0`$). The simulation domain, then, consists of an $`r`$$`z`$ slice through a single quadrant, with the stellar core specified by a circular, inner BC.

For simplicity, we want to ignore effects of a wind from the core, so we have chosen a “solid” BC such that no matter can flow from or onto the core. The core behaves as a perfect electric conductor and initially is rotating as a solid body and is threaded by a rotation-axis aligned dipole magnetic field. The computational domain (existing between the core boundary and the outer boundary) is initialized with a hydrostatic envelope whose density falls off as $`R^2`$ (where $`R`$ is the spherical radius). The envelope is initially stagnant (no rotation or motion) and is threaded everywhere by the dipole magnetic field that is anchored to the core. We also assume that the mass of the core is much greater than that of the envelope, so that the envelope’s self-gravity can be ignored.

The upper left panel of figure 1 shows the initial configuration in the inner part of the computational domain (from the center, out to 13 core radii, $`R_{\mathrm{core}}`$). The initial configuration is static and in equilibrium. At $`t=0`$ the core begins to rotate. Differential rotation between the core and stagnant envelope leads to an azimuthal twisting of the dipole field lines just above the core. The azimuthal magnetic field $`B_\varphi `$ increases linearly in time, until it becomes strong enough to be dynamically important. At that time, magnetic pressure forces (i.e., $`B_\varphi ^2`$), begin to drive material outward, away from the core.

In our chosen configuration, the key parameters are the relative values of the characteristic speeds in the system—the equatorial rotation speed of the core ($`v_{\mathrm{rot}}`$), the initial Alfvén speed at the core equator ($`v_\mathrm{A}`$), and the escape speed from the surface of the core ($`v_{\mathrm{esc}}`$)—which are related to the relative values of various energies in the system. Note that in hydrostatic equilibrium, the sound speed is proportional to the escape speed, so the sound speed (or thermal energy) is not a key, independent parameter. Here, we present a case with $`v_{\mathrm{rot}}/v_{\mathrm{esc}}=0.1`$ (14% of breakup rotation) and $`v_\mathrm{A}/v_{\mathrm{esc}}=1.0`$.

The upper right, lower left, and lower right panels of figure 1 show the evolution after 1.5, 3.0, and 4.5 rotations of the core, respectively. The expansion of newly generated azimuthal magnetic field is evident in the figure, as it drives the envelope material outward. This also leads to an expansion of the poloidal magnetic field (lines), though some field lines near the equator remain closed, where material is forced and held in corotation with the core. Due to the coupling of the rotation to the dipole field, $`B_\varphi `$ has a maximum near the core at mid latitudes. Thus, magnetic pressure forces point (from strong to weak field), not only radially away from the core, but also toward the rotation axis and toward the magnetic equator. This produces structures in the outburst such that the expansion is faster toward the equatorial and polar directions. Finally, notice that the density near the core decreases with time. This is because that region is drained of material by the expansion of the magnetic field (our “solid” core BC prevents the region from being replenished by mass flow from the core). This process is truly a transient phenomenon.

The four panels of figure 2 show the evolution of the system at the same times and for the same simulation as in figure 1. However, the data in figure 2 is shown on a scale that is 10 times larger than in figure 1 (by taking advantage of nested simulation grids), and the data has been reflected about the rotation axis and equator to better illustrate the flow.

As the simulation proceeds, one sees a swept up shell of envelope material moving outward from the core. In the upper right panel, the outer edge of this shell is quite spherical, and it represents the initial pressure wave expanding outward to make room for the expanding magnetic field. As the field expansion accelerates and exceeds the sound speed, it overtakes the spherical pressure wave and begins to protrude outward along the rotation axis and magnetic equator. This is evident in the bottom two panels, where the outer dense region is swept up envelope material, and the “hollow” inner region is actually filled with magnetic (mostly $`B_\varphi `$) energy. Note the strongly quadrupolar shape of this inflating bubble. Again, the asymmetry is due to the fact that a spinning dipole produces a $`B_\varphi `$ that is strongest at mid latitudes, so magnetic pressure forces direct outward from the core, and also toward the pole and equator. The magnetic field responsible for inflating the bubble is generated by the rotation of the core. Thus we are seeing the spin energy of the core being converted into linear kinetic energy of the envelope (and the core will spin down, as a result).

We have run several simulations with different relative values of $`v_{\mathrm{rot}}`$, $`v_\mathrm{A}`$, and $`v_{\mathrm{esc}}`$. We find that the envelope will be fully ejected from the system, as long as the dimensionless parameter $`v_{\mathrm{rot}}v_\mathrm{A}v_{\mathrm{esc}}^2>0.1`$. For smaller values, the outburst may occur only at selected latitudes (e.g., the pole and equator), or, for very small values, there is no envelope ejection. For cases where the envelope is ejected, the expansion speed of the swept up shell is approximately twice the rotation speed, $`v_{\mathrm{rot}}`$. Also, various shell shapes are produced, with different relative speeds in the polar and equatorial regions.

The enhanced equatorial flow is a feature of the spinning dipole field. Different field geometries will not necessarily have an equatorial enhancement, though there will always be a polar enhancement. To demonstrate this more clearly, we ran a simulation with the same parameters as the data in figures 1 and 2, but with a split monopole field geometry (radial magnetic field, but with a direction reversal at the equator). The right panel of figure 3 shows the split monopole case after three rotations of the core. For comparison, the left panel shows the dipole case (same as the lower right panel of fig. 2).

## 3. Discussion

The magnetic explosion mechanism presented here, if triggered at the end of the AGB phase, may solve all of the problems discussed in section 1. Triggering requires that the field and the rotation of the core must either be sustained throughout the lifetime of the AGB phase or be generated near the end of the AGB phase, perhaps stimulated by binary interaction. To compare our simulations to real pPNe, we can scale our parameters to real values. For example, we can assume a core mass of $`0.5M_{}`$, a swept up shell mass and size of $`0.062M_{}`$ and $`10^4`$ AU (for CRL 2688, Bujarrabal et al., 2001, see also Kastner in these proceedings), and $`v_{\mathrm{esc}}=1000`$ km s-1. In that case, the simulation presented in figures 1 and 2 applies to a core with a radius of $`0.2R_{}`$, $`v_{\mathrm{rot}}=100`$ km s-1, and a dipole field strength of $`2\times 10^5`$ Gauss. If instead we assume $`v_{\mathrm{esc}}=3000`$ km s-1, the same simulation applies to a core with a radius of $`10^9`$ cm, $`v_{\mathrm{rot}}=300`$ km s-1, and a dipole field strength of $`10^8`$ Gauss.

The time taken to extract all of the rotational kinetic energy from the core (i.e., the spin-down time) is $`100`$ years, comparable to the typical pPNe lifetime. This may explain the short acceleration times apparent in pPNe outflows and the relative lack of quickly rotating white dwarfs (Koester et al., 1998).

### Acknowledgments.

We thank the organizers for a productive and enjoyable meeting. This research was supported by NSERC, McMaster University, and CITA through a CITA National Fellowship.