text:
# Entanglement capabilities of non–local Hamiltonians

## Abstract

We quantify the capability of creating entanglement for a general physical interaction acting on two qubits. We give a procedure for optimizing the generation of entanglement. We also show that a Hamiltonian can create more entanglement if one uses auxiliary systems.

In the last 2-3 years there has been very considerable increase in experimental activity aiming to create entangled quantum states. One reason is the potential applications of entanglement to quantum information processing. Creating entanglement has been possible in quantum optics for more than a decade, however now many new communities, working in a variety of experimental areas (for example NMR, condensed matter physics) are also joining the field . In general, entanglement between two systems can be generated if they interact in a controlled way. However, in most experiments these interactions are weak which make the production of entanglement a very difficult task. Thus, it would be very convenient to have a theory which would provide us with the best way of exploiting interactions to produce entanglement.

In this Letter we analyze the entanglement capabilities of Hamiltonians. In particular, we would like to answer questions like: given an interaction (Hamiltonian), what is the most efficient way of entangling particles? Can we make the process more efficient by supplementing the action of the Hamiltonian with some local unitary operations? Can we increase the entanglement more efficiently by using some ancillas?

So far, much of the theoretical effort in quantum information theory has been devoted to the characterization and quantification of the entanglement of a given state. Very recently, it has been realised that there is a parallel notion of the entanglement in the dynamics of a system . In , the authors consider the situation that one has a given unitary transformation and ask, for example, how much state entanglement is needed to produce it. Here we focus on a different issue: given an interaction (i.e. a Hamiltonian) how can we make the most effective use of it. What we propose here is to define and determine the entanglement capabilities of physical processes, in particular, of unitary evolutions . This is a very relevant problem not only from the theoretical point of view, but also from the experimental one. Of course, this problem is even more difficult than the one of quantifying the entanglement of states. In any case, in the present work we give the first steps in this direction by considering the case in which the physical process is acting on two qubits.

From our results it turns out that: (i) It is more efficient to produce entanglement if one initially has already some; (ii) The best initial entanglement is universal, i.e. independent of the physical process; (iii) One can improve the performance of a physical process by complementing it with fast local operations; (iv) One can also improve it (in certain cases) by using auxiliary systems; (v) All entangling Hamiltonians can simulate each other and are thus qualitatively equivalent; we also provide an upper bound on the time required for one Hamiltonian to simulate another.

We consider two qubits interacting via a non–local Hamiltonian $`H`$. We want to determine the most efficient way in which we can use such an interaction to produce entanglement. We will characterize the entanglement of a state of the qubits at a given time $`t`$, $`|\mathrm{\Psi }(t)`$, by some entanglement measure $`E`$. In order to quantify the entanglement production, we define the entanglement rate $`\mathrm{\Gamma }`$ at a particular time $`t`$ of the interaction as follows:

$$\mathrm{\Gamma }(t)\frac{dE(t)}{dt}.$$
(1)
This quantity depends on $`|\mathrm{\Psi }(t)`$ not only through its entanglement $`E`$. The goal is then to find the conditions which must be satisfied in order to obtain a maximal entanglement rate. In particular, we will be interested in determining the following:

For any initial entanglement $`E`$ of the two-qubit system, what is the state $`|\mathrm{\Psi }`$, say $`|\mathrm{\Psi }_E`$, for which the interaction produces the maximal rate $`\mathrm{\Gamma }_E`$.

The maximal achievable entanglement rate $`\mathrm{\Gamma }_{\mathrm{max}}`$ , $`\mathrm{\Gamma }_{\mathrm{max}}\mathrm{max}_E\mathrm{\Gamma }_E`$ and the state $`|\mathrm{\Psi }_{\mathrm{max}}`$ for which $`\mathrm{\Gamma }=\mathrm{\Gamma }_{\mathrm{max}}`$.

These quantities are interesting because the knowledge of the state $`|\mathrm{\Psi }_E`$ will allow us to find out the most efficient way of entangling the qubits. The idea is to supplement the interaction Hamiltonian $`H`$ with appropriate local unitary operations in such a way that the state of the qubits at any time $`t`$ is precisely $`|\mathrm{\Psi }_{E(t)}`$, for which the increase of entanglement is optimal. In order to show how this can be achieved, let us consider that the evolution given by $`H`$ proceeds in very small time steps $`\delta t`$. Let us also assume that the qubits are initially disentangled. Using local operations, we can always prepare the state $`|\mathrm{\Psi }_0`$ —that is, the product state which most efficiently becomes entangled under the action of $`H`$. After a time step $`\delta t`$, the state will change and its entanglement will increase to $`\delta E`$. Then, we use (fast) local unitary operations to transform the new state of the qubits into the state $`|\mathrm{\Psi }_{\delta E}`$ for which $`\mathrm{\Gamma }`$ is optimal. Note that this is always possible, since for qubits all states with the same value of $`E`$, say $`\delta E`$, are connected by local unitary transformations. By proceeding in the same way after every time step, and taking the continuous time limit $`\delta t0`$ we obtain that the state of the qubits at time $`t`$ is always the optimal one, $`|\mathrm{\Psi }_{E(t)}`$. Obviously, in an experimental realization, this procedure requires that we can apply the appropriate local transformations in times which are short compared to the typical time scale $`\tau _H`$ associated to $`H`$, $`\tau _H=(e_{\mathrm{max}}e_{\mathrm{min}})^1,`$ where $`e_{\mathrm{max}}`$ and $`e_{\mathrm{min}}`$ are the maximum and minimum eigenvalues of $`H`$, and we have set $`\mathrm{}=1`$.

Knowledge of $`\mathrm{\Gamma }_E`$ also permits us to determine the maximum amount of entanglement $`E_{\mathrm{max}}`$ produced as a function of time. We just have to express $`\mathrm{\Gamma }_E`$ as an explicit function of $`E`$, substitute it in (1) and solve that differential equation to determine $`E_{\mathrm{max}}(t)`$. Note that the optimal procedure described above will precisely reach the entanglement $`E_{\mathrm{max}}(t)`$.

The state $`|\mathrm{\Psi }_{\mathrm{max}}`$ is important since it gives rise to the maximal increase of entanglement, and therefore corresponds to the best operational point. After reaching the state $`|\mathrm{\Psi }_{\mathrm{max}}`$ with the procedure described in the previous paragraph, the entanglement would be produced in a very efficient way, if one could transfer the entanglement that is gained after each time step $`\delta t`$ to other qubits (using entanglement dilution  or some other means). In particular, it would increase proportionally to the time, $`\mathrm{\Gamma }_{\mathrm{max}}`$ being the proportionality constant.

In the following, we will show how to determine $`|\mathrm{\Psi }_E`$, $`\mathrm{\Gamma }_E`$, $`|\mathrm{\Psi }_{\mathrm{max}}`$, and $`\mathrm{\Gamma }_{\mathrm{max}}`$ for an arbitrary Hamiltonian $`H`$. To this end, it is convenient to use the Schmidt decomposition of the state of the qubits $`|\mathrm{\Psi }(t)`$ to write

$$|\mathrm{\Psi }=\sqrt{P}|\phi ,\chi +e^{i\alpha }\sqrt{1P}|\phi ^{},\chi ^{},$$
(2)
where for the sake of short–hand notation we have omitted the time dependence of all these quantities. Here, $`\phi |\phi ^{}=\chi |\chi ^{}=0`$ and $`P1/2`$. Note that $`E`$ must only depend on the Schmidt coefficient $`P`$, given the fact that it must be invariant under local unitary operations. For example, if we choose as entanglement measure the entropy of entanglement  — the entropy of the reduced density operator of one of the qubits—, we will have

$$E(P)=P\mathrm{log}_2(P)(1P)\mathrm{log}_2(1P).$$
(3)
Note that the entropy of entanglement quantifies the amount of EPR entanglement contained asymptotically in a pure state $`|\mathrm{\Psi }`$. That is $`E(P)`$ gives the ratio of maximally entangled EPR states $`|\mathrm{\Psi }^{}=1/\sqrt{2}(|01|10)`$ which can be distilled from \[are needed to create\] $`|\mathrm{\Psi }`$ respectively. Thus, we can write

$$\mathrm{\Gamma }(t)=\frac{dE}{dP}\frac{dP}{dt}.$$
(4)
In (4), given a particular entanglement measure $`E(P)`$, we just have to determine $`dP/dt`$. In order to do that, we need to find the (infinitesimal) time evolution of the Schmidt coefficients of the state of the qubits. After a time $`\delta t`$ we will have $`|\mathrm{\Psi }(t+\delta t)=\mathrm{exp}(iH\delta t)|\mathrm{\Psi }(t)(1iH\delta t)|\mathrm{\Psi }(t)`$. The corresponding reduced density operator $`\rho _A(t+\delta t)`$, where $`\rho _{A,B}=\mathrm{Tr}_{B,A}(|\mathrm{\Psi }\mathrm{\Psi }|)`$, can then be written as $`\rho _A(t+\delta t)=\rho _A(t)i\delta t\mathrm{Tr}_B\{[H,|\mathrm{\Psi }(t)\mathrm{\Psi }(t)|]\}.`$ The eigenvalues (Schmidt coefficients) of this operator can be easily determined starting from $`\rho _A|\phi =P|\phi `$ and using standard perturbation theory. We find

$$\frac{dP}{dt}=2\sqrt{P(1P)}\times \mathrm{Im}[e^{i\alpha }\phi ,\chi |H|\phi ^{},\chi ^{}],$$
(5)
where we have omitted the time–dependence. Upon substitution in (4) we obtain the entanglement rate. Since we are interested in maximizing $`\mathrm{\Gamma }`$, it is clear that we can always choose $`\alpha `$ such that

$$\mathrm{\Gamma }=f(P)|h(H,\phi ,\chi )|.$$
(6)
where

$`f(P)`$ $`=`$ $`2\sqrt{P(1P)}E^{}(P),`$ (8)
$`h(H,\phi ,\chi )`$ $`=`$ $`\phi ,\chi |H|\phi ^{},\chi ^{}.`$ (9)
By analyzing Eq. (6) we can extract some interesting conclusions, even before determining the maximum value of $`\mathrm{\Gamma }`$ explicitly. Given the fact that $`f`$ and $`h`$ depend on different parameters, in order to determine the quantities mentioned in (i–ii) we can maximize the functions $`f`$ and $`|h|`$ independently. First, if we want to determine the quantities mentioned in (i), we have to fix the value of $`E`$. In that case, $`P`$ is also fixed and therefore the maximum of the entanglement rate will correspond to a state of the form (2) with some fixed $`|\phi `$, $`|\chi `$, and $`\alpha `$ (which maximize $`|h|`$). That is, for any value $`E`$ of the entanglement, the states $`|\phi `$ and $`|\chi `$ for which the maximal entanglement rate $`\mathrm{\Gamma }_E`$ is obtained do not depend on $`E`$, but only on the form of the Hamiltonian $`H`$. Let us denote by $`h_{\mathrm{max}}`$ the maximum value of $`|h|`$; that is,

$$h_{\mathrm{max}}=\underset{\phi ,\chi =1}{\mathrm{max}}|\phi ,\chi |H|\phi ^{},\chi ^{}|.$$
(10)
Then, we can easily determine how the entanglement would evolve with time if we always drive the qubits with local operations so that at each time their state corresponds to the optimal one. We can simply solve the differential equation (5), obtaining $`P(t)=\mathrm{sin}^2[h_{\mathrm{max}}t+\varphi _0],`$ with $`P(0)=\mathrm{sin}^2(\varphi _0)`$. Using the explicit dependence of $`E`$ on $`P`$, we can directly then calculate $`E(t)`$. The evolution of the entanglement is fully characterized by $`h_{\mathrm{max}}`$, which is a quantity that only depends on the interaction Hamiltonian. That is, for a given $`H`$, $`h_{\mathrm{max}}`$ measures the capability of creating entanglement. In the following we give a simple way of determining $`h_{\mathrm{max}}`$, which allows us to classify the entanglement capability of any Hamiltonian. On the other hand, once the entanglement measure is specified, we can calculate the value $`P_0`$ of $`P`$ for which we obtain the maximal rate by simply considering the function $`f(P)`$. For example, choosing the expression (3) for the entanglement, we find that $`P_0`$ solves the equation $`\mathrm{ln}\frac{1P_0}{P_0}=\frac{2}{12P_0},`$ i.e. $`P_00.0832`$ which gives $`E(P_0)0.413`$. This shows that, in order to increase the entanglement of a two-qubit system in an optimal way, it is better to start with some initially entangled state rather than a product state . Note that the optimal initial entanglement $`E(P_0)`$ is independent of $`H`$.

In the following we will show how to determine the entanglement capability $`h_{\mathrm{max}}`$ of a general Hamiltonian $`H`$ acting on the qubits. First, we will show how, if we allow to supplement the evolution of $`H`$ by local unitary operations, we can express $`H`$ in a standard form that only depends on three parameters. Then we will derive an expression for $`h_{\mathrm{max}}`$ in terms of those parameters.

Except for a trivial constant, we can always express a general Hamiltonian as

$`H`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{3}{}}}\alpha _i\sigma _i^A1𝐥_B+{\displaystyle \underset{j=1}{\overset{3}{}}}\beta _j1𝐥_A\sigma _j^B+`$ (12)
$`{\displaystyle \underset{i,j=1}{\overset{3}{}}}\gamma _{i,j}\sigma _i^A\sigma _j^B.`$
Here, $`\sigma _i`$ are the Pauli operators, and $`\stackrel{}{\alpha }`$, $`\stackrel{}{\beta }`$, and $`\gamma `$ are two real vectors and a real matrix, respectively. We now show that by supplementing the evolution operator with local unitary operations we can obtain an effective Hamiltonian which has one of the two standard forms

$$\widehat{H}^\pm =\mu _1\sigma _1^A\sigma _1^B\pm \mu _2\sigma _2^A\sigma _2^B+\mu _3\sigma _3^A\sigma _3^B,$$
(13)
where $`\mu _1\mu _2\mu _30`$ are the (sorted) singular values of the matrix $`\gamma `$ . We first note that the terms corresponding to $`\stackrel{}{\alpha },\stackrel{}{\beta }`$ in (12) give no contribution to $`h_{\mathrm{max}}`$ (10) and can therefore be neglected. Second, we apply the local operations $`U`$ ($`V`$) and $`U^{}`$ ($`V^{}`$) to the first (second) qubit at the beginning and end of the evolution process, respectively. We select them such that $`U^{}\sigma _i^AU=_{k=1}^3O_{k,i}^A\sigma _k^A,V^{}\sigma _j^BV=_{l=1}^3O_{j,l}^B\sigma _l^B,`$ where $`O^{A,B}`$ are orthogonal matrices of determinant one, each being plus or minus the orthogonal matrices in a singular value decomposition of $`\gamma `$. Thus the total (non-local) effect of the evolution for a time $`t`$ is equivalent to the one obtained with the Hamiltonian $`\widehat{H}^+`$ ($`\widehat{H}^{}`$) for the same time if det$`\gamma 0`$ (det$`\gamma <0`$). Without loss of generality, we may take $`H`$ of the form $`\widehat{H}^+`$ (13) in what follows .

Now let us determine $`h_{\mathrm{max}}`$ in terms of $`\mu _{1,2,3}`$. We can write $`h(H,\phi ,\chi )=_{k=1}^3\mu _k\phi |\sigma _k^A|\phi ^{}\chi |\sigma _k^B|\chi ^{}.`$ Using the Cauchy–Schwarz inequality, it can be checked that the maximum of (the absolute value of) this function is reached for $`|\chi =|\phi ^{}`$. In this case, using the fact that $`|\phi \phi |+|\phi ^{}\phi ^{}|=1𝐥`$ we obtain $`h(H,\phi ,\phi )=_{k=1}^3\mu _k_{k=1}^3\mu _k\phi |\sigma _k|\phi ^2.`$ Taking into account that $`\mu _1\mu _2\mu _3`$, we see that the maximum value occurs when $`|\phi =|0`$ or $`|\phi =|1`$, i.e. an eigenstate of $`\sigma _3`$. For that choice we obtain

$$h_{\mathrm{max}}=\mu _1+\mu _2.$$
(14)
Summarizing, once we have transformed the Hamiltonian $`H`$ to the standard form (13) we obtain that for a given value of $`E`$ (and therefore of $`P`$),

$`|\mathrm{\Psi }_E`$ $`=`$ $`\sqrt{P}|0,1+i\sqrt{1P}|1,0,`$ (16)
$`\mathrm{\Gamma }_E`$ $`=`$ $`f(P)h_{\mathrm{max}}`$ (17)
where $`h_{\mathrm{max}}=\mu _1+\mu _2`$. The maximum rate $`\mathrm{\Gamma }_{\mathrm{max}}`$ is obtained for $`P=P_0`$, where $`P_0`$ is the value that maximizes $`f(P)`$. Thus, $`|\mathrm{\Psi }_{\mathrm{max}}`$ and $`\mathrm{\Gamma }_{\mathrm{max}}`$ are given by (14) with $`P=P_0`$. For example, for the entanglement measure (3), $`P_00.0832`$ which leads to $`f(P_0)1.9123`$.

So far, we have calculated the most efficient way of entangling two qubits if we can use local unitary operations acting on each of the qubits. We have not allowed, however, local operations which entangle each of the qubits with local ancillas. We will now show that this possibility permits us to increase the maximum entanglement rate $`\mathrm{\Gamma }_{\mathrm{max}}`$ for certain kind of Hamiltonians. We will first generalize the formulas derived above to the case of multilevel systems, given that the system qubit–plus–ancilla is of this sort. We consider a state $`|\mathrm{\Psi }`$ with Schmidt decomposition $`|\mathrm{\Psi }=_{n=1}^N\sqrt{\lambda _n}|\phi _n,\chi _n.`$ As before, any entanglement measure $`E`$ will only depend on the Schmidt coefficients $`\lambda _n0`$. In particular, in the following we will use the entropy of entanglement, $`E(\stackrel{}{\lambda })=_{n=1}^N\lambda _n\mathrm{log}_2(\lambda _n).`$ Using the definition (1) of entanglement rate, we have

$$\stackrel{~}{\mathrm{\Gamma }}=\underset{n=1}{\overset{N}{}}\frac{E}{\lambda _n}\frac{d\lambda _n}{dt}=\frac{1}{N}\underset{n,m=1}{\overset{N}{}}\left[\frac{E}{\lambda _n}\frac{E}{\lambda _m}\right]\frac{d\lambda _n}{dt},$$
(18)
where we have used the fact that the sum of all the Schmidt coefficients is constant. Using perturbation theory as before, we find $`\frac{d\lambda _n}{dt}=2\sqrt{\lambda _n}\mathrm{Im}[\phi _n,\chi _n|H|\mathrm{\Psi }]=2_{m=1}^N\sqrt{\lambda _n\lambda _m}\mathrm{Im}[\phi _n,\chi _n|H|\phi _m,\chi _m].`$

Rather than proceeding in complete generality we now consider an example which demonstrates that adding ancillas may allow one to increase entanglement more efficiently than is possible without the use of ancillas. We will consider the case in which the ancillas are also qubits. We write $`P=\lambda _1`$ and concentrate on the case in which $`\lambda _2=\lambda _3=\lambda _4=(1P)/3`$. In that case, Eq. (18) simplifies to

$$\stackrel{~}{\mathrm{\Gamma }}=\stackrel{~}{f}(P)\stackrel{~}{h}(H,\phi _n,\chi _n)$$
(19)
where now

$`\stackrel{~}{f}(P)=2\sqrt{P(1P)/3}\mathrm{log}_2[(1P)/(3P)],`$ (21)
$`\stackrel{~}{h}(H,\phi _n,\chi _n)={\displaystyle \underset{n=2}{\overset{4}{}}}\mathrm{Im}[\phi _1,\chi _1|H|\phi _n,\chi _n].`$ (22)
We can always choose the phase of the states $`|\phi _n`$ such that all the terms on the sum add with the same sign. We can therefore replace the imaginary parts of the terms in the above expression by their absolute values, and in (19) we can replace $`\stackrel{~}{f}(P)`$ by $`|\stackrel{~}{f}(P)|`$. We find that $`\stackrel{~}{P}_00.8515`$ (which corresponds to an entropy of entanglement $`E(\stackrel{~}{P}_0)0.8415`$) maximizes $`|\stackrel{~}{f}(P)|`$ (21) and leads to $`|\stackrel{~}{f}(\stackrel{~}{P}_0)|1.6853`$. Proceeding as before, we can easily maximize $`\stackrel{~}{h}`$. We obtain that the maximum value is $`\stackrel{~}{h}_{\mathrm{max}}=\mu _1+\mu _2+\mu _3,`$ which occurs when $`|\phi _n=|\chi _n`$ are orthogonal maximally entangled states between the qubit and the ancilla. For example in the case that $`det\gamma 0`$, the choice $`|\phi _1=|\varphi ^+,|\phi _2=i^{\frac{3}{2}}|\psi ^+,|\phi _3=i^{\frac{1}{2}}|\psi ^{},|\phi _4=i^{\frac{3}{2}}|\varphi ^{},`$ where $`\{|\varphi ^\pm ,|\psi ^\pm \}`$ are Bell states , together with $`P=P_0=0.8515`$ leads to a maximal \[under the previous assumptions on the $`\lambda _i`$’s\] entanglement rate $`\stackrel{~}{\mathrm{\Gamma }}=\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{max}}`$.

Let us compare the cases in which we use ancillas and the one in which we do not use them. On the one hand, we have that $`|\stackrel{~}{f}(\stackrel{~}{P}_0)|<|f(P_0)|`$. But on the other, $`\stackrel{~}{h}_{\mathrm{max}}h_{\mathrm{max}}`$. Thus, if $`\mu _30`$ it may be the case that the use of ancillas can help to increase the maximum rate of entanglement $`\mathrm{\Gamma }_{\mathrm{max}}`$ as well as the rate $`\mathrm{\Gamma }_E`$ for a given entanglement $`E`$ of state $`|\mathrm{\Psi }`$. This is in fact the case, if we have, for example, $`\mu _1=\mu _2=\mu _3`$ (i.e. $`\stackrel{~}{h}_{\mathrm{max}}=3/2h_{\mathrm{max}}`$). In this case we obtain $`\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{max}}1.3220\mathrm{\Gamma }_{\mathrm{max}}`$. In a similar way, one can check for this specific Hamiltonian that $`\stackrel{~}{\mathrm{\Gamma }}_E\mathrm{\Gamma }_E`$ if the initial entanglement satisfies $`E0.08`$.

Finally, it is easy to show that all entangling Hamiltonians are qualitatively equivalent when assisted by local operations. In particular, given two Hamiltonians $`H`$ and $`H^{}`$ with either $`h_{\mathrm{max}}=\alpha h_{\mathrm{max}}^{}`$ or $`\stackrel{~}{h}_{\mathrm{max}}=\alpha \stackrel{~}{h}_{\mathrm{max}}^{}`$, one can simulate the action of $`H^{}`$ for any time $`t`$ by applying $`H`$ for at most $`3\alpha ^1t`$. This can be seen as follows: Applying $`H`$ of the form (13) for $`\delta t/2`$ followed by a local unitary operation $`\sigma _1`$ in $`A`$ before and after another application of $`H`$ for $`\delta t/2`$ is equivalent to the application of the Hamiltonian $`H_1=\mu _1\sigma _1\sigma _1`$ for the time $`\delta t`$, provided that $`\delta t`$ is infinitessimally small. Since $`H_1`$ is locally equivalent to $`H_k=\mu _1\sigma _k\sigma _k`$, applying sequentially $`H_k`$ for $`\delta t\mu _k^{}/\mu _1`$ is equivalent to the application of $`H^{}`$ for the time $`\delta t`$. Using the restrictions on $`\mu _k,\mu _k^{}`$ the claim readily follows. The question of efficient simulation of another Hamiltonian in the same time $`t`$ will be addressed in future work.

In summary, we have found the optimal way of using any non–local interaction to entangle a pair of qubits. The idea is to use local operators to drive the instantaneous state to the one that maximizes the entanglement rate, at each moment of the evolution. We have found that the entanglement capacity of any given Hamiltonian is determined by the sum of the two largest singular values of the matrix $`\gamma `$ defined in (12). Finally, we have shown that for certain Hamiltonians one can overcome this maximal entanglement rate by using ancillas prepared in maximally entangled states with the qubits.

We are grateful to C. H. Bennett for very valuable input, particularly concerning the issue of the qualitative equivalence of Hamiltonians. This work was supported by the Austrian SF, the European Community (TMR network ERB–FMRX–CT96–0087 ; project EQUIP (contract IST-1999-11053)), the ESF, and the Institute for Quantum Information GmbH. G.V also acknowledges funding from EC through grant HPMF-CT-1999-00200.

Note added: The idea of universal simulation of Hamiltonians as discussed in this manuscript originated in a larger collaboration, including several other people and is developed in detail in . After completion of this manuscript we have also learned that the notion of universal simulation of Hamiltonians has been independetly addressed by Dodd et al. .


format: text

text:
# A spiral structure in the disk of EX Draconis on the rise to outburst maximum

## 1 Introduction

Dwarf novae are close binaries in which a late type star (the secondary) overfills its Roche lobe and transfers matter to a companion white dwarf via an accretion disk. They show recurrent outbursts of 2–5 magnitudes on timescales of days to several months, resulting either from an instability in the mass transfer from the secondary (the MTI model), or from a thermal instability in the accretion disk (the DI model) which switches the disk from a low to a high-viscosity regime (Warner 1995 and references therein). Observations of dwarf novae at the early phases of the outburst offer good prospects to discriminate between the two competing models. The MTI model predicts that the accretion disk should shrink at the beginning of the outburst due to the sudden accretion of material with low angular momentum while the luminosity of the bright spot should increase in response to the higher mass-transfer rate (e.g. Warner 1995). No such effects are expected in the DI model. Most of the few current observations seem to support the DI model (e.g., Vogt 1983; Rutten et al. 1992a), although recent results suggest an MTI origin for the outbursts in the intermediate polar EX Hya (Hellier et al. 2000).

Spiral shocks have long been proposed as a possible mechanism for transport of angular momentum in accretion disks (e.g., Sawada, Matsuda & Hachisu 1986; Savonije, Papaloizou & Lin 1994) and may be the key, together with magnetic viscosity (Hawley, Balbus & Winters, 1999), in understanding the anomalous viscosity mechanism responsible for the spiraling inwards of the disk material. The recent discovery of spiral shocks in the accretion disk of the dwarf nova IP Pegasi in outburst (e.g., Steeghs, Harlaftis & Horne 1997; Baptista, Harlaftis & Steeghs 2000a) confirmed the results of hydrodynamical simulations (e.g., Stehle 1999) and boosted the research in this topic. Two-armed spiral shocks are excited in the outer regions of the disk by the tides raised by the secondary star when the outbursting accretion disk extends far enough into the Roche Lobe.

EX Draconis is a long period ($`P_{orb}=5`$ hr) eclipsing dwarf nova, the outbursts of which have typical amplitudes of 2.0 mag, duration of $`10`$ days, and recurrence timescale of about 20 days (Baptista, Catalán & Costa 2000b; hereafter BCC). The analysis of a set of light curves along the outburst cycle led BCC to conclude that the outbursts do not start in the outer disk regions. (i.e., they are inside-out outbursts). Recently, Joergens, Spruit & Rutten (2000) found evidence of spiral shocks in the accretion disk of EX Dra from Doppler tomography close to outburst maximum.

In this letter we report on the comparative analysis of the light curves of BCC during the early rise to outburst maximum with those of the quiescent state of the binary. Our results show that a one-armed spiral structure develops in the accretion disk of EX Dra on the rise to maximum. Section 2 describes the data analysis. The results are presented in section 3 and discussed in section 4.

## 2 Data Analysis

Our data are $`R`$ band light curves obtained with the 0.9-m James Gregory Telescope at the University of St. Andrews in 1995. The reader is referred to BCC for a detailed description of the data set and of the reduction procedures. We collected three eclipse light curves of EX Dra on the rising branch of the 1995 November outburst (eclipse cycles 7812, 7813 and 7817), two days after the onset of the outburst and two days before the outburst maximum (outburst D of BCC; see their Fig. 1). The quiescent state is represented by four light curves (eclipse cycles 8007, 8008, 8011 and 8013), obtained on 1995 December 27-28, about 10 days after the end of the subsequent outburst (outburst E of BCC). In order to improve the signal-to-noise ratio and to reduce the influence of flickering, the individual light curves were combined to produce average eclipse light curves for the two observed states. For each light curve, we divided the data into phase bins of 0.003 cycle and computed the median for each bin. The median of the absolute deviations with respect to the median for each bin was taken as the corresponding uncertainty. The light curves were phase folded according to the sinusoidal ephemeris of BCC,

$$T_{mid}=\mathrm{HJD}\mathrm{\hspace{0.33em}2\hspace{0.17em}448\hspace{0.17em}398.4530}(\pm 1)+\mathrm{0.209\hspace{0.17em}936\hspace{0.17em}98}(\pm 4)E+$$

$$+(8.2\pm 1.5)\times 10^4\mathrm{sin}\left[2\pi \frac{(E968)}{7045}\right]d.$$
(1)
Out-of-eclipse brightness changes are not accounted for by the eclipse mapping method, which assumes that all variations in the eclipse light curve are due to the changing occultation of the emitting region by the secondary star. Orbital variations were therefore removed from the light curves by fitting a spline function to the phases outside eclipse, dividing the light curve by the fitted spline, and scaling the result to the spline function value at phase zero. This procedure removes orbital modulations with only minor effects on the eclipse shape itself.

We applied eclipse mapping techniques (Horne 1985; Rutten, van Paradijs & Tinbergen 1992b; Baptista & Steiner 1993) to our light curves to solve for a map of the disk brightness distribution and for the flux of an additional uneclipsed component in each case. The uneclipsed component accounts for all light that is not contained in the eclipse map in the orbital plane (i.e., light from the secondary star and/or a vertically extended disk wind). The reader is referred to Rutten et al. (1992b) and Baptista, Steiner & Horne (1996) for a detailed description of and tests with the uneclipsed component. For the reconstructions we adopted the default of limited azimuthal smearing of Rutten et al. (1992b), which is better suited for recovering asymmetric structures than the original default of full azimuthal smearing (see Baptista et al. 1996). Simulations that show the ability of the eclipse mapping method to reconstruct asymmetric structures such as spiral arms are presented and discussed by Baptista et al. (2000a). As our eclipse map we adopted a grid of $`75\times 75`$ pixels centered on the primary star with side 2 RL1, where RL1 is the distance from the disk center to the inner Lagrangian point. The eclipse geometry is defined by the mass ratio $`q`$ and the inclination $`i`$. We adopted the parameters of BCC, $`q=0.72`$ and $`i=85\mathrm{°}`$, which correspond to an eclipse width of the disk center of $`\mathrm{\Delta }\varphi =0.1085`$. These combination of parameters ensure that the white dwarf is at the center of the map.

The statistical uncertainties of the eclipse maps were estimated with a Monte Carlo procedure (e.g., Rutten et al. 1992). A set of 20 artificial light curves is generated, in which the data points are independently and randomly varied according to a Gaussian distribution with standard deviation equal to the uncertainty at that point. The light curves are fitted with the eclipse mapping algorithm to produce a set of randomized eclipse maps. These are combined to produce an average map and a map of the residuals with respect to the average, which yields the statistical uncertainty at each pixel.

## 3 Results

The data and the model light curves are shown in the left-hand panels of Fig. A spiral structure in the disk of EX Draconis on the rise to outburst maximum. Horizontal dashed lines indicate the uneclipsed component in each case. The middle panels show the corresponding eclipse maps in the same logarithmic grayscale. The maps in the right-hand panels show the asymmetric part of the maps in the middle panels and are obtained by calculating and subtracting azimuthally-averaged intensities at each radius.

The uneclipsed fluxes are 1.45 and 1.6 mJy, respectively, for the light curves in quiescence and during the rise. The uneclipsed component in quiescence is mostly due to the contribution of the secondary star to the $`R`$ band flux, and is comparable to the quiescent mid-eclipse flux. The larger uneclipsed flux during the rise is possibly the result of additional emission from a (variable) vertically-extended disk wind component, which becomes non-negligible during the outbursts (see Baptista & Catalán 2000).

The light curve in quiescence shows a flat-bottomed asymmetric eclipse that is the result of the almost total occultation of a bright source at disk center and of a faint compact bright spot at disk rim. The spot in the eclipse map lies along the gas stream trajectory and yields a quiescent disk radius of $`R_q=0.51\pm 0.02R_{L1}`$, in good agreement with the value measured by BCC. The width of the eclipse in the rise light curve is the same as in quiescence but the eclipse is skewed towards later phases, indicating the occultation of an asymmetric brightness distribution in the trailing side of the disk. The resulting eclipse map has a symmetrical, centered brightness distribution, and an asymmetric structure in the form of a spiral arm that is clearly visible in the right-hand panel of Fig. A spiral structure in the disk of EX Draconis on the rise to outburst maximum. The eclipse maps obtained with the Monte Carlo procedure consistently show a similar spiral structure. The asymmetric arc contributes about 22 per cent of the total flux of the eclipse map.

In order to trace the distribution in radius and azimuth of the asymmetric structures, we divided the eclipse maps in azimuthal slices (i.e., ‘slices of pizza’) and computed the distance at which the intensity is a maximum for each azimuth. We restricted the search to radii $`R>0.2R_{L1}`$ in order to avoid being affected by the central brightness source. The results are plotted in Fig. A spiral structure in the disk of EX Draconis on the rise to outburst maximum as a function of binary phase (binary phases increase clockwise in the eclipse maps of Fig. A spiral structure in the disk of EX Draconis on the rise to outburst maximum and phase zero coincides with the inner Lagrangian point L1). The lower panel shows the dependency of the maximum intensity with binary phase, while the upper panel gives the radial position of the maximum intensity as a function of binary phase (only the regions for which the radius of maximum intensity is larger than $`0.2R_{L1}`$ are shown). The error bars were derived via Monte Carlo simulations with the eclipse light curves (section 2).

There are clear differences between the asymmetric structures in quiescence and on the rise to maximum. The quiescent bright spot is a well defined, compact source with narrow azimuthal ($`\mathrm{\Delta }\varphi 30\mathrm{°}`$) and radial ($`\mathrm{\Delta }R=0.04R_{L1}`$) extents. On the other hand, the asymmetric structure on the rise map is extended both in azimuth ($`\mathrm{\Delta }\varphi 180\mathrm{°}`$) and radius (from $`R0.2`$ to $`0.43R_{\mathrm{L1}}`$), showing a spiral pattern in which the radius of the maximum intensity continuously increases with binary phase. An estimate of the disk radius at the rising stage can be obtained from the outer radius of the spiral arm, $`R_r=0.43R_{L1}`$.

We computed equivalent Keplerian velocities along the spiral structure by assuming $`M_1=0.75\pm 0.15M_{}`$ and $`R_{\mathrm{L1}}=0.85R_{}`$ (BCC). We obtain velocities in the range $`600800\mathrm{km}\mathrm{s}^1`$, larger than the velocities ($`400600\mathrm{km}\mathrm{s}^1`$) inferred for the spirals seen close to outburst maximum (Joergens et al. 2000). Since the eclipse map and the Doppler tomogram correspond to different outburst stages, we are not able to distinguish whether the different velocities are the consequence of different radii for the spirals, or if they are the signature of sub-Keplerian velocities in spiral shocks – as seen in IP Peg (Baptista et al. 2000a).

## 4 Discussion

Figures A spiral structure in the disk of EX Draconis on the rise to outburst maximum and A spiral structure in the disk of EX Draconis on the rise to outburst maximum show that the spiral arm is brighter and closer to disk center than the quiescent bright spot. The ratio of the intensity of the outer parts of the spiral arm and that of the quiescent bright spot ($`3.5`$, lower panel of Fig. A spiral structure in the disk of EX Draconis on the rise to outburst maximum) is much larger than the expected increase in accretion luminosity owing to the reduction in radius (for a constant mass accretion rate), $`L_{acr}(R_r)/L_{acr}(R_q)=R_q/R_r=1.2`$. This discrepancy indicates that the mass accretion rate on the rise to maximum is larger than the mass inflow at the quiescent bright spot by a factor of $`3`$111 In fact, this is a lower limit, since the bulk of the accretion luminosity at the smaller $`R_r`$ radius should move to shorter wavelengths due to the higher effective temperature, making the ratio of intensities in the R-band $`I_\nu (R_r)/I_\nu (R_q)<L_{acr}(R_r)/L_{acr}(R_q)`$.. Moreover, if the outer radius of the spiral arm can be assigned to the disk radius, the comparison of $`R_q`$ and $`R_r`$ indicates that the accretion disk shrinks at the onset of the outburst. This result is in agreement with the predictions of the MTI model (Warner 1995), and suggests that the outburst in EX Dra is driven by a burst of enhanced mass-transfer from the secondary star.

We now turn our attention to the nature of the observed asymmetric structure. It is hard to understand the observed one-armed spiral structure in terms of tidally-induced spiral shocks. First, because the accretion disk radius at the corresponding outburst stage (again assumed to be of the same order of the outer radius of the spiral arm, $`R_r=0.43R_{L1}`$) is significantly smaller than the range in radius ($`0.560.75R_{L1}`$) required to excite spiral shocks strong enough to be observed (Steeghs & Steele 1999). Secondly, because the tidal effect of the secondary star should, in principle, lead to the formation of a two-armed spiral structure (e.g., Yukawa, Boffin & Matsuda 1997; Armitage & Murray 1998), instead of the observed one-armed spiral.

An alternative explanation for the spiral structure in terms of enhanced gas stream emission might be advocated, although the position of the spiral arm is clearly not consistent with the location of the gas stream trajectory (see Fig.A spiral structure in the disk of EX Draconis on the rise to outburst maximum). Inconsistently low mass ratios ($`q\stackrel{<}{}0.1`$) would be required to bring the gas stream trajectory in closer agreement with the observed spiral pattern. If valid, this alternative would also indicate an MTI origin for the outbursts of EX Dra.

A third explanation may be devised if we note that outbursting accretion disks might depart from the geometrically thin disk approximation. At the high inclination of EX Dra, this may have important effects on the eclipse maps. It may be possible that the observed spiral pattern is the combination of emission from a bright spot/stream running along the edge of the disk and the (apparently) enhanced emission from the far (with respect to the secondary star) side of a flared accretion disk. This scenario makes the bright spot on the rise a factor of 3 brighter than the quiescent bright spot, also in agreement with the predictions of the MTI model.

The observed spiral structure bears resemblance with the transient arc-shaped asymmetric structures seen in 2-D and 3-D numerical simulations of accretion flows, at the onset of mass transfer (Makita, Miyawaki & Matsuda 1998; Makita & Matsuda 1998) or at the start of an instability-driven outburst (Różyczka & Spruit 1993). In this regard, useful constraints on future numerical simulations attempting to reproduce the observed spiral pattern may be obtained by noting that the individual light curves combined to produce the average rise light curve cover a time interval of 5 binary cycles. Their eclipse shape is the same under the uncertainties. Hence, the spiral structure seems to be stationary on a reference frame co-rotating with the binary and stable for a time length of, at least, 5 binary orbits.

Finally, we note that, contrary to the usual assumption, inside-out outbursts are not necessarily inconsistent with the MTI model, provided that the burst of mass transfer leads to the formation of a spiral structure well inside the primary lobe, that expands thereafter both inwards and outwards under the effects of mass and angular momentum redistribution to form a fully developed accretion disk.

We thank Emilios Harlaftis for helpful discussions and comments which stimulated the writing of this paper. This work was partially supported by the PRONEX/Brazil program through the research grant FAURGS/FINEP 7697.1003.00. RB acknowledges financial support from CNPq/Brazil through grant no. 300 354/96-7. MSC acknowledges financial support from a PPARC post-doctoral grant during part of this work.


text:
# 1 Introduction

## 1 Introduction

$`O(N)`$ models are of general relevance to condensed matter physics and to quantum field theory, because many physical systems exhibit a second-order phase transition with universal properties, which may be derived analytically from $`O(N)`$ symmetric effective field theories . A further aspect of these theories is remarkable. Due to the existence of massless Goldstone modes in $`O(N)`$ models with $`N>1`$ and dimension $`2<d4`$  singularities are expected on the whole coexistence line $`T<T_c,H=0`$, in addition to the known critical behaviour at $`T_c`$. For $`N=2`$ and $`d=3`$ there exists a rigorous proof . Numerically these theoretical predictions have been recently confirmed by simulations of the $`3d`$ $`O(4)`$ model . In the same article a parametrization of the equation of state including the Goldstone effect has been worked out and compared to perturbative predictions . Here we want to extend these numerical studies to the case of the $`O(2)`$ model. For the latter approximate representations of the equation of state have been derived, starting from high-temperature expansions .

Our interest in the $`O(2)`$ model is moreover motivated by its relation to the staggered formulation of quantum chromodynamics (QCD) on the lattice. At finite temperature QCD undergoes a chiral phase transition. For two degenerate light-quark flavors this transition is supposed to be of second order in the continuum limit and to belong to the same universality class as the $`3d`$ $`O(4)`$ model -. QCD lattice data have therefore been compared to the $`O(4)`$ scaling function, as determined numerically in . For staggered fermions the comparison to $`O(4)`$ is still not conclusive -, but results for Wilson fermions  seem to agree quite well with the predictions, though for the Wilson action the chiral symmetry is only restored in the continuum limit. In the staggered formulation a part of the chiral symmetry is remaining even for finite lattice spacing, and that is $`O(2)`$. For the test of lattice QCD with staggered fermions it is therefore also important to know the corresponding universal $`O(2)`$ scaling function or equation of state.

The $`O(2)`$-invariant nonlinear $`\sigma `$-model (or $`XY`$ model), which we want to investigate in the following is defined by

$$\beta =J\underset{<i,j>}{}𝐒_i𝐒_j𝐇\underset{i}{}𝐒_i,$$
(1)
where $`i`$ and $`j`$ are nearest-neighbour sites on a $`d`$dimensional hypercubic lattice, and $`𝐒_i`$ is an 2-component unit vector at site $`i`$. It is convenient to decompose the spin vector $`𝐒_i`$ into a longitudinal (parallel to the magnetic field $`𝐇`$) and a transverse component

$$𝐒_i=S_i^{}\widehat{𝐇}+𝐒_i^{}.$$
(2)
The order parameter of the system, the magnetization $`M`$, is then the expectation value of the lattice average $`S^{}`$ of the longitudinal spin component

$$M=<\frac{1}{V}\underset{i}{}S_i^{}>=<S^{}>.$$
(3)
On finite lattices and $`H=0`$ system flips occur and lead to $`<S^{}>=0`$. Therefore one usually resorts to approximate order parameter definitions, as e.g. $`<|S^{}|>`$ \- see the discussion in Tapalov and Blöte . Here this is unnecessary since we always work at finite $`H`$.

There are two types of susceptibilities. The longitudinal susceptibility is defined as usual by the derivative of the magnetization, whereas the transverse susceptibility corresponds to the fluctuation of the lattice average $`𝐒^{}`$ of the transverse spin component

$`\chi _L`$ $`=`$ $`{\displaystyle \frac{M}{H}}=V(<S_{}^{}{}_{}{}^{2}>M^2),`$ (4)
$`\chi _T`$ $`=`$ $`V<𝐒_{}^{}{}_{}{}^{2}>.`$ (5)
Both susceptibilities are predicted to diverge on the coexistence line $`T<T_c,H=0`$. The Goldstone singularities lead for all temperatures below $`T_c`$ for small $`H`$ to strong finite-size effects, which have been studied in the context of chiral perturbation theory for the three-dimensional $`O(2)`$ model already in . An explicit check on the $`H`$-dependence of the magnetization for $`V\mathrm{}`$ close to $`H=0`$ is however lacking. In comparison to the $`O(4)`$ model, we expect the Goldstone effects to be weaker, because we have only one transverse spin component in $`O(2)`$, not three like in $`O(4)`$. Another difference to $`O(4)`$ is that sizeable corrections to scaling appear in the $`XY`$ model . The determination of the universal equation of state of the three-dimensional $`O(2)`$ model, therefore requires a careful consideration of the corrections. We could have avoided the problems with these corrections by working, as it was done in Ref.  at $`\lambda =2.1`$, instead of $`\lambda =\mathrm{}`$. For a better comparison to the $`O(4)`$ case , we chose however to use the same model definition. Also there is then no need for an additional modulus update.

The plan of the paper is as follows. In the next section we briefly discuss the perturbative predictions for the magnetization and the susceptibilities at low temperatures, the form and the analytic results for the magnetic equation of state. A more general review can be found in . Our numerical results are presented in Section 3, the equation of state is determined and parametrized in Section 4. We close with a summary and our conclusions in Section 5.

## 2 Perturbative Predictions and Equation of State

For $`T<T_c`$ the system is in a broken phase, i.e. the magnetization $`M(T,H)`$ attains a finite value $`M(T,0)`$ at $`H=0`$. We explicitly assume here $`H>0`$. As a consequence the transverse susceptibility, which is directly related to the fluctuation of the Goldstone modes, diverges as $`H^1`$ when $`H0`$ for all $`T<T_c`$. This is immediately clear from the identity

$$\chi _T=\frac{M(T,H)}{H},$$
(6)
which is valid for all values of $`T`$ and $`H`$. It is non-trivial that also the longitudinal susceptibility is diverging on the coexistence line for $`2<d4`$. The leading term in the perturbative expansion for $`2<d<4`$ is $`H^{d/22}`$ . The predicted divergence in $`d=3`$ is thus

$$\chi _L(T<T_c,H)H^{1/2}.$$
(7)
This is equivalent to an $`H^{1/2}`$-behaviour of the magnetization near the coexistence curve

$$M(T<T_c,H)=M(T,0)+cH^{1/2}.$$
(8)
An interesting question is whether this form of the magnetization is compatible with the general Widom-Griffiths form of the equation of state  describing the critical behaviour in the vicinity of $`T_c`$. It is given by

$$y=f(x),$$
(9)
where

$$yh/M^\delta ,xt/M^{1/\beta }.$$
(10)
The variables $`t`$ and $`h`$ are the normalized reduced temperature $`t=(TT_c)/T_0`$ and magnetic field $`h=H/H_0`$. We take the usual normalization conditions

$$f(0)=1,f(1)=0.$$
(11)
The critical exponents $`\delta `$ and $`\beta `$ appearing in Eqs. 9 and 10 specify all the other critical exponents

$$d\nu =\beta (1+\delta ),\gamma =\beta (\delta 1),\nu _c=\nu /\beta \delta .$$
(12)
Possible dependencies on irrelevant scaling fields and exponents are however not taken into account in Eq. 9, the function $`f(x)`$ is universal. Another way to express the dependence of the magnetization on $`t`$ and $`h`$ is

$$M=h^{1/\delta }f_G(t/h^{1/\beta \delta }),$$
(13)
where $`f_G`$ is a scaling function. This type of scaling equation is used for comparison to QCD lattice data. The scaling forms in Eqs. (9) and (13) are clearly equivalent.

The equation of state (9) has been derived by Brézin et al.  to order $`ϵ^2`$ in the $`ϵ`$-expansion, where $`ϵ=4d`$. The resulting approximation has been considered by Wallace and Zia  in the limit $`x1`$, i.e. at $`T<T_c`$ and close to the coexistence curve. In this limit the function is inverted to give $`x+1`$ as a double expansion in powers of $`y`$ and $`y^{d/21}`$

$$x+1=\stackrel{~}{c}_1y+\stackrel{~}{c}_2y^{d/21}+\stackrel{~}{d}_1y^2+\stackrel{~}{d}_2y^{d/2}+\stackrel{~}{d}_3y^{d2}+\mathrm{}.$$
(14)
The coefficients $`\stackrel{~}{c}_1`$, $`\stackrel{~}{c}_2`$ and $`\stackrel{~}{d}_3`$ are then obtained from the general expression of . The above form is motivated by the $`H`$-dependence in the $`ϵ`$-expansion of $`\chi _L`$ at low temperatures . In Section 4 we propose a fit of our Monte Carlo data to the perturbative form of the equation of state, using (14).

As for the large-$`x`$ limit (corresponding to $`T>T_c`$ and small $`H`$), the expected behaviour is given by Griffiths’s analyticity condition

$$f(x)=\underset{n=1}{\overset{\mathrm{}}{}}a_nx^{\gamma 2(n1)\beta }.$$
(15)
## 3 Numerical Results

Our simulations are done on three-dimensional lattices with periodic boundary conditions and linear extensions $`L=8160`$ using the cluster algorithm of Ref. . The value $`J_c=0.454165(4)`$, obtained in simulations of the zero-field model , is used in the following. We have simulated at increasingly larger values of $`L`$ at fixed values of $`J=1/T`$ (i.e. at fixed temperature $`T`$) and $`H`$ in order to eliminate finite-size effects.

As an example we show in Fig. 1 part of our data for the magnetization for low temperatures, at $`J=0.5,0.47`$, and at $`J_c`$, plotted versus $`H^{1/2}`$. The picture is rather similar to the one obtained in $`O(4)`$ : strong finite-size effects appear for small $`H`$ and persist as one moves away from $`T_c`$, the results from the largest lattices are at first sight linear in $`H^{1/2}`$, as predicted by Eq. 8, below but not at $`T_c`$. A closer look at the curves reveals however subtle differences compared to the $`O(4)`$ case, which are confirmed by the attempt to determine the value of $`M(T,H=0)`$ by a linear fit in $`H^{1/2}`$. Very close to $`H=0`$ the fixed temperature curves become slighty flatter, leading to a higher value for $`M(T,0)`$ than expected from the data at larger $`H`$ values. This behaviour is more pronounced close to $`T_c`$ than at lower temperatures. The determination of the magnetization in the thermodynamic limit ($`V\mathrm{}`$) on the coexistence line requires therefore more precise data and smaller values of $`H`$ than for the $`O(4)`$ model. In order to extrapolate the data to $`H0`$ and $`V\mathrm{}`$ we apply two different strategies. The first is to extend the linear form in $`H^{1/2}`$, Eq. 8, to a quadratic one

$$M(T<T_c,H)=M(T,0)+c_1H^{1/2}+c_2H,$$
(16)
and to fit the data from the largest lattices, which we assume to represent data on an infinite volume lattice, to this form. The second way to find $`M(T,0)`$ is just opposite to the first. Here we exploit the $`L`$ or volume dependence at fixed $`J`$ and fixed small $`H`$ to determine via chiral perturbation theory (CPT) the magnetization $`\mathrm{\Sigma }`$ of the continuum theory for $`V\mathrm{},H=0`$, which is related to

$`M(T,0)`$ by

$$M(T,0)=\frac{\mathrm{\Sigma }}{\sqrt{J}}.$$
(17)
According to the $`ϵ`$-expansion of CPT the quantity $`\mathrm{\Sigma }`$ and the Goldstone boson decay constant $`F`$ are the only parameters which determine the finite-size corrections to the partition function to order $`1/L^2`$ . We summarize the relevant formulae for the three-dimensional $`XY`$ model from Ref. . In our notation the magnetization is

$$M=\frac{\mathrm{\Sigma }}{\sqrt{J}}u[\rho _1\eta +2\rho _2\alpha ^2],$$
(18)
where $`u=\rho _1\mathrm{\Sigma }HV/\sqrt{J}`$, $`\alpha =1/(F^2L)`$ and

$$\rho _1=1+\frac{1}{2}\beta _1\alpha +\frac{1}{8}(\beta _1^22\beta _2)\alpha ^2;\rho _2=\frac{1}{4}\beta _2.$$
(19)
For a symmetric three-dimensional box $`\beta _1=0.225785`$ and $`\beta _2=0.010608`$. The parameter $`\eta `$ is given in terms of modified Bessel functions as

$$\eta =\frac{I_1(u)}{uI_2(u)}.$$
(20)
By construction the $`ϵ`$-expansion is only applicable in a range where $`m_\pi L<1`$ which translates into the condition

$$H\frac{\mathrm{\Sigma }}{\sqrt{J}}<\left(\frac{F}{L}\right)^2.$$
(21)
In Fig. 2 we show the data for the magnetization as a function of $`H^{1/2}`$ for six fixed values of $`J`$ in the low-temperature phase. We observe a remarkable coincidence of the fits according to Eq. 16 with the CPT results at $`H=0`$. Details of the fits are presented in Table 1. The $`\chi ^2`$ per degree of freedom for the first fit type is in the range 0.01-0.5 . The corresponding numbers for the CPT fits are of the order of 1. Like for the $`O(4)`$ model we see in Fig. 2 that the region where $`M`$ is linear in $`H^{1/2}`$ shrinks when $`T`$ approaches $`T_c`$, the value of $`c_2`$ in Table 1 increases correspondingly. At $`J=0.460`$ the fit parabola coincides only on a small piece with the largest volume data (here we simulated even on $`L=160`$ lattices), the errors on $`c_1`$ and $`c_2`$ become large, though the value at $`H=0`$ agrees again nicely with the CPT result. We have further convinced ourselves that $`M(T,0)`$ can also be obtained from calculations at $`H=0`$ with the modulus definition. The approach to the thermodynamic limit is here from above.

In the neighbourhood of the critical temperature the results for $`M(T,0)`$ should show the usual critical behaviour. Contrary to the $`O(4)`$ case we expect here a sizeable correction to the leading scaling behaviour . In order to determine the critical amplitude $`B`$ we therefore make the following ansatz

$$M(T\mathrm{}<T_c,H=0)=B(T_cT)^\beta [1+b_1(T_cT)^{\omega \nu }+b_2(T_cT)].$$
(22)
As an input we take the critical exponents from Ref.

$$\beta =\mathrm{\hspace{0.33em}0.3490}(6),\nu =\mathrm{\hspace{0.33em}0.6723}(11),\omega =\mathrm{\hspace{0.33em}0.79}(2).$$
(23)
A fit of all points of Table 1, apart from the one at $`J=0.46`$, using the form (22) leads to the result $`B=0.945(5),b_1=0.053(23)`$ and $`b_2=0.098(23)`$. In Fig. 3a we show this fit and also the leading term separately.

As the critical point is reached the $`H`$-dependence of the magnetization changes to satisfy critical scaling. We therefore fit the data from the largest lattice sizes at $`T_c`$ to the form

$$M(T_c,H)=d_cH^{1/\delta }[1+d_c^1H^{\omega \nu _c}].$$
(24)
A further term proportional to $`H`$ is unnecessary, because the corrections to scaling are much smaller here, than on the coexistence line. The largest $`L`$ data can be fitted very well with the ansatz (24) and the critical exponents as input, as can be seen in Fig. 3b. We find $`d_c=0.978(2)`$ and $`d_c^1=0.075(5)`$.

We have performed in addition some simulations in the high temperature phase of the model. Here the finite size effects are rather small. In Fig. 4 we show the results. One observes that with increasing $`T`$ (decreasing $`J`$) the region of linear dependence on $`H`$ of the magnetization becomes larger.

## 4 Determination of the Equation of State

In the last section we derived from our Monte Carlo data the magnetization in the thermodynamic limit. At temperatures below $`T_c`$ this was achieved for the larger $`H`$ values by simply performing simulations on lattices with increasing size until all volume dependence was gone. For the small $`H`$ region down to $`H=0`$ we took then advantage of the Goldstone effect to extrapolate our data and confirmed this by the CPT results for $`\mathrm{\Sigma }`$. In the high-temperature phase we reached the thermodynamic limit already on lattices with $`L72`$ for all $`H`$ values.

The two critical amplitudes $`B`$ and $`d_c`$ can now be utilized to normalize temperature and magnetic field according to Eq. (11). We find from our fits (22) and (24)

$$T_0=B^{1/\beta }=1.18(2),\mathrm{and}H_0=d_c^\delta =1.11(1).$$
(25)
As mentioned already, the equation of state does not account for possible corrections to scaling. A more general form of Eq. (13) is

$$Mh^{1/\delta }=\mathrm{\Psi }(th^{1/\beta \delta },h^{\omega \nu _c}).$$
(26)
Expanding the function $`\mathrm{\Psi }`$ in $`h^{\omega \nu _c}`$ leads to

$$Mh^{1/\delta }=f_G(th^{1/\beta \delta })+h^{\omega \nu _c}f_G^{(1)}(th^{1/\beta \delta })+\mathrm{}.$$
(27)
In order to obtain the scaling function $`f_G`$ we therefore perform quadratic fits to our data in $`h^{\omega \nu _c}`$ at fixed values of $`th^{1/\beta \delta }`$ in the low-temperature region, where the corrections are strong. Fortunately, the corrections are very small in the high-temperature region and the data scale directly. Altogether data with $`H0.0075`$ and $`0.43J0.55`$ were used for this purpose. In Fig. 5 we have plotted the uncorrected and the final results for the scaling function $`f_G`$.

Like in the case of $`O(4)`$  we want to describe the low-temperature part of the equation of state by its perturbative form as discussed in Section 2, but with nonperturbative coefficients. We then want, again as in Ref. , to interpolate this result with a fit to the large-$`x`$ form (15). The variables $`x`$ and $`y`$ are simply related to the scaling function $`f_G`$ and its argument by

$$y=f_G^\delta ,x=(t/h^{1/\beta \delta })f_G^{1/\beta }.$$
(28)
We first perform a fit to $`x`$ in the interval $`[1,\mathrm{\hspace{0.17em}1.4}]`$ (which corresponds to $`y[0,3]`$) with the three leading terms in (14)

$$x_1(y)+1=(\stackrel{~}{c}_1+\stackrel{~}{d}_3)y+\stackrel{~}{c}_2y^{1/2}+\stackrel{~}{d}_2y^{3/2}.$$
(29)
In the fit we require $`\stackrel{~}{d}_2=1(\stackrel{~}{c}_1+\stackrel{~}{d}_3+\stackrel{~}{c}_2)`$ to fix the normalization $`y(0)=1`$. We obtain

$$\stackrel{~}{c}_1+\stackrel{~}{d}_3=\mathrm{\hspace{0.17em}0.352}(30),\stackrel{~}{c}_2=\mathrm{\hspace{0.17em}0.592}(10).$$
(30)
The fit describes the corrected scaling function at $`T<T_c`$ and and the direct data in the high temperature phase up to $`x1.7`$. This confirms, as in $`O(4)`$, that the expression (14) is valid also away from $`x1`$. Our coefficients can be compared to those calculated perturbatively for $`N=2`$ in Ref.

$$\stackrel{~}{c}_1+\stackrel{~}{d}_3=\mathrm{\hspace{0.17em}0.818},\stackrel{~}{c}_2=\mathrm{\hspace{0.17em}0.229}.$$
(31)
The result (30) is much closer to the coefficients found in $`O(4)`$  than to the $`ϵ`$-expansion results (31), though the Goldstone effect is decreasing somewhat with decreasing $`N`$: $`\stackrel{~}{c}_2`$ and the $`x`$ region where the ansatz (29) is valid are a little smaller.

For large $`x`$ we have done a 2-parameter fit of the behaviour (15), in the corresponding form for $`x`$ in terms of $`y`$

$$x_2(y)=ay^{1/\gamma }+by^{(12\beta )/\gamma }.$$
(32)
Considering data points with $`y`$ in the interval  we find

$$a=\mathrm{\hspace{0.17em}1.2595}(30),b=1.163(20).$$
(33)
Expression (32) describes the data for $`x>1.5`$. The small- and large-$`x`$ curves cover the whole range of values of $`x`$ remarkably well. The two curves overlap around $`y3.5`$. This is shown in Fig. 6a. We therefore interpolate them smoothly with the ansatz

$$x(y)=x_1(y)\frac{y_0^6}{y_0^6+y^6}+x_2(y)\frac{y^6}{y_0^6+y^6},$$
(34)
where $`y_0=3.5`$. In Fig. 6b we compare our interpolation (34) of the equation of state to the high-temperature data and to two parametric forms for $`f(x)`$ obtained from the high-temperature expansion in Ref. . According to the authors , the difference between these two curves gives an idea of the uncertainty of their approach. For large $`x`$ this difference can be traced back to the uncertainty in the universal amplitude ratio $`R_\chi `$ (Eq. (49) of Ref. ), which was determined in  to $`R_\chi =1.4(1)`$, that is with a 7% error. In fact, $`R_\chi `$ is related to the quantity $`a`$ of our large-$`x`$ parametrization by

$$R_\chi =a^\gamma =\mathrm{\hspace{0.33em}1.356}(4),$$
(35)
which is in agreement with the result from Ref. . Their low temperature parameter $`c_f`$ may be calculated from our value of $`\stackrel{~}{c}_2`$ and amounts to $`c_f=2.85(7)`$. Finally we show in Fig. 7 the scaling function $`f_G`$ obtained parametrically from $`x(y)`$ in (34) and the corresponding scaling function for $`O(4)`$. As can be seen from Fig. 7, the dependence on $`t/h^{1/\delta \beta }`$ is similar, but the $`O(2)`$ curve is flatter than the one for $`O(4)`$.

## 5 Summary and Conclusions

We have simulated the three-dimensional $`O(2)`$ model on cubic lattices as a function of the magnetic field $`H`$ and the temperature $`T`$. From the behaviour of the measured magnetization $`M`$ below the critical temperature $`T_c`$ and close to the coexistence line ($`H=0`$), we could clearly verify the Goldstone-mode effects. This was done in two independent ways which led to the same values $`M(T<T_c,H0)`$ for the magnetization in the thermodynamic limit $`V\mathrm{}`$ : on one hand we were able to profit by the observed Goldstone behaviour to extrapolate our data to $`H0`$ and on the other hand this result was confirmed from the finite-size dependence induced by the Goldstone modes, using chiral perturbation theory.

The $`M`$ values on the coexistence line were subsequently used to calculate the critical amplitude $`B`$ of the magnetization. Not unexpectedly, we found here strong corrections to scaling. On the critical line $`T=T_c`$ we computed then the second critical amplitude $`d_c`$ of the magnetization. Here, the corrections to scaling are less pronounced ; in the high-temperature phase they are unimportant.

The problem with corrrections to scaling appeared again in the determination of the equation of state in the low-temperature phase. By generalizing the scaling equation to include possible corrections, we were able to derive the universal form of the equation of state for the $`O(2)`$ model from our data. This form was then analyzed in a similar manner as it was done for the $`O(4)`$ model in Ref. . In particular, we found again an efficient parametrization of the equation of state in the low $`x`$ region with the perturbative form (29), which is based on expansion (14) by Wallace and Zia . The effect due to the Goldstone modes is nearly as large as in the $`O(4)`$ case, as can be seen from the coefficient $`\stackrel{~}{c}_2`$ and much larger than the prediction from the $`ϵ`$-expansion.

Like for $`O(4)`$ we found a very good large-$`x`$ fit to the high-temperature data. The coefficient $`a`$ obtained from the fit implies a value of $`R_\chi =1.356(4)`$, in agreement with parametrizations of the equation of state by Campostrini et al. . A further, indirect check of the equation of state is the computation of the universal ratio $`A_+/A_{}`$. This requires however an integration of the magnetization with respect to $`h`$ and derivatives with respect to $`t`$. We shall consider this in the near future.

Upon interpolation with the low-$`x`$ curve, a complete description for the equation of state is obtained, which can be plotted parametrically also for the scaling function. In Ref.  the lattice QCD data of the MILC Collaboration for $`N_\tau =4`$ were compared to the $`O(4)`$ scaling function. The test failed because the data were indicating a steeper scaling function. Since the $`O(2)`$ scaling function is flatter than the one for $`O(4)`$, the situation will be worse there. A way out may be the comparison to finite-size-scaling functions, since lattice QCD is presumably far from the thermodynamic limit. We are investigating this currently.

Acknowledgements

We thank Attilio Cucchieri for helpful suggestions and comments. This work was supported by the Deutsche Forschungsgemeinschaft under Grant No. Ka 1198/4-1.


text:
# Role of Bose enhancement in photoassociation11footnote 1Paper submitted to the special issue of J. Mod. Opt. on “Seminar on Fundamentals of Quantum Optics V”, Kühtai, Austria, January 16-21, 2000; ed. F. Ehlotzky.

## 1 Introduction

The theoretical realization is now emerging that it may be possible to effect analogs of coherent optical transients, such as Rabi flopping , rapid adiabatic passage , and stimulated Raman adiabatic passage (STIRAP) , in photoassociation of a Bose-Einstein condensate (BEC). In contrast, the feasibility of transients in PA of a thermal nondegenerate gas has been controversial. There are predictions of STIRAP in ordinary gases , whereas our results suggest otherwise .

The purpose of the present paper is to clarify the status of coherent optical transients in both nondegenerate and degenerate thermal gases. Our key concept is Bose enhancement. Suppose we have a large number of atoms $`N`$ in one quantum state, as in a BEC . Due to Bose-Einstein statistics, the transition matrix element referring to the BEC will pick up a multiplicative factor $`\sqrt{N}`$ as compared to the matrix element for a single atom. It turns out that Bose enhancement will make a difference. In the thermodynamic limit, when the volume and particle number of a gas tend to infinity while the density remains constant, coherent optical transients will vanish in photoassociation of a nondegenerate gas, but Bose enhancement will facilitate such transients in a condensate. In the rest of the present paper we outline the details of this argument using STIRAP as our explicit example.

## 2 Ordinary STIRAP

As a prelude to our development we briefly reiterate the salient features of STIRAP in an ordinary three-level $`\mathrm{\Lambda }`$ system , as in Fig. 1. Two laser fields are tuned to exact two-photon resonance between the two lower states $`a`$ and $`g`$ of the $`\mathrm{\Lambda }`$ system. First the laser intensities are arranged in such a way that the coupling between states $`g`$ and $`b`$ is much stronger than the coupling between $`a`$ and $`b`$. One of the dressed states of the three-level system then coincides with the bare state $`a`$. When the laser intensities are adjusted in such a way that the coupling becomes much stronger in the transition from $`a`$ to $`b`$, the dressed state that initially coincided with the bare state $`a`$ turns instead into the bare state $`g`$.

Suppose now that, before the adjustment of laser intensities, the system started in the bare state $`a`$. It then also started in the dressed state whose fate we follow. If the lasers were adjusted slowly enough, adiabaticity guarantees that the system stays in the dressed state all along. In the end, the system therefore winds up in the bare state $`g`$. Moreover, a detailed analysis shows that, in the ideal adiabatic limit, there is never any population at all in the intermediate state $`b`$.

In short, the intensities of the two light pulses are switched in such a way that initially the coupling between the unoccupied states is much stronger than the coupling to the occupied state, and the same situation prevails when the coupling strengths have reversed. This is called counter-intuitive pulse order. The result is that the system is transferred from the initial state $`a`$ to the final state $`g`$. Besides, ideally, the system never visits the intermediate state $`b`$, which is a major virtue if $`b`$ is plagued by dissipation.

## 3 Quasicontinuum photoassociation

In a process of photoassociation two atoms and a photon combine to make the corresponding diatomic molecule. As the two atoms may be considered to be a dissociated state of the molecule, photoassociation is really about free-bound transitions. There are internal atomic and molecular states involved in photoassociation as well, but we assume that, by properly choosing the laser frequency, one has selected a transition from the initial state of the atoms to a unique rovibrational state of the molecule.

Quantum optics normally deals with bound-bound transitions responsible for characteristic spectral lines, and with bound-free transition that manifest themselves in decay processes such as spontaneous emission. From this angle, photoassociation is irreversible decay in reverse, and may seem to violate the second law of thermodynamics. Nonetheless, photoassociation spectroscopy is alive and well, and is the source of the most accurate molecular structure data available at this time .

To model free-bound transitions, especially in a thermal sample, it is desirable to start with two atoms whose relative motion is in an eigenstate of energy. Unfortunately, as energy eigenstates of the relative motion for unconfined atoms are not normalizable, there is no mathematically sound way to write down such a quantum state. The cure is straight from textbooks of quantum mechanics: Assume that the relative motion is restricted to a quantization volume $`V`$, and at the end of the calculations take the limit $`V\mathrm{}`$.

This is the stratagem of our quasicontinuum (QC) approach . The dissociation continuum is broken up into discrete states, a QC, whereupon problems with the normalization of the states vanish. In an unexpected windfall, the mathematics also turns out to work out in such a way that the experience about few-level systems gained over decades of quantum optics and laser spectroscopy is directly transferable to understanding photoassociation.

As it comes to the theme of the present paper, the observation of most immediate relevance is that the dipole matrix element between any free (two-atom) state and the bound (rovibrational molecular) state scales with the quantization volume as $`dV^{1/2}`$. There is no mystery to this. A unit-normalized dissociated state fills the entire volume $`V`$, so its square is proportional to $`1/V`$ and the normalization constant is $`V^{1/2}`$. On the other hand, the normalization constant of the bound state is independent of $`V`$. Eventually the $`V^{1/2}`$ makes its way to the free-bound dipole matrix element.

The immediate consequence is that in the limit of an infinitely large quantization volume, the dipole matrix element between the bound state and any single QC state vanishes. Accordingly, when the volume is increased and Bose enhancement is absent, the photoassociation coupling can always be treated with perturbation theory .

## 4 Nondegenerate thermal gas

### 4.1 Theoretical method

The main reason why we have insisted on eigenstates of energy is the basic random-phase postulate of statistical mechanics, which states that the thermal density operator is diagonal in eigenstates of energy. It is therefore always permissible to do whatever analysis one is aiming at by first assuming that the system starts in a given eigenstate of energy $`m`$, and at the end of the calculations averaging the results over the thermal probability distribution of the states $`m`$.

When discussing a nondegenerate thermal gas, we make another assumption as well. We analyze photoassociation for just two atoms, and calculate the free-bound transition rate, $`R`$. However, a typical experiment involves $`N2`$ atoms. When thinking of photoassociation of any given ‘probe’ atom, we add the transition rates due to all colliders, so that the photoassociation rate per atom becomes $`NR`$.

In the limit of infinite quantization volume, the rate of photoassociation for two atoms vanishes with the coupling matrix element as $`Rd^21/V`$. We of course expect as much since two atoms cease to collide in an infinite volume. Nonetheless, for $`N`$ atoms the photoassociation rate per atom is $`NRN/V=\varrho `$, proportional to the density of the gas. In the thermodynamic limit when both $`N`$ and $`V`$ tend to infinity in such a way that $`\varrho `$ remains constant, the photoassociation rate per atom has a finite limit proportional to density. The result is reasonable, and our QC approach  in fact exactly reproduces the free-bound transition rate obtained from collision theory .

### 4.2 The demise of STIRAP

In our model, and in current practice , the primary photoassociated molecular state is reached by absorption of a photon. This invariably means that the molecule is subject to spontaneous emission, and decays away. Moreover, the reverse of PA, photodissociation, tends to break molecules back into atoms. It is then a natural idea to add a second laser field tuned between the initially photoassociated state and another (more) stable molecular state, and to attempt to utilize STIRAP to avoid losses from the primary photoassociated state . We sketch such a scheme in Fig. 2.

Unfortunately, this idea hits a roadblock. According to the random-phase postulate, for a thermal sample we may assume that the system starts in a given QC state $`m`$, as denoted in Fig. 2. But the free-bound matrix element $`d`$ scales with the quantization volume $`V`$ as $`dV^{1/2}`$, and tends to zero with $`V\mathrm{}`$. On the contrary, the bound-bound matrix element is volume independent. The free-bound coupling is therefore always small compared to the coupling for bound-bound transitions, and it is impossible to effect the counter-intuitive reversal of the coupling strengths needed for STIRAP . Of course, the final average over a distribution of initial states $`m`$ is not expected to create STIRAP either.

We are not arguing that there would be no two-color photoassociation; see Ref. . However, in the absence of any evidence to the contrary, we do not believe that the advantages of counter-intuitive pulse order and STIRAP, protection from decay of the intermediate state and the ensuing improvement in transfer efficiency, will materialize in free-bound-bound photoassociation of a thermal nondegenerate gas.

## 5 Bose-Einstein condensate

### 5.1 Theoretical method

Calculating the free-bound transition rate for one pair of atoms and then multiplying by the number of available colliders is a process which implicitly assumes that we can distinguish between the atoms. Such an approach is fundamentally flawed in the case of a BEC. Instead, we have adopted a phenomenological second-quantized Hamiltonian for photoassociation .

Our basic premises are that one may treat atoms and molecules as bosons in their own right, and that photoassociation conserves momentum. Given the momentum representation for atoms and molecules, annihilation operators $`a_𝐤`$ and $`b_𝐤`$, the part of the Hamiltonian responsible for photoassociation reads,

$$H=\mathrm{}\frac{1}{2}\underset{𝐤,𝐩,𝐪}{}𝐝(𝐤𝐩)𝐄_𝐪b_{𝐤+𝐩+𝐪}^{}a_𝐤a_𝐩+\mathrm{}.$$
(1)
Here $`𝐄_𝐪`$ is the Fourier component $`𝐪/\mathrm{}`$ of the positive frequency part of the electric field driving photoassociation. By translational symmetry, the dipole matrix element may depend only on the difference of the momenta of the atoms, $`𝐤𝐩`$, and in the dipole approximation it cannot depend on the photon momentum $`𝐪`$. The term written down is simply a sum of processes in which two atoms with momenta $`𝐤`$, $`𝐩`$ and a photon with momentum $`𝐪`$ are combined into a molecule with momentum $`𝐤+𝐩+𝐪`$.

It remains to determine the values of the dipole matrix elements $`𝐝`$ in the Hamiltonian. We do this by demanding that for the nondegenerate thermal gas the results from the Hamiltonian (1) be the same as we obtain from our QC approach. In the process a number of subtleties come up having to do, e.g. with Bose-Einstein statistics and the Wigner threshold law for photodissociation . The bottom line, though, is that we know how to deduce the matrix elements from considerations such as the standard molecular-structure calculations, or measurements of the photodissociation rate.

Consider now an ideal zero-momentum condensate photoassociated by a plane wave of light, where all photons have momentum $`𝐪`$. The molecules made by photon absorption all have the momentum $`𝐪`$. The converse is not trivially true. By momentum conservation alone, the induced emission of a molecule with momentum $`𝐪`$ need not return two atoms into the condensate; but, in a process we call rogue photodissociation , the two atoms may emerge with any opposite nonzero momenta. However, the photodissociation processes that return the atoms back to the condensate are favored by Bose enhancement, and rogue photodissociation is further suppressed by energy conservation. Thus we adopt a two-mode model, only taking into account atoms with zero momentum and molecules with momentum $`𝐪`$. The corresponding creation and annihilation operators are denoted by $`a`$ and $`b`$.

### 5.2 The return of STIRAP

Turning now towards two-color photoassociation of a degenerate gas , we reconsider the possibility of STIRAP . It is thus assumed that a further laser beam couples the primarily photoassociated molecule to another bound molecular state, whose annihilation operator is denoted by $`g`$. The three-mode Hamiltonian reads

$$\frac{H}{\mathrm{}}=\mathrm{\Delta }g^{}g\delta b^{}b\frac{1}{2}\kappa (b^{}aa+ba^{}a^{})\frac{1}{2}\mathrm{\Omega }(b^{}g+bg^{}).$$
(2)
Here $`\mathrm{\Delta }`$ and $`\delta `$ are the two-photon and intermediate detunings, including the proper photon recoil energies, the free-bound QC Rabi frequency is

$$\kappa =\frac{𝐝𝐄}{2\mathrm{}}$$
(3)
and $`\mathrm{\Omega }`$ is the bound-bound Rabi frequency.

Given the Hamiltonian (2), the Heisenberg equations of motion for the boson operators read

$`\dot{a}`$ $`=`$ $`i\kappa a^{}b,`$ (4)
$`\dot{b}`$ $`=`$ $`i\delta b+\frac{1}{2}i(\kappa aa+\mathrm{\Omega }g),`$ (5)
$`\dot{g}`$ $`=`$ $`i\mathrm{\Delta }g+\frac{1}{2}i\mathrm{\Omega }b.`$ (6)
Suppose now that, if all molecules were dissociated to atoms, there were $`N`$ atoms. The boson operators in the system are then of the order $`\sqrt{N}`$. We thus introduce the rescaled boson operators $`\alpha =a/\sqrt{N}`$, …. These operators are of the order of unity, and may roughly be interpreted as the second-quantized counterparts of the probability amplitudes that an atom is in the atomic condensate ($`\alpha `$) or in one of the two molecular condensates ($`\beta `$, $`\gamma `$). The scaled operators obey the equation of motion

$`\dot{\alpha }`$ $`=`$ $`i\chi \alpha ^{}\beta ,`$ (7)
$`\dot{\beta }`$ $`=`$ $`i\delta \beta +\frac{1}{2}i(\chi \alpha \alpha +\mathrm{\Omega }\gamma ),`$ (8)
$`\dot{\gamma }`$ $`=`$ $`i\mathrm{\Delta }\gamma +\frac{1}{2}i\mathrm{\Omega }\beta .`$ (9)
The key point of our argument emerges from an inspection of the new Rabi frequency after scaling,

$$\chi =\frac{\sqrt{N}𝐝𝐄}{2\mathrm{}}.$$
(10)
The $`\sqrt{N}`$ is nothing but Bose enhancement in the present context. In the thermodynamic limit $`𝐝\sqrt{1/V}`$, so that $`\chi \sqrt{N/V}=\sqrt{\rho }`$ remains finite. More precisely, suppose that the same laser field with amplitude $`𝐄`$ were tuned in such a way that the the photodissociation of bound molecules would produce two atoms with the reduced mass $`\mu `$ and relative velocity $`v`$ at the rate $`\mathrm{\Gamma }(v)`$, then we have

$$\chi =\underset{v0}{lim}\sqrt{\frac{2\pi \mathrm{}^2\mathrm{\Gamma }(v)\varrho }{\mu ^2v}}.$$
(11)
The limit is finite and nonzero by virtue of the Wigner threshold law.

The couplings $`\chi `$ and $`\mathrm{\Omega }`$ in Eqs. (7)-(9) do not depend on the quantization volume anymore; the dependence on $`V`$ is replaced by a dependence on $`\varrho `$, the density of atoms if all molecules were to dissociate. In a BEC there is no longer any intrinsic restriction on the relative size of the couplings, and photoassociative STIRAP is feasible even in the thermodynamic limit.

We have constructed explicit examples of STIRAP by solving Eqs. (7)-(9) in a semiclassical or mean-field  approximation, treating $`\alpha `$, $`\beta `$ and $`\gamma `$ as $`c`$-numbers instead of quantum operators . Unlike the equations for probability amplitudes for the ordinary $`\mathrm{\Lambda }`$ system, these equations are nonlinear. Nonetheless, the basic character of STIRAP remains intact. In particular, for completely adiabatic switching of pulse strengths (for which Ref.  gives quantitative criteria), there is never any probability in the primary photoassociated state. Spontaneous-emission losses and rogue photodissociation are then shut off.

We cannot think of any matter of principle that could go catastrophically wrong with STIRAP. In practical experiments the atoms and possibly also the molecules are trapped and do not form infinite homogeneous condensates. Trapped particles come with a time scale that is of the order of tens of milliseconds in the magnetic case. If the STIRAP pulses are faster, trapping should not matter a whole lot. We have also ignored collisions between atoms, between molecules, and between atoms and molecules. But again, the corresponding time scales could easily be much longer than the time scale of the laser pulses, so that collisions are also negligible. The worst practical enemy of STIRAP might be light shifts of the two-photon resonance, which come about as a result of virtual transitions accompanying rogue photodissociation . However, in principle one can always compensate for light shifts by chirping the laser pulses.

## 6 Ruminations

Our renunciation of photoassociative STIRAP in a thermal nondegenerate gas most likely applies to other coherent optical transients as well . In the limit $`V\mathrm{}`$ the free-bound Rabi frequency $`\kappa 1/\sqrt{V}`$ tends to zero, and the corresponding time scale for the transient $`1/\kappa `$ tends to infinity. At least in a large enough sample, something else is likely to happen before a coherent transient can run its course.

Similarly, coherent transients should work in photoassociation of a BEC. Because of Bose enhancement, the coupling matrix element picks up a the factor $`\sqrt{N}`$, giving the effective Rabi frequency $`\chi =\sqrt{N}\kappa \sqrt{N/V}=\sqrt{\varrho }`$. In the thermodynamic limit $`\chi `$ remains finite, and the time scale of transients actually scales with atom density as $`\rho ^{1/2}`$. Such transients have already been discussed in the photoassociation literature , and equivalent mathematics is involved in a wide variety of topics ranging from second-harmonic generation  to the Feshbach resonance .

One final point concludes our discussion of transients in photoassociation. Bose enhancement oftentimes makes little practical difference even when it is formally present. For instance, suppose that we were to calculate a transition rate for a condensate under a normal bilinear coupling of the form $`b^{}a`$. Bose enhancement gives the factor $`\sqrt{N}`$ in the matrix element, and the transition rate for the condensate in terms of the one-atom rate $`R`$ reads $`NR`$. However, in the second quantized formulation all $`N`$ condensate atoms are treated at once, so that the transition rate per atom is still $`R`$. On the other hand, the coupling for photoassociation, $`b^{}aa`$, is trilinear. The ensuing Bose enhancement $`(\sqrt{N})^2`$ is thus stronger than for a bilinear coupling; in fact, just strong enough to offset the decrease of the coupling matrix element in the thermodynamic limit. As a result of the confluence of BEC and cubic coupling, a condensate responds to photoassociation as a unit on a time scale that depends on density, but oddly enough not directly on atom number or size of the condensate.

## 7 Acknowledgements

One of us (JJ) specifically thanks for Prof. Fritz Ehlotky for the opportunity to participate in the Obergurgl/Kühtai meetings already for the fifth time — and hopefully running…. This work is supported in part by NSF, Grant No. PHY-9801888, and by NASA, Grant No. NAG8-1428.


text:
# Improvements of the Discrete Dipole Approximation method

## Abstract

We report improvements in complex conjugate gradient algorithms applied to the discrete dipole approximation (DDA). It is shown that computational time is reduced by using the Bi-CGSTAB version of the CG algorithm, with diagonal left preconditioning.

Key words: scattering, non-spherical particles, discrete dipole approximation.

Optics Letters 1997, volume 22, number 16, 1205-1207.
© Optical Society of America, 1997.

The discrete-dipole approximation (DDA) is a flexible technique for computing scattering and absorption by targets of arbitrary geometry. In  the discrete dipole approximation (DDA) for scattering calculations is reviewed. Rather than “direct” methods for solving linear system of equations arising in DDA problem iterative methods for finding the solution have proven effective and efficient.

In this paper we perform systematic study of various non-stationary iterative (conjugate gradient) methods in search for the most efficient one. We document implementation of these methods in our public domain code DDSCAT.5a code

Numerical aspects of the discrete dipole approximation continue to be of great interest. Yung  applied a conjugate gradient method to the in DDA approach. Hoekstra  identifies Yung’s scheme as the conjugate gradient (CG) algorithm proposed by Hestenes . Rahola  discusses solution of dense systems of linear equations in the discrete-dipole approximation and choice of of the best iterative method in this application. Draine  implemented a conjugate gradient method based on work of Petravic and Kuo-Petravic.  This implementation is quite robust and has been used for many years.  However, Lumme and Rahola  applied the quasi-minimal residual (QMR) conjugate gradient algorithm to the system of linear equations arising in the DDA applications. They claim that the QMR method is approximately 3 times faster in comparison to the one employed in the DDSCAT code.  It was this work which prompted us to perform the analysis reported here.

PIM is a collection of Fortran 77 routines designed to solve systems of linear equations on parallel and sequential computers using a variety of iterative methods.

PIM contains implementations of various methods: conjugate-gradient (CG); Conjugate-Gradients for normal equations with minimization of the residual norm (CGNR); Conjugate-Gradients for normal equations with minimization of the error norm (CGNE); Bi-Conjugate-Gradients (Bi-CG); Conjugate-Gradients squared (CGS); the stabilised version of Bi-Conjugate-Gradients (Bi-CGSTAB); the restarted, stabilised version of Bi-Conjugate-Gradients (RBi-CGSTAB); the restarted, generalized minimal residual (RGMRES); the restarted, generalized conjugate residual (RGCR), the quasi-minimal residual with coupled two-term recurrences (QMR); the transpose-free quasi-minimal residual (TFQMR); and Chebyshev acceleration. The routines allow the use of preconditioners; the user may choose to use left-, right- or symmetric-preconditioning.

The convergence rate of iterative methods depends on the coefficient matrix. Hence one may attempt to transform the linear system into one that is equivalent (in the sense that it has the same solution) but is easier to solve. A preconditioner is a matrix $`M`$ that effects such a transformation. It is possible to introduce left- and right preconditioners. The simplest preconditioner consists of just the diagonal of the coefficient matrix. This is known as the (point) Jacobi preconditioner.

To compare these different algorithms we have used them to find solutions to the problem of scattering by a homogeneous sphere. The scattering problem is specified by the usual size parameter $`x=2\pi a/\lambda `$, where a is the radius.

Tables 1 and 2 presents the number of iterations and CPU time for size parameter $`x=0.1`$ and $`x=1`$ and for several values of refractive index. The conjugate gradient methods are defined as above. Label (L) indicates left Jacobi preconditioning. For example CGNE(L) is the conjugate gradient method for normal equations with minimization of the error norm and left Jacobi preconditioning. Similarly, (R) indicates right Jacobi preconditioning. CPU time (sequential Silicon Graphics workstation) is normalized to the “best” method. Star indicates that the method did not converge in the maximum allowed number of iterations or that the method failed to converge. Fractional error $`10^5`$ was used as the stopping criterion. The DDSCAT.5a code with the newly implemented GPFA fast Fourier transform method was used. For Bi-CGSTAB and CGNE we used left and right Neumann polynomial preconditioner truncated after the first term. Thus, Bi-CGSTAB(N)(L) indicates the stabilised version of Bi-Conjugate-Gradients method with left Neumann polynomial preconditioner.

Table 1 presents results for size parameter $`x=0.1`$ and real refractive index $`\mathrm{n}=1.33,2,3,5`$ as well as one case with small complex part of refractive index $`\mathrm{n}=(5,0.0001)`$ and size parameter $`\mathrm{x}=0.1`$. In Table 1 the CPU times are normalized to the CG(L) method, which was found to be the best method. For example it is 4.0 times faster in comparison with the CGNE for $`\mathrm{n}=(1.33,0)`$. For larger values of real refractive index the CGNE is almost an order of magnitude slower in comparison to CG. This is because more iterations are needed for the same convergence and because cost of one CG iteration is less than cost of one CGNE iteration. The QMR algorithm is never competitive and actually fails to converge for large real refractive indices. For small refractive index the Bi-CGSTAB algorithm is comparable to the CG and requires less iterations. However, the cost per iteration is larger in comparison to CG which offsets the advantage of lesser number of iterations. The Petravic and Kuo-Petravic  algorithm used by us for many years  is similar to CGNR and CGNE. However, we observed on occasion slightly different convergence rates due to stabilization of Petravic and Kuo-Petravic algorithm every 10th time step.  This is true for all other cases. The storage requirements of CG, CGNE, CGNR is $`6\times N`$, for BiCG it is $`8\times N`$, for CGS, Bi-CGSTAB, TFQMR it is $`10\times N`$, QMR requires $`11\times N`$. Thus, for pure real refractive index, the CG is not only the fastest method but also it requires the least amount of temporary storage. It can be seen that left preconditioning by the inverse of diagonal of the DDA matrix  reduces the number of iterations needed. The added time needed for division by diagonal elements is generally negligible in comparison to the time saved by smaller amount of iterations. It can be seen that for Bi-CGSTAB, Bi-CGSTAB(L), and Bi-CGSTAB(R) the left Jacobi preconditioning is the only method converging for larger refractive index. Restarted methods (RBi-CBSTAB and RGCR) appear to be not competitive but further study may be needed (we used the orthogonal base of 10 vectors for all restarted methods). The CG method is also competitive in cases with small absorption (see last column of Table 1). We have also calculated (not presented) results for size parameter of $`x=0.1`$ and increasing complex part of refractive index $`n=(1.33,0),(1.33,0.01),(1.33,0.1),(1.33,1),(1.33,2),(1.33,3)`$. The BiCGSTAB(L), which proved to be the most robust method. However the CGS(L) is competitive and faster for $`n=(1.33,3)`$. Both CGS(L) and BiCGSTAB(L) require the same amount of iteration for convergence and their cost is similar. These methods are between 2.9 and 1.6 times faster in comparison to CGNR — the method used in DDSCAT code. The QMR and TFQMR which Lumme and Rahola  claim to be faster in comparison to CGNR and the DDSCAT implementations do not converge on occasion and when they work they are only slightly better in this case. As before, left Jacobi preconditioning is almost always beneficial. The CG(L) algorithm is faster than BiCGSTAB(L) for refractive index $`n=(1.33,0),(1.33,0.01),(1.33,0.1)`$.

Table 2 is for size parameter $`x=1`$. All the results are normalized to Bi-CGSTAB(L). This method is clearly superior to the CGNR method and it is 2-4.3 faster. It can be seen that CGNR converges slowly, and has not satisfied the stopping criterion in 140 iterations for $`n=(3,0.0001)`$. For this larger value of size parameter the QMR algorithm doesn’t converge well but its smooth version TFQMR does. However, TFQMR is slower in comparison to Bi-CGSTAB(L) and comparable to CGNR. The CG(L) method for refractive index $`n=(1.33,0)`$ and $`n=1.33,0.01`$ is faster than the reference scheme Bi-CGSTAB(L). It can be seen that the Neumann polynomial preconditioning Bi-CGSTAB(N)(L) or Bi-CGSTAB(N)(R) does reduce the number of iterations needed for certain cases of refractive index. However the cost associated with the additional calculations always offsets this improved convergence rate. As before, the left Jacobi preconditioner is superior to right or no-preconditioner cases. CG(L) works well for small refractive index but is comparable to Bi-CGSTAB(L). The QMR algorithm fails to converge but the transpose-free quasi-minimal residual (TFQMR) algorithm converges well and is comparable to CGNR. The CG method is theoretically valid for Hermitian positive definite matrices. The matrix arising in the DDA is not Hermitian but symmetric. Therefore, strictly speaking, the CG method is not valid for use in the DDA. The users are advised to test the CG method when extrapolating results presented here to different size parameters, particle sizes, and refractive index values.

We recommend use of the stabilized version of the Bi-conjugate gradient algorithm with left Jacobi preconditioning \[Bi-CGSTAB(L)\]. This algorithms requires 67% greater storage than the CGNR algorithm, but is typically 2-3 times faster.

The recent version of Discrete Dipole Approximation code DDSCAT5a developed by Draine and Flatau contains improvements documented in this paper. The code is available via anonymous ftp from the `ftp.astro.princeton.edu` site or from the Light Scattering and Radiative Transfer Codes Library — `SCATTERLIB` (`http://atol.ucsd.edu/~pflatau`).

I have been supported in part by the Office of Naval Research Young Investigator Program and in part by DuPont Corporate Educational Assistance. I would like to thank Drs M. J. Wolff and A. E. Ilin who helped with computer tests. Bruce Draine checked the manuscript. Dr. R. J. Riegert of Du Pont if acknowledged for his continuing interest in DDSCAT developments.