text:
# Fast Convolution on the Sphere

## I Introduction

A major near-term objective in the field of Cosmology today is to gain a detailed measurement and statistical understanding of the anisotropies of the cosmic microwave background (CMB). While the theory of primary CMB anisotropy is well-developed (see  for a review) and we are facing a veritable flood of data from a new generation of instruments and missions, perhaps the single most limiting factor for interpreting these data is the exorbitant computational cost involved in realistic mission simulation and careful analysis of the data products.

Important and computationally expensive tasks for both simulation and analysis of microwave data are to simulate and to correct for the systematic errors due to imperfections of realistic microwave telescopes, such as beam asymmetries and far side lobes. The effect of an asymmetric “point spread function” is to distort the shapes of the detected anisotropies. What makes far side lobes an important issue is the fact that the CMB anisotropy signal has an amplitude of one in $`10^5`$ relative to the $`2.7`$K background. In regions of low galactic latitude, foregrounds from galactic synchrotron radiation and dust emission are expected to exceed this signal by many orders of magnitude over a wide range of frequencies . Even though CMB experiments will obviously not target these regions to obtain measurements of the background anisotropy, the large amplitudes of these galactic sources may induce systematic errors even when “looking” in directions far away from the galactic plane if the instrument allows diffraction of stray light into the detectors. Solar system bodies, including the earth, are other possible sources of stray light.

To assess these problems and formulate solutions we must be able to compute the detector response at every pointing of the telescope. The inputs are a physical model of the “beam” over $`4\pi `$ steradian and a model of the “sky” containing both simulated signal as well as foreground sources possibly including ground emission. Note that in the general case not just the direction of the pointing is important but also the orientation of the beam about the pointing axis. The detector response is then the solution to a quadrature problem at each orientation.

Analysis methods of CMB data have neglected this difficulty by assuming azimuthal symmetry of the beam which greatly simplifies the calculation . Simulation work which did include an asymmetric beam and far side lobes using pixel based methods  ran up against computational challenges for angular scales smaller than one degree, running for hundreds of hours even with optimized adaptive mesh algorithms. Such algorithms are clearly inadequate for modern high resolution experiments which achieve resolutions of a few minutes of arc.

In this paper we describe a numerical method which greatly accelerates the computations which are necessary to correctly account for realistic beam profiles in simulation and analysis of directional data on the sphere. This is achieved by rewriting the problem in such a way that we can take advantage of the Cooley-Tukey Fast Fourier Transform (FFT) algorithm.

The following section of this paper defines the general problem in terms of rotations of the beam with respect to the sky. We then introduce a geometrically motivated split of the rotation operator in section three. This enables us, in section four, to derive the general solution for the detector response for all possible relative orientations of the beam and the sky within a given section on the sphere. Section five then discusses the solution and derives special cases from it, amongst others the well-known algorithm for convolution with azimuthally symmetric kernels. We conclude in section six with a timing example. An appendix applies the same ideas to accelerating the computation of the transverse convolution, an operation which becomes important in the inverse problem of map estimation.

While we are motivated by the goal of achieving and interpreting precision measurements of the anisotropies of the cosmic microwave background, the methods we present are general and apply to the convolution or filtering of any scalar field on the sphere with an arbitrary, asymmetric but constant kernel. We generalise our methods to tensor fields on the sphere in reference .

## II Statement of the Problem

Consider two bandlimited functions on the sphere $`b(\stackrel{}{\gamma })`$ and $`s(\stackrel{}{\gamma })`$. For definiteness and to aid the imagination we will refer to them in the following as the beam and the sky, respectively, but they could be completely general bandlimited functions — in particular neither of them is constrained to be positive definite or even real.

The task is to compute the scalar product of the beam and the sky at a set of beam orientations. To describe these orientations, we use the Euler angles $`\mathrm{\Phi }_1,\mathrm{\Theta }`$ and $`\mathrm{\Phi }_2`$***Our Euler angle convention refers to active right handed rotations of a physical body in a fixed coordinate system. The coordinate axes stay in place under all rotations and the object rotates around the $`z`$, $`y`$ and $`z`$ axes by $`\mathrm{\Phi }_1`$, $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }_2`$, respectively, according to the right handed screw rule.. The convolved signal for each beam orientation $`(\mathrm{\Phi }_1,\mathrm{\Theta },\mathrm{\Phi }_2)`$ can then be written as

$$T(\mathrm{\Phi }_2,\mathrm{\Theta },\mathrm{\Phi }_1)=𝑑\mathrm{\Omega }_\stackrel{}{\gamma }[\widehat{D}(\mathrm{\Phi }_2,\mathrm{\Theta },\mathrm{\Phi }_1)b](\stackrel{}{\gamma })^{}s(\stackrel{}{\gamma }).$$
(1)
Here the integration is over all solid angles, $`\widehat{D}`$ is the operator of finite rotations such that $`\widehat{D}b`$ is the rotated beam, and the asterisk denotes complex conjugation.

If $`(\mathrm{\Phi }_1,\mathrm{\Theta },\mathrm{\Phi }_2)`$ can be written as a continuous function of a parameter $`t[0,T]`$, say, then we call the ordered set of tuples $`(\mathrm{\Phi }_1(t),\mathrm{\Theta }(t),\mathrm{\Phi }_2(t))`$ a scan path. Note that Eq. (1) assumes that time varying signals in the sky vary either on time scales much longer than the duration of the scan or much smaller than the integration time per sample. In the context of CMB missions this is a good approximation with the exceptions of planets (for long duration missions), time varying point sources, and atmospheric foregrounds. Of these only atmospheric foregrounds present a problem for the convolution, because they are extended - convolution with a point source is a simple operation in position space (the pixel basis) and can be computed separately. Linearity allows us to then add the results of the $`r\pi `$ convolution of extended sources to the point source convolution.

In the most general case, the bandlimits (see Eq. (5) below for a definition) of the beam and the sky are $`𝖫_b`$ and $`𝖫_s`$, respectively. Define $`𝖫\mathrm{min}(𝖫_b,𝖫_s)`$. Note that we actually only need one of $`s`$, $`b`$ to be bandlimited as long as the multipoles of the other are bounded as $`l\mathrm{}`$. Then the numerical evaluation of the integral in Eq. (1) takes $`𝒪\left(𝖫^2\right)`$ operations for each tuple $`(\mathrm{\Phi }_1,\mathrm{\Theta },\mathrm{\Phi }_2)`$. These integrals need to be evaluated for a grid of beam locations that has to contain at least $`𝒪\left(𝖫^3\right)`$ grid points to allow subsequent interpolation at arbitrary locations. Therefore the total computational cost for evaluating the convolution using Eq. (1) scales as $`𝒪\left(𝖫^5\right)`$ .

## III Factorizing the rotation

It is possible to simplify the evaluation of Eq. (1) significantly by factorizing the rotation into two auxiliary rotations such as

$$\widehat{D}(\mathrm{\Phi }_2,\mathrm{\Theta },\mathrm{\Phi }_1)\widehat{D}(\varphi _E,\theta _E,0)\widehat{D}(\varphi ,\theta ,\omega ).$$
(2)
We will define the various angles and motivate this split in the following. Figure 1 is intended to illustrate this discussion.

To introduce these coordinates let us first consider basic scan paths. Imagine a scan path where the beam sweeps over the sky by scanning on rings of angular radius $`\theta [0,\pi /2)`$. The centers of these scanning circles lie on a ring of constant latitude, a polar angle $`\theta _E[0,\pi )`$ away from the north pole. The angle $`\varphi _E[0,2\pi )`$ selects a given scanning circle and is defined as the longitude of its center, while $`\varphi [0,2\pi )`$ measures the angle along each scanning ring defined as increasing in a right-handed way about the outward normal at the center, starting from zero at the southernmost point on the ring. Hence, for such a basic scan path we can write the convolution as set of scalar products $`T(\varphi _E,\varphi )`$. The angles $`\theta `$ and $`\theta _E`$ are thought of as parameters which fixed to define the scanning geometry.

As a generalization of basic scan paths, we allow as a further degree of freedom an additional right handed rotation of the beam about its outward axis by an angle $`\omega [0,\pi /2)`$.

Now we can see geometrically that all beam orientations on generalized basic scan paths can be arrived at by successively applying the two rotations in Eq. (2). Define as the null position of the beam when it is oriented along the $`z`$-axis $`\theta =\theta _E=\varphi _E=\varphi =\omega =0`$. Starting from the null position, acting on it with $`\widehat{D}(\varphi ,\theta ,\omega )`$ rotates it about its axis by $`\omega `$ and moves it out onto a ring with opening angle $`\theta `$ at an azimuthal angle $`\varphi `$. Then acting with $`\widehat{D}(\varphi _E,\theta _E,0)`$ moves the beam into position.

Using the factorization Eq. (2) we can re-write Eq. (1) as

$`T(\varphi _E,\varphi ,\omega )={\displaystyle 𝑑\mathrm{\Omega }_\stackrel{}{\gamma }[\widehat{D}(\varphi _E,\theta _E,0)\widehat{D}(\varphi ,\theta ,\omega )b](\stackrel{}{\gamma })^{}s(\stackrel{}{\gamma })},`$ (3)
$`(\theta ,\theta _E\mathrm{fixed}).`$ (4)
The function $`T(\varphi _E,\varphi ,\omega )`$ contains all possible integrals for a given scanning geometry. In fact, for the special case $`\theta =\theta _E=\pi /2`$, these angles parameterize all possible orientations of the beam on the sky, i.e. $`(\varphi _E,\varphi ,\omega )`$ parameterize the group of rotations in three dimensions. It is well known that in this case these coordinates cover SO(3) twice, but this can be easily remedied by restricting the range of one of the angles to half its range. We defer removing this redundancy until the end of our calculation.

## IV Solution

To exploit the form of Eq. (4), it is expedient to represent the functions $`s`$ and $`b`$ as well as the rotation operators in the spherical harmonic basis. A bandlimited function $`f(\stackrel{}{\gamma })`$ can be expanded in spherical harmonics as

$$f(\stackrel{}{\gamma })=\underset{l=0}{\overset{l=𝖫_f}{}}\underset{m=l}{\overset{m=l}{}}f_{lm}Y_{lm}(\stackrel{}{\gamma }),$$
(5)
where $`\stackrel{}{\gamma }`$ denotes a unit vector. For practical applications the bandlimit $`𝖫_f`$ is chosen such that higher terms contribute insignificantly. We use the notation where all quantities carrying both an $`l`$ and an $`m`$ index vanish for $`m>l`$. This saves having to write explicit limits for sums over azimuthal quantum numbers.

Invariance of the scalar product under a change of basis then allows us to re-write Eq. (4) as

$$T(\varphi _E,\varphi ,\omega )=$$

$$\underset{lm}{}s_{lm}\left[\widehat{D}(\varphi _E,\theta _E,0)\widehat{D}(\varphi ,\theta ,\omega )b\right]_{lm}^{}=$$
(6)
$$\underset{lmMM^{}}{}s_{lm}D_{mM}^l(\varphi _E,\theta _E,0)D_{MM^{}}^l(\varphi ,\theta ,\omega )b_{lM^{}}^{}.$$

A simple explicit expression for the matrix elements $`D_{mm^{}}^l(\varphi _2,\theta ,\varphi _1)`$ can be given. One can define a real function $`d_{mm^{}}^l(\theta )`$ such that

$$D_{mm^{}}^l(\varphi _2,\theta ,\varphi _1)=e^{im\varphi _2}d_{mm^{}}^l(\theta )e^{im^{}\varphi _1}$$
(7)
Thus the dependence of $`D`$ on the Euler angles $`\varphi _1`$ and $`\varphi _2`$ is only in terms of complex exponentials. While explicit formulas for the $`d`$-functions exist , they are more conveniently computed numerically using their recursion properties .

Substituting into Eq. (6) and defining the three–dimensional Fourier transform of $`T(\varphi _E,\varphi ,\omega )`$ as

$$T_{mm^{}m^{\prime \prime }}=$$

$$\frac{1}{(2\pi )^3}_0^{2\pi }𝑑\varphi _E𝑑\varphi 𝑑\omega T(\varphi _E,\varphi ,\omega )e^{im\varphi _Eim^{}\varphi im^{\prime \prime }\omega }.$$
(8)
we obtain

$$T_{mm^{}m^{\prime \prime }}=\underset{l}{}s_{lm}d_{mM}^l(\theta _E)d_{MM^{}}^l(\theta )b_{lM^{}}^{}.$$
(9)
This equation is the main result of this paper, in effect generalising fast 2D Fourier transfrom convolution from the plane to the sphere. Its properties and specialisations will be discussed in the next section. Here we give a geometrical interpretation. We have arrived at Eq. (9) by writing convolution problems in such a way that the results are fields on 3-tori instead of subsets of the 3-sphere, which is the group manifold of rotations in 3 dimensions. Convolutions over azimuthally symmetric and connected sections of the 2-sphere (such as polar caps or annuli) can be parameterised by $`\theta `$ and $`\theta _E`$ and hence can be extended to 3-tori as shown. Since exponentials are a complete and orthonormal basis on the 3-torus and because we assumed that $`s`$ and $`b`$ are band-limited, the $`T_{mm^{}m^{\prime \prime }}`$ contain all information about the inverse transform, Eq. (4),

$$T(\varphi _E,\varphi ,\omega )=$$

$$\underset{m,m^{},m^{\prime \prime }=𝖫}{\overset{𝖫}{}}T_{mm^{}m^{\prime \prime }}e^{im\varphi _E+im^{}\varphi +im^{\prime \prime }\omega }.$$
(10)
Not all tuples $`(\varphi _E,\varphi ,\omega )`$ correspond to distinct beam orientations but this redundancy is more than compensated for by the efficiency of the method.

## V Discussion

Several remarks about Eq. (9) are in order.

#### 1 Computational cost

Computing the $`T_{mm^{}m^{\prime \prime }}`$ in Eq. (9) costs $`𝒪\left(𝖫^4|\mathrm{sin}\theta _E|2\theta /\pi \right)`$ operations. The factors $`\theta `$ and $`|\mathrm{sin}\theta _E|`$ come from the fact that the bandlimit $`𝖫`$ for the sky implies a bandlimit $`2\theta 𝖫/\pi `$ on a ring of radius $`\theta `$ and hence the ranges of $`m`$ and $`m^{}`$ can be reduced by factors of $`|\mathrm{sin}\theta _E|`$ and $`2\theta /\pi `$, respectively, if the rings and the ring of ring centers are not great circles. Using the Fast Fourier Transform (FFT) algorithm, the inverse Fourier transform takes $`𝒪\left(𝖫^3\mathrm{log}𝖫|\mathrm{sin}\theta _E|2\theta /\pi \right)`$ operations. If the convolved sky is assumed to be real, we have

$$T_{mm^{}m^{\prime \prime }}=T_{mm^{}m^{\prime \prime }}^{},$$
(11)
reducing memory and processor requirements by a factor 2.

#### 2 Quadrature and interpolation

In pixel space each evaluation of $`T(\varphi _E,\varphi ,\omega )`$ is an explicit quadrature problem and hence necessarily approximate. In our approach, all sums have a finite number of terms and the results are exact as long as $`s`$ and $`b`$ are band-limited. Quadrature issues only have to be dealt with if $`b`$ or $`s`$ are given in pixel space and we have to evaluate the beam and sky multipole coefficients $`b_{lm}`$ and $`s_{lm}`$. The details of which pixelization to choose on the sphere and how to solve this generalized quadrature problem for the multipole coefficients are outside of the scope of this work but an efficient and practical approach to the quadrature problem is implemented in the HEALPix package  and will be discussed in a future publication .

An interesting property of Eq. (9) is that as long as $`𝖫`$ was chosen appropriately one is guaranteed to have the convolved sky sampled sufficiently densely for worry-free interpolation on either of the three indices.

### A Special cases

We will now discuss certain special cases of Eq. (9).

#### 1 Total convolution

Let us obtain the convolved sky at all possible beam orientations $`\omega `$ on an equidistant coordinate grid in $`\varphi `$ (corresponding to the polar angle) and $`\varphi _E`$ (corresponding to the azimuthal angle). We will refer to this case as the total convolution. This can be achieved by evaluating Eq. (9) setting $`\theta =\theta _E=\frac{\pi }{2}`$. In this case we only need to know $`d_{m^{}m}^l(\frac{\pi }{2})`$. This means we only have to evaluate a single recursion relation to evaluate the sum on $`l`$, which simplifies the computation. The inverse FFT gives the desired result.

A further simplification arises in this case from the fact that if $`\theta _E=\frac{\pi }{2}`$, not all components of $`T_{mm^{}m^{\prime \prime }}`$ are independent. The redundancy in the parametrization where the polar angle $`\varphi [0,2\pi )`$ leads to the symmetry

$$T(\varphi _E,\varphi ,\omega )T(\pi +\varphi _E,2\pi \varphi ,\pi +\omega ).$$
(12)
This translates into the identity

$$T_{mm^{}m^{\prime \prime }}(1)^{m+m^{\prime \prime }}T_{mm^{}m^{\prime \prime }},$$
(13)
which cuts the required memory and computation time by a factor 2.

#### 2 Exact or approximate azimuthal symmetry of the beam

In many practical situations the “beam” represents the response function of an optical system with only mild imperfections. If this is the case, the beam has only slowly varying azimuthal structure, implying a cutoff wavenumber $`𝖬`$ such that $`b_{lm}0`$ for $`m𝖬`$. in this case the computational cost for a total convolution scales as $`𝒪\left(𝖫^3𝖬\theta \mathrm{sin}\theta _E\right)`$ .

In the limit of an azimuthally symmetric beam, $`𝖬=0`$, we obtain an $`𝒪\left(𝖫^3\theta \mathrm{sin}\theta _E\right)`$ method. However, it is known  that at least in principle there exist faster methods for convolution of a function on the full sphere ($`\theta =\theta _E=\pi /2`$) with an azimuthally symmetric beam which scale as $`𝒪\left(𝖫^2(\mathrm{log}(𝖫))^2\right)`$ . We can show how this limit is obtained from Eq. (9) by using the facts that the $`d_{mm^{}}^l(\frac{\pi }{2})`$ are the Fourier coefficients of the $`d_{mm^{}}^l(\theta )`$ and that $`d_{m0}^l(\theta )=P_{lm}(\theta )`$. Then Eq. (9) reduces to the form

$$T(\varphi _E,\varphi )=\underset{lm}{}Y_{lm}(\pi \varphi ,\varphi _E+\pi /2)b_{l0}s_{lm},$$
(14)
where the arguments of $`Y_{lm}`$ are the polar angle and the azimuthal angle, respectively. The algorithm by  succeeds precisely in reducing the computational cost of evaluating this expression to $`𝒪\left(𝖫^2(\mathrm{log}^2𝖫)\right)`$ under the proviso of the technical difficulties there outlined.

We note here for completeness, that by choosing a delta function beam (and hence $`b_{lm}=const`$), we recover the Fourier summation method for the spherical harmonic transform, described in equations (5.2) to (5.4) in . This computes Eq. (5) on an equidistant coordinate grid by doing Fourier transforms on latitudinal and longitudinal lines. The forward transform is obtained by simply working all steps in reverse.

#### 3 Basic scan paths

Consider an application where the convolution is required only along a “basic scan path”. This is one of the proposed scan strategies for the Planck satellite mission. From our definition of basic scan paths in section III we see that they correspond to setting $`\omega =0`$ in Eq. (9).

Computing the inverse Fourier transform of Eq. (9) with $`\omega =0`$ just amounts to summing over $`m^{\prime \prime }`$. Then only the two–dimensional Fourier transform

$$T_{mm^{}}(\omega =0)=\underset{l}{}s_{lm}d_{mm^{}}^l(\theta _E)X_{lm^{}}$$
(15)
remains to be evaluated. The quantity

$$X_{lm}\underset{M}{}d_{mM}^l(\theta )b_{lM}^{}$$
(16)
can be precomputed. All in all the computational time needed for evaluating these expressions is $`𝒪(𝖫^3\theta \mathrm{sin}\theta _E)`$. Storage requirements scale only as $`𝒪\left(𝖫^2\right)`$ .

Note that in this case azimuthal symmetry does not necessarily imply reduced computational cost. If the beam is concentrated at the north pole into a region of size $`\sigma `$ then $`X_{lm}`$ will have $`m`$ modes populated up to $`𝖬\theta /\sigma `$.

Geometrically, the basic scan path corresponds to a 2-torus which is the section of the 3-sphere of rotations at constant $`\omega `$. Note that in this case there is no redundancy in the parametrization — every tuple $`(\varphi _E,\varphi ,0)`$ corresponds to a distinct beam orientation.

#### 4 Perturbations about basic scan paths

A slight generalization of the previous case are scan paths which are close to basic but include a variation in $`\theta _E`$ from scanning circle to scanning circle. Such paths result for example from precessing or “wobbling” the spin axis of a scanning satellite.

Such scan paths can be composed by computing several convolutions along basic scan paths for different angles $`\theta _E`$ and then choosing scanning circles at will from among these basic ones. This method suggests itself if the precession angle is small and hence a small number of convolutions is sufficient to sample the variation in $`\theta _E`$. Convolutions at points which do not coincide with sampling points can then be determined by interpolation.

Another approach to this type of problem and further generalisations are discussed in the next paragraphs.

#### 5 Other special cases

Other potentially interesting special cases of Eq. (9) can be worked out by fixing any of the parameters to special values and evaluating the inverse transform, analogous to the calculation for basic scan paths. For example one obtains expressions for

* All possible beam orientations along a circle of constant latitude $`\theta _E`$. In this case $`\varphi `$ and $`\omega `$ have the same meaning; formally, $`T_{mm^{}m^{\prime \prime }}=T_{mm^{}}\delta _{m^{}m^{\prime \prime }}`$ and we obtain an $`𝒪\left(𝖫^2𝖬\mathrm{sin}\theta _E\right)`$ method:

$$T_{mm^{}}=\underset{l}{}s_{lm}d_{mm^{}}^l(\theta _E)b_{lm^{}}^{}.$$
(17)
* Individual scanning rings of a basic scan path. Here, $`\omega =0`$, and the only free parameter is $`\varphi `$.

The details of the calculations for this and similar cases are now easy exercises.

#### 6 Generalizations

Further, it is clear from the derivation that more general types of paths can be constructed by factorizing the rotation operator more than twice, so as to generate for example a ring of ring of rings, etc. For particular applications some of these may be advantageous, for example if they simplify the interpolation problem on the output ring set. A specific example is the precessing scan path mentioned in the previous subsection. Inserting another rotation operator $`\widehat{D}(0,\theta _P,\varphi _P)`$ between the two operators in Eq. (2), and setting $`\omega =0`$ produces a set of rings whose centers lie on circles of radius $`\theta _P`$ about $`theta_E`$. This may simplify the interpolation problem. The rotation corresponding to $`\varphi _P`$ can be sampled sparsely if $`\theta _P`$ is small ($`𝖬_P\mathrm{sin}\theta _P\mathrm{sin}\theta _E𝖫`$ with obvious notation) and the interpolation problem becomes simpler.

## VI Conclusions

This paper presents a general algorithm which greatly reduces the computational cost of convolving two bandlimited but otherwise arbitrary functions on the sphere. The speedup increases linearly with the smallest angular scale of the smoother of the two functions in the problem. The scalings of the necessary operation counts are discussed in detail in section five.

We quote in the appendix formulas showing how the ideas presented in this paper can be applied to the inverse problem of “deconvolution” by speeding up the iterative solution of the normal equation in an analogous way.

This paper focuses on the convolution of scalar valued functions on the sphere such as temperature, elevation, etc. In order to be able to deal with the polarization of the cosmic microwave background we extend the methods presented here to tensor valued functions on the sphere in reference .

The algorithms which are presented here are already being used as a core component of the prototype simulation pipeline of the Planck satellite. To give an example for the timing gains one makes by applying this method, we computed the following case: both sky and beam were interpolated and pixelised very densely, with millions of pixels each to resolve the steep variations over many orders of magnitude. The bandlimit was somewhat generously chosen as $`𝖫=1024`$. Then the convolution of the sky with the beam of a single detector for a whole year of mission data, consisting of $`(2049)^24\times 10^6`$ convolved samples along a basic scan path was generated in less than 15 minutes on a single Silicon Graphics R10000 processor. This compares with several days of computation on a severely coarsened sampling grid with several hundred times fewer samples on the same machine, using the adaptive mesh method . For the same resolution which we achieved with our methods, the adaptive mesh code would have run for months.

Due to our methods, future CMB missions can go beyond having to approximate the treatment of realistic beams. Our methods lend themselves to being used in conjunction with iterative map–making methods to remove from the data artefacts which are due to beam distortions and far side lobes (see appendix).

Lastly, we feel that the geometric constructions, analogies to group properties and algebraic results we introduce in this article may be useful more generally for CMB data analysis and plan to explore these issues in future work.

###### Acknowledgements.

We would like to thank F. Hansen for coding help in the early stages of this project, F. Bouchet for suggesting consideration of the transpose operation and R. Caldwell for reading the manuscript before submission. Part of this work was supported by the Dansk Grundforskningsfond through its funding for TAC. BDW is supported by the NASA MAP/MIDEX program.

## A The convolution transpose

We sketch here how to set up the inverse problem of reconstructing the true sky from the convolved observations. If we start with a noise free set of convolutions, then the equation to be inverted in order to estimate the true underlying sky is, schematically,

$$𝐀𝐬=𝐝.$$
(A1)
Here, $`𝐬`$ is the true sky, $`𝐝`$ is the vector containing the time-ordered data after observation. The convolution operator is represented by $`𝐀`$.

The least-square estimator for the true sky, $`\widehat{𝐬}`$ then satisfies the normal equation

$$𝐀^𝐓𝐀\widehat{𝐬}=𝐀^𝐓𝐝.$$
(A2)
For a perfect observation with a delta function beam, $`𝐀^𝐓𝐀1\mathrm{I}`$. So it may be reasonable to expect that we can make progress by considering a mildly imperfect optical system and consider iterative techniques for solving the normal equation. In this case the ability to solve for $`\widehat{s}`$ iteratively (e.g. using a Conjugate Gradient technique) relies on convergence (which is assured up to numerical effects because the normal matrix $`𝐀^𝐓𝐀`$ is positive definite and being able to compute the matrix products in (A2) quickly. The application of $`𝐀`$ can be computed efficiently using the formulae set out in sections four and five. We now find an algorithm for the efficient application of $`𝐀^𝐓`$, the transpose convolution.

### 1 Applying the transpose convolution

We can write down the expression for $`𝐀^𝐓`$ in a similar way to Eq. (4)

$$y(\stackrel{}{\gamma })=$$

$$𝑑\varphi _E𝑑\varphi 𝑑\omega [\widehat{D}(\varphi _E,\theta _E,0)\widehat{D}(\varphi ,\theta ,\omega )b](\stackrel{}{\gamma })^{}T(\varphi _E,\varphi ,\omega ).$$
(A3)
Now the derivation is analogous to the one preceding Eq. (9), and $`y(\stackrel{}{\gamma })`$ is given in terms of the Wigner $`d`$ functions as

$$y(\stackrel{}{\gamma })=\underset{lmm^{}m^{\prime \prime }}{}Y_{lm}^{}(\stackrel{}{\gamma })d_{mm^{}}^l(\theta _E)d_{m^{}m^{\prime \prime }}^l(\theta _E)b_{lm^{\prime \prime }}^{}T_{mm^{}m^{\prime \prime }}$$
(A4)
This formula can be generalised or applied to special cases just as we showed in section V for Eq. (9).


format: text

text:
# Discovery of a Close Pair of 𝑧 = 4.25 Quasars from the Sloan Digital Sky Survey 1footnote 11footnote 1Based on observations obtained with the Sloan Digital Sky Survey, which is owned and operated by the Astrophysical Research Consortium.,2footnote 22footnote 2Based on observations obtained with the Hobby-Eberly Telescope, which is a joint project of the University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximillians-Universität München, and Georg-August-Universität Göttingen.

## 1 Introduction

One of the first observations of the $`z>4`$ universe was also one of the most startling: the serendipitous discovery of the $`z=4.4`$ quasar Q 2203+29 by McCarthy et al. (1988). While obtaining a spectrum of the $`z=0.707`$ radio galaxy 3C 441, the spectrograph slit was aligned along the radio jet. The data included the spectrum of a faint red source 51′′ from the nucleus; the secondary source’s redshift of 4.41 made it the second most distant known object at the time! The authors showed that this was a very low probability event unless the number density of high-redshift quasars had been underestimated by orders of magnitude (a slight possibility at the time given that only a few $`z>4`$ quasars were known and thus the number density poorly constrained).

Several years later PC 0027+0521, a faint $`z=4.21`$ quasar, was found by chance in a high-redshift quasar survey (Schneider, Schmidt, & Gunn 1994). The object happened to lie in the slit ($`i.e.,`$ the slit was placed at the parallactic angle, not rotated to acquire additional objects) while the spectrum of a high-redshift quasar candidate (which turned out to be a star!) was obtained.

Post facto calculations showed that the chance of either event occurring was on the order of 0.1%. In this paper we describe the third serendipitous discovery of a $`z>4`$ quasar; it was discovered, like PC 0027+0521, when by chance it was included in a spectroscopic observation of a high-redshift quasar candidate. In this case, however, the candidate was indeed a quasar, and it has the same redshift as the serendipitous source.

## 2 Observations

### 2.1 Sloan Digital Sky Survey

The Sloan Digital Sky Survey (SDSS; York et al. (2000)) uses a CCD camera (Gunn et al. (1998)) on a dedicated 2.5-m telescope (Siegmund et al. (2001)) at Apache Point Observatory, New Mexico, to obtain images in five broad optical bands over 10,000 deg2 of the high Galactic latitude sky centered approximately on the North Galactic Pole. The five filters (designated $`u^{}`$, $`g^{}`$, $`r^{}`$, $`i^{}`$, and $`z^{}`$) cover the entire wavelength range of the CCD response (Fukugita et al. (1996)). Photometric calibration is provided by simultaneous observations with a 20-inch telescope at the same site. The survey data processing software measures the properties of each detected object in the imaging data, and determines and applies astrometric and photometric calibrations (Pier et al. (2001); Lupton et al. (2001)).

The high photometric accuracy of the SDSS images and the information provided by the $`z^{}`$ filter (central wavelength of 9130 Å) make the SDSS data an excellent source for identification of high-redshift quasar candidates. In the past two years the SDSS has identified many $`z>3.5`$ quasars (Fan et al. 1999a, 2000a,c; Schneider et al. 2000; Zheng et al. 2000), culminating with the discovery of a $`z=5.80`$ quasar (Fan et al. 2000b ).

We have started a survey of faint, high-redshift quasars using the SDSS imaging data. This survey uses a multicolor selection technique similar to that of Fan et al. (1999a,2000a); the major modification is a change from $`i^{}`$ = 20.0 to $`i^{}=21.0`$ for the limiting magnitude. The $`i^{}20.6`$ object SDSSp J143952.58$``$003359.2 was flagged as a probable $`z4.3`$ quasar, primarily based on its location in the ($`g^{}r^{}`$), ($`r^{}i^{}`$) diagram, from imaging data taken on 22 March 1999 (SDSS imaging run 756). For brevity, we will sometimes refer to this object as “A” in this paper.

The area containing the quasar had also been observed on 20 March 1999 (SDSS imaging run 745). The SDSS photometric measurements of the quasar on both nights are presented in Table 1, which shows that the photometry/errors are internally consistent.

Notes on SDSS nomenclature: The names for sources have the format SDSSp Jhhmmss.ss+ddmmss.s, where the coordinate equinox is J2000, and the “p” refers to preliminary. The reported magnitudes are based on a preliminary photometric calibration; to indicate this the filters have an asterisk instead of a prime superscript. The estimated astrometric accuracies in each coordinate are 0.1′′ and the calibration of the photometric measurements are accurate to 0.05 magnitudes. Note that the $`u^{}`$ measurements indicate essentially zero flux in this band.

### 2.2 Spectroscopy of Quasar Candidates

A spectrum of SDSSp J143952.58$``$003359.2 was obtained on 2 June 2000 with the Marcario Low Resolution Spectrograph (LRS; Hill et al. 1998a,1998b) of the Hobby-Eberly Telescope (HET; Ramsey et al. (1998)) 333A summary of the commissioning of the LRS is given in Schneider et al. (2000).. The LRS is mounted in the Prime Focus Instrument Package, which rides on the HET tracker. The 900 s observation was acquired with the 2′′ wide slit, rotated to the parallactic angle. The dispersive element was a 300 line mm-1 grism blazed at 5500 Å. An OG515 blocking filter was installed to permit calibration of the spectra out to 1 micron. The detector is a thinned, antireflection-coated 3072 $`\times `$ 1024 Ford Aerospace CCD, and was binned $`2\times 2`$ during readout; this produced an image scale of 0.50′′ pixel-1 and a dispersion of $``$ 4.5 Å pixel-1. The spectra covered the range from 5100–10,200 Å at a resolution of approximately 20 Å.

The wavelength calibration was provided by Ne, Cd, and Ar comparison lamps, and the relative flux and atmospheric band absorption calibration was provided by the observation of the spectrophotometric standard BD +26 2606 (Oke & Gunn 1983). The observing conditions were marginal; the seeing was 2′′ with moderate cloud cover (the signal from the standard indicated more than a magnitude of extinction). Absolute spectrophotometric calibration was carried out by scaling each spectrum so that $`i^{}`$ magnitudes synthesized from the spectra matched the SDSS photometric measurements.

Figure 1 shows the SDSS $`i^{}`$ image of the field of SDSSp J143952.58$``$003359.2, which is denoted by an “A” in the figure; this object was placed in the center of the slit. The slit was oriented at a position angle of 26; the spectra of the bright star and object “B”, both located southwest of A, were obtained as they fell into the slit. A single 900 s exposure was made.

Figure 2 displays the region centered at 6680 Å of the spectroscopic exposure. The spectra run vertically from 5790 Å (bottom) to 7575 Å (top). The data are dominated by the spectrum of the star slightly to the right of center. The spectrum of A is clearly visible in the center of the frame, and the spectrum possesses the characteristic signal of a high-redshift quasar (strong, broad emission line, and a significant drop in the continuum as one moves from red to blue across the line). The initial visual inspection of the frame revealed that the spectrum of B, located midway between the star and the right of the frame in Figure 2, appeared to have a very similar spectrum to A.

The flux and wavelength calibrated spectra of the two objects are presented in Figure 3 (the brighter object, A, is the top panel). Although the signal-to-noise ratio is not large (6 per binned pixel at 7000 Å in A, a factor of two lower in B), there is little doubt that both objects are $`z4.25`$ quasars: the broad, strong, asymmetric Lyman $`\alpha `$ emission line and the continuum depression due to Lyman $`\alpha `$ forest absorption are the most obvious features. The Lyman $`\alpha `$ forest absorption appears to be significantly different in the two quasars; both the photometry and the spectroscopy indicate that the spectrum of A suffers considerably less absorption than B. We are unable to search for any Lyman limit absorption because the OG515 filter does not transmit light below a wavelength of $``$ 980 Å in the rest frame of the quasars.

The redshifts, determined from the Si+O and C IV emission lines in A and the C IV line in B, are identical to within the measurement errors. The two objects, however, are almost certainly a pair of quasars and not a lens, given the relatively large separation and markedly different profiles of the Lyman $`\alpha `$ emission line and the contrast in the strengths of the Si+O feature.

## 3 Discussion

Object “B” was identified with SDSSp J143951.60$``$003429.2 in the SDSS database; its photometric measurements are presented in Table 1. Figure 4 shows the locations of B and A in the $`g^{}r^{}`$,$`r^{}i^{}`$ diagram, which is the SDSS color-color plot in which $`z4`$ are most clearly separated from the stellar locus (the Lyman $`\alpha `$ emission line falls in $`r^{}`$, the Lyman forest occurs in $`g^{}`$; see Fan 1999). Object A satisfies the high-redshift selection criteria in this diagram (shaded area) and is approximately 0.4 magnitudes brighter than the $`i^{}=21.0`$ limit. Object B lies in the same direction as A from the stellar locus, but is slightly outside the selection region and is also nearly half a magnitude too faint in $`i^{}`$ to be included in the sample. For both objects the $`g^{}`$ measurement has a very low signal-to-noise ratio.

A summary of the properties of the two quasars is given in Table 2; the absolute magnitudes were calculated assuming $`H_0`$ = 50 km s-1 Mpc-1 and $`q_0`$ = 0.5. (We will use this cosmological model throughout this section.) We do not quote a continuum depression value ($`D_A`$) because of the limited signal-to-noise ratio of the data. The quasars are relatively faint (indeed, the fainter of the pair is quite similar in redshift and luminosity to the second serendipitous $`z>4`$ object, PC 0027+0521), but the luminosity of each is well above the canonical “quasar/Seyfert” dividing line of $`M_B=23`$. The observed radio emission from each quasar must be less than 1 mJy based on the lack of detection in the FIRST (Becker, White & Helfand 1995) survey data; the quasars are radio-quiet, as are the vast majority of $`z>4`$ quasars (Schneider et al. 1992; Schmidt et al. 1995; Stern et al. 2000; Carilli et al. 2000).

The quasars are separated by 33.4′′ and B is at a position angle of 206 relative to A. In our adopted cosmology, this corresponds to a projected spatial separation of 208 kpc (1.1 Mpc comoving), much closer than any other published pair of $`z>4`$ quasars. We have measured the redshift difference in the spectra both from individual line redshifts and cross-correlating the spectra; both techniques yield a redshift difference consistent with zero but with an uncertainty of 0.01.

If the redshift difference is 0.01 and this is due to cosmological expansion rather than peculiar velocities, then the radial separation is 950 kpc (5 Mpc comoving). It is quite possible, however, that the redshift difference is much smaller, in which case assigning a distance to the redshift difference is suspect. If we assume that we have no radial information and are viewing the pair at a random angle, the mean deprojected separation is 327 kpc (1.7 Mpc comoving).

The only other known pair of $`z>4`$ quasars with a comoving separation of less than 100 Mpc was found in observations for the Palomar Scan Grism Survey (PSGS; Schneider, Schmidt, & Gunn 1999): PC 0027+0525 ($`z=4.099`$, $`M_B=24.8`$) and the previously mentioned serendipitously discovered PC 0027+0521 ($`z=4.210`$, $`M_B=24.0`$). This pair is separated on the sky by 319′′ and has a comoving separation of approximately 58 Mpc.

Both these pairs are chance discoveries, so calculation of the probability of finding them will necessarily involve post facto arguments. With this caveat in mind, we investigate this probability, first assuming that quasars are unclustered, then taking clustering into account. What are the chances of finding pairs of quasars of a given spatial separation between redshifts of 4.0 and 4.5 where both members have $`M_B<24.0`$? In the following discussion, we will examine three examples: a “minimum” separation model (330 kpc; 1.7 Mpc comoving) and a “maximum” separation model (1 Mpc; 5.3 Mpc comoving) for the SDSS pair, and a 58 Mpc separation for the PSGS objects.

Fan et al. (2000d) have investigated the luminosity function of high-redshift quasars from a sample of 39 $`z>3.6`$ objects identified by the SDSS (18 quasars with redshifts larger than four). They find the comoving number density of quasars more luminous than $`M_B26`$ at a redshift of 4.25 to be $`2.8\times 10^8`$ Mpc-3, and an $`nL^{2.5}`$ differential luminosity function. This numbers are in excellent agreement with the previous studies by Schmidt, Schneider, & Gunn (1995) and Kennefick, Djorgovski, & de Carvalho (1995), which were based on samples that contained 9 and 10 $`z>4`$ quasars, respectively. This luminosity function produces a number density of $`4.4\times 10^7`$ Mpc-3 for quasars with $`M_B<24.0`$. This calculation requires an extrapolation of more than 1.5 magnitudes from the Schmidt et al. and Fan et al. data, and should be viewed as an upper limit on the number density given the evidence ($`e.g.,`$ the PSGS) that the quasar luminosity function flattens at luminosities below $`M_B=26`$.

For the three examples mentioned above, the volumes enclosed by spheres whose radii are equal to the separations are 20.6 Mpc3, 624 Mpc3, and $`8.2\times 10^5`$ Mpc3; the chances of more than one quasar occurring in each volume are $`4.1\times 10^{11}`$, $`3.9\times 10^8`$, and 0.051. The volume enclosed by the shell bounded by redshifts of 4.0 and 4.5 is $`1.4\times 10^{11}`$ Mpc3 so one would expect 0.28, 8.9, and 8700 such pairs over the entire sky. Given that the PSGS/SDSS have examined only 1% of the celestial sphere for high-redshift quasars, the likelihood of finding a 1.7 Mpc pair is approximately one part in 350 and there is only a nine percent chance of identifying a 5.3 Mpc pair. The existence of only two, not 87, known pairs with smaller than 58 Mpc separation is merely a reflection of the fact that the current sample of $`z>4`$ quasars is woefully incomplete to $`M_B=24`$, which comes as no surprise. This implied incompleteness, much more than a factor of ten, makes the discovery of the SDSS pair even more striking.

This calculation suggests that given the existence of the SDSS pair, it is very unlikely that high-redshift quasars are not clustered. Kundić (1997) and Stephens et al. (1997) found that luminous $`z>2.7`$ quasars are clustered based on the three pairs with less than 30 Mpc (comoving) separation found in the 90 high-redshift quasars of the Palomar Grism Transit Survey (Schneider, Schmidt, and Gunn 1994). Their estimate of the comoving correlation length was 35 Mpc, with a relatively large uncertainty. Bahcall & Chokshi (1991) discussed the origin of low-redshift quasar correlations, and two recent papers (Martini & Weinberg 2000; Haiman & Hui 2000) predict that high-redshift quasars should be clustered based on theoretical arguments.

How strong must the quasar-quasar correlation function be to produce the observed SDSS pair in only 1% of the sky? If we assume a correlation function of the form seen in the local universe and that for $`z3`$ galaxies (Giavalisco et al. 1998)

$$\xi (r)=\left(\frac{r}{r_0}\right)^\gamma $$

the ratio of the number of objects in the clustered case to the unclustered case within a sphere of radius $`R`$ is

$$\frac{n_{\mathrm{corr}}}{n_{\mathrm{field}}}|_R=\frac{_0^R4\pi r^2\xi (r)𝑑r}{_0^R4\pi r^2𝑑r}=\frac{3}{3\gamma }\left(\frac{R}{r_0}\right)^\gamma $$

Adopting a representative value of $`\gamma `$ =$`1.8`$, this expression becomes

$$\frac{n_{\mathrm{corr}}}{n_{\mathrm{field}}}|_R=2.5\left(\frac{R}{r_0}\right)^{1.8}$$

If this expression describes the high-redshift quasar-quasar correlation function, then the scale length $`r_0`$ must approach 30 Mpc for the minimum model and 12 Mpc for the maximum model for there to be a reasonable chance of finding a pair in the SDSS database to date. (Remember that this calculation assumes that the quasar sample is complete to $`M_B=24`$!)

One could argue that the above discussion is invalid because current surveys are failing to find most of the quasars at this redshift due to effects such as reddening ($`e.g.,`$ Ostriker & Heisler 1984); therefore the estimate of the mean number density of quasars is grossly underestimated and hence “chance” pairs much more prevalent the calculation suggests. This may indeed be the state of affairs, but even if true would not invalidate the above reasoning. The quasar surveys mentioned above have detection algorithms that rely on the strength of the Lyman $`\alpha `$ absorption and strong Lyman $`\alpha `$ emission (although the SDSS has demonstrated the ability to find high-redshift objects without Lyman $`\alpha `$ emission; see Fan et al. 1999b). Perhaps we can identify but a tiny fraction of the luminous objects at redshift greater than four, but we can state that luminous objects with strong Lyman $`\alpha `$ emission are clearly strongly clustered.

These numbers indicate that it is highly likely that a significant number of $`z>4`$ quasar pairs will be discovered in current large area surveys. Recent investigations of the environments of high-redshift quasars ($`e.g.,`$ Djorgovski 1999) demonstrate that these distant beacons do not occur in isolation. The SDSS spectroscopic survey should discover several thousand quasars with redshifts larger than four, and the SDSS image database will contain useful color information for objects brighter than $`i^{}21.5`$, similar to the brightness of the fainter of the SDSS pair in this study, in all of the fields. While the efficiency of identification of high-redshift quasars in the primary SDSS survey at this limit will fall well below the 50%–70% levels seen in the SDSS to date ($`e.g.,`$ Fan et al. 1999a,2000a), the number of such candidates within an arcminute or two of known high-redshift quasars should be sufficiently small for spectroscopic investigation. The southern SDSS survey, consisting of an area of 225 sq deg that will reach a magnitude or more deeper than the primary SDSS survey, should be ideal for studies of the high-redshift quasar correlation function.

We would like to thank Gary Hill, Marsha Wolf, Grant Hill, and Matt Shetrone for aid with the HET data acquisition, and John Booth and his team for their efforts that have markedly improved the HET’s image quality during the month preceding the observation. The Sloan Digital Sky Survey (SDSS) is a joint project of the University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, the Johns Hopkins University, the Max-Planck-Institute for Astronomy, Princeton University, the United States Naval Observatory, and the University of Washington. Apache Point Observatory, site of the SDSS, is operated by the Astrophysical Research Consortium. Funding for the project has been provided by the Alfred P. Sloan Foundation, the SDSS member institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, and Monbusho. The SDSS Web site is http://www.sdss.org/. The Hobby-Eberly Telescope (HET) is a joint project of the University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximillians-Universität München, and Georg-August-Universität Göttingen. The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly. The LRS is named for Mike Marcario of High Lonesome Optics who fabricated several optics for the instrument but died before its completion. This work was supported in part by National Science Foundation grants AST99-00703 (DPS) and AST96-18503 (MAS). MAS and XF acknowledge additional support from Research Corporation, the Princeton University Research Board, and a Porter O. Jacobus Fellowship.

| TABLE 1. SDSS Photometrya | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
|  | SDSS |  |  |  |  |  |
| Objectb | Run | $`u^{}`$ | $`g^{}`$ | $`r^{}`$ | $`i^{}`$ | $`z^{}`$ |
| SDSSp J143952.58$``$003359.2 | 745 | 23.17 | 23.41 | 20.98 | 20.56 | 20.20 |
|  |  | $`\pm `$0.51 | $`\pm `$0.24 | $`\pm `$0.05 | $`\pm `$0.06 | $`\pm `$0.16 |
| SDSSp J143952.58$``$003359.2 | 756 | 24.20 | 23.35 | 20.95 | 20.57 | 20.27 |
|  |  | $`\pm `$0.59 | $`\pm `$0.20 | $`\pm `$0.04 | $`\pm `$0.05 | $`\pm `$0.12 |
| SDSSp J143951.60$``$003429.2 | 745 | 23.58 | 23.11 | 21.69 | 21.67 | 21.39 |
|  |  | $`\pm `$0.61 | $`\pm `$0.19 | $`\pm `$0.09 | $`\pm `$0.15 | $`\pm `$0.42 |
| SDSSp J143951.60$``$003429.2 | 756 | 22.95 | 23.48 | 21.89 | 21.46 | 20.95 |
|  |  | $`\pm `$0.38 | $`\pm `$0.22 | $`\pm `$0.10 | $`\pm `$0.10 | $`\pm `$0.21 |

a Photometry is reported in terms of asinh magnitudes; see Lupton, Gunn, & Szalay (1999) for details. In this system, zero flux corresponds to 23.24, 24.91, 24.53, 23.89, and 22.47 in $`u^{}`$, $`g^{}`$, $`r^{}`$, $`i^{}`$, and $`z^{}`$, respectively.

b Coordinate equinox is J2000.

| TABLE 2. Properties of Quasarsa | | | |
| --- | --- | --- | --- |
| Object | $`z`$ | $`AB_{1450}`$ | $`M_B^a`$ |
| SDSSp J143952.58$``$003359.2 | 4.255 $`\pm `$ 0.010 | 20.67 | $`24.4`$ |
| SDSSp J143951.60$``$003429.2 | 4.258 $`\pm `$ 0.010 | 22.03 | $`25.8`$ |

a Calculated assuming $`H_0`$ = 50, $`q_0`$ = 0.5, and $`\alpha `$ = $`0.5`$, and standard Galactic reddening law with $`E(BV)`$ = 0.039 (Schlegel, Finkbeiner, & Davis (1998)).


text:
# Magnetotransport measurements on freely suspended two-dimensional electron gases

## Abstract

We present magnetotransport measurements on freely suspended two-dimensional electron gases from AlGaAs/GaAs heterostructures. The technique to realize such devices relies on a specially MBE grown GaAs/AlGaAs-heterostructure including a sacrificial layer. We fabricated simple Hall-bars as well as quantum cavities and quantum dot systems. We find well-pronounced Shubnikov-de Haas oscillations and observe commensurability resonances, allowing characterization of the electron gas in these 100 nm thin membranes.

Recently the realization of suspended, monocrystalline, GaAs nanostructures containing a three-dimensional electron gas has been demonstrated . Previous work has shown the possibility of using a GaAs high-electron mobility transistor to achieve sensitive piezoelectric detection of strain . We are combining these techniques to enable new means of detecting motion in nanomechanical systems and to study interactions in coupled electron-phonon systems of reduced dimensionality for both the electrons and the phonons. Free-standing structures incorporating a high-mobility electron gas have interesting device applications: A high mobility 2DEG system provides a unique approach to implement wideband, extremely sensitive displacement detection. These systems also offer prospects for very sensitive bolometers and represent model systems for high sensitivity calorimetry . One very interesting potential application is in heat capacity measurements on two-dimensional electron gases (2DEGs). The sensitivity of these types of thermal devices is enhanced by their small volume, which can be of order of $`3\mu `$m3, significantly smaller than the usual dimensions of $`10^8\mu `$m3 .

In this work we discuss the processing technique of such devices and present first results on magnetotransport measurements on suspended two- and zero-dimensional systems. The processing employs a specially-designed MBE-grown 2DEG heterostructure in which a sacrificial layer is included. The layer structure is shown in Fig. 1(a): The structural layer stack, from which the devices are formed, comprises seven individual layers having a total thickness of 100 nm. Top and bottom are formed by thin GaAs cap layers preventing oxidation of the AlGaAs:Si donor layers which follow beneath. The central 15 nm thick GaAs layer forms a quantum well sustaining a high mobility 2DEG located 37 nm below the top surface and is surrounded by two AlGaAs spacer layers. Below the structural layer stack is a 400 nm Al0.8Ga0.2As sacrificial layer.

The heterostructure was designed by modelling the conduction band lineup numerically. These calculations employed a Poisson-Schrödinger solver program, written by Laux and Kumar, which neglects many-body effects . In contrast to common 2DEG heterostructure configurations, here an additional GaAs layer is included to avoid deleterious carrier depletion effects from the lower surface once it becomes exposed, i.e. after the sacrificial etch. The calculations indicate the donor regions are particularly susceptible to illumination, and that parallel conduction might be possible. However, in our measurements no obvious contribution from such channels even after continuous illumination was evident as will be discussed below. Hence, we only have to consider a single 2DEG layer in the evaluation below.

To realize these suspended nanostructures with three-dimensional relief we employ multiple steps of optical and electron beam lithography, followed by a combination of pattern transfer steps. The latter involve both anisotropic ion-beam (dry) and chemically-selective (wet) etching techniques. The first step of the processing procedure is fabrication of standard Au/Ni/Ge ohmic contacts to the GaAs. These require an annealing step after their deposition. The second step is definition of contacts by optical and electron beam lithography and patterning with Ni. These contacts serve as a metallic etch mask during the subsequent ion-beam etch step. The Ni etch mask is selectively removed by another brief wet etch step employing FeCl3 after the dry etch. The device geometry is machined by anisotropic, chemically assisted ion beam etching (CAIBE) with chlorine gas. The final step then is the selective chemical etch to remove the sacrificial layer. Depending on the stochiometry of the structural layer, either a 2% solution of HF or concentrated HCl was employed .

A typical freely suspended 2DEG samples is shown in Fig. 1(b): In this case the device has a shape similar to a conventional Hall bar with dimensions of $`1\times 2\mu `$m2 (width $`\times `$ length). In Fig. 2(a) a scanning electron microscope (SEM) micrograph of the same device as in Fig. 1(b) is shown, but taken at a steeper angle; the sacrificial layer is clearly seen to be completely removed from beneath the structure. In comparison to usual Hall-bars with a wide 2DEG region and small contact areas we had to choose a simplified design, since the underetching effectively limits the lateral extension of the system. For larger membranes the mechanical supports are commonly undercut as well, leading to a collapse of the whole structure. Although the contact areas are fairly large compared to the device size, the mini-Hall-bar allows us to monitor the longitudinal voltage drop properly.

Our quantum devices, also exemplified in Fig. 2(b), represent the first truly suspended quantum dots. The single dot has a final diameter of 800 nm, while the diameter of the two coupled dots is of order of 400 nm. As discussed below, the actual electronic diameter of these devices is reduced by edge depletion . In the present case, coupling between the dots is mediated by the constriction regions connecting them. Careful design of the geometry of these regions allows to control the carrier depletion and hence the degree of coupling to the leads. Fabrication of additional gate contacts in these regions, not attempted in these preliminary investigations, would permit definition of variable tunneling barriers.

At first, low temperature magnetotransport measurements were carried out to characterize the electron gas in the mini-Hall-bar and the quantum devices. Fig. 2(c) displays Shubnikov-de Haas (SdH) oscillations in the longitudinal resistance for one of the suspended Hall bar samples after illumination with a red light emitting diode for 60 sec. Prior to illumination the resistance commonly was a factor of five larger, as seen in the values obtained for the cavities A and B in Fig. 3. The reason for this deviations might be also found in a varying surface tension of the different suspended 2DEGs. The electron density is finally obtained from the SdH-oscillations shown in Fig. 2(c) through the standard relation $`\mathrm{\Delta }(1/B)=g_se/(hn_s)`$, where $`g_s`$ is the spin degeneracy factor (here $`g_s=2`$), which yields $`n_s5.5\times 10^{15}`$ m-2. Note that the Landau levels are well developed even at fairly low fields. In the inset of Fig. 2(c) the $`(1/B)`$ spectrum of the oscillations for the illuminated shows only a single period. Hence, a parallel conduction channel appears to be absent or can at least be neglected for the transport data, despite illumination and the model calculations. The most likely explanation for this is the fact that lattice strain effects are not included in the model calculations, but are expected to play a significant role for these suspended devices.

For the mini-Hall-bar the sheet resistance is of the order of $`\rho =\rho _{xx}4`$ k$`\mathrm{\Omega }/\mathrm{}`$ without illumination, and $`700\mathrm{\Omega }/\mathrm{}800\mathrm{\Omega }/\mathrm{}`$ after illumination. The values are obtained by extrapolating the resistance values from $`B=1.5`$ T to $`B0`$ T, in order to circumvent contributions by scattering around $`B=0`$ (seen in the strong negative magnetoresistance). As noted before, the length/width ratio of the Hall bars is $`l/w=2\mu `$m$`/1\mu `$m$`=2`$. Combining this zero-field resistance of approximately $`700\mathrm{\Omega }/\mathrm{}`$ with the measured carrier density, we obtain an electron mobility after illumination of $`1.41.6`$ m2/Vs. This is a factor of 27 below the value for the non-suspended 2DEG. It is likely that this strong reduction of electron mobility results from the combined effects of damage imparted by ion bombardment during the dry etch step, and from inhomogeneous surface tension in the suspended structure (see Fig 2(a)). This is currently under investigation. For future samples the etching process will be optimized in order to minimize damage to the device . At low magnetic fields negative magnetoresistance is manifested, indicating the prevalence of weak localization or interaction phenomena.

In Fig. 3(a,b) are shown the magnetoconductance traces of two different cavities with open leads as depicted in the SEM micrographs in Fig. 2(b). The devices constitute either quantum cavities or quantum dots, depending on the coupling to the leads, i.e. number of transmitting channels $`N`$ in the leads with a final conductance of $`G=N2e^2/h`$. We denote the two different structures measured as cavity A and B, since the effective total resistance is far below the resistance quantum $`h/e^2`$. The measurements on the coupled dot device are discussed below (see Fig. 5). The density obtained from the SdH-resonances at larger magnetic fields is $`n_s5.5\times 10^{15}`$ m-2 (cavity A) and $`n_s6.8\times 10^{15}`$ m-2 (cavity B). Here we do not take into account the spin degeneracy factor, since the spin splitting only emerges at $`B4`$ T (cf. arrows in Fig. 3(a,b)). The resulting mobilities are $`\mu _{cavA,cavB}=`$ 0.4 m2/Vs and 0.3 m2/Vs, strongly reduced by the barriers introduced (compared to the ’free’ 2DEG). As seen, the data at larger fields exhibit the common signature of a 2DEG, while at lower fields an intricate resonance structure is found.

In focusing on this fine structure around $`B=0`$ T, shown in Fig. 4, we are able to extract some of the properties of the electron gas in the suspended membrane: The low temperature measurement at 4.2 K yields two distinct features, namely the coarse modulation of the resistance and a fine structure with a modulation period of $`\delta B_{cavA}=(26.1\pm 1.0)`$ mT. This fine structure vanishes rapidly upon increasing the temperature (see traces for 10 K in Fig. 4(a,b)). Under illumination the carrier density varies and the $`B`$-periodicity disappears, even at 4.2 K (Fig. 4(a)).

Evaluating the periodicity of the fine structure by treating it as Aharonov-Bohm (AB) oscillations, we obtain a radius of the first cavity of $`r_{AB}^{cavA}=\sqrt{\varphi _0/(\pi \times \delta B_{cavA})}=`$ $`\sqrt{4.14\times 10^{15}\text{Tm}^2/(\pi \times 20\text{mT})}225`$ nm and with $`\delta B_{cavB}=(26.3\pm 0.5)`$ mT for the second one $`r_{AB}^{cavA}=220`$ nm. With this value we obtain an ’electronic’ diameter of the cavity of 450 nm for the first ($`A`$) and 440 nm for the second cavity ($`B`$). Compared to the lithographically defined size of 800 nm, we thus find a depletion depth of $``$ 175 nm from the edges. The fact that AB-like oscillations are found without an opening in the center could possibly be explained by some local inhomogeneity within the dot, e.g. local strain in the 2DEG, which might result in an effective depletion of the electrons in the center of the cavity. This also agrees with the disappearance of AB oscillations after illumination (higher density). Apart from the pecularity of the existence of AB-oscillations in such a sample it has to be noted that dissipation phenomena in free-standing 2DEGs seem to be less pronounced than in comparable 2DEGs connected to the bulk GaAs-crystal.

The electronic radius found from the AB like oscillations is confirmed by the coarse modulation, which can be explained by geometric resonances found earlier in strongly modulated one-dimensional electron systems . From this model we obtain the classical cyclotron radius $`r_c=\mathrm{}k_F/eB`$ with $`k_F=\sqrt{2\pi n_s}`$ being the Fermi wave vector. For the two quantum cavities measured (Fig. 3) we find a variety of different radii: For the first cavity A (Fig. 4(a)) we obtained $`r_c=652`$ nm, $`r_c^d=680`$ nm, $`r_c^c=484`$ nm, $`r_c^b=357`$ nm, $`r_c^a=265`$ nm after illumination. Prior to illumination the only cyclotron orbit found is $`r_c=r_c^c`$. Similar commensurability oscillations and an almost identical fine structure are found for cavity B (Fig. 4(b)): $`r_c^c=874`$ nm, $`r_c^b=539`$ nm, $`r_c^a=210`$ nm. The insets in Fig. 4(a,b) show some of the possible orbits: The upper sketch $`\alpha `$ shows the ideal case, while in the lower one ($`\beta `$) we allowed a finite contact width, demonstrating the occurence of the different cyclotron radii. In the third case $`\gamma `$ the electron follows a trajectory with maximum cyclotron radius (with non-ideal contacts), elucidating the large values of $`600800`$ nm measured. Considering the different radii and the corresponding orbits we obtain for the smallest radii A: $`r_c^a=265`$ nm and B: $`r_c^a=210`$ nm a reasonable agreement with the values obtained before under the assumption of AB oscillations. Since we find similar radii for the commensurability oscillations, being caused by the coarse modulation, and under the assumption of AB-oscillations, caused by the fine structure, we are confident that the depletion depth is correctly determined.

Finally, we want to present some of the first transport data on suspended quantum dots as shown in Fig. 2(b), i.e. for the case when $`G<e^2/h`$. Already strong indications were found that the relaxation mechanisms of single electrons tunneling through coupled dots is influenced by discrete phonon modes . However, these quantum dots were still embedded in the bulk GaAs crystal and hence the phonon coupling is not well controlled. The conductance measurements on a quantum dot under bias voltage strongly decoupled from the leads are shown in Fig. 5. For these measurements we employed the coupled dot system shown in Fig. 2(b), since the smaller total size and the narrower contacts result in an increased resistance. As seen we find a transition into the tunneling regime for $`B=0`$ T in the low drain/source bias ($`V_{ds}<\pm 75`$ mV) regime where the conductance $`G=dI/dV`$ is found to be below $`e^2/h`$. The conductance drops below $`e^2/h`$ for bias voltages around +20 mV to +50 mV. This drop is found to occur stepwise. Including the magnetic field dependence we are able to identify three different step widths, which are denoted $`E_C^0`$, $`E_C^1`$, and $`E_C^2`$. Interpreted as charging states of the quantum dots these steps correspond to energies approximately 40 meV, 10 meV, and 30 meV, respectively. A crude estimation of the charging energy expected from the electronic dot radius of roughly $`r_e200`$ nm (using the values found above for the depletion depth) we obtain $`E_C=e^2/2C_\mathrm{\Sigma }1`$ meV. We assumed a dielectric disk with $`C_\mathrm{\Sigma }=8ϵ_0ϵ_rr_e`$, where $`C_\mathrm{\Sigma }`$ is the total capacitance of the dot, $`ϵ_r`$ the dielectric constant of AlGaAs and $`ϵ_0=8.85\times 10^{12}`$ F/m. Although this value is far below the stepsize found in Fig. 5, it is still sufficient to allow Coulomb blockade (CB) at 4.2 K. Moreover, it is possible that the electronic diameter of the coupled dots is smaller than we estimated here.

Unfortunately it was not possible to verify charging of the quantum dot device by an external gate directly. Thus no clear CB signature can be identified for this device. Instead we applied a magnetic field perpendicular to the plane of the dot. Increasing the field, we observe a suppression and an additional modulation of the conductance around $`V_{ds}=0`$, as seen for $`B=4`$ T and 8 T. Discrete states in the quantum dot must have formed, indicated by the structure in the $`IV`$-characteristic of Fig. 5. This is evidenced by the modulation at finite magnetic field. Additionally, we find an enhanced background noise when keeping the dots under forward bias, which can be attributed to an effective mechanical vibration of the whole sample, due the resulting Lorentz force in the magnetic field.

In summary, we have demonstrated a technique of how to fabricate freely suspended two-dimensional electron gases. In transport measurements we observe the characteristic features of 2DEGs. Moreover, we find at low temperatures a strong negative magnetoresistance effect, indicating the lower mobilites of the suspended electron gas compared to the bulk values. In further measurements on suspended quantum dots and quantum cavities we found geometric resonances and a fine structure from which we estimate the lateral depletion length to be 175 nm. Finally, we presented first measurements of transport spectroscopy on a suspended quantum dot structure revealing the existence of discrete states. Although the properties of the suspended electron gas are not yet optimized in comparison to conventional AlGaAs/GaAs-heterostructures, we clearly demonstrated functioning of these devices. Moreover, this technique offers a broad variety of possible applications, including high-frequency strain sensing in nano-mechanical devices, studies on the single electron/phonon interaction, measurements of electronic specific heat and ultra-sensitive bolometry.

We gratefully acknowledge funding for this work from DARPA MTO/MEMS under grant DABT-63-98-1-0012. We thank A.N. Cleland for expert technical help and S. Ulloa for critically reading the manuscript. RHB thanks the Max-Planck-Gesellschaft (Otto-Hahn stipend), the Max-Planck-Institute for Solid State Physics, Stuttgart, and the Alexander-von-Humboldt Stiftung (Feodor-Lynen stipend) for support.

current address: Center for NanoScience and Sektion Physik, Ludwig-Maximilians-Universität München, Geschwister-Scholl-Platz 1, 80539 Munich, Germany.
e-mail: robert.blick@physik.uni-muenchen.de

new address: Institut für Angewandte und Experimentelle Physik, Universität Regensburg, 93040 Regensburg, Germany.

IBM Research Staff Member Emeritius.

Fig. 1:

(a) Calculation of the band-structure of the freely suspended two-dimensional electron gas. As indicated the filled squares mark the non-illuminated trace and the triangles the illuminated one. The Fermi level is pinned at mid-gap for top and bottom surfaces. The temperature assumed in the calculations is $`T=`$ 10 K. Top and bottom of the heterostructure are formed by GaAs cap layers, two AlGaAs:Si donor layers and two AlGaAs spacer layers around the two-dimensional electron gas.

(b) Scanning electron beam micrograph of the suspended mini-Hall-bar used in the measurements. The dimensions of this structure are length $`\times `$ width $`=2\times 1\mu `$m2, while the total thickness of the whole structure is 100 nm.

Fig. 2:

(a) Side-view of the device in Fig. 1(b), demonstrating the clear undercut; for all samples the same wafer material was used.

(b) Aerial view of the suspended single and double quantum dots. The lithographical diameter of the single cavity is 800 nm, while the double dots are on the order of 400 nm. We fabricated several of these samples; data shown in Fig. 3,4 is obtained from two different cavities A and B, while the data from the coupled dot is presented in Fig. 5.

(c) Shown is the longitudinal magnetoresistance $`R(B)`$ of one of the suspended mini-Hall-bars at $`T=4.2`$ K prior and after illumination as indicated (second sample’s data not shown). The inset gives a $`(1/B)`$-plot demonstrating the periodicity of the oscillations after illumination (electron density for samples used is obtained as: $`n_s5.5\times 10^{15}`$ m-2).

Fig. 3:

(a) Magnetoresistance trace of the suspended ’open’ quantum dot (cavity A), prior to illumination (note the large resistance). The leads contain several transport channels, resulting in a strong coupling, i.e. $`G>2e^2/h`$. In the high field region the classical Shubnikov-de Haas oscillations appear, while at low fields we find commensurability oscillations (circle). The two arrows mark the onset of the spin splitting. Inset gives a clear $`1/B`$-periodicity of the oscillations, we obtain a density of $`n_s^{cavA}=5.65\times 10^{15}`$ m-2.

(b) Corresponding Shubnikov-de Haas data of cavity B (evaluation of the data in the inset gives $`n_s^{cavB}=6.38\times 10^{15}`$ m-2). For this cavity we find even better pronounced commensurability oscillations (circle). Inset depicts the experimental setup for the transport measurements. Also marked by arrows is the occurence of the first spin splitting.

Fig. 4:

(a) At the lowest temperatures (4.2 K-plot) we find commensurability resonances and superimposed $`B`$-periodic Aharonov-Bohm type oscillations from which we derive an ’electronic’ diameter of the cavity of $``$ 450 nm for cavity A. The oscillations clearly vanish at 10 K and under illumination. The geometric resonances are also modified under illumination. The commensurability resonances equally allow the determination of the cavity diameter: Calculating the cyclotron radii we find $`r_c=652`$ nm, $`r_c^d=680`$ nm, $`r_c^c=484`$ nm, $`r_c^b=357`$ nm, $`r_c^a=265`$ nm (see text for details). Insets give some of the possible orbits: The upper sketch $`\alpha `$ shows the ideal case, while in the lower one $`\beta `$ we allowed a finite contact diameter, demonstrating the occurence of the different cyclotron radii.

(b) Similar commensurability oscillations and almost identical fine structure for cavity B. Calculating the radii of the cyclotron orbits yields: $`r_c^c=874`$ nm, $`r_c^b=539`$ nm, $`r_c^a=210`$ nm. Inset ($`\gamma `$) sketches a possible electron trajectory with maximum cyclotron radius.

Fig. 5:

Drain/source conductance of the coupled dot in the tunneling regime for different magnetic fields at 4.2 K. Marked with $`E_C^0`$, $`E_C^1`$, and $`E_C^2`$ are the different possible charging states, changing upon increasing the field (see text for details). Inset indicates sequential electron tunneling through the coupled dot.


text:
# Secure and linear cryptosystems using error-correcting codes

## Abstract

A public-key cryptosystem, digital signature and authentication procedures based on a Gallager-type parity-check error-correcting code are presented. The complexity of the encryption and the decryption processes scale linearly with the size of the plaintext Alice sends to Bob. The public-key is pre-corrupted by Bob, whereas a private-noise added by Alice to a given fraction of the ciphertext of each encrypted plaintext serves to increase the secure channel and is the cornerstone for digital signatures and authentication. Various scenarios are discussed including the possible actions of the opponent Oscar as an eavesdropper or as a disruptor.

The goal of cryptography is to enable two people, usually referred to as Alice and Bob, to communicate over an insecure channel in such a way that the opponent Oscar cannot understand and decrypt the transmitted message. A block message is called a plaintext, and a long message is a sequence of plaintexts. In a general scenario, the plaintext is encrypted by Alice through the key $`E_k`$ and the result, ciphertext, is sent over the channel. A third party, eavesdropping on the channel, cannot determine what the plaintext was. However, Bob, who knows the encryption key, can decrypt the ciphertext using the key $`D_k`$ and recover the plaintext.

In a private-key system, the keys $`E_k`$ and $`D_k`$ are known only to Alice and Bob, and it obviously increases the security of the channel. However, a private-key system requires communication between Alice and Bob prior to the transmission of any plaintext. This prerequisite makes the private-key communication impractical in modern communication, especially in such areas as electronic commerce and Internet-based communication. The goal of public-key systems is to devise a cryptosystem where it is computationally infeasible to determine $`D_k`$ given $`E_k`$, and hence the encryption rule $`E_k`$ can be made public.

The secure channel and the efficiency of a public-key cryptosystem depends on many parameters, among them: (a) the complexity to determine $`D_k`$ given $`E_k`$; (b) the complexity of the encryption/decryption processes; (c) the length of the ciphertext and the public-key in comparison to the length of the plaintext.

The commonly used RSA cryptosystem is based on the difficulty of factorizing large integers. Its main drawback is that the complexity of the encryption/decryption processes is of $`O(N^2)/O(N^3)`$, where $`N`$ is the length of the plaintext (see figure 1). For small $`N`$ these complexities are also small, but then the complexity to determine $`D_k`$ given $`E_k`$ may also be accessible to Oscar. It was recently found that even $`N=512`$ may be too small to ensure a secure channel. Hence the complexity of the encryption/decryption becomes the bottleneck of public-key cryptosystems as well as for other tasks of the secure channel (digital signature, authentication, etc.) based on such methods.

From the known cryptosystems it appears that there is a trade off between the secure channel and the complexity of the encryption/decryption processes. In this work, we propose a new secure cryptosystem, based on the preliminary bridge built between error-correcting codes and cryptosystems by McEliece with the following features and ingredients: (a) The complexity of the encryption/decryption processes scale linearly with the size of the plaintext $`N`$. These complexities can be easily reduced even further under parallel dynamics. (b) The method is based on boolean operations between two sparse matrices, in contrast of factorizing large integers in cryptosystems based on number theory. (c) Our method consists of many stochastic ingredients; Bob adds noise to the public-key whereas Alice adds noise to the ciphertext. (d) The method is applicable as a public-key cryptosystem, as well as for digital signatures and authentication. A digital signature is used to specify the person responsible for the message, and an authentication ensures the integrity of the plaintexts constructing the message. It is a challenge to have a secure public-key cryptosystem operating with low complexity which can serve also for all the different tasks of the secure channel.

Before describing the details of our method, let us first categorize the possible capabilities of Oscar: (a) An eavesdropper: Oscar may try to reveal the plaintext Alice is sending to Bob from the transmitted ciphertext and/or the digital signature. (b) A distruptor: The message Alice is sending to Bob can be repeated, replaced or corrupted during transmission by Oscar. We may distinguish between a forged meaningful/meaningless signed plaintext. Note that the ability to forge many meaningless but legally signed messages could cause a disastrous effect in the event of real-time procedures. It may take some critical time for Bob to realize that legally signed messages are forged messages rather than noisy ones. Note that in cryptosystems such as RSA, it is easy to forge a meaningless signed message or to repeat the transmission of the same message or previously legally signed messages to Bob.

Our cryptosystem is based on an error correcting method known as the Gallager method or its MN version. It comprises two sparse boolean matrices, $`A`$ and $`B`$, of dimensionalities $`M\times N`$ and $`M\times M`$ respectively, and the rate $`RN/M1`$. Note that all operations are in (mod 2). The encryption public-key is an $`M\times N`$ matrix

$$E_k=B^1A\text{(mod 2)}$$
(1)
and the encrypted ciphertext is

$$C=E_ks\text{(mod 2)}$$
(2)
Before the transmission, Alice adds noise (bits flipped with probability $`f`$) to $`C`$, such that the received ciphertext is

$$r=C+n\text{(mod 2)}$$
(3)
where $`n`$ represents the noise. The decryption key

$$D_k=[A,B]$$
(4)
is known only to Bob, who can find $`s`$ by multiplying $`r`$ with the matrix $`B`$ to obtain $`z=B(E_ks+n)=As+Bn`$. It requires solving the equation

$$[A,B]\left[\begin{array}{c}s^{}\\ n^{}\end{array}\right]=z,$$
(5)
where $`s^{}`$ and $`n^{}`$ are the unknowns, but their statistics (for instance, unbiased message for $`s^{}`$, and flip rate $`f`$ for $`n^{}`$) are known. This may be carried out using standard methods such as that of belief network decoding. In this method, representing a special case of parity-check codes, each bit of the ciphertext $`C`$ is derived from the parity of a sum of certain plaintext’s bits. It has recently been shown, based on insight gained from the study of diluted spin systems that specific choices of cascading sub-matrices $`A`$ and $`B`$ can nearly saturate Shannon’s bound.

With the lack of noise and invertible $`E_k`$ (otherwise decryption cannot be terminated successfully with probability one), Oscar is able to easily recover $`s=E_k^1r`$. In order to make Oscar’s task more difficult, we follow the line of McEliese where noise is added to $`C`$, namely, bits are flipped ($`01`$ or $`10`$) with probability $`f`$. In this event $`E_k^1r`$ results in an approximated plaintext, $`s_{app}`$, where a fraction of the bits are wrong. For a given rate $`R`$ and large $`N`$, the maximal noise $`f`$ (for which the decryption could terminate successfully without error bits in the decrypted plaintext) is given by the maximal channel capacity $`R_c=1H_2(f)`$ where $`H_2(f)=f\mathrm{log}_2(f)(1f)\mathrm{log}_2(1f)`$. Oscar’s task to recover $`s`$ is difficult, since he has to decompose $`E_k`$ into the matrices $`B^1`$ and $`A`$, which is known to be an NP-complete problem.

The main drawbacks of this cryptosystem are: (a) Strong finite size effects which are visible even for large $`N=O(10^4)`$. The decryption typically terminates with some percentages of error bits, which is catastrophic from a practical point of view. (b) For large $`N`$, the encryption evolves a product of a matrix $`(M\times N)`$ $`E_k`$ by the plaintext $`s`$, hence its complexity is $`O(N^2)`$. Similarly, the complexity of each step of the decryption is $`O(N^2)`$, eq. (5). Clearly it is less than the cubic complexity of the decryption in RSA. However, if in practice one has to work with $`N=10^4`$, the reduction in the complexity becomes a drawback. Furthermore, the size of the public-key which has to be downloaded by Alice diverges as $`O(N^2)`$.

To overcome these drawbacks, which prevent the usefulness of error-correcting codes as a cryptosystem, we used the following three observations. In error-correcting codes, the nature of the channel typically has a statistical nature; a probability for a bit to flip or a width of the Gaussian noise. Each transmitted bit has the same probability to flip from a given distribution. The first observation is that in a case where an error-correcting method serves as a cryptosystem, Bob can pre-corrupt the public-key $`E_k=B^1A`$ in the following sense. In a fraction $`p_q`$ of rows, part (all) of the elements are flipped at random. The location of these $`p_qM`$ rows is known only to Bob. Hence, a fraction $`p_q`$ of the ciphertext is corrupted with an average probability $`1/2`$. Since this pre-corruption of $`E_k`$ is common to all ciphertexts, we denote it as a quenched noise, $`n_q`$. The main purpose of the quenched noise is to make Oscar’s task of decomposing $`E_k`$ to $`B^1`$ and $`A`$ more difficult. Note that our cryptosystem works properly for some sub-classes of matrices $`A`$ and $`B`$, but for a random choice our cryptosystem fails with probability one. Hence the task of Oscar is to find a decomposition which works properly as a cryptosystem. The second observation: in addition to the pre-corruption process, Bob publicizes a given fraction, $`p`$, of the ciphertext where Alice’s private-noise, $`n_a`$, can be added. This localized private-noise consists of a flip rate $`f`$ of given $`pM`$ bits of the ciphertext. The resulting ciphertext then comprises of frozen (non-flipped) bits, randomly flipped bits and flipped bits with probability $`f`$. The presence of frozen/flipped bits in the plaintext serves to increase the secure channel and to suppress finite size effects. Similar to Shannon’s bound, one can show that for a given rate $`R`$ the maximal fraction of flipped bits with probability $`f`$ is

$$p_c=\frac{1p_qR}{H_2(f)}$$
(6)
We assume that a fraction $`p_q`$ of the bits are flipped with probability $`1/2`$, however, $`p_c`$ might even be further improved for the following reason. In an error-correction scenario only statistical properties of the plaintext and the flip rate are known, hence any decoded state obeying these statistical features is valid. In contrast, Bob knows the manner in which $`E_k`$ was corrupted and hence the error in the $`p_qM`$ corrupted bits should be consistent with the decrypted plaintext.

The most striking observation comes from various simulations on different random constructions of the matrix $`B`$, indicating the following rule. As long as the average connectivity, number of non-zero elements per column, of $`B`$ is smaller than $`2`$, $`B^1`$ is sparse. A random construction means that the elements of each row are chosen at a random position, with no spatial structure. Since $`A`$ is a sparse matrix with random positions of the non-zero elements, it is clear that $`E_k`$ is also sparse. Hence, for such matrices $`B`$, the size of the public-key $`E_k`$ scales linearly with the size of the plaintext. Furthermore, the complexity of the decryption process also scales linearly with the size of the plaintext, as the number of iterations is of $`O(1)`$ (see details below). A comparison with encryption/decryption complexities of the RSA system is presented in figure 1 for various sizes of plaintexts. A sparse public-key is a necessary requisite for an efficient encryption process of large plaintexts, which are of great practical importance. For an average connectivity greater than $`2`$, $`B^1`$ is heavily dense, and the number of non-zero elements in $`E_k`$ is around $`MN/2`$.

The sparseness of $`B`$ for a connectivity less than $`2`$ can be supported by the following theoretical argument. Assume that the matrix $`B`$ is constructed such that $`B_{ij}=\delta _{i,j}+\delta _{i+c,j}`$ for $`i\rho M`$ and $`B_{ij}=\delta _{i,j}`$ for $`i>\rho M`$, where $`\rho <1`$ and the average connectivity is $`1+\rho <2`$. This matrix can easily be inverted and one can show that for $`c=O(1)`$ $`B^1`$ is dense (with a small prefactor), but for $`c=O(N)`$, $`B^1`$ is sparse. In a random construction the typical distance between two non-zero elements belonging to the same row is of $`O(N)`$. With respect to this property it is similar to $`c=O(N)`$. Furthermore, for connectivity below $`2`$ the generic graph represented by $`B`$ is below the percolation threshold.

We perform simulations on $`R=1/2`$ and $`256N2048`$ with a few different classes of matrices $`A`$ and $`B`$ and here we report only limited results. Each parameterization of the matrices $`A`$ and $`B`$ was averaged over at least $`10^5`$ plaintexts and $`50`$ realizations. The construction follows the spirit of the constructions for error-correcting codes of the Gaussian channel and $`R=1/2`$. The structure of the matrix $`B`$ is such that for $`i\rho M`$ there are two non-zero elements at random positions where for $`i>\rho M`$ there is only one non-zero element. The matrix $`A`$ is constructed such that in the first $`\rho ^{}M`$ rows there are two non-zero elements and in the remaining rows there are six non-zero elements chosen at random positions. In order to break the inversion symmetry, where each row of $`A`$ consists of an even number of non-zero elements, a small number of rows were changed from $`21`$ and $`76`$ non-zero elements. Note that the spatial separation between different rows of the matrices was done only for demonstration, and to increase the security of the channel one can mix their locations. The performance depends on the success rate of the decryption as a function of ($`N,p,p_q,f`$), where the private-noise was added to the first $`pM`$ bits of the ciphertext. Let us present a few examples among many where the decryption terminates successfully over at least $`10^5`$ plaintexts in a finite fraction of the realizations: (a) $`\rho =1/2`$, $`\rho ^{}=7/8`$ and $`(512,0.53,00.04,0.04)`$, (b) $`\rho =\rho ^{}=3/4`$ and $`(1024,0.53,00.04,0.075)`$ and (c) $`\rho =3/4`$, $`\rho ^{}=7/8`$ and $`(768,0.53,00.04,0.088)`$. These results indicate that the probability for a wrongly decrypted block (plaintext) is $`P_B<10^5`$. In all the above-mentioned classes, the number of iterations of the belief algorithm is typically $`10`$ steps, where the complexity of each step of the algorithm is of the order of the number of non-zero elements in matrices $`A`$ and $`B`$, $`O(N)`$. No long tail in the distribution of the convergence time was observed. Note that each belief iteration can be implemented in parallel such that the time complexity can be reduced by $`O(1/N)`$. Results indicate that finite size effects are dramatically suppressed by the frozen bits (in contrast to homogeneous noise), and can be improved even further by increasing $`N`$.

The location of frozen (non-flipped) bits of the ciphertext $`(1p)M`$ is known also to Oscar. Hence, the secure channel forces the number of frozen bits $`(1p)M<N`$. Otherwise Oscar may try to solve $`E_k^{1p}s_{1p}=C_{1p}`$, where $`1p`$ indicates the relevant part of the matrix/vector corresponding to the frozen bits. In such a case where $`E_k^{1p}`$ (of dimentionalities $`(1p)M\times N`$) is invertible, Oscar can easily find the plaintext $`s`$.

Let us now discuss a possible attack on our cryptosystem. Oscar’s goal is to find a partial public-key, $`E_k^{part}`$, obeying the following constraints: (a) The dimentionalies of $`E_k^{part}`$ is $`M^{}\times N`$ where $`NM^{}M`$. (b) The corresponding $`M^{}`$ bits of the ciphertext are the correct ones. (c) $`E_k^{part}`$ is invertible. In such an event Oscar can easily find the plaintext $`s`$, and the question is, what is the probability of such an event? The number of frozen bits of the ciphertext is $`(1p)M`$ which was chosen to be less than $`N`$. Assuming these $`(1p)M`$ rows are linearly independent, Oscar has to guess additional $`N(1p)M=N(R+p1)/R`$ correct rows in order to construct a plausible invertible $`E_k^{part}`$. The probability of such an event is $`(1f)^{NM(1p)}`$ and it becomes negligible as we increase the size of our plaintext. Furthermore, in simulations we realized that the rank of the $`(1p)M`$ correct rows is is $`0.9(1p)M`$. Hence Oscar has to guess additional correct rows and the probability of such an event decreases even further. Last but not least, how does Oscar know that he chose the correct rows? A plausible answer to this question is that any set of additional correct rows results in the same plaintext $`s`$. Hence, Oscar may repeat the above attack many times and deduce that the most probable outcome is the desired plaintext. The first difficulty with this scenario is that $`s`$ has to appear many times as an outcome, since we would like to distinguish between the signal, $`s`$, and the noise of other possible outcomes. Secondly, the most probable (exponentially dominated) $`E_k^{part}`$ consists of $`(NM(1p))f>0`$ wrong rows. It is true that all these $`E_k^{part}`$ do not necessarily result in the same plaintext, but their distribution is still in question.

The secret information of Bob is the decomposition of the public-key $`E_k`$ into $`B^1`$ and $`A`$ and the quenched noise ,$`n_q`$, used to corrupt the public-key. The secret information of Alice is the plaintext $`s`$ and the private noise, $`n_a`$. The key point of our signature scheme is that after the decryption process terminates successfully Bob recovers not only the plaintext $`s`$ but also the private noise, $`n_a`$, added to the ciphertext. More precisely, on one hand side, from the decryption of the plaintext $`s`$ Bob knows the corrupted ciphertext by using the corrupted public-key, $`E_ks`$. On the other hand, Bob has in his hand the received ciphertext, $`E_ks+n_a`$. From the difference between these two pieces of information Bob can easily find $`n_a`$. The ability of Bob to reveal $`n_a`$ besides the plaintext is at the center of the discussion below. It shows how Alice can use the additional information to sign and to keep the integrity of the message.

A simple signature is based on the following two ingredients: (a) Alice constructs an additional plaintext comprising of a linear combination of $`s`$ and $`n_a`$, $`X(s,n_a)`$. The new plaintext $`X`$ is encrypted by Alice to a new ciphertext, $`t_1`$, using $`E_k`$ and a new private noise $`n_{a1}`$. Alice transmits to Bob both ciphertexts, $`t`$ and $`t_1`$. (b) For verification, Alice constructs by a known procedure a verifiable vector, $`V=V(s,n_a,n_{a1})`$. After the decryption of both ciphertexts Bob knows all the ingredients of $`V`$ and the verification can be carried out. Note that for a one-time signature scheme where Oscar is functioning as an eavesdropper only, our channel is secure. The usefulness of these signature schemes is twofold: (a) The signature/verification procedure is very easy for Alice/Bob to implement with complexities of $`O(N)`$. (b) A plaintext repeated twice has in each transmission a different signature due to the different private-noise. The main drawback of the above signature scheme is that Oscar can easily forge a legal plaintext. There are exponentially many plaintexts $`s`$ and private-noise $`n_a`$ and $`n_{a1}`$ which give the same verifiable vector $`V`$

It appears that in order to have a secret personal signature Alice has (a) to send additional information besides the plaintext; and (b) to use her own known signature scheme based on her cryptosystem, similar to using RSA for both the encryption and signature. Surprisingly, we demonstrate below that Alice can construct a secure signature without the transmission of any additional information besides the encrypted plaintext.

A simple scenario for an advanced secure signature is one in which Alice first generates a vector $`V`$ of rank $`N^{}<N`$ using $`s`$ and $`n_a`$ following her public protocol. Next, the number of $`1`$’s in $`V`$ is truncated to a fixed number $`K`$ (or $`K`$) following Alice’s public prescription. (For rare events where there are no $`1`$’s in $`V`$, Alice provides a special procedure). The signature, $`E_k^{part}V`$ is left publicized by Alice, where $`E_k^{part}`$ stands for the relevant rows corresponding to the rank $`N^{}`$. Determining $`V`$ from the knowledge of $`E_k`$ and the signature is known to be a NP-complete problem (pg. 280 in ref. (14)). Bob, who knows $`s`$ and $`n_a`$, can easily verify the signature. One can easily make the situation more complex when Alice creates a new cryptosystem in an on-line manner. For instance, based on the public-key, $`E_k`$, Alice implements some permutations among the rows as a function of the detailed structure of $`s`$ (and/or $`n_a`$). The permutation scheme is publicized by Alice, which may be fixed for all plaintexts/users or can be chosen time dependently. Let us denote the permute public-key by $`E_k^P`$, hence the signature of $`s`$ is $`t_1=E_k^PV`$, publicized by Alice but does not have to be transmitted over the channel. Bob first decrypts $`t`$ and obtains $`s`$ and $`n_a`$. Next, Bob finds his permuted public-key, $`E_k^P`$, using the public permutation prescription chosen by Alice, and then easily verifies $`t_1`$ from $`E_k^PV`$. Since the signature depends on $`s`$ and $`n_a`$ as well as on $`E_k`$, the same plaintext transmitted to different addresses or at different times (different $`n_a`$) is characterized by different signatures. As an eavesdropping, Oscar does not know $`s`$, $`n_a`$ and also $`E_k^P`$. As a disruptor, Oscar may try to replace the ciphertext $`t`$ by $`t^{}`$ consisting of $`s^{}`$ and $`n_a^{}`$ such that the signature of $`V=V(n_a^{})`$ using the permuted public-key $`E_k^P(s^{})`$ has the same signature as that publicized by Alice. The lack of an independent permuted public-key as a function of the plaintext seems to make the work of a disruptor even harder. In principle Alice can corrupt her signature $`E_k^PV`$ either by flipping some bits among a small fraction of the ciphertext or by adding a noise consisting of $`n_a`$ and $`s`$. In such a scenario Bob has first to decrypt $`s`$ and $`n_a`$ and then to verify the decrypted signature with the public protocol of Alice. For a permuted public key the decryption follows the permuted matrices $`A_{per}=A`$, and $`B_{per}`$ is nothing more than the same permutation as made for the rows of $`E_k`$, but now for the columns of $`B`$.

The aim of the authentication procedure is to keep the integrity of the message constructed from a sequence of plaintexts, such that Oscar cannot forge (add/delete) ciphertexts. Using error correcting codes as a cryptosystem this goal can be achieved by using correlated noise for successive ciphertexts. Let us describe some possible scenarios. The private-noise for the first ciphertext is chosen as explained above, where the noise for the next ciphertext is related to the previous one by some permutations. The permutations may depend on (a) the private-noise of the previous ciphertext; (b) the previous plaintext; (c) both the previous ciphertext, and the private-noise. One may think that the permuted noise has to be bounded to the allowed regime by Bob in order to ensure a successful decryption (which can be easily achieved). However, there is no necessity for such a restriction, since after the decryption of the first plaintext, the private-noise of the next ciphertext is uniquely determined and hence also the plaintext.

The advantage of such an authentication scheme is that Bob has only to decrypt the first plaintext, whereas the rest of the message is uniquely defined, since the noise is known. On the other hand, Oscar knows the authentication scheme and may concentrate only on the decryption of the first ciphertext, or alternatively on an intermediate ciphertext (the easy one) which reveals all successive plaintexts. In order to ensure the same security of (almost) all plaintexts, one can use accumulated permutations. The private-noise for the current ciphertext depends on all previous plaintexts/private-noise by a publicized procedure.

The tasks of our cryptosystem can be extended to other functions of the secure channel, such as an undeniable signature. The private-noise is added out of the allowed range such that the decryption cannot terminate successfully without Alice partially revealing her noise.

Alice has to keep as public information all previous signatures. The list of the signatures may load Alice’s resources, and furthermore it may take a long time for Bob to find the appropriate signature among many. This drawback can be alleviated by removing the signature into an archive after verification by Bob.

In conclusion, let us briefly discuss a few of the advantages of our cryptosystem over methods based on numbers theory, such as an RSA cryptosystem. First, the matrix operations/belief network decoding in the decryption/encryption process can be carried out and implemented in parallel. Secondly, a one-time success by Oscar to reveal a plaintext does not automatically help or ensure the recovery of other plaintexts that Alice sent to the same Bob. Finally, in the RSA method, for instance, Oscar’s task requires a check of many possible trails, where each trail can be examined by the same algorithm. Hence, the task of Oscar can be easily split among many resources. In contrast, our cryptosystem is based on many stochastic ingredients with time dependent features of Alice and Bob. Hence the strategy of Oscar may need to vary between different messages and users of the channel. Even for a given channel, the challenge for future research is to find how to parallelize and to simplify the task of Oscar.

I. K. acknowledges fruitful discussions and comments by W. Kinzel. The content of this work is patent pending.