format: text

text:
# Efficient Synthesis of Linear Reversible Circuits

## 1 Introduction

A reversible or information lossless circuit is one that implements a bijective function, or loosely, a circuit where the inputs can be recovered from the outputs and all output values are achievable. A major motivation for studying reversible circuits is the emerging field of quantum computation . A quantum circuit implements a unitary function, and is therefore reversible. Circuit synthesis for reversible computations is an active area of research . The goal in circuit synthesis is, given a gate library, to synthesize an efficient circuit performing a desired computation. In the quantum context, the individual gates correspond to physical operations on quantum states called qubits, and therefore reducing the number of gates in the synthesis generally leads to a more efficient implementation.

Linear reversible classical circuits form an important sub-class of quantum circuits, which can be generated by a single type of gate called a C-NOT gate (see Figure 1c). This gate is an important primitive for quantum computation because it forms a universal gate set when augmented with single qubit rotations . Moreover, current quantum circuit synthesis algorithms can generate circuits with strings of C-NOT gates, and therefore more efficient synthesis for these classical linear reversible sub-circuits would imply more efficient synthesis for the overall quantum computation.

In this paper we consider the problem of efficiently synthesizing an arbitrary linear reversible circuit on $`n`$ wires using C-NOT gates. This problem can be mapped to the problem of row reducing a $`n\times n`$ binary matrix. Until now the best synthesis methods have been based on standard row reduction methods such as Gaussian elimination and LU-decomposition, which yield circuits with $`O(n^2)`$ gates . However, the best lower bound leaves open the possibility that synthesis with as few as $`O\left(n^2/\mathrm{log}n\right)`$ gates in the worst case may exist .

We present a new synthesis algorithm that meets the lower bound, and is therefore asymptotically optimal up to a multiplicative constant. Furthermore, our algorithm is also asymptotically faster than previous methods. Simulation results show that the proposed algorithm outperforms previous methods even for relatively small $`n`$. Generically our algorithm can be interpreted as a matrix decomposition algorithm, that yields an asymptotically efficient elementary matrix decomposition of a binary matrix. Generalizations to matrices over larger finite fields are straightforward.

## 2 Background

We can represent the action of an $`n`$-input $`m`$-output logic gate as a function mapping the values of the inputs to those of the outputs: $`f:𝔽_2^n𝔽_2^m`$, where $`f`$ maps each element of $`𝔽_2^n`$ to an element in $`𝔽_2^m`$. Here $`𝔽_2`$ is the two-element field, and $`𝔽_2^n`$ is the set of all $`n`$-dimensional vectors over this field. A gate is reversible if this function is bijective, that is, $`f`$ is one-to-one and onto. Intuitively, this means that the inputs can be uniquely determined from the outputs and all output values are achievable. For example, the AND gate (Figure 1a) is not reversible since it maps three input values to the same output value. The NOT gate (Figure 1b), on the other hand, is reversible since both possible input values yield unique output values, and both possible output values are achievable. The controlled-NOT or C-NOT gate, shown in Figure 1c, is another important reversible gate. This gate passes the first input, called the control, through unchanged and inverts the second, called the target, if the control is a one. As its truth table shows, this gate is reversible since it maps each input vector to a unique output vector and all output vectors are achievable.

A reversible circuit is a directed acyclic combinatorial logic circuit where all gates are reversible and are interconnected without fanout . An example of a reversible circuit consisting of C-NOT gates is shown in Figure 2. Note that, as is the case for reversible gates, the function computed by a reversible circuit is bijective.

We say a circuit or gate, computing the function $`f`$, is linear if $`f(x_1x_2)=f(x_1)f(x_2)`$ for all $`x_1,x_2𝔽_2^n`$, where $``$ is the bitwise XOR operation. The C-NOT gate is an example of a linear gate:

$$\begin{array}{ccc}f\left(\left[\mathrm{0\; 0}\right]\right)f\left(x\right)=f\left(x\right)& f\left(\left[\mathrm{0\; 1}\right]\right)f\left(\left[\mathrm{1\; 0}\right]\right)=f\left(\left[\mathrm{1\; 1}\right]\right)& f\left(\left[\mathrm{1\; 0}\right]\right)f\left(\left[\mathrm{1\; 1}\right]\right)=f\left(\left[\mathrm{0\; 1}\right]\right)\\ f\left(x\right)f\left(x\right)=f\left([\mathrm{0\; 0}]\right)& f\left(\left[\mathrm{0\; 1}\right]\right)f\left(\left[\mathrm{1\; 1}\right]\right)=f\left(\left[\mathrm{1\; 0}\right]\right).& \end{array}$$

The action of any linear reversible circuit on $`n`$ wires can be represented by a linear transformation over $`𝔽_2`$. In particular, we can represent the action of the circuit as multiplication by a non-singular $`n\times n`$ matrix $`A`$ with elements in $`𝔽_2`$:

$$Ax=y,$$

where $`x`$ and $`y`$ are $`n`$-dimensional vectors representing the values of the input and output bits respectively. Using this representation, the action of a C-NOT gate corresponds to multiplication by an elementary matrix, which is the identity matrix with one off-diagonal entry set to one. Multiplication by an elementary matrix performs a row operation, the addition of one row of a matrix or vector to another. Applying a series of C-NOT gates corresponds to performing a series of these row operations on the input vector or equivalently to multiplying it by a series of elementary matrices. For example, the linear transform computed by the circuit in Figure 2 is given by

$$A=\stackrel{G_6}{\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 1& 1\end{array}\right]}\stackrel{G_5}{\left[\begin{array}{cccc}1& 1& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]}\stackrel{G_4}{\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 1& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]}\stackrel{G_3}{\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 1& 1& 0\\ 0& 0& 0& 1\end{array}\right]}\stackrel{G_2}{\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 1& 1\end{array}\right]}\stackrel{G_1}{\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]}=\left[\begin{array}{cccc}1& 0& 1& 0\\ 0& 0& 1& 0\\ 1& 1& 1& 0\\ 1& 1& 0& 1\end{array}\right]$$

We can use the matrix notation to count the number of different $`n`$-input linear reversible transformations. In order for the transformation to be reversible, its matrix must be non-singular, in other words, all nontrivial sum of the rows should be non-zero. There are $`2^n1`$ possible choices for the first row, all vectors except for the all zeros vector. There are $`2^n2`$ possible choices for the second row, since it cannot be the equal to the first row or the all zeros vector. In general, there are $`2^n2^{i1}`$ possible choices for the $`i`$th row, since it cannot be any of the $`2^{i1}`$ linear combinations of the previous $`i1`$ rows (otherwise the matrix would be singular). Therefore there are

$$\underset{i=0}{\overset{n1}{}}\left(2^n2^i\right)$$

unique $`n`$-input linear reversible transformations.

Since any non-singular matrix $`A`$ can be reduced to the identity matrix using row operations, we can write $`A`$ as a product of elementary matrices. Therefore, any linear reversible function can be be synthesized from C-NOT gates. Moreover, the problem of C-NOT circuit synthesis is equivalent to the problem of row reduction of a matrix $`A`$ representing the linear reversible function: any synthesis of the circuit can be written as a product of elementary matrices equal to $`A`$ and any such product yields a synthesis. The length of the synthesized circuit is given by the number of elementary matrices in the product. Standard Gaussian elimination and LU-decomposition based methods requires $`O(n^2)`$ gates in the worst-case . However, the best lower bound is only $`\mathrm{\Omega }\left(n^2/\mathrm{log}n\right)`$ gates .

###### Lemma 1 (Lower Bound)

There are $`n`$-bit linear reversible transformation that cannot be synthesized using fewer than $`\mathrm{\Omega }(n^2/\mathrm{log}n)`$ C-NOT gates.

Proof Let $`d`$ be the maximum number of C-NOT gates needed to synthesize any linear reversible function on $`n`$ wires. The number of different C-NOT gates which can act on $`n`$ wires is $`n(n1)`$. Therefore the number of unique C-NOT circuit with no more than $`d`$ gates must be no more than $`\left(n^2n+1\right)^d`$, where we have included a do-nothing NOP gate in addition to the $`n^2n`$ C-NOT gates to account for circuits with fewer than $`d`$ gates. Since the number of circuits with no more than $`d`$ C-NOT gates must be greater than the number of unique linear reversible function on $`n`$ wires, we have the inequality

$$\left(n^2n+1\right)^d\underset{i=0}{\overset{n1}{}}\left(2^n2^i\right)2^{n(n1)}.$$
(1)
Taking the $`\mathrm{log}`$ of both the left and right sides of the equations gives

$$d\frac{n(n1)\mathrm{log}2}{\mathrm{log}(n^2n+1)}=\frac{n^2n}{\mathrm{log}_2(n^2n+1)}=\mathrm{\Omega }\left(\frac{n^2}{\mathrm{log}n}\right).$$
(2)
$`\mathrm{}`$

This lemma suggests a synthesis more efficient than standard Gaussian elimination may be possible. The multiplicative constant in this lower bound is $`1/2`$ (assuming logs are taken base 2).

## 3 Efficient Synthesis

In this section we present our synthesis algorithm, which achieves the lower bound given in the previous section. In Gaussian elimination, row operations are used to place ones on the diagonal of the matrix and to eliminate any remaining ones. One row operation is required for each entry in the matrix that is targeted. Since there are $`n^2`$ matrix entries, $`O(n^2)`$ row operation are required in the worst case. If instead we group entries together and use single row operations to change these groups, we can reduce the number of row operation required, and therefore the number of gates needed to synthesize the circuit.

The basic idea is as follows. We first partition the columns of the $`n\times n`$ matrix into sections of no more than $`m`$ columns each. We call the entries in a particular row and section a sub-row. For each section we use row operations to eliminate sub-row patterns that repeat in that section. This leaves relatively few ($`<2^m`$) non-zero sub-rows in the section. These remaining entries are handled using Gaussian elimination. If $`m`$ is small enough ($`<\mathrm{log}_2n`$), most of the row operations result from the first step, which requires a factor of $`m`$ fewer row operations than full Gaussian elimination. As with the Gaussian elimination based method, our algorithm is applied in two steps; first the matrix is reduced to an upper triangular matrix, the resulting matrix is transposed, and then the process is repeated to reduce it to the identity. Detailed pseudo-code for our algorithm is shown on the next page.

The following example illustrates our algorithm for a $`6`$-wire linear reversible circuit.

1) Choose $`m=2`$ and partition matrix.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 1& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}1& 1\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \end{array}\right]$$

2) (Step A - section 1) Eliminate duplicate sub-rows.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 1& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 1\\ \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}& \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}\colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}& \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \end{array}\right]\stackrel{\begin{array}{c}14\\ 15\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 1& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 1\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \end{array}\right]$$

3) (Step B - section 1, column 1) One already on diagonal.

4) (Step C - section 1, column 1) Remove remaining ones in column below diagonal.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ \text{1}& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 1\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \end{array}\right]\stackrel{\begin{array}{c}12\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 1\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \end{array}\right]$$

3) (Step B - section 1, column 1) One already on diagonal.

4) (Step C - section 1, column 1) Remove remaining ones in column below diagonal.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \end{array}\right]\stackrel{\begin{array}{c}23\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \end{array}\right]$$

5) (Step A - section 2) Eliminate duplicate sub-rows below row 2.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}\colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{0}$}& \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}\\ \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}& \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \end{array}\right]\stackrel{\begin{array}{c}35\\ 46\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 1& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 0& 1\end{array}& \end{array}\right]$$

6) (Step B - section 2, column 3) Place one on diagonal.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}\text{0}& 1\\ 1& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 0& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}43\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \begin{array}{cc}1& 1\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 0& 1\end{array}& \end{array}\right]$$

7) (Step C - section 2, column 3) Remove remaining ones in column below diagonal.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ \text{1}& 1\end{array}& \begin{array}{cc}1& 1\\ 1& 1\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 0& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}34\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}1& 1\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 0& 1\end{array}& \end{array}\right]$$

8) Matrix is now upper triangular. Transpose and continue.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 1\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 1& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}1& 1\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 0& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}\text{transpose}\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 1& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 1\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]$$

9) (Step A - section 1) Eliminate duplicate sub-rows.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 1& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}\colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{0}$}& \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}45\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 1& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]$$

| 10) | (Step B - section 1, column 1) Because matrix is triangular and non-singular there |
| --- | --- |
|  | will always be ones on the diagonal. |

11) (Step C - section 1, column 1) Remove remaining ones in column.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ \text{1}& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}12\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 1\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]$$

12) (Step C - section 1, column 2) Remove remaining ones in column.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}24\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 1& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]$$

13) (Step A - section 2) Eliminate duplicate sub-rows.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}& \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{0}$}\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}36\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 1\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]$$

14) (Step C - section 2, column 1) Remove remaining ones in column.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}\colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}& 1\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}35\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 1\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]$$

15) (Step C - section 2, column 2) Remove remaining ones in column.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& \colorbox[rgb]{0.749019607843137,0.749019607843137,0.749019607843137}{$\text{1}$}\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}45\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 1& 1\end{array}& \end{array}\right]$$

16) (Step C - section 3, column 1) Remove remaining ones in column.

$$\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ \text{1}& 1\end{array}& \end{array}\right]\stackrel{\begin{array}{c}56\end{array}}{}\left[\begin{array}{cccc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}& \end{array}\right]$$

The synthesized circuit is specified by the row operations and is shown in Figure 3.

In general, the length of the synthesized circuit is given by the number of row operations used in the algorithm. By accounting for the maximum number of row operations in each step, we can calculate an upper bound on the maximum number of gates that could be required in synthesizing an $`n`$-wire linear reversible circuit. C-NOT gates are added in the steps marked Step A-C in the algorithm. Step A is used to eliminate the duplicates in the subsections. It is called fewer than $`n+m`$ times per section (combined for the upper/lower triangular stages of the algorithm), giving a total of no more than $`(n+m)n/m`$ gates. Step B is used to place ones on the diagonal. It can be called no more than $`n`$ times. Step C is used to remove the ones remaining after all duplicate sub-rows have been cleared. Since there are only $`2^m`$ $`m`$-bit words, there can be at most as many non-zero sub-rows below the $`m\times m`$ sub-matrix on the diagonal. Therefore, Step C is called fewer than $`m(2^m+m)`$ times per section, or fewer than $`2n/mm(2^m+m)`$ times in all. Adding these up we have

total row ops $``$ $`(n+m){\displaystyle \frac{n}{m}}+n+2{\displaystyle \frac{n}{m}}m\left(2^m+m\right)`$ (3)
$``$ $`{\displaystyle \frac{n^2}{m}}+n+n+m+n+2n2^m+2nm+2m2^m+2m^2.`$ (4)
If $`m=\alpha \mathrm{log}_2n`$,

$`\text{total row ops}{\displaystyle \frac{n^2}{\alpha \mathrm{log}_2n}}`$ $`+`$ $`3n+\alpha \mathrm{log}_2n+2n^{1+\alpha }+2n\alpha \mathrm{log}_2n`$ (5)
$`+`$ $`2\alpha \mathrm{log}_2nn^\alpha +2\left(\alpha \mathrm{log}_2n\right)^2.`$
If $`\alpha <1`$, the first term dominates as $`n`$ gets large. Therefore the number of row operations is $`O(n^2/\mathrm{log}n)`$. Combining this result with Lemma 1, we have the following theorem.

###### Theorem 3.1

The worst-case length of an $`n`$-wire C-NOT circuit is $`\mathrm{\Theta }(n^2/\mathrm{log}n)`$ gates.

In Equation 5, $`\alpha `$ can be chosen to be arbitrarily close to 1. In the limit, the multiplicative constant in the $`O(n^2/\mathrm{log}n)`$ expression becomes $`1`$ (assuming logs are taken base 2). By contrast, the multiplicative constant in the lower bound in Lemma 1 is $`1/2`$.

This algorithm, in addition to generating more efficient circuits than the standard method, is also asymptotically more efficient in terms of run time. The execution time of the algorithm is dominated by the row operations on the matrix, which are each $`O(n)`$. Therefore the overall execution time is $`O(n^3/\mathrm{log}n)`$ compared to $`O(n^3)`$ for standard Gaussian elimination \[8, p. 42\].

Our algorithm is closely related to Kronrod’s Algorithm (also known as “The Four Russians’ Algorithm”) for construction of the transitive closure of a graph . One important difference between the two is that in their case the goal was a fast algorithm for their application, which is only of secondary concern for our application. Our primary goal is an algorithm that produces an efficient circuit synthesis. Generically, our algorithm can be interpreted as producing an efficient elementary matrix decomposition of a binary matrix.

## 4 Simulation Results

Though Algorithm 1 is asymptotically optimal, it would be of interest to know how large $`n`$ must be before the algorithm begins to outperform standard Gaussian elimination. For this purpose we have synthesized linear reversible circuits using both our method and Gaussian elimination for randomly generated non-singular 0-1 matrices. The results of these simulations are summarized in Figure 4. Our algorithm shows an improvement over Gaussian elimination for $`n`$ as small as $`8`$. The length of the circuit synthesis produced by Algorithm 1 is dependent on the choice $`m`$, the size of the column sections. Here we have somewhat arbitrarily chosen $`m=\text{round}((\mathrm{log}_2n)/2)`$. The performance for some values of $`n`$ could be significantly improved by optimizing this choice. This would also smooth out the performance curve in Figure 4 for Algorithm 1.

## 5 Conclusions and Future Work

We have given an algorithm for linear reversible circuit synthesis that is asymptotically optimal in the worst-case. We show that the algorithm is also asymptotically faster than current methods. While our results are primarily asymptotic, simulation results show that even in the finite case our algorithm outperforms the current synthesis method. Applications of our work include circuit synthesis for quantum circuits.

While the primary motivations for the synthesis method we have given here are to provide an asymptotic bound on circuit complexity and a practical method to synthesize efficient circuits, another application is to bounds on circuit complexity for the finite case. In particular, we can use our method to determine an upper bound on the maximum number of gates required to synthesize any $`n`$ wire C-NOT circuit. For this application the particular partitioning of the columns can be very important. For example, much better bounds can be determined if the size of the sections are a function of the location of the section in the matrix. Sections to the left have more rows below the diagonal and therefore should be larger than sections towards the right of the matrix which have fewer rows below the diagonal. An ongoing area of work is determining optimal column partitioning methods.

Our algorithm basically yields an efficient decomposition for matrices with elements in $`𝔽_2`$, and can be generalized in a straightforward manner for matrices over any finite field. The asymptotic size of the generalized decomposition is $`O(n^2/\mathrm{log}_{|F|}n)`$, where $`|F|`$ is the order of the finite field. Our algorithm, particularly in this generalized form, is quite generic and may lend itself to a wide range of other applications. Related algorithms  have applications in finding the transitive closure of a graph, binary matrix multiplication, and pattern matching.

A major area of future work is extending our results to more general reversible circuits, particularly quantum circuits. Currently, there is an asymptotic gap between the best upper and lower bounds on the worst-case circuit complexity both for general classical reversible circuits and quantum circuits. The gap for classical reversible circuits is the same logarithmic factor that previously existed for linear reversible circuits , which suggests it may be possible to extend our methods to this problem.


text:
# Four Numerical Approaches for Solving the Radiative Transfer Equation in Magnetized White-Dwarf Atmospheres

## 1. The Radiative Transfer Equations

The polarization properties of light are described by the four Stokes parameters $`I`$, $`Q`$, $`V`$, and $`U`$. $`I`$ is the intensity. Linear polarization is described by $`Q`$ and $`U`$, and circular polarization by $`V`$. In general, the absorption coefficients $`\kappa _l`$ and $`\kappa _r`$ for left- and right-handed circularly polarized light, and $`\kappa _p`$ for linearly polarized radiation traveling perpendicular to the magnetic field, are not equal at a given wavelength. $`\kappa _l`$, $`\kappa _p`$, and $`\kappa _r`$ correspond to transitions where the magnetic quantum number changes by $`\mathrm{\Delta }m=1,0,`$ and $`+1`$. The absorption coefficients, normalized to the Rosseland mean opacity ($`\kappa _{\mathrm{Ross}}`$), are defined by $`\eta _p=\kappa _p/\kappa _{\mathrm{Ross}}`$, $`\eta _l=\kappa _l/\kappa _{\mathrm{Ross}}`$, and $`\eta _r=\kappa _r/\kappa _{\mathrm{Ross}}`$. These absorption coefficients are combined into $`\eta _I=\frac{1}{2}\eta _p\mathrm{sin}^2\psi +\frac{1}{4}(\eta _l+\eta _r)(1+\mathrm{cos}^2\psi )`$, $`\eta _Q=\frac{1}{2}\eta _p\frac{1}{4}(\eta _l+\eta _r)\mathrm{sin}^2\psi `$, and $`\eta _V=\frac{1}{2}(\eta _r\eta _l)\mathrm{cos}\psi `$. Here $`\psi `$ denotes the angle between the magnetic field direction and the line of sight. When the magneto-optical parameters $`\rho _R`$ (Faraday rotation) and $`\rho _W`$ (Voigt effect) originating from spectral lines and free electrons in a magnetic field (see Jordan et al. 1992) are taken into account, the three radiative transport equations of Unno (1956) expand into four equations (Beckers 1969):

$$\mu \frac{\text{d}I}{\text{d}\tau }=\eta _I(IB)+\eta _QQ+\eta _VV,$$
(1)
$$\mu \frac{\text{d}Q}{\text{d}\tau }=\eta _Q(IB)+\eta _IQ+\rho _RU,$$
(2)
$$\mu \frac{\text{d}U}{\text{d}\tau }=\rho _RQ+\eta _IU\rho _WV,$$
(3)
$$\mu \frac{\text{d}V}{\text{d}\tau }=\eta _V(IB)+\rho _WU+\eta _VV.$$
(4)
The optical depth scale is $`d\tau \kappa _{\mathrm{Ross}}dz^{}`$, and $`\mu \mathrm{cos}\vartheta `$, where $`\vartheta `$ is the angle between the normal to the surface and the line of sight (along which the $`z^{}`$ axis of the local coordinate system on each surface element points). The $`x^{}`$ axis is defined by the direction of the projection of the magnetic field vector onto the plane of the sky. $`y`$ completes the right-handed coordinate system. $`B`$ is the source function (the Planck formula in the case of LTE).

In the case of large Faraday rotation, radiative transfer can be described in the framework of generalized Stokes parameters, with $`I_+`$ being the ordinary mode and $`I_{}`$ being the extraordinary mode. For the ordinary mode the electric field vector, the magnetic field, and the direction of propagation are in one plane; for the extraordinary mode the plane of propagation and the electric field vector is perpendicular to the magnetic field. The phase relation between the two modes is given by $`I_c=\sqrt{I_+I_{}}\mathrm{cos}\delta `$, and $`I_s=\sqrt{I_+I_{}}\mathrm{sin}\delta `$, which correspond to two polarization ellipses, the so-called normal modes.

We have used three different methods to solve the radiative transfer equations formally: (a) The slightly modified Martin und Wickramasinghe (M&W, 1979), (b) an accelerated lambda iteration scheme (ALI), (c) an approximation for large Faraday rotation (APPROX), or (d) an analytical procedure using matrix exponential function (MATEXP).

### 1.1. Algorithm by M&W

The solution is based on the assumption that the source function is linear in the optical depth and that between two successive depth points the Stokes parameters can be described by

$$\left(\begin{array}{c}I\\ Q\\ U\\ V\end{array}\right)=\left(\begin{array}{c}I_a\\ Q_a\\ U_a\\ V_a\end{array}\right)+\left(\begin{array}{c}I_b\\ Q_b\\ U_b\\ V_b\end{array}\right)\tau +\underset{i=4}{\overset{4}{}}\left(\begin{array}{c}I_{ci}\\ Q_{ci}\\ U_{ci}\\ V_{ci}\end{array}\right)\mathrm{exp}(a_i\tau ).$$
(5)
This assumption is inserted into the four radiative transfer equations, using an Unno (1956) condition at the inner boundary. Then the coefficients $`I_a,\mathrm{}V_{ci}`$, as well as $`a_i`$, are evaluated by comparing the constant, linear, and exponential terms on both sides of the equations. One obtains finally a relatively simple scheme, by which the transport equations can be solved from the inner depth points to the outside.

### 1.2. ALI method

This powerful and very flexible method is described in more detail by Deetjen et al. in these proceedings.

### 1.3. APPROX method

This method uses generalized Stokes parameters (Ramaty 1969). Of particular interest is the limiting case of large Faraday rotation, which is always valid in the case of magnetic white dwarfs. This means that the rotation between two depth points is larger than $`2\pi `$, and the phase relation between the polarization ellipses $`I_s`$ and $`I_c`$ can assume arbitrary values; $`I_s=I_c=0`$ on the average. Then the transfer equations decouple (Ramaty 1969):

$$\mu \frac{dI_\pm }{d\tau }=\alpha _\pm \left(\frac{B}{2}I_\pm \right).$$
(6)
According to Väth & Chanmugam (1995) the formal solution (using a linear approximation of the source function) of this equation can be written as

$`I_+(\tau _i)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\alpha _B(1e^{\alpha _+\mathrm{\Delta }\tau /\mu }){\displaystyle \frac{1}{2}}{\displaystyle \frac{\beta _B}{\alpha _+}}\left[\alpha _+(\tau _i\tau _{i+1}e^{\alpha _+\mathrm{\Delta }\tau /\mu })(1e^{\alpha _+\mathrm{\Delta }\tau /\mu })\right]`$
$`+I_+(\tau _{i\pm 1})e^{\alpha _+\mathrm{\Delta }\tau /\mu }`$
$`I_{}(\tau _i)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\alpha _B(1e^{\alpha _{}\mathrm{\Delta }\tau /\mu }){\displaystyle \frac{1}{2}}{\displaystyle \frac{\beta _B}{\alpha _{}}}\left[\alpha _+(\tau _i\tau _{i+1}e^{\alpha _{}\mathrm{\Delta }\tau /\mu })(1e^{\alpha _{}\mathrm{\Delta }\tau /\mu })\right]`$
$`+I_{}(\tau _{i\pm 1})e^{\alpha _{}\mathrm{\Delta }\tau /\mu }`$
### 1.4. MATEXP method

A direct analytical way to solve the radiative transfer equations uses matrix exponential functions. After Dittmann (1995) the radiative transfer equation $`d\stackrel{}{I}/d\tau =𝜼(\stackrel{}{I}\stackrel{}{S})`$, where $`\stackrel{}{I}(I,Q,V,U)`$, can be written as

$$\stackrel{}{I}(0)=\mathrm{exp}(𝜼\tau )\left[\stackrel{}{I}(\tau )𝜼_\tau ^0\mathrm{exp}(𝜼s)\stackrel{}{S}(s)𝑑s\right]$$
(7)
where the matrix exponential function is defined as $`e^𝜼_{n=0}^{\mathrm{}}𝜼^n/n!`$. The infinite series can be rewritten as $`e^{𝜼\tau }=\underset{i=0}{\overset{3}{}}c_i(\tau )𝐌^i`$. The matrix $`𝐌`$ is given by $`𝐌=𝜼\eta _I𝐄`$. We assume a linear source function between successive depth points: $`S=(\alpha +\beta \tau ,0,0,0)`$. If we define $`x=\frac{1}{2}\sqrt{2(\eta _5^2+\rho _3^2)}`$, $`y=\frac{1}{2}\sqrt{2(\eta _5^2\rho _3^2)}`$, $`z=iy`$, $`\rho _3^2=\sqrt{\eta _5^4+4(\eta _Q\rho _W+\eta _V\rho _R)^2}`$, $`\eta _5=\sqrt{\eta _Q^2+\eta _V^2\rho _R^2\rho _W^2}`$ and take into account that $`y`$ can be become imaginary if magneto-optical parameters are present (a case that has been neglected by Dittmann 1995), then the coefficients are given by:

$`c_0`$ $`=`$ $`{\displaystyle \frac{e^{\eta _I\tau }}{\rho _3^2}}\left[x^2\mathrm{cos}(z\tau )y^2\mathrm{cosh}(x\tau )\right]`$ (8)
$`c_1`$ $`=`$ $`{\displaystyle \frac{e^{\eta _I\tau }}{\rho _3^2}}\left[{\displaystyle \frac{x^2}{z}}\mathrm{sin}(z\tau ){\displaystyle \frac{y^2}{x}}\mathrm{sinh}(x\tau )\right]`$ (9)
$`c_2`$ $`=`$ $`{\displaystyle \frac{e^{\eta _I\tau }}{\rho _3^2}}\left[\mathrm{cosh}(x\tau )\mathrm{cos}(z\tau )\right]`$ (10)
$`c_3`$ $`=`$ $`{\displaystyle \frac{e^{\eta _I\tau }}{\rho _3^2}}\left[{\displaystyle \frac{1}{x}}\mathrm{sinh}(x\tau ){\displaystyle \frac{1}{z}}\mathrm{sin}(z\tau )\right]`$ (11)
For large magneto-optical parameters $`lim_{x0}\mathrm{sinh}(x\tau )/x=\tau `$, and the expressions can be somewhat simplified.

## 2. Conclusion

In all tests, the resulting spectra and polarization turned out to be very similar for the four different approaches. Small numerical instabilities in the case of the M&W and the MATEXP algorithms lead to small oscillations of the Stokes parameters (particularly $`Q`$ and $`V`$) in the vicinity of spectral lines; thsese are caused by the terms containing $`\mathrm{sin}`$, $`\mathrm{cos}`$, $`\mathrm{sinh}`$, or $`\mathrm{cosh}`$. The ALI method, having the potential to account for NLTE and scattering is numerically very stable. However, the APPROX method is about a factor of 11 faster and very stable and should therefore be preferred for practical purposes. For final fits to observations it is always possible to use the ALI method for comparison.

## References

Beckers, J. M. 1969, Sol.Phys., 9, 372

Dittmann, O. 1995, PhD thesis, Heidelberg

Jordan S., 1992, A&A, 265, 570

Jordan, S. 2001, in ASP Conf. Ser. Vol. 226 (San Francisco: ASP), 269

Martin, B., Wickramasinghe D., 1979, MNRAS, 189, 883

Ramaty, R. 1969, ApJ, 158, 753

Unno, W. 1956, PASJ, 8, 108

Väth, H., Chanmugam, G. 1995, ApJS, 98, 295

Wickramasinghe, D. T. & Ferrario, L. 2000, PASP, 112, 873


text:
# Noncommuting Coordinates in the Landau problem

## 1 Introduction

Because of its relevance to string theory, the idea of noncommuting spatial coordinates has gained much attention recently, but the idea actually predates string theory. Coordinate noncommutativity, defined by the equation

$$[x^i,x^j]=i\theta ^{ij}$$
(1)
where $`\theta ^{ij}`$ is a constant, anti-symmetric two-index object, implies a coordinate uncertainty relation, resulting in a discretization of space itself (similar to the familiar quantum phase space) and the elimination of spatial singularities. The idea was suggested by Heisenberg in the 1930s as a way to eliminate divergences in quantum field theory arising from the assumption of point interactions between fields and matter, and the first paper on the subject appeared some years later . Since then, much attention has been given to the study of quantum field theory on noncommutative spaces. Because quantum mechanics is just the one-particle nonrelativistic sector of quantum field theory, it is also relevant to study the quantum mechanics of particles in noncommutative spaces, and to understand the transition between the commutative and noncommutative regimes.

There exists a well known phenomenological realization of noncommuting coordinates in the realm of quantum mechanics: a charged particle in an external magnetic field so strong that projection to the lowest Landau level is justified.A charged particle in an external magnetic field is effectively confined to a two-dimensional space perpendicular to the field, which becomes noncommuting when motion is projected onto the lowest Landau level. It is interesting to note that this example is similar to the instances of noncommuting coordinates which arise in string theory in that both are characterized by the presence of a strong background magnetic-like field. Because the lowest Landau level is the only physically realized example of noncommuting coordinates, it should be useful to our understanding of noncommutative spaces in physics to know precisely how noncommutativity arises in this system. To this end, a quantum solution to the Landau problem is reviewed, and then the coordinate commutator is calculated after a projection to a truncated space of Landau levels.

## 2 The Landau Problem

We consider a charged ($`e`$) and massive ($`m`$) particle in a constant magnetic field ($`B`$), which we chose to point along the $`z`$ axis without loss of generality. This is called the Landau problem because the eigenstates and eigenvalues were investigated for the first time by Landau. The Lagrangian describing the particle’s motion is

$$L=\frac{1}{2}m(\dot{x}^2+\dot{y}^2+\dot{z}^2)+\frac{e}{c}(\dot{x}A_x+\dot{y}A_y+\dot{z}A_z).$$
(2)
¿From this, we determine the Hamiltonian for the particle in the usual way, with the result

$$H=\frac{1}{2m}(\stackrel{}{p}\frac{e}{c}\stackrel{}{A})^2$$
(3)
where the canonical momentum $`\stackrel{}{p}`$ of the particle is no longer the usual $`m\stackrel{}{v}`$ but is equal to $`m\stackrel{}{v}+\frac{e}{c}\stackrel{}{A}`$. The Hamiltonian is thus simply $`H=\frac{1}{2}m\stackrel{}{v}^2`$, which is what we expect since the magnetic field does no work and thus cannot contribute to the energy.

To solve the eigenstates and eigenenergies of this problem, we must prescribe a vector potential $`\stackrel{}{A}`$, satisfying $`\times \stackrel{}{A}=B\widehat{z}`$. We chose $`\stackrel{}{A}=(0,xB,0)`$ Substituting this gauge choice into the Hamiltonian, we have

$$H=\frac{1}{2m}(p_x^2+p_z^2+(p_y\frac{exB}{c})^2).$$
(4)
Because $`H`$ commutes with both $`p_z`$ and with $`p_y`$, we can write the eigenvalue equations $`p_z|\mathrm{\Psi }=\mathrm{}k_z|\mathrm{\Psi }`$ and $`p_y|\mathrm{\Psi }=\mathrm{}k_y|\mathrm{\Psi }`$, where $`|\mathrm{\Psi }`$ denotes the eigenstates of the full Hamiltonian, and $`\mathrm{}k_z`$ and $`\mathrm{}k_y`$ are the eigenvalues of their respective momentum operators. Using these eigenvalues, the Hamiltonian can be written as

$$H=\frac{1}{2m}(\mathrm{}k_z)^2+\frac{1}{2m}\left(p_x^2+\left(\frac{eB}{c}\right)^2\left(x\frac{c\mathrm{}k_y}{eB}\right)^2\right).$$
(5)
¿From now on, motion in the $`z`$ direction is suppressed since it is not quantized, and only planar motion in the $`x`$-$`y`$ plane is retained. The second term in the Hamiltonian is just a shifted harmonic oscillator with angular frequency $`\omega =\frac{eB}{cm}`$. The eigenstates $`|\mathrm{\Psi }`$ are thus labeled by the number of oscillator quanta $`n`$ and by $`k_y`$. The associated eigenenergies are

$$E_n=\mathrm{}\frac{eB}{mc}\left(n+\frac{1}{2}\right).$$
(6)
$`k_y`$ is free, so each $`n`$ indexes an energy eigenstate, called a Landau level, which is infinitely degeneracy in $`k_y`$. Note that adjacent Landau levels are separated by energy $`\mathrm{}\frac{eB}{mc}`$. The eigenfuctions in coordinate space are

$$x,y|n,k_y=\frac{1}{\sqrt{2\pi \mathrm{}}}e^{ik_yy}\varphi _n(x\frac{c\mathrm{}k_y}{eB})$$
(7)
where $`\varphi _n`$ are the normalized harmonic oscillator wavefunctions. ¿From now on, we let $`k\mathrm{}k_y`$ and label the eigenstates $`|n,k`$.

## 3 Projection to the lowest Landau level

Since the separation between the states $`|n,k`$ is $`𝒪(B/m)`$, if the magnetic field is strong, only the lowest Landau level $`|0,k`$ is relevant. The higher states are essentially decoupled to infinity. The large $`B`$ limit is equivalent to the small $`m`$ limit. So the Lagrangian can be modified to describe only the lowest Landau level by setting $`m`$ to zero in (2). With this modification, and the addition of a potential $`V(x,y)`$ to represent impurities in the plane, the Lagrangian becomes

$$L_{lLl}=\frac{e}{c}Bx\dot{y}V(x,y).$$
(8)
This has the same form as $`L=p\dot{q}H(p,q)`$, and thus we recognize $`\frac{e}{c}Bx`$ and $`y`$ as canonical conjugates, with the corresponding commutator:

$`[{\displaystyle \frac{eB}{c}}x,y]`$ $`=`$ $`i\mathrm{}`$ (9)
$`[x,y]`$ $`=`$ $`i{\displaystyle \frac{\mathrm{}c}{eB}}`$ (10)
So we see that by restricting the particle to the first Landau level, the space it moves in no longer obeys the standard Heisenberg algebra. This is called the “Peierls substitution” .

We can verify this result in a less heuristic fashion, by directly calculating the matrix elements of the coordinate commutator .

$`n,k|[x,y]|n^{},k^{}`$ $`=`$ $`n,k|xy|n^{},k^{}n,k|yx|n^{},k^{}`$ (11)
$`=`$ $`n,k|xy|n^{},k^{}n^{},k^{}|yx|n,k^{}`$ (12)
$`=`$ $`f(nk,n^{}k^{})f^{}(n^{}k^{},nk)`$ (13)
where

$$f(nk,n^{}k^{})=n,k|xy|n^{},k^{}$$
(14)
We evaluate $`f`$ by inserting a complete set of states in the product $`xy`$:

$$f(nk,n^{}k^{})=\underset{m,q}{}n,k|x|m,qm,q|y|n^{},k^{}.$$
(15)
If we restrict the particle’s world to include only the lowest Landau level, then we should only include the lowest Landau level (with its infinite degeneracy) in the intermediate state sum of this calculation. Additionally, since we’re pretending that the world is confined to the lowest Landau level, it only makes sense to calculate this matrix element for $`n=n^{}=0`$. To evaluate that element, we must calculate

$$f(0k,0k^{})=𝑑q0,k|x|0,q0,q|y|0,k^{}.$$
(16)
To evaluate each matrix element, we explicitly write the plane wave functions and the harmonic oscillator wavefunctions. These expressions are readily evaluated by integration over the $`x`$-$`y`$ plane.

$`0,k|x|0,q`$ $`=`$ $`{\displaystyle 𝑑x𝑑y\frac{1}{\sqrt{2\pi \mathrm{}}}e^{iky/\mathrm{}}(\frac{m\omega }{\pi \mathrm{}})^{1/4}e^{\frac{m\omega }{2\mathrm{}}(x\frac{ck}{eB})^2}x}`$ (17)
$`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{}}}}e^{iqy/\mathrm{}}({\displaystyle \frac{m\omega }{\pi \mathrm{}}})^{1/4}e^{\frac{m\omega }{2\mathrm{}}(x\frac{cq}{eB})^2}`$
$`=`$ $`\delta (kq){\displaystyle \frac{cq}{eB}}`$ (18)
$`0,q|y|0,k^{}`$ $`=`$ $`{\displaystyle 𝑑x𝑑y\frac{1}{\sqrt{2\pi \mathrm{}}}e^{iqy/\mathrm{}}(\frac{m\omega }{\pi \mathrm{}})^{1/4}e^{\frac{m\omega }{2\mathrm{}}(x\frac{cq}{eB})^2}y}`$ (19)
$`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{}}}}e^{ik^{}y/\mathrm{}}({\displaystyle \frac{m\omega }{\pi \mathrm{}}})^{1/4}e^{\frac{m\omega }{2\mathrm{}}(x\frac{ck^{}}{eB})^2}`$
$`=`$ $`i\mathrm{}\delta ^{}(qk^{})`$ (20)
where the prime denotes differentiation with respect to the argument.

The commutator can then be obtained:

$$0,k|[x,y]|0,k^{}=i\frac{\mathrm{}c}{eB}0,k|0,k^{}$$
(21)
which is consistent with (10).

## 4 Projection to higher Landau levels

One could now pretend that the world is restricted to the lowest $`N+1`$ Landau levels ($`n=0,1,\mathrm{},N`$). In this case, we cannot obtain the commutator $`[x,y]`$ heuristically by modifying the Lagrangian and reading off a canonical pair, but our method of explicitly calculating the relevant matrix elements is still valid. To evaluate $`f`$ we include the lowest $`N+1`$ levels (each with its infinite degeneracy) in the intermediate state sum.

$$f(nk,n^{}k^{})=\underset{m=0}{\overset{N}{}}𝑑qn,k|x|m,qm,q|y|n^{},k^{}.$$
(22)
This is defined for $`n`$ and $`n^{}`$ less than or equal to $`N`$. Each matrix element is evaluated by explicitly writing the coordinate space plane wavefunctions, and representing the harmonic oscillator wavefunctions as $`\varphi _n(x)`$. The expressions are simplified by exploiting the orthonormality of the harmonic oscillator wavefunctions and integrating over the plane.

$`n,k|x|m,q`$ $`=`$ $`{\displaystyle 𝑑x𝑑y\frac{1}{\sqrt{2\pi \mathrm{}}}e^{iky/\mathrm{}}\varphi _n(x\frac{ck}{eB})x}`$ (23)
$`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{}}}}e^{iqy/\mathrm{}}\varphi _m(x{\displaystyle \frac{cq}{eB}})`$
$`=`$ $`\delta (qk)\delta _{nm}{\displaystyle \frac{ck}{eB}}+\delta (qk)n|x|m`$ (24)
$`m,q|y|n^{},k^{}`$ $`=`$ $`{\displaystyle 𝑑x𝑑y\frac{1}{\sqrt{2\pi \mathrm{}}}e^{iqy/\mathrm{}}\varphi _n(x\frac{cq}{eB})y}`$ (25)
$`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{}}}}e^{ik^{}y/\mathrm{}}\varphi _n^{}(x{\displaystyle \frac{ck^{}}{eB}})`$
$`=`$ $`\mathrm{}i\delta ^{}(qk^{})\delta _{mn^{}}{\displaystyle \frac{c}{eB}}\delta (qk^{})n^{}|p|m`$ (26)
By evaluating the sums and integrals of (22), one can find that the matrix elements of the coordinate commutator vanish unless $`n=n^{}=N`$. For this case, we get the result

$$N,k|[x,y]|N,k^{}=i\mathrm{}\frac{c}{eB}(N+1)\delta (kk^{}).$$
(27)
So, for example, the matrix representation of the coordinate commutator for a particle confined to the two lowest Landau levels is

$$[x,y]=\left(\begin{array}{cc}0& 0\\ 0& 2i\mathrm{}\frac{c}{eB}\end{array}\right)$$
(28)
If we expanded our world to include another Landau level, we would find that the matrix element $`2,k|[x,y]|2,k^{}`$ vanishes because of the additional $`m=2`$ term in the sum in (22). Thus we have shown that the phenomena of noncommutative space is not specific to the lowest Landau level but can be obtained by projecting to an arbitrary finite number of Landau levels.

## 5 Motion in the symmetric gauge

By choosing a different gauge, we can verify this result, as well as come upon it much more elegantly. Substituting the symmetric gauge $`\stackrel{}{A}=\frac{B}{2}(y,x)`$ into the Hamiltonian describing the planar motion, we find

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2m}}\left((p_x+{\displaystyle \frac{eB}{2c}}y)^2+(p_y{\displaystyle \frac{eB}{2c}}x)^2\right)`$ (29)
$`=`$ $`{\displaystyle \frac{1}{2m}}(p_x^2+p_y^2)+{\displaystyle \frac{1}{2}}m\left({\displaystyle \frac{eB}{2mc}}\right)^2(x^2+y^2){\displaystyle \frac{eB}{2mc}}L,`$ (30)
where $`L=xp_yyp_x`$ is the angular momentum in the $`x`$-$`y`$ plane. In this gauge, the system looks like a two-dimensional harmonic oscillator with an additional interaction $`\frac{eB}{2mc}L`$. Since we are dealing here only with harmonic oscillator wavefunctions and not plane waves, it is clear that this problem will most easily solved by introducing harmonic oscillator creation and annihilation operators. To this end, we define:

$`a`$ $`=`$ $`{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{eB}{2c\mathrm{}}}}(xiy)+{\displaystyle \frac{i}{2}}\sqrt{{\displaystyle \frac{2c}{eB\mathrm{}}}}(p_xip_y)`$ (31)
$`b`$ $`=`$ $`{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{eB}{2c\mathrm{}}}}(x+iy)+{\displaystyle \frac{i}{2}}\sqrt{{\displaystyle \frac{2c}{eB\mathrm{}}}}(p_x+ip_y).`$ (32)
satisfying the following commutation relations

$`[a,a^{}]=[b,b^{}]=1`$ (33)
with all other commutators vanishing. ($`a,a^{}`$) and ($`b,b^{}`$) are two pairs of independent harmonic oscillator operators. We can now write $`L`$ and $`H`$ in terms of these operators:

$`L`$ $`=`$ $`\mathrm{}(a^{}ab^{}b)`$ (34)
$`H`$ $`=`$ $`\mathrm{}{\displaystyle \frac{eB}{2mc}}(a^{}a+b^{}b+1)\mathrm{}{\displaystyle \frac{eB}{2mc}}(a^{}ab^{}b)`$ (35)
$`=`$ $`\mathrm{}{\displaystyle \frac{eB}{mc}}(b^{}b+{\displaystyle \frac{1}{2}}).`$ (36)
Eigenstates of the Hamiltonian are labeled by the number $`j`$ of excitation quanta of the oscillator $`a`$, and the number $`n`$ of excitation quanta of the oscillator $`b`$:

$`a^{}a|n,j=j|n,j,`$ (37)
$`b^{}b|n,j=n|n,j.`$ (38)
Both $`j`$ and $`n`$ can take on any non-negative integer value. Only the $`n`$ quanta contribute to the energy, so the Landau levels, labeled by $`n`$, are each infinitely degenerate in $`j`$.

To evaluate the commutator $`[x,y]`$ using this basis, we first rewrite $`x`$ and $`y`$ in terms of $`\alpha `$ and $`\alpha ^{}`$, where $`\alpha =a+b^{}`$

$`x`$ $`=`$ $`\sqrt{{\displaystyle \frac{\mathrm{}c}{2eB}}}(\alpha +\alpha ^{}),`$ (39)
$`y`$ $`=`$ $`i\sqrt{{\displaystyle \frac{\mathrm{}c}{2eB}}}(\alpha \alpha ^{}).`$ (40)
Then the coordinate commutator is

$`[x,y]`$ $`=`$ $`{\displaystyle \frac{i}{2}}{\displaystyle \frac{\mathrm{}c}{eB}}[\alpha +\alpha ^{},\alpha \alpha ^{}]`$ (41)
$`=`$ $`i{\displaystyle \frac{\mathrm{}c}{eB}}[\alpha ,\alpha ^{}].`$ (42)
It we evaluate this exactly, we find that the commutator vanishes as expected.

Now using this new basis we can explicitly evaluate matrix elements for an incomplete state space. We will denote the eigenstates $`|n,j`$ explicitly as product states: $`|n|j`$. To determine the commutator, we need to evaluate the matrix elements:

$$n|j|\alpha \alpha ^{}|j^{}|n^{}n|j|\alpha ^{}\alpha |j^{}|n^{}(1)(2).$$
(43)
We evaluate each term separately by inserting intermediate states.

$`(1)`$ $`=`$ $`{\displaystyle \underset{m,l=0}{\overset{m=N,l=\mathrm{}}{}}}n|j|\alpha |l|mm|l|\alpha ^{}|j^{}|n^{}`$ (44)
$`(2)`$ $`=`$ $`{\displaystyle \underset{m,l=0}{\overset{m=N,l=\mathrm{}}{}}}n|j|\alpha ^{}|l|mm|l|\alpha |j^{}|n^{}`$ (45)
Summing over $`l`$, one readily finds

$`(1)(2)=\delta _{nn^{}}\delta _{jj^{}}(1+{\displaystyle \underset{m=1}{\overset{N+1}{}}}m\delta _{nm}{\displaystyle \underset{m=0}{\overset{N1}{}}}(m+1)\delta _{nm})`$ (46)
Since $`m`$ never gets larger than $`N`$, we can write this as

$`=`$ $`\delta _{nn^{}}\delta _{jj^{}}(1+{\displaystyle \underset{m=1}{\overset{N}{}}}m\delta _{nm}{\displaystyle \underset{m=0}{\overset{N1}{}}}(m+1)\delta _{nm})`$ (47)
$`=`$ $`\delta _{nn^{}}\delta _{jj^{}}(1+N\delta _{nN}{\displaystyle \underset{m=0}{\overset{N1}{}}}\delta _{nm})`$ (48)
Let us analyze this expression. It clearly vanishes for off diagonal elements. For $`n=n^{}N`$, it also vanishes because the sum will yield unity since $`n`$ is necessarily less than $`N`$. For $`n=n^{}=N`$, we obtain our earlier result:

$`N|j|[x,y]|j^{}|N=i{\displaystyle \frac{\mathrm{}c}{eB}}(N+1)j|j^{}.`$ (49)
Suppose we had let $`N\mathrm{}`$, then (46) would have been

$`=`$ $`\delta _{nn^{}}\delta _{jj^{}}(1+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}m\delta _{nm}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(m+1)\delta _{nm})`$ (50)
$`=`$ $`\delta _{nn^{}}\delta _{jj^{}}(1\delta _{n0}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\delta _{nm})`$ (51)
which clearly vanishes for any choice of $`n`$. Noncommutativity in the Landau problem is clearly associated with a truncated state space.

## 6 Acknowledgements

Thanks to Roman Jackiw for giving me this problem to solve.


text:
# Dust emission from young outflows: the case of L 1157

## 1 Introduction

Mm and submm mapping of the dust thermal continuum emission has proven to be a powerful tool to study the density structure of low-mass protostellar condensations at different evolutionary stages (e.g. Motte et al. 1998; Chandler & Richer 2000; Hogerheijde & Sandell 2000; Shirley et al. 2000; Motte & André 2000). However, virtually all protostars are driving powerful molecular outflows (Bachiller 1996), which are much more extended than the protostellar envelope size (e.g. André et al. 1999). Strong interaction between the outflow and the envelope is observed in a few cases (e.g. L 1527, Motte & André 2000). But besides such local perturbations, there is increasing observational evidence that the outflow itself can be responsible for (sub)mm continuum emission also at long distances from the driving source (Shirley et al. 2000; Chini et al. 2001). Since such continuum emission can be a valuable tracer of the outflow material, its detailed origin and structure clearly deserves further study.

The outflow located in the L 1157 molecular cloud, at a distance of 440 pc, is one of the few which have been detected so far in the continuum submm and mm emission (Shirley et al. 2000; Chini et al. 2001). This outflow is driven by an extremely young Class 0 source (Gueth et al. 1997). The southern lobe consists of two large overlapping cavities, each one ended by a strong bow shock, which are misaligned due to the probable precession of the ejection axis (Gueth et al. 1996, 1998; Bachiller et al. 2001). The two shocks exhibit a rich molecular composition (Tafalla & Bachiller 1995; Zhang et al. 1995; Avery & Chiao 1996; Bachiller & Pérez-Gutiérrez 1997) which makes L 1157 the prototype of chemically active outflows (Bachiller et al. 2001). In this letter, we report new high-sensitivity bolometer observations of L 1157 at mm wavelengths, which, in combination with previously published maps, provides basic information on the continuum emitting regions.

## 2 Observations

The observations were carried out in 1997 February and 1999 December, with the IRAM 30–m telescope. The detector was the MPIfR 1.3 mm bolometer array (MAMBO), which had 19 channels in 1997 and 37 channels in 1999. Each channel has an equivalent bandwidth of $``$70 GHz. 8 individual on-the-fly maps were obtained, with a different scanning direction for each coverage, from nearly perpendicular to almost parallel to the flow axis. The antenna beam width was $`12^{\prime \prime }`$. Skydips were observed every one to two hours to monitor the atmospheric zenith optical depth, which was found to be 0.1 to 0.4. Flux calibration was achieved by observing Uranus. The final flux accuracy is estimated to be $``$30%. The data were reduced with the IRAM package (NIC), using the EKH restoration algorithm. Fig. 1a presents the resulting image. The rms at the map center is $``$3.5 mJy/beam.

Beltrán et al. (2003) obtained IRAM Plateau de Bure interferometric observations of the 1.3 mm continuum emission of the L 1157 protostar. The bolometer map presented in this paper was used to derive the short-spacing information and thereby complement the interferometric dataset. The resulting map, which covers only the central region of the bolometer image, is shown in Fig. 1b.

## 3 Results

### 3.1 Four continuum components

The data presented in Fig. 1 show that the mm continuum emission in L 1157 arises from four different components. Their properties are summarized in Table 1.

Central source. The interferometric observations (Fig. 1b) reveal the existence of a compact ($`<1^{\prime \prime }`$) source, of $``$80 mJy, located at 20h39m06s.24, 680215′′.6 (J2000). Obviously, this structure is closely associated to the protostar itself, and it is tempting to identify it with the thermal dust emission of a circumstellar disk. We refer to Beltrán et al. (2003) for a more detailed discussion.

Envelope. The central source is surrounded by an extended emission, that originates from the infalling (Gueth et al. 1997; Mardones et al. 1997) protostellar envelope. On the bolometer map, the peak flux density is 500 mJy/beam, and the flux integrated in a circle of 40′′ diameter is 960 mJy. Assuming a typical dust temperature of 30 K, the corresponding mass would be $``$2.1 $`M`$$`_{\mathrm{}}`$ (Table 1). At the resolution of these observations ($`12^{\prime \prime }`$), the envelope is smooth and symmetric (the contour at half intensity is nearly round and has a diameter of $``$13′′, cf. Fig. 1a). However, the interferometric map (Fig. 1b) shows that the perturbation of the envelope by the outflow affects the internal structure of the envelope, which is no longer isotropic at that resolution, but extended along the flow axis.

Extended envelope. The L 1157 protostellar envelope is surrounded by an even more diffuse emission extended roughly east-west, i.e. approximately perpendicular to the flow axis (Fig. 1a). As shown in Fig. 2, this feature is also seen in the 850 $`\mu `$m continuum map obtained by Shirley et al. (2000), as well as in the lower resolution NH3 (1,1) map obtained by Bachiller et al. (1993). It is also detected in quiescent emission of 13CO (Gueth et al. 1997), DCO+, and N2H+ (Bachiller et al. 2001). Very likely this structure is a dense remnant of the molecular cloud in which the protostar was formed. The size of this envelope in the mm map is about 140$`{}_{}{}^{\prime \prime }\times `$60′′, its flux density is roughly constant, $``$12 mJy/beam, and its total flux is 320$`\pm `$50 mJy (after removing the 960 mJy corresponding to the central source and its envelope). If the dust temperature is 13 K (as suggested by the NH3 observations of Bachiller et al. 1993), the mass of this extended structure is $``$2.1 $`M`$$`_{\mathrm{}}`$.

Outflow. Besides the emission associated with the protostar vicinity, Fig. 1a shows that there is also an emission extending toward the south-east. A comparison with the CO (2–1) map of Bachiller et al. (2001), which has the same angular resolution, clearly shows that this emission is associated to the outflow (see Fig. 2): the continuum emission presents two main structures (labeled B1 and B2), whose positions and shapes correspond very well to those of the brightest CO emission. This is further demonstrated by the morphological agreement that exists between the mm continuum map and the NH3 (3,3) emission (Fig. 2) and the high-resolution SiO (2–1) maps of Gueth et al. (1998). The NH3 and SiO emission trace shocks in molecular outflows. B1 does also coincide with the HH 375 Herbig-Haro object (Chini et al. 2001). We thus conclude that significant emission arises from the shocked regions of the L 1157 outflow. The peak flux densities of B1 and B2 are $``$50 and $``$25 mJy/beam, respectively, and their integrated fluxes are 245 and 90 mJy.

### 3.2 Origin of the outflow emission

The broad bolometer bandwidth includes many molecular lines, which are quite strong in L 1157, and could thus create a spurious continuum emission. At the B1 position, the CO (2–1) line-integrated emission is 155 K km-1 (Bachiller et al. 2001), which corresponds to $``$8 mJy/beam in the bolometer band, i.e. 16% of the detected emission. Other lines are weaker, but their sum can however be important: adding the integrated intensities of all lines (but CO) detected by Bachiller & Pérez-Gutiérrez (1997) at B1 yields 110 K km-1. This puts the lines contribution to about 30% of the emission detected by the bolometer. This value could possibly be somewhat larger, since the line survey was not complete and not covering the full bolometer band.

Fig. 2 also shows the SCUBA 850 $`\mu `$m map obtained by Shirley et al. (2000). This map looks very similar to the 1.3 mm observations. From Hirano & Taniguchi (2001), we can roughly estimate the CO (3–2) line-integrated emission to be 200 K km-1 at the B1 position, which corresponds to $``$20% of the bolometer measurement. Other lines must contribute and, as in the 1.3 mm case, could match the CO intensity and thus raise the line contribution to about 40%.

We conclude that the lines contributions to the 1.3 mm and 850 $`\mu `$m bolometer maps are significant, but are unlikely to account for the whole emission, which is thus dominated by continuum emission. A comparison of both maps reveals a spectral index in the range 3 to 5.5 (see below), which is typical for dust thermal emission, and thus suggests that any free-free emission is negligible in these images. Dust is actually thought to play a key role in the chemical processes taking place in shocked outflows (e.g. Schilke et al. 1997).

It is difficult to estimate the temperature of the dust in the outflow: the IRAS measurements (that would be needed to build the SED and then derive a temperature) do not have the angular resolution required to resolve the different components contributing to the L 1157 emission. Tafalla & Bachiller (1995) used multiline NH3 analysis to derive outflowing gas temperatures of 60 to 80 K. These values are upper limits to the temperature of the dust, which has a very short post-shock cooling time (e.g. Clark et al. 1986). Chini et al. (2001) suggested a dust temperature of 27 K for HH 375 (B1). From the measured mm fluxes, assuming optically thin emission and a temperature between 25 and 60 K, we estimate that the mass of the blueshifted shocked regions is $``$0.3 to 0.9 $`M`$$`_{\mathrm{}}`$.

### 3.3 Spectral index

Fig. 3 presents a map of the spectral index $`\alpha `$ computed between 850 $`\mu `$m and 1.3 mm: $`\alpha `$ = log($`S_{850\mu \mathrm{m}}/S_{1.3\mathrm{mm}}`$)$`/`$log($`1300/850`$). This image shows systematic variations across the mapped region. The lowest values of the index, in the range 3 to 4, are obtained around the protostellar condensation. They are in good agreement with the average value $`\alpha `$ = 3.4$`\pm `$0.3 found by Shirley et al. (2000) in a sample of Class 0 sources. The highest values are in the range from 4.5 to 5.5, and are measured towards the outflow region.

Possible artifacts. A calibration error in one or both maps would translate into an incorrect scale for the spectral indexes, but would not introduce a position-dependent effect. However, line contamination of the bolometer measurements in the outflow region can affect the apparent spectral index by introducing an offset to the actual index of the continuum emission. Simple algebra shows that the values discussed in the previous section result in $`\alpha `$ being increased by 0.35 at B1, the position where the line contamination reaches its highest value. The variations of $`\alpha `$ at other positions are expected to be smaller. However, uncertainties on the line contamination, as well as possible calibration errors, make that this estimate of the $`\alpha `$ variation is only a rough approximation. In fact, the offset could be smaller (or even negative), if e.g. line contamination at 1.3 mm was underestimated. Inversely, if the line contribution at 850 $`\mu `$m is higher, the offset can be larger: in order to introduce the observed 1.5 difference between the outflow and the protostellar envelope, a line/continuum ratio at 850 $`\mu `$m of $``$2 is required. This would imply that the 850 $`\mu `$m map is essentially line emission, a possibility which –from the fluxes of the lines observed in the band– seems unlikely.

Dust properties. Besides the possible effects of the line contamination, a significant fraction of the variation in the spectral index has to be traced back to intrinsic properties of the continuum emission. For optically thin thermal dust emission, in the Rayleigh-Jeans limit, $`\alpha =2+\beta `$, where $`\beta `$ is the exponent in the dust opacity law $`\kappa \nu ^\beta `$. Hence, our observations would suggest that $`\beta `$ varies substantially from $``$1.5 around the protostar to $``$3 in the shocked regions111 Note however that the factor $`h\nu `$/$`k`$ is $``$11 K at 1.3 mm and $``$17 K at 850 $`\mu `$m, so the Rayleigh-Jeans approximation (T$``$$`h\nu `$/$`k`$) fails e.g. in the extended envelope, where Tdust$``$13 K (if the dust is at the temperature indicated by the NH3 observations). In that case, $`\alpha =2+\beta +\gamma (T_{\mathrm{dust}})`$, where $`\gamma `$ can easily be derived from the Planck function (see e.g. Visser et al. 1998). At 13 K, $`\gamma 0.6`$, but this correction is only 0.2 at 30 K, and $`0.1`$ for $`T_{\mathrm{dust}}60`$ K. Hence, this temperature effect cannot explain the observed variations of $`\alpha `$ across the L 1157 area, and one has to invoke variations of $`\beta `$.. The frequency variation of dust opacity is not well known, and depends on dust characteristics (see e.g. Visser et al. 1998, and references therein). Interestingly, a larger $`\beta `$ – as suggested in the flow – points toward smaller dust grains (e.g. Krügel & Siebenmorgen 1994) and/or different chemical properties (e.g. Pollack et al. 1994). Since the mm continuum emission in L 1157 arises from shocked regions, as shown by the morphological coincidences with shock tracers such as NH3 and SiO, this result would be consistent with models predicting mantle evaporation (hence the rich chemistry observed in the shocks of L 1157) and even grain destruction (which may be the only way to release silicon and produce strong SiO emission, Schilke et al. 1997).

## 4 Conclusion

The observations presented in this letter demonstrate that the mm continuum emission in L 1157 originates from four distinct components: the circumstellar disk, the infalling envelope, the dense remnant of the cloud in which the protostar has formed, and the outflow. The emission from the outflow represents $``$20% of the total flux at 1.3 mm, and is dominated by dust thermal emission.

The spectral index computed between 1.3 mm and 850 $`\mu `$m shows a significant difference between the protostellar envelope and the outflow emission. This could be related to an important line contamination of the 850 $`\mu `$m map, and/or to variations of the dust opacity law index. In the latter case, this may be an indication of different dust properties, possibly of smaller grains in the outflow. Observations of additional bipolar flows are needed to check if such variations in the spectral index are a common phenomenon, and to investigate their origin.

Finally, we caution that large-scale mm and submm continuum mapping of molecular clouds (e.g. Motte et al. 1998) or of the vicinity of protostellar objects (e.g. Shirley et al. 2000) can potentially be confused by the presence of outflow emission, such as that from L 1157. Some emission knots could easily be misinterpreted as, e.g., starless clumps. The observations reported here demonstrate that emission from bipolar outflows can be important, and must be taken into account when interpreting the structure of mm maps of star forming regions.

###### Acknowledgements.

The authors are grateful to the IRAM staff at Pico Veleta, to R. Neri for valuable help with the observations, and to Y.L. Shirley and N.J. Evans for providing a digital version of the 850 $`\mu `$m map. This research have been partially supported by Spanish DGES grant AYA2000-927.