text:
# Bose-Einstein condensation in a simple Microtrap

Website: \]http://www.bec.uni-hd.de

## Abstract

A Bose-Einstein condensate is created in a simple and robust miniature Ioffe-Pritchard trap, the so-called Z trap. This trap follows from the mere combination of a Z-shaped current carrying wire and a homogeneous bias field. The experimental procedure allows condensation of typically $`3\times 10^5`$ 87Rb atoms in the $`|F=2,m_F=2`$ state close to any mirroring surface, irrespective of its structure, thus it is ideally suited as a source for cold atom physics near surfaces.

Bose-Einstein condensation, Magnetic Traps, Cold Atom Physics, Atom Chip

Bose-Einstein condensates are at the heart of research in quantum optics. Since the first realization Anderson et al. (1995); Davis et al. (1995), various methods of creating a condensate have been shown, ranging from large-scale magnetic traps to miniature surface traps Ott et al. (2001); Hänsel et al. (2001) as well as optical dipole traps Barrett et al. (2001). Many fascinating properties of this form of matter have been studied Ketterle et al. (1999); K. Ziemelis, Ed. (2002); GSU .

Manipulation of ultra cold atoms using surface mounted micro-structures, so-called atom chips Folman et al. (2000, 2002), promises accuracy and versatility in manipulating atoms at the quantum level, and perhaps even the implementation of quantum information processing Bouwmeester et al. (2000). For a comprehensive review on surface mounted micro-traps, see Ref. Folman et al. (2002). One of the main experimental challenges is to create a simple source of ultra-cold (ground state) atoms for the experiments, irrespective of the structures on the surface.

In this Rapid Communication we report on the creation of a sizable ($`3\times 10^5`$ atoms) Bose-Einstein condensate in a Z-wire trap, which yields a simple, easy to handle and robust small-scale Ioffe-Pritchard trap Haase et al. (2001). By virtue of the trapping wire being part of the atom chip mounting, both the experimental procedure leading to a condensate, and its alignment to the atom chip are independent of the actual structure of the atom chip, or any other reflecting surface. Thus it provides the means to combine the versatile atom chip assembly and the stringent UHV demands for Bose-Einstein condensation with a high number of atoms.

The most basic microtrap uses the magnetic field minimum created by superposing the magnetic field of a current in a straight wire with a homogeneous magnetic bias field perpendicular to it Frisch and Segre (1933); Folman et al. (2002). The two fields add up to a 2D quadrupole field along the wire, which can guide atoms. The trap depth is given by the homogeneous bias field, the field gradient is inversely proportional to the wire current. Such a geometry lends itself to miniaturization of the wire size using micro-fabrication techniques. In this way traps can be made even steeper, while maintaining mechanical stability and robustness. A typical example is our atom chip Folman et al. (2000), where the 1 - 200 $`\mu `$m wide wires are micro fabricated in a 2 - 5 $`\mu `$m thick gold layer on a silicon substrate.

By bending the ends of the wire one creates slightly more complex structures that provide a three dimensional confinement Denschlag et al. (1999); Haase et al. (2001); Reichel et al. (1999).
a) The U-trap: by bending a wire in a U-shape, the two fields from the leads close the guide along the base. The result is a three dimensional quadrupole field, with a trap minimum $`B_0=0`$.
b) The Z-trap: if the leads are pointing in opposite directions (Z-shaped wire) a field component parallel to the base remains, giving a Ioffe-Pritchard type trap with a non-zero minimum $`B_0>0`$.

A common property of these traps is the scaling of the strong linear confinement in the transverse direction, inversly proportional to the wire current $`I_W`$ and the distance $`z_0`$ between the wire center and the trap minimum.

$$z_0\frac{I_W}{B_x},B_{}^{}=\frac{B}{x,z}\frac{B_x^2}{I_W}.$$
(1)
Here $`B_x`$ is the homogeneous bias field perpendicular to the wire. In the case of a Z-trap, the additional longitudinal field arising from the leads turns the bottom of the quadrupole trap into a harmonic potential. The transverse angular trap frequencies $`\omega _{x,z}`$ in this potential scale as:

$$\omega _{x,z}\frac{B_{}^{}}{\sqrt{B_0}}.$$
(2)
Such wire traps provide a linear field gradient given by Eq. (1). This enables a large flexibility to first trap a large number of atoms and second to efficiently compress the trap to small volumes and high trap frequencies for effective evaporative cooling. For instance, using a wire current of $`I_W=50`$ A and properly tuned bias fields, thermal atoms can be trapped up to 3 mm from the wire center. By only varying the bias fields, the trap can be compressed to field gradients in excess of 400 G/cm at a trap distance of 1.5 mm from the wire center, reaching trap frequencies of $`\omega _{x,z}>600`$ Hz. The large distance allows a micro-fabricated structure to be placed between the wire and the trap. This is the principle advantage of the experimental setup. An additional benefit of such small-scale circuits is that their low self-inductance enables rapid changing and switching of the potentials.

As presented in Fig. 1, the experimental apparatus mainly consists of a vertically aligned double magneto-optical trap (MOT) setup with a mirror MOT Lee et al. (1996); Pfau and Mlynek (1996); Reichel et al. (1999); Folman et al. (2000) in the upper ultra-high vacuum chamber. The atom chip gold layer serves as the mirror Folman et al. (2000). The magnetic quadrupole field for the MOT is supplied by external coils or by the field of the U-trap (U-MOT). For the mirror MOT to work, the magnetic axis of the quadrupole must be tilted by $`45^{}`$ with respect to the mirror. To maintain the optical access to the trap the coils are wound in a conical shape. Two counter-propagating laser beams overlap the quadrupole axis, being reflected on the mirror at $`45^{}`$. Two additional counter-propagating laser beams run horizontally, parallel to the mirror and perpendicular to the plane of the $`45^{}`$-reflected beams. In the case of the U-MOT, the $`45^{}`$ angle of the field axis follows naturally from the geometry of the wire and the bias field.

The fields for the Z-trap and the U-MOT are derived from currents through the adequate ports of a monolithic H-shaped copper structure, see the inset of Fig. 1. The central bar of the H-structure forms the base wire of both the U- and the Z-trap. The base wire length for the U-MOT is 14.5 mm (outer ports), while for the Z-trap it is 7.25 mm (inner ports). The device is tightly fitted in an insulating ceramic, allowing any structure to be mounted directly against it. The distance between the center of the base wire and the surface of an atom chip is 1.2 mm. The $`1.2\times 1`$ mm2 cross-section of the base wire withstands currents of up to 50 A for over a minute without significant heating. Thus it enables trapping of large atom numbers both near and far away from the chip-surface, while preserving the ultra high vacuum in the upper chamber.

The base of the trap mounting is made from stainless steel tubes welded to a CF63 UHV-flange, allowing a coolant flow through the mounting. The flange also carries four high-power current feedthroughs for operating the H-structure and a 35-pin connector to access the atom chip. The ceramic block containing the H-structure is fastened to the stainless steel part of the mounting by a copper rod, acting as a heat sink. The mounting is hung from its flange in the upper, ultra high vacuum, chamber of the setup, which is machined from a single piece and is fitted with optical quality windows. The distance between the atom trap center and the closest window is $`40`$ mm.

A regular vapour cell MOT in the lower part of the setup traps $`10^9`$ 87Rb atoms from the background gas obtained from a heated Rb-dispenser ($`P10^8`$ mbar). A continuous laser beam, the push beam, transfers the atoms through a differential pumping stage and over a 40 cm vertical distance from the vapour cell MOT chamber($`P10^8`$ mbar) to the UHV chamber($`P10^{11}`$ mbar) Bartenstein (2001); Wohlleben et al. (2001).

The experimental cycle starts by loading the mirror-MOT by the push beam for 20 s, resulting in a cloud of $`10^9`$ atoms located $`3`$ mm below the atom chip. Next, the quadrupole field from the mirror-MOT is replaced in 300 ms by that of the U-MOT. Using a wire current of $`I_W=25.4`$ A and bias field of $`B_x=6.5`$ G, a field gradient of 14 G/cm is reached. The brief transfer to the U-MOT compresses the cloud and moves it to 2 mm from the surface. During this process we observe no atom losses.

After establishing the U-MOT, all magnetic fields are switched off for 15 ms to allow optical molasses cooling to $`50\mu `$K. It should be noted that the atom chip surface is not a perfect mirror, due to the structures on the chip. However, the diffraction from these edges does not harm in any aspect the functionality of the mirror-MOT, nor does it inhibit molasses cooling. A 200 $`\mu `$s optical pumping pulse to the $`|F=2,m_F=2`$ state prepares the atoms for the transfer to the Z-trap. Optical pumping increases the number of trapped atoms by a factor of 3. With careful matching of the trap parameters, more than $`10^8`$ atoms can be loaded into the magnetic trap.

The trap frequencies for the initial Z-trap, with a Z-current of 49.7 A and a bias field of $`B_x=25.4`$ G (dotted line in Fig. 2), are $`\omega _y=2\pi \times 18`$ Hz and $`\omega _{x,z}=2\pi \times 51`$ Hz. Immediately after loading, the offset field of the Z-trap is reduced by the addition of a bias field $`B_y`$, parallel to the base, opposing the fields of the leads. This compresses the trap transversely and deepens it. After 100 ms this bias field reaches its final value of $`B_y=5.5`$ G (dashed line in Fig. 2). The trap is further compressed in all three directions by linearly increasing the bias field $`B_x`$ from 25.4 G to 58.3 G, which also moves the trap center from initially 1.8 mm to a distance of $`300\mu `$m from the chip surface (solid line in Fig. 2). During the 19 s compression stage, forced evaporative cooling is applied by RF-radiation from a coil-antenna outside the vacuum chamber. The frequency is linearly ramped down from 19 MHz to typically 600 kHz. The simultaneous processes of compression and cooling avoid crushing the initially large atom cloud into the chip surface. For the final trapping potential, the angular oscillation frequencies are $`\omega _{x,z}=2\pi \times 600`$ Hz and $`\omega _y=2\pi \times 70`$ Hz (solid line in Fig. 2). The compression increases the transverse field gradient by more than a factor of 5, from 80 G/cm to 440 G/cm. The wire current is kept constant throughout the procedure.

The atom cloud is studied by absorption imaging with a weak ($`120\mu `$W/cm2) circularly polarised laser beam, resonant to the $`\mathrm{F}=2\mathrm{F}^{}=3`$ transition. The imaging system uses a 8-Bit CCD camera and achromatic optics with a resolution of $`3\mu `$m. The direction of view is parallel to the plane of the atom chip and perpendicular to the central bar of the H-structure. The laser beam grazes the chip, enabling detailed pictures down to very short distances to the chip. Care has to be taken in aligning the absorption laser beam, because slight misalignment can cause severe diffraction stripes to appear in the images.

A helpful feature of the setup is that the distance of the cloud to the chip surface may be inferred from slightly tilting the imaging laser beam; one observes two images of the condensate, where one results from a partially reflected absorption beam (see Fig. 3).

Modeling the trap geometry with a straightforward approach using infinitely long leads and infinitely thin wires, neglecting the finite size of the structure, gives a qualitative estimate for the parameters of the trap. Detailed simulations of the current distribution in the Z-structure, employing finite element algorithms, reveal that the effective current path length $`L_{\mathrm{eff}}`$ in the Z-structure is shorter than the geometrical value, $`L_{\mathrm{geom}}`$, defined by the separation of the leads along the base wire (see inset in Fig. 2). Using an effective baselength of $`L_{\mathrm{Eff}}6.7`$ mm, instead of $`L_{\mathrm{Geom}}=7.25`$ mm, in the model with infinitely thin wires, we reach good quantitative agreement with the experimental results.

The procedure described above results in a Bose-Einstein condensate of approx. $`3\times 10^5`$ atoms, at a density of $`9\times 10^{13}`$ cm-3 and with a transition temperature of $`600`$ nK, see Fig. 4. The temperature difference of at least 8 orders of magnitude between the condensate and the uncooled atom chip can be separated by less than 250 $`\mu `$m. The condensate survives for over a second when using RF-shielding at the final frequency of the cooling ramp.

Since in this work the atom chip is only used as a reflecting surface for the mirror-MOT, we have shown that a single current carrying wire in a homogeneous bias field suffices to create and store a Bose-Einstein condensate near a hot reflective surface, irrespective of its structure. This makes the setup very well suited for studying the interaction between degenerate atom clouds and surfaces. Since the large volume of the macroscopic Z-trap allows trapping, storing and cooling of large atom numbers, the technique can easily be used to load pre-cooled ensembles into small-scale traps of a different nature, such as miniature optical traps, traps made of permanent magnets, electric traps, etc., which may have shallow depths or small trapping volumes, thus lacking the ability to store sufficient many atoms for compression and evaporative cooling.

We especially like to thank Toni Schönherr for the fine mechanical work. This work was supported by the EU project ACQUIRE, Austrian Science Fund (FWF), and the Deutsche Forschungsgemeinschaft (Grant Schm1599/2-1). L.F. acknowledges the support of the Alexander von Humboldt Foundation.


format: text

text:
# Quantum Computing of Quantum Chaos in the Kicked Rotator Model

## I Introduction

It is only recently that it was realized that quantum mechanics can be used to process information in fundamentally new ways. In particular, Feynman  emphasized that the massive parallelism due to the superposition principle may allow to simulate efficiently some problems intractable on classical computers, the most obvious being many-body quantum systems. Since that time, a model of quantum computer has been set up, viewed as an ensemble of n qubits, i.e. two-level systems, with a Hilbert space of dimension $`2^n`$ (see reviews ). Computation is performed through unitary transformations applied to the quantum wave functions of this many-body system. In fact, it has been shown that any unitary transformation on this $`2^n`$ dimensional space can be written in terms of a set of universal gates, for example one- and two-qubit transformations. Also, important quantum algorithms have been developed, such as Shor’s algorithm for factoring large numbers  which is exponentially faster than any known classical method, and Grover’s algorithm to search a database , where the gain is polynomial.

Motivated by these developments, many experimental implementations for actual realization of such a quantum computer were proposed (see  and references therein). Recent results include for example the NMR implementation of factorization algorithm with seven qubits made from nuclear spins in a molecule , and the simulation of the quantum baker map . Thus small quantum computers with a few qubits are already available experimentally, and systems of larger size can be envisioned at relatively short term.

Still, algorithms such as the one of Shor require large number of qubits and the use of many gates. It is therefore important to develop algorithms which need a smaller number of qubits and gates and still can yield interesting quantities. In particular, algorithms enabling to simulate quantum mechanical systems, as originally envisioned by Feynman, can be implemented relatively easily and solve problems inaccessible to classical computers with less expanse in number of qubits and gates. Several such algorithms have been developed for various systems, including many-body Hamiltonians  or spin lattices . An especially interesting class of systems corresponds to chaotic quantum maps. Such systems can have a very complex dynamics while their Hamiltonians keep a relatively simple form. Algorithms for fast simulation on a quantum computer were built for the quantum baker map , the kicked rotator , and the sawtooth map . We note that recently the quantum baker map was implemented on a NMR quantum computer . The kicked rotator is an especially rich and generic system, which has been a cornerstone for the study of quantum chaos . In the classical limit it reduces to the Chirikov standard map which has been also extensively studied in the field of classical chaos . Implementation of this model can be done on a small quantum computer with a few tens of qubits, and classical supercomputers will be outperformed with a few hundreds of qubits. Still, real quantum computers will not be free of imperfections and errors, and this will affect the results of the computation. It is therefore important to understand the effects of different sources of errors on the results of such an algorithm. For example, first numerical simulations of the quantum computation of this model  have shown that errors affect in a different way the various physical quantities characterizing the model, and that for some of them the effect of errors can be exponentially strong.

In this paper, after presenting in more detail the physics of the kicked rotator, we study the effects of errors on several physical quantities. We focus on random unitary errors, which may arise when imperfect gates are applied, and study first how global quantities such as second moment or fidelity are affected by errors. Our results confirm and extend those obtained in  showing a marked contrast in the behavior of these two quantities in presence of errors. We also investigate how well the whole wave function is reproduced by an imperfect quantum computer. A particularly interesting way to display wave functions is to express them through phase space distributions, such as Wigner and Husimi functions. These distributions display the same information as the wave functions, but in a form which allows direct comparisons between classical and quantum dynamics, a property especially interesting to probe the classical limit of quantum mechanics. They have been extensively used in many fields, and recently a method has been devised  to measure such distribution for quantum simulations on quantum computers. The effects of errors on Wigner and Husimi functions will be investigated in details, showing how imperfections affect the different parts of phase space, and discussing how information can be retrieved through quantum measurement. A separate section is focused on how a localized distribution may escape from an island of integrability, showing an especially large effect of quantum errors on a quantity which is directly relevant to quantum tunneling.

## II The kicked rotator

The classical kicked rotator is described by the Chirikov standard map :

$$\overline{n}=n+k\mathrm{sin}\theta ;\overline{\theta }=\theta +T\overline{n}$$
(1)
where $`(n,\theta )`$ is the pair of conjugated momentum (action) and angle variables, and the bars denote the resulting variables after one iteration of the map. It describes a free angle rotation and a kick in momentum. This area-preserving map has been extensively studied during the past decades, and has been applied to problems such as particle confinement in magnetic traps, beam dynamics in accelerators, comet trajectories and many others .

The dynamics of this map takes place on a cylinder (periodicity in $`\theta `$), and is controlled by a single parameter $`K=kT`$. For $`K=0`$ the system is integrable and all trajectories lie on one-dimensional tori (lines $`n=`$constant). For $`K>0`$, the system undergoes a transition to chaos, which follows the Kolmogorov-Arnold-Moser (KAM) theorem. Periodic orbits corresponding to rational frequencies are transformed into chains of integrable islands mixed with chaotic region. On the contrary, tori with irrational frequencies are deformed but survive, forming invariant curves which separate zones in phase space. As $`K`$ is increased, these surviving tori become Cantor sets (cantori) and disappear. The most robust torus corresponds to the golden number $`(1+\sqrt{5})/2`$, and disappears for $`K=K_g0.9716\mathrm{}`$. Thus, for $`K>K_g`$ global chaos sets in, with appearance of an extended chaotic region in phase space, and with dynamics characterized by a positive Kolmogorov-Sinai entropy $`h\mathrm{ln}(K/2)>0`$ (for $`K6`$). In this régime, a typical trajectory shows diffusive growth of momentum which statistically can be described by the Fokker-Planck equation, with diffusion rate $`D=n^2/tk^2/2`$ where $`t`$ is measured in number of iterations (kicks) . For lower values of $`K`$, the phase space displays a complex hierarchical structure with integrable islands surrounded by chaotic zones at smaller and smaller scales.

The map (1) is periodic in $`n`$ with period $`2\pi /T`$, so the phase space structures repeat themselves on each cell of size $`2\pi /T`$. Such a cell is shown on Fig.1 for $`K=K_g`$, displaying the complex hierarchical structures which appear in the phase space.

The quantization of (1) yields a Hamiltonian which after integration over one period gives a unitary evolution operator acting on the wave function $`\psi `$

$`\overline{\psi }=\widehat{U}\psi =e^{ik\mathrm{cos}\widehat{\theta }}e^{iT\widehat{n}^2/2}\psi ,`$ (2)
where $`\widehat{n}=i/\theta `$, $`\mathrm{}=1`$, and $`\psi (\theta +2\pi )=\psi (\theta )`$. The quantum dynamics depends on two parameters $`k`$ and $`T`$ (instead of the single parameter $`K=kT`$ for the classical one). The classical limit corresponds to $`k\mathrm{}`$, $`T0`$ while keeping $`K=kT=`$constant . In a sense, $`T`$ plays the role of an effective $`\mathrm{}`$.

Depending on the values of theses parameters, the system follows different régimes, from regularity to quantum chaos. Due to this variety of behaviors, the quantum kicked rotator has been intensively studied (see  and Refs. therein). Indeed, most of the phenomena characteristic of quantum chaos are present, such as quantum ergodicity, Random Matrix Theory statistics, chaos assisted tunneling, and others. In particular, for $`K>K_g`$, the phenomenon of dynamical localization appears. Although in this régime a typical classical trajectory diffuses in momentum, the eigenstates $`\chi _m(n)`$ of the operator $`\widehat{U}`$ in momentum space are exponentially localized for typical values of $`k`$ and $`T`$. Their envelopes obey the law $`\chi _m(n)\mathrm{exp}(|nm|/l)/\sqrt{l}`$ where $`m`$ marks the center of the eigenstate and $`l`$ is the localization length. For $`kK1`$ this length is determined by the classical diffusion rate $`l=D/2k^2/4`$ . This phenomenon has close relationship with the Anderson localization of electrons in disordered solids , and investigation of the kicked rotator gives information on this important solid-state problem still under intensive investigation nowadays. The quantum kicked rotator describes also the properties of microwave ionization of Rydberg atoms . It has been realized experimentally with cold atoms, and the effects of dynamical localization, external noise and decoherence have been studied experimentally .

For numerical studies of the quantum evolution (2) it is convenient to choose the case of quantum resonance with $`T/(4\pi )=M/N`$ where $`M,N`$ are integers . In this way the quantum dynamics takes place on a torus with $`N`$ levels. For $`lN`$ the eigenstates of evolution operator become ergodic and the level spacing statistics is described by random matrix theory .

The algorithm for the quantum simulation of the kicked rotator was presented in . The evolution (2) consists of the product of two unitary operators, which are diagonal in the angle and momentum bases respectively. The most efficient classical algorithm available consists in changing back and forth between the angle and momentum representation by Fast Fourier Transforms (FFT). The operator which is diagonal in the basis is then implemented by direct multiplication of the coefficients of the wave function. In total, one iteration of (2) on a Hilbert space of dimension $`N=2^{n_q}`$ requires $`O(N\text{log}N)`$ classical operations, the limiting steps being the FFT. The quantum algorithm follows the classical one, and speeds up all parts of it to obtain exponential increase of computation rate. First an initial distribution is built, in a polynomial number of operations (in $`n_q`$). Various initial wave functions can be built in such a way. In the following, we will use as initial state $`|\mathrm{\Psi }_0`$ a wave function localized at a precise value of momentum $`n`$, which can be built in $`n_q`$ single-qubit rotations starting from the ground state. The general state of the system can be written as $`_{n=0}^{N1}a_n|n`$, where $`a_n`$ are the amplitudes of the wave function on the $`|n`$ basis state. Then the first unitary operator is applied. In the $`n`$ representation it is diagonal and can be written $`\mathrm{exp}(iTn^2/2)`$. This operator can be implemented efficiently by using the binary decomposition of $`n`$: if $`n=_{j=0}^{n_q1}\alpha _j2^j`$, then $`n^2=_{j_1,j_2}\alpha _{j_1}\alpha _{j_2}2^{j_1+j_2}`$. Therefore $`\mathrm{exp}(iTn^2/2)=\mathrm{\Pi }_{j_1,j_2}\mathrm{exp}(iT\alpha _{j_1}\alpha _{j_2}2^{j_1+j_21})`$ with $`\alpha _{j_{1,2}}=0`$ or $`1`$. Thus one needs to implement the two-qubit gate applied to each qubit pair $`(j_1,j_2)`$ which keeps the states $`|00,|01,|10`$ unchanged while $`|11`$ is transformed to $`\mathrm{exp}(iT2^{j_1+j_21})|11`$. $`O(n_q^2)`$ applications of this gate are sufficient to simulate $`\mathrm{exp}(iTn^2/2)`$.

Then a quantum Fourier transform (QFT) (see e.g. ) is performed to shift from $`n`$ to $`\theta `$ representation, yielding $`_{i=0}^{N1}b_i|\theta _i`$. This transformation needs only $`O(n_q^2)`$ one and two-qubit gates, and yields the wave function in $`\theta `$ representation. In this representation, the second operator $`\mathrm{exp}(ik\mathrm{cos}\widehat{\theta })`$ is diagonal. Direct (sequential) multiplication by $`\mathrm{exp}(ik\mathrm{cos}\theta _i)`$ for each $`\theta _i`$ will require exponentially many operations, so a parallel way to apply this operator has to be devised. In , it was proposed to use supplementary registers on which the values of $`\mathrm{cos}(\theta _i)`$ will be computed in parallel. The procedure transforms $`_{i=0}^{N1}b_i|\theta _i|0`$ into $`_{i=0}^{N1}b_i|\theta _i|\mathrm{cos}\theta _i`$ with $`\mathrm{cos}(\theta _i)`$ computed up to a fixed precision using a recursive method based on Moivre’s formula . This is actually the slowest step of the algorithm, requiring $`O(n_q^3)`$ elementary operations. From the state $`_{i=0}^{N1}b_i|\theta _i|\mathrm{cos}\theta _i`$, it is easy by using $`n_q`$ one-qubit operations to build the state $`_{i=0}^{N1}b_i\mathrm{exp}(ik\mathrm{cos}\theta _i)|\theta _i|\mathrm{cos}\theta _i`$. Then the cosines in the last register are reversibly erased by running backward the sequence of gates that constructed them, and one ends up with the state $`_{i=0}^{N1}b_i\mathrm{exp}(ik\mathrm{cos}\theta _i)|\theta _i|0`$, which is the result of the action of the unitary operator $`\mathrm{exp}(ik\mathrm{cos}\widehat{\theta })`$. Another QFT (requiring $`O(n_q^2)`$ operations) takes the wave function back to the $`n`$ representation.

Increasing $`n_q`$, which exponentially increases the dimension of the Hilbert space available, enables to probe various physical limits in the system. If $`K=kT`$ is kept constant, the classical mechanics remains the same. If $`T`$ is kept constant, the effective $`\mathrm{}`$ is fixed, and increasing $`n_q`$ will increase exponentially the size of the phase space of the system (number of cells). In contrast, if $`T=2\pi /N`$, with $`N=2^{n_q}`$, the size of the phase space remains the same, all $`N`$ momentum states corresponding to the same cell of size $`2\pi /T`$. In this case, increasing $`n_q`$ increases the number of quantum levels corresponding to the same classical structure, and is equivalent to decreasing $`\mathrm{}`$ toward the classical limit.

The whole quantum algorithm described above requires $`O(n_q^3)`$ gate operations to perform one iteration of the quantum map (2), exponentially less than the classical algorithm. Still, a physical quantum computer will not be an ideal perfect machine, and there will be imperfections, which may hamper the computation. In the following sections, we will investigate the effects of noise and imperfections on the physical quantities that are simulated, and estimate the accuracy of the quantum computation of the quantum map (2). The numerical simulation of many qubits is very resource-consuming on a classical computer. Due to that we took in all numerical computations the action of $`\mathrm{exp}(ik\mathrm{cos}\widehat{\theta })`$ as exact, and performed by direct multiplication in $`\theta `$ representation, all other operations being made with errors. We think that this approximation does not alter the qualitative features of the results, although the number of quantum gates is reduced from $`O(n_q^3)`$ to $`O(n_q^2)`$. Also in this approximation all supplementary registers required for the computation of $`|\mathrm{cos}\theta _i`$ are eliminated and the quantum evolution on $`N=2^{n_q}`$ levels is performed with only $`n_q`$ qubits.

## III Global quantities

We first study the effects of imperfections and errors for the global quantities of the system.

To model these imperfections, we introduce a random unitary error during the operation of elementary quantum gates. These errors are present for each gate performing the quantum Fourier transform and the action of the unitary operator $`e^{iT\widehat{n}^2/2}`$. Two elementary gates are used: single-qubit Hadamard gates $`H=\text{diag}(1,1,1,1)`$ and the two-qubit gate $`B=\text{diag}(1,1,1,\mathrm{exp}(i\alpha ))`$ where $`\alpha `$ is a phase. The transformation $`H`$ can be written $`H=\stackrel{}{u_0}\stackrel{}{\sigma }`$ where $`\stackrel{}{u_0}=(1/\sqrt{2},0,1/\sqrt{2})`$ and $`\stackrel{}{\sigma }=(\sigma _x,\sigma _y,\sigma _z)`$. It is replaced by an imperfect gate $`H^{}=\stackrel{}{u}\stackrel{}{\sigma }`$, where $`\stackrel{}{u}`$ is a unit vector with a random angle $`\beta `$ from $`\stackrel{}{u_0}`$. In a similar way, each $`B`$ is replaced by $`B^{}=\text{diag}(1,1,1,\mathrm{exp}(i\alpha +i\gamma ))`$ where $`\gamma `$ is again a random angle. At a given strength $`\epsilon >0`$ of noise, each gate is implemented with a $`\beta `$ or $`\gamma `$ randomly selected from a uniform distribution such that $`|\beta |<\pi \epsilon `$ or $`|\gamma |<\pi \epsilon `$ . As explained in Section II, we made the approximation of taking the action of $`\mathrm{exp}(ik\mathrm{cos}\widehat{\theta })`$ as exact, all other operations being made with errors. The use of $`n_q`$ qubits gives a Hilbert space for wave functions of the kicked rotator on $`N`$ levels, with $`N=2^{n_q}`$, i. e. values of momentum range from $`n=1`$ to $`n=N`$. In all numerical computations, the initial state $`|\mathrm{\Psi }_0`$ was chosen as localized on a precise value of the momentum $`n_0`$, i. e. $`|\mathrm{\Psi }_0=|n_0`$, with $`n_0=1`$ (lowest value of momentum) or $`n_0=N/2`$. The rotation is computed as $`\mathrm{exp}(iT(n\overline{n})^2/2)`$ with $`\overline{n}=N/2`$.

Depending on the choice of parameters in (2), increasing the number of qubits $`n_q`$ will increase the number of values of momentum in each phase space cell of size $`\mathrm{\Delta }n=2\pi /T`$, or increase the number of cells, or both. In , it was shown that if $`T`$ is constant while $`n_q`$ increase, errors in the QFT may lead to an exponential growth of errors with $`n_q`$ for the second moment $`n^2`$ of the probability distribution. In this case, the size of phase space grows exponentially with $`n_q`$ but $`K`$ and the effective $`\mathrm{}`$ are kept fixed. Due to quantum localization, exact wave functions cannot spread beyond a region of size given by the localization length, which remains fixed when $`n_q`$ increases. Therefore for all values of $`n_q`$, the second moment of a distribution initially located at $`n_0=N/2`$ will saturate with time at a value independent of $`n_q`$ (full line in Fig.2) if (2) is exactly simulated. On the contrary, errors in the QFT lead to small transfer of probabilities to regions of phase space which are exponentially far away from where the exact wave function is localized. This induces the exponential increase of the second moment with $`n_q`$. We confirm here this effect in Fig.2 for different parameters with the more complete set of errors used in this paper, and with simulations up to larger number of qubits.

To be more quantitative, Fig.3 shows the time scale $`t_q`$ on which the presence of errors leads to a doubling of the value of the second moment $`n^2`$ as a function of $`n_q`$ and error strength $`\epsilon `$. In  the formula:

$$t_qC_qk^4/(\epsilon ^2n_q2^{2n_q})$$
(3)
was proposed and checked numerically with up to 13 qubits. It stems from the fact that each imperfect gate operation transfers on average a probability of $`\epsilon ^2`$ equally divided among $`n_q`$ spurious peaks located at integer powers of $`2`$. Thus due to imperfections $`n^2n_q\epsilon ^22^{2n_q}t`$ (each time step involves $`n_q^2`$ gate operations), whereas for the exact wave function $`n^2D^24l^2k^4/4`$. Both expressions become comparable at the time $`t_q`$ given by (3). Fig.3 confirms this formula by extensive numerical computations, with up to 20 qubits, and for two different values of $`K`$. This enables to get the numerical constant $`C_q0.23`$.

Although the time scale $`t_q`$ drops exponentially with $`n_q`$, there are other observables which show only polynomial sensitivity to errors. A standard quantity used to characterize the global influence of errors is the fidelity defined by the projection of the wave function with errors $`\psi _\epsilon (t)`$ on the perfect one $`\psi _0(t)`$: $`f(t)=|\psi _\epsilon (t)|\psi _0(t)|^2`$. The dependence of this fidelity on time in presence of errors is shown in Fig.4, showing that it decreases slowly with $`t`$ and amplitude of noise $`\epsilon `$. One can define a time scale $`t_f`$ such that $`f(t_f)=0.5`$. Fig.5 presents the variation of $`t_f`$ with system parameters in two different régimes. It shows that the relation:

$$t_fC_f/(\epsilon ^2n_q^2)$$
(4)
holds with the numerical constant $`C_f0.35`$. Fig.4 and Fig.5 are consistent with a fidelity decay $`f(t)\mathrm{exp}(\mathrm{\Gamma }t)`$ where $`\mathrm{\Gamma }\epsilon ^2n_q^2`$.

The relation (4) can be understood from the following physical considerations. Each imperfect unitary gate is rotated by a random angle of order $`\epsilon `$ from the exact one. Therefore a probability of order $`\epsilon ^2`$ is transferred from the exact state at each gate operation. Each time step of the map (2) takes $`O(n_q^2)`$ operations, in the approximation which we have taken where the building of the cosines is supposed exact. This implies that $`t_f`$, which is in units of time steps of (2), should vary as $`1/(\epsilon ^2n_q^2)`$. We expect that if the full algorithm was implemented, with the cosines computed following the procedure explained in Section II, a time step of (2) should take $`O(n_q^3)`$ operations, and accordingly $`t_f`$ should vary as $`1/(\epsilon ^2n_q^3)`$.

The data shown in this section exemplify the sharp contrast in the behavior of the different observables in presence of errors. The fidelity shows only a polynomial decrease with respect to both $`\epsilon `$ and $`n_q`$, whereas the second moment of the wave function grows exponentially with $`n_q`$, but polynomially with $`\epsilon `$. The resolution of this apparent paradox is related to the fact that the second moment is sensitive to the size of the Hilbert space, which grows exponentially with $`n_q`$. Small spurious peaks due to imperfections do not spoil the fidelity, but modify strongly the variance $`n^2`$ if they appear very far away from the exact location of the wave function .

## IV Wigner and Husimi distributions

In the previous section we focused mainly on the case where $`T`$ (effective $`\mathrm{}`$) is fixed but the phase space size grows exponentially with $`n_q`$. In contrast, at $`T=2\pi /N`$ and $`N=2^{n_q}`$ the system size in classical momentum (number of $`2\pi /T`$ cells in $`n`$) remains fixed when $`n_q`$ increases. In this way the effective $`\mathrm{}`$ drops exponentially with $`n_q`$ and going to larger number of qubits means approaching the classical limit (exponentially fast). Smaller and smaller details of the classical structure will be visible in the quantum wave functions. In this régime, data presented in Fig.4 and Fig.5 have already shown that the fidelity follows the law (4) as in the case $`T=`$constant. However, the fidelity characterizes in one number the accuracy of the whole wave function, and does not tell how well the local properties are reproduced. To study the local properties of wave functions, one can express it in $`\theta `$ or $`n`$ representation. However, a very useful representation corresponds to phase space distributions, such as Wigner or Husimi distributions. They are especially used in the field of quantum chaos, since such representations permit a direct comparison with classical Hamiltonian mechanics which takes place in phase space. They also enable to probe the classical/quantum border when $`\mathrm{}`$ is decreased compare to other parameters of the system. Plotting such quantities in presence of errors allows to probe how local properties of the wave functions are sensitive to imperfections in the quantum algorithm.

An additional motivation to study such phase space representations stems from the fact that recently an algorithm was proposed  which enables to compute the Wigner function on a chosen point in phase space by the use of an ancilla qubit.

For a continuous system with two conjugate variables $`p`$ and $`q`$ the Wigner transform  of a wave function $`\psi `$ is defined by:

$$W(p,q)=\frac{e^{\frac{i}{\mathrm{}}p.q^{}}}{\sqrt{2\pi \mathrm{}}}\psi (q+\frac{q^{}}{2})^{}\psi (q\frac{q^{}}{2})𝑑q^{}$$
(5)
In a discrete system with $`N`$-dimensional Hilbert space, one is led to define the Wigner function on a lattice of $`2N\times 2N`$ points (see e. g. ). In the case of the kicked rotator, the formula becomes:

$$W(\theta ,n)=\underset{m=0}{\overset{N1}{}}\frac{e^{\frac{2i\pi }{N}n(m\mathrm{\Theta }/2)}}{2N}\psi (\mathrm{\Theta }m)^{}\psi (m),$$
(6)
with $`\mathrm{\Theta }=\frac{N\theta }{2\pi }`$. The Wigner function is always real, but contrary to classical Liouville phase space distributions it can take negative values. It verifies $`_iW(\theta _i,n)=|\psi (n)|^2`$ and $`_iW(\theta ,n_i)=|\psi (\theta )|^2`$.

The Wigner transform has the drawback of being negative or positive. Nevertheless, coarse-graining this function over cells of size $`\mathrm{}`$ gives non negative values. Such a procedure gives the Husimi distribution (see e.g. ) which corresponds to a Gaussian smoothing of the Wigner function. In the case of the kicked rotator, the Husimi distribution can be computed through:

$$h(\theta ,n)=\underset{m=nN/2}{\overset{n+N/2}{}}(\frac{T}{\pi })^{\frac{1}{4}}\frac{\psi (m)}{\sqrt{N}}e^{\frac{T}{2}(mn)^2}e^{im\theta }$$
(7)
where the Gaussian for simplicity is truncated for values larger than $`N/2`$, and $`\psi (m)`$ is the wave function in momentum representation. The Husimi distribution is always non negative, and allows a direct comparison between classical Liouville density distributions and quantum wave functions.

Wigner and Husimi distributions of wave functions of the quantum kicked rotator simulated on a quantum computer are shown on Fig.6 for different level of errors. Both functions have similar patterns, although as expected the Wigner function displays interference structures absent in the Husimi distribution. In the regime of parameters studied, classical invariant curves are still present in phase space and prevent the exact wave function to enter the large elliptical island in the middle. In the presence of moderate level of noise, main structures are still present and distinguishable.

Fig. 7 confirms this result, showing the Husimi distribution for larger number of qubits, together with the classical phase space distribution. The Husimi distributions in phase space show features mimicking the classical phase space distributions, in accordance with the correspondence principle. Fig.7 (left) shows that when $`n_q`$ is changed, finer and finer details of the classical structures are visible in the exact quantum wave function, in accordance with the fact that increasing $`n_q`$ amounts to reduce $`\mathrm{}`$ and approach the classical limit. The same figure shows that the wave function is spread over a larger domain of phase space as $`n_q`$ increases. This can be explained by the following effect. In this mixed régime between integrability and strong chaos at $`K=1.3`$ the invariant classical curves which prevent any transport are no longer present since the last one is destroyed at $`K=K_g=0.97..`$. But cantori are present, which are remnants of the disappeared invariant curves. They have a fractal structure, and a wave packet can cross them only if the holes are large enough. These holes scale as $`(KK_g)^3`$ and become comparable with the minimal area scale of the Husimi distribution determined by the effective $`\mathrm{}`$ given by $`T`$. Hence, for $`KK_g1`$, the wave function is prevented to cross the cantorus for $`(KK_g)^3<T`$. Due to that quantum interference prevents the transport via cantori .

The quantum Husimi distributions shown in Fig. 6 and Fig.7 display structures of increasing complexity with larger $`n_q`$. Still, with moderate level of noise, the quantum computer is able to reproduce the exact distributions with reasonable accuracy. For larger errors in gate operations, significant probability is present at wrong phase space locations, and phase space structures become blurred. The comparison with the effect of classical noise visible in Fig.7 shows that in this case the quantum errors enable the wave function to enter classically forbidden zones much faster, a fact which will be analyzed in more details in the next section.

It is interesting to evaluate the effects of noise and imperfections not only on the broad features of the full Wigner function, but also on individual values. In Fig.8 and Fig.9, the behavior of individual values of the Wigner function in presence of noise in the gates is investigated. Fig.8 shows that the relative error (i.e. the error $`|WW_\epsilon |`$ divided by the average individual value of the exact Wigner function $`|W|`$) increases slowly with the growth of $`t`$ and $`\epsilon `$ even in the chaotic zone. Similar results can be observed in the integrable zone and in the localized régime (data not shown). In a more quantitative way, Fig.9 shows the behavior of the time scale $`t_W`$ when the error on the Wigner function become comparable to its mean value in the régime chosen ($`|W(t_W)W_\epsilon (t_W)|=|W|/2`$). In all three cases considered, one obtains

$$t_WC_W/(n_q^\alpha \epsilon ^2),$$
(8)
with $`\alpha =1`$ or $`\alpha =1.5`$. Thus individual values of the Wigner function in the kicked rotator model are robust quantities with respect to noise, even in the chaotic régime. These results are interesting also in view of the recent discussion on the effects of decoherence on Wigner functions . Our results clearly show that in the framework of quantum computation, the errors on the Wigner function are polynomial and not exponential.

As noted previously, a recent algorithm  enables to measure the value of the Wigner function of a system of density matrix $`\rho `$ on a selected point in phase space, with the help of an ancilla qubit $`a`$. First $`H`$ (Hadamard gate) is applied on $`a`$, followed by a controlled-$`U`$ operation ($`U`$ is applied to the system to be measured depending on the state of $`a`$) and again $`H`$ is applied on $`a`$. Then the expectation value of $`a`$ is $`\sigma ^z=Re[Tr(U\rho )]`$. The use of a particular operator $`U`$, which can be implemented efficiently , enables to get $`W(p,q)=\sigma ^z/2N`$ (where $`N=2^{n_q}`$).

However, we should note that even if $`U`$ can be implemented efficiently, $`\sigma ^z`$ can be evaluated only by iterating the procedure enough times to get a good estimate. Therefore the amplitude of the signal is crucial to make the whole process efficient. Thus it is interesting to study the amplitude of peaks in the Wigner function, in order to know if strong peaks are present which can be detected reliably through this method. This can be investigated through a quantity which we call inverse participation ratio of the Wigner function, in analogy with the inverse participation ratio for wave functions used in quantum chaos and systems with Anderson localization . For a wave function with $`N`$ projections $`\psi _i`$ on some basis, the inverse participation ratio $`|\psi _i|^2/(|\psi _i|^4)`$ measures the number of significant components in this basis. For the Wigner function, one has the additional sum rules $`W_i=1`$ and $`W_i^2=\frac{1}{N}`$. To define an inverse participation ratio for the Wigner function, we therefore use the formula $`\xi =1/(N^2W_i^4)`$. If $`N`$ peaks of approximately equal weights $`1/N`$ are present, then $`\xi =N`$, whereas $`N^2`$ components of equal weights (in absolute value) $`1/N^{3/2}`$ give $`\xi =N^2`$. This quantity $`\xi `$ therefore permits to estimate the number of main components of the Wigner function. Fig.10 and Fig.11 show the scaling of this quantity with $`n_q`$ for different values of parameters. In all the cases where $`T=2\pi /N`$ ($`N=2^{n_q}`$) the ratio $`\xi /N^2`$ reaches a saturation value. This implies that asymptotically $`\sigma _z=NW(p,q)1/\sqrt{N}`$, a value which requires $`N`$ iterations followed by measurements to be reliably estimated. In this case, the asymptotic gain in number of operations compared to the classical algorithm is only $`O(\mathrm{log}(N))`$, although the resources needed are exponentially smaller ($`n_q`$ qubits instead of $`2^{n_q}`$ classical registers). This should be contrasted with the case where the number of cells increases ($`T`$ constant) where Fig.11 shows that $`\xi N`$. This gives $`\sigma _z=NW(p,q)1`$, which means that in this regime with localization, any of the $`N`$ components of the Wigner function which are important can be estimated reliably and efficiently through this method (provided one knows beforehand the approximate position of the localized state). The results presented in Fig.8 and Fig.9 show that despite the different scaling laws of $`|W|`$, the relative errors grow only polynomially in all cases considered, thus enabling such measurements of individual values of $`W`$ to be reliable for moderate amounts of noise.

## V Quantum tunneling through invariant curves

In the previous section, it was shown that the classical and quantum errors affect the dynamics in a rather different way. This difference is particularly striking in the régime where classical invariant curves are present (integrable or mixed systems, which correspond to moderate values of $`K`$ here, as in Fig.7). Such invariant curves cannot be crossed classically, and only quantum tunneling can transfer probability inside integrable islands from chaotic regions. However, whereas small classical errors enable to cross only neighboring invariant curves, small quantum errors may lead to long-distance “jumps” of probability deep into integrable island (see Fig.7, last column).

To study the effects of errors on quantum tunneling we show in Fig.12 the dependence of probability of Husimi distribution $`h(\theta ,n)`$ inside the classically forbidden region on time $`t`$. The quantity $`I(t)=_Dh(\theta ,n)𝑑\theta 𝑑n`$, where $`D`$ is the domain enclosed by the circle in Fig.10 (inset), shows a linear growth with $`t`$. This can be understood by a physical argument similar to the one justifying (4). Indeed, imperfect gates transfer on average a probability of order $`\epsilon ^2`$ from the exact wave function to wrong phase space positions. However, not all gates will transfer probability inside $`D`$ but only a subset of them. This predicts that $`I(t)n_q^\alpha \epsilon ^2t`$. Data from Fig.12 and Fig.13 and additional data (not shown) confirm this prediction, with $`\alpha 1.3`$.

To exemplify the effect of quantum errors, Fig.13 shows $`I(t)`$ at fixed time $`t`$ as a function of number of qubits $`n_q`$ for zero and nonzero noise in the gates. In the case of zero noise, there is an exponential decrease with $`n_q`$. Indeed, the only process which allows to enter the island for the wave packet is quantum tunneling. In general, the probability of such a transition scales like $`\mathrm{exp}(S/\mathrm{})`$ where $`S`$ is a classical action. Increase of $`n_q`$ amounts to decrease the effective $`\mathrm{}`$ and leads to the exponential drop of $`I`$ obtained numerically at $`\epsilon =0`$. In sharp contrast, the presence of imperfections in the gates ($`\epsilon >0`$) leads to direct jumps inside the island that gives an increase of $`I`$ with $`n_q`$ according to the estimate of the previous paragraph. Thus for this specific process, the effect of noise in the gates results in a qualitative change of the dependence of tunneling probability $`I`$ on $`n_q`$.

## VI Conclusion

The results presented in this paper show that it is possible to simulate efficiently the quantum kicked rotator on a quantum computer. For the quantum algorithm simulating the dynamics of kicked rotator we investigated the effects of gate errors and showed that certain quantities like fidelity, Wigner and Husimi distributions are sufficiently robust against noise in the gates. Thus for small amplitude of noise these quantities can be computed reliably without application of quantum error corrections. At the same time we found that there exist other characteristics, e.g. variance of probability distribution and tunneling probability inside stability islands, which are very sensitive to errors in quantum gates. In addition, the study of the Wigner function shows that individual values of this function are robust with respect to quantum errors and can be reliably estimated. However, the computation of the Wigner function at specific points meets certain readout problems in deep quasiclassical regime where generally a large number of measurements is required.

On the basis of obtained results we believe that the quantum algorithms simulating quantum chaotic maps will provide important grounds for testing the accuracy of the next generation of experimental implementations of quantum computers.

We thank the IDRIS in Orsay and CalMiP in Toulouse for access to their supercomputers. This work was supported in part by the NSA and ARDA under ARO contract No. DAAD19-01-1-0553, by the EC RTN contract HPRN-CT-2000-0156 and by the project EDIQIP of the IST-FET programme of the EC.


text:
# Observations of possible symbiotic star V335 Vul

## 1. Introduction

V335 Vul = (AS 356 = GSC 02128.00676 = LD 120 = IRAS 19211+2421) was first noticed by Merril and Burwell (1950) as a star with an apparent H-$`\alpha `$ emission. Nassau and Blanco (1957) resolved V335 Vul as a cool carbon star. This was confirmed by Stephenson (1971).

The star was found to be variable by Collins (1991) on photographic plates with maximum magnitude at 10.1 mag and minimum under the limiting magnitude of the plates at 12.7 mag. Collins also determined the type of variability to be SR with ephemeris like Max = 2447767 + 336 x $`E`$. This information caused the appearance of the star in the 71st name-list of variable stars (Kazarovets et al. 1993). The independent discovery of variability came from Dahlmark (1993). From his photographic measurements between 1985 and 1992, he could derive following ephemeris: 2446740 + 342.0 x $`E`$. He also noted that the period of light changes was stable during 1967 and 1992 and found $`BV`$ = 3.3 mag. The type of variability was determined to be on the edge between Mira and SRa.

Munari, Tomov and Rejkuba (1999) argued (using echellete spectra) that V335 Vul is probably a symbiotic binary. Munari et al. (1999) discovered an apparent outburst of V335 Vul with significant change of colour index $`BV`$ from 5.01 to 3.06 mag. Belczynski et al. (2000) attributed V335 Vul to stars suspected from symbiotic nature. Munari and Jurdana-Sepic (2002) examined Asiago photographic plates for data on V335 Vul, but insufficient coverage of two years by only 26 points in $`B`$ or $`V`$ filters did not permit detailed analysis.

## 2. Observations

Since the spring of 2000, MEDUZA group has been conducting extensive observing campaign on V335 Vul. This include long-term VRI CCD measurements from N. Copernicus Observatory and planetarium in Brno by P. Sobotka and O. Pejcha, long-term $`V`$ band CCD measurements conducted by L. Smelcer at Valasske Mezirici Observatory, high-speed $`V`$ band CCD photometry on three nights made at Hradec Kralove Observatory by P. Marek, and visual observations carried out by amateur astronomers in the Czech and Slovak Republics, mostly by P. A. Dubovsky, L. Brat and O. Pejcha. As a result of the campaign, 1086 CCD long-term measurements in all filters and 192 visual estimates have been made.

CCD observations from Brno and Hradec Kralove observatories were reduced by the means of Munidos software (Hroch, Novak and Kral, 2002); observations from Valasske Mezirici were reduced using CCDOPS software bundled with SBIG CCD cameras. CCD data from Brno Observatory have been transformed to the standard Johnson-Cousins system, CCD data from Valasske Mezirici Observatory have not. Systematical error of the transformations should not exceed 0.1 mag level.

All observations were made with respect to the precisely measured comparison stars of Henden and Munari (2000). The only exceptions are the two brightest stars used only for visual observations - their magnitudes were taken from the Tycho and Tycho-2 catalogues.

Before doing any analysis, visual observations have been averaged by moving averages (with 50 days width and 1 day step) using MedPrum procedure written by M. Haltuf. For plot clarity, weighted averages of CCD observations made on the same night with the same filter and by one observer were done. Times of minima and maxima were derived by Kwee and Van Woerden method using program AVE by R. Barbera.

## 3. Results

High-speed photometry did not reveal anything unusual above 0.02 mag level. This is not unexpected - amplitudes of flickering in the $`V`$ band tend to be very small among symbiotic binaries. This is not the best way how to prove or disprove the symbiotic nature of the star.

In Figure 1 we present averaged VRI light curves of V335 Vul as obtained from the multi-colour observations by P. Sobotka and O. Pejcha (and also L. Smelcer in the case of $`V`$ filter). All our observations suggest ordinary Mira or SRa type variations with slightly different depths of minima. The observed maximum is symmetric in all passbands. Maximum and minimum brightness vary in time, see Table 1.

We conclude that the outburst announced by Munari et al. (1999) is identical with ordinary Mira-like pulsations. If V335 Vul is a symbiotic binary, it is similar e.g. to R Aqr.

In Figure 2 we present changes of $`RI`$ and $`VR`$ colour indices compared with $`V`$ band light changes. As expected, the star is redder near minimum light. Changes of colour indices are presented in Table 2.

Times of maximum depends on wavelength: $`V=2452169\pm 1`$, $`R=2452173\pm 1`$ and $`I=2452183\pm 1`$. Maxima of $`R`$ and $`I`$ band curves occurred 4 and 14 days after the $`V`$ band maximum, respectively.

Observed $`V`$ band maxima occurred at JD $`2451822.2\pm 0.5`$ and $`2452169\pm 1`$ leading to the period of approx. 347 days, which is similar to previous works.

All visual observations as well as the chart with sequence can be downloaded via MEDUZA webpage at http://www.meduza.info.

### Acknowledgments.

We are thankful to all visual and CCD observers who dedicated their time to V335 Vul and L. Kral and M. Haltuf for software support. This research has made use of NASA’s Astrophysics Data System.

## References

Belczynski, K., et al. 2000, A&A, 146, 407

Collins, M., 1991, The Astronomer, 28, 182

Dahlmark, L., 1993, IBVS 3855

Henden, A., Munari, U., 2000, A&A, 143, 343

Hroch, F., Novak, R., & Kral, L., 2002, http://munipack.astronomy.cz

Kazarovets, E. V., Samus, N. N., & Goranskij, V. P., 1993, IBVS 3840

Merrill, P. W., Burwell, C. G., 1950, ApJ, 112, 72

Munari, U., Jurdana-Sepic, R., 2002, www.arxiv.org/astro-ph/0202118

Munari, U., Tomov, T., Moro, D., Henden, A., 1999, IAUC 7329

Munari, U., Tomov, T., Rejkuba, M., 1999, IBVS 4668

Nassau, J. J., Blanco, V. M., 1957, ApJ, 125, 195

Stephenson, C. B., 1971, PWSO, 1, 4


text:
# Electromagnetic focusing through a tilted dielectric surface

## I Introduction

Focusing of light has drawn considerable attention for many centuries due to its importance in many fields of physics. Microscopy, spectroscopy and optical data storage are only a few of the fields utilizing focused light. The earliest treatment of electromagnetic focusing problems are due to IgnatowskyIgnatowsky . However, the structure of the focused field was not clarified until the 1950’s and 60’s, when Wolf and coworkers did the first detailed studies of aplanatic electromagnetic focusing systemsWolf ; Richards ; Boivin . Here the socalled Debye approximation was adopted, where only the plane waves with propagation vectors that fall inside the geometrical cone whose apex is at the focal point contribute to the field in the focal regionWolf ; Debye . The Debye approximation predicts that the electromagnetic field is symmetrical about the focal plane. However, later it was found that the Debye theory is valid only when the focal point is located many wavelengths away from the aperture, the angular aperture is sufficiently large, and the wavefront deviation kept reasonably lowWolf1 ; Stamnes ; Stamnes1 ; Sheppard ; Li . In the general case the electromagnetic field is not symmetrical about the geometrical focal plane, and focal shifts may occur.

Focusing through dielectric interfaces is of broad interest since this geometry is often used in optical data storage and biology (see e.g. Refs. Stamnes2 ; Dhayalan ; Sheppard1 ; Wilson ; Torok ; Ando ; Hendriks ; Helseth  and references therein). To date, most studies have concentrated about focusing through a dielelectric interface with its normal along the optical axis. However, in some applications it is of importance to understand what happens when the interface is tilted. For example, consider an optical disk which rotates around its axis in a DVD player. If the disk is slightly tilted (e.g. by shaking the DVD player), then the wavefront is aberrated, and the focus is distorted. To find out how much the focus is distorted, it is most common to perform optical ray tracing or use scalar wave theory (see e.g. Refs. Ando ; Hendriks  and references therein). However, also in the case of electromagnetic focusing it is of interest to obtain an explicit expression for the intensity distribution in the focal region, since this may give a clearer understanding of the physical problem. The aim of this Letter is to treat the problem of electromagnetic focusing through a tilted dielectric interface in more detail. First I will formulate the problem, and then look at the small angle limit in the Debye approximation. It should be emphasized that no numerical analysis is given here.

## II General formulation

In general, the diffracted field near the focal plane can be calculated in the Kirchoff approximation as Stamnes2

$$𝑬=\frac{ik_i}{2\pi }_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑻(𝒔_i)\mathrm{exp}[ik_i(s_{ix}x+s_{iy}y+s_{iz}z)]𝑑s_{ix}𝑑s_{iy},$$
(1)
where $`k_i=2\pi n_i/\lambda `$ is the wavenumber in medium i, $`𝒔_i=(s_{ix},s_{iy},s_{iz})`$ is the unit vector along a typical ray, and $`𝑻`$ the vector pupil distribution which accounts for the polarization, phase and amplitude distributions at the exit pupil. Let us consider focusing through a tilted isotropic dielectric interface when the optical axis crosses the interface at a distance $`z_i=d`$ (see Fig. 1). Without loss of generality, we let the angle $`\beta `$ be the tilting angle about the y-plane. Thus, the z-coordinate at the interface is given by $`z_i=d+xtan\beta `$. When $`\beta =0`$ we have $`𝒔_1=(s_{1x},s_{1y},s_{1z})`$. The unit vector corresponding to a finite tilt angle is then given by a rotation of the coordinate system $`𝒔_1^{}=𝑲𝒔_1`$, where

$$𝐊=\left[\begin{array}{ccc}cos\beta & 0& sin\beta \\ 0& 1& 0\\ sin\beta & 0& cos\beta \end{array}\right],$$

and $`𝒔_1^{}=(s_{1x}cos\beta s_{1z}sin\beta ,s_{1y},s_{1x}sin\beta +s_{1z}cos\beta )`$. Here we have used the equivalence between tilting the optical axis and the dielectric interface for a single ray. We will see that this change of coordinate system is useful for expressing the electric field in the second medium. Moreover, Snell’s law is now expressed as; $`n_1s_{1y}=n_2s_{2y}`$ and $`n_1(s_{1x}cos\beta s_{1z}sin\beta )=n_2s_{2x}`$. We note that these equations are not valid when $`n_1=n_2`$ and $`\beta 0`$, since in this case we can not immediately assume equivalence between tilting the (absent) substrate and the optical axis.

It is convenient to write the electric fields in the first and second media as (see also Ref. Torok )

$$𝑬_1=\frac{ik_1}{2\pi }_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\frac{𝑻(𝒔_1)}{s_{1z}}\mathrm{exp}[ik_1(s_{1x}x+s_{1y}y+s_{1z}z)]𝑑s_{1x}𝑑s_{1y},$$
(2)
and

$$𝑬_2=\frac{ik_2}{2\pi }_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑻(𝒔_2)\mathrm{exp}[ik_2(s_{2x}x+s_{2y}y+s_{2z}z)]𝑑s_{2x}𝑑s_{2y}.$$
(3)
To express the electric field in medium 2 in terms of ’untilted’ coordinates, $`(s_{1x},s_{1y})`$, we note that

$$ds_{2x}ds_{2y}=Jds_{1x}ds_{1y},$$
(4)
where the Jacobian determinant is given by

$$J=\left(\frac{k_1}{k_2}\right)^2\left(cos\beta +\frac{s_{1x}}{s_{1z}}sin\beta \right).$$
(5)
Moreover, we must match the electric field at the interface, which gives

$$𝑻_2=𝑭(𝒔_1,\beta )\frac{k_1𝑻(𝒔_1)}{Jk_2s_{1z}}\mathrm{exp}\left[iz_i(k_1s_{1z}k_2s_{2z})\right],$$
(6)
where $`𝑭(𝒔_1,\beta )`$ is the transmission factor (Fresnel coefficients) through the tilted surface. Now we may write

$$𝑻_2=\stackrel{\mathbf{~}}{𝑻}(𝒔_1,\beta )\mathrm{exp}(i\mathrm{\Psi }_M+i\mathrm{\Psi }_T),$$
(7)
where

$$\stackrel{\mathbf{~}}{𝑻}(𝒔_1,\beta )=\frac{k_2}{k_1}\frac{𝑭(𝒔_1,\beta )𝑻(𝒔_1)}{\left(cos\beta +\frac{s_{1x}}{s_{1z}}sin\beta \right)s_{1z}}.$$
(8)
$`\mathrm{\Psi }_M`$ and $`\mathrm{\Psi }_T`$ are given by

$$\mathrm{\Psi }_M=d(k_1s_{1z}k_2s_{2z}),$$
(9)
and

$$\mathrm{\Psi }_T=xtan\beta (k_1s_{1z}k_2s_{2z}).$$
(10)
Naturally, $`\mathrm{\Psi }_T`$ vanishes when there is no index mismatch. Note also that $`s_{2z}=\sqrt{1s_{2x}^2s_{2y}^2}`$, where $`s_{2x}`$ and $`s_{2y}`$ are given by Snells law. Finally, the expression for the electric field in medium 2 is

$$𝑬_2=\frac{ik_2^2}{2\pi k_1}_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑻_2\mathrm{exp}\left[ik_1(s_{1x}cos\beta s_{1z}sin\beta )x+ik_1s_{1y}y+ik_2s_{2z}z\right]𝑑s_{1x}𝑑s_{1y}.$$
(11)
This is expression enable us to calculate the electric field for reasonable tilt angles in the Kirchoff approximation. In the next section I will study the special case of high angular apertures and small tilt angles, and it will be shown that only a simple aberration term remains.

## III Small tilt angles in the Debye approximation

We now consider focusing with high angular aperture assuming that the Debye approximation can be used. Then only the plane waves with propagation vectors that fall inside the geometrical cone whose apex is at the focal point contribute to the field in the focal region. Moreover, we assume that $`\beta 1`$ and $`z=0`$. Then the influence of $`\beta `$ on $`\stackrel{~}{𝑻}(𝒔_1,\beta )`$ can be neglected, and only the lowest order contribution in the phase remains. Equation (11) can therefore be written as

$$𝑬_2=\frac{ik_2^2}{2\pi k_1}_\mathrm{\Omega }\stackrel{~}{𝑻}(𝒔_1)\mathrm{exp}(i\mathrm{\Psi }_T^{}+i\mathrm{\Psi }_M^{})\mathrm{exp}\left[ik_1s_{1x}(xu)+ik_1s_{1y}y\right]𝑑s_{1x}𝑑s_{1y},$$
(12)
where $`\mathrm{\Omega }`$ is the solid angle formed by all the geometrical rays,

$$\mathrm{\Psi }_T^{}=x\beta k_2s_{2z0},$$
(13)
$$\mathrm{\Psi }_M^{}=dk_1s_{1z}+dk_2s_{2z0},$$
(14)
$`u=d\beta (k_1/k_2)(s_{1z}/s_{2z0})`$ and $`s_{2z0}=\sqrt{1(n_1/n_2)^2(s_{1x}^2+s_{1y}^2)}`$. These expression were found by expanding $`s_{2z}`$ and only keeping the lowest order in $`\beta `$. Note that $`\mathrm{\Psi }_M^{}`$ is the usual aberration introduced by focusing through a planar, nontilted dielectric surface, see Ref. Torok . This aberration term will be neglected here, which is a reasonable approximation for systems corrected for the index mismatch (in absence of tilting). It is seen that the tilt introduces a coordinate shift, u, which depends on the tilt and index mismatch. Note in particular that $`ud\beta `$ when $`n_1n_2`$, which is just the shift expected from a tilt in the coordinate system. Moreover, u is clearly altered when the index mismatch and numerical aperture change. In the further studies I will neglect this shift, concentrating on the remaining aberration term $`\mathrm{\Psi }_T^{}`$. This term is also nonzero when $`n_1n_2`$. However, I argue that it has the correct symmetry, and is a real aberration term due to the finite tilt $`\beta `$.

I now assume that the aperture is circular symmetric, which means that it is most convenient to adopt spherical coordinates:

$$𝒔_i=(sin\theta _icos\varphi ,sin\theta _isin\varphi ,cos\theta _i),$$
(15)
and

$$𝒓_c=(r_csin\theta _ccos\varphi _c,r_csin\theta _csin\varphi _c,z).$$
(16)
Equation (12) can now be written as

$`𝑬_2`$ $`=`$ $`{\displaystyle \frac{ik_2^2}{2\pi k_1}}{\displaystyle _0^\alpha }{\displaystyle _0^{2\pi }}\stackrel{~}{𝑻}(\theta _1,\varphi )\mathrm{exp}(ik_2s_{2z0}\beta r_csin\theta _ccos\varphi _c)`$ (18)
$`\mathrm{exp}[ik_1r_csin\theta _1sin\theta _csin\theta _1cos(\varphi \varphi _c)]d\varphi d\theta _1,`$
In order to be able to evaluate this integral, it remains to find the vector pupil function $`\stackrel{~}{𝑻}(\theta _1,\varphi )`$. This can be done following the guidelines given in e.g. Ref. Helseth , which results in

$$\stackrel{~}{𝐓}(\theta _1,\varphi )=A(\theta _1)\left[\begin{array}{ccc}a[t_pcos\theta _2cos^2\varphi +t_ssin^2\varphi ]+b[t_pcos\theta _2sin\varphi cos\varphi t_ssin\varphi cos\varphi ]& & \\ a[t_pcos\theta _2cos\varphi sin\varphi t_ssin\varphi cos\varphi ]+b[t_pcos\theta _2sin^2\varphi +t_scos^2\varphi ]& & \\ t_psin\theta _2[acos\varphi +bsin\varphi ]& & \end{array}\right],$$

where $`A(\theta _1)`$ is an apodization factor, $`a(\theta _1,\varphi )`$ is the strength of the incident x polarized light, $`b(\theta _1,\varphi )`$ the strength of the incident y polarized light and $`t_{p,s}`$ the Fresnel transmission coefficients. In the case of linearly polarized light (a=1, b=0) we obtain the following electric field components:

$$E_xi\left(I_0+I_2cos2\varphi _c\right),$$
(19)
$$E_yiI_2sin2\varphi _c,$$
(20)
$$E_z2I_1cos\varphi _c,$$
(21)
where

$$I_0=_0^\alpha A(\theta _1)(t_s+t_pcos\theta _2)sin\theta _1J_0(kr_csin\theta _1sin\theta _c)\mathrm{exp}(ik_2s_{2z0}\beta r_csin\theta _ccos\varphi _c)𝑑\theta _1,$$
(22)
$$I_1=_0^\alpha A(\theta _1)t_psin\theta _1sin\theta _2J_1(kr_csin\theta _1sin\theta _c)\mathrm{exp}(ik_2s_{2z0}\beta r_csin\theta _ccos\varphi _c)𝑑\theta _1,$$
(23)
$$I_2=_0^\alpha A(\theta _1)(t_st_pcos\theta _2)sin\theta _1J_2(kr_csin\theta _1sin\theta _c)\mathrm{exp}(ik_2s_{2z0}\beta r_csin\theta _ccos\varphi _c)𝑑\theta _1.$$
(24)
We see that now the only effect of the tilt is to introduce an aberration term which depends on the azimuthal angle $`\varphi _c`$ as well as $`\theta _1`$. When $`\beta =0`$, these equations reduces to the ones found in e.g. Refs. Torok ; Helseth .

In conclusion, focusing through a tilted dielectric interface has been investigated. It is found that the tilt introduces additional aberration terms. In the small angle limit and Debye approximation, only a simple aberration term remains.