text:
# Microarcsecond Radio Imaging using Earth Orbit Synthesis

## 1 Introduction

Flux density variability on intra-day timescales has recently been detected in a few tens of extragalactic radio sources (Heeschen 1984; Witzel et al. 1986, Kedziora-Chudczer et al. 2001). There is now overwhelming evidence that the principal mechanism responsible for intra-day variability (IDV) at centimeter wavelength is interstellar scintillation (ISS), rather than intrinsic source variability (Kedziora-Chudczer et al. 1997; Dennett-Thorpe & de Bruyn 2001a,b; Jauncey et al. 2001; Rickett et al. 2001; Jauncey & Macquart 2001).

ISS is observed in compact ($`100\mu `$as) radio sources due to the scattering of radiation in the Galactic interstellar medium (ISM) as it propagates toward earth. Inhomogeneities in the ISM impart phase structure on the wavefront, which leads in turn to a pattern of flux density deviations on the observer’s plane. Flux density variability is observed in a scattered source if an observer moves relative to the source’s scintillation pattern, ie. due to relative motion between the ISM, the source and the earth. The presence of such radio IDV implies extremely small source sizes ($`5100\mu `$as), and high brightness temperatures ($`T_B10^{13}10^{15}`$ K, e.g. Kedziora-Chudczer et al. 1997). This is well above the inverse Compton limit, and it is not clear whether such emission is explainable in terms of Doppler boosted synchrotron emission (Begelman, Rees & Sikora 1994). The properties — including the radio structures — of IDV sources are therefore of considerable interest.

The scintillation pattern of an IDV source contains information on the structure of the source, which can be derived from the statistical properties of the scintillation pattern given the statistical properties of the scattering medium. This has already been demonstrated for sources undergoing interplanetary scintillation; the experiment by Cornwell, Anantharamaiah & Narayan (1989) is an elegant illustration of this technique.

Some IDV sources exhibit large annual modulations in their variability timescales (Dennett-Thorpe & de Bruyn 2001a; Rickett et al. 2001; Jauncey & Macquart 2001). This is explained in terms of the earth’s orbital motion; as the earth’s velocity changes during its orbit about the sun, so does an earth-based observer’s velocity relative to the ISM. Both the speed and direction change greatly if the speed of the scintillation pattern relative to the heliocentre, due to the intrinsic motion of the ISM, is comparable to the earth’s speed. This effect is expected to be prevalent for IDV sources because the ISM velocity is comparable to earth’s, but is not generally observed in the scintillations of pulsars because of their high intrinsic speeds. Their speeds, typically $`150`$ km/s are much greater than that of the earth or the ISM, and usually dominate the scintillation speed.

The large variation in scintillation speed and direction over the course of a year makes it possible to observe two-dimensional structure in the scintillation pattern. Hence ISS can be used to “image” a source on $`5100\mu `$as scales. In this paper we show how to use changes in the earth’s velocity around the sun to probe the scintillation pattern in two dimensions, and how to invert this information to obtain source structure. We refer to this technique as Earth Orbit Synthesis in light of its similarity to the conventional radio interferometric imaging technique of Earth Rotation Synthesis.

As for the radio imaging technique of Earth Rotation Synthesis, Earth Orbit Synthesis assumes that the source structure remains relatively stable over the course of the observations. This appears to be a good assumption for several IDV sources that have been well-studied over a number of years. B0917$`+`$624 appears to have been scintillating over the decade and a half that this source has been followed in detail (e.g. Jauncey & Macquart 2001), while J1819$`+`$3845 has shown a clear annual cycle that repeats closely over at least the last two years (Dennett-Thorpe & de Bruyn 2001a), and PKS 1519$``$273 has shown IDV each time it has been observed in sufficient detail over the past decade (Kedziora-Chudczer et al. 2001). However, it should be remembered that this technique is not applicable to those sources that exhibit episodic IDV such as PKS 0405$``$385 (Kedziora-Chudczer et al. 2000).

In §2 the apparent scintillation velocity of an extragalactic source is derived in terms of the velocities of the ISM and the earth. In §3 we relate the statistics of intensity and polarization fluctuations to source structure. A model for the scintillation velocity is combined with statistics of the intensity fluctuations to derive source structure in §4, and the technique is demonstrated with some specific examples. The conclusions are presented in §5.

## 2 Earth’s orbital motion

The velocity of an earth-based observer across the scintillation pattern of a scattered source varies annually due to the earth’s orbit about the sun. Here we derive the observed scintillation velocity in terms of the position of a distant scintillating source and the velocity of the interstellar medium relative to the heliocentre (see Fig. 1).

We assume that the internal motions of the scattering material are negligible compared to the bulk motion of the material across the line of sight. We also assume that the scattering occurs predominantly in one layer of the ISM only, at a distance $`L`$ from earth, so that motion of the scattering material can be characterized by a single velocity. This assumption appears to be valid for a few IDV radio sources whose scattering screens appear to be nearby and localized (e.g. Dennett-Thorpe & de Bruyn 2000), and is a good approximation for other IDV sources.

The observed velocity of the line of sight to an extragalactic source through a scattering screen is given by

$`𝐯_{\mathrm{obs}}=𝐯_{\mathrm{ISM}}𝐯_{},`$ (2-1)
where $`𝐯_{\mathrm{ISM}}=(v_{\mathrm{ISM}}^{}{}_{x}{}^{},v_{\mathrm{ISM}}^{}{}_{y}{}^{},v_{\mathrm{ISM}}^{}{}_{z}{}^{})`$ and $`𝐯_{}`$ are, respectively, the velocities of the ISM and the earth relative to the heliocentre. For extragalactic sources, the proper motions of the source across the line of sight can be ignored. The co-ordinate system is aligned so that $`𝐞_z`$ points normal to the ecliptic plane and $`𝐞_x`$ points along ecliptic longitude $`0^{}`$ in the ecliptic plane.

We ignore the $`0.01^{}`$ inclination of the earth’s orbit to the ecliptic, and approximate the earth’s orbital motion as occurring in the ecliptic plane. The eccentricity of the earth’s orbit ($`e=0.0167`$) results in a negligible variation in orbital speed. We therefore approximate the earth’s orbital velocity relative to the heliocentre as:

$`𝐯_{}`$ $`=`$ $`v_{}(\mathrm{sin}2\pi t,\mathrm{cos}2\pi t,0),`$ (2-2)
where $`v_{}=29.78`$km/s is the mean orbital speed of the earth, $`t`$ is the time in years measured from the vernal equinox.

As shown in Figure 1, and choosing the earth to be at the origin of the co-ordinate system we write location of the scintillation screen in terms of ecliptic co-ordinates $`𝐚=a(\mathrm{cos}\varphi \mathrm{cos}\theta ,\mathrm{sin}\varphi \mathrm{cos}\theta ,\mathrm{sin}\theta )`$, where $`\varphi `$ and $`\theta `$ are, respectively, the ecliptic longitude and latitude of the scintillating source. For typical screen distances $`20`$pc, changes in $`\widehat{𝐚}=𝐚/a`$ caused by orbital changes in the earth’s position relative to the screen are negligible.

The observed motion of the scintillation pattern across the line of sight is

$`𝐯_{\mathrm{scint}}=[V_{}+v_{}\mathrm{cos}(2\pi t\varphi )]𝐞_{}+[V_{}+v_{}\mathrm{sin}\theta \mathrm{sin}(2\pi t\varphi )]𝐞_{},`$ (2-3)
where $`𝐞_{}=(\mathrm{sin}\varphi ,\mathrm{cos}\varphi ,0)`$ and $`𝐞_{}=(\mathrm{cos}\varphi \mathrm{sin}\theta ,\mathrm{sin}\varphi \mathrm{sin}\theta ,\mathrm{cos}\theta )`$ are, respectively, the unit vectors parallel and perpendicular to the ecliptic along the line of sight to the scintillating source. The intrinsic speed of the ISM relative to the heliocentre and tangential to the line of sight is $`V_{}`$ along $`𝐞_{}`$ and $`V_{}`$ along $`𝐞_{}`$, where

$`V_{}`$ $`=`$ $`v_{\mathrm{ISM}}^{}{}_{y}{}^{}\mathrm{cos}\varphi v_{\mathrm{ISM}}^{}{}_{x}{}^{}\mathrm{sin}\varphi ,`$ (2-4)
$`V_{}`$ $`=`$ $`v_{\mathrm{ISM}}^{}{}_{z}{}^{}\mathrm{cos}\theta \mathrm{sin}\theta (v_{\mathrm{ISM}}^{}{}_{x}{}^{}\mathrm{cos}\varphi +v_{\mathrm{ISM}}^{}{}_{y}{}^{}\mathrm{sin}\varphi ).`$ (2-5)
Equation (2-3) implies that the velocity of the scintillation pattern relative to the earth changes in both magnitude and direction on an annual cycle.

## 3 Source structure

The scintillation pattern of a scattered source is related to both its intrinsic structure and the inhomogeneities in the ISM. Source structure can be derived if the statistical properties of the scintillation pattern due to a point-source are known. In §3.1 we express for the scintillation pattern due to a point source in terms of the scattering medium’s statistical properties. In §3.2 we show how the scintillation pattern is modified when the scintillating source has extended structure.

### 3.1 The scintillation pattern of a point source

The scintillation pattern of a point source depends upon the phase fluctuations on the scattering screen which are characterized by the phase structure function

$`D_\varphi (𝐫)=[\varphi (𝐫^{}+𝐫)\varphi (𝐫^{})]^2.`$ (3-1)
The angular brackets denote an ensemble average over all possible realisations of the scattering medium, and $`𝐫`$ is a two-dimensional vector on the scattering screen.

For most lines of sight through the ISM the phase structure function reflects an underlying power law spectrum of phase inhomogeneities between some inner and outer scales $`l_0`$ and $`L_0`$ respectively, and one has

$`D_\varphi (𝐫)=\left({\displaystyle \frac{r}{r_{\mathrm{diff}}}}\right)^{\beta 2},`$ (3-2)
where $`\beta `$ is close to $`11/3`$ (Armstrong, Rickett & Spangler 1995), the value expected for Kolmogorov turbulence. The strength of the phase fluctuations is measured by $`r_{\mathrm{diff}}`$, the length scale over which the root-mean-square phase difference on the phase screen is one radian.

The strength of the scattering is determined by the ratio $`r_\mathrm{F}/r_{\mathrm{diff}}`$; the quantity $`r_\mathrm{F}=\sqrt{cL/(2\nu \pi )}`$ is the Fresnel scale. The scattering is weak when phase irregularities on the scattering screen are small, when $`r_{\mathrm{diff}}>r_\mathrm{F}`$ and strong when they are large, when $`r_{\mathrm{diff}}<r_\mathrm{F}`$. The scattering strength varies as a function of frequency since both $`r_{\mathrm{diff}}\nu ^{2/(\beta 2)}`$ and $`r_\mathrm{F}`$ vary with frequency. It is convenient to define $`\nu _t`$ as the frequency at which the scattering changes from strong to weak. The scattering is weak at frequencies above $`\nu _t`$, while strong scattering occurs at frequencies below $`\nu _t`$.

The qualitative properties of the scintillation depend heavily on the scattering strength, and hence on $`\nu _t`$. In general, the value of $`\nu _t`$ toward a scintillating source is unknown, but can be determined from observations of its variability over a range of frequencies (e.g. Kedziora-Chudczer et al. 1997; Quirrenbach et al. 2000; Kedziora-Chudczer et al. 2001), or can be inferred from models of the ISM (Walker 1998). Along lines of sight out of the Galactic plane the transition frequency is in the range $`\nu _t38`$ GHz (Walker 1998).

A useful statistic of the intensity fluctuations of a point source due to scattering by the ISM is the intensity autocorrelation function

$`C_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐫)=\delta I_{\mathrm{pt}}(𝐫^{}+𝐫)\delta I_{\mathrm{pt}}(𝐫^{}),`$ (3-3)
where $`\delta I_{\mathrm{pt}}(𝐫)=[I(𝐫)I]/I`$ is the fractional deviation in the intensity at the point $`𝐫`$ on the observer’s plane. This function describes the correlation between the intensity fluctuations between points $`𝐫^{}`$ and $`𝐫^{}+𝐫`$ on the observer’s plane. An alternative but equivalent representation of the intensity fluctuations is provided by the power spectrum, $`W_{II}^{}{}_{\mathrm{pt}}{}^{}`$, which is related to the intensity autocorrelation function by a Fourier transform

$`W_{II}^{}{}_{\mathrm{pt}}{}^{}={\displaystyle d^2𝐫\mathrm{exp}[i𝐪𝐫]C_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐫)}.`$ (3-4)
We adopt the power spectrum representation in the discussion below.

The functional dependence of $`W_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐫)`$ depends on the scattering strength. In the regime of weak scattering, one has (Salpeter 1967, Jokipii & Hollweg 1970, Narayan 1992)

$`W_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐪)`$ $`=`$ $`Cq^2\left({\displaystyle \frac{1}{qr_{\mathrm{diff}}}}\right)^{\beta 2}\mathrm{sin}^2\left({\displaystyle \frac{1}{2}}r_\mathrm{F}^2q^2\right),`$ (3-5)
with $`C=11.2`$ for $`\beta =11/3`$. The power spectrum peaks at $`q1/r_\mathrm{F}`$, showing that the typical length scale of scintillations in the regime of weak scattering is $`rr_\mathrm{F}`$. An object undergoing weak scintillations would exhibit flux density variations on a timescale $`t_\mathrm{F}=r_\mathrm{F}/v_{\mathrm{scint}}`$.

In the regime of strong scattering the power spectrum of intensity fluctuations for $`\beta <4`$ is (Gochelashvily & Shishov 1975)

$`W_{II}^{}{}_{\mathrm{pt}}{}^{}(q)=\{\begin{array}{cc}2^{\beta 3}\frac{\pi (\beta 2)(4\beta )\mathrm{\Gamma }(\beta /2)}{\mathrm{\Gamma }\left(\frac{6\beta }{2}\right)}q^2\left(\frac{r_\mathrm{F}}{r_{\mathrm{diff}}}\right)^{\beta 2}(qr_\mathrm{F})^{6\beta }\mathrm{exp}[(qr_{\mathrm{ref}})^{\beta 2}],\hfill & qr_{\mathrm{ref}}^1\hfill \\ \frac{2\pi }{\beta 2}\mathrm{\Gamma }\left(\frac{1}{\beta 2}\right)r_{\mathrm{diff}}^2,\hfill & r_{\mathrm{ref}}^1qr_{\mathrm{diff}}^1\hfill \\ \frac{2^\beta \pi (\beta 2)\mathrm{\Gamma }(\beta /2)}{\mathrm{\Gamma }\left(\frac{4\beta }{2}\right)}q^2(qr_{\mathrm{diff}})^{2\beta },\hfill & r_{\mathrm{diff}}^1q\hfill \end{array},`$ (3-9)
where $`r_{\mathrm{ref}}=r_\mathrm{F}^2/r_{\mathrm{diff}}`$ is known as the refractive scale. Intensity fluctuations can occur predominantly on two timescales in this regime, corresponding to the two peaks in $`W_{II}^{}{}_{\mathrm{pt}}{}^{}(q)`$. The outermost peak, corresponding to the effect of diffractive scintillation, occurs at $`q1/r_{\mathrm{diff}}`$ and is associated with fast-timescale, narrowband, large amplitude (up to 100% modulated) variability.

The other peak in $`W_{II}^{}{}_{\mathrm{pt}}{}^{}(q)`$ occurs at $`q1/r_{\mathrm{ref}}`$, and is associated with refractive scintillation. Refractive scintillations occur on a longer timescale than weak scintillations, are broadband and have a typically low ($`10`$%) modulation index. The timescale associated with refractive scintillation is $`t_{\mathrm{ref}}=r_{\mathrm{ref}}/v_{\mathrm{scint}}`$. This is a factor $`r_\mathrm{F}/r_{\mathrm{diff}}>1`$ slower than for weak scattering.

Intensity fluctuations of IDV sources scintillating in the strong scattering regime are seen to occur on timescales $`110`$ days with amplitudes up to $`10`$% of the total source intensity (e.g. Kedziora-Chudczer et al. 2001). This is attributed to refractive scintillation which is exhibited by any source with angular diameter $`\theta _Sr_{\mathrm{ref}}/L`$. No IDV source has yet been observed to exhibit diffractive scintillation. The stringent limitations on source size required for a source to exhibit diffractive scintillation, $`\theta _Sr_{\mathrm{diff}}/L`$ are harder to satisfy for IDV sources. For example, one requires $`\theta _S6\mu `$as at $`\nu =2`$ GHz to exhibit fully modulated diffractive scintillation scattering, where we assume $`\nu _t=5`$ GHz and $`D=100`$ pc. These limitations are less stringent if $`D`$ is smaller, as has been suggested for at least one source (Dennett-Thorpe & de Bruyn 2000).

Expressions (3-5) and (3-9) are only valid when the statistical properties of the phase fluctuations on the scattering screen are isotropic, in the sense that the length scale $`r_{\mathrm{diff}}`$ does not vary with orientation on the screen. Anisotropy introduces a preferred direction in the scattering medium so that $`C_{II}^{}{}_{\mathrm{pt}}{}^{}`$ is no longer a function of $`r`$ only, but depends on the orientation. Typical elongations of the scattering pattern in the ISM can exceed 2:1 (see Chandran & Backer submitted for a full discussion). Anisotropic screen structure can therefore cause elongations in the scintillation pattern similar to apparent elongations caused by intrinsic source structure. One way to distinguish between source structure and screen anisotropy involves imaging the scattered source at low frequency, where the scattering is strong (ie. $`\nu <\nu _t`$) and the condition $`\theta _S<r_{\mathrm{ref}}/D`$ is more likely to be satisfied. Then the shape of the scattered image is entirely dominated by the scattering and any anisotropy is readily apparent as an elongation of the scattering disk along the axis corresponding to small $`r_{\mathrm{diff}}`$. This permits a direct measurement of the effects of screen anisotropy. Of course, the measurement can only be attempted at frequencies low enough that the scattering disk is sufficiently large to be resolved; this may be an important limitation in practice.

We see in §4.1 that the effect of screen anisotropy is not important when the angular size of a source exceeds the angular scale of the scintillation pattern for a point source (ie. $`r_\mathrm{F}/D`$ for weak scattering), and is readily distinguished from source structure by its characteristic structure in $`W_{II}^{}{}_{\mathrm{pt}}{}^{}`$. Thus, to first order, effects of screen anisotropy are diminished in many IDV sources. Hereafter we restrict our discussion to isotropic phase screens for the purposes of simplicity.

It possible to determine which of equations (3-5) and (3-9) is applicable to a set of observations provided that the characteristics of the intensity fluctuations are known at several frequencies. For weak scattering one has $`t_\mathrm{F}\nu ^{1/2}`$, whereas $`t_{\mathrm{ref}}\nu ^{2.2}`$ (assuming $`\beta =11/3`$) for strong scattering. The large difference in the timescales of the intensity variations allows a determination of whether the variability is occurring in the regime of weak or strong scattering.

Many sources are observed to scintillate at frequencies that straddle the regime between weak and strong scattering. It is difficult to provide analytic expressions for $`W_{II}^{}{}_{\mathrm{pt}}{}^{}`$ in this transition regime, and instead, equation (3-5) may be used on the understanding that it is only approximately correct, or $`W_{II}^{}{}_{\mathrm{pt}}{}^{}`$ may be derived numerically from scattering simulations.

### 3.2 The scintillation pattern of an extended source

We now discuss the effect of source structure on the observed scintillation pattern of an extended source. Direct measurements of the scattered wavefield from an extended source may be inverted to completely determine the source structure. The intensity scintillation pattern of a source also contains information on source structure. However, analysis of the intensity fluctuations alone only provides partial information on source brightness distribution.

#### 3.2.1 The wavefield

The wavefield of a scattered source is directly related to its structure. The mutual coherence function is the simplest nontrivial moment of the wavefield. If the wavefield measured by a receiver at the point $`𝐫`$ on the observer’s plane is $`u(z;𝐫)`$, the normalized two-point correlation function is then

$`\mathrm{\Gamma }_{\alpha \alpha ^{}}(z;𝐫)={\displaystyle \frac{u_\alpha (z;𝐫^{})u_\alpha ^{}^{}(z;𝐫^{}+𝐫)}{I_0}},`$ (3-10)
where $`\alpha ,\alpha ^{}=[x,y]`$ corresponds to the component of the wavefield oscillating along the $`x`$ and $`y`$ axes orthogonal to the direction of propagation, along the $`z`$-axis, and $`I_0`$ is the mean source total intensity. The normalized Stokes visibilities are defined by

$`\mathrm{\Gamma }_I(z;𝐫)`$ $`=`$ $`(\mathrm{\Gamma }_{xx}(z;𝐫)+\mathrm{\Gamma }_{yy}(z;𝐫))/2\mathrm{\Gamma }_Q(z;𝐫)=(\mathrm{\Gamma }_{xx}(z;𝐫)\mathrm{\Gamma }_{yy}(z;𝐫))/2`$
$`\mathrm{\Gamma }_U(z;𝐫)`$ $`=`$ $`(\mathrm{\Gamma }_{xy}(z;𝐫)+\mathrm{\Gamma }_{yx}(z;𝐫))/2\mathrm{\Gamma }_V(z;𝐫)=i(\mathrm{\Gamma }_{xy}(z;𝐫)\mathrm{\Gamma }_{yx}(z;𝐫))/2.`$ (3-11)
It is assumed that the source is spatially incoherent and that the wavefront from the IDV source incident upon the scattering screen is planar, which is a very good approximation for extragalactic sources. (The planar approximation is not valid for Galactic objects, but it is possible to make a simple correction for the curvature of the wavefront in this case (e.g., Goodman & Narayan 1989).) The mutual coherence function measured by an observer is then related to the mutual coherence function of the source $`\mathrm{\Gamma }_X(0;𝐫),X=[I,Q,U,V]`$ by

$`\left(\begin{array}{c}\mathrm{\Gamma }_I(z;𝐫)\hfill \\ \mathrm{\Gamma }_Q(z;𝐫)\hfill \\ \mathrm{\Gamma }_U(z;𝐫)\hfill \\ \mathrm{\Gamma }_V(z;𝐫)\hfill \end{array}\right)=e^{D_\varphi (𝐫)/2}\left(\begin{array}{c}\mathrm{\Gamma }_I(0;𝐫)\\ \mathrm{\Gamma }_Q(0;𝐫)\mathrm{cos}2\varphi _V+\mathrm{\Gamma }_U(0;𝐫)\mathrm{sin}2\varphi _V\\ \mathrm{\Gamma }_Q(0;𝐫)\mathrm{sin}2\varphi _V+\mathrm{\Gamma }_U(0;𝐫)\mathrm{cos}2\varphi _V\\ \mathrm{\Gamma }_V(0;𝐫)\end{array}\right),`$ (3-20)
where $`\varphi _V=\lambda ^2\mathrm{RM}`$ is the mean phase change induced by Faraday rotation along the ray path, and we assume no inhomogeneity in the Faraday rotation transverse to the line of sight.

The normalized mutual coherence of the source, $`\mathrm{\Gamma }_X(0;𝐫)`$ is related to its angular brightness distribution normalized by the mean intensity via a Fourier transform

$`\mathrm{\Gamma }_X(0;𝐫)={\displaystyle d^2\text{}b_X(\text{})\mathrm{exp}\left[2\pi i\frac{\text{}𝐫}{\lambda }\right]}.`$ (3-21)
The behaviour of the mutual coherence function depends on the angular size of the source relative to the angular scale of the scintillation pattern as seen by an observer on the earth. Equations (3-20) and (3-21) imply the behaviour of $`\mathrm{\Gamma }(z;𝐫)`$ is dominated by the source angular structure if its angular size exceeds that of the scintillation pattern seen by an observer, $`\theta _Sk/r_{\mathrm{diff}}`$. In this instance the source is said to be resolved by the scattering medium.

Conversely, the mutual coherence function of a source with angular size $`\theta _Sk/r_{\mathrm{diff}}`$ is dominated by the scattering, and $`\mathrm{\Gamma }(z;𝐫)\mathrm{\Gamma }(z;0)\mathrm{exp}[D_\varphi (𝐫)/2]`$. It is difficult to derive source structure when the source is unresolved by the scattering medium because the scattering properties dominate the mutual coherence function. This limitation is important in the regime of strong scattering where $`k/r_{\mathrm{diff}}`$ is may be large relative to the size of an IDV source.

The quantities $`\mathrm{\Gamma }_X(𝐫)`$ can be measured over short coherent integrations (see, e.g. Goodman & Narayan 1989); typically two or three scintles are sufficient to obtain a good estimate.

#### 3.2.2 Intensity fluctuations

Measurement of the intensity fluctuations contains partial information on the structure of a source. In this section we derive the relationship between the brightness distribution of a scattered source and its intensity fluctuations. The pattern of intensity fluctuations from an extended source is derived by dividing the source into small elements and computing the individual intensity scintillation pattern from each element (again assuming a spatially incoherent source). The scintillation pattern of the extended source is therefore the sum of the scintillation patterns from all the elements of the source.

Figure 2 shows that the wavefront incident on the screen at an angle \`  to the line of sight, and heading towards the point $`𝐫+L\text{}`$, intersects a thin scattering screen at the same point as a normally incident wavefront, heading towards $`𝐫`$. Each of these wavefronts experiences identical scattering, and the amplification caused by the scattering screen is thus identical for these two wavefronts. It follows that scintillation pattern of an extended source is (Little & Hewish 1966),

$`\delta S_X(𝐫)={\displaystyle d^2\text{}b_X(\text{})\delta I_{\mathrm{pt}}(𝐫L\text{})},`$ (3-22)
where $`\delta S_X=(S_XS_X)/S_X`$ is the fractional flux density deviation in Stokes parameter $`X=[I,Q,U,V]`$, and $`b_X(\text{})`$ is the angular brightness distribution of the source as observed at the scattering screen, normalized such that $`b_I(\text{})d^2\text{}=1`$ and $`b_Y(\text{})d^2\text{}=Y/I_0,Y=[Q,U,V]`$. The brightness distributions $`b_Q`$ and $`b_U`$ at the screen differ from their distributions at the source if Galactic Faraday rotation is important. The distributions at the source are found by determining the Faraday rotation along the line of sight to the source. A magnetized ISM may also induce fluctuations in the amount of radiation transferred between Stokes $`Q`$, $`U`$ and $`V`$, however the fluctuations in the degree of birefringence are of order $`\mathrm{\Omega }_e/\omega _p10^8`$ smaller than the phase fluctuations considered here, and are probably unimportant for the scattering considered here (see Macquart & Melrose 2000).

The flux density deviations observed in an extended source are connected to the scintillation statistics of a point source by

$`C_{XY}(𝐫)`$ $`=`$ $`\delta S_X(𝐫^{})\delta S_Y(𝐫+𝐫^{})`$ (3-23)
$`=`$ $`{\displaystyle d^2\text{}d^2\text{}^{}b_X(\text{})b_Y(\text{}^{})\delta I_{\mathrm{pt}}(𝐫^{}L\text{})\delta I_{\mathrm{pt}}(𝐫^{}+𝐫L\text{}^{})}.`$
The Fourier-transformed equivalent of equation (3-23) illustrates clearly the connection between source structure and the scintillation statistics via the power spectrum of intensity fluctuations,

$`W_{XY}(𝐪)=\stackrel{~}{C_{XY}}(𝐪)`$ $`=`$ $`\stackrel{~}{S}_X(𝐪)\stackrel{~}{S}_Y^{}(𝐪)`$ (3-24)
$`=`$ $`{\displaystyle \frac{1}{L^4}}[\stackrel{~}{B}_X(𝐪)\stackrel{~}{B}_Y^{}(𝐪)W_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐪)],`$
where $`B_X(\text{})=b_X(\text{}/L)`$, the Fourier transform of a function is denoted with a tilde, and we have made use of equation (3-4). The quantity $`B_X(\text{})`$ is identical to the visibility $`\mathrm{\Gamma }_X(0;𝐫)`$ of the previous section. Equation (3-24) may be inverted to derive information on the structure of a scintillating source if the scattering properties of the ISM, embodied in $`W_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐪)`$, are known.

Equation (3-24) may be interpreted as follows. The quantity $`W_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐪)`$ is the power spectrum of intensity fluctuations in the scintillation pattern that would be observed for a point source. The quantity $`\stackrel{~}{B}_X(𝐪)\stackrel{~}{B}_Y^{}(𝐪)`$ is the power spectrum of intensity structure in the source. The resultant power spectrum of intensity fluctuations, $`\stackrel{~}{C_{XY}}(𝐪)`$, is that of a point source if the angular structure contained in source is smaller than the scale of the scintillation pattern for a point source (ie. the source is unresolved by the ISM). On the other hand, the observed scintillation pattern is dominated by source angular structure if the source is larger than the scale of the scintillation pattern for a point source. For example, scintillations in the weak scattering regime contain power on length scales $`rr_\mathrm{F}`$ (equivalent to $`q1/r_\mathrm{F}`$). Source structure has no measurable effect if the source angular size is less than the angular scale of the scintillation pattern (ie. when $`\theta _S<r_\mathrm{F}/D`$).

If the angular size of the source is smaller than the scale size of the scintillation pattern for a point source, the term in square brackets in (3-24) is approximated by $`W_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐪)`$. This is because $`W_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐪)`$ decreases much faster with $`q`$ than $`\stackrel{~}{B_I}(𝐪)\stackrel{~}{B_I}^{}(𝐪)`$. On the other hand, if the source angular size exceeds that of the scintillation pattern for a point source, we can use the approximation $`W_{II}^{}{}_{\mathrm{pt}}{}^{}2.8q^{4\beta }r_\mathrm{F}^4/r_{\mathrm{diff}}^{\beta 2}`$ for small $`q`$ to obtain (see equation (3-23))

$`C_{XY}(𝐫){\displaystyle \frac{2.8r_\mathrm{F}^4r_{\mathrm{diff}}^{2\beta }}{(2\pi )^2L^4}}{\displaystyle d^2𝐪e^{i𝐪𝐫}q^{4\beta }\stackrel{~}{B_X}(𝐪)\stackrel{~}{B_Y}^{}(𝐪)}.`$ (3-25)
Sources sufficiently large to be resolved by the scattering medium exhibit smaller and slower intensity fluctuations relative to unresolved sources. The intensity modulation index, $`C_{XY}(0)`$, of an resolved source is therefore smaller than that of an unresolved source. This is evident from equation (3-24). The amplitude of $`C_{XY}(0)`$ decreases if the terms due to the source cut off $`C_{XY}(𝐪)`$ before $`W_{II}^{}{}_{\mathrm{pt}}{}^{}`$ reaches its maximum value (e.g. at $`q1/r_\mathrm{F}`$ for weak scattering).

## 4 Earth Orbit Synthesis

The earth’s orbital motion allows the scintillation pattern of a scattered source to be mapped in two dimensions. This is because the velocity of an earth-based observer relative to the ISM changes as the earth’s orbital velocity changes during the course of a year. An observer moving at a velocity $`𝐯`$ relative to the ISM samples the scintillation pattern along the direction parallel to $`𝐯`$. It is possible to measure the correlations $`\mathrm{\Gamma }_X(𝐫=𝐯t)`$ or $`C_{XY}(𝐫=𝐯t)`$ for observations a time $`t`$ apart. The quantity $`\mathrm{\Gamma }_X(𝐫)`$ can be measured over short coherent integrations (see, e.g. Goodman & Narayan 1989). A good approximation to the quantity $`C_{XY}(𝐫)`$ can typically be derived from observations of $`5`$ scintles, typically requiring one to two days flux density monitoring for typical IDV sources.

Over the course of the year observations at different epochs sample the scintillation pattern in different directions, as the velocity of the earth changes around the solar circle. In particular, if the velocity of the ISM is comparable to the earth’s, then a strong “annual cycle” effect will appear in the observed IDV pattern, as has been observed already in J1819+3845 (Dennett-Thorpe & de Bruyn 2001a) and 0917+624 (Jauncey & Macquart 2001; Rickett et al. 2001). This enables measurement of $`\mathrm{\Gamma }_X(𝐫)`$ or $`C_{XY}(𝐫)`$ for a large range of orientations of $`𝐫`$, and the information may be inverted to obtain source structure provided that it is stable on a timescale long enough for it to be mapped. The main limitation to the technique is that one does not always know $`𝐯`$ because the velocity of the ISM is not generally known. Thus one must solve simultaneously for $`𝐯`$ and the source structure, or measure $`𝐯`$ directly from observations of the time delay between the scintillation pattern as observed at two widely separated ($`3000`$ km) receivers as the scintillation pattern passes across the earth (Jauncey et al. 2001, Dennett-Thorpe & de Bruyn 2002).

We demonstrate the technique with several examples, concentrating on Earth Orbit Synthesis using intensity fluctuations.

### 4.1 An example using Intensity Fluctuations

An intensity image of the structure in a scintillating source requires measurement of the autocorrelation function $`C_{II}(𝐫)`$. It is not possible to measure the spatial structure of the scintillation pattern of a source, $`C_{II}(𝐫)`$, directly. The scale of the scintillation pattern is usually larger than that of earth so that it is impossible to measure $`C_{II}(𝐫)`$ on sufficiently long baselines $`𝐫`$. This problem is overcome by measuring the correlations of the scintillation signal as a function of time, rather than baseline. One measures the temporal correlation function

$`C_{II}^{}{}_{t}{}^{}(t)={\displaystyle \frac{I(t^{})I(t+t^{})}{I^2}}1.`$ (4-1)
The correlation between two scintillation signals measured a time $`t`$ apart at the same receiver may be derived from the correlation in the scintillation pattern measured by two receivers a distance $`𝐫`$ apart. The relation between the temporal and spatial intensity autocorrelation functions is thus

$`C_{II}^{}{}_{t}{}^{}(t)=C_{II}(𝐫=𝐯t)={\displaystyle \frac{1}{(2\pi )^2L^4}}{\displaystyle d^2𝐪e^{i𝐪𝐫}\left[\stackrel{~}{B_I}(𝐪)\stackrel{~}{B_I}^{}(𝐪)W_{II}^{}{}_{\mathrm{pt}}{}^{}(𝐪)\right]}.`$ (4-2)
An important point in connection with (4-2) is that the modulation depth of the scintillation pattern, $`C_{II}(0)`$ is independent of the velocity of the scattering medium relative to earth. For instance, the depth of the intensity modulations does not depend on whether the scintillation pattern moves parallel or orthogonal to the long axis of an elongated source. However, we see that the characteristic timescale of the scintillations depends strongly on the direction of motion if $`𝐯`$ varies during the course of a year.

The effect of source structure is particularly simple when a source undergoes weak scattering and its angular size is larger than that of the scintillation pattern of a point source (ie. $`\theta _S>r_\mathrm{F}/D`$). Equation (3-25) is applicable in this case, and we assume $`\beta =4`$ instead of the usual Kolmogorov value of $`\beta =11/3`$ to simplify the algebra and derive analytic results. We obtain

$`C_{II}(𝐫)={\displaystyle \frac{2.8r_\mathrm{F}^4}{r_{\mathrm{diff}}^2L^4}}{\displaystyle d^2𝐪e^{i𝐪𝐫}|\stackrel{~}{B_I}(𝐪)|^2},L\theta _S>r_\mathrm{F}.`$ (4-3)
This expression, and the assumptions made in deriving it, are employed in the following examples.

#### 4.1.1 Elongated source

We illustrate Earth Orbit Synthesis by examining the combined effect of source structure and velocity variations on a source with elliptical brightness distribution:

$`b_I(\text{})={\displaystyle \frac{1}{\alpha _x\alpha _y\pi }}\mathrm{exp}\left[\left({\displaystyle \frac{\theta _x}{\alpha _x}}\right)^2\left({\displaystyle \frac{\theta _y}{\alpha _y}}\right)^2\right],`$ (4-4)
where $`\alpha _x`$ and $`\alpha _y`$ are chosen to be the source angular scales parallel to the $`𝐞_{}`$ and $`𝐞_{}`$ axes respectively. This choice is made for convenience only, and it is simple to re-orient the source at an angle $`\xi `$ to the axes by making the transformation $`(x,y)(x\mathrm{cos}\xi +y\mathrm{sin}\xi ,x\mathrm{sin}\xi +y\mathrm{cos}\xi )`$. Combining equations (4-3) and (4-4) yields the following intensity autocorrelation function

$`C_{II}(𝐫)={\displaystyle \frac{2.8r_\mathrm{F}^4}{2\pi \alpha _x\alpha _yL^2r_{\mathrm{diff}}^2}}\mathrm{exp}\left[\left({\displaystyle \frac{x^2}{2L^2\alpha _x^2}}+{\displaystyle \frac{y^2}{2L^2\alpha _y^2}}\right)\right],L\alpha _{x,y}>r_\mathrm{F}`$ (4-5)
It is possible to determine the temporal decorrelation of the scintillation signal, for a given scintillation speed and direction. In the following we use $`V_{}=15`$ km/s, $`V_{}=10`$ km/s and we assume the source has ecliptic co-ordinates $`\theta =30^{}`$ and $`\varphi =0`$. The variation in the scintillation speed over the course of a year is shown in Fig. 3. Figures 4(a) and 5(a) show how the intensity temporal autocorrelation function $`C_{II}(t)`$ varies throughout the year for an elliptical source with these screen parameters.

The effect of source structure is evident in the timescale of variability. It is convenient to define the decorrelation timescale, $`t_{\mathrm{dc}}`$ as the time in which the intensity correlation falls to $`1/e`$ of its maximum value. The decorrelation timescale for the elliptical source described by equation (4-5) is

$`t_{\mathrm{dc}}=\left[{\displaystyle \frac{(𝐯_{\mathrm{scint}}𝐞_{})^2}{2L^2\alpha _x^2}}+{\displaystyle \frac{(𝐯_{\mathrm{scint}}𝐞_{})^2}{2L^2\alpha _y^2}}\right]^{1/2}.`$ (4-6)
Changes in the magnitude of $`𝐯_{\mathrm{scint}}`$ cause changes in the scintillation timescale even if the source is circularly symmetric, whereas changes in the direction of $`𝐯_{\mathrm{scint}}`$ only result in changes in $`t_{\mathrm{ISS}}`$ if the source is elongated, when $`\alpha _x\alpha _y`$.

Figures 4(b) and 5(b) plot the variation in decorrelation timescale as a function of epoch for a specific scintillation velocity for the elliptical source with two different source parameters, $`L\alpha _x=10^9`$ m and $`L\alpha _y=10^9`$ m and $`5\times 10^9`$ m respectively. Fig. 4(b) shows the expected scintillation timescale for a circularly symmetric source, and Fig. 5(b) the timescale for an elongated source. The peak in the scintillation timescale at $`T0.9`$ yr is common to both plots and occurs because the scintillation speed reaches a minimum. The emergence of another peak in Fig. 5(b) occurs due to the presence of structure in the source, here an elongation along the $`y`$-axis. The peak occurs at the time of the year that the scintillation direction is tangential to the long axis of the source.

#### 4.1.2 Two-component source

A source with two components shows richer structure in the statistics of the intensity variations. We consider two gaussian patches of emission, the central peak having an intensity $`a_0I_0`$ ($`0<a_0<1`$) and the other, located at an angular separation of $`\text{}_0`$, having an intensity $`(1a_0)I_0`$. The two subsources have identical angular scales $`\alpha `$ and we again assume that the angular diameters of the subsources exceed the angular diameter of the scintillation pattern due to a point source, so that equation (3-25) applies in the regime of weak scattering. The source angular brightness distribution is

$`b_I(\text{})={\displaystyle \frac{a_0}{\pi \alpha ^2}}\mathrm{exp}\left[{\displaystyle \frac{\text{}^2}{\alpha ^2}}\right]+{\displaystyle \frac{1a_0}{\pi \alpha ^2}}\mathrm{exp}\left[{\displaystyle \frac{(\text{}\text{}_0)^2}{\alpha ^2}}\right],`$ (4-7)
for which the resulting spatial intensity autocorrelation function is

$`C_{II}(𝐫)={\displaystyle \frac{r_\mathrm{F}^4}{4\alpha ^2L^2\pi ^2r_{\mathrm{diff}}^2}}\mathrm{exp}\left[{\displaystyle \frac{r^2+L^2\text{}_0^2}{2\alpha ^2L^2}}\right]\left((12a_0+2a_0^2)\mathrm{exp}\left[{\displaystyle \frac{\text{}_0^2}{2\alpha ^2}}\right]+2a_0(1a_0)\mathrm{cosh}\left[{\displaystyle \frac{𝐫\text{}_0}{\alpha ^2L}}\right]\right).`$ (4-8)
Figure 6 plots the spatial intensity autocorrelation function of this source. This function has three peaks. The lightcurve of a two-component source therefore exhibits interesting behaviour when the scintillation direction is parallel to the line joining the two substructures. An observer then measures the spatial autocorrelation function along a line that joins the three peaks in the autocorrelation function.

A specific example of such behaviour is illustrated in Figure 7, which demonstrates how the temporal autocorrelation function of a two component source changes throughout the course of a year. At $`T0.00.2`$ yr the scintillation velocity is parallel to the direction of the two separated peaks, and a second bump is visible in the temporal intensity autocorrelation function.

#### 4.1.3 Polarization reference mapping

Only limited information about source structure is available from Earth Revolution Synthesis using intensity fluctuations. It is impossible to determine $`b_I(\text{})`$ directly from $`C_{II}(𝐫)`$. Only the correlation of the angular brightness distribution, $`b_I(\text{}+\text{}^{})b_I(\text{}^{})`$, can be derived.

However, Earth Orbit Synthesis is useful even when the brightness distribution $`b_I(\text{})`$ cannot be ascertained exactly. IDV sources often exhibit fluctuations in Stokes $`Q`$, $`U`$ and $`V`$ that are dissimilar to those observed in $`I`$. This occurs if the polarized structure of the source differs from its structure in total intensity. This substructure can be mapped relative to the total intensity structure, using the correlations $`C_{IQ}(𝐫)`$, $`C_{IU}(𝐫)`$ and $`C_{IV}(𝐫)`$ to find the quantities $`b_I(\text{}^{}+\text{})b_Q(\text{}^{})`$, etc.. It is thus possible to map the brightness distribution of a source in one Stokes parameter relative to another Stokes parameter, as has been shown for both the linear and circular polarized components in PKS 1519$``$273 (Macquart et al. 2000).

## 5 Conclusion

Measurement of the scintillation characteristics of IDV sources allows their structure to be investigated with microarcsecond resolution. The radio structures of IDV sources are currently of interest because of their high brightness temperatures. Earth Orbit Synthesis has the potential to determine the origin of these compact structures; it can determine, for example, whether the observed emission emanates from the central cores of these objects or from the bases of jets.

Short timescale polarization variability in IDV sources indicates the presence of polarization structure on similarly fine scales and — in many cases — on even smaller scales than the structure in total intensity. The advantage of Earth Orbit Synthesis over traditional polarimetric imaging is its high resolution. Such a high resolution probe is less subject to beam-depolarization. Indeed, this mapping technique is already discovering remarkable detail about the polarized sub-structure of IDV sources (e.g. Macquart et al. 2000).

The observed two-dimensional scintillation pattern of a scattered radio source is highly influenced by source structure. This pattern can be measured in two dimensions if the apparent velocity of the scintillation across the line of sight to a scattered source changes; changes in direction allow an observer to probe source structure along different directions. Direction changes occur when the velocity of the scattering material is comparable to the earth’s velocity, and the vector sum of these two velocities changes during the course of a year.

The derivation of source structure requires measurement of the scintillation fluctuations at several epochs over the course of a year. This is achieved by measuring the temporal decorrelation of the intensity $`I(t)I(t+t^{})`$. Many IDV sources also exhibit polarization variability, which is often not well correlated with the total intensity. This occurs when the source structure in polarized and total intensity emission differs. Polarized substructure information is available by measuring the temporal cross-correlations $`X(t)Y(t+t^{})`$ between Stokes parameters $`X`$ and $`Y`$. It is possible to measure the polarized substructure of a source with reference to the intensity emission to high precision.

An advantage of the Earth Orbit Synthesis technique is the capability to compare source structure between frequencies. For example, one could compare the angular separation between source components as a function to yield information relating to opacity changes with frequency. Analysis of the data on PKS 1519$``$273 (Macquart et al. 2000) demonstrates not only the power to compare source structure at two frequencies, but even to compare the relative source properties in different Stokes parameters.

It should be emphasised that the use of Earth Orbit Synthesis as an imaging technique is subject to limitations. Firstly, the scintillation process only probes source structure over a limited angular size range. Consequently, while Earth Orbit Synthesis probes microarcsecond structures, it provides no information on the relation between microarcsecond and larger scale structure. Secondly, the use of autocorrelation functions in the imaging process is a necessary but severe limitation to any imaging that can be achieved. This technique yields the power spectrum of the source brightness distribution and, as such discards phase information that is present in normal rotational synthesis. The stochastic nature of the scintillation process forces us to use ensemble average quantities such as autocorrelation functions.

The angular resolution achievable using scintillation-based imaging depends on the properties of the scattering medium. In the regime of weak scattering, generally applicable at frequencies $`\nu 5`$ GHz, the intensity scintillations are sensitive to angular scales $`\theta r_\mathrm{F}/D\nu ^{1/2}`$. The scattering is strong at frequencies $`\nu 5`$ GHz. The refractive scintillation exhibited by IDV sources in this regime is less sensitive to small angular scales because the angular resolution achievable is $`\theta r_{\mathrm{ref}}/D`$ which is large, and scales as $`\nu ^{2.2}`$ for a Kolomogorov spectrum of inhomogeneities. Strong scattering is typically sensitive to sizes $`0.1`$ mas, and hence can be found in many more flat-spectrum sources, particularly at low Galactic latitudes. Thus, while less useful competing against VLBI for milliarcsecond imaging, strong scattering here remains an excellent probe of the ISM close to the Galactic plane.

A program to use Earth Orbit Synthesis to map structure in IDV sources is currently underway using the Australia Telescope Compact Array (Bignall et al. 2001).

The authors thank Barney Rickett, Simon Johnston and Lucyna Kedziora-Chudczer for valuable discussions and comments on the manuscript. JPM is particularly indebted to Don Melrose, from whose Festschrift the idea for this paper originated. The Australia Telescope is funded by the Commonwealth Government for operation as a national facility by the CSIRO.


format: text

text:
# Molecular Quantum Computing by an Optimal Control Algorithm for Unitary Transformations

## Abstract

Quantum computation is based on implementing selected unitary transformations which represent algorithms. A generalized optimal control theory is used to find the driving field that generates a prespecified unitary transformation. The approach is illustrated in the implementation of one and two qubits gates in model molecular systems.

PACS number(s): 82.53.Kp 33.90.+h 32.80.Qk 03.67.LX

A universal model of a quantum computer can be constructed from an array of two level systems used as registers (qubits). Any general quantum gate can be decomposed to a one and two qubit unitary transformation nielsen00 . Thus a physical realization of a quantum computer should be able control such unitary transformations by use of external driving fields. Realizations based on nuclear magnetic resonance techniques already have been demonstrated for non-trivial quantum gates. For example, a three qubit quantum Fourier transform has been performed by a sequence of pulses corresponding to one and two qubit transformations weinstein01 . This implementation was based on addressing each qubit independently by their spectral separation. However, when many qubits become involved in a computation, the spectrum becomes congested so that it becomes more difficult to address each qubit individually. The total fidelity of the algorithm will depend on the accumulation of errors at each step. Another difficulty results from decoherence processes which are unavoidable due to coupling to the external environment. These processes will degrade the performance in proportion to the time required to carry out the computation task. Thus it is desirable to minimize the number of computation steps and the total computation time.

An alternative possibility for implementing quantum computing is the use of molecules driven by shaped light pulses. A hypothetical model can be based on using the vibrational and rotational states as registers. A shaped light pulse can write data as amplitudes on these states. The computing algorithm is then implemented by performing a unitary transformation employing a second shaped pulse. In the end the output can be read by a probe pulse. This approach implies changing the computation model by tying together many single and two qubit operations into one resulting in a combined unitary transformation.

Such an approach has been used for the experimental implementation of elements of coherent computation in Li2 vala02 ; amitay02 . For example, two states of the ro-vibrational manifold of an electronic state can form a single qubit. By taking advantage of multiple electronic states, the evolution of the molecular system can be controlled using electromagnetic fields in the optical region realized by shaped laser pulses. The fast development of pulse shaping technology can make possible the implementation of simple quantum gates in molecules. This program has to overcome interference from the large number of other molecular levels coupled to the field but not assigned to the qubits. The task becomes therefore to implement the quantum algorithm on the molecular levels used as registers while avoiding the intervention of other states from the same system.

Tesh et. al. tesch01 , proposed the use of optimal control theory (OCT) to calculate the field which can induce a specific transformation used for quantum computation. Originally, OCT was designed as a method to obtain a light field which could induce a specific state to state transformation pierce88 ; kosloff89 . In the quantum computing context, the goal is to obtain the optimal pulse that induces a given unitary transformation, irrespective of the initial state of the system. The present approach generalized OCT to obtain directly the driving field that induces a target unitary transformation in the system. The approach is based on the equation of motion of the unitary transformation by making use of the main ingredients of OCT. This allows the utilization of the large number of tools that have been developed in this field.

The model system consists of a free Hamiltonian $`\widehat{𝐇}_0`$ controlled by an external field $`ϵ(t)`$,

$$\widehat{𝐇}(t)=\widehat{𝐇}_0\widehat{𝝁}ϵ(t).$$
(1)
where $`\widehat{𝝁}`$ is a system operator. In the molecular system, $`\widehat{𝝁}`$ is the transition dipole operator and $`ϵ(t)`$ describes a shaped short light pulse. For simplicity, the field $`ϵ(t)`$ is assumed real but the generalization to an electromagnetic field with two independently controlled polarizations gerber2001  is straightforward.

The algorithm is represented in the target time $`T`$ by a unitary transformation $`\widehat{𝐔}(T)`$ generated from the Hamiltonian Eq. (1);

$$\frac{\widehat{𝐔}(t)}{t}=\frac{i}{\mathrm{}}\widehat{𝐇}(t)\widehat{𝐔}(t),$$
(2)
with the initial condition $`\widehat{𝐔}(t=0)=𝟙`$, where $`𝟙`$ denotes the identity operator. The objective is to obtain the optimal driving field $`ϵ(t)`$ that induces a given unitary transformation $`\widehat{𝐎}`$ at $`t=T`$, i.e., $`\widehat{𝐔}(T)=e^{i\varphi }\widehat{𝐎}`$. $`\varphi `$ denotes a physically irrelevant global phase that is related to the energy origin and is therefore uncontrollable by the field.

The task is a typical inversion problem which can be solved by employing a variational procedure maximizing the projection of the generated operator on the target operator,

$$|\tau |=|\mathrm{Tr}\{\widehat{𝐎}^{}\widehat{𝐔}(T)\}|,$$
(3)
where the projection $`(\widehat{𝐀}\widehat{𝐁})`$ is defined by $`\mathrm{Tr}\{\widehat{𝐀}^{}\widehat{𝐁}\}`$. The functional $`\tau `$ is a complex number inside the circle $`|\tau |N_H`$, where $`N_H`$ is the dimension of the Hilbert space of the system. Equality is reached only when the argument of the trace is $`e^{i\varphi }𝟙`$ and then a maximum of $`|\tau |`$ is equivalent to $`\widehat{𝐔}(T)=e^{i\varphi }\widehat{𝐎}`$. The optimal solution is then found by maximizing the functional Eq. (3), with respect to the control field $`ϵ(t)`$. Since a direct algorithm to maximize $`|\tau |`$ was not found, a working alternative is used based on formulating the problem as the optimization of $`\mathrm{Re}[\tau ]`$, or of $`\mathrm{Im}[\tau ]`$, or of a linear combination of both. For simplicity, the optimization of the real part represented by the functional $`J=\mathrm{Re}[\mathrm{Tr}\{\widehat{𝐎}^{}\widehat{𝐔}(T)\}]`$ is considered. Two constraints are introduced pierce88 ; kosloff89 , the first restricts the dynamics to obey the Schrödinger equation Eq. (2), and the second restricts the total field energy. Using Lagrange multipliers a modified functional is obtained,

$$\overline{J}=\mathrm{Re}\left[\tau _0^T\mathrm{Tr}\left\{\left(\frac{\widehat{𝐔}(t)}{t}+\frac{i}{\mathrm{}}\widehat{𝐇}(t)\widehat{𝐔}(t)\right)\widehat{𝐁}(t)\right\}𝑑t\right]\lambda _0^T\frac{1}{s(t)}|ϵ(t)|^2𝑑t,$$
(4)
where $`\widehat{𝐁}(t)`$ is an operator Lagrange multiplier, $`\lambda `$ is a scalar Lagrange multiplier and $`s(t)`$ is a shape function which turns the pulse on and off sundermann99 . The use of more elaborate constraints and choices of $`\lambda `$, allows a higher degree of control on the shape of the optimal pulse hornung01 ; hornung02 .

Applying the calculus of variations, $`\delta \overline{J}=0`$ with respect to $`\widehat{𝐁}`$, $`\widehat{𝐔}`$, and $`ϵ`$, a set of equations is obtained: a) The Schrödinger equation Eq. (2) with the initial condition $`\widehat{𝐔}(t=0)=𝟙`$ for $`\widehat{𝐔}`$; b) The inverse Schrödinger equation

$$\frac{\widehat{𝐁}(t)}{t}=\frac{i}{\mathrm{}}\widehat{𝐁}(t)\widehat{𝐇}(t)$$
(5)
with the condition $`\widehat{𝐁}(t=T)=\widehat{𝐎}^{}`$ for $`\widehat{𝐁}`$; c) The field equation:

$$ϵ(t)=\frac{s(t)}{2\lambda \mathrm{}}\mathrm{Im}[\mathrm{Tr}\{\widehat{𝐁}(t)\widehat{𝝁}\widehat{𝐔}(t)\}],$$
(6)
Eq. (2) and (5) represent two counter currents with information from the initial condition and the target unitary transformation respectively. The equations are solved iteratively, the Krotov method tannor92 , similar to the methods described in Ref. zhu98 , was found to be the most efficient. The input is a “guess” field, $`ϵ^{(0)}(t)`$, so that in the $`k`$ iteration ($`k=1,2\mathrm{}`$): (i) $`\widehat{𝐁}^{(k1)}(t)`$ is propagated backwards from $`t=T`$ to $`t=0`$ using Eq. (5) and $`ϵ^{(k1)}`$; and (ii) $`\widehat{𝐔}^{(k)}(t)`$ is propagated forward using Eq. (2) and $`ϵ^{(k)}`$ is evaluated using

$$ϵ^{(k)}(t)=\frac{s(t)}{2\lambda \mathrm{}}\mathrm{Im}[\mathrm{Tr}\{\widehat{𝐁}^{(k1)}(t)\widehat{𝝁}\widehat{𝐔}^{(k)}(t)\}].$$
(7)
The procedure is repeated until the desired convergence has been reached. The hard numerical task is the propagation of the operators $`\widehat{𝐔}`$ and $`\widehat{𝐁}`$ with the time dependent Hamiltonian for which a second order Newton polynomial integrator kosloff94  was used.

A direct use of Eq. (7) in the algorithm leads to saturation. This is because the constraint related with the field energy becomes more important than the original objective. A remedy is to interpret the right hand side of Eq.(7), denoted by $`\mathrm{\Delta }ϵ^{(k)}`$, as a correction of the field in the previous interaction bartana97 . The field after the $`k`$ iteration is then given by $`ϵ^{(k)}(t)=ϵ^{(k1)}(t)+\mathrm{\Delta }ϵ^{(k)}(t)`$. The OCT procedure was first applied to a 1 qubit operation the Hadamard rotation (Eq. 10). In a similar fashion, 2-qubit operation such as the “controlled not” were obtained. The optimal fields generating unitary operations with eight (3-qubits) were also attempted.

It was found that the number of iterations required to converge the results increased approximately exponentially (or factorially) with the number of states in the problem. A possible reason for this scaling is that not only a specific state to state transition has to be forced but this has to be carried out without disturbing the other state to state transitions in $`\widehat{𝐔}`$. As a result the scaling becomes $`O(n!)`$ which is consistent with numerical experience.

In a molecular environment, obtaining the optimal field to carry out an algorithm is more involved. The register levels which are used to write the input and output are only part of a much larger manifold of molecular energy levels. Considering the advances in pulse shaping techniques in the visible region of the spectrum the transitions of choice are electronic. For such a molecular construction, imposing a unitary transformation on the total set of levels which are addressed by the field is too restrictive. Relying on the experience that the convergence is close to factorial, an extremely large number of iterations would be required to converge.

The strategy is therefore to restrict the target objective to only the states used directly as registers keeping the other states in the system as passive observers. The reduced objective is obtained by changing the previous expressions $`\mathrm{Tr}\{\widehat{𝐀}\widehat{𝐁}\}`$ to $`_ii|_R\widehat{𝐀}\widehat{𝐁}|i_R`$, where $`\{|i_R\}`$ is a basis of the subspace of registers. The subindex $`R`$ is used to denote the operators in that subspace, for example $`\widehat{𝐎}_R`$ denotes the target unitary transformation. The substitution of the restricted condition instead of the trace in Eq. (3) keeps the maximum condition when the unitary transformation $`\widehat{𝐔}_R(T)`$ is equal to $`\widehat{𝐎}_R`$. The maximum value becomes equal to the dimension of the subspace $`N_R`$. In the condition for $`\widehat{𝐁}`$ at time $`T`$, a particular dependence must be specified for all the levels. For simplicity the identity in the passive subspace has been chosen.

The approach is illustrated using two models. For the first model, the implementation of a unitary transformation in one qubit while minimizing the population transfer to other levels in the molecule is studied. The model consists of a molecule with two electronic surfaces described by the Hamiltonian,

$$\widehat{𝐇}=\left(\begin{array}{cc}\widehat{𝐇}_g& \widehat{𝝁}ϵ(t)\\ \widehat{𝝁}ϵ(t)& \widehat{𝐇}_e\end{array}\right),$$
(8)
where $`\widehat{𝐇}_g`$ and $`\widehat{𝐇}_e`$ are the ground and excited surface Hamiltonians. The electronic surfaces are coupled by the transition dipole operator $`\widehat{𝝁}`$, controlled by the shaped field $`ϵ(t)`$. The two first levels of the ground electronic surface are chosen as the registers representing the qubit. The model includes $`15`$ rovibroinc levels in the ground electronic state and $`5`$ in the excited state described by the Hamiltonians,

$$\widehat{𝐇}_g=\underset{i=1}{\overset{15}{}}E_{gi}|g_ig_i|;\widehat{𝐇}_e=\underset{j=1}{\overset{5}{}}E_{ej}|e_je_j|.$$
(9)
shown in Fig. 1. Next, a transition dipole operator with equal coupling strength between levels was chosen, $`\widehat{𝝁}=\mu _0_{i=1}^{15}_{j=1}^5|g_ie_j|+\mathrm{H}.\mathrm{c}.`$, where $`\mathrm{H}.\mathrm{c}.`$ denotes the Hermitian conjugate.

In Fig. 1 the optimized field is shown when the target unitary transformation is a Hadamard rotation given by,

$$\widehat{𝐎}_R^H=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),$$
(10)
restricted to the $`2\times 2`$ qubit subspace. The shape function was chosen as $`s(t)=\mathrm{sin}^2(2\pi t/T)`$ and the initial guess was $`ϵ^{(0)}(t)=s(t)\mathrm{cos}(\mathrm{\Omega }t)`$, where $`\mathrm{\Omega }`$ is the frequency of the 00-line of the electronic transition. Fig. 1 shows that the dominant frequencies are the ones related to the known vibronic molecular transitions. In the time domain the field is split into a symmetric sequence of sub pulses. A phase relation correlating the dominant frequencies is observed in a Wigner plot (not shown). These phase relations guarantee that no population is lost to the excited levels.

In the second model the two electronic surfaces, Eq. (8), have two orthogonal vibrational modes, denoted as $`\alpha `$ and $`\beta `$. The first two levels of each vibrational mode in the ground electronic surface are chosen as the the physical implementation of the two qubits. The goal is to induce an operation that involves entanglement between them. In this model each vibrational mode of the ground surface is coupled by the field to the corresponding mode in the excited electronic surface. The modes in the excited surface are coupled by a static term modeling Duschinsky rotation. This last term which is not controlled by the field can generate entanglement between the qubits. The model consists of two levels for each vibrational normal mode denoted by $`\alpha `$ and $`\beta `$. The electronic surface Hamiltonians are $`\widehat{𝐇}_g=\widehat{𝐇}_{g\alpha }𝟙_\beta +𝟙_\alpha \widehat{𝐇}_{𝕘\beta }`$ and $`\widehat{𝐇}_e=\widehat{𝐇}_{e\alpha }𝟙_\beta +𝟙_\alpha \widehat{𝐇}_{𝕖\beta }+\widehat{𝐕}_{\alpha \beta }`$ with,

$`\widehat{𝐇}_{g\nu }`$ $`=`$ $`\left(E_{g\nu _0}|g_0_\nu g_0|_\nu +E_{g\nu _1}|g_1_\nu g_1|_\nu \right),`$
$`\widehat{𝐇}_{e\nu }`$ $`=`$ $`\left(E_{e\nu _0}|e_0_\nu e_0|_\nu +E_{e\nu _1}|e_1_\nu e_1|_\nu \right),`$
$`\widehat{𝐕}_{\alpha \beta }`$ $`=`$ $`\delta _{\alpha \beta }\left(|e_0_\alpha |e_1_\beta e_1|_\alpha e_0|_\beta +|e_1_\alpha |e_0_\beta e_0|_\alpha e_1|_\beta \right),`$ (11)
with $`\widehat{𝐕}_{\alpha \beta }`$ the Duschinsky term that couples the vibrational modes in the excited surface. A schematic representation of the levels is given in Fig. 2. The transition dipole operator is chosen as $`\widehat{𝝁}=\widehat{𝝁_𝜶}𝟙_\beta +𝟙_\alpha \widehat{𝝁_𝜷}`$, with $`\widehat{𝝁}_\nu =\mu _{0\nu }(|e_0_\nu +|e_1_\nu )(g_0|_\nu +g_1|_\nu )+\mathrm{H}.\mathrm{c}.`$ Operators in the combined Hilbert space of this system are represented by $`16\times 16`$ matrices.

The target unitary transformation is a two qubits quantum Fourier transform weinstein01 ,

$$\widehat{𝐎}_R^{FT}=\frac{1}{2}\left(\begin{array}{cccc}1& 1& 1& 1\\ 1& i& 1& i\\ 1& 1& 1& 1\\ 1& i& 1& i\end{array}\right)$$
(12)
where the operator is represented in the basis $`|g_i_\alpha |g_j_\beta `$. Notice that operators in the subspace of interest are represented by $`4\times 4`$ matrices. The spectra of the field resulting from the optimization is shown in Fig. 2. The frequencies are confined in the region of the vibronic transitions. The structure of the spectra is much more complex than in the previous model. An obvious reason is the larger dimension of the subspace of interest. In addition the complicated mechanism required to create entanglement is indirect. During the pulse significant population is transfered to the excited electronic states. The resulting more complex field is reflected by also in the Wigner distribution, where the field shows a very high degree of correlation between time and frequency.

In conclusion, Optimal Control Theory has been generalized to obtain the driving field that generates any target unitary transformation. In principle the scheme allows the implementation of any quantum gate using a single step. Compared to the implementation of the gate using a sequence of known simple pulses, the loss by the possible complexity of the shaped pulse can be more than compensated by the faster implementation. Convergence is efficient for a small number of levels. Besides, with the right choice of the guess field, the method can obtain the driving field for more complex systems where a large number of levels are involved and the pulses for simple operations are unknown, as in molecular systems. The advantage of such systems is the short time in which these algorithms can be executed. Using Li2 as an example, the Hadamard rotation can be executed in $``$ 1.5 psec. The 2-qubit Fourier transform could be executed in a molecule for example OCS in approximately $`6`$ psec. These fast timescales give hope that the quantum computation can be carried out before decoherence processes take place.

J.P.P. acknowledges the Golda Maier Fund of the Hebrew University. This work was supported by the Spanish MCT BMF2001-3349 and the Israel Science Foundation. The Fritz Haber Center is supported by the Minerva Gesellschaft für die Forschung, GmbH München, Germany. The authors thank David Tannor and Zohar Amity for helpful discussions.