format: text

text:
# Time Evolution of Two-Level Systems Driven by Periodic Fields

## I Introduction

Periodically (or more generically quasi-periodically) driven quantum two-level systems are of basic importance in many physical applications an have been widely studied since the pioneering works of Rabi , of Bloch and Siegert  and of Autler and Townes  (see also  for more recent discussions and  for other general references on the subject). The time evolution of such systems is governed by the Schrödinger equation (we henceforth adopt $`\mathrm{}=1`$)

$$i\frac{d}{dt}|\mathrm{\Psi }=H(t)|\mathrm{\Psi },$$
(1)
where $`|\mathrm{\Psi }=|\mathrm{\Psi }(t)=\left(\genfrac{}{}{0pt}{}{\psi _+(t)}{\psi _{}(t)}\right)`$ and $`H(t)`$ is the Hamiltonian of the system. We may consider, for instance, a spin-1/2 system under the influence of a time-dependent (periodic) magnetic field $`\stackrel{}{B}(t)`$. In this situation, $`H(t)`$ takes the usual form $`H(t)=\frac{1}{2}\stackrel{}{B}(t)\stackrel{}{\sigma }`$, where $`\stackrel{}{\sigma }=(\sigma _1,\sigma _2,\sigma _3)`$ are the Pauli matrices. The interest in the solution of (1) in this case is not restricted to the investigation of the quantum system. As first pointed by Feynman, Vernon and Hellwarth  (see also the discussion in ), the quantum system is equivalent to the classical Hamiltonian system describing a classical gyromagnet precessing in a magnetic field: $`\frac{d}{dt}\stackrel{}{𝒮}=\stackrel{}{B}(t)\times \stackrel{}{𝒮}`$, where $`\stackrel{}{𝒮}`$ is a (three dimensional) unit vector. In fact, the methods described below can be directly applied to the analysis of this system as well, since the components of $`\stackrel{}{𝒮}`$ can be written in terms of the components $`\psi _\pm (t)`$ of $`|\mathrm{\Psi }`$ .

The evolution of the systems governed by (1) is determined by the time evolution operator $`U(t,s)`$ (also known as the propagator) which connects the state $`|\mathrm{\Psi }(s)`$ at time $`s`$ to the state $`|\mathrm{\Psi }(t)`$ at time $`t`$: $`|\mathrm{\Psi }(t)=U(t,s)|\mathrm{\Psi }(s)`$. Defining $`U(t)=U(t,0)`$ one has $`U(t,s)=U(t)U(s)^{}`$. For a time-dependent Hamiltonian the propagator $`U(t)`$ can be computed via the Dyson expansion:

$$U(t)=\widehat{1}+\underset{n=1}{\overset{\mathrm{}}{}}(i)^n_0^tH(t_1)𝑑t_1\mathrm{}_0^{t_{n1}}H(t_n)𝑑t_n.$$
(2)
Although (2) gives a straightforward manner to compute $`U(t)`$, the series in the r.h.s. is not generally uniformly convergent in time. For practical purposes this gives rise to difficulties when one is interested in the large-time behaviour of the system. For instance, if one considers a periodic Hamiltonian of the form $`H(t)=_mH_me^{im\omega t}`$, with $`H_00`$, two successive integrations in (2) would produce a linear term in $`t`$. Higher order terms in $`t`$ would appear with further integrations. These polynomial terms are known as secular terms and they plague the expansion of $`U(t)`$ in such a way that its uniform convergence is spoiled.

Of particular interest is the situation where the Schrödinger equation (1) takes the form

$$i\frac{d}{dt}|\mathrm{\Psi }=H_1(t)|\mathrm{\Psi },\text{with}H_1(t):=ϵ\sigma _3f(t)\sigma _1,$$
(3)
where $`f(t)`$ is a periodic function of time $`t`$ and $`ϵ`$ is a real constant. By a rotation of $`\pi /2`$ around the 2-axis, we get the equivalent system

$$i\frac{d}{dt}|\mathrm{\Phi }=H_2(t)|\mathrm{\Phi },\text{with}H_2(t):=ϵ\sigma _1+f(t)\sigma _3,$$
(4)
where $`|\mathrm{\Phi }=|\mathrm{\Phi }(t)=\mathrm{exp}(i\pi \sigma _2/4)|\mathrm{\Psi }(t)`$. The Hamiltonian in (3) can be interpreted as describing a system with an unperturbed diagonal Hamiltonian $`H_0:=ϵ\sigma _3`$, subjected to a (periodic) time-dependent perturbation $`H_I(t):=f(t)\sigma _1`$, inducing a time-depending transition between the unperturbed eigenstates of $`H_0`$. Of course, (3) can be also interpreted as a spin-1/2 system under a magnetic field $`\stackrel{}{B}(t)=(2f(t),0,2ϵ)`$.

In  we investigated the system described by (3) or (4) in the situation where $`f`$ is a periodic or quasi-periodic function of time and $`ϵ`$ is “small” and a special perturbative expansion (power series in $`ϵ`$) was developed, whose main virtue is to be free of secular terms. The algorithm employed involves an inductive “renormalisation” of a sort of effective field introduced through an exponential Ansatz (the function $`g`$, to be introduced below). Moreover, in the periodic case, we have been able in  to prove rigorously the uniform convergence of our expansions for all $`t`$, for $`ϵ`$ sufficiently small, provided the function $`f`$ satisfies the special conditions I, II and III we describe below.

A feature of our method is the fact that we are able to present the complete (and convergent) $`ϵ`$-expansion for the secular frequency $`\mathrm{\Omega }`$ (also known as Rabi frequency) associated to the solution of (3)-(4). This is particularly important for the qualitative study of the large-time behaviour of that solutions. Another relevant conclusion is that, under conditions I-III, our method provides a complete perturbative solution of (3)-(4) for the monochromatic field (ad-dc field), except perhaps for spurious situations .

Since our series are uniformly convergent in time, it is possible to use them to study the long-time behaviour of observable quantities in a controlled way. It turns out that our expansions are also very practical for numerical calculations, providing very small errors even with few terms.

In this paper we describe the algorithm employed in the numerical computation of our perturbative solution and show the results obtained in the particular case of the monochromatic interactions. We begin in Section II by given a quick review of the perturbative method developed in . Section III presents a description of the numerical algorithm developed to compute the expansions. Next, in Section IV, we show and discuss the results obtained in the case of monochromatic interactions. Finally, in Section V, we draw some conclusions and final remarks.

## II Description of the Method

In , it was shown that the propagator $`U(t)`$ associated with the system (4) can be written as

$$U(t)=\left(\begin{array}{cc}R(t)(1+ig_0S(t))& iϵR(t)S(t)\\ iϵ\overline{R(t)}\overline{S(t)}& \overline{R(t)}(1i\overline{g_0}\overline{S(t)})\end{array}\right),$$
(5)
where

$$R(t):=\mathrm{exp}\left(i_0^t(f(\tau )+g(\tau ))𝑑\tau \right)$$
(6)
and

$$S(t):=_0^tR(\tau )^2𝑑\tau ,$$
(7)
where $`g`$ is a particular solution of the generalised Riccati equation

$$g(t)^{}ig(t)^22if(t)g(t)+iϵ^2=0,$$
(8)
and $`g_0g(0)`$. To solve (8) we pose

$$g(t)=\underset{n=1}{\overset{\mathrm{}}{}}G^{(n)}(t)ϵ^n,$$
(9)
where $`G^{(n)}(t):=q(t)c_n(t)`$, with

$$q(t):=\mathrm{exp}\left(i_0^tf(\tau )𝑑\tau \right).$$
(10)
Inserting (9) into (8) yields a sequence of recursive equations for the coefficients $`c_n`$, whose solutions are

$`c_1(t)`$ $`=`$ $`\alpha _1q(t),`$ (11)
$`c_2(t)`$ $`=`$ $`q(t)\left[i{\displaystyle _0^t}(\alpha _1^2q(\tau )^2q(\tau )^2)𝑑\tau +\alpha _2\right],`$ (12)
$`c_n(t)`$ $`=`$ $`q(t)\left[i\left({\displaystyle _0^t}{\displaystyle \underset{p=1}{\overset{n1}{}}}c_p(\tau )c_{np}(\tau )d\tau \right)+\alpha _n\right],`$ (13)
for $`n3`$, where the $`\alpha _n`$ are arbitrary integration constants. Our strategy consists in fixing these constants inductively to cancel the secular terms. For instance, to cancel the secular term in $`c_2`$, the integrand in (12) cannot contain a constant term. Defining the mean value of a quasi-periodic function $`h`$ by

$$M(h):=\underset{T\mathrm{}}{lim}\frac{1}{2T}_T^Th(t)𝑑t,$$
(14)
we must require that $`M(\alpha _1^2q^2q^2)=0`$. This leads to $`\alpha _1^2=\overline{M(q^2)}/M(q^2)`$, provided $`M(q^2)0`$. We remark that the mean value of a quasi-periodic function $`h`$ equals its zero order Fourier coefficient. In  we identified three classes of quasi-periodic functions $`f`$ for which this procedure of elimination of secular terms can be applied to all orders (i.e., to all functions $`c_n`$). These classes are defined by the following conditions:

$`\mathrm{I}.`$ $`M(𝒬_0)0,`$ (15)
$`\mathrm{II}.`$ $`M(𝒬_0)=0\text{but}M(𝒬_1)0,`$ (16)
$`\mathrm{III}.`$ $`M(𝒬_0)=M(𝒬_1)=0\text{but}M(𝒬_3)0,`$ (17)
where $`𝒬_0(t):=q(t)^2`$,

$`𝒬_1(t)`$ $`:=`$ $`𝒬_0(t){\displaystyle _0^t}\left(𝒬_0(\tau )^1M\left(𝒬_0^1\right)\right)𝑑\tau ,`$ (18)
$`𝒬_3(t)`$ $`:=`$ $`𝒬_0(t){\displaystyle _0^t}(𝒬_1(\tau )M(𝒬_1))𝑑\tau .`$ (19)
Under condition II, the Ansatz (9) has to be slightly modified to

$$g(t)=\underset{n=1}{\overset{\mathrm{}}{}}\stackrel{~}{G}^{(n)}(t)\lambda ^n,\text{with}\lambda :=ϵ^2$$
(20)
and $`\stackrel{~}{G}^{(n)}(t):=q(t)e_n(t)`$. The solutions for $`e_n`$ are

$`e_1(t)`$ $`=`$ $`q(t)\left(i{\displaystyle _0^t}q(\tau )^2𝑑\tau +\beta _1\right),`$ (21)
$`e_n(t)`$ $`=`$ $`q(t)\left[i\left({\displaystyle _0^t}{\displaystyle \underset{p=1}{\overset{n1}{}}}e_p(\tau )e_{np}(\tau )d\tau \right)+\beta _n\right],`$ (22)
for $`n2`$. Under condition III, we found in  the same solution for $`g`$ as in I, namely (9). The difference between the solutions obtained in I and III lie in the set of constants $`\alpha _n`$, which is conventionally fixed to remove the secular terms in each case.

In  we have shown that the $`ϵ`$-expansion (9) or (20) converges uniformly for $`|ϵ|`$ small and for all $`t`$, provided $`f`$ is periodic and satisfies I or II. In  we extended this result for case III. As we will discuss soon, case III is particularly important for the investigation of the dynamical localisation effect (also called, less properly, coherent destruction of tunnelling) for monochromatic interactions.

The secular frequency $`\mathrm{\Omega }`$, in conditions I or III, is given by

$$\mathrm{\Omega }=M(f)+M(g)=F_0+\underset{n=1}{\overset{\mathrm{}}{}}M\left(G^{(n)}\right)ϵ^n$$
(23)
and, in condition II, by

$$\mathrm{\Omega }=F_0+\underset{n=1}{\overset{\mathrm{}}{}}M\left(\stackrel{~}{G}^{(n)}\right)\lambda ^n.$$
(24)
Hence, with our previous definitions, for condition I, we have

$`\mathrm{\Omega }`$ $`=`$ $`F_0+ϵ\alpha _1M(𝒬_0)`$ (25)
$`+`$ $`ϵ^2\left[i\alpha _1^2M(𝒬_2)iM(𝒬_1)+\alpha _2M(𝒬_0)\right]`$ (26)
$`+`$ $`ϵ^3\left[2\alpha _1M(𝒬_3)+\alpha _3M(𝒬_0)\right]+𝒪(ϵ^4),`$ (27)
where

$$𝒬_2(t):=𝒬_0(t)_0^t(𝒬_0(\tau )M(𝒬_0))𝑑\tau .$$
(28)
One easily sees, by computing the zero order Fourier coefficient of $`𝒬_2`$ (details in ), that $`M(𝒬_2)=0`$ whenever $`M(𝒬_0)=0`$. Hence, for condition II, we get

$$\mathrm{\Omega }=F_0iϵ^2M(𝒬_1)+𝒪(ϵ^4)$$
(29)
and, for condition III,

$$\mathrm{\Omega }=F_0+2\alpha _1ϵ^3M(𝒬_3)+𝒪(ϵ^4).$$
(30)
## III Numerical Computation of the Solution

We present here a step-by-step algorithm which enables us to compute the propagator (5) solely from the Fourier coefficients of the (periodic) interaction $`f`$. This algorithm is based on the convergent expansions presented in , which expresses the propagator $`U(t)`$ in terms of its Floquet form.

Let us first begin with some conventions. We suppose that the interaction function $`f`$ is periodic with period $`T_\omega =2\pi /\omega `$ and that its Fourier decomposition $`f=_nF_ne^{in\omega t}`$ has only a finite number of terms. Since $`f`$ is real, excluding the constant term $`F_0`$, $`f`$ must have an even number of non-vanishing Fourier coefficients, say $`2J`$, with $`J1`$. Denoting the set of integers $`\{n,n0|F_n0\}`$ by $`\{n_1,\mathrm{},n_{2J}\}`$, we may write

$$f(t)=F_0+\underset{a=1}{\overset{2J}{}}f_ae^{in_a\omega t},$$
(31)
with the convention that $`n_a=n_{2Ja+1}`$, for all $`1aJ`$, and with $`f_aF_{n_a}`$. Clearly $`\overline{f_a}=f_{2Ja+1}`$, $`1aJ`$.

### A Interactions with $`F_0=0`$

In this case, the Fourier decomposition of the functions $`q`$ and $`q^2`$ can be written as

$$q(t)=\underset{m}{}Q_me^{im\omega t}\text{and}q(t)^2=\underset{m}{}Q_m^{(2)}e^{im\omega t}.$$
(32)
The coefficients $`Q_m`$ and $`Q_m^{(2)}`$ are of basic importance in the numerical computation of the propagator (5). They can be obtained in a closed form from the coefficients $`F_m`$ of $`f`$. Explicitly, we have

$$Q_m=e^{i\gamma _f}\underset{p_1,\mathrm{},p_{2J}=0}{\overset{\mathrm{}}{}}\delta (P,m)\underset{a=1}{\overset{2J}{}}\left[\frac{1}{p_a!}\left(\frac{f_a}{n_a\omega }\right)^{p_a}\right],$$
(33)
for all $`m`$, where

$$PP(p_1,\mathrm{},p_{2J},n_1,\mathrm{},n_{2J}):=\underset{b=1}{\overset{2J}{}}p_bn_b,$$
(34)
and

$$\gamma _f:=i\underset{a=1}{\overset{2J}{}}\frac{f_a}{n_a\omega }.$$
(35)
The symbol $`\delta (P,m)`$ denotes the Krönecker delta: $`\delta (P,m)=1`$, if $`P=m`$ and $`\delta (P,m)=0`$, otherwise. To compute the Fourier coefficients of $`q^2`$, we simply note that $`q^2`$ is obtained from $`q`$ by the substitution $`f2f`$. Hence, for all $`m`$,

$$Q_m^{(2)}=e^{2i\gamma _f}\underset{p_1,\mathrm{},p_{2J}=0}{\overset{\mathrm{}}{}}\delta (P,m)\underset{a=1}{\overset{2J}{}}\left[\frac{1}{p_a!}\left(\frac{2f_a}{n_a\omega }\right)^{p_a}\right].$$
(36)
Formulas (33) and (36) can be computed either analytically (in some cases) or numerically. For the monochromatic interactions, a closed form in terms of Bessel functions of first kind is obtained. The exact result will be presented in Section IV.

Once we know the coefficients $`Q_m`$ and $`Q_m^{(2)}`$, we proceed to compute the particular solution (9) or (20) of the generalised Riccati equation (8). To decide whether we use (9) or (20), we must check which condition I, II or III is satisfied by $`f`$. Hence, we have to look at the mean values of $`𝒬_0`$, $`𝒬_1`$ and $`𝒬_3`$. One obviously has $`M(𝒬_0)=M(q^2)=Q_0^{(2)}`$. To obtain $`M(𝒬_1)`$ and $`M(𝒬_3)`$, we need the Fourier decompositions of $`𝒬_1`$ and $`𝒬_3`$. From the definitions (18) and (19), after some simple computations, we get

$$M(𝒬_1)=\frac{i}{\omega }\underset{\genfrac{}{}{0pt}{}{m}{m0}}{}\frac{\overline{Q_m^{(2)}}\left(Q_0^{(2)}Q_m^{(2)}\right)}{m}$$
(37)
and

$`M(𝒬_3)`$ $`=`$ $`{\displaystyle \frac{1}{\omega ^2}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{n,m}{n0,m0}}{}}{\displaystyle \frac{\overline{Q_m^{(2)}}}{nm}}(Q_0^{(2)}Q_n^{(2)}`$ (38)
$``$ $`Q_0^{(2)}Q_{nm}^{(2)}+Q_n^{(2)}Q_{nm}^{(2)}).`$ (39)
The numerical value of $`M(𝒬_1)`$ and $`M(𝒬_3)`$ can be calculated trivially from the above expressions once $`Q_m^{(2)}`$ are known.

#### 1 Computing $`g(t)`$ in Cases I and III

We remember that condition I applies whenever $`f`$ satisfies $`M(𝒬_0)0`$. In this case, by properly fixing the constants $`\alpha _n`$, we can completely eliminate the secular terms from the functions $`c_n`$ . For this reason, we may write

$$c_n(t)=\underset{m}{}C_m^{(n)}e^{im\omega t}.$$
(40)
The Fourier coefficients $`C_m^{(n)}`$ are obtained from equations (11)-(13). Their inductive structure is given by the relations

$`C_m^{(1)}`$ $`=`$ $`\alpha _1Q_m,`$ (41)
$`C_m^{(2)}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{n_1}{n_10}}{}}{\displaystyle \frac{Q_{mn_1}\left(\alpha _1^2Q_{n_1}^{(2)}\overline{Q_{n_1}^{(2)}}\right)}{n_1\omega }}+\alpha _2Q_m,`$ (42)
$`C_m^{(n)}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{n_1,n_2}{n_1+n_20}}{}}{\displaystyle \frac{Q_{m(n_1+n_2)}}{(n_1+n_2)\omega }}\left({\displaystyle \underset{p=1}{\overset{n1}{}}}C_{n_1}^{(p)}C_{n_2}^{(np)}\right)`$ (43)
$`+`$ $`\alpha _nQ_m,\text{for }n3.`$ (44)
The constants $`\alpha _n`$ have closed forms in terms of the coefficients $`Q_m^{(2)}`$ and $`C_m^{(p)}`$, for $`pn1`$. Since they involve somewhat large expressions, we refrain from writing them here. The complete expressions can be found in . We may see that the whole inductive structure of the coefficients $`C_m^{(n)}`$ is known and, therefore, the computation of $`c_n`$ is just a matter of numerically evaluating the above expressions.

Since in case I $`g`$ is given by (9), we may write

$$g(t)=q(t)\underset{n=1}{\overset{\mathrm{}}{}}c_n(t)ϵ^n=:\underset{n=1}{\overset{\mathrm{}}{}}\left(\underset{m}{}G_m^{(n)}e^{im\omega t}\right)ϵ^n,$$
(45)
where $`G_m^{(n)}`$ is given by the convolution

$$G_m^{(n)}=\underset{p}{}Q_{mp}C_p^{(n)},$$
(46)
whose numerical value can be easily computed since we already know $`Q_m`$ and $`C_m^{(n)}`$, for all $`m`$ and $`n`$. This gives $`g`$ for all $`t`$ under condition I. For condition III ($`M(𝒬_0)=M(𝒬_1)=0`$, but $`M(𝒬_3)0`$) one has essentially the same solution except for the constants $`\alpha _n`$ which are calculated differently from condition I. Their formulas can be found in .

#### 2 Computing $`g(t)`$ in Case II

Condition II applies whenever $`f`$ satisfies $`M(𝒬_0)=0`$, but $`M(𝒬_1)0`$. In this case, the perturbative solution of the generalised Riccati equation (8) is given by (20). The constants $`\beta _n`$ which appear in the functions $`e_n`$ (equations (21)-(22)) can be chosen such that no secular terms emerge . Hence, with this particular choice of the constants $`\beta _n`$, we may write

$$e_n(t)=\underset{m}{}E_m^{(n)}e^{im\omega t}.$$
(47)
The recursive structure of the coefficients $`E_m^{(n)}`$ is given by

$`E_m^{(1)}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{n}{n0}}{}}{\displaystyle \frac{Q_{mn}\overline{Q_n^{(2)}}}{n\omega }}+Q_m\left(\beta _1+{\displaystyle \underset{\genfrac{}{}{0pt}{}{n}{n0}}{}}{\displaystyle \frac{\overline{Q_n^{(2)}}}{n\omega }}\right),`$ (48)
$`E_m^{(n)}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{n_1,n_2}{n_1+n_20}}{}}{\displaystyle \frac{Q_{mn_1n_2}}{(n_1+n_2)\omega }}\left({\displaystyle \underset{p=1}{\overset{n1}{}}}E_{n_1}^{(p)}E_{n_2}^{(np)}\right)`$ (49)
$``$ $`Q_m{\displaystyle \underset{\genfrac{}{}{0pt}{}{n_1,n_2}{n_1+n_20}}{}}\left({\displaystyle \underset{p=1}{\overset{n1}{}}}E_{n_1}^{(p)}E_{n_2}^{(np)}\right){\displaystyle \frac{1}{(n_1+n_2)\omega }}`$ (50)
$`+`$ $`\beta _nQ_m,\text{for }n2.`$ (51)
The constants $`\beta _n`$ have also closed forms  in terms of the known $`Q_m^{(2)}`$ and $`E_m^{(p)}`$, for all $`m`$ and $`pn1`$. Thus, the solution (20) can be evaluated numerically in any order for all $`t`$. For future convenience, we write

$$g(t)=q(t)\underset{n=1}{\overset{\mathrm{}}{}}e_n(t)\lambda ^n=:\underset{n=1}{\overset{\mathrm{}}{}}\left(\underset{m}{}\stackrel{~}{G}_m^{(n)}e^{im\omega t}\right)\lambda ^n,$$
(52)
where $`\stackrel{~}{G}_m^{(n)}`$ is given by

$$\stackrel{~}{G}_m^{(n)}=\underset{p}{}Q_{mp}E_p^{(n)},$$
(53)
which completely specifies $`g`$ for all $`t`$ under condition II.

#### 3 Computing the Propagator in Cases I, II and III

Once the coefficients $`G_m^{(n)}`$ (for cases I and III) or $`\stackrel{~}{G}_m^{(n)}`$ (for case II) are known, the propagator $`U(t)`$ expressed in (5) can be computed in a straightforward manner, as will be shown now. We illustrate our procedure with the coefficients $`G_m^{(n)}`$ of conditions I and III. For condition II, where the coefficients of $`g`$ are $`\stackrel{~}{G}_m^{(n)}`$, the discussed procedure has to be adapted with self-evident modifications.

We begin by defining

$$G_m(ϵ):=\underset{n=1}{\overset{\mathrm{}}{}}G_m^{(n)}ϵ^n.$$
(54)
The secular frequency (23) is clearly given by (see (45))

$$\mathrm{\Omega }\mathrm{\Omega }(ϵ)=\underset{n=1}{\overset{\mathrm{}}{}}G_0^{(n)}ϵ^n,$$
(55)
since we are supposing that $`F_0=0`$. Hence, by (45),

$$g(t)=\mathrm{\Omega }+\underset{\genfrac{}{}{0pt}{}{m}{m0}}{}G_m(ϵ)e^{im\omega t}.$$
(56)
Looking at expression (5) for the propagator, we see that the Fourier series of $`R`$ (see (6)) can be computed if we first find the Fourier decomposition of

$$W(t):=\mathrm{exp}\left(i_0^tg(\tau )𝑑\tau \right).$$
(57)
Indeed, since $`R=\overline{q}W`$, we obtain the Fourier series of $`R`$ by taking a convolution of the coefficients of $`\overline{q}`$ and $`W`$.

It is easy to see that

$$W(t)=e^{i\gamma _f(ϵ)}e^{i\mathrm{\Omega }t}\mathrm{exp}\left(\underset{m}{}H_me^{im\omega t}\right),$$
(58)
with

$$H_mH_m(ϵ):=\frac{G_m(ϵ)}{m\omega },\text{for }m0$$
(59)
and $`H_0=0`$. Moreover,

$$\gamma _f(ϵ):=i\underset{m}{}H_m.$$
(60)
Writing

$$W(t)=e^{i\mathrm{\Omega }t}\underset{m}{}W_me^{im\omega t},$$
(61)
and using (58), we find

$`W_mW_m(ϵ)=e^{i\gamma _f(ϵ)}(H_m+{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^{p+1}}{(p+1)!}}`$ (62)
$`\times {\displaystyle \underset{n_1,\mathrm{},n_p}{}}H_{n_1}\mathrm{}H_{n_p}H_{mN_p}),\text{for }m0`$ (63)
and

$`W_0W_0(ϵ)=e^{i\gamma _f(ϵ)}(1+{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^{p+1}}{(p+1)!}}`$ (64)
$`\times {\displaystyle \underset{n_1,\mathrm{},n_p}{}}H_{n_1}\mathrm{}H_{n_p}H_{N_p}),`$ (65)
with $`N_p:=_{a=1}^pn_a`$. Using now (61) and the fact that $`R=\overline{q}W`$, we conclude that $`R`$ can be written as

$$R(t)=e^{i\mathrm{\Omega }t}\underset{m}{}R_me^{im\omega t},$$
(66)
with the coefficients $`R_m`$ given by the convolution

$$R_m=\underset{p}{}\overline{Q_{pm}}W_p.$$
(67)
This finishes with the computation of $`R`$ in terms of its Fourier series (66). We note from formulas (54)-(67) that once the coefficients $`G_m^{(n)}`$ of $`g`$ are given, we can numerically evaluate $`G_m`$, $`H_m`$, $`W_m`$ and, hence, $`R_m`$ with a trivial computer code.

Next we proceed to compute the Fourier series of $`S`$ (see (7)). First we find the Fourier coefficients of $`R^2`$. This is an easy task since $`R^2`$ is obtained from $`R`$ by replacing $`(f+g)2(f+g)`$. Hence, we must replace $`H_m2H_m`$ and $`\overline{q}q^2`$. Consequently, we get

$$R(t)^2=e^{2i\mathrm{\Omega }t}\underset{m}{}R_m^{(2)}e^{im\omega t},$$
(68)
where

$$R_m^{(2)}=\underset{p}{}Q_{mp}^{(2)}W_p^{(2)}.$$
(69)
with

$`W_m^{(2)}W_m^{(2)}(ϵ)=e^{2i\gamma _f(ϵ)}(2H_m+{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{2^{p+1}}{(p+1)!}}`$ (70)
$`\times {\displaystyle \underset{n_1,\mathrm{},n_p}{}}H_{n_1}\mathrm{}H_{n_p}H_{mN_p}),\text{for }m0`$ (71)
and

$`W_0^{(2)}W_0^{(2)}(ϵ)=e^{2i\gamma _f(ϵ)}(1+{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{2^{p+1}}{(p+1)!}}`$ (72)
$`\times {\displaystyle \underset{n_1,\mathrm{},n_p}{}}H_{n_1}\mathrm{}H_{n_p}H_{N_p}).`$ (73)
Now $`S`$ is obtained by a simple integration of $`R^2`$. A trivial computation from (68), gives

$$S(t)=\sigma _0+e^{2i\mathrm{\Omega }t}\underset{m}{}S_me^{im\omega t},$$
(74)
with

$$S_m:=i\frac{R_m^{(2)}}{m\omega +2\mathrm{\Omega }}\text{and}\sigma _0:=\underset{m}{}S_m.$$
(75)
We assume that $`m\omega +2\mathrm{\Omega }0`$ for all $`m`$ (see the discussion of crossings in ).

We have found expressions for $`R`$ and $`S`$ in terms of its Fourier series. This series converge absolutely and uniformly as we showed in . To compute the propagator given in (5), we still need $`g_0`$, which can be easily obtained from (56):

$$g_0g(0)=\mathrm{\Omega }+\underset{\genfrac{}{}{0pt}{}{m}{m0}}{}G_m(ϵ).$$
(76)
Formulas (66), (74) and (76) can now be used to evaluate $`U(t)`$ for all times. We stress that these formulas depend essentially on the Fourier coefficients $`G_m^{(n)}`$ of $`g`$. These, in turn, depend on $`C_m^{(n)}`$ which are direct linked to $`Q_m`$ and $`Q_m^{(2)}`$, derived from the Fourier coefficients of the interaction function $`f`$. In short, to help us visualise the necessary steps towards the computation of $`U(t)`$, we may draw the following “chain”:

$$\begin{array}{ccccccc}f(t)& & F_m& & Q_m,Q_m^{(2)}& & \alpha _n\\ & & & & & & \\ H_m& & \mathrm{\Omega }& & G_m^{(n)},g_0& & C_m^{(n)}\\ & & & & & & \\ W_m,W_m^{(2)}& & R_m,S_m& & R(t),S(t)& & U(t).\end{array}$$
(77)
### B Interactions with $`F_00`$

When $`F_00`$, one automatically has $`M(𝒬_0)=0`$, except perhaps when $`2F_0=k\omega `$, for some integer $`k`$. These facts were shown in . Since the situation where $`2F_0=k\omega `$, for some integer $`k`$, was nowhere investigated in our previous works , we will ignore this possibility by henceforth assuming that $`2F_0k\omega `$, for all $`k`$. Since $`M(𝒬_0)=0`$, condition I is never satisfied when $`F_00`$. The convergence of our expansions for $`F_00`$ in condition III, however, has not yet been studied. Hence, we will only consider condition II. In this case , the function $`q`$ turns to be

$$q(t)=e^{iF_0t}\underset{m}{}Q_me^{im\omega t}$$
(78)
and the functions $`e_n`$ (see (21)-(22)),

$$e_n(t)=e^{iF_0t}\underset{m}{}E_m^{(n)}e^{im\omega t}.$$
(79)
When $`F_00`$, the coefficients $`E_m^{(n)}`$ in the above expression assume a simple form :

$`E_m^{(1)}`$ $`=`$ $`{\displaystyle \underset{a}{}}{\displaystyle \frac{Q_{m+a}\overline{Q_a^{(2)}}}{a\omega +2F_0}},`$ (80)
$`E_m^{(n)}`$ $`=`$ $`{\displaystyle \underset{p=1}{\overset{n1}{}}}{\displaystyle \underset{a,b}{}}{\displaystyle \frac{Q_{mab}E_a^{(p)}E_b^{(np)}}{(a+b)\omega 2F_0}},\text{for }n2.`$ (81)
Finally, from (78) and (79), we conclude that $`g`$ given in (52) assumes the form

$$g(t)=\underset{m}{}\stackrel{~}{G}_me^{im\omega t},$$
(82)
where

$$\stackrel{~}{G}_m\stackrel{~}{G}_m(ϵ)=\underset{n=1}{\overset{\mathrm{}}{}}\stackrel{~}{G}_m^{(n)}\lambda ^n,$$
(83)
with $`\lambda =ϵ^2`$ and

$$\stackrel{~}{G}_m^{(n)}=\underset{p}{}Q_{mp}E_p^{(n)}.$$
(84)
This gives $`g`$ in terms of its Fourier series.

To compute the propagator $`U(t)`$ accordantly to (5), all we have to do is follow the procedure detailed in Section III A 3 with the replacement of $`G_m^{(n)}\stackrel{~}{G}_m^{(n)}`$ and of $`ϵ\lambda `$.

## IV Numerical Results for the Monochromatic Interactions

We now apply the algorithm discussed in Section III to study the system (3)-(4) under the influence of monochromatic interactions (ac-dc field):

$$f(t)=F_0+\phi \mathrm{cos}(\omega t),$$
(85)
which are relevance for many physical applications. With the conventions introduced in Section III, we have $`J=1`$, $`f_1=f_2=\phi /2`$, $`n_1=n_2=1`$. A simple application of formula (36) gives

$$𝒬_0(t)=\underset{n}{}J_n(\chi _1)e^{i(n+\chi _2)\omega t},$$
(86)
where $`J_n`$ is the Bessel function of first kind and order $`n`$ and where we defined $`\chi _1:=2\phi /\omega `$ and $`\chi _2:=2F_0/\omega `$. Depending on the parameters $`\chi _1`$ and $`\chi _2`$, the function $`f`$ given in (85) satisfies one of conditions I, II or III, except, perhaps, for spurious situations. The detailed analysis of these facts can be found in . Table I summarises the conclusions presented in  and gives a classification of the conditions satisfied by $`f`$ as a function of the parameters $`\chi _1`$ and $`\chi _2`$.

Although not indicated in Table I there is in case (C), on each interval $`(k,k+1),k=1,2,\mathrm{}`$, a special value $`\chi _2^s`$ of $`\chi _2`$ (depending on $`\chi _1`$) for which $`M(𝒬_1)=0`$, and we would be out of condition II . We refrain from studying this rather spurious situation here. A more detailed analysis of this case can be found in .

We next show some graphical results of transition probabilities calculated via the algorithm described in section III for the situations (A), (B) and (C) described in Table I. Let $`|\mathrm{\Phi }_+=\left(\genfrac{}{}{0pt}{}{1}{0}\right)`$ and $`|\mathrm{\Phi }_{}=\left(\genfrac{}{}{0pt}{}{0}{1}\right)`$ be two orthogonal states of a system described by (4) (the eigenstates of the unperturbed Hamiltonian $`H_0`$ of $`H_1`$ in (3)). The probability for the transition from the initial state $`|\mathrm{\Phi }_+`$ to the final state $`|\mathrm{\Phi }_{}`$ at time $`t`$ is given by

$$P(t):=|\mathrm{\Phi }_+\left|U(t)\right|\mathrm{\Phi }_{}|^2=|U_{12}(t)|^2.$$
(87)
We can evaluate $`P(t)`$ numerically using the methods described in Section III to compute $`U(t)`$. To estimate the accuracy of our calculations, we tested the unitarity of the time evolution operator, $`U(t)^{}U(t)=\widehat{1}`$, and considered the quantity

$$N(t):=|U_{11}(t)|^2+|U_{12}(t)|^21,$$
(88)
which should be identically equal to 0 for unitary $`U(t)`$.

Let us first consider case (A) of Table I with $`\omega =1.0`$, $`\chi _1=2/\omega `$ and $`\chi _2=F_0=0`$. Figure 1 shows graphs of $`P(t)`$ and $`N(t)`$ for $`ϵ=0.01`$, $`ϵ=0.10`$ and $`ϵ=0.40`$, plotted from $`t=0`$ to $`t=T_\mathrm{\Omega }=2\pi /\mathrm{\Omega }`$, in units of $`T_\omega =2\pi /\omega `$, the cycle of the external field. The calculations were performed using an expansion for $`g`$ up to $`𝒪(ϵ^{25})`$. We took all the Fourier coefficients (generically called $`_m`$) involved in the computations of $`U(t)`$ within the range $`m=40,\mathrm{},40`$. From the deviations of $`N(t)`$ from 0, we can infer very small errors in the calculations, leading to very accurate values of $`P(t)`$. For $`ϵ=0.01`$ and $`ϵ=0.10`$, we have errors of the order of only $`4.0\times 10^5\%`$. For $`ϵ=0.40`$, the errors jump to $`6.0\times 10^2\%`$, indicating that the parameter $`ϵ`$ is coming close to the radius of convergence of our expansions. Since in case (A) we are under condition I, the secular frequency $`\mathrm{\Omega }`$ is given by (27), hence $`\mathrm{\Omega }=𝒪(ϵ)`$ (we choose $`F_0=\chi _2=0`$). For $`ϵ=0.01`$, from the full expansion of $`\mathrm{\Omega }`$ (see (23)), we get $`T_\mathrm{\Omega }450T_\omega `$. For $`ϵ=0.10`$, $`T_\mathrm{\Omega }45T_\omega `$ and for $`ϵ=0.40`$, $`T_\mathrm{\Omega }25T_\omega `$. We may notice from the graphs of Figure 1, that the transition probability $`P(t)`$ behaves like a Rabi oscillation with a frequency $`\mathrm{\Omega }`$, i.e. $`P(t)\mathrm{sin}^2(\mathrm{\Omega }t)`$. Small oscillations of frequency $`\omega `$, with amplitudes of order $`|ϵ|`$, are superposed with this Rabi oscillation, leading to a quasi-periodic evolution for the system (as ensured by the Floquet theorem).

We now consider case (B) of Table I with $`\omega =10.0`$, $`\chi _1=x_1`$, $`x_1`$ being the first positive zero of $`J_0`$, and $`\chi _2=F_0=0`$. Since we are under condition III and $`F_0=0`$ the secular frequency becomes, according to (30), $`\mathrm{\Omega }=𝒪(ϵ^3)`$, a fact first pointed in . This weak dependence on $`ϵ`$ indicates long transition times for the probability amplitude $`P(t)`$. This phenomenon is known as the dynamical localisation effect  (see  for a general criterion and for a more complete list of references). To compute the solution in this case, $`g`$ was expanded up to $`ϵ^6`$. This contrasts with the expansion for $`g`$ in the previous situation, where $`g`$ has been computed up to $`𝒪(ϵ^{25})`$. The reason for this lies in the fact that under condition III the expressions for the constants $`\alpha _n`$ (see (11)-(13)) are somewhat more intricate than those of condition I. Fortunately, however, the expressions for the constants $`\alpha _1`$, $`\alpha _2`$ and $`\alpha _3`$ are not so complicated and this suffices to give an expansion of $`g`$ free of secular terms up to the $`ϵ^6`$ term (see  for details). If we limit ourself to study situations of small $`ϵ`$, i.e. well inside of the radius of convergence of our expansions, excellent results can be achieved. To improve the accuracy of our results, it is better to use a larger value of $`\omega `$.

Figure 2 shows graphs of the numerical values of $`P(t)`$ and $`N(t)`$ for $`ϵ=0.01`$, $`ϵ=0.10`$ and $`ϵ=0.20`$, plotted from $`t=0`$ to $`t=T_\mathrm{\Omega }=2\pi /\mathrm{\Omega }`$, in units of $`T_\omega =2\pi /\omega `$. The errors in Figure 2 (measured out of the deviations of $`N(t)`$ from zero) are greater than those presented in Figure 1 due to the lower order $`ϵ`$-expansion of $`g`$ in the former case. However, for small $`ϵ`$, the errors are still small, being of the order of only $`3.0\times 10^3\%`$ for $`ϵ=0.01`$ and of the order of $`3.0\times 10^1\%`$ for $`ϵ=0.10`$. For larger $`ϵ`$ ($`=0.20`$), the errors reach the value of about $`1.0\%`$, indicating that a higher order $`ϵ`$-expansion is needed to improve accuracy. We observe that the plots of $`P(t)`$ in Figure 2 have the same qualitative aspect of those in Figure 1. There is a predominant Rabi oscillation (with frequency $`\mathrm{\Omega }`$) aspect superposed by minor oscillations of frequency $`\omega `$ whose amplitudes are possible bounded by $`|ϵ|^3`$ terms. This fact, however, has yet not been proven directly from our expansions, so it has to be faced more as a qualitative analysis rather than a rigorous quantitative one.

The main distinction between the graphs of Figure 1 and those of Figure 2 is, undoubtedly, the Rabi oscillation period $`T_\mathrm{\Omega }`$. In the later, the long time needed for the system to transit from the initial state $`|\mathrm{\Phi }_+`$ to the final state $`|\mathrm{\Phi }_{}`$, compared to the basic cycle $`T_\omega `$ of the external perturbation, is the (approximate) dynamical localisation effect. The secular period $`T_\mathrm{\Omega }`$ obtained for the situations studied in Figure 2 were $`T_\mathrm{\Omega }1.6\times 10^9T_\omega `$ for $`ϵ=0.01`$, $`T_\mathrm{\Omega }1.6\times 10^6T_\omega `$ for $`ϵ=0.10`$ and $`T_\mathrm{\Omega }2.2\times 10^5T_\omega `$ for $`ϵ=0.20`$. These are somewhat much larger values of $`T_\mathrm{\Omega }`$ than those of Figure 1.

As we have mentioned, the graphs of Figure 2 were computed with $`\chi _1=x_1`$, where $`x_1`$ is the first positive zero of $`J_0`$. It is interesting to test our solution considering other possible zeros of $`J_0`$. In Figure 3 we show plots of $`P(t)`$ and $`N(t)`$ calculated with $`\chi _2=x_2`$, the second positive zero of $`J_0`$ (we also took $`\omega =10.0`$ and, of course, $`\chi _2=0`$). The three situations presented in Figure 3 correspond to $`ϵ=0.10`$, $`ϵ=0.20`$ and $`ϵ=0.30`$. We may note that the qualitative behaviour of $`P(t)`$ presented in Figure 3 is the same of Figure 2. In particular, the effect of dynamical localisation is preserved since we still have $`\mathrm{\Omega }=𝒪(ϵ^3)`$ when $`\chi _1=x_2`$. We may note from the deviation of unitarity $`N(t)`$ shown in Figure 3 that the errors, compared with the ones in Figure 2 for the same values of $`ϵ`$, are smaller by a factor $`1/10`$. We may understand this fact as follows: since $`x_2>x_1`$, the strength $`\phi `$ of the interaction presented in Figure 3 is greater than in Figure 2, hence, for the same values of $`ϵ`$, the effective perturbation $`ϵ/\phi `$ is smaller in the former case, leading to a more precise perturbative computation.

We now investigate case (C) of Table I. We consider $`\omega =1.0`$, $`\chi _1=1.0`$ and $`\chi _2=0.3`$. This implies condition II and we have to follow the prescriptions of Section III B. Figure 4 shows the results obtained for $`P(t)`$ and $`N(t)`$ for three values of $`ϵ`$: 0.05, 0.10 and 0.20. The various plots were calculated using an expansion of $`g`$ up to $`ϵ^{20}`$. It was not necessary, thus, to use large $`\omega `$ to ensure convergence of the expansions, as we did in the situations presented in Figure 2 (condition III). As usual, $`P(t)`$ and $`N(t)`$ were computed from $`t=0`$ to $`t=T_\mathrm{\Omega }`$ in units of $`T_\omega `$. The qualitative behaviour of the transition probability is significantly different from the previous results (Figures 1 and 2). This is a consequence of the non-vanishing constant field $`F_0`$ presented in the interaction $`f`$. Indeed, the secular frequency $`\mathrm{\Omega }`$ is now given by (29) and, hence, is of order of $`F_0`$. For $`ϵ=0.05`$, we obtained $`T_\mathrm{\Omega }6.3T_\omega `$, for $`ϵ=0.15`$, $`T_\mathrm{\Omega }6.0T_\omega `$ and for $`ϵ=0.20`$, $`T_\mathrm{\Omega }5.0T_\omega `$. Since we choose $`F_0`$ as the same order of $`\omega `$, there is a strong competition between the Rabi oscillation (governed by $`\mathrm{\Omega }`$) and the external field oscillations. This leads to the patterns shown in Figure 4, which do not behave purely like $`\mathrm{sin}^2(\mathrm{\Omega }t)`$. In particular, $`P(t)<1`$ for all times, leading to the conclusion that the state $`|\mathrm{\Phi }_+`$ never transits completely to $`|\mathrm{\Phi }_{}`$. One sees, moreover, that the transition amplitude $`P(t)`$ tends to zero as $`ϵ0`$ (c.p. (91), below), much in contrast to the cases pictured in Figs. 1, 2 and 3. The leading $`ϵ`$-dependence of $`P(t)`$ can be algebraically determined, in principle, but this was not yet performed due to the complexity of our expansions.

It is interesting to note from the graphs of Figure 4 that the transition probability gets closer to 1 as $`ϵ`$ increases. It would seem that for $`ϵ`$ large enough, we could have $`P(t)>1`$. However, our expansions would not converge in this case, since $`ϵ`$ would be greater than the radius of convergence. To illustrate this situation, let us consider the trivial case where $`\chi _1=0`$ and $`\chi _2=2F_0/\omega `$ is not an integer. Since now $`f(t)=F_0=\omega \chi _2/2`$, the Schrödinger equation (4) becomes time-independent, with a simple Hamiltonian given by

$$H_2=ϵ\sigma _1+F_0\sigma _3=\left(\begin{array}{cc}F_0& ϵ\\ ϵ& F_0\end{array}\right).$$
(89)
The propagator $`U(t)`$ can be computed by elementary methods (f.i., by diagonalising the Hamiltonian), leading to

$$U(t)=\mathrm{cos}(\omega _0t)\widehat{1}\frac{\mathrm{sin}(\omega _0t)}{\omega _0}\left(iF_0\sigma _3+ϵ\sigma _1\right),$$
(90)
where $`\omega _0:=\sqrt{F_0^2+ϵ^2}`$. Hence,

$$P(t)=\frac{ϵ^2}{F_0^2+ϵ^2}\mathrm{sin}^2(\omega _0t).$$
(91)
Note that $`P(t)<1`$ for all $`ϵ`$ and all $`F_0`$. For $`f(t)=F_0`$, the generalised Riccati equation (8) admits a particular solution given by the constant $`g_0=F_0+\sqrt{F_0^2+ϵ^2}`$. It was shown in  that the $`ϵ`$-expansion (20) for $`f(t)=F_0`$ coincides, as expected, with the Taylor expansion (centred at $`ϵ=0`$) of $`g_0`$. Thus, our method has a clearly restricted region of convergence defined by $`|ϵ|<|F_0|`$.

Let us see what these last considerations mean numerically. First we set $`F_0=0.4`$. Figure 5 shows plots of $`P(t)`$ and $`N(t)`$ for three critical values of $`ϵ`$: $`ϵ=0.30<F_0`$, $`ϵ=0.40=F_0`$ and $`ϵ=0.43>F_0`$, where our expansions are not supposed to converge. We also show in Figure 5 (dashed lines) plots of $`P(t)`$ calculated according to (91). Looking at $`N(t)`$ and at the deviation of the perturbatively computed transition probability from the one calculated via (91), we conclude that the region of convergence of our expansions is restricted to $`|ϵ|<|F_0|`$, as expected. It is important to note that when $`ϵ<F_0`$, our computation of the propagator matches exactly with the non-perturbative solution (90). This can be seen from comparing $`P(t)`$ calculated perturbatively and non-perturbatively via (91), as shown in Figure 5.

## V Final Remarks

We stress that, at least for the case of monochromatic interactions examined in Section IV, the errors obtained are very small and bounded as time increases (see the behaviour of $`N(t)`$ in figures of Section IV). This is due to the absence of secular terms in our perturbative expansions and its uniform convergence in time. As a consequence, one can study the long-time behaviour of the quantum system (3)-(4) in a controlled way.

Another important feature of our method is that it can be numerically implemented with relatively simple computer codes. Indeed, the algorithm described in Section III to calculate the unitary propagator $`U(t)`$ consists, essentially, in simple computations involving the Fourier coefficients of the functions $`q`$ and $`q^2`$, which are known in a closed form (see (33) and (36)). This is an important advantage against other perturbative approaches, based f.i. on the Dyson expansion (2) which, in general, cannot be evaluated in a simple manner (not to mention the fact that such expansion is not uniform convergent in time, as we have stressed in the introduction of this paper). To sum up, our method is not only mathematically rigorous, but also very useful for practical purposes, where it can be applied with great generality, leading to very accurate results.

###### Acknowledgements.

We are grateful to A. Sacchetti for discussions. J. C. A. Barata was partially supported by CNPq. D. A. Cortez was supported by FAPESP.


text:
# The ARAUCARIA project. Discovery of Cepheid Variables in NGC 300 from a Wide-Field Imaging Survey 1footnote 11footnote 1Based on observations obtained with the 2.2 m ESO/MPI telescope at the European Southern Observatory

## 1 Introduction

One of the most fascinating and challenging chapters of modern astrophysics has been the effort to establish the extragalactic distance scale with an ever increasing accuracy. Key to this effort has been the use of the Hubble Space Telescope over the last years which has allowed to detect and measure the brightness of Cepheid variables in some 20 nearby spiral galaxies and determine their distances from the Cepheid period-luminosity relation. In fact, the HST Key Project team on the extragalactic distance scale has just recently presented their final results for the Cepheid distances to their selected galaxies (Freedman et al. 2001), from which secondary methods of distance measurement have yielded an improved determination of the Hubble constant. However, in the HST Key Project procedure to determine the Hubble constant there are still a number of systematic uncertainties which prevent, at the present time, a truly accurate determination of this constant which is so important to cosmology. The most nagging problem continues to be the distance to the LMC, whose Cepheids are used as a fiducial to determine the relative distances to the other galaxies. The distance to the LMC is currently the largest single contributor to the systematic uncertainty on $`\mathrm{H}_0`$ (Mould et al. 2000). Another problem, affecting the calibrating Cepheid samples in different galaxies, is the dependence of Cepheid absolute magnitudes on metallicity. Currently there are very controversial claims about the importance of this effect for the determination of Cepheid distances to galaxies, and a truly accurate empirical determination of the effect is urgently needed to further reduce the uncertainty on the Hubble constant.

In spite of the very important progress achieved over the last years, not only by HST but by many other teams of researchers using ground-based telescopes, it seems clear that the remaining uncertainty in the distance scale can only be significantly reduced if we succeed to make true progress in determining the distance to the LMC, and a number of other nearby galaxies. To achieve this, we must learn more about the systematic uncertainties in the principal methods currently used to find these distances. The most promising stellar methods, apart from Cepheid variables, use red clump giants, the tip of the red giant branch, eclipsing binaries, and planetary nebulae as distance indicators. Another, relatively new tool to determine extragalactic distances from individual stars is the wind momentum-luminosity relationship (WLR) for blue supergiants (Kudritzki et al. 1999; Bresolin et al. 2001). Only when the true capabilities and systematics of these various methods of distance determination are fully understood, the techniques properly calibrated, and agreement on the distances to a number of nearby galaxies with different environmental properties from all these methods is achieved, can we trust that the local calibration of the distance scale has been finally settled. A working group at the Universidad de Concepción (including several of the authors of this paper) has recently been established to work on the improvement of stellar distance indicators. We call our project the ARAUCARIA project (making reference to a famous tree at our southern location), and we hope to make a significant contribution to the improved calibration of the local extragalactic distance scale, over the next years.

Our target to start with is the Sculptor Group spiral galaxy NGC 300. At a distance of about 2.0 Mpc (Freedman et al. 2001), this galaxy is sufficiently close to resolve its stellar populations. Since NGC 300 shows clear signs of recent, massive star formation, it seemed an excellent target to discover a significant number of Cepheid variables. Indeed, a relatively small number of Cepheids in NGC 300 had already been previously discovered by Graham (1984) with photographic techniques at the CTIO 4 m telescope, and a subsample of these variables was later observed with a ground-based CCD camera by Freedman et al. (1992). However, taking advantage of new wide-field imaging capabilities, we hoped to find many more Cepheids and determine their properties with a higher accuracy than what had been achieved in former work, a hope which is borne out by the results presented in this paper. In a parallel project, we are conducting a survey for blue supergiant stars in the disc of NGC 300 with the ESO VLT telescope; first results have recently been presented in Bresolin et al. (2001). The blue supergiants will be used to establish the (stellar) metallicity gradient in the disc of NGC 300 to perform a new and accurate empirical calibration of the effect of metallicity on Cepheid luminosities, one of the principal goals of our current project. They will also be used to improve the WLR calibration for blue supergiants. In addition, we have also secured narrow-band images of the galaxy to conduct a complete survey for its planetary nebulae population and improve the calibration of the PNLF method of distance determination.

The present paper is organized as follows. In section 2, we describe the observations from which the current survey for Cepheids was conducted. In section 3 we describe the reduction techniques and calibration of the images. In section 4, we present the catalog of newly discovered Cepheid variables in NGC 300 (variable stars of different types will be presented in a forthcoming paper). In section 5, we compare our new data to the ones for those Cepheids for which previous photometry exists from the work of Graham, and of Freedman et al. (1992). A discussion of the results and conclusions are presented in section 6.

## 2 Observations

All observations presented in this paper have been collected with the ESO/MPG 2.2 m telescope at the La Silla observatory in Chile. The telescope was equipped with the ESO mosaic Wide Field Camera (WFI) consisting of eight 2048 x 4098 pixel arrays. The CCDs were separated by gaps of 23.8 and 14.3 arcsec in right ascension and declination directions, respectively. The total field of view was about 34 x 33 arcmin with a scale of 0.238 arcsec / pixel.

Our observations started on July 31, 1999 and lasted until January 8, 2001. During this period, we secured 150 B, 138 V, 49 R and 145 I-band images centered on NGC 300, on 29 different nights allocated to our programme. Most of these nights were photometric, and the seeing ranged from 0.7 to 1.5 arcsec during the vast majority of the exposures taken. In order to obtain photometry for all stars from the whole observed field including those located in the gaps between CCDs, the observations were carried out in a ”dithering” mode. During a given night 5 consecutive 360 s exposures, slightly shifted in right ascension and/or declination with respect to each other, were made in a given filter. This technique also allowed us to obtain a better accuracy of the photometric measurements by evaluating a mean from all observations taken during a given night through a given filter, and help to recognize bad measurements caused by weather, blending, or presence of bad pixels. On about half of the nights we obtained full observing sequences through BVRI filters, in the other nights we could only observe BVR or BVI sequences. On a few nights, only BV sequences could be obtained. Effectively, 23-30 epochs, depending on the brightness of a given star and filter, were obtained. The journal of B and V observations which were used to conduct the search for Cepheids is presented in Table 1.

## 3 Reductions and Calibrations

Preliminary reductions (i.e. debiasing and flatfielding) were carried out in the standard way with the IRAF222IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the NSF. package. In order to avoid possible variations of the PSF across the image each of the 8 chips was divided into 2 slightly overlapping subframes. Then profile photometry for all stars detected on all subframes was performed with the $`DAOPHOT`$ and $`ALLSTAR`$ programs. The PSF model was derived iteratively. First, 10 relatively bright and isolated stars were selected and the first approximation of the PSF model was calculated. In the next step we subtracted all neighbor stars with the $`ALLSTAR`$ program and derived the PSF model again. The PSF model obtained after three such loops was finally adopted. The photometric data from each two subframes were finally tied together using stars located in the overlapping regions.

In order to prepare an unique system of identification of stars in the observed field, the best (i.e. obtained under the best seeing conditions) set of five V images was selected and stacked together. The resulting ”template” included all stars located in the gaps between chips and covered a region of about 34 x 33 arcmin (see Fig. 1). The reference list with coordinates of objects from the template image was prepared with the $`DAOPHOT`$ program. After rejecting evidently spurious detections due to saturated stars or wild pixels, the total number of objects we detected was about 32000. The pixel coordinates of all stars from a given field were transformed to the template coordinate system, and a unique number corresponding to the numbering system from the reference list was assigned to each object.

A detailed description of the adopted photometric calibration procedure, completeness tests, a comparison with previous results together with a presentation and discussion of CMDs will be presented in a forthcoming paper. Briefly, to calibrate the photometric magnitudes of all stars on the wide field images an extensive sequence of secondary standard stars, distributed over almost the whole observed area and spanning a broad range in magnitudes and colors, was set up (Pietrzyński et al. 2001). The WFI photometry was transformed to the standard system using about 100 stars from this list. Derived residuals usually did not exceed 0.04 mag and did not show any dependence on color or brightness. The accuracy of the zero point of our WFI photometry is about 0.03 mag in the standard V and B bands. We found good agreement (the mean differences were $`\mathrm{\Delta }V=0.008\pm 0.028`$ mag, $`\mathrm{\Delta }(BV)=0.018\pm 0.031`$ mag) between our photometry and that obtained by Walker (1995), for 22 common stars in both samples.

Taking into account the periods of the potential Cepheids, all consecutive observations obtained in a given filter, during a given night (amounting to typically 30 minutes) still define practically the same phase in the light curve, even for the shortest-period Cepheids we detected (see section 4). For that reason we decided to improve the accuracy of our Cepheid measurements, for a given night and in a given filter, by adopting the mean magnitude from the individual dithered observations.

Equatorial coordinates were calculated with the algorithm developed and used by the OGLE team (Udalski et al. 1998). Shortly, a FITS file slightly larger than our template was extracted from the Digital Sky Survey (DSS) images and all stars having more than 200 counts above the sky level were detected on it. The pixel (x,y) coordinates of the detected stars were transformed to equatorial coordinates, and then to (x’,y’) pixel coordinates on the plane tangential to the celestial sphere at the center of our field. The resulting two sets of pixel coordinates were tied together using third order polynomials. The internal accuracy of our transformation was about 0.3 arcsec.

## 4 The Cepheid Catalog

As a first step, all stars were subjected to a period search using the AoV algorithm (Schwarzenberg-Czerny 1989). We searched for periods between 0.2 and 90 days, the latter value corresponding to about a half of the time span of our observations. For the detected Cepheids, the accuracy of the derived periods is about 0.5-1 $`\times 10^2\times \mathrm{P}`$, or roughly 1 percent of the period lengths. However, for some of the Cepheids, the ones with previous observations from Freedman et al. (1992) and/or Graham (1984), we were able to determine much more accurate periods by combining the different sets of photometric measurements (see chapter 5 for more details).

In order to distinguish Cepheid variables from other types of variable stars we used the following selection criteria:

$$\begin{array}{c}\text{1. Stars must exhibit the typical asymmetrical Cepheid-like shape of the light curves}\hfill \\ \text{2. }\mathrm{A}_\mathrm{B}>0.4\text{ mag}\hfill \\ \text{3. }\mathrm{A}_\mathrm{B}>\mathrm{A}_\mathrm{V}\hfill \\ \text{4. 0.4 }<\text{ B-V }<\text{ 1.5}\hfill \end{array}$$

Here, $`\mathrm{A}_\mathrm{B}`$ and $`\mathrm{A}_\mathrm{V}`$ stand for the amplitudes in the B and V bands, respectively. Regarding the color selection criterion, we note that the foreground reddening in the direction to NGC 300 is only E(B-V)=0.02 (Burstein & Heiles 1984) so that the observed colors correspond closely to the intrinsic colors of the stars, supposing that intrinsic absorption within NGC 300 is small.

To derive the mean magnitudes for those variables passing the selection test for a Cepheid, the light curves were approximated by Fourier series of an order ranging from 2 to 5, depending on the accuracy of the photometry and/or the phase coverage. Fig. 2 shows such fits for two exemplary B light curves: one of a bright, well observed Cepheid, and one of a much fainter Cepheid having a noisier light curve and less points because the star could only be detected and measured on nights of good seeing. However, even in such cases mean magnitudes could be determined rather accurately. They were obtained by fitting the light curves, converted to intensity units, and transforming the mean intensities back to a magnitude scale. The statistical accuracy of our intensity mean B and V magnitudes was typically better than 0.01 mag (i.e. much smaller than the uncertainty of the zero point of our photometry).

Alltogether, 117 stars satisfied our selection criteria and entered our catalog of Cepheids presented in this paper. Table 2 contains their description. The first column is the star identification number. In the next columns the equatorial coordinates, derived periods, epochs of zero phase corresponding to maximum light in V, and intensity mean brightnesses in the B and V bands are given. The last column contains remarks on the objects. V2, V3, etc correspond to the numbering system introduced by Graham (1984) for the 34 variable star candidates in NGC 300 discovered by him. The digits of the periods given reflect their uncertainties. The phased B and V light curves for all Cepheids are presented in Fig. 3. Note that each observational point on the light curves is the mean from all observations (up to 5) obtained through a given filter during a particular night.

In addition to the 117 Cepheids in Table 2, we present another 12 Cepheid candidates in Table 3. These objects do not satisfy all criteria described above, but the shapes of their (noisy) light curves suggest that they may be Cepheids, too. The phased B and V light curves of these objects are displayed in Fig. 4. For all Cepheids, the individual B and V observations are given in Table 4.

Finding charts for all Cepheid variables, prepared based on V band image are presented in Figure 5 (available only in the electronic edition of the Journal). They can be also be obtained from the authors upon request.

## 5 Previously Known Cepheids

As mentioned before, a previous survey for variable stars in NGC 300 has already been conducted by Graham (1984). His observations, based on photographic plates taken at the prime focus of the CTIO 4 m telescope between November 1974 and October 1981, yielded a list of 34 variables or variable star candidates which included 18 Cepheids with derived periods, and another 5 objects suspected to be Cepheids, but without enough data for a period determination. Later Freedman et al. (1992) obtained CCD BVRI photometry at the same telescope for 16 of the Cepheids discovered by Graham . The baseline of their new observations was about 3 years, and as expected the quality of the CCD data was much higher than the previous photographic observations. However, in many cases the phase coverage of the resulting light curves was poor, making the determination of the mean brightnesses rather uncertain. Combining the photographic and their own CCD data, Freedman et al. (1992) had already obtained refined periods for some of the Cepheids.

From our new WFI data we rediscovered all 18 Cepheids reported by Graham (1984) (see Table 2). We also confirm 3 of his Cepheid candidates to be truly Cepheids (V7, V20, and V31). All the Graham variables were among the brightest objects in our database and have the most accurate photometry and the best coverage of their light curves. We also found that the object designated by Graham as V4 and classified by him as an eclipsing binary is another Cepheid with P = 15.37 d. Most probably the spurious classification was due to poor resolution of the photographic data. We did not confirm the variability of the objects marked as V14 and V19. Regarding the periods, our values generally agree, within the quoted uncertainties, with those obtained in the past. However, in the case of 5 variables (V4, V18, V21, V24, V25) we found a significantly different value for the period. The discrepancy in the period determination for V24 may be the result of a genuine period change for this variable during the last 26 years. Unfortunately the previous observations are too scarce and uncertain to confirm this beyond doubt. For V21 we obtained a period of about 24.0 days instead of 9.667 days derived previously by Freedman et al. However, the Freedman et al. observations for V21 do fit much better with our period than with that derived by themselves. The same happens in the case of star V25. For both of these Cepheids the previous data look much better when phased with our periods than with those obtained before, and we were able to refine the periods (see below) combining our new with the previous data. The photographic data obtained for stars V4 and V18 are not accurate enough to be useful to improve the periods from the current new data.

Including our new observations, the time baseline since the first epoch observations is now 26 years, corresponding to about 76-1000 pulsation cycles for the Cepheids, which makes further improvement of the periods of these stars possible. Before doing this we need to make sure that all the data sets share the same photometric zero point. Walker (1988) had already shown and corrected a zero point and scale error in Graham’s photometric data. Using this correction, we tied the Graham data to Freedman’s CCD system using the relations derived by Freedman et al. (1992). Finally, we need to check if the Freedman et al. (1992) photometry is not significantly different from ours. We did this by comparing the intensity mean magnitudes derived from our measurements to those given by Freedman et al., for the common stars. Fig. 6 presents this comparison for 8 Cepheids with reasonably good photometry from Freedman et al. It is seen that both data sets yield mean magnitudes which are in very good agreement. We also checked the agreement between individual observations from these two data sets. An example is given in Fig. 7 for the variable V3. It is seen that a systematic zero point offset, if any, is very small, and that the individual observations generally agree very well. It is also worth noting that apart from their much larger scatter, the photographic data fit also quite well our observations.

After making sure that there are no significant shifts in photometric zero points of the three available data sets, we phased all the data for common objects with the best available period (usually that from Freedman et al., (1992)), looking for phase shifts to our present new data. In 15 cases we found clear shifts in the phase indicating the need to adjust the period. The periods of these stars given in Table 1 are already the improved periods resulting from a combination of the available data sets. For 7 remaining Cepheids, the long gaps between the individual data sets, the accuracy of the available periods, the large scatter in the photographic data, and a very limited number of previous CCD observations (if available) did not allow us to refine the periods that we derived from our data alone.

## 6 Summary and Discussion

Results are presented of a wide-field multi-color photometric monitoring of NGC 300. Based on observations obtained on 29 nights during a 5.3 month period, 117 Cepheids and 12 Cepheid candidates were discovered from images in B and V bands. Their periods range between 115 and 5.42 days. For each object we provide equatorial coordinates, period, intensity mean B and V magnitudes and phased light curves. The accuracy of the derived periods is generally about $`0.5110^2`$ P. We have re-discovered all previously known Cepheid variables, and for many of these the periods could be improved by combining the available data sets. The photometric zero points in both B and V are secure to within 0.03 mag.

The spatial distribution of the detected Cepheids is shown in Fig. 1. It can be appreciated that we were able to detect Cepheids even in the very dense regions near the center of NGC 300. Since the Cepheid distribution traces the spiral arms of the galaxy, most of our Cepheids lie in dense fields and are most probably blended, to some degree, with unresolved stars. This blending problem will be discussed in more detail in a forthcoming paper in which we will derive the distance of NGC 300 based on Cepheid PL relations.

In Fig. 8, we show the location of the detected Cepheids on the observed V, B-V color-magnitude diagram constructed from our data. The Cepheids delineate the typical Cepheid instability strip. Fig. 7 also indicates that the completeness of our Cepheid search drops suddenly for objects fainter than about 22.5 mag in V. For magnitudes brighter than this, we expect the survey to be near-complete. The corresponding period cutoff is close to 10 days, in the sense that our survey should have discovered a very large fraction of the Cepheids in NGC 300 with periods larger than 10 days.

Our Cepheid catalog contains 22 entries common to the list of variable stars prepared by Graham (1984). Sixteen of the Cepheids had also accurate (albeit sometimes sparse) previous CCD photometry (Freedman et al. 1992). We found very good agreement between the zero points in B and V of our data and those presented by Freedman et al. Using combined data sets, the periods of common objects were examined and improved whenever possible.

Fig. 9 presents the period distribution for all the 117 Cepheids in our catalog. There are three evident peaks, corresponding to periods of about 8, 14 and 24 days. Their existence and locations do not depend on the size of the binning we are using. A precise interpretation of the physical meaning of these peaks is difficult (and beyond the scope of this paper) but one possibility is that they reflect the fact that the Cepheids we are observing in NGC 300 have been formed during different bursts of star formation, in different places in the galaxy.

Finally, we show the period-luminosity relations in B and V derived from our data in Fig. 10. For periods below about 10 days, the data are clearly affected by a Malmquist bias-only the brighter Cepheids, at any given period below 10 days, have been detected in our survey. A detailed analysis of the PL relations to determine the distance of NGC 300 will be the subject of a forthcoming paper, which will include the presentation and analysis of photometric data obtained at longer wavelengths.

We are very grateful to the European Southern Observatory for providing large amounts of observing time with the Wide Field Camera on the 2.2 m ESO/MPI telescope to this survey. It is a great pleasure to thank the 2p2 team for expert support. Among the astronomers who have been especially helpful with the observations we would like to thank H. Boehnhardt, P. Leisy, J. Brewer, E. Pompei, D. Bersier, H. Jones, P. Francois and T. Augusteijn. Pierre Leisy and L. Rizzi were very helpful with the data reductions. Furthermore, we would like to thank Drs Andrzej Udalski, Michał Szymański, and Grzegorz Pojmański for making their computer programs available and valuable discussion. WG would like to acknowledge financial support from Fondecyt grant 8000002, and express his thanks to the Centrum fuer internationale Migration und Entwicklung (CIM) in Germany for providing funds for a Sun Ultra 60 workstation on which most of the present data analysis was carried out. We also acknowledge use of the Digitized Sky Survey which was produced at the Space Telescope Science Institute based on photographic data obtained with the UK Schmidt Telescope, operated by the Royal Observatory Edinburgh.


text:
# Time–optimal Hamiltonian simulation and gate synthesis using homogeneous local unitaries

## I Introduction

### A Simulation of non–local Hamiltonians and synthesis of non–local gates

The simulation of Hamiltonians has attracted considerable attention in several areas of quantum physics, including those of nuclear magnetic resonance  (NMR), quantum control  and quantum information . Let us consider a quantum system that naturally evolves according to some Hamiltonian $`H`$ for a time $`t`$. By interspersing its evolution $`\mathrm{exp}(iHt)`$ with a sequence of fast unitary operations $`W_k`$ , another evolution

$$e^{iH^{}t^{}}=e^{iHt_n}W_n\mathrm{}e^{iHt_2}W_2e^{iHt_1}W_1$$
(1)
is obtained which, when $`|Ht_k|<<1`$, approximately corresponds to having the Hamiltonian

$$H^{}=\frac{c}{t}\underset{k}{}t_kV_kHV_k^{},$$
(2)
$`V_k^{}W_k\mathrm{}W_1`$ \[we assume $`V_n=I`$\], acting for a time $`t^{}t/c`$, where $`c0`$ is the time overhead of the simulation. This means that we can effectively generate, given an available Hamiltonian $`H`$ and some control operations $`W_k`$, a new Hamiltonian $`H^{}`$, which may otherwise not be accessible experimentally. Hamiltonian simulation is therefore a powerful technique to control the evolution of a quantum system. As such, it has a broad range of potential applications, and may eventually become commonplace in experimental implementations of quantum information and quantum computation.

The synthesis of a gate $`U\mathrm{exp}(iH^{}T)`$ by using $`H`$ and control operations $`W_k`$ is a closely related issue, for $`U`$ can be achieved by simulating $`H^{}`$ for a time $`T`$. We notice, however, that inequivalent sequences of control operations $`W_k`$ may also lead to the same gate $`U`$ —the only requirement is that at the end of the protocol the gate has been performed— and thus simulating $`H^{}`$ may not be the most convenient way to proceed.

A series of recent contributions  has recently addressed, in the context of multipartite quantum systems, the problem of simulating one non–local Hamiltonian (or interaction) by using another non–local Hamiltonian and local unitary transformations (LU) as control operations. In particular, special emphasis has been placed on assessing the universality of this model , as well as on determining time-optimal simulation protocols and developing a quantitative theory of interactions based on their simulation capabilities . Similarly, time-optimal protocols for non–local gate synthesis have been thoroughly analyzed in the context of quantum control , where the problem has been reduced, in the case of a two-qubit system, to a problem of non–local Hamiltonian simulation. In  the corresponding optimal non–local Hamiltonian, and thereby the minimal time for two–qubit gate synthesis, have been presented.

### B Hamiltonian simulation in solid state physics: <br>homogeneous local manipulation

In most of previous work on non-local Hamiltonian simulation, arbitrary LU are assumed to be available as control operations. This model is motivated by the fact that in a number of experimental schemes for quantum information processing (using e.g. ion traps, neutral atoms or photons ) independent LU operations can be performed on each of the interacting systems supporting the qubits. In solid state implementations, however, the relevant systems are frequently too close to each other to be addressed individually. A good example of this is a set of spins in a NMR molecule or in a lattice of quantum dots, where an extraordinarily focused (and thus unfeasible) magnetic field would be required to independently address each spin. Therefore, and due to experimental limitations, the current model for interaction simulation does not apply to a number of solid state schemes.

Most relatedly, an alternative, solid-state oriented model for quantum computation has recently been put forward . It overcomes the above lack of controllability by encoding logical qubits in convenient subspaces of few physical qubits, in a way that universal computation can be accomplished by just tunning the intensity of a specific non–local Hamiltonian —such as the isotropic exchange interaction, naturally available in many systems. Notice, however, that in solid state not all forms of local manipulation are equally inaccessible. In particular, a uniform magnetic field can be easily applied on all spins, inducing the same local evolution on each of them and thereby producing what we shall call an homogeneous local unitary transformation (HLU). In this sense, HLU appear as a rather restricted, but definitely feasible, set of control operations in solid state physics, that may be worth considering.

### C Results

In this paper we address the problems of non–local Hamiltonian simulation and gate synthesis using a given non–local Hamiltonian and fast HLU. Our results, mostly concerned with two-body interactions in multi-qubit systems, imply that an interaction, when only supplemented with HLU, does not constitute a universal resource for simulating other interactions or synthesizing arbitrary gates, in sharp contrast with the inhomogeneous LU model . Nevertheless, the HLU model still allows for a rich variety of Hamiltonian simulations and gate synthesis. In particular, HLU can be used to eliminate unwished anisotropic terms (as coming e.g. from a spin-orbit interaction) in a system of spins with imperfect isotropic exchange interaction. This task is of practical relevance, for instance, in the context of encoded universality, and has been considered in . Therefore the HLU model for non-local Hamiltonian simulation provides us with a complementary or alternative technique to that discussed in .

We shall present the following results.

($`i`$) Necessary and sufficient conditions for two-qubit Hamiltonian simulation under HLU, both for symmetric and antisymmetric interactions. These results also apply to some multiple-qubit scenarios with two-body interactions.

($`ii`$) Time-optimal synthesis of two-qubit gates using a symmetric interaction.

($`iii`$) Recovery of universality in gate synthesis using HLU and just one initial and one final non-homogeneous LU per two-qubit gate.

A remarkable feature of the HLU model of simulation will be that it is ruled by the usual majorization relation , instead of the special majorization relation that rules the non–homogeneous LU model . As we shall explain, this gives a neat picture of which two-qubit Hamiltonian simulations are possible in the symmetric case: $`H`$ can simulate $`H^{}`$ if and only if this implies a gain of isotropy in the interaction.

The organization of this paper goes as follows. In sect. II we address the problem of two-qubit Hamiltonian simulation and present the optimal solution in the case of symmetric and antisymmetric interactions. Using these results, the time-optimal synthesis of symmetric gates using a symmetric interaction is analized in sect. III. In sect. IV we discuss the synthesis of symmetric gates using an arbitrary interaction. In sect. V we also show how to achieve any two-qubit gate by allowing for two non–homogeneous interventions. The paper finishes with some conclusions.

## II Optimal Hamiltonian simulation using HLU

Consider a collection of $`N`$ interacting systems, say $`N`$ qubits for concreteness, that evolve according to a Hamiltonian $`H_{12\mathrm{}N}`$. If only fast HLU can be performed on the qubits in order to control their evolution, i.e. $`W_k=w_k^N`$ in Eq. (1), then, according to Eq. (2), any simulated Hamiltonian $`H_{12\mathrm{}N}^{}`$ must be of the form

$$H_{12\mathrm{}N}^{}=c\underset{k}{}p_k(v_k^N)H_{12\mathrm{}N}(v_k^N)^{},$$
(3)
where $`p_kt_k/t`$ sum up to one. It already follows from this expression that, in general, it is not possible to simulate an arbitrary Hamiltonian $`H_{12\mathrm{}N}^{}`$ by $`H_{12\mathrm{}N}`$ and HLU alone. As an example, suppose that $`H_{12\mathrm{}N}`$ is symmetric under exchange of qubits. Then $`H_{12\mathrm{}N}^{}`$ necessarily has the same symmetry , so that we cannot simulate any inhomogeneous evolution. Thus, the HLU model for simulation is rather restricted. Nevertheless, it is still a useful and —most important— available tool to control the evolution of the qubits.

We start our analysis by considering $`N=2`$ qubits, for which we present most of our results. Later in this section we shall present some generalizations to the general $`N`$ case.

### A Two-qubit Hamiltonians

Consider two qubits, $`A`$ and $`B`$, that naturally evolve according to some Hamiltonian $`H_{AB}`$. The Hamiltonians achievable by interspersing the evolution of $`H_{AB}`$ with control operations $`w_kw_k`$, $`w_kSU(2)`$, are now

$$H_{AB}^{}=c\underset{k}{}p_kv_kv_kH_{AB}v_k^{}v_k^{},c0,$$
(4)
where we have followed the same steps than in . Because the set $`\{I,\sigma _1,\sigma _2,\sigma _3\}`$ forms a basis for the self-adjoint operators acting on $`𝒞^2`$, we can always expand $`H_{AB}`$ as

$$H_{AB}=\alpha II+h_AI+Ih_B+\underset{i,j=1}{\overset{3}{}}(M)_{ij}\sigma _i\sigma _j,$$
(5)
where $`M`$ is a real $`3\times 3`$ matrix called the “Pauli representation” of $`H_{AB}`$, and $`h_A`$ and $`h_B`$ are some local terms, that we set to $`h_A=h_Bh`$ for simplicity . The term $`\alpha II`$ is an irrelevant phase and we ignore it. We can always add or eliminate equal local terms by applying an appropriate HLU, $`e^{iht}e^{iht}`$, because

$$\left(e^{iht}e^{iht}\right)e^{iHt}=e^{i(hI+Ih+H)t}+𝒪(t^2).$$
(6)
Therefore, we can cancel the local terms of $`H_{AB}`$, simulate the interaction part of $`H_{AB}^{}`$ using the interaction part of $`H_{AB}`$ and, then, add the local terms of $`H_{AB}^{}`$ with the corresponding HLU. Hence, we need consider Hamiltonians with only interaction terms. For any unitary $`uSU(2)`$ there always exists a rotation $`RSO(3)`$ such that

$$u\sigma _iu^{}=\underset{j}{}R_{ij}^T\sigma _j.$$
(7)
It is worth mentioning that this argument only works for qubits . Using this, we write (4) as

$$M^{}=c\underset{k}{}p_kR_kMR_k^T.$$
(8)
The resulting $`M^{}`$ is thus proportional to an orthogonal (real unitary) mixing of $`cM`$. Every term $`R_kMR_k^T`$ in Eq. (8) conserves the trace and the symmetry of $`M`$. This implies that

$$tr(M^{})=ctr(M).$$
(9)
Thus, we have no freedom in choosing the time overhead $`c`$ (inverse of the efficiency $`s`$ of ) unless $`tr(M)=tr(M^{})=0`$. Another consequence of Eq. (8), is that with a symmetric (antisymmetric) $`M`$, $`M=M^T`$ ($`M=M^T`$), we can only simulate a symmetric (antisymmetric) $`M^{}`$. Let us study separately these two cases.

#### 1 Symmetric M

Let us start by considering a Hamiltonian symmetric under the exchange of the subsystems, $`SH_{AB}S^{}=H_{AB}`$, where the swap operator $`S`$ is defined as to exchange the states of qubit $`A`$ and $`B`$, $`S(|\psi _A|\varphi _B)=|\varphi _A|\psi _B`$. We note that this is the case of many relevant interactions, such as the Ising interaction $`\sigma _x\sigma _x`$, the $`XY`$ model or anisotropic exchange interaction $`\sigma _x\sigma _x+\sigma _y\sigma _y`$, and the Heisenberg or isotropic exchange interaction $`\sigma _x\sigma _x+\sigma _y\sigma _y+\sigma _z\sigma _z`$.

Symmetric real matrices are the only ones that can be diagonalized with $`SO(3)`$ transformations. Calling $`\stackrel{}{\lambda }`$ ($`\stackrel{}{\lambda }^{}`$) the vector of eigenvalues of $`M`$ ($`M^{}`$) we can rewrite (8) as

$$\lambda _i^{}=c\underset{k}{}p_k\underset{j}{}(R_k)_{ij}^2\lambda _j,$$
(10)
where we have absorbed the orthogonal matrices that diagonalize $`M`$ and $`M^{}`$ in redefinitions of the $`R_k`$ matrices. Doubly stochastic matrices form a convex set with the permutation matrices as the generators (see Ref. ). The squared elements of an orthogonal matrix, $`(R)_{ij}^2`$, always form a doubly stochastic matrix. Thus, when Eq. (10) is satisfied there exists a doubly stochastic matrix, $`D`$, such that

$$\stackrel{}{\lambda }^{}=cD\stackrel{}{\lambda }.$$
(11)
On the other hand, for every permutation matrix $`P`$ there exist orthogonal matrices $`R`$ fulfilling $`(R)_{ij}^2=(P)_{ij}`$. Therefore, given a doubly stochastic matrix $`D`$ we always can find an orthogonal mixing such that $`(D)_{ij}=_kp_k(R_k)_{ij}^2`$. Concluding, $`M`$ simulates $`M^{}`$ iff there exists one doubly stochastic matrix $`D`$ fulfilling Eq. (11). The last is equivalent to say that $`c\stackrel{}{\lambda }`$ majorizes  $`\stackrel{}{\lambda }^{}`$ (the proof is in ).

Result 1: *Let $`H_{AB}`$ and $`H_{AB}^{}`$ be two symmetric interactions. $`H_{AB}`$ simulates $`H_{AB}^{}`$ by using fast HLU if, and only if, the eigenvalues of $`M^{}`$ are majorized by the eigenvalues of $`cM`$, with the appropriate time overhead $`c`$ given by Eq. (9)*.

Using the notation of , we write Result 1 as

$$H_{AB}^{}_{\mathrm{HLU}}cH_{AB}\stackrel{}{\lambda }^{}c\stackrel{}{\lambda }.$$
(12)
Notice that the above simulation strategy is automatically optimal because the time overhead $`c`$ is fixed.

The above result can be understood in physical terms as follows. The more isotropic a Hamiltonian is, the less simulating power it carries. External homogeneous manipulations can only produce mixing, thus isotropy. More precisely, the isotropy of $`H_{AB}`$ —as measured by the degree of mixing of $`\stackrel{}{\lambda }`$— sets the arrow of allowed simulations by HLU. Accordingly, highly anisotropic Hamiltonians, e.g. the Ising model type $`H_{AB}=\sigma _z\sigma _z`$, are excellent resources to simulate other Hamiltonians and, later on, to synthesize efficiently all sort of non-local gates. Instead, the isotropic exchange interaction $`H_{AB}=\sigma _x\sigma _x+\sigma _y\sigma _y+\sigma _z\sigma _z`$ is a fixed point under HLU transformations, and can not simulate any other interaction.

#### 2 Time–optimal simulation with antisymmetric M

We consider now the case of an interaction fulfilling $`SH_{AB}S^{}=H_{AB}`$. There is an isomorphism between $`3\times 3`$ antisymmetric matrices and 3-dimensional vectors via $`v(M)_i=ϵ_{ijk}M_{jk}`$. Due to the invariance of the Levi-Civita tensor $`ϵ_{ijk}`$ under rotations, we obtain $`\stackrel{}{v}(RMR^T)=R\stackrel{}{v}(M)`$. Then, we can write Eq. (8) in vectorial form:

$$\stackrel{}{v}(M^{})=c\underset{k}{}p_kR_k\stackrel{}{v}(M).$$
(13)
By the triangular inequality we have the following relation among the two moduli

$$v(M^{})cv(M).$$
(14)
Result 2: *Let $`H_{AB}`$ and $`H_{AB}^{}`$ be two antisymmetric interactions. $`H_{AB}`$ can always simulate $`H_{AB}^{}`$ using fast $`HLU`$ for any time overhead $`c`$ fulfilling (14)*, i.e.,

$$H_{AB}^{}_{\mathrm{HLU}}cH_{AB}v^{}cv.$$
(15)
The minimal time overhead or optimal simulation is achieved when the equality (14) holds, $`c=v^{}/v`$. An interesting consequence of Result 2 is that any Hamiltonian with antisymmetric $`M`$ can simulate itself backwards in time with overhead $`1`$, also that it can simulate the zero Hamiltonian. In other words, we can always stop its action. We shall use this property later on.

We finish our analysis for two qubits by noticing that the above cases correspond to the general decomposition of a rank two tensor in irreducible representations of the rotation group, $`J=0`$, $`J=1`$ and $`J=2`$. The scalar pieces produces the Eq. (9) restriction, the antisymmetric case comes from $`J=1`$ and gives Eq. (13) and the standard majorization is related to the $`J=2`$ irreducible representation.

### B Multi-qubit systems with two-body interactions

Let us now consider a system of $`N`$ qubits with a Hamiltonian $`H_{12\mathrm{}N}`$ made out of two-qubit interactions,

$$H_{12\mathrm{}N}=\underset{i,j;j>i}{}H_{ij}.$$
(16)
If interactions $`H_{ij}`$ are arbitrary, the characterization of allowed simulations using $`H_{12\mathrm{}N}`$ and HLU becomes excessively involved. However, in those systems where, e.g., $`H_{ij}=H`$ for some fixed interaction $`H`$, we automatically recover the Results 1 and 2 discussed above. This also happens if some couples of qubits do not interact at all.

As an example, we may consider an array of spins with first-neighbour interaction $`\sigma _z\sigma _z`$. By applying convenient HLU, we can make the system evolve according to first-neighbour anisotropic (or isotropic) exchange interactions. More generally, we can achieve any first-neighbour interaction compatible with Eq. (12).

## III  Time-optimal synthesis of non-local gates using HLU

The problem of producing a two-qubit gate $`U`$ in optimal time using a given interaction $`H_{AB}`$ and fast LU was throughoutly explored in  and finally solved in . In this section we study this problem in the case where the given interaction $`H_{AB}`$ is symmetric ($`SH_{AB}S^{}=H_{AB}`$) and only HLU are allowed. Notice that symmetric interactions are the most common ones in nature. Because HLU are also symmetric, the net effect of interspersing the action of $`H_{AB}`$ with HLU must give a symmetric $`U`$. Thus, we only consider unitary gates $`U`$ fulfilling

$$SUS^{}=U.$$
(17)
### A Synthesis of symmetric two-qubit gates

In what follows we take advantage of a result of , which says that any unitary matrix $`U`$ acting on the Hilbert space $`𝒞^2𝒞^2`$ can be written as

$$U=u_Au_B\mathrm{exp}\left(i\underset{k=1}{\overset{3}{}}\lambda _k^{}\sigma _k\sigma _k\right)v_Av_B,$$
(18)
the decomposition being unique if we make the restriction

$$\pi /4\lambda _1^{}\lambda _2^{}|\lambda _3^{}|.$$
(19)
The recipe for getting $`\stackrel{}{\lambda }^{}`$ is also given in . When imposing Eq. (17) we get that the one-qubit unitaries $`u_A`$ and $`v_A`$ must be equal to $`u_B`$ and $`v_B`$ respectively, but since now we cannot modify the sign of the coefficients $`\lambda _k`$, a unique decomposition for $`U`$ arises only if we require

$$\pi /4|\lambda _1^{}||\lambda _2^{}||\lambda _3^{}|.$$
(20)
Because we assume HLU can be performed arbitrarily fast and we are only concerned with simulating times, simulating $`U`$ is equivalent to just simulating

$$U_\stackrel{}{\lambda }^{}\mathrm{exp}\left(i\underset{k=1}{\overset{3}{}}\lambda _k^{}\sigma _k\sigma _k\right).$$
(21)
The construction $`U_\stackrel{}{\lambda }^{}`$ is called the canonical form of $`U`$. Note that all commutators $`[\sigma _j\sigma _j,\sigma _k\sigma _k]`$ vanish and that $`\mathrm{exp}(i\pi \sigma _k\sigma _k/2)=i\sigma _k\sigma _k`$ is a HLU. This implies that, for any vector $`\stackrel{}{n}=(n_1,n_2,n_3)`$ with integer components $`n_k`$, $`U_{\stackrel{}{\lambda }^{}+\pi \stackrel{}{n}/2}=U_\stackrel{}{\lambda }^{}U_{\pi \stackrel{}{n}/2}`$, which is equivalent to $`U_\stackrel{}{\lambda }^{}`$ up to HLU. Accordingly, we can associate to $`U`$ any one of the Hamiltonians

$$H_{\stackrel{}{\lambda }^{}+\pi \stackrel{}{n}/2}\underset{k=1}{\overset{3}{}}\left(\stackrel{}{\lambda }_k^{}+\frac{\pi }{2}\stackrel{}{n}\right)\sigma _k\sigma _k.$$
(22)
Now, if we are able to simulate the action of one of these Hamiltonians for a time $`t^{}=1`$, we can implement the gate $`U`$ on our two-qubit system, since $`U`$ and $`U_{\stackrel{}{\lambda }^{}+\pi \stackrel{}{n}/2}`$ are HLU–equivalent. In addition, because $`H_{AB}`$ and $`H_{\stackrel{}{\lambda }^{}+\frac{\pi }{2}\stackrel{}{n}}`$ are both symmetric, Result 1 sets the simulation condition:

$$\stackrel{}{n}|(\stackrel{}{\lambda }^{}+\frac{\pi }{2}\stackrel{}{n})c\stackrel{}{\lambda }.$$
(23)
That is, even if $`H_{AB}`$ is not able to simulate $`H_\stackrel{}{\lambda }^{}`$ with HLU, it may still be able to simulate $`H_{\stackrel{}{\lambda }^{}+\pi \stackrel{}{n}/2}`$ for some vector $`\stackrel{}{n}`$, and thus perform gate $`U`$. Notice, however, that in the case when $`\stackrel{}{\lambda }`$ is proportional to $`(1,1,1)`$ (isotropic interaction) there is no integer vector $`\stackrel{}{n}`$ that makes the majorization possible unless $`\stackrel{}{\lambda }^{}`$ is itself proportional to $`(1,1,1)`$. In the next paragraph we explore for which $`\stackrel{}{\lambda }`$s the condition (23) is satisfied whatever the value of $`\stackrel{}{\lambda }^{}`$ is.

The sum of the components of $`\stackrel{}{\lambda }`$, $`\mathrm{\Sigma }_k\lambda _k`$, can be positive, negative or zero. If $`\mathrm{\Sigma }_k\lambda _k`$ is positive, by choosing $`n`$ large enough, we can make $`\stackrel{}{\lambda }^{}+\frac{\pi }{2}(n,n,n)`$ as close as we want to a vector proportional to $`(1,1,1)`$, which is always majorized by any $`c\stackrel{}{\lambda }`$ not proportional to $`(1,1,1)`$. The same happens if $`\mathrm{\Sigma }_k\lambda _k`$ is negative, but in this case, $`n`$ must be large and negative. If $`\mathrm{\Sigma }_k\lambda _k=0`$, $`c\stackrel{}{\lambda }`$ only majorizes vectors having also the sum of their components equal to zero, and, not always exists an integer vector $`\stackrel{}{n}`$ making $`\mathrm{\Sigma }_k(\lambda _k^{}+\pi n_k/2)=0`$. We conclude that, interacting hamiltonians with $`\mathrm{\Sigma }_k\lambda _k`$ different from zero and $`\stackrel{}{\lambda }`$ not proportional to $`(1,1,1)`$ can be used to synthesize any symmetric gate.

Notice, from Eq. (23) that increasing the modulus of $`\stackrel{}{\lambda ^{}}+\pi \stackrel{}{n}/2`$ by changing $`\stackrel{}{n}`$ carries an increase of the time-overhead $`c`$ (equivalently, decrease of the efficiency $`s`$ of ). Thus, we can conlcude that the more isotropic $`H_{AB}`$ is (i.e. the closer $`\stackrel{}{\lambda }`$ is to a vector proportional to $`(1,1,1)`$) the less efficient $`H_{AB}`$ is to synthesize gates. As a limiting case, if $`\stackrel{}{\lambda }`$ is precisely proportional to $`(1,1,1)`$, it will never majorize any $`\stackrel{}{\lambda }^{}`$ not proportional to $`(1,1,1)`$.

### B Time-optimal gate synthesis

So far we have seen that any symmetric gate can be synthesized by almost all interactions. A priori, the best synthesis protocol ought not to make use of a unique simulated interaction since there is freedom to change it along the protocol. Yet, the following result states that, indeed, time-optimal gate synthesis is achieved maintaining a unique simulated Hamiltonian. As before, $`\stackrel{}{\lambda }`$ and $`\stackrel{}{\lambda }^{}`$ are the vectors associated to $`H_{AB}`$ and $`U`$, respectively.

Result 3 (Theorem 10 in ): *The time-optimal way to synthesize a symmetric gate $`U`$ using the interaction $`H_{AB}`$ and fast HLU consists of simulating, for a time $`t^{}=1`$, among all Hamiltonians $`H_{\stackrel{}{\lambda }^{}+\frac{\pi }{2}\stackrel{}{n}}`$ such that*

$$(\stackrel{}{\lambda }^{}+\frac{\pi }{2}\stackrel{}{n})c\stackrel{}{\lambda }$$
(24)
*the one with smallest time overhead $`c`$. The minimal interaction time (interaction cost of ) is given by $`c`$.*

The proof of this result follows from particularizing theorem 10 in Ref.  as indicated in , when taking into consideration that, up to permutations of the coefficients $`\lambda _k`$ , Eq. (22) contains all Hamiltonians $`H_{AB}^{}`$ such that $`U`$ is HLU–equivalent to $`\mathrm{exp}(iH_{AB}^{})`$.

Let us illustrate the above construction with a specific example. We would like to synthesize the non-local part of the CNOT gate, $`|i_A|j_B|i_A|ij_B`$, by using a physical system which interacts according to the anisotropic exchange interaction: $`H=\sigma _x\sigma _x+\sigma _y\sigma _y`$, that is $`\stackrel{}{\lambda }=(1,1,0)`$. We first note that the CNOT gate is equal to:

$`U_{\mathrm{CNOT}}=\left(e^{i\frac{\pi }{4}\sigma _z}U_{\mathrm{Hadamard}}e^{i\frac{\pi }{4}\sigma _z}\right)e^{i\frac{\pi }{4}\sigma _z\sigma _z}`$ (25)
$`\left(IU_{\mathrm{Hadamard}}\right),`$ (26)
that is, its interaction part reduces, up to HLU, to

$$U_{(\pi /4,0,0)}=\mathrm{exp}(i\frac{\pi }{4}\sigma _x\sigma _x).$$
(27)
A direct simulation of this Hamiltonian cannot be achieved since the condition $`(\pi /4,0,0)c(1,1,0)`$ cannot be fulfilled for any $`c`$ (vector $`(1,1,0)`$ is already more isotropic than $`(\pi /4,0,0)`$). Instead, we try to simulate the Hamiltonian $`H_{(\pi /4,\pi /2,\pi /2)}`$, by which we would also achieve the goal according to (22). Now the simulation is possible, $`(\pi /4,\pi /2,\pi /2)c(1,1,0)`$, with time overhead $`c=5\pi /8`$. It can be seen that this is the time-optimal way to do this. A possible mixing of permutations of the components of $`\stackrel{}{\lambda }=(1,1,0)`$ leading to $`(\pi /4,\pi /2,\pi /2)`$ is

$$\left(\begin{array}{c}\pi /4\\ \pi /2\\ \pi /2\end{array}\right)=c\left[\frac{1}{5}\left(\begin{array}{c}1\\ 1\\ 0\end{array}\right)+\frac{1}{5}\left(\begin{array}{c}1\\ 0\\ 1\end{array}\right)+\frac{3}{5}\left(\begin{array}{c}0\\ 1\\ 1\end{array}\right)\right].$$
(28)
Note that to permute the $`x`$-$`z`$($`y`$-$`z`$) axes we can perform a rotation of $`\pi /2`$ along the $`y`$($`x`$) axis. Summing up, in order to synthesize the non-local part of the CNOT gate we have to divide the time $`t=5\pi /8`$ into small equal intervals $`ϵ`$ and, in each one of those intervals, we have to perform the following sequence:

$`\left(e^{i\frac{\pi }{4}\sigma _x}e^{i\frac{\pi }{4}\sigma _x}\right)e^{iH\frac{1}{5}ϵ}\left(e^{i\frac{\pi }{4}\sigma _x}e^{i\frac{\pi }{4}\sigma _x}\right)`$ (29)
$`\left(e^{i\frac{\pi }{4}\sigma _y}e^{i\frac{\pi }{4}\sigma _y}\right)e^{iH\frac{1}{5}ϵ}\left(e^{i\frac{\pi }{4}\sigma _y}e^{i\frac{\pi }{4}\sigma _y}\right)`$ (30)
$`e^{iH\frac{3}{5}ϵ}`$ (31)
The gate is synthesized once this protocol is repeated for a time $`t=5\pi /8`$. A more detailed explanation of the type of manipulations needed above is given in .

## IV  Hamiltonian simulation and gate synthesis with an Asymmetric interaction and HLU

Let us now consider how to simulate symmetric interactions and synthesize symmetric gates by using HLU and a two-qubit interaction $`H_{AB}`$ which is neither symmetric nor antisymmetric under exchange of qubits.

Recall that we can associate to $`H_{AB}`$ a real matrix $`M`$ (its Pauli representation) through Eq. (5). In what follows we give a method which gets rid of the antisymmetric part of $`H_{AB}`$, $`M_a(MM^\mathrm{T})/2`$, and the local terms, $`h_AI`$ and $`Ih_B`$, without affecting the symmetric part $`M_s(M+M^\mathrm{T})/2`$. To achieve this goal we use the following property:

$$\underset{k=0}{\overset{3}{}}\frac{1}{4}\sigma _k\sigma _k\left(\sigma _n\sigma _m\right)\sigma _k^{}\sigma _k^{}=\delta _{nm}\sigma _n\sigma _n,$$
(32)
where $`n`$ and $`m`$ take the values $`0,1,2,3`$ and $`\sigma _0`$ is the identity matrix. Performing this unitary mixing on a general Hamiltonian $`H`$, we get

$$\underset{k=0}{\overset{3}{}}\frac{1}{4}\sigma _k\sigma _kH\sigma _k^{}\sigma _k^{}=\underset{k=1}{\overset{3}{}}(M)_{kk}\sigma _k\sigma _k.$$
(33)
That is, only the diagonal part of $`M`$ survives. Suppose that $`Q`$ is the $`SO(3)`$ matrix that diagonalizes $`M_s`$, and $`u`$ is the unitary matrix corresponding to $`Q`$ via the morphism (7). Then, if instead of (33), we perform the mixing

$`{\displaystyle \underset{k=0}{\overset{3}{}}}{\displaystyle \frac{1}{4}}(u^{}\sigma _ku)(u^{}\sigma _ku)H(u^{}\sigma _ku)^{}(u^{}\sigma _ku)^{}`$ (34)
$`={\displaystyle \underset{i,j=1}{\overset{3}{}}}(M_s)_{ij}\sigma _i\sigma _j,`$ (35)
we get the projection onto the symmetric interaction part of the Hamiltonian and no local terms, with no loss of efficiency. In order to make a simulation using only the symmetric interaction part, we divide every infinitesimal time between the action of two HLUs into four equal intervals and perform the mixing (35). Therefore,

Result 4: *Results 1 and 3 also hold in the case when the simulating Hamiltonian $`H_{AB}`$ is not symmetric.*

## V  Universal Quantum computation with an hybrid model of simulation

We have already argued that an interaction and HLU are not sufficient to perform universal quantum computation, since not all gates can be achieved using only these resources. In this last section we analize which extra resources are needed so that the HLU model recovers universality. As extra resources we consider inhomogeneous LU. A context where this hybrid model may apply is the case of two physical qubits which are brought very close to each other so that they can interact. When the qubits are close, only HLU can be performed in order to control the evolution, and the HLU model applies. But when the systems are brought far apart and no longer interact, they can be addressed individually and therefore arbitrary LU can be applied.

We have seen in sect. III that any symmetric interaction (a part from the isotropic interaction and interactions with tr$`M=0`$) can be used to perform any symmetric gate. This result has been extended in sect. IV also to asymmetric two-qubit interactions (with $`M_s`$ not proportional to the identity and tr$`M0`$). Bearing this in mind, we announce the following result.

Result 5: *Any two-qubit gate can be performed using HLU and any given interaction $`H_{AB}`$ (such that $`M_s`$ is neither trace-less nor proportional to the identity), provided that we allow for two non-homogeneous LU*.

To justify the last assertion, we recall that any two-qubit gate $`U`$ can be written as in Eq. (18). Thus, we can use the available interaction $`H_{AB}`$ to perform the non-local part of $`U`$, which between two non-homogeneous LU interventions leads to $`U`$.

This method is very inefficient when $`M_s`$ is almost proportional to the identity matrix or has very small trace, but alternative methods can be developed in these specific cases, and the hybrid model turns out to be fully universal. For instance, the isotropic interaction $`\stackrel{}{\lambda }=(1,1,1)`$ can be used to perform the gate given by $`\stackrel{}{\lambda }^{}=\pi /8(1,1,1)`$. Then, we can compose two of these gates with a local unitary to achieve a $`(\pi /4,0,0)`$ gate, that is the non–local part of a CNOT gate, which together with LU is known to be sufficient for universal quantum computation.

## VI  Conclusions

We have studied the simulation of non–local Hamiltonians and the synthesis of non-local gates when only HLU can be used to control the available interaction. This model is motivated by the difficulties that often arise in solid state physics, where individual systems can not be addressed independently because they are very close to each other.

Although notoriously less powerful than inhomogeneous manipulation, homogeneous local unitary transformations are sufficient to establish simulations within broad classes of Hamiltonians arranged by their symmetry. In particular, HLU can be used to dispose of the anisotropic terms of any interaction, as required in . Symmetric gate synthesis can be performed using symmetric interaction Hamiltonians and LHU. Moreover, these resources can be supplemented with just a few extra non-homogeneous transformations to fully simulate any two-qubit unitary transformation, and thereby allow for universal quantum computation.

We thank E. Jané and A. Acín for their useful comments. We thank financial support from the following projects: AEN99-0766, 1999SGR-00097, IST-1999-11053. G.V. acknowledges grant HPMF-CT-1999-00200 (Marie Curie fellowship) by the European Community. This work was supported in part by the National Science Foundation (of U.S.A.) under Grant No. EIA-0086038.