text:
# Controlling a Telescope Chopping Secondary Mirror Assembly Using a Signal Deconvolution Technique

## I Introduction

Chopping scans are widely used in radioastronomy as they provide an efficient way to reduce the adverse effects that any instabilities present either in the sky signal or some telescope equipment can have on the detection of weak signals. A chopping scan is defined as a mode of observation where the telescope’s secondary mirror is rotated back and forth through some angle and where the signals from both “end” positions are integrated separately and subtracted from each other. This mode is to be compared with the so-called ON-OFF position (beam switching) scan where the telescope actually moves back and forth from one end position to the other. Because of the much greater speed at which the secondary mirror can move compared to the telescope, the signal subtraction happens much faster and thus an increase in the ability to detect weak signals. By moving or chopping the secondary mirror even at a relatively low frequency (e.g., 1 Hz) one can obtain a significant improvement in baseline quality when compared to a typical beam switch. In what will follow, the secondary mirror displacements are in units of arcseconds as measured on the sky.

At the Caltech Submillimeter Observatory (CSO) a chopping secondary mirror assembly was installed in 1994 and has since been used both for heterodyne receivers and bolometer cameras (e.g., SHARC and HERTZ) observations. It is composed, in part, of a carbon fiber mirror mounted on a DC brushless motor along with a system of counterweights which greatly reduces the amount of vibration noise transmitted to the observatory’s receivers or cameras. The huge advantage that this vibration suppression technique brings, for the detection of weak signals comes, however, at the cost of an increase in the inertia of the chopper assembly which causes a reduction in the speed and an increase in the settling time in the response of the system.

We show in Figure 1 a block diagram of the chopping secondary system as it has been used since its installation at the CSO. Once the user of the telescope has selected a chop throw amplitude (in arcseconds) and frequency, a square wave is sent to the input of a typical Proportional-Integral-Derivative (PID) electronic controller (Ogata, 1970) where it is compared to the position signal of the mirror (obtained through a Linear Variable Differential Transformer or LVDT). The processed error signal is then sent to a power amplifier which feeds the motor and thus continuously re-positions the mirror while the PID controller acts to minimize the error signal.

Because of the relatively slow response of the chopper assembly, and the non-linearities inherent to the system (see section IV), the parameters of the PID controller cannot be held fixed at a given set of values but have to be adjusted by the user of the telescope for different chop throw amplitudes and frequencies. Although this does not present a problem in principle it has been the experience that the tuning of the controller’s parameters can sometimes be a time consuming effort that reduces the efficiency of the observatory. Also, since, as will be shown later, the response time of the assembly is of the order of the chopping period (or more) it is often quite difficult to find the appropriate set of parameters that will give optimum results. Too often, the outcome of such exercise is a reduction in the performance of the chopping assembly; both in its settling time and positioning accuracy.

In the following sections of this paper we will demonstrate how a signal processing method based on the Lucy deconvolution technique (Lucy, 1974) was implemented at the CSO to solve this problem and provide a system that requires no intervention from the telescope user, while keeping hardware changes to a minimum. We will start in the next section with a brief exposition of the set of equations that govern the Lucy deconvolution technique followed by a presentation of the new chopping secondary assembly (section III). We will finish by showing how the deconvolution technique was implemented, along with the needed modifications, and by presenting some results obtained so far.

## II Lucy’s Deconvolution Technique

An iterative method for signal deconvolution based on the Bayes rule for conditional probabilities was introduced by Lucy (1974) and has been successfully used in astronomy for the processing and extraction of precise photometric information from originally blurred images taken under average seeing conditions (see for example Houde & Racine (1994)).

Limiting ourselves to a one-dimensional problem, the set of equations governing Lucy’s technique is relatively simple. Denoting by $`r\left(t\right)`$ and $`s\left(t\right)`$ the input and output signals of a linear system, respectively, we know that they are related to each other through the transfer function $`h\left(t\right)`$ of this same system by the following convolution integral

$$s\left(t\right)=r\left(\tau \right)h\left(t\tau \right)𝑑\tau ,$$
(1)
where the limits of integration in equation (1), and in all of the integrals that will follow, are from $`\mathrm{}`$ to $`+\mathrm{}`$.

The goal of a deconvolution technique is to invert equation (1) and express $`r\left(\tau \right)`$ as a function of $`s\left(t\right)`$ using a new function $`g\left(\tau t\right)`$ as follows

$$r\left(\tau \right)=s\left(t\right)g\left(\tau t\right)𝑑t.$$
(2)
Lucy’s idea was to liken the (reversed) time shifted transfer function $`h\left(t\tau \right)`$ to a Bayes density function of conditional probability. In doing so, the new function $`g\left(\tau t\right)`$ can be interpreted as a new density function and readily determined using the Bayes rule by (Haykin, 1978)

$$g\left(\tau t\right)=\frac{r\left(\tau \right)h\left(t\tau \right)}{s\left(t\right)},$$
(3)
or alternatively

$$g\left(\tau t\right)=\frac{r\left(\tau \right)h\left(t\tau \right)}{r\left(\lambda \right)h\left(t\lambda \right)𝑑\lambda }.$$

Evidently, it is impossible to directly determine $`g\left(\tau t\right)`$ with equation (3) since it is expressed as a function of $`r\left(\tau \right)`$ which is the unknown that we are trying to evaluate. But the form of equations (1), (2) and (3) suggests a simple iterative method that can be used to solve the problem.

If we supply an initial “guess” $`r_0\left(\tau \right)`$ for $`r\left(\tau \right)`$ and insert it in equation (1), we find a first approximate solution $`s_0\left(t\right)`$ to $`s\left(t\right)`$. We then in turn insert $`s_0\left(t\right)`$ along with $`r_0\left(\tau \right)`$ in equation (3) to get an approximation $`g_0\left(\tau t\right)`$ for $`g\left(\tau t\right)`$. Finally, $`g_0\left(\tau t\right)`$ is used in equation (2) to get a new function $`r_1\left(\tau \right)`$, and so on. This process can be repeated as often as desired or until convergence is attained.

The final set of equations that define this iterative algorithm can then be written as follows

$`s_i\left(t\right)`$ $`=`$ $`{\displaystyle r_i\left(\tau \right)h\left(t\tau \right)𝑑\tau }`$ (4)
$`r_{i+1}\left(\tau \right)`$ $`=`$ $`r_i\left(\tau \right){\displaystyle \frac{s\left(t\right)}{s_i\left(t\right)}h\left(t\tau \right)𝑑t}.`$ (5)
for $`i=0,1,2,\mathrm{}`$

Finally, two comments to end this section:

* the applicability of the solution to the problem given by equations (4) and (5) is based on the implied assumption that the transfer function of the system $`h\left(t\right)`$ can be measured independently or is known *a priori*. This is true for the problem of the chopping secondary that will be addressed starting in the next section
* it will be noted that the integral in equation (5) is actually a correlation. It follows that the algorithm dictated by the final set of equations can easily be programmed (i.e., computer coded) using subroutines based on the so-called Fast-Fourier-Transform (FFT) methods for convolution and correlation integrals. This is what we have done in the implementation of our technique where we have used Fortran routines presented by Press, Teukolsky, Vetterling & Flannery (1996).

## III The new CSO Chopping Secondary Mirror Assembly

In a simple implementation of equations (4) and (5) one needs a way to generate an input signal to be applied to a given system, measure the output of the system when subjected to this input and finally evaluate the transfer function of the system. In order to accomplish this with our chopping secondary system we modified our assembly from that of Figure 1 to that of Figure 2. We have replaced the square wave generator of our original system by a Real Time Linux (RT Linux) computer which is equipped with the necessary input/output devices (i.e., Analog-to-Digital and Digital-to-Analog converters) to achieve these tasks.

The RT Linux computer also serves as host to the program that performs the necessary calculations and measurements that will allow for the determination of the optimum input to the chopper assembly.

Ideally the sequence of operations would go like this

1. Calibration of the system: signals of constant level are sent to the input of the assembly and the corresponding output levels are measured. In this manner, the “gain” and “offset” of the system are determined and applied to all subsequent input/output operations.
2. Evaluation of the transfer function: this is done by sending a step input signal of a given amplitude to the chopper assembly and calculating the normalized time derivative of the corresponding output signal. This is a very simple way to evaluate a transfer function since the convolution of an arbitrary function with a unit step function is equivalent to the primitive of the original function. This is the technique we use although it should be noted that we also smooth the resulting time derivative with a Savitzky-Golay filter (Press, Teukolsky, Vetterling & Flannery, 1996) to reduce the impact of noise in the application of the Lucy deconvolution technique. We will show some examples of measured transfer functions in the next section.
3. Determination of the desired or targeted output signal $`s\left(t\right)`$.
4. Determination of the optimum input signal: to do so one would *i)* choose an arbitrary waveform as a hypothetical input of the assembly ($`r_0\left(\tau \right)`$ in equations (4) and (5)), *ii)* calculate the corresponding output response $`s_0\left(t\right)`$ of the system using equation (4) and *iii)* calculate a new input $`r_1\left(\tau \right)`$ using equation (5). Repeat *ii)* and *iii)* (using $`r_i\left(\tau \right)`$ and $`s_i\left(t\right)`$, with $`i=1,2,\mathrm{}`$, instead of $`r_0\left(\tau \right)`$ and $`s_0\left(t\right)`$) until convergence to the best input $`r_f\left(\tau \right)`$ signal is attained.
5. Finally, $`r_f\left(\tau \right)`$ is applied to the input of the assembly to produce the output $`s_f\left(t\right)`$ that most resemble the desired output $`s\left(t\right)`$.

We have tested this technique on simple linear systems (e.g., electrical RC filters) with very good results. However, when applied to our chopping secondary mirror assembly the technique did not work in general. It was determined that the non-linearities in the system’s response were the cause of this failure and forced us to bring some changes to the algorithm discussed here.

## IV Implementation of the Method

### IV.1 Non-linearities

Since we have a DC motor as one of the main components of the chopper assembly, it is not surprising that the system should include some non-linearities in its response. As one should expect, the magnetic core of the motor is inherently non-linear as it will experience different amounts of saturation depending on the amplitude of the excitation it is subjected to. That is to say that the transfer function of the system changes with the input signal and that the system reacts differently to different chop throw amplitudes. Moreover, it is also the case that the sign of chop throw will affect the shape of the transfer function. Simply stated, the system has hysteresis and, therefore, does not go “up” the same way it goes “down”.

This will be made clear with the results presented in Figure 3. In this figure we can see the effect that the non-linearities have on the transfer function of the system. The transfer functions shown were measured using the method discussed earlier using rising (“up”) and falling (“down”) step functions (90 arcseconds in amplitude) at two different rest positions (0 and 180 arcseconds for the top and bottom graphs, respectively).

From this it is clear why the algorithm defined by equations (4) and (5) would be doomed to failure. The same two equations can, however, be easily adapted to the problem at hand and make it possible to use the Lucy deconvolution method (albeit somewhat modified) for this kind of non-linear system.

### IV.2 Modifications to the Lucy Deconvolution Method

As was mentioned in the last section, the fact there does not exist a single transfer function that defines the system does not imply that we cannot use the Lucy deconvolution technique to achieve our goal, but we must acknowledge the existence of a family of transfer functions that are dependent on the input signal to the system. That is to say, we should replace $`h\left(t\right)`$ by $`h(t;r)`$, the aforementioned dependence on the input signal $`r\left(t\right)`$ now being made explicit. In practice this means that we now have to measure the transfer functions of the system along a sufficiently refined two-dimensional grid of different step amplitudes (positive and negative) and rest positions. Four examples of such measurements were shown in Figure 3. For the results that will be presented later in this section, we have used a grid where the step amplitude ranges from $`240`$ to 240 arcseconds with a resolution of 30 arcseconds and the rest position spans a similar domain with half the resolution (i.e., 60 arcseconds). It should be noted that this measurement of the transfer functions requires a fair amount of time (as much as 15 to 20 minutes for the grid defined above). But, on the other hand, it needs to be done only once and does not have to be repeated for different chop throw amplitudes and frequencies.

Another important thing to realize is that, contrary to instances where one uses the Lucy technique to deconvolve astronomical images (Houde & Racine, 1994), we are here free to use the system to measure its response to any given input signal and not forced to calculate it through equation (4). This means that the original set of equations (4) and (5) can be reduced to only one equation, namely

$$r_{i+1}\left(\tau \right)=r_i\left(\tau \right)\frac{s\left(t\right)}{s_i\left(t\right)}h(t\tau ;r_i)𝑑t.$$
(6)
With these modifications, the sequence of operations defined in section III now becomes

1. Calibration of the system: signals of constant level are sent to the input of the assembly and the corresponding output levels are measured. In this manner, the “gain” and “offset” of the system are determined and applied to all subsequent input/output operations.
2. Evaluation of the transfer functions: a set of step input signals of differing amplitudes and rest positions are sequentially sent to the chopper assembly and the transfer functions are measured by calculating the normalized time derivative of the corresponding output signals. A Savitzky-Golay filter (Press, Teukolsky, Vetterling & Flannery, 1996) is applied to the functions to reduce the impact of noise on the deconvolution.
3. Determination of the desired or targeted output signal $`s\left(t\right)`$.
4. Determination of the optimum input signal: to do so one would *i)* send an arbitrary wave form $`r_0\left(\tau \right)`$ to the input to the assembly, *ii)* measure the corresponding output response $`s_0\left(t\right)`$ of the system and *iii)* use equation (6) to determine a new input signal $`r_1\left(\tau \right)`$. Repeat *i)*, *ii)* and *iii)* (using $`r_i\left(\tau \right)`$ and $`s_i\left(t\right)`$, with $`i=1,2,\mathrm{}`$, instead of $`r_0\left(\tau \right)`$ and $`s_0\left(t\right)`$) until convergence to the best input $`r_f\left(\tau \right)`$ signal is attained.
5. Finally, $`r_f\left(\tau \right)`$ is applied to the input of the assembly to produce the output $`s_f\left(t\right)`$ that most resemble the desired output $`s\left(t\right)`$.

We have applied this technique to our chopping secondary mirror assembly at the CSO with success. We show typical results in Figures 4 and 5 for chop throws of 60 arcseconds and 300 arcseconds, respectively, at a frequency of 1 Hertz. For this, we chose the initial input signal $`r_0\left(\tau \right)`$ to be a square wave with corresponding amplitudes and frequency, the system’s response to this input is labeled “uncorrected output” in the legends. The desired or “targeted output” signals corresponding to $`s\left(t\right)`$ (also shown on the graphs) in equation (6) in both cases rise (or fall) at the same rate of 3 arcseconds per millisecond when not constant. A comparison of the “uncorrected output” ($`s_0\left(t\right)`$) with the “corrected output” ($`s_f\left(t\right)`$) shows the power of this deconvolution method when applied to this type of problems. In both cases the improvement is significant. Furthermore, it would have been next to impossible to guess which form should the final input signal $`r_f\left(\tau \right)`$ (shown by the “applied input” curves in the graphs) take to obtained the desired outcome. The residual error signal is also plotted in the bottom graph of each figure. As can be seen, the RMS error calculated on flatter parts of the curves are in both cases small ($`1.1`$ arcsecond).

Referring to Figure 3 we see that the transfer function of the system settles down in about 0.5 seconds, which is exactly equal to half of the period of the signals displayed in Figures 4 and 5. This means that the assembly would have just enough time to settle into steady state during half of a cycle when subjected to a square wave of a frequency of 1 Hertz. It would be interesting to see how our technique fares when the period of the input signal is reduced to a value that is significantly less than that the system’s settling time. To test this we have subjected the chopper assembly to a signal of a frequency of 4 Hz and tried to obtained an output of 60 arcseconds in amplitude. This is shown on Figure 6 where now the “targeted output” rises and falls at a rate of 6 arcseconds per millisecond when not constant. Although as could be expected the overall shape of the resulting output signal is somewhat more “rounded” when compared to the results shown in Figure 4, the improvement obtained in going from the “uncorrected output” to the final output signal (i.e., “corrected output” on the graph) is rather significant. In fact, we can see from the bottom graph that for about 52% of a period the response is at most within a few arcseconds from the desired position; the RMS value of the error on that portion of the signal is 1.4 arcseconds.

###### Acknowledgements.

The Caltech Submillimeter Observatory is funded by the NSF through contract AST 9980846. The chopping secondary mirror assembly was designed and built by R. H. Hildebrand’s group at the University of Chicago and was supported by NSF Grant # AST 8917950.

## Figure Captions

Figure 1: The existing chopping secondary mirror assembly at the CSO. A square wave signal is sent to the input of the PID controller and compared with the mirror output position signal (from a LVDT). The resulting processed error signal is sent to a power amplifier which feeds the positioning motor.

Figure 2: The new chopping secondary controller. The new mirror assembly is the same as that of Figure 1, but the input signal is now fed to the PID controller from a RT Linux computer, which hosts the deconvolution program that determines the needed input signal.

Figure 3: Effects of the non-linearities as seen through transfer functions obtained with rising (“up”) and falling (“down”) step functions (90 arcseconds in amplitude) at two different rest positions (0 and 180 arcseconds for the top and bottom graphs, respectively).

Figure 4: Results obtained with our deconvolution technique for a throw of 60 arcseconds at a frequency of 1 Hz (top). The residual error signal is plotted in the bottom graph and its RMS value (0.6 arcsecond) was calculated using data points located between the vertical lines (on the flatter parts of the curve which represent about 84% of a period).

Figure 5: Results obtained with our deconvolution technique for a throw of 300 arcseconds at a frequency of 1 Hz (top). The residual error signal is plotted in the bottom graph and its RMS value (1.1 arcsecond) was calculated using data points located between the vertical lines (on the flatter parts of the curve which represent about 60% of a period).

Figure 6: Results obtained with our deconvolution technique for a throw of 60 arcseconds at a frequency of 4 Hz (top). The residual error signal is plotted in the bottom graph and its RMS value (1.4 arcseconds) was calculated using data points located between the vertical lines (on the flatter parts of the curve which represent about 52% of a period).


format: text

text:
# Probing Decoherence with Electromagnetically Induced Transparency in Superconductive Quantum Circuits

## Abstract

Superconductive quantum circuits (SQCs) comprise quantized energy levels that may be coupled via microwave electromagnetic fields. Described in this way, one may draw a close analogy to atoms with internal (electronic) levels coupled by laser light fields. In this Letter, we present a superconductive analog to electromagnetically induced transparency (S-EIT) that utilizes SQC designs of present day experimental consideration. We discuss how S-EIT can be used to establish macroscopic coherence in such systems and, thereby, utilized as a sensitive probe of decoherence.

Superconductive quantum circuits (SQCs) comprising mesoscopic Josephson junctions can exhibit quantum coherence amongst their macroscopically large degrees of freedom Leggett85a . They exhibit quantized flux and/or charge states depending on their fabrication parameters, and the resultant quantized energy levels are analogous to the quantized internal levels of an atom. Spectroscopy, Rabi oscillation, and Ramsey interferometry experiments have demonstrated that SQCs behave as “artificial atoms” under carefully controlled conditions Nakamura99a ; Friedman00a ; Wal00a ; Vion02a ; Yu02a ; Martinis02a ; Chiorescu03a ; Pashkin03a . This Letter extends the SQC-atom analogy to another quantum optical effect associated with atoms: electromagnetically induced transparency (EIT) Boller91a ; Harris97a . We propose the demonstration of microwave transparency using a superconductive analog to EIT (denoted S-EIT) in a superconductive circuit exhibiting two meta-stable states (e.g., a qubit) and a third, shorter-lived state (e.g., the readout state). We show that driving coherent microwave transitions between the qubit states and the readout state is a demonstration of S-EIT. We further propose a means to use S-EIT to experimentally probe the qubit decoherence rate in a sensitive manner. The philosophy is similar to that in Ref. Ruostekoski99a, , where it was proposed to use EIT to measure phase diffusion in atomic Bose-Einstein condensates.

The three-level $`\mathrm{\Lambda }`$ system illustrated in Fig. 1a is a standard energy level structure utilized in EIT Boller91a ; Harris97a . It comprises two meta-stable states $`|1`$ and $`|2`$, each of which may be coupled to a third excited state $`|3`$. In atoms, the meta-stable states are typically hyperfine or Zeeman levels, while state $`|3`$ is an excited electronic state that may spontaneously decay at a relatively fast rate $`\mathrm{\Gamma }_3`$. In an atomic EIT scheme, a resonant “probe” laser couples the $`|1|3`$ transition, and a “control” laser couples the $`|2|3`$ transition. The transition coupling strengths are characterized by their Rabi frequencies $`\mathrm{\Omega }_{j3}𝐝_{j3}𝐄_{j3}`$ for $`j=1,2`$ respectively, where $`𝐝_{j3}`$ are the dipole matrix elements and $`𝐄_{j3}`$ are the slowly varying envelopes of the electric fields. For particular Rabi frequencies $`\mathrm{\Omega }_{j3}`$, the probe and control fields are effectively decoupled from the atoms by a destructive quantum interference between the states of the two driven transitions. The result is probe and control field transparency Boller91a ; Harris97a . In more recent experiments, ultra-slow light propagation due to EIT-based refractive index modifications in atomic clouds have also been demonstrated Hau99a ; Kash99a ; Budker99a .

SQCs have also been demonstrated to exhibit $`\mathrm{\Lambda }`$-like energy level structures Martinis02a ; Segall03a ; Crankshaw03a ; Yu03a . One example is the persistent-current (PC) qubit, a superconductive loop interrupted by two Josephson junctions of equal size and a third junction scaled smaller in area by the factor $`\alpha <1`$ (Fig. 1b) Orlando99a ; Mooij99a . Its dynamics are described by the Hamiltonian

$$\begin{array}{c}_{pc}=\frac{1}{2}C\left(\frac{\mathrm{\Phi }_0}{2\pi }\right)^2(\dot{\phi }_p^2+(1+2\alpha )\dot{\phi }_m^2)\hfill \\ \hfill +E_j\left[2+\alpha 2\mathrm{cos}\phi _p\mathrm{cos}\phi _m\alpha \mathrm{cos}(2\pi f+2\phi _m)\right],\end{array}$$
(1)
in which $`C`$ is the capacitance of the larger junctions, $`\phi _{p,m}(\phi _1\pm \phi _2)/2`$, $`\phi _i`$ is the gauge-invariant phase across the larger junctions $`i=\{1,2\}`$, $`E_J`$ is the Josephson coupling energy, and $`f`$ is the magnetic flux through the loop in units of the flux quantum $`\mathrm{\Phi }_0`$. The second term in Eq. (1) defines the magnetic-flux-dependent qubit potential landscape. For flux biases near one-half of a flux quantum, $`f1/2`$, the potential may be approximated by a double-well potential, with each well corresponding to a distinct, stable, classical state of the electric current, i.e., left or right circulation through the loop. In turn, each current state has a net magnetization of opposite direction that is measurable using a dc SQUID Orlando99a . As a quantum object, the potential wells exhibit quantized energy levels corresponding to the quantum states of the macroscopic circulating current Segall03a ; Crankshaw03a . These levels may be coupled using microwave radiation Wal00a ; Yu03a , and their quantum coherence has been experimentally demonstrated Chiorescu03a . Note that for this system, the terms “population” and “occupation probability” are used synonymously.

Tuning the flux bias away from $`f=1/2`$ results in the asymmetric double-well potential illustrated in Fig. 1c. The three states in the left well constitute the superconductive analog to the atomic $`\mathrm{\Lambda }`$ system. States $`|1`$ and $`|2`$ are “meta-stable states,” with a tunneling and coherence time much longer than the excited “readout” state $`|3`$. State $`|3`$ has weakly-coupled intra-well transitions, but has a strong inter-well transition when tuned on resonance with state $`|4`$ Segall03a ; Crankshaw03a ; Yu03a . Using tight-binding models with experimental PC qubit parameters Orlando99a ; Yu03a ; Crankshaw03a  at a flux bias f=0.5041, we estimate the tunneling times from states $`|1`$, $`|2`$, and $`|3`$ to the right well are $`1/\mathrm{\Gamma }_1`$ 1 ms, $`1/\mathrm{\Gamma }_2`$$`\mu `$s, and $`1/\mathrm{\Gamma }_3`$ 1 ns respectively. Thus, a particle reaching state $`|3`$ will tend to tunnel quickly to state $`|4`$, and this event results in a switching of the circulating current that may be detected using a fast-measurement scheme. Alternatively, one may detune states $`|3`$ and $`|4`$, and then apply a resonant $`\pi `$-pulse to transfer the population from state $`|3`$ to $`|4`$; since the states are now off-resonance, the relaxation rate back to the left well is reduced and a slower detection scheme may be used. We note that a single-junction qubit Martinis02a  shares this property, since the right well in Fig. 1c is effectively replaced by a quasi-continuum of states, and transitions out of the left well will not return.

Transitions between the quantized levels are driven by resonant microwave-frequency magnetic fields. Assuming the Rabi frequencies $`\mathrm{\Omega }_{ij}`$ to be much smaller than all level spacings $`|\omega _{kl}||\omega _k\omega _l|`$, the system-field interaction may be written within the rotating wave approximation (RWA)Scully97a ,

$`_{int}^{(\mathrm{RWA})}={\displaystyle \frac{\mathrm{}}{2}}\left[\begin{array}{ccc}0& \mathrm{\Omega }_{12}^{}& \mathrm{\Omega }_{13}^{}\\ \mathrm{\Omega }_{12}& 0& \mathrm{\Omega }_{23}^{}\\ \mathrm{\Omega }_{13}& \mathrm{\Omega }_{23}& i\mathrm{\Gamma }_3\end{array}\right],`$ (2)
in which the decay from state $`|3`$ is treated phenomenologically as a non-Hermitian matrix element Scully97a ; Dutton02a . For small microwave perturbations of amplitude $`f_\mathrm{\Delta }`$, the associated Rabi frequencies are given by $`\mathrm{\Omega }_{pq}=f_\mathrm{\Delta }p|\mathrm{sin}(2\pi f+2\varphi _m)|q`$; numerical simulations of the matrix elements using PC qubit parameters are consistent with recent experimental results (see caption Fig. 1c) Segall03a ; Crankshaw03a ; Yu03a . In general, all three intra-well transitions are allowed in SQCs. For example, consider states $`|1`$ and $`|2`$ to be a qubit that is prepared in an arbitrary superposition state $`|\mathrm{\Psi }=c_1|1+c_2|2`$ by temporarily driving the $`\mathrm{\Omega }_{12}`$ transition. Then, by applying a probe field, the population of state $`|1`$ may be read out through a transition to state $`|3`$ followed by a rapid escape to the right well Martinis02a . In this case, both the preparation and readout transitions were allowed and absorptive.

One may achieve S-EIT in a superconductive $`\mathrm{\Lambda }`$ system that is prepared in state $`|\mathrm{\Psi }=c_1|1+c_2|2`$, by simultaneously and solely applying the microwave fields $`\mathrm{\Omega }_{13}`$ and $`\mathrm{\Omega }_{23}`$ such that

$`{\displaystyle \frac{\mathrm{\Omega }_{13}}{\mathrm{\Omega }_{23}}}={\displaystyle \frac{c_2}{c_1}}.`$ (3)
Under this condition (with $`\mathrm{\Omega }_{12}=0`$), the state $`|\mathrm{\Psi }`$ is an eigenstate of $`_{int}^{(\mathrm{RWA})}`$ in Eq. (2) with eigenvalue zero, and the SQC becomes transparent to the microwave fields. As in conventional EIT, the amplitudes for the two absorption transitions into $`|3`$ have equal and opposite probability amplitudes, leading to a destructive quantum interference. Thus, in the absence of decoherence, preparing the qubit in state $`|\mathrm{\Psi }`$ with an ideal preparation ($`\mathrm{\Omega }_{12}`$) field and subsequently applying ideal probe ($`\mathrm{\Omega }_{13}`$) and control ($`\mathrm{\Omega }_{23}`$) microwave fields which satisfy Eq. (3) would result in no population loss through the readout state $`|3`$. In this way, S-EIT would confirm, without disturbing the system, that we had indeed prepared the qubit in the desired state.

In a practical SQC, there will be decoherence of the state $`|\mathrm{\Psi }`$, and this must be measured, characterized, and minimized for quantum information applications. S-EIT is one such sensitive decoherence probe, since deviations in the amplitude and/or relative phase of the complex coefficients $`c_i`$ from the condition established in Eq. (3) result in a small probability $`|(c_1\mathrm{\Omega }_{13}+c_2\mathrm{\Omega }_{23})/\mathrm{\Omega }|^2`$ of the SQC being driven into the readout state $`|3`$ on a time scale $`\mathrm{\Gamma }_3/\mathrm{\Omega }^2`$. In general, there are two categories of decoherence: loss and dephasing. Loss refers to population losses from the metastable states $`|1,|2`$, and it is present in an SQC due to, for example, the finite loss rate of level $`|2`$, $`\mathrm{\Gamma }_21/\mu \mathrm{s}=(2\pi )0.2\mathrm{MHz}`$. Dephasing refers to interactions of the SQC with other degrees of the freedom in the system that cause the relative phase between $`c_1`$ and $`c_2`$ to diffuse. The incorporation of dephasing is facilitated by the use of a density matrix formalism.

We describe the system with a $`3\times 3`$ density matrix with diagonal elements $`\rho _{ii}`$ describing the populations, and $`\rho _{ij},ij`$ describing the coherences between levels. In the presence of the EIT fields $`\mathrm{\Omega }_{13}`$ and $`\mathrm{\Omega }_{23}`$ with no direct coupling ($`\mathrm{\Omega }_{12}=0`$), the Bloch equations govern the evolution of the density matrix Scully97a :

$`\dot{\rho }_{11}`$ $`=\mathrm{\Gamma }_1\rho _{11}{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{13}^{}\rho _{31}+{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{13}\rho _{13},`$ (4)
$`\dot{\rho }_{22}`$ $`=\mathrm{\Gamma }_2\rho _{22}{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{23}^{}\rho _{32}+{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{23}\rho _{23},`$ (5)
$`\dot{\rho }_{33}`$ $`=\mathrm{\Gamma }_3\rho _{33}+{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{13}^{}\rho _{31}{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{13}\rho _{13}`$
$`+{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{23}^{}\rho _{32}{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{23}\rho _{23},`$ (6)
$`\dot{\rho }_{12}`$ $`=\gamma _{12}\rho _{12}{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{13}^{}\rho _{32}+{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{23}\rho _{13},`$ (7)
$`\dot{\rho }_{13}`$ $`=\gamma _{13}\rho _{13}+{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{13}^{}(\rho _{11}\rho _{33})+{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{23}^{}\rho _{12},`$ (8)
$`\dot{\rho }_{23}`$ $`=\gamma _{23}\rho _{23}+{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{23}^{}(\rho _{22}\rho _{33})+{\displaystyle \frac{i}{2}}\mathrm{\Omega }_{13}^{}\rho _{21}.`$ (9)
The remaining three elements’ equations are determined by $`\rho _{ij}^{}=\rho _{ji}`$. The decoherence rates $`\gamma _{ij}=(\mathrm{\Gamma }_i+\mathrm{\Gamma }_j)/2+\gamma _{ij}^{(\mathrm{deph})}`$ include both loss and dephasing contributions. We concentrate on the regime in which the readout state escape rate $`\mathrm{\Gamma }_3=1\text{ ns}^1=(2\pi )130\mathrm{MHz}`$ dominates all other loss and dephasing rates, thus $`\gamma _{13}\gamma _{23}\mathrm{\Gamma }_3/2`$. Furthermore, we ignore the meta-stable state losses $`\mathrm{\Gamma }_1,\mathrm{\Gamma }_2`$ relative to the dephasing $`\gamma _{12}^{(\mathrm{deph})}`$ and set $`\gamma _{12}\gamma _{12}^{(\mathrm{deph})}`$. Theoretical estimates of dephasing rates, such as $`\gamma _{12}^{(\mathrm{deph})}`$, in multi-level systems were recently obtained in Ref. Burkhard03a, .

We illustrate an S-EIT decoherence probe example by applying EIT fields $`\mathrm{\Omega }_{13}=\mathrm{\Omega }_{23}=(2\pi )150\mathrm{MHz}`$ to the dark state $`|\mathrm{\Psi }=(|1|2)/\sqrt{2}`$ with a dephasing rate $`\gamma _{12}=(2\pi )5\mathrm{MHz}`$ and numerically integrating Eqs. (4)-(9). The dephasing $`\gamma _{12}`$ causes a small population transfer to the excited state $`\rho _{33}`$ (Fig. 2a). The excited state population initially exhibits a rapid rise (see inset Fig. 2a) with transitory oscillations, reaching its maximum value $`\rho _{33}^{(\mathrm{max})}`$ within about 4 ns. This is followed by a smooth decay with a $`1/e`$ time of about 80 ns. The solid curve in Fig. 2b traces the total population $`P=\rho _{11}+\rho _{22}+\rho _{33}`$ remaining in the system as a function of time. When the excited state maximum $`\rho _{33}^{(\mathrm{max})}`$ is reached, the total remaining population is $`P(4\text{ ns})=0.973`$. In contrast, the dashed line in Fig. 2b illustrates the rapid population loss expected when the same fields are applied to the state $`|\mathrm{\Psi }=(|1+|2)/\sqrt{2}`$ \[$`\pi `$ out of phase with the dark state in Eq. (3)\]. In the absence of S-EIT quantum interference, the entire population is lost on a time scale $`\mathrm{\Gamma }_3/\mathrm{\Omega }^21\mathrm{ns}`$. The general behavior presented in Figs. 2a and 2b is observed over a wide parameter regime of experimental interest.

We now use Eqs. (4)-(9) to show how measuring the slow population loss in S-EIT can be used to extract the decoherence rate $`\gamma _{12}`$. The elements $`\rho _{33}`$, $`\rho _{13}`$, and $`\rho _{23}`$ in Eqs. Probing Decoherence with Electromagnetically Induced Transparency in Superconductive Quantum Circuits8, and 9 are damped at a rapid rate $`\mathrm{\Gamma }_3`$, allowing their adiabatic elimination Javanainen95a ; Dutton02a ; we solve for their quasi-steady state values by setting $`\dot{\rho }_{33}=\dot{\rho }_{13}=\dot{\rho }_{23}=0`$. This approximation is accurate once initial transients have passed and the plateau value $`\rho _{33}^{(\mathrm{max})}`$ has been reached. Using these results in Eq. (7) yields an equation for $`\dot{\rho }_{12}`$ with a strong damping term $`\mathrm{\Omega }^2/\mathrm{\Gamma }_3`$, and it too can be solved for its quasi-steady state value. In the limit $`\gamma _{12}\mathrm{\Gamma }_3/\mathrm{\Omega }^21`$ we get Murali03a

$`\rho _{12}(t){\displaystyle \frac{\mathrm{\Omega }_{13}\mathrm{\Omega }_{23}}{\mathrm{\Omega }^2}}\left(1{\displaystyle \frac{2\gamma _{12}\mathrm{\Gamma }_3}{\mathrm{\Omega }^2}}\right)\left(\rho _{11}(t)+\rho _{22}(t)\right).`$ (10)
The ratio $`2\gamma _{12}\mathrm{\Gamma }_3/\mathrm{\Omega }^2`$ represents the small fractional deviation of $`\rho _{12}`$ from its dark state value. There is a competition between the “preparation rate” $`\mathrm{\Omega }^2/\mathrm{\Gamma }_3`$ (which constantly acts to drive the system into the dark state) and the decoherence rate $`\gamma _{12}`$ (which attempts to drive it back out).

We now consider the conditions under which one can use S-EIT to estimate the decoherence rate with little population loss in the system. Eqs. (4), (5), and (10) reveal that deviations from the dark state cause population loss through $`|3`$ at a rate $`R=2\gamma _{12}(\mathrm{\Omega }_{13}^2\mathrm{\Omega }_{23}^2/\mathrm{\Omega }^4)`$; this ultimately leads to the exponential decay of $`P`$ seen in Fig. 2b. By assumption, population escapes the system only through the decay term $`\rho _{33}\mathrm{\Gamma }_3`$ in Eq. (Probing Decoherence with Electromagnetically Induced Transparency in Superconductive Quantum Circuits), yielding $`\rho _{33}=(R/\mathrm{\Gamma }_3)P`$. At the early time when $`\rho _{33}^{(\mathrm{max})}`$ has been reached, the population $`P`$ remains close to unity, and the excited state population reaches

$`\rho _{33}^{(\mathrm{max})}2{\displaystyle \frac{\mathrm{\Omega }_{13}^2\mathrm{\Omega }_{23}^2}{\mathrm{\Omega }^4}}{\displaystyle \frac{\gamma _{12}}{\mathrm{\Gamma }_3}}.`$ (11)
The time $`T_{ss}`$ to reach $`\rho _{33}^{(\mathrm{max})}`$ is generally the smaller of the preparation time $`\mathrm{\Gamma }_3/\mathrm{\Omega }^2`$ and the inverse of the decay rate $`1/\mathrm{\Gamma }_3`$. At this time, the total population loss will be $`T_{ss}\rho _{33}^{(\mathrm{max})}(2\mathrm{\Omega }_{13}^2\mathrm{\Omega }_{23}^2/\mathrm{\Omega }^4)\mathrm{Max}(\gamma _{12}/\mathrm{\Gamma }_3,\gamma _{12}\mathrm{\Gamma }_3/\mathrm{\Omega }^2)`$. So long as the loss during this initial transient time is small, the population will follow a simple exponential decay $`P(t)=\mathrm{exp}(\rho _{33}^{\text{(max)}}\mathrm{\Gamma }_3t)`$, and the dephasing rate $`\gamma _{12}`$ can be easily extracted. To keep this loss small, we require both ratios in the $`\mathrm{Max}(\mathrm{})`$ argument to be small ($`\mathrm{\Omega }\sqrt{2\gamma _{12}\mathrm{\Gamma }_3}`$ and $`\mathrm{\Gamma }_3\gamma _{12}`$) in order to use this approach to accurately estimate the decoherence rate while causing little loss from the system. Since $`\mathrm{\Omega }`$ is experimentally controllable, it can be chosen to satisfy the first constraint. If $`\mathrm{\Gamma }_3`$ is comparable or smaller than $`\gamma _{12}`$, then S-EIT remains a decoherence probe, although the strong damping assumption leading to Eq. (10) no longer holds and so the analysis is different.

We have performed a series of numerical simulations, varying $`\gamma _{12}`$ to test the validity of the above approach. The results are presented in Fig. 3. Fig. 3a indicates $`\rho _{33}^{(\mathrm{max})}`$ versus $`\gamma _{12}`$ and compares the results with the analytic estimate \[Eq. (11)\]. The agreement is quite good for $`\gamma _{12}<(2\pi )4\mathrm{MHz}`$, which corresponds to the $`2\gamma _{12}\mathrm{\Gamma }_3/\mathrm{\Omega }^2<0.056`$. Higher dephasing rates compete more with the preparation rate, making the adiabatic elimination approach less valid; this leads to deviations from our analytic prediction \[Eq. (11)\]. In such cases, one observes a significant loss by the time $`\rho _{33}^{(\mathrm{max})}`$ is reached, as illustrated in Fig. 3b

The RWA invoked in the calculation ignores far off-resonance couplings induced by other applied fields at other transitions (e.g., in our case, the $`\mathrm{\Omega }_{23}`$ field drives $`|1|2`$ at $`\mathrm{\Delta }4\mathrm{GHz}`$ off resonance). We have performed calculations including all such off-resonant couplings. Generally, they lead to small shifts of the energy levels (analogous to A.C. Stark shifts) and loss rates. These losses scale as $`\mathrm{\Gamma }_3(\mathrm{\Omega }^2/\mathrm{\Delta }^2)`$ and $`\mathrm{\Omega }^4/\mathrm{\Gamma }_3^2\mathrm{\Delta }`$, putting a limit on the field strengths $`\mathrm{\Omega }^2`$ which can be used for a given $`\mathrm{\Delta }`$. For the parameters considered here, we found shifts of $`(2\pi )8\mathrm{MHz}`$ and loss rates totalling $`(2\pi )19\mathrm{kHz}`$; this should not effect measurements in the regime $`\gamma _{12}(2\pi )19\mathrm{kHz}`$ Murali03a .

We have proposed using the superconductive analog to EIT (S-EIT) to demonstrate macroscopic quantum interference in superconductive quantum circuits. S-EIT provides an accurate and sensitive means to probe the accuracy and phase coherence of qubit preparation, and we have calculated analytic expressions for the field strengths required for this purpose.

This work was supported in part by the AFOSR grant No. F49620-01-1-0457 under the Department of Defense University Research Initiative in Nanotechnology (DURINT). The work at Lincoln Laboratory was sponsored by the AFOSR under Air Force Contract No. F19628-00-C-0002.


text:
# Evidence for grain growth in T Tauri disksBased on observations made with the ESO 3.6m Telescope at the La Silla Observatory under program ID 68.D-0537(A), 69.C-0268 and 70.C-0544.

## 1 Introduction

T Tauri stars are low mass young stellar objects with ages of $`10^610^7`$ yr. One of the characteristics of so-called classical T Tauri stars is a distinctive IR excess. The origin of the IR radiation is the accretion disk which plays a crucial role in star formation by transporting material to the central star while angular momentum is transferred outwards. Different models were proposed to bring the observed facts into a general framework (Lynden-Bell & Pringle (1974); Adams, Lada, & Shu (1987); Chiang & Goldreich (1997)). Yet there remain many open questions concerning the details of configuration, chemistry and evolution of disks. A particular feature of the infrared spectra of such disks is the silicate band in the region from 8 to 12 $`\mu `$m. The band originates from the stretching mode of the Si-O bond of silicate minerals like olivine (\[Mg,Fe\]2SiO4), forsterite (Mg2SiO4), enstatite (MgSiO3) and silica (SiO2). The silicate feature observed in T Tauri stars appears in emission as well as in absorption (Cohen & Witteborn (1985)). Silicate emission is assumed to emerge from a warm, optically thin disk layer (disk atmosphere) which is heated by the radiation of the central star (Chiang & Goldreich (1997), Natta, Meyer, & Beckwith (2000)). Honda et al. (2003) showed that crystalline silicates are present in T Tauri stars, indicating substantial grain processing. In this Paper we present studies on the silicate emission from a sample of 14 T Tauri stars based on mid-infrared spectroscopy performed during three observation campaigns in February, June and December 2002 on the ESO 3.6 m telescope at La Silla. We analyse the spectra in particular with respect to a possible correlation between shape and strength of the silicate feature. Such a correlation has been recently observed for Herbig Ae/Be stars (van Boekel et al. (2003)) and has been interpreted as evidence of grain processing in the circumstellar disks of those stars.

## 2 Observations and data reduction

The observations of the T Tauri stars were performed with TIMMI2, the Thermal Infrared Multi Mode Instrument 2 (Reimann et al. (1998)) installed at the 3.6 m telescope at ESO’s observatory at La Silla, Chile. We observed each object in two modes: first we used the imaging mode to obtain photometry at 11.9 $`\mu `$m (band width 1.2 $`\mu `$m). In a number of cases in addition we used the 8.9 $`\mu `$m filter (band width 0.8 $`\mu `$m). Second we performed spectroscopy using the low-resolution grism mode (R=160) with a slit width of 1.2″. In both modes the system was working with chopping and nodding with an amplitude of 10″. In our data reduction, the measurements of the spectrophotometric standards were used to determine the shape of the spectra, and the photometric data at 11.9 $`\mu `$m were used to establish the absolute flux level. The spectrum of the atmospheric extinction, which is needed to perform the first of these steps, was determined from observations of standard stars at different airmass. We regularly observed standard stars to monitor variations in atmospheric transmission. The flux data for the spectrophotometric standard stars, (based on models by Cohen (1998)), were taken from ESO’s TIMMI2 webpage.

In total we observed 21 T Tauri stars with the silicate feature in emission. We selected the 14 objects for which the silicate feature was measured with a signal-to-noise ratio better than 2 (see Tab. 1). Several of these are known to be binary or multiple systems. In two cases we were able to obtain separate spectra of the components. The remaining 7 sources (DO Tau, DQ Tau, FU Ori, BBW 76, VZ Cha, VW Cha, Sz 82) were degraded by noise to the extent that a meaningful spectral analysis was not possible. For some stars the observations were repeated to verify the results.

For our studies we assumed that the continuum outside of the silicate band is actually reached at the edges of our measured spectral range (at 8.2 $`\mu `$m and at 13.0 $`\mu `$m). This assumption is supported by comparison of our spectra with ISOPHOT spectra (Ábrahám et. al 2003 in prep, see also Natta, Meyer, & Beckwith (2000)). Figure 1 illustrates this for two examples. Here, the continuum was estimated by fitting a second order polynomial to all data points shortwards of 8 $`\mu `$m and beyond 12.8 $`\mu `$m. Note that the fitted curve is matching the ends of the TIMMI2 spectra and has an almost linear characteristic over the range of the silicate feature. Since we were not able to find such additional data of good quality for all objects of our sample, we estimated the continuum level within the silicate feature by a linear interpolation between the end points of the TIMMI2 spectra at 8.2 $`\mu `$m and at 13.0 $`\mu `$m, where the values for these end points were obtained by averaging over the data over an interval of 0.2 $`\mu `$m. The continuum level determined this way agreed with the second order fit, where available, to typically $`\pm `$ 10 % . For all sources we consistently used the derived linearly varying continuum level to determine the continuum normalized111to preserve the shape of the emission feature we define:
$`F_{\text{norm}}(\lambda )=1+[F_{\text{total}}(\lambda )F_{\text{continuum}}(\lambda )]/F_{\text{continuum}}`$ spectra; these are shown in Fig. 2.

## 3 Analysis of the silicate feature

The spectra in Fig. 2 show a variety in strength and shape. In some spectra (AK Sco, Haro 1-16, CR Cha) a strong silicate emission with a peak near 9.8 $`\mu `$m is visible. In other cases (Ru Lup, DR Tau, HBC 639), the silicate emission is weaker and the profile looks more like a plateau. It is known that the shape of the band is strongly affected by the chemical composition and in particular by the size of the dust grains (Henning et al. (1995), Bouwman et al. (2001)). A correlation between strength and shape of the silicate emission could be expected if there are systematic changes in particle size within the group of objects.

To estimate the strength of the silicate emission, we use the peak of the continuum normalised spectrum. For the characterisation of the band shape we measured the continuum subtracted flux at two wavelengths, 11.3 and 9.8 $`\mu `$m. The flux ratio for these wavelengths can be interpreted as a indicator for the grain size composition following Bouwman et al. (2001). They point out that (see Fig. 3)

* the mid-infrared absorption coefficient of small amorphous olivine dust grains (with a size of 0.1 $`\mu `$m) has a triangular shape with the maximum at 9.8 $`\mu `$m. This spectrum is compatible with observations of the silicate emission of unprocessed dust of the ISM (see also Bowey, Adamson, & Whittet (1998)).
* the mid-infrared absorption coefficient of larger amorphous olivine dust grains shows a plateau-like shape in the range from 9.5 to 12 $`\mu `$m.

Assuming a dust temperature of around 300 K, the blackbody thermal emission will not vary very much over our wavelength band, and the silicate emission profiles would have a shape quite close to the shape of the absorption coefficients. Since the absorption cross section per unit mass of large grains is in general smaller than that of small grains of similar shape, the emission of large grains is in general expected to be weaker than that of small grains.

Our observations suggest that most observed mid-infrared spectra of young low or intermediate mass stars can be characterised simply by emission corresponding to the two grain sizes covered in Fig. 4. We note that the ratio of absorption coefficients at the wavelengths of 11.3 $`\mu `$m and 9.8 $`\mu `$m as function of grain size varies from 0.43 (for 0.1 $`\mu `$m grains) to about 1.0 (2.0 $`\mu `$m grains). Therefore, the ratio of the (continuum subtracted) flux at 11.3 $`\mu `$m over that at 9.8 $`\mu `$m should be a measure of the grain size. In Fig. 4 we plotted this ratio against the strength of the silicate feature. The distribution of the data points shows a clear correlation. Strong silicate features have 11.3/9.8 flux ratios down to 0.53 (for CR Cha), while weak features exhibit higher ratios up to about 1.0 (in case of the S CrA NW-component even 1.42). The remaining points are distributed over the range between these two extremes. The error bars quoted in Fig. 4 take into account the signal to noise of the observations, and the uncertainty in the continuum estimate. A weighted linear fit to the data points with the relation $`y=abx`$, where x is the strength of the feature and y the flux ratio, gives the coefficients: $`a=1.11\pm 0.11`$ and $`b=0.16\pm 0.04`$. A very similar correlation between strength and shape of the silicate emission was found by van Boekel et al. (2003) for isolated Herbig Ae/Be stars, with coefficients of $`a=1.48\pm 0.09`$, $`b=0.28\pm 0.04`$.

Emission by Polycyclic Aromatic Hydrocarbons (PAHs) at 11.3 $`\mu `$m does not appear to be a problem for our data set. There was no evidence for PAH emission in our observed spectra nor in the available ISOPHOT-S spectra. This absence of PAH emission is not unexpected, since one can assume that the UV radiation of T Tauri stars normally is too weak to cause a significant excitation of PAH (Natta & Kruegel (1995)).

## 4 Discussion

The large difference in geometric cross section per unit mass between large and small dust grains implies that if the silicate emission band is dominated by small grains, a substantial amount of large grains could be present without being noticed. On the other hand, the spectra with plateau-like emission are in geneal weak, likely due to a much lower number of small grains present in the upper disk regions. In some spectra (S CrA both components, RU Lup) the flux at 11.3 $`\mu `$m is even higher than the flux at 9.8 $`\mu `$m, which results to a ratio greater than 1. This cannot be explained by the assumption of amorphous olivine dust with the absorption spectra shown above. It may be that crystalline forsterite (Mg2SiO4) is responsible for this effect, since its absorption spectrum shows a peak at 11.3 $`\mu `$m (see Bouwman et al. (2001), Honda et al. (2003)). Possible indications for crystalline forsterite are also seen in some other objects with stronger features (SU Aur, GG Tau).

We also note that the silicate emission emerges from warm, optically thin material, assumed here to be located in the upper layer of a flared accretion disk, which is irradiated by the central star. The dust in these upper layers of different T Tauri disks appears to vary continously between two extreme types: a) small amorphous dust grains or a mixture of small and large grains, with a spectrum very similar to the dust of the ISM and b) large amorphous dust grains with a possible addition of crystalline silicate. If we think of these differences in terms of evolution, then the reason for the difference between these two dust types could be either the removal of the small grains from the upper disk layers or the coagulation of small grains to larger ones. A removal of small dust grains by pure settling to the mid-plane of the disk is unlikely to explain the observed kind of differentiation since the settling time for larger grains is shorter than for the smaller ones. A further mechanism which could be responsible for the vanishing of small grains is by the radiation pressure of the central star. A star of solar luminosity would be quite effective in removing tenth micron sized particles. But in this case, a permanent diffusion of particles from the mid-plane on shorter time scales should avoid a significant decrease of the amount of small particles in the upper layer (Takeuchi & Lin (2003)). The effect of coagulation is a possible explanation, if constructive collisions to larger grains dominate over destructive collisions which produce again small grains. This should be the case in T Tauri disks where the SED can be well fit by gas-rich disks in hydrostatic equilibrium.

An interesting detail of our studies is the availability of separate spectra from the components of two binary T Tauri systems (AS 205 and S CrA). We observed significant differences in strength and also differences in the shape of the silicate band of the components, with the secondary having the - relatively - stronger emission. Since the components presumably have the same age, the difference must have other reasons than the evolutionary state of the stars.

###### Acknowledgements.

This work is partly supported by the Hungarian Research Fund (OTKA No. T037508) and by the Bolyai Fellowship (P.Àbrahàm). S.Y. Melnikov thanks the German Academic Exchange Service (DAAD) for the DAAD fellowship grant.


text:
# Neural cryptography with feedback

## I Introduction

Neural networks learn from examples. When a system of interacting neurons adjusts its couplings to a set of externally produced examples, this network is able to estimate the rule which produced the examples. The properties of such networks have successfully been investigated using models and methods of statistical physics Hertz et al. (1991); Engel and Van den Broeck (2001).

Recently this research program has been extended to study the properties of interacting networks Metzler et al. (2000); Kinzel et al. (2000). Two networks which learn the examples produced by their partner are able to synchronize. This means that after a training period the two networks achieve identical time dependent couplings (synaptic weights). Synchronization by mutual learning is a phenomenon which has been applied to cryptography Kanter et al. (2002); Kinzel and Kanter (2003).

To send a secret message over a public channel one needs a secret key, either for encryption, decryption, or both. In 1976, Diffie and Hellmann have shown how to generate a secret key over a public channel without exchanging any secret message before. This method is based on the fact that—up to now—no algorithm is known which finds the discrete logarithm of large numbers by feasible computer power Stinson (1995).

Recently it has been shown how to use synchronization of neural networks to generate secret keys over public channels Kanter et al. (2002). This algorithm, called neural cryptography, is not based on number theory but it contains a physical mechanism: The competition between stochastic attractive and repulsive forces. When this competition is carefully balanced, two partners $`A`$ and $`B`$ are able to synchronize whereas an attacking network $`E`$ has only a very low probability to find the common state of the communicating partners.

The security of neural cryptography is still being debated and investigated Klimov et al. (2003); Kinzel and Kanter (2002a, b); Mislovaty et al. (2002); Kanter and Kinzel (2003). In this paper we introduce a mechanism which is based on the generation of inputs by feedback. This feedback mechanism increases the repulsive forces between the participating networks, and the amount of the feedback, the strength of this force, is controlled by an additional parameter of our model.

A measure of the security of the system is the probability $`P_E`$ that an attacking network is successful. We calculate $`P_E`$ obtained from the best known attack Klimov et al. (2003) for different model parameters and search for scaling properties of the synchronization time as well as for the security measure. It turns out that feedback improves the security significantly, but it also increases the effort to find the common key. When this effort is kept constant, feedback only yields a small improvement of security.

## II Repulsive and attractive stochastic forces

The mathematical model used in this paper is called a tree parity machine (TPM), sketched in Fig. 1. It consists of $`K`$ different hidden units, each of them being a perceptron with an $`N`$-dimensional weight vector $`𝐰_k`$. When a hidden unit $`k`$ receives an $`N`$-dimensional input vector $`𝐱_k`$ it produces the output bit

$$\sigma _k=\text{sgn}(𝐰_k𝐱_k).$$
(1)
The $`K`$ hidden units $`\sigma _k`$ define a common output bit $`\tau `$ of the total network by

$$\tau =\underset{k=1}{\overset{K}{}}\sigma _k.$$
(2)
In this paper we consider binary input values $`x_{k,j}\{1,+1\}`$ and discrete weights $`w_{k,j}\{L,L+1,\mathrm{},L1,L\}`$, where the index $`j`$ denotes the $`N`$ components and $`k`$ the $`K`$ hidden units.

Each of the two communicating partners $`A`$ and $`B`$ has its own network with an identical TPM architecture. Each partner selects random initial weight vectors $`𝐰_k^A(t=0)`$ and $`𝐰_k^B(t=0)`$.

Both of the networks are trained by their mutual output bits $`\tau ^A`$ and $`\tau ^B`$. At each training step, the two networks receive common input vectors $`𝐱_k`$ and the corresponding output bit $`\tau `$ of its partner. We use the following learning rule.

1. If the output bits are different, $`\tau ^A\tau ^B`$, nothing is changed.
2. If $`\tau ^A=\tau ^B\tau `$ only the hidden units are trained which have an output bit identical to the common output, $`\sigma _k^{A/B}=\tau ^{A/B}`$.
3. To adjust the weights we consider three different learning rules.

1. Anti-Hebbian learning

    $$𝐰_k^+=𝐰_k\tau 𝐱_k\mathrm{\Theta }(\sigma _k\tau )\mathrm{\Theta }(\tau ^A\tau ^B).$$
    (3)
    2. Hebbian learning

    $$𝐰_k^+=𝐰_k+\tau 𝐱_k\mathrm{\Theta }(\sigma _k\tau )\mathrm{\Theta }(\tau ^A\tau ^B).$$
    (4)
    3. Random walk

    $$𝐰_k^+=𝐰_k+𝐱_k\mathrm{\Theta }(\sigma _k\tau )\mathrm{\Theta }(\tau ^A\tau ^B).$$
    (5)

If any component $`w_{k,j}`$ moves out of the interval $`L,\mathrm{},L`$, it is replaced by $`\mathrm{sgn}(w_{k,j})L`$.

Note that for the last rule, the dynamics of each component is identical to a random walk with reflecting boundaries. The only difference to usual random walks is that the dynamics is controlled by the $`2K`$ global signals $`\sigma _k^{A/B}`$ which, in turn, are determined by the ensemble of random walks. Two corresponding components of the weights of $`A`$ and $`B`$ receive an identical input $`x_{k,j}`$, hence they move into the same direction if the control signal allows both of them the move. As soon as one of the two corresponding components hits the boundary their mutual distance decreases. This mechanism finally leads to complete synchronization, $`𝐰_k^A(t)=𝐰_k^B(t)`$ for all $`tt_{sync}`$.

On average, a common step leads to an attractive force between the corresponding weight vectors. If, however, only the weight vector of one of the two partners is changed the distance between corresponding vectors increases, on average. This may be considered as a repulsive force between the corresponding hidden units.

A learning step in at least one of the $`K`$ hidden units occurs if the two output bits are identical, $`\tau ^A=\tau ^B`$. In this case, there are three possibilities for a given pair of hidden units:

1. an attractive move for $`\sigma _k^A=\sigma _k^B=\tau ^{A/B}`$;
2. a repulsive move for $`\sigma _k^A\sigma _k^B`$;
3. and no move at all for $`\sigma _k^A=\sigma _k^B\tau ^{A/B}`$.

We want to calculate the probabilities for repulsive and attractive steps Klimov et al. (2003); Rosen-Zvi et al. (2002a). The distance between two hidden units can be defined by their mutual overlap

$$\rho _k=\frac{𝐰_k^A𝐰_k^B}{\sqrt{𝐰_k^A𝐰_k^A}\sqrt{𝐰_k^B𝐰_k^B}}.$$
(6)
The probability $`ϵ_k`$ that a common randomly chosen input $`𝐱_k`$ leads to a different output bit $`\sigma _k^A\sigma _k^B`$ of the hidden unit is given by Engel and Van den Broeck (2001)

$$ϵ_k=\frac{1}{\pi }\mathrm{arccos}\rho _k.$$
(7)
The quantity $`ϵ_k`$ is a measure of the distance between the weight vectors of the corresponding hidden units. Since different hidden units are independent, the values $`ϵ_k`$ determine also the conditional probability $`P_r`$ for a repulsive step between two hidden units given identical output bits of the two TPMs. In the case of identical distances, $`ϵ_k=ϵ`$, one finds for $`K=3`$

$`P_r`$ $`=`$ $`P(\sigma _k^A\sigma _k^B|\tau ^A=\tau ^B)`$ (8)
$`=`$ $`{\displaystyle \frac{2(1ϵ)ϵ^2}{(1ϵ)^3+3(1ϵ)ϵ^2}}.`$
On the other side, an attacker $`E`$ may use the same algorithm as the two partners $`A`$ and $`B`$. Obviously, it will move its weights only if the output bits of the two partners are identical. In this case, a repulsive step between $`E`$ and $`A`$ occurs with probability $`P_r=ϵ`$ where now $`ϵ`$ is the distance between the hidden units of $`E`$ and $`A`$.

Note that for both the partners and the attacker one has the important property that the networks remain identical after synchronization. When one has achieved $`ϵ=0`$ at some time step, the distance remains zero forever, according to the previous equations for $`P_r`$. However, although the attacker uses the same algorithm as the two partners, there is an important difference: $`E`$ can only listen but it cannot influence $`A`$ or $`B`$. This fact leads to the difference in the probabilities of repulsive steps; the attacker has always more repulsive steps than the two partners. For small distances $`ϵ1`$, the probability $`P_r`$ increases linear with the distance for the attacker but quadratic for the two partners. This difference between learning and listening leads to a tiny advantage of the partners over an attacker. The subtle competition between repulsive and attractive steps makes cryptography feasible.

On the other side, there is always a nonzero probability $`P_E`$ that an attacker will synchronize, too Mislovaty et al. (2002). For neural cryptography, $`P_E`$ should be as small as possible. Therefore it is useful to investigate synchronization for different models and to calculate their properties as a function of the model parameters.

Here we investigate a mechanism which decreases $`P_E`$, namely we include feedback in the neural networks. The input vectors $`𝐱_k`$ are no longer common random numbers, but they are produced by the bits of the corresponding hidden units. Therefore the hidden units of the two partners no longer receive an identical input, but two corresponding input vectors separate with the number of training steps. To allow synchronization, one has to reset the two inputs to common values after some time interval.

For nonzero distance $`ϵ>0`$, this feedback mechanism creates a sort of noise and increases the number of repulsive steps. After synchronization $`ϵ=0`$, feedback will produce only identical input vectors and the networks move with zero distance forever com .

Before we discuss synchronization and several attacking scenarios, we consider the properties of the bit sequence generated by a TPM with feedback.

## III Bit Generator

We consider a single TPM network with $`K`$ hidden units, as in the preceding section. We start with $`K`$ random input vectors $`𝐱_k`$. But now, for each hidden unit $`k`$ and for each time step $`t`$, the input vector is shifted and the output bit $`\sigma _k(t)`$ is added to its first component Eisenstein et al. (1995). Simultaneously, the weight vector $`𝐰_k`$ is trained according to the anti-Hebbian rule, Eq. (3). Consequently, the bit sequence $`\tau (t)`$ generated by the TPM is given by the equation

$$\tau (t)=\underset{k=1}{\overset{K}{}}\text{sgn}\left(\underset{j=1}{\overset{N}{}}w_{k,j}(t)\sigma _k(tj)\right).$$
(9)
Similar bit generators were introduced in Ref. Zhu and Kinzel (1998) and the statistical properties of their generated sequences were investigated Metzler et al. (2001). Here we study the corresponding properties for our TPM with discrete weights.

The TPM network has $`2^{KN}`$ possible input and $`(2L+1)^{KN}`$ weight vectors. Therefore our deterministic finite state machine can only generate a periodic bit sequence whose length $`l`$ is limited by $`(4L+2)^{KN}`$.

Our numerical simulations show that the average length $`l`$ of the period indeed increases exponentially fast with the size $`KN`$ of the network, but it is much smaller than the upper bound. For $`K=3`$ and $`L>N`$ we find $`l(2.69)^{3N}`$, independent of the number $`L`$ of weight values.

The network takes some time before it generates the periodic part of the sequence. We find that this transient time also scales exponentially with the system size $`KN`$. This means that, for sufficiently large values of $`N`$, say $`N100`$, any simulation of the bit sequence remains in the transient part and will never enter the cycle.

The bit sequence generated by a TPM with $`K>2`$ cannot be distinguished from a random bit sequence. For $`K=L=3`$ we have numerically calculated its entropy and found the value $`\mathrm{ln}2`$ as expected from a truly random bit sequence. In addition, we have performed several tests on randomness as described by Knuth Knuth (1981). We did not find any correlations between consecutive bits; the bit sequence passed all tests on randomness within strict confidence levels.

Although the bit sequence passed many known tests on random numbers we know that it is generated by a neural network. Does this knowledge help to estimate correlations of the sequence and to predict it? In fact, for a sequence generated by a perceptron (TPM with $`K=1`$), another perceptron trained on the sequence could achieve an overlap to the generator Metzler et al. (2000).

Consider a bit sequence generated by a TPM with the anti-Hebbian rule. Another TPM (the “student”) is trained on this sequence using the same rule. In addition, if the output bit disagrees with the corresponding bit of the sequence, we use the geometric method of Ref. Klimov et al. (2003) to perform a training step.

Figures 2 and 3 show that for $`K=3`$ hidden units, it is not possible to obtain an overlap to the generating TPM by learning the sequence. Only if the initial overlap between the generator and the student is very large there is a nonzero probability $`P_{sync}`$ that the student will synchronize with the generator. If it does not synchronize, the overlap between student and generator decays to zero.

Summarizing, a TPM network generates a pseudorandom bit sequences which cannot be predicted from part of the sequence. As a consequence, for cryptographic applications, the TPM can be used to encrypt and decrypt a secret message after it has generated a secret key.

## IV  Synchronization

As shown in the preceding section, a TPM cannot learn the bit sequence generated by another TPM since the two input vectors are completely separated by the feedback mechanism. This also holds for synchronization by mutual learning: With feedback, two networks cannot be attracted to an identical time dependent state. Hence, to achieve synchronization, we have to introduce an additional mechanism which occasionally resets the two inputs to a common vector. This reset occurs whenever the system has produced $`R`$ different output bits, $`\tau ^A(t)\tau ^B(t)`$. For $`R=0`$ we obtain synchronization without feedback, which has been studied previously, and for large values of $`R`$ the system does not synchronize. Accordingly, we have added a new parameter in our algorithm which increases the synchronization time as well as the difficulty to attack the system. In the following two sections, we investigate synchronization and security of the TPM with feedback quantitatively.

We consider two TPMs $`A`$ and $`B`$ which start with different random weights and common random inputs. The feedback mechanism is defined as follows.

1. After each step $`t`$ the input is shifted, $`x_{k,j}(t+1)=x_{k,j1}(t)`$ for $`j>1`$.
2. If the output bits agree, $`\tau ^A(t)=\tau ^B(t)`$, the output of each hidden unit is used as a new input bit, $`x_{k,1}(t+1)=\sigma _k(t)`$, otherwise all $`K`$ pairs of input bits $`x_{k,1}(t)`$ are set to common public random values.
3. After $`R`$ steps with different output, $`\tau ^A(t)\tau ^B(t)`$, all input vectors are reset to public common random vectors, $`x_{k,j}^A(t+1)=x_{k,j}^B(t+1)`$.

Feedback creates correlations between the weights and the inputs. Therefore the system becomes sensitive to the learning rule. We find that only for the anti-Hebbian rule, Eq. (3), the components of the weights have a broad distribution. The entropy per component is larger than 99% of the maximal value $`\mathrm{ln}(2L+1)`$. For the Hebbian or random walk rule, the entropy is much smaller, because the values of the weights are pushed to the boundary values $`\pm L`$. Therefore the network with the anti-Hebbian rule offers less information to an attack than the two other rules.

In Fig. 4 we have numerically calculated the average synchronization time as a function of the number $`L`$ of components for the anti-Hebbian rule. Obviously, there is a large deviation from the scaling law $`t_{sync}L^2`$ as observed for $`R=0`$. Our simulations for larger values of $`N`$, which are not included here, show that there exist strong finite size effects which do not allow to derive a reliable scaling law from the numerical data.

Fortunately, the limit $`N\mathrm{}`$ can be performed analytically. The simulation of the $`KN`$ weights is replaced by a simulation of an $`(2L+1)\times (2L+1)`$ overlap matrix $`f_{a,b}^k`$ for each hidden unit $`k`$ which measures the fraction of weights which are in state $`a`$ for the TPM $`A`$ and in state $`b`$ for $`B`$ Rosen-Zvi et al. (2002a, b).

We have extended this theory to the case of feedback. A new variable $`\lambda _k(t)`$ is introduced which is defined as the fraction of input components $`x_{k,j}`$ which are different between the corresponding hidden units of $`A`$ and $`B`$. This variable changes with time, and it influences the equation of motion for the overlap matrix $`f_{a,b}^k(t)`$. Details are described in the Appendix.

Figure 5 shows the results of this semianalytic theory. Now, in the limit of $`N\mathrm{}`$, the average synchronization time can be fitted to increase with a power of $`L`$, roughly proportional to $`L^2`$. The data indicate that only the prefactor but not the exponent depends on the strength $`R`$ of the feedback; the prefactor seems to increase linearly with $`R`$.

Hence, if the network is large enough, feedback has only a small effect on synchronization. In the following section we investigate the effect of feedback on the security of the network: How does the probability that an attacker is successful depend on the feedback parameter $`R`$?

## V Ensemble of attackers

Up to now, the most successful attack on neural cryptography is the geometric attack Klimov et al. (2003); Mislovaty et al. (2002). The attacker $`E`$ uses the same TPM with an identical training step as the two partners. That means, only for $`\tau ^A=\tau ^B`$ the attacker performs a training step. When its output bit $`\tau ^E`$ agrees with the two partners, the attacker trains the hidden units which agree with the common output. For $`\tau ^E\tau ^{A/B}`$, however, the attacker first inverts the output bit $`\sigma _k`$ for the hidden unit with the smallest absolute value of the internal field and then performs the usual training step.

For the geometric attack the probability $`P_E`$ that an attacker synchronizes with $`A`$ and $`B`$ is nonzero. Consequently, if the attacker uses an ensemble of sufficiently many networks there is a good chance that at least one of them will find the secret key.

We have simulated an ensemble of attackers using the geometric attack for the two TPMs with feedback and anti-Hebbian learning rule. Of course, each attacking network uses the same feedback algorithm as the two partner networks. Figure 6 shows the results of our numerical simulations. The success probability $`P_E`$ decreases with the feedback parameter $`R`$. For the model parameters shown in Fig. 6 we find that $`P_E`$ can be fitted to an exponential decrease with $`L^2`$,

$$P_Ee^{uL^2}.$$
(10)
The coefficient $`u`$ increases linearly with $`R`$, as shown in Fig. 7. The scaling \[Eq. (10)\], however, is a finite size effect. For large system sizes $`N`$, the success probability decreases exponentially with $`L`$ instead of $`L^2`$,

$$P_Ee^{yL}.$$
(11)
This can be seen from the limit $`N\mathrm{}`$ which can be performed with the analytic approach of the preceding section. Now the dynamics of the system is described by a tensor $`f_{a,b,e}^k`$ for the three networks $`A`$, $`B`$, and $`E`$ and corresponding variables $`\lambda _k^A,\lambda _k^B,\lambda _k^E`$. Details are given in the Appendix.

Figure 8 indicates the exponential scaling behavior \[Eq. (11)\] for several values of $`R`$. The coefficient $`y`$ increases linearly with $`R`$, as shown in Fig. 9.

These results show that feedback improves the security of neural cryptography. The synchronization time, on the other side, increases, too. Does the security of the system improve for constant effort of the two partners?

This question is answered in Fig. 10 which shows the probability $`P_E`$ as a function of the average synchronization time, again for several values of the feedback parameter $`R`$. On the logarithmic scale shown for $`P_E`$, the security does not depend much on the feedback. For constant effort to find the secret key, feedback yields a small improvement of security, only.

## VI Conclusions

Neural cryptography is based on a delicate competition between repulsive and attractive stochastic forces. A feedback mechanism has been introduced which amplifies the repulsive part of these forces. We find that feedback increases the synchronization time of two networks and decreases the probability of a successful attack.

The numerical simulations up to $`N=10^5`$ do not allow to derive reliable scaling laws, neither for the synchronization time nor for the success probability. But the limit $`N\mathrm{}`$ which can be performed analytically indicates that the scaling laws with respect to the number $`L`$ of component values are not changed by the feedback, only the respective coefficients are modified. The average synchronization time increases with $`L^2`$ while the probability $`P_E`$ of a successful attack decreases exponentially with $`L`$, for huge system sizes $`N`$.

Accordingly, the security of neural cryptography is improved by including feedback in the training algorithm. But simultaneously the effort to find the common key rises. We find that for a fixed synchronization time, feedback yields a small improvement of security, only.

After synchronization, the system is generating a pseudorandom bit sequence which passed all tests on random numbers applied so far. Even if another network is trained on this bit sequence it is not able to extract some information on the statistical properties of the sequence. Consequently, the neural cryptography cannot only generate a secret key, but the same system can be used to encrypt and decrypt a secret message, as well.

*

## Appendix A Semianalytical calculation for synchronization with feedback

In this appendix we describe our extension of the semianalytic calculation Rosen-Zvi et al. (2002a, b) to the case of feedback.

The effect of the feedback mechanism depends on the fraction $`\mathrm{\Lambda }`$ of newly generated input elements $`x_{k,j}`$ per step and hidden unit. In the numerical simulations presented in this paper $`\mathrm{\Lambda }`$ is equal to $`N^1`$. In this case the effect of the feedback mechanism vanishes in the limit $`N\mathrm{}`$. But it is also possible to generate several input elements $`x_{k,j}`$ per hidden unit and step. For that purpose one can multiply the output bit $`\sigma _k`$ with $`\mathrm{\Lambda }N`$ random numbers $`z\{1,+1\}`$. As we want to compare the results of the semianalytical approach with simulations for $`N=1000`$, we set $`\mathrm{\Lambda }=10^3`$ in the following calculations.

In the case of two TPMs the development of the input noise $`\lambda _k`$ is given by

$$\lambda _k^+=(1\mathrm{\Lambda })\lambda _k+\mathrm{\Lambda }\mathrm{\Theta }(\sigma _k^A\sigma _k^B)\mathrm{\Theta }(\tau ^A\tau ^B).$$
(12)
At the beginning and after $`R`$ steps with $`\tau ^A\tau ^B`$ all variables $`\lambda _k`$ are set to zero (according to the algorithm described in Sec. IV).

The input noise generated by the feedback mechanism affects the output of the hidden units. An input element with $`x_{k,j}^B=x_{k,j}^B`$ causes the same output $`\sigma _k^B`$ as a change of sign in $`w_{k,j}^B`$ together with equal inputs for both $`A`$ and $`B`$. Therefore the probability $`ϵ_{k,\mathrm{eff}}`$ that two hidden units with overlap $`\rho _k`$ and input error $`\lambda _k`$ disagree on the output bit is given by

$$ϵ_{k,\mathrm{eff}}=\frac{1}{\pi }\mathrm{arccos}(12\lambda _k)\rho _k.$$
(13)
The distance $`ϵ_{k,\mathrm{eff}}`$ between the hidden units of $`A`$ and $`B`$ is used to choose the output bits $`\sigma _k^A`$ and $`\sigma _k^B`$ with the correct probabilities in each step Rosen-Zvi et al. (2002a, b).

The feedback mechanism influences the equation of motion for the overlap matrix $`f_{a,b}^k`$, too. Here we use additional variables $`\mathrm{\Delta }_k^m=\mathrm{\Theta }(\sigma _k^m\tau ^m)\mathrm{\Theta }(\tau ^A\tau ^B)`$ to determine if the weights of hidden unit $`k`$ in the TPM of $`m\{A,B,E\}`$ change ($`\mathrm{\Delta }_k^m=1`$) or not ($`\mathrm{\Delta }_k^m=0`$). Therefore we are able to describe the update of elements $`f_{a,b}^k`$ away from the boundary ($`L<a,b<L`$) in only one equation:

$`f_{a,b}^{k+}`$ $`=`$ $`{\displaystyle \frac{1\lambda _k}{2}}(f_{a+\mathrm{\Delta }_k^A,b+\mathrm{\Delta }_k^B}^k+f_{a\mathrm{\Delta }_k^A,b\mathrm{\Delta }_k^B}^k)`$ (14)
$`+`$ $`{\displaystyle \frac{1}{2}}\lambda _k(f_{a+\mathrm{\Delta }_k^A,b\mathrm{\Delta }_k^B}^k+f_{a\mathrm{\Delta }_k^A,b+\mathrm{\Delta }_k^B}^k).`$
The second term in Eq. (14) which is proportional to $`\lambda _k`$ shows the repulsive effect of the feedback mechanism. Similar equations can be derived for elements on the boundary.

In the limit $`N\mathrm{}`$ the number of steps required to achieve full synchronization diverges Kinzel and Kanter (2002a). Because of that one has to define a criterion which determines synchronization in order to analyze the scaling of $`t_{sync}`$ using semianalytic calculations. As in Ref. Rosen-Zvi et al. (2002a) we choose the synchronization criterion $`\overline{\rho }^{AB}=\frac{1}{3}_{k=1}^K\rho _k^{AB}0.9`$.

In order to analyze the geometric attack in the limit $`N\mathrm{}`$ one needs to extend the semianalytical calculation to three TPMs. In this case the development of the input noise is given by the following equations:

$`\lambda _k^{A+}`$ $`=`$ $`\mathrm{\Lambda }\mathrm{\Theta }(\sigma _k^A\sigma _k^B)\mathrm{\Theta }(\sigma _k^A\sigma _k^E)\mathrm{\Theta }(\tau ^A\tau ^B)`$ (15)
$`+`$ $`(1\mathrm{\Lambda })\lambda _k^A,`$
$`\lambda _k^{B+}`$ $`=`$ $`\mathrm{\Lambda }\mathrm{\Theta }(\sigma _k^B\sigma _k^A)\mathrm{\Theta }(\sigma _k^B\sigma _k^E)\mathrm{\Theta }(\tau ^A\tau ^B)`$ (16)
$`+`$ $`(1\mathrm{\Lambda })\lambda _k^B,`$
$`\lambda _k^{E+}`$ $`=`$ $`\mathrm{\Lambda }\mathrm{\Theta }(\sigma _k^E\sigma _k^A)\mathrm{\Theta }(\sigma _k^E\sigma _k^B)\mathrm{\Theta }(\tau ^A\tau ^B)`$ (17)
$`+`$ $`(1\mathrm{\Lambda })\lambda _k^E.`$
Analogical to Eq. (13) the distance $`ϵ_{k,\mathrm{eff}}^{mn}`$ between two hidden units can be calculated from the overlap $`\rho _k^{mn}`$ and the variables $`\lambda _k^m`$ and $`\lambda _k^n`$:

$$ϵ_{k,\mathrm{eff}}^{mn}=\frac{1}{\pi }\mathrm{arccos}(12\lambda _k^m2\lambda _k^n)\rho _k^{mn}.$$
(18)
But for the geometric attack the attacker $`E`$ needs to know the local fields $`h_k^E`$. The joint probability distribution of $`h_k^A`$, $`h_k^B`$ and $`h_k^E`$ is given by Rosen-Zvi et al. (2002a)

$$P(h_k^A,h_k^B,h_k^E)=\frac{e^{(1/2)(h_k^A,h_k^B,h_k^E)𝒞_k^1(h_k^A,h_k^B,h_k^E)^T}}{\sqrt{(2\pi )^3det𝒞_k}}.$$
(19)
The covariance matrix in this equation describes the correlations between the three neural networks:

$$𝒞_k=\left(\begin{array}{ccc}Q_k^A& R_{k,\mathrm{eff}}^{AB}& R_{k,\mathrm{eff}}^{AE}\\ R_{k,\mathrm{eff}}^{AB}& Q_k^B& R_{k,\mathrm{eff}}^{BE}\\ R_{k,\mathrm{eff}}^{AE}& R_{k,\mathrm{eff}}^{BE}& Q_k^E\end{array}\right).$$
(20)
From the tensor $`f_{a,b,e}^k`$ and the variables $`\lambda _k^m`$ one can easily calculate the elements of $`𝒞_k`$:

$`Q_k^A`$ $`=`$ $`{\displaystyle \underset{a,b,e=L}{\overset{L}{}}}a^2f_{a,b,e}^k,`$ (21)
$`Q_k^B`$ $`=`$ $`{\displaystyle \underset{a,b,e=L}{\overset{L}{}}}b^2f_{a,b,e}^k,`$ (22)
$`Q_k^E`$ $`=`$ $`{\displaystyle \underset{a,b,e=L}{\overset{L}{}}}e^2f_{a,b,e}^k,`$ (23)
$`R_{k,\mathrm{eff}}^{AB}`$ $`=`$ $`(12\lambda _k^A2\lambda _k^B){\displaystyle \underset{a,b,e=L}{\overset{L}{}}}abf_{a,b,e}^k,`$ (24)
$`R_{k,\mathrm{eff}}^{AE}`$ $`=`$ $`(12\lambda _k^A2\lambda _k^E){\displaystyle \underset{a,b,e=L}{\overset{L}{}}}aef_{a,b,e}^k,`$ (25)
$`R_{k,\mathrm{eff}}^{BE}`$ $`=`$ $`(12\lambda _k^B2\lambda _k^E){\displaystyle \underset{a,b,e=L}{\overset{L}{}}}bef_{a,b,e}^k.`$ (26)
We use a pseudorandom number generator to determine the values of $`h_k^A`$, $`h_k^B`$, and $`h_k^E`$ in each step. The application of the *rejection method* Hartmann and Rieger (2002) ensures that the local fields have the right joint probability distribution $`P(h_k^A,h_k^B,h_k^E)`$. Then the output bits $`\sigma _k^m`$ of the hidden units are given by $`\sigma _k^m=\mathrm{sgn}(h_k^m)`$. If $`\tau ^A=\tau ^B\tau ^E`$ the hidden unit $`k`$ with the smallest absolute local field $`|h_k^E|`$ is searched and its output $`\sigma _k^E`$ is inverted (geometric attack). Afterwards the usual training of the neural networks takes place.

The equation of motion for tensor elements $`f_{a,b,e}^k`$ away from the boundary ($`L<a,b,e<L`$) is given by

$`f_{a,b,e}^{k+}`$ $`=`$ $`{\displaystyle \frac{1\lambda _k^A\lambda _k^B\lambda _k^E}{2}}f_{a+\mathrm{\Delta }_k^A,b+\mathrm{\Delta }_k^B,e+\mathrm{\Delta }_k^E}^k`$ (27)
$`+`$ $`{\displaystyle \frac{1\lambda _k^A\lambda _k^B\lambda _k^E}{2}}f_{a\mathrm{\Delta }_k^A,b\mathrm{\Delta }_k^B,e\mathrm{\Delta }_k^E}^k`$
$`+`$ $`{\displaystyle \frac{1}{2}}\lambda _k^Af_{a\mathrm{\Delta }_k^A,b+\mathrm{\Delta }_k^B,e+\mathrm{\Delta }_k^E}^k`$
$`+`$ $`{\displaystyle \frac{1}{2}}\lambda _k^Af_{a+\mathrm{\Delta }_k^A,b\mathrm{\Delta }_k^B,e\mathrm{\Delta }_k^E}^k`$
$`+`$ $`{\displaystyle \frac{1}{2}}\lambda _k^Bf_{a+\mathrm{\Delta }_k^A,b\mathrm{\Delta }_k^B,e+\mathrm{\Delta }_k^E}^k`$
$`+`$ $`{\displaystyle \frac{1}{2}}\lambda _k^Bf_{a\mathrm{\Delta }_k^A,b+\mathrm{\Delta }_k^B,e\mathrm{\Delta }_k^E}^k`$
$`+`$ $`{\displaystyle \frac{1}{2}}\lambda _k^Ef_{a+\mathrm{\Delta }_k^A,b+\mathrm{\Delta }_k^B,e\mathrm{\Delta }_k^E}^k`$
$`+`$ $`{\displaystyle \frac{1}{2}}\lambda _k^Ef_{a\mathrm{\Delta }_k^A,b\mathrm{\Delta }_k^B,e+\mathrm{\Delta }_k^E}^k.`$
Similar equations can be derived for elements on the boundary. An attacker is considered successful if one of the conditions $`\overline{\rho }^{AE}0.9`$ or $`\overline{\rho }^{BE}0.9`$ is achieved earlier than the synchronization criterion $`\overline{\rho }^{AB}0.9`$.