text:
# Center vortex model for the infrared sector of Yang-Mills theoryTalk presented by M. Engelhardt

## Abstract

A model for the infrared sector of $`SU(2)`$ Yang-Mills theory, based on magnetic vortices represented by (closed) random surfaces, is presented. The model quantitatively describes both confinement (including the finite-temperature transition to a deconfined phase) and the topological susceptibility of the Yang-Mills ensemble. A first (quenched) study of the spectrum of the Dirac operator furthermore yields a behavior for the chiral condensate which is compatible with results obtained in lattice gauge theory.

Diverse nonperturbative effects characterize strong interaction physics. Color charge is confined, chiral symmetry is spontaneously broken, and the axial $`U(1)`$ part of the flavor symmetry exhibits an anomaly. Various model explanations for these phenomena have been advanced; to name but two widely accepted ones, the dual superconductor mechanism of confinement, and instanton models, which describe the $`U_A(1)`$ anomaly and spontaneous chiral symmetry breaking. However, no clear picture has emerged which comprehensively describes infrared strong interaction physics within one common framework. The vortex model presented here  aims to bridge this gap. On the basis of a simple effective dynamics, it simultaneously reproduces the confinement properties of $`SU(2)`$ Yang-Mills theory (including the finite-temperature deconfinement transition), as well as the topological susceptibility, which encodes the $`U_A(1)`$ anomaly. Furthermore, a first (quenched) study of the spectrum of the Dirac operator yields a behavior for the chiral condensate which is compatible with results obtained in lattice gauge theory; this indicates that the model also correctly reproduces the spontaneous breaking of chiral symmetry.

Center vortices are closed chromomagnetic flux lines in three-dimensional space; thus, they are described by closed two-dimensional world-surfaces in four-dimensional space-time. In the $`SU(2)`$ case, their magnetic flux is quantized such that they modify any Wilson loop by a phase factor $`(1)`$ when they pierce an area spanned by the loop. To arrive at a tractable vortex model, it is useful to compose the vortex world-surfaces out of plaquettes on a hypercubic lattice. The spacing of this lattice is a fixed physical quantity (related to a thickness of the vortex fluxes), and represents the ultraviolet cutoff inherent in any infrared effective framework. The model vortex surfaces are regarded as random surfaces, and an ensemble of them is generated using Monte Carlo methods. The corresponding weight function penalizes curvature by associating an action increment $`c`$ with every instance of two plaquettes which are part of a vortex surface, but which do not lie in the same plane, sharing a link (note that several such pairs of plaquettes can occur for any given link).

Via the definition given above, Wilson loops (and, in complete analogy, Polyakov loop correlators) can be evaluated in the vortex ensemble, and string tensions extracted. For sufficiently small curvature coefficient $`c`$, one finds a confined phase (non-zero string tension) at low temperatures, and a transition to a high-temperature deconfined phase. For $`c=0.24`$, the $`SU(2)`$ Yang-Mills relation between the deconfinement temperature and the zero-temperature string tension, $`T_C/\sqrt{\sigma _0}=0.69`$, is reproduced. When furthermore setting $`\sigma _0=(440\text{MeV})^2`$ to fix the scale, measurement of $`\sigma _0a^2`$ yields the lattice spacing $`a=0.39\text{fm}`$. The full temperature dependence of the string tensions is displayed in Fig. 1.

Note that the confined and deconfined phases can alternatively be characterized by certain percolation properties  of the vortex clusters.

Complementarily, also the topological properties of the Yang-Mills ensemble encode important nonperturbative effects. The topological charge $`Q`$ of a vortex surface configuration is carried by its singular points , i.e. points at which the set of tangent vectors to the surface configuration spans all four space-time directions (a simple example are surface self-intersection points). Since a vortex surface carries a field strength characterized by a nonvanishing tensor component associated with the two space-time directions locally orthogonal to the surface , these singular points are precisely the points at which the topological charge density $`ϵ_{\mu \nu \lambda \tau }\text{Tr}F_{\mu \nu }F_{\lambda \tau }`$ is non-vanishing. Note that, in order to carry nontrivial global topological charge, surfaces must be non-oriented; lines along which the orientation flips can be associated with Abelian magnetic monopoles . Generic vortex surfaces in the ensemble studied here are indeed non-orientable. In practice, identifying all singular points of the hypercubic lattice surfaces used in the present model involves resolving ambiguities  reminiscent of those contained in lattice Yang-Mills link configurations. The resulting topological susceptibility $`\chi =Q^2/V`$, where $`V`$ denotes the space-time volume under consideration, is exhibited in Fig. 2.

A comprehensive description of the nonperturbative phenomena which determine strong interaction physics must furthermore include the coupling of the vortices to quark degrees of freedom and the associated spontaneous breaking of chiral symmetry. The latter can be quantified via the chiral condensate $`\overline{\psi }\psi `$, which is related to the spectral density $`\rho (\lambda )`$ of the Dirac operator in a vortex background via the Casher-Banks formula $`\overline{\psi }\psi =\pi \rho (0)`$. One can locally associate the chromomagnetic flux represented by an arbitrary vortex surface with a continuum gauge field ; this allows to construct the Dirac operator directly in the continuum, preserving exact chiral symmetry. Globally, one must allow for the gauge field to be defined on different space-time patches, in order to accomodate the non-orientability of the vortex surfaces. While the gauge fields on each patch can be chosen Abelian (i.e. in 3-direction in color space), the transition functions in the overlap regions in general have to be non-Abelian. Since the spectrum of the Dirac operator here is obtained using the finite element method, it is natural to use the different finite elements as the space-time patches. A sample Dirac spectrum generated using this method is displayed in Fig. 3.

Note the anomalous enhancement of $`\rho (\lambda )`$ near $`\lambda =0`$, in accordance with the divergence expected for a Dirac operator with good chiral properties . Also the linear extrapolation of the bulk of the spectrum to $`\lambda =0`$, i.e. the spectral density obtained by truncating the $`\lambda =0`$ divergence, quantitatively lies in the range obtained in full $`SU(2)`$ Yang-Mills theory . Of course, the chiral condensate by itself does not directly represent a renormalization group invariant, physical quantity, but can only interpreted in conjunction with a definition of the current quark masses. These can only be reliably fixed in the present model after more detailed (ultimately unquenched) measurements of hadronic properties.

It is also instructive to consider a change in the construction of the gauge field which trivializes the topology. This can be done by replacing the monopole loops on the vortex surfaces by double vortex sheets spanning areas bounded by the loops. Then, the surfaces can be oriented, the global topological charge vanishes, and presumably also many local clusters of nontrivial topology are trivialized (although the topological density as such does not necessarily vanish). In the Dirac spectrum, this leads to a truncation of the $`\lambda =0`$ divergence, and instead the finite-volume gap near $`\lambda =0`$ familiar e.g. from chiral random matrix theory opens, without the bulk of the spectrum being qualitatively changed. While the $`\lambda =0`$ divergence thus seems intimately related to topological properties , chiral symmetry breaking as such merely appears to require randomness of the gauge fields in a very general sense, without specific reference to topology.

Several improvements of the present model are envisaged. For one, the treatment must be generalized from $`SU(2)`$ to $`SU(3)`$ color. Also the Dirac operator can still be developed further to explicitly take into account a finite thickness of the vortices; the above construction, while embodying the correct topology, still deals with magnetic flux localized on infinitely thin world-surfaces. Eventually, it is hoped that this model will become a useful tool for phenomenological considerations.


format: text

text:
# The application of computational mechanics to the analysis of geomagnetic data

## I Introduction

Computational mechanics (CM) Badii1999  has a formalism Crutchfield1989  that has been proved to construct the minimal model capable of statistically reproducing all the resolvable causal structure of any infinite sequence of discrete measurements (be they scalar, vector, tensor, or descriptive)Shalizi2001 . The size of a model so defined, measured by a quantity termed statistical complexity, $`C_\mu `$, Crutchfield1989  is a reliable and falsifiable indication of the amount of structure the data containCrutchfield1994 ; Shalizi2000 .

The particular strengths of this approach are that it enables the complexities and structures of different sets of data to be quantifiably compared and that it directly discovers detailed causal structure within those data. By examining data in this way it is possible to appreciate, in a well-defined abstract sense, how a system actually functions and what scales are most important to it. This information can then be used to optimise the efficiency of physically plausible modelsPalmer2000 .

As with all other analytical tools, CM has some limitations in the face of certain real-world problems that affect the information content of the signal under study. These problems may include:

1. Gaps in the data
2. Noise
3. Restricted sequence length
4. Correlations at a very wide range of scales.

The problem of correlations at a wide range of scales is particularly interesting and relevant to geophysical and other natural time series because of their typically power law (coloured noise) Fourier spectra pt1999 . Theoretically, the minimum resolvable scale will be constrained by the data sampling interval and the maximum resolvable scale by the length of the data series. In practice, the range of resolvable scales will also be set by the available computational resources. Thus, the range of resolvable scales may be less than those necessary to evaluate correlations on all relevant scales. Consequently, it is important to understand how structural analysis is affected by unresolved structure due to correlation.

In this paper, we address these issues in detail. In section II, we discuss how the ideal formalism of CM can be adapted to apply to a non-infinite series of corrupted and correlated data. In particular, three concepts are defined and discussed: (1) A tolerance parameter Crutchfield1989  to account for the statistical uncertainty introduced by a non-infinite series that destroys the exact equivalence of different causal states sharing the same outcome. (2) A new expletive filter that removes signal corruption by assuming that corruption creates rare causal states or words that are not in the dictionary of the true signal. (3) The new concept of effective soficity in which a data series has a finite set of equivalent causal states that is stable to small changes in the effective memory of those states.

The latter concept distinguishes between unresolvable “random” structure and resolvable structure whose discovery is only prevented by the effective memory being used and by the length of the data series.

In section III, we apply the CM algorithm with these additional concepts to the analysis of structure in four simulated time series: (1) Uncorrelated, white noise. (2) Periodic signal with white noise corruption. (3) A biased Poisson switch (i.e., a sequence of pulses whose pulse durations and inter-pulse intervals are determined by stationary Poisson processes). (4) A sequence of bursts similar to (3) but with fixed pulse duration. The structure of the time series is analysed by searching for regions of effective soficity in maps of statistical complexity over the parameter space of the CM model.

The simulated time series represent four types of signal thought to be present in time series measurements of the geomagnetic field. In section IV, we use the CM algorithm to examine a real geomagnetic time series measured at Halley, Antarctica, in which deflections of the earth’s magnetic field are due mainly to electrical currents in the ionosphere. The CM analysis yields a structural model that comprises a diurnal component corresponding to the oscillation of the measuring apparatus with the rotation of the earth and a Poisson-switched, fixed-duration, pulse component that is likely associated with the magnetospheric substorm Borovsky1993 .

In section V, we discuss some general principles that have been learnt in applying CM to the analysis of structure in real data, and draw conclusions in section VI.

## II Method

Here we give an introduction to the practical use of CM in the analysis of real data. We concentrate only on describing in detail the formalism for the parsing structure that we have used in the analyses. For a fuller description of the potential intricacies of the method see referenceShalizi2001 . Defining some new terminology, we highlight the difficulties associated with analysing experimental data in this way, and explain solutions to these problems.

To start with, one has a set of measurements - either a spatial or temporal series where the separation between each point is known. The total time or length for which data exist is their span, $`S`$. After coarse-graining at a fixed scale $`s`$, the series has $`N=\frac{S}{s}`$ equally spaced measurements. Next, we digitise the signal amplitude. For reasons that will be apparent later, the number of possible digits should be low unless the series length is extremely large. The digitised sequence is then a concatenation of $`N`$ letters $`𝐥^N=\{l_0,l_1,\mathrm{},l_i,l_{N1}\}`$, where there are $`L`$ types of such letters, ranging from $`0,1,2,\mathrm{}`$ up to $`L1`$. In order to maximize the prior probable information content of the processed sequence, digitisation should normally be performed such that there are equal numbers of each letter present. For example, in the case of binarisation (where $`L=2`$), this would mean that the threshold for letter $`1`$ would be the median value of the data. It should be noted, though, that the best way to digitise the sequence is that which actually maximises the information content of the result; but that this cannot usually be guessed. Another approach that has been suggestedPalmer2000  is to use the formalism of Maximum Entropy.

The next step is to parse the sequence. One begins by composing words from each group of $`n`$ consecutive letters; the $`i`$th word, $`𝐖_i`$, is defined by:

$$𝐖_i=𝐥_i^n=\{l_i,l_{i+1},\mathrm{},l_{i+n1}\}$$
(1)
Thus there are $`L^n`$ possible words, each represented by a unique scalar, $`W_i`$:

$$W_i=\underset{j=0}{\overset{n1}{}}L^{(n1)j}\times l_{i+j}$$
(2)
The total number of words generated from the sample is $`N(n1)`$. We now introduce some terminology; any word $`𝐖_i`$ may be called a proword, $`𝐖_P`$ when followed by any word $`𝐖_{i+1}`$. This latter is called the epiword, $`𝐖_E`$. For this sentence we digress slightly to note that it may sometimes be beneficial to perform the initial digitisation on each separate block of data $`2n`$ letters long rather than the entire dataset.

We now proceed to capture causal structure in the word sequence by compiling a tally of epiwords following each proword. This means going through the sequence incrementing an array $`𝐓_{(W_P,W_E)}`$ accordingly. Representing summation over an index by its omission, we see that the total tally is $`T=N(2n1)`$. Thus, contracting over epiwords gives a tally of prowords only:

$$𝐓_{(W_P)}=\underset{W_E=0}{\overset{L^n1}{}}𝐓_{(W_P,W_E)}$$
(3)
and the fractional prevalence of each proword in the sequence is therefore contained in the vector

$$𝐏_{(W_P)}=\frac{𝐓_{(W_P)}}{T}$$
(4)
Finally, the fractional profile of each proword by epiword is given by the array

$$𝐏_{(W_EW_P)}=\frac{𝐓_{(W_P,W_E)}}{𝐓_{(W_P)}}$$
(5)
where the repeated indices in the division are not summed over. Given a particular proword, this tells us the likelihoods of transitions to the various epiwords.

The crux of the technique now lies in identifying prowords with equivalent epiword profiles. Such prowords are said to belong to the same “equivalence class” or “causal state” - i.e. they share statistically equivalent probabilistic futures (at the level of analysis one has been pursuing). The identification is made via an equivalence relation, denoted by $``$. For an infinite sequence, $``$ can demand exact correspondence between profiles, in which case it is always transitive (meaning $`AB,BCAC`$). In a practical situation, where even the finite length of the sequence introduces fluctuations in the calculated profilesfootnote1 , it is not possible to be so exact. We therefore introduce a tolerance parameter, $`\tau `$, within the bounds of which the profiles of words in the same equivalence class are allowed to vary: Two prowords, $`A`$ and $`B`$, are in the same equivalence class if, $``$ $`W_E`$;

$$\left|𝐏_{(W_EW_P=A)}𝐏_{(W_EW_P=B)}\right|\tau $$
(6)
where the large vertical bars signify absolute magnitude.

Although this appears to destroy the formal transitive property of $``$, it is postulated here that after some finite sequence length is surpassed, there always exists a definite range of values for $`\tau `$ within which the transitive property is observed correctly. The transitive property can be enforced where $``$ is not transitive — by grouping equivalence classes defined by $``$ that share at least one word.

Having identified the words lying within each equivalence class, a model which outputs a series of letters statistically equivalent to the original can be constructed. It is a particular strength of the technique that the model generated is always a minimal representation of the data’s statistical structure for the amount of memory the analysis employsShalizi2000 . The model is easiest to describe in terms of its representation as a labelled “diagraph”. Two very simple labelled diagraphs, with extracts from their outputs, are presented in Figure 1. A more complicated labelled diagraph, also representing a minimal model, is shown in Figure 2. Each diagraph comprises a node or nodes, indicated by a circle with a number in it, and lines joining one node to another or to itself. Each numbered node of a diagraph represents a causal state corresponding to each of the model’s equivalence classes, while each line (unidirectionally) joining two nodes is labelled with the string of letters (the word) that is output when that line is followed. In addition, each line is associated with a probability. The word output on going from one causal state to another is in the equivalence class of the future state. The probability of each word’s output may therefore be trivially given by $`𝐏_{(W_EW_P)}`$. It should be held in mind that only a subset of all possible labelled diagraphs represent minimal models. Even so, an arbitrary diagraph’s output can naturally be used to construct the appropriate minimal model.
Models like this are useful for three reasons:

1. Their minimality allows the structure of two sets of data to be directly compared.
2. Once a model has been synchronized with current data it optimizes one’s ability to forecast the behaviour of the system in the future.
3. The information concerning scales of causal structure in the data can be used to optimize the performance of more physically plausible models.

If a recursive decomposition is employedZurek1990 ; LamNaroditsky1992 ; Crutchfield1994 ; Hanson1993 ; Perry1999 ; Shalizi2000 ; Young1991 , diagraphs labelled with words to output can be reformed into equivalent diagraphs labelled with single letters. This is perhaps a mathematically aesthetic thing to do. However, because any real analysis is performed with an effective memory of only $`n`$ symbols, only the last $`n`$ symbols are of any use for prediction. This fact is not manifest in the single-letter diagraph obtained through the decomposition of the transition matrix, $`P_{(W_EW_P)}`$. Moreover, it is necessary to synchronise a single-letter diagraph to an input stream of data before it is of any use for the purpose of prediction. In this paper we concentrate only upon the identification of proword equivalence classes, because it is a powerful tool for pattern discovery in its own right.

A measure of the structure of such models is given by the statistical complexityCrutchfield1989 :

$$C_\mu \underset{i}{}P(C_i)\mathrm{log}_2P(C_i)$$
(7)
where logarithms are canonically taken to base $`2`$ and the prevalence $`P_{(C_i)}`$ of equivalence class $`i`$ is given by the sum of the prevalences of the words in that class. For example, the models represented by the labelled diagraphs in Figure 1 both have a statistical complexity of zero because they only have one causal state (and therefore one equivalence class) each. This is sensible because they both output noise. The model represented in Figure 2, though, has four causal states with equal prevalences and a correspondingly higher statistical complexity of two bits:

$$C_\mu =\underset{i=1}{\overset{4}{}}P(C_i)\mathrm{log}_2P(C_i)=4\times \frac{1}{4}\times \frac{\mathrm{ln}(\frac{1}{4})}{\mathrm{ln}(2)}=2$$
(8)
$`C_\mu `$ is extremely important, not only because it reflects the complexity of the system, but also because it does not converge until the data have been fully characterized. It is a hard fact that if the sequence length $`N`$ is too small, full characterization will not be possible. This is because fluctuations in the proword profiles will corrupt the identification of equivalence classes. In this text the resultant unresolvable structure is called noise. In contrast, resolvable but as yet unresolved structure is described as random Poirier2001 . Such random structure is likely to be encountered in analysing data sets with correlation lengths comparable to or exceeding the maximum wordlength. Making this distinction is very important, even though it is not possible to discern whether unresolved structure is noisy or random until further computation has resolved it. In other words, the data appear to be noisy until the series is found to be effectively sofic, at which point $`C_\mu `$ attains its correct value and the model is complete. We note that sofic sequences are those which still have a finite number of equivalence classes when $`N`$ is infinite and $`n`$ is semi-infinite (see Badii and PolitiBadii1999 , page 80 for a longer explanation). Effective soficity is here defined to mean that a sequence has equivalence classes that are stable to an increase in wordlength. Thus, a sequence could be effectively sofic at one range of wordlengths but not at another where either more or less structure is in the process of being identified.

Both random structure and noise will redistribute the original tally from what would be expected if only resolved structure was present, raising $`C_\mu `$ from the value corresponding to resolved structure alone and increasing the complexity of its model. A simple one-parameter model for the corruption process is to assume that the probability that any letter is corrupted to any other letter is $`\chi `$. Then the probability any letter stays as it is is $`\sigma =1\chi `$ and the corruption will have been governed by the redistribution function

$$𝐓_{(W_P^C,W_E^C)}^{\text{corrupt}}=\underset{W_P=0}{\overset{L^n1}{}}\underset{W_E=0}{\overset{L^n1}{}}P(\$𝐖_P=𝐖_P^C,\$𝐖_E=𝐖_E^C)𝐓_{(W_P,W_E)}^{\text{pure}}$$
(9)
where $`\$𝐀=𝐁^C`$ reads, the pure word $`𝐀`$, when corrupted by noise in a certain way, is identical to the corrupt word $`𝐁^C`$. The label “pure” implies effective soficity. Thus, assuming the corruption of prowords and epiwords are independent we have

$$P(\$𝐖_P=𝐖_P^C,\$𝐖_E=𝐖_E^C)=P(\$𝐖_P=𝐖_P^C)P(\$𝐖_E=𝐖_E^C)$$
(10)
where

$$P(\$𝐖=𝐖^C)=\underset{i=1}{\overset{n}{}}\left\{\sigma \delta (w_i=w_i^c)+\chi \delta (w_iw_i^c)\right\}$$
(11)
where $`\delta `$s is a Kronecker delta and the corruption of letters is assumed to be independent.

It happens that arbitrarily corrupted distributions can be uniquely deconvolved as long as one knows $`\chi `$, but this is not usually the case in an experimental situation. We have two alternative options. The first is to scan through $`\chi `$, deconvolving the proword prevalences each time. This will produce a drastic decrease in the statistical complexity at some point, signifying correct parameterization of $`\chi `$. A good guess for $`\chi `$ might be the first value which results in a single proword having a prevalence of zero.

Whilst the assumption of independent corruption of letters is likely to be a good model of noise, it is unlikely to be a good model of the uncharacterised correlated structure that we call random. Consequently, a second option is to ignore the details of any corruption and simply assume that the prevalence of any expletive (corrupted word) is below a certain expletive prevalence, $`x`$. We scan through $`x`$, eradicating any prowords whose prevalence is less than $`x`$, and recalculate $`C_\mu `$ each time. This procedure can alternatively be performed after the identification of preliminary equivalence classes to eradicate expletive equivalence classes. In any case, the approach can only work when the actual structure-to-noise ratio (SNR) is high enough to ensure that expletives are eradicated before meaningful words are. If the pure proword prevalence distribution is very uneven this method cannot work. In general, a combined method would probably be most successful - that is, where one first attempts the deconvolution after making some bold assumptions and then removes the resulting low-prevalence words completely. It is always possible to determine all the resolvable structure of a sequence for which the SNR is arbitrarily small, so $`C_\mu `$ is independent of SNR. Of course though, if the SNR is zero, so is $`C_\mu `$, because the model suddenly collapses to a single equivalence class.

Note that deconvolution can always be achieved by inversion of an assumed convolution matrix, but that this is not always easy. In particular, if one knew the actual matrix then the “noise” would not be noise at all, but resolved structure. The only deconvolution that is strictly necessary is that which removes the noise (unresolvable structure) from the signal. It should therefore assume that the redistribution is Gaussian. In practice though, some random (resolvable) structure may be so computationally difficult to identify that a messy deconvolution is required to remove it, allowing the analysis of more easily resolvable structure to proceed. It is admissible to remove expletives from the prevalence distribution because they destroy the effective soficity of the data.

It is instructive at this point to go through the uncertainties present in the profile and prevalence distributions. When the sequence length is large compared to $`L^n`$ the probability that any individual word has been corrupted is approximately $`\mathrm{\Delta }=n\chi `$. Following the definition of the prevalence distribution, we find that the uncertainty in the prevalence of a proword $`\mathrm{\Delta }P_{(W_P)}`$ is governed by an inequality:

$$\frac{\mathrm{\Delta }}{\sqrt{T}}<\mathrm{\Delta }P_{(W_P)}<\mathrm{\Delta }$$
(12)
where the lower limit corresponds to uncorrelated errors and the upper limit to systematic errors. We indeed expect the uncertainty to be somewhere in this range because the errors are due to unresolved structure. The uncertainty in the prevalence of a single epiword within a particular proword’s profile is expected to be greater:

$$\frac{\mathrm{\Delta }}{\sqrt{𝐓_{(W_P)}}}<\mathrm{\Delta }P_{(W_EW_P)}<\mathrm{\Delta }$$
(13)
These inequalities go some way to justifying the use of the blanket tolerance $`\tau `$ to identify the equivalence classes, because we know nothing about the nature of the errors. In some cases it is conceivable that $`\tau `$ would have to be scaled by $`1/𝐓_{(W_P)}^{}{}_{}{}^{m}`$, where $`0<m<\frac{1}{2}`$ in order to correctly identify equivalence relations between profiles. In such cases $`m`$ is an extra parameter.

## III Examples

We now turn to the analysis of test sets of data by the algorithm described in detail above. The test data represent signal types thought to be present in time series measurements of the geomagnetic field that we shall study in the next section.

### III.1 White noise.

A white noise (temporally uncorrelated) signal was generated by a sequence of 5000 independent samples from a uniform distribution and converted to binary by setting those values above the median to unity and those below the median to zero. Figure 3 shows the variation of statistical complexity, $`C_\mu `$, versus word length, $`n`$, and tolerance, $`\tau `$, for this signal. The absence of a plateau in this graph indicates that, for the range of memories (wordlengths) tested, the analysis does not discern any structure at all in the signal. The linear variation of $`C_\mu `$ with $`n`$ for $`\tau 0`$ represents models with as much arbitrariness as is possible at each level of memory used in the analysis. These models collapse to a single equivalence class as the tolerance parameter is increased. An increasing amount of tolerance is required for this collapse for increasing wordlength, as expected from equation 13. Thus, no complex models at all were constructed for this noisy sequence at any time during this analysis. This was expected; we would have been disappointed with the random number generator that was used to construct the sequence (the IDL “randomu” function, see alsoPress1996 ) if we had easily found correlations.

### III.2 Periodic signal with white noise corruption.

Figure 4 shows the result of the analysis on a binary period four sequence (i.e., 00110011…) of length 5000, where 10% of the bits have been randomly flipped. This graph has a stable, but rather jagged, plateau at $`C_\mu 2.8`$ which begins at wordlength 4 for tolerances in the range $`0.1\tau 0.24`$. This plateau corresponds to a group of models that capture the essential structure in the signal. In the absence of noise the statistical complexity of a binary period four signal should be $`C_\mu =2`$. The apparently anomalously high level of the plateau is caused by both the noise and the finite sequence length corrupting the identification of the equivalence classes. It is not entirely flat because the corruption is different at each value of wordlength and tolerance. In fact, there is a gentle downward trend which would converge to $`C_\mu =2`$ in the limit of the extra, spurious, states decreasing in prevalence at longer and longer wordlengths, if the sequence was long enough. Note that the gradient of the increase of statistical complexity with wordlength changes at a memory equal to half the period of the structure in this signal. This is the point at which the structure is first discovered: It is important to note that the convergence of $`C_\mu `$ is not immediate, suggesting an analogue of the Nyquist sampling theorem for CM. Note also steep drops in $`C_\mu `$ where previously distinguishable equivalence classes have suddenly collapsed together as the tolerance parameter $`\tau `$ exceeds some critical value, supporting our earlier postulate that there is a definite range of $`\tau `$ within which a true transitive property is observed.

Figure 5 shows results of a similar analysis on the same sequence, excepting that this time, words of prevalence less than $`x`$ (the expletive prevalence parameter) were eradicated from the probability distributions. $`x`$ was chosen to be 0.075 for this graph. As we can see, this approach was entirely successful in the respect that $`C_\mu `$ converges to a plateau for a broad range of the tolerance parameter $`\tau `$. A minimal model that was capable of outputting sequences with statistical structure identical to that characterized from the input was effectively constructed at every point on this plateau. The value of $`C_\mu `$ for periodic sequences was always found to directly reflect the amount of memory required by the system to produce such data: A sequence with a sole period, $`Q`$, has a statistical complexity of $`\mathrm{log}_2(Q)\text{ bits}`$ (in this case the period four signal has a statistical complexity of $`2.0\text{ bits}`$). Moreover, we can appreciate that the analysis only yields a convergent value after the wordlength has exceeded at least half the period of the sequence. More generally; convergence begins when an analysis first has greater memory than a system. If the system has certain structure with greater memory than it may be feasible to analyse, for example, a red noise signal that consists of many Fourier modes with a power law distribution of amplitudes and random phases, $`C_\mu `$ will not ever truly converge. However, there may be stages where the analysis has enough memory to identify some structure, and this is indicated by approximately flat regions, or, at the very least, dips in the gradient of $`C_\mu `$ with increasing wordlength.

### III.3 Biased Poisson switch.

We next turned to more detailed analyses of two other illustratively important diagraph’s outputs. The first we considered was the biased Poisson switchBendat1958 , represented as a labelled diagraph in Figure 6. The circled states, $`0`$ and $`1`$ may each generate either a one or a zero with the probabilities shown.

In the figure, $`\overline{\alpha }=1\alpha `$ and $`\overline{\beta }=1\beta `$. Note that when $`\alpha =\beta `$ the output sequence is no longer biased. It turns out that the values of $`C_\mu `$ we can derive for different values of $`\alpha `$ and $`\beta `$ provide some nice insights into the nature of information and the optimization of measurement processes. The measure has two distinct regimes: where $`\alpha +\beta =1`$, and where they do not. Since the diagraph only has two states, it is clear that as far prediction of the next epiword is concerned, only the last bit of any proword can ever matter. Therefore, all words usually separate into two equivalence classes (corresponding to odd and even words). If, however, $`\alpha +\beta =1`$ then $`\alpha =1\beta =\overline{\beta }`$ and $`\beta =\overline{\alpha }`$. This always results in the two equivalence classes collapsing into one, giving a statistical complexity of zero, corresponding to pure noise. This is appropriate because in this degenerate situation the possible outcomes of node 0 in Figure 6 are identical to those of node 1 and the diagraph collapses to a single state too (see inset), and can only produce noise anyway. If the diagraph does not collapse in this way there will always be two equivalence classes. Their prevalences are found to be $`1/(\frac{\alpha }{\beta }+1)`$ and $`1/(\frac{\beta }{\alpha }+1)`$. Thus, when $`\alpha +\beta 1`$, we have:

$`C_\mu `$ $`=`$ $`{\displaystyle \underset{i}{}}P_i\mathrm{log}_2\left(P_i\right)`$ (14)
$`=`$ $`{\displaystyle \frac{1}{\left[\frac{\alpha }{\beta }+1\right]}}\mathrm{log}_2\left({\displaystyle \frac{\alpha }{\beta }}+1\right)+{\displaystyle \frac{1}{\left[\frac{\beta }{\alpha }+1\right]}}\mathrm{log}_2\left({\displaystyle \frac{\beta }{\alpha }}+1\right)`$
and if $`\alpha +\beta =1`$, $`C_\mu `$ is always zero. A graph of this function is shown in Figure 7. Note that it always evaluates to unity when $`\alpha =\beta `$, except when $`\alpha `$ = $`\beta =\frac{1}{2}`$. If $`\alpha `$ does not equal $`\beta `$ (and $`\alpha +\beta 1`$) then it is less than unity. In fact, as the switch becomes more and more biased the statistical complexity goes down and down, reaching zero when only one digit is ever output. This is to be expected because a biased data set (e.g. more ones than zeros) is a symptom of an inefficient measurement apparatus: If one symbol is more prevalent than any other then the system is under-characterised by the alphabet in use. In the parlance of Shannon’s theory of communication this statistical complexity is equivalent to the maximum rate of information.

Given that the collapse discussed above takes a slice out of the graph in Figure 7, we would expect sequences generated by certain Poisson switches to be more difficult to characterize. For example, sequences produced by a switch with $`\alpha =\beta =0.49`$ have a statistical complexity of unity, but it is difficult to distinguish them from noise (where $`C_\mu =0`$) because they are so close to the collapse at $`\alpha =\beta =\frac{1}{2}`$. Such a sequence, 5000 symbols long, was analysed up to a wordlength of $`n=7`$ at 100 equal intervals between $`\tau =0.00`$ and $`\tau =0.01`$, without assuming any noise was present (i.e. $`\chi =0`$ and $`x=0`$). The results are presented in Figure 8. The plateau corresponding to the optimal model is that which has a statistical complexity of unity. We can see that it is difficult to construct this model because the plateau is radically constricted at higher wordlengths. On one side of it $`\tau `$ is too small to identify the equivalence classes, so every proword occupies its own equivalence class and $`C_\mu =n`$, its maximum value at any wordlength $`n`$. On the other side, $`\tau `$ is too large, so the two equivalence classes collapse together, producing degenerate models that would output noise. The $`C_\mu =1`$ plateau has a distinct end at $`n=5`$ because the sequence is not long enough to support analysis at a wordlength of $`n=6`$. At the latter wordlength statistical fluctuations in every proword profile mean that the correct classification of equivalence classes is no longer possible at any range of $`\tau `$. We are not too concerned about this here because we have already identified the optimal model which was stable from $`n=1`$ to $`n=5`$. In fact, a “more optimal” model would be able to predict the flipping of the switch itself to some extent. The construction of such a model would probably need a lot of computation and would probably require $`N`$ to be very large. These things depend on how random or noisy the switch is.

If a set of data is very complicated, no stable model might be identified before the wordlength becomes too large to be statistically supportable by the sequence length. The only solution is to gather more data. The alternative is to settle with models that are either inadequate or arbitrarily complicated. Although the latter models reproduce structure well (and are therefore most useful to engineers), studying them can reveal little about underlying processes. They are scientifically unaesthetic. In contrast, one can tell a lot about the intricacies of a system from the minimal adequate model associated with it at a certain level of analysis. This is the concern of scientists.

### III.4 Fixed pulse duration Poisson switch.

The next class of labelled diagraphs we consider produce binary sequences that are simple models of a process with bursts. These sequences have the structure of sustained switches - that is, when the switch is down it has a constant probability of switching up, and when up, it stays up for a fixed count, $`U`$. When the sequence is unbiased the up-switching probability is $`\frac{1}{U+1}`$. See Figure 9 for an example of this kind of labelled diagraph. The exact statistical complexities of such unbiased sustained switches are given by

$`C_\mu `$ $`=`$ $`\left[\left({\displaystyle \frac{U+1}{2U}}\right)\mathrm{log}_2\left({\displaystyle \frac{U+1}{2U}}\right)+(U1)\left({\displaystyle \frac{1}{2U}}\right)\mathrm{log}_2\left({\displaystyle \frac{1}{2U}}\right)\right]`$ (15)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{ln}(2)}}\left[\mathrm{ln}(2U)\left({\displaystyle \frac{U+1}{2U}}\right)\mathrm{ln}(U+1)\right]`$
We now investigate the practical analysis of a sequence one million binary symbols long that was produced by a sustained switch with $`U=4`$. The statistical complexities of the models constructed by the analysis are shown in Figure 10. It can be seen that the first convergent values are at wordlengths one greater than $`U`$. That is to say, good models can be constructed when the analysis first has a greater memory than the system. The plateau identifiable with a model of the form shown in Figure 9 begins at a wordlength of $`4`$ and extends laterally from $`\tau 0.02`$ to $`\tau 0.06`$. The remarkable thing about this plateau is that, although it is very flat, it is not entirely flat. It begins at $`C_\mu ^51.77069`$ which is significantly higher than the theoretical statistical complexity of: $`C_\mu =\frac{1}{\mathrm{ln}(2)}\left[\mathrm{ln}(8)\left(\frac{5}{8}\right)\mathrm{ln}(5)\right]1.54879`$ and subsequently oscillates around this value while it converges to it (e.g. $`C_\mu ^{10}1.52160`$). This behaviour is caused by the phase ambiguity due to the absence of information concerning the synchronisation of a burst when a word is composed entirely of “up” symbols. For example, at wordlength six, the profile of word $`63`$, (i.e., $`111111`$ in binary), is a superposition of the profiles of sequences like $`10[111111]`$, $`01[111111]`$ and $`11[111111]`$, from each of which it cannot be distinguished at that level of analysis. Therefore, in this case, the profile of word $`63`$ does not match that of any other word, and is allocated its own equivalence class. Although the prevalence of this word, and thence its class, is very low, it is sufficient to distort the statistical complexity.

In an analysis with recourse to infinite memory, the prevalence of an infinite sequence of “up” symbols is zero. Thus, the $`U`$ causal states of such a sequence would be correctly identified, and the statistical complexity of the model constructed would match exactly with the theoretical value. Of course, in practice no analysis can have infinite memory. If one wishes to retain optimal predictability of future data then it is necessary to accept whatever model is actually constructed by an analysis with finite memory.

## IV Analysing geomagnetic data

The test data examples analysed in the previous section represent signal types thought to be present in time series measurements of the geomagnetic field. If this is true, we may expect to see similar structure emerging from a CM analysis of a real geomagnetic time series.

The CM analysis detailed in section II was performed on 3-hour averaged measurements of the variation of the east-west component of the geomagnetic field at Halley, Antarctica, from three separate years: 24 February \- 16 December, 1995, 26 January - 28 December, 1998, and 2 January - 30 December 2000. A graph of the data from 26 January to 28 December 1998 is shown in Figure 11. It can be seen that the magnetic deflections have both a linear trend and a high frequency signal with an annual amplitude modulation that maximises in the austral summer. The linear trend is caused by the movement of the ice shelf upon which Halley is situated and was removed by subtracting the result of a linear regression for each of the three years. The detrended time series was then binarised with respect to the median, giving three sequences of 2352, 2688, and 2896 symbols, respectively. These series were then analysed up to a wordlength of 10 and with tolerances varying in 80 equal steps from 0.05 to 0.25. Words with prevalences less than $`x=0.004`$ were eradicated. The graph of statistical complexity, $`C_\mu `$, is shown in Figure 12. Two plateaus are evident, one at $`C_\mu 0.9`$, covering a wide range of tolerances and between word lengths of 1 and 3, and the other plateau at $`C_\mu 5.0`$, at the top left-hand corner of the graph, between tolerances of about 0.05 and 0.07 and at word lengths of 8 or more.

The convergence of statistical complexity at a word length of 8 corresponds to a time scale of 8 x 3 = 24 hours. Such a diurnal variation is well known and is primarily caused by the rotation of the observing station with the earth under the so-called SQ ionospheric current system that is driven by pressure gradients caused by solar heating and is thus fixed in the Sun-Earth frame Akasofu1972 . The variation can be seen in the raw data, as illustrated by plotting a typical month of Halley geomagnetic data in Figure 13. The associated ground magnetic variation has neither a pure sinusoidal shape nor a fixed period of exactly 24 hours, that is likely to account for the higher observed statistical complexity of $`C_\mu 5.0`$ compared to the expected $`C_\mu =3`$ for a pure binary period 8 signal.

The other plateau in Figure 12 at word lengths of 1 to 3 indicates the presence of significant structure at 3 to 9 h time scales. This plateau has a statistical complexity of approximately 0.9 and an overall structure similar to that of Figure 8, suggesting the possibility of some random pulse-like process. Such a possibility is intriguing because pulse-like geomagnetic perturbations on hour time scales (known as magnetic bays) are particularly prominent during the nighttime at high (auroral zone) latitudes and are associated with magnetospheric substorms Smith1999  whose occurrence has been argued to be a stationary Poisson process with mean recurrence time of 5 h Borovsky1993 . Figure 14 shows a single day of Halley geomagnetic data that illustrates the presence of such pulse-like disturbances on hour time scales sitting on top of the diurnal variation.

To investigate this further, an analysis was made of a 40-minute averaged time series of the east-west component of the geomagnetic field at Halley from 00:00 UT, 25 January, 1998 to 00:00 UT, 26 December, 1998. After removing the linear trend in the data due to the movement of the ice shelf, the time series was binarised with respect to the median, giving a sequence of 12011 symbols. The series was analysed up to a word length of 11 for tolerances in 60 equal steps between 0.00 and 0.15. Words with prevalences less than $`x=0.015`$ were eradicated. The graph obtained for $`C_\mu `$ is shown in Figure 15. The plateaus in this graph are stable to variation of $`x`$. The higher plateaus have models that are more useful for prediction of future data, if they are stable to an increase in the amount of data available to the analysis. The lower plateaus have models that show the most dominant structures – and are easier to understand and interpret physically.

It can be seen from the graph that, at a tolerance between $`\tau =0.12`$ and $`\tau =0.15`$, more structure is identified between wordlengths six and eight than it was possible to resolve with a memory of only five symbols. The model which corresponds to this plateau is represented, for a wordlength of seven, in Figure 16. The details of the model are in the Appendix. Comparing with Figure 6, the transitions between states 0 and 1 of this diagraph are an approximately Poisson-switched process with a timescale of about five hours. This value is given by the range of wordlengths capable of resolving this structure from the sequence within this range of $`\tau `$ ($`n=6,7\text{and}8`$); at $`n=7`$ the characteristic timescale is $`7\times 40\text{ minutes}4\text{ hours and }40\text{ minutes}`$.

It was thought that the other states and transitions in Figure 16 would be caused by the diurnal variation of the data alone. This was investigated by analysing, in exactly the same way, a pure binary sequence with a period of 36 symbols – corresponding to a period of one day if each symbol were to represent a 40-minute average. The principal transitions of the model constructed for this sequence are shown in the diagraph drawn in Figure 17. The structural similarities and differences between this diagraph and the one in Figure 16 are obvious, and support the idea that the transitions between states 0 and 1 of Figure 16 are due to substorm activity, rather than merely being an artifact of a partially characterised 24-hour period.

## V Discussion

In the previous sections, we have demonstrated how CM can measure the statistical complexity of linear data sequences and construct the minimal model necessary to describe the data. The reader may have noticed that there are seven degrees of freedom in making such a model:

1. digitisation method (binary, trinary, etc)
2. coarse-graining scale, $`s`$
3. sequence length, $`L`$
4. wordlength, $`n`$
5. tolerance, $`\tau `$
6. corruption frequency, $`\chi `$
7. expletive frequency, $`x`$

These degrees of freedom express the level of information in the data and the depth of knowledge with which the model is probing the system from which the data are measured. For example, increasing the sequence length, $`L`$, reducing the coarse-graining scale, $`s`$, or increasing the digitisation from binary to trinary, all provide increased information and thereby increased knowledge of the system that the data represent. Conversely, increasing the tolerance or the expletive frequency reduces information by admitting different states to be equivalent or to be omitted, respectively, thereby reducing knowledge of the system. Consequently, we might anticipate that the best model of the system is the model corresponding to the region of the multi-dimensional parameter space in which information is maximised. Whilst such a model is the most accurate description of the data sequence with the greatest information content, it is not necessarily the optimal model of the system. This is because any data sequence is not a complete representation of the system it is measured from. In particular, it is limited in two important respects: First, there is structure in a data sequence, that we have termed noise, that cannot be resolved under any amount of computation. This will create differences in the profiles of words that are statistically insignificant and should be ignored by allowing some non-zero value of tolerance, corruption frequency or expletive frequency. Second, there is structure in a data sequence, that we have termed random, that has not been resolved at a certain level memory or wordlength but that is resolvable at a greater wordlength. Recognising these sources of structure, we advance a hypothesis about model construction:

> Given enough data relevant to a system, there exists at least one connected region of finite size in the space of parameters used to construct models from those data within which an optimal model of the system is constructed.

In other words, meaningful models of the data can only be found within certain, usually finite, zones of the parameter spacefootnote2 . Within each zone, $`C_\mu `$ is constant and the model is both stable and minimal. Outside this zone, the model is either too degenerate or overly complicated. For example, it will be degenerate (and $`C_\mu `$ will be too low) if $`\tau `$ is set too large. This is because equivalence classes will collapse into one another. Similarly, the model will be unnecessarily complicated (and $`C_\mu `$ will be too high) if $`\tau `$ is set too small. This is because distinctions will be made between words on the basis of insignificant differences in their profiles.

An analogy is the construction of a vocabulary for the structure of speciation of feline animals. If one is too fussy about the tail, Manx cats cannot be classed as domestic cats. If one’s sole criterion is purring or a meow, a lion cub may be misclassed as a domestic cat. The correct classification of feline animals needs a finite amount of information to fall within the boundaries of a finite number of provisos.

In the case of Computational Mechanics we interpret effectively sofic models to be optimal. Thus we seek plateaus in the multi-dimensional parameter space. Generally, this space can contain many plateaus, the heights of which are the corresponding models’ statistical complexities. If we want to forecast the data most accurately, we are looking for the highest plateau, which has the most stringent conditionsfootnote3 . More physically understandable models may exist on some lower plateaus where only the more dominant causal structures are preserved.

Thus, in the end, the success of the analysis depends upon the existence of effectively sofic plateaus of statistical complexity in the multi-dimensional parameter space and our ability to discover them. This is contingent upon the data that are supplied and how much computing power is available. It is important to bear in mind that the data are not only a function of the physical system’s behaviour, but also of the measurement apparatus and any pre-processing. There are four main pitfalls (represented by corresponding model parameters):

1. Mischaracterization of the system by the measurement apparatus ($`\mathrm{}\chi ,x`$)
2. Degradation of data prior to the analysis by processing ($`s`$)
3. Insufficient data to resolve all structure present ($`L`$)
4. Insufficient computing power to resolve random structure ($`n,\tau `$)

The apparatus may easily mischaracterize the system, either by introducing structure to the data which is foreign to the system’s behaviour, or by neglecting to transcribe structure that should be present. This situation is most apparent when the apparatus is clearly only taking measurements from a cross-section of the system. Nevertheless, if it is reasonable to assume in a particular case that the apparatus is capable of providing a good representation, then identified structure can be attributed to the system. In such cases we would also expect the statistical complexity to scale with the system’s true complexity. Naturally, this is not valid when the cross-section happens to be an exact sub-system.

Although all processing degrades data, it may still be possible to correctly characterize all the structure present. This is because the degradation will usually produce noise, which can be ignored. A graver problem is when (uncharacterizable) noise represents some of the system’s structure. The only solution may be to collect more data, but other preliminary approaches are to use a finer scale when coarse-graining and/or to digitise more finely. However, it is always necessary to choose sensible margins for the parameter search because some regions of the parameter space are computationally very costly to explore. For example, a trinary sequence is about seven hundred times as hard to fully analyse at a word length of eight than a binary sequence. You must have a good reason not to use binary.

An alternative approach may be useful when the data have a number of widely separated scales with structure; it may be more computationally efficient to construct higher-level equivalence classes than to persist with using longer and longer words. If this is the case there will be a change in the constant increase of $`C_\mu `$ with wordlength at the wordlength where the lower scale of structure is found, but the gradient will thereafter remain constant (until the next scale of structure is reached).

Classes on the next highest level are found by applying the same analysis method to the sequence expressed in terms of a set of primary level causal states for which $`C_\mu `$ has not yet converged. All information between the scales $`n_1s`$ and $`n_2n_1s`$ is lost in this process. Even so, it is a more preferable approach than simply further coarse-graining the data to intervals of $`n_1s`$ if one has reason to believe that the system’s degrees of freedom at the two scales are coupled. The total statistical complexity is the sum of those calculated at each level, so it is in fact possible to test for such coupling by comparing the coarse-grained $`C_\mu `$ with the hierarchical value.

There is actually no reason why the prowords and epiwords should not come from different sequences, enabling the direct causal correlation of two systems, such as the solar wind and the magnetosphere.

## VI Conclusion

Computational Mechanics is an intuitive and powerful way to study complicated linear outputs from physical systems. This is because the analysis identifies causal structure from data presented to it and constructs the minimal adequate model that fits these data. The information about this structure, and in particular its scales, can then be used to optimise more physically plausible models. In this paper we have discussed in detail how the original formalism has to be used when applied to non-infinite sequences. The main conclusion is that models constructed by Computational Mechanics are good if, and only if, they are stable to the variation of the parameters used to construct them from the data. In addition, various more general postulates and definitions are made. These concern the general constructibility of models from a set of observations:

1. Structure which cannot be resolved from a set of data under any amount of computation is most usefully called noise
2. Structure which has not been resolved at a certain level of computation or memory, but which is resolvable from the set of data is usefully called random
3. Given enough data relevant to a system, there exists at least one connected region of finite size in the space of parameters used to construct models from those data within which an optimal model of the system is constructed.

The prior undecidability of whether unresolved structure is noise or randomness is a direct parallelism of Gödel’s famous theorem. For a proof relating the two fields, but in a slightly different context, see G.J. ChaitinChaitin1974 .

The method developed in this paper was applied to magnetometer measurements of ionospheric currents for the years 1995, 1998 and 2000. The technique successfully constructed models, the simplest of which comprised a diurnal component and a Poisson-switched process with a timescale of about five hours that likely relates to the occurrence of magnetic substorms. The most complicated model could be used to forecast space weather.

A similar method was also proposed to characterize the causal relationship of any two systems, such as the solar wind and the magnetosphere.

*

## Appendix A Details of the geomagnetic data model

Details of the simple stable model at wordlength seven, $`\tau 0.14`$. $`x=0.015`$, at which value 114 of 128 words are cut.

| Class | Equivalent surviving words | | | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 0 | 0 | 64 | 96 | 112 | 120 | 124 | 126 |
| 1 | 1 | 3 | 7 | 15 | 31 | 63 |  |
| 2 | 127 |  |  |  |  |  |  |

| Transitions from Class 0 | | |
| --- | --- | --- |
| Word Number (Class) | Word | Probability |
| 0 (0) | 0000000 | 0.481382 |
| 1 (1) | 0000001 | 0.0759726 |
| 3 (1) | 0000011 | 0.0729505 |
| 7 (1) | 0000111 | 0.0641128 |
| 15 (1) | 0001111 | 0.0540634 |
| 31 (1) | 0011111 | 0.0546018 |
| 63 (1) | 0111111 | 0.0532191 |
| 64 (0) | 1000000 | 0.0338069 |
| 96 (0) | 1100000 | 0.0234980 |
| 127 (2) | 1111111 | 0.0514929 |

| Transitions from Class 1 | | |
| --- | --- | --- |
| Word Number (Class) | Word | Probability |
| 0 (0) | 0000000 | 0.0423272 |
| 64 (0) | 1000000 | 0.0291209 |
| 96 (0) | 1100000 | 0.0209377 |
| 112 (0) | 1110000 | 0.0209523 |
| 120 (0) | 1111000 | 0.0293505 |
| 124 (0) | 1111100 | 0.0442261 |
| 126 (0) | 1111110 | 0.0822315 |
| 127 (2) | 1111111 | 0.684994 |

| Transitions from Class 2 | | |
| --- | --- | --- |
| Word Number (Class) | Word | Probability |
| 0 (0) | 0000000 | 0.0789801 |
| 64 (0) | 1000000 | 0.0907960 |
| 96 (0) | 1100000 | 0.0945274 |
| 112 (0) | 1110000 | 0.103234 |
| 120 (0) | 1111000 | 0.108831 |
| 124 (0) | 1111100 | 0.105721 |
| 126 (0) | 1111110 | 0.101990 |
| 127 (2) | 1111111 | 0.268035 |

| Average Transitions from Class 0 | | |
| --- | --- | --- |
| Class | Average Word | Probability |
| 0 | \[ 0.161, 0.102, 0.061, 0.040, 0.027, 0.009, 0.000 \] | 0.573587 |
| 1 | \[ 0.000, 0.142, 0.288, 0.432, 0.603, 0.797, 1.000 \] | 0.374920 |
| 2 | \[ 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000 \] | 0.051493 |

| Average Transitions from Class 1 | | |
| --- | --- | --- |
| Class | Average Word | Probability |
| 0 | \[ 0.843, 0.735, 0.657, 0.579, 0.470, 0.306, 0.000 \] | 0.269146 |
| 1 | \[ 0.000, 0.309, 0.444, 0.590, 0.714, 0.906, 1.000 \] | 0.045860 |
| 2 | \[ 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000 \] | 0.684994 |

| Average Transitions from Class 2 | | |
| --- | --- | --- |
| Class | Average Word | Probability |
| 0 | \[ 0.885, 0.752, 0.614, 0.463, 0.304, 0.149, 0.000 \] | 0.684080 |
| 1 | \[ 0.000, 0.117, 0.221, 0.325, 0.532, 0.701, 1.000 \] | 0.047886 |
| 2 | \[ 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000 \] | 0.268035 |

| Class | Prevalence |
| --- | --- |
| 0 | 0.483737 |
| 1 | 0.269561 |
| 2 | 0.246701 |

Statistical Complexity $`=`$ 1.51477 bits


text:
# Bang-Bang Operations from a Geometric Perspective

## I Introduction

Recently, controlling the evolution of a system by using strong, short pulses has been introduced as a new means for quantum error correction/prevention Viola and Lloyd (1998); Vitali and Tombesi (1999); Zanardi (1999); Duan and Guo (1999). These operations have been termed “bang-bang” (BB) pulses Viola and Lloyd (1998) (a name derived from classical control theory Velimir Jurdjevic (1997)), “parity kicks” Vitali and Tombesi (1999); D. Vitali and P. Tombesi (2001) (for the special case of a sign changing operation), decoupling operations L. Viola, S. Lloyd and E. Knill (1999); L. Viola, E. Knill, and S. Lloyd (2000) (since they can serve to decouple the system from environmental degrees of freedom), and symmetrization procedures Zanardi (1999) (which are associated with a group symmetrization/averaging). The advantage they have over the active and passive error correction procedures associated with quantum error correcting codes (QECCs) (see P.W. Shor (1995); Gottesman (1997); A.M. Steane (1999) and references therein) and decoherence free subspaces (DFSs) (or noiseless subsystems; see J. Kempe, D. Bacon, D.A. Lidar, and K.B. Whaley (2001); E. Knill, R. Laflamme and L. Viola (2000); Zanardi (2001) and references therein) is the use of external pulses rather than requiring several physical qubits to encode one logical qubit. Since today’s experiments use $`<10`$ qubits, this may, for the time being, make BB controls a method of choice for small-scale quantum computer implementations. However, it is clear that the time constraints imposed by bang-bang operations on the system are great L. Viola, S. Lloyd and E. Knill (1999) and may not be practical for eliminating noise altogether. Even if this technique cannot completely eliminate the noise, it can still be used to reduce noise Viola and Lloyd (1998); Duan and Guo (1999); D. Vitali and P. Tombesi (2001). This is important for possibly reducing error rates and thus extending computing time and/or the utility of QECCs and/or DFSs.

In order to take full advantage of the BB technique, the symmetrization operations and their effects must be made clear so that the benefit from the implementation can be readily determined. The work put forth here will aid in the analysis of the results of BB operations by providing an explicit geometric representation for the group-algebraic elements describing such interactions. This geometric picture also has the advantage of clarifying the error between a desired and a modified evolution. In addition, for two-state systems (qubits), we recover a familiar Bloch sphere representation and this provides us with an intuitive understanding of BB/symmetrization operations.

## II Decoupling by Symmetrization

The process of decoupling by symmetrization, counteracts decoherence by applying sequences of frequent pulses Vitali and Tombesi (1999); Viola and Lloyd (1998); Duan and Guo (1999). The time scales are crucial: roughly speaking, one needs to perform a complete cycle of symmetrization operations in a time shorter than the bath correlation time. An elegant group theoretical treatment shows that if the applied pulses are unitary transformations forming a finite-dimensional group, then the application of that series of pulses amounts to an average (symmetrization) over this group L. Viola, E. Knill and S. Lloyd (1999); L. Viola, S. Lloyd and E. Knill (1999); L. Viola, E. Knill, and S. Lloyd (2000); Zanardi (1999, 2001). We briefly review this theory.

The general evolution of a system and a bath coupled to it can be written in the form

$$H=H_S1\mathrm{l}_B+1\mathrm{l}_SH_B+\underset{\gamma }{}S_\gamma B_\gamma ,$$
(1)
where $`H_S`$ acts on the system alone, $`H_B`$ acts on the bath alone, and $`H_I_\gamma S_\gamma B_\gamma `$ is the interaction part of the Hamiltonian composed of traceless operators $`S_\gamma `$ ($`B_\gamma `$) which act on the system (bath). The objective of the BB procedure is to modify this evolution.

A set of symmetrization or BB operations can be chosen such that they form a discrete (finite order) subgroup of the full unitary group of operations on the Hilbert space of the system. Denote this subgroup $`𝒢`$ and its elements $`g_k`$, $`k=0,1,\mathrm{},|𝒢|1`$, where $`|𝒢|`$ is the order of the group. The cycle time is $`T_c=|𝒢|\mathrm{\Delta }t`$, where $`|𝒢|`$ is now also the number of symmetrization operations, and $`\mathrm{\Delta }t`$ is the time that the system evolves freely between operations under $`U_0=\mathrm{exp}(iHt)`$. The symmetrized evolution is given by

$$U(T_c)=\underset{k=0}{\overset{|𝒢|1}{}}g_k^{}U_0(\mathrm{\Delta }t)g_ke^{iH_{eff}T_c}.$$
(2)
$`H_{eff}`$ denotes the resulting effective Hamiltonian. Since the approximation requires very strong, short pulses to be implemented in a sequence, they have been termed bang-bang (BB) operations (we will use symmetrization and BB operations interchangeably). In this (BB) limit

$$HH_{eff}=\frac{1}{|𝒢|}\underset{k=0}{\overset{|𝒢|1}{}}g_k^{}Hg_k\mathrm{\Pi }_𝒢(H),$$
(3)
where $`H_{eff}`$ is the desired Hamiltonian (without noise). The map $`\mathrm{\Pi }_𝒢`$ is the projector into the centralizer, $`Z(𝒢)`$, defined as

$$Z(𝒢)=\{X|[X,g_k]=0,g_k𝒢\}.$$
(4)
It is clear that $`\mathrm{\Pi }_𝒢`$ commutes with all $`g_k`$ so that, if our group is generated by $`\{1\mathrm{l},H_S,S_\gamma \}`$, the system is effectively decoupled from its environment. The *control algebra* is the algebra generated by the set $`\{g_k\}`$. Even if the symmetrization is performed under less than ideal conditions, it can still reduce the noise in the system Viola and Lloyd (1998); Duan and Guo (1999); D. Vitali and P. Tombesi (2001).

## III Geometric Interpretation of the Effect of BB Operations

Now consider a set of unitary operators $`\{U_k\}`$, $`U_01\mathrm{l}_S`$, as an explicit realization of the subgroup $`𝒢`$ and choice of our set of BB operations. Then the following condition must be satisfied for an evolution generated by the effective Hamiltonian:

$$H_{eff}=\frac{1}{|𝒢|}\underset{k=0}{\overset{|𝒢|1}{}}U_k^{}HU_k.$$
(5)
Note that $`H_{eff}=0`$ is the case of *storage*. Considering Eq. (1), we can always include the terms $`H_S1\mathrm{l}_B`$ and $`1\mathrm{l}_SH_B`$ in $`H_I`$. We do not include the identity component $`1\mathrm{l}_S1\mathrm{l}_B`$ since it only gives rise to an overall phase. Thus $`H`$ and $`H_{eff}`$ are traceless. Let us now introduce $`Nn^21`$ traceless, Hermitian generators $`\{\lambda _i\}_{i=1}^N`$ of $`SU(n)`$. These generators are closed under commutation and span the space of traceless Hermitian matrices. For $`SU(2)`$, the Pauli matrices are commonly used; for $`SU(3)`$, the Gell-Mann matrices, and for higher dimensions, one may use a direct generalization of the Gell-Mann matrices. For dimensions that are a power of two it is often convenient to use the Pauli group (tensor products of Pauli matrices). The $`\{\lambda _i\}`$ satisfy trace-orthogonality,

$$\text{Tr}(\lambda _i\lambda _j)=M\delta _{ij},$$
(6)
where $`M`$ is a normalization constant (often taken to be $`2`$ for Lie algebras or $`n`$ for $`n\times n`$ matrices). Expanding the system operators in terms of the $`\{\lambda _i\}`$ yields:

$$S_\gamma =\underset{i}{}a_{i\gamma }\lambda _i$$
(7)
where the expansion coefficients are

$$a_{i\gamma }=\frac{1}{M}\text{Tr}(\lambda _iS_\gamma ).$$
(8)
Using this, $`H`$ can be written as as follows:

$$H=\underset{\gamma }{}S_\gamma B_\gamma =\underset{\gamma }{}\underset{i=1}{\overset{N}{}}a_{i\gamma }\lambda _iB_\gamma \underset{\gamma }{}(\stackrel{}{a}_\gamma \stackrel{}{\lambda })B_\gamma .$$
(9)
Here $`\stackrel{}{a}_\gamma `$ and $`\stackrel{}{\lambda }`$ are vectors of length $`N`$. In this representation, used extensively in G. Mahler and V.A. Weberruss (1998), an $`n\times n`$ Hamiltonian, $`H`$, is a vector with coordinates $`\stackrel{}{a}_\gamma `$ for each error $`\gamma `$ in an $`N`$-dimensional vector space spanned by the $`\{\lambda _i\}`$ as basis vectors, with ordinary vector addition and scalar multiplication.

As is well-known, there is a homomorphic mapping between the Lie groups $`SU(2)`$ and $`SO(3)`$ Michael Tinkham (1964). This mapping is generalized as follows for $`SU(n)`$ and a subgroup of the rotation group $`SO(N)`$:

$$U_k^{}\lambda _iU_k=\underset{j=1}{\overset{N}{}}R_{ij}^{(k)}\lambda _j,$$
(10)
where the matrix $`R^{(k)}SO(N)`$, the adjoint representation of $`SU(n)`$.

The BB operation \[Eq. (5)\] may now be viewed as a weighted sum of rotations of the (adjoint) vectors $`\stackrel{}{a}_\gamma `$. To see this, first let

$$\stackrel{}{a}_\gamma ^{(k)}=R^{(k)}\stackrel{}{a}_\gamma .$$
(11)
This represents the rotation by $`R^{(k)}`$ of the coordinate vector $`\stackrel{}{a}_\gamma `$. Next average over all rotations:

$$\stackrel{}{a}_\gamma ^{}=\frac{1}{|𝒢|}\underset{k=0}{\overset{|𝒢|1}{}}\stackrel{}{a}_\gamma ^{(k)}.$$
(12)
Finally, note that the effective Hamiltonian, after the BB operations, can be rewritten as:

$$H_{eff}=\frac{1}{|𝒢|}\underset{k=0}{\overset{|𝒢|1}{}}U_k^{}HU_k=\underset{\gamma }{}(\stackrel{}{a}_\gamma ^{}\stackrel{}{\lambda })B_\gamma .$$
(13)
Eq. (13) \[compare to Eq. (9)\] is our desired geometric representation of BB operations. Their effect is to simply transform, for each error $`\gamma `$, the coordinates $`\stackrel{}{a}_\gamma `$ to $`\stackrel{}{a}_\gamma ^{}`$. It is simplest to interpret this in the case of storage, where we seek BB operations such that $`H_{eff}=0`$. Since the errors can be decomposed in the linearly independent basis set indexed by $`\gamma `$, each term $`\stackrel{}{a}_\gamma ^{}\stackrel{}{\lambda }`$ must vanish separately. Furthermore, since the $`\lambda _i`$ are independent this can only be satisfied if $`\stackrel{}{a}_\gamma ^{}=\stackrel{}{0}`$ for each $`\gamma `$. This means that

$$\stackrel{}{a}_\gamma ^{}=\left(\frac{1}{|𝒢|}\underset{k}{}R^{(k)}\right)\stackrel{}{a}_\gamma =\stackrel{}{0},$$
(14)
i.e., the sum of all rotations applied to the original coordinate vector $`\stackrel{}{a}_\gamma `$ must vanish.

Similarly, to obtain a modified evolution corresponding to a target Hamiltonian $`H_{eff}^t=_\gamma (\stackrel{}{a}_\gamma ^t\stackrel{}{\lambda })B_\gamma `$, we require the weighted sum of rotations applied to the original coordinate vector to be equal to the corresponding target coordinate vector $`\stackrel{}{a}_\gamma ^t`$. I.e., for $`H_{eff}0`$, the following condition should be satisfied to obtain the desired evolution:

$$\stackrel{}{a}_\gamma ^{}=\stackrel{}{a}_\gamma ^t$$
(15)
This may require a combination of switching strategies for the BB pulses L. Viola, E. Knill, and S. Lloyd (2000).

It should be noted that our geometrical picture is an explicit representation of a subset of the group algebra $`𝒢`$ using the set of traceless Hermitian matrices and the identity as the basis. When the coefficients of the adjoint vector are real, the resulting matrix $`H_{eff}`$ is Hermitian. When they are complex, the resulting matrix is not Hermitian and the evolution is not unitary but may still be treated empirically M.S. Byrd and D.A. Lidar (2001).

## IV Errors

The picture developed above also gives an intuitive way in which to evaluate the error that remains in the system evolution after the application of the BB pulses. Let $`\stackrel{}{a}_\gamma ^t`$ be the coordinates vector corresponding to the desired Hamiltonian evolution and $`\stackrel{}{a}_\gamma ^{}`$ the actual vector after BB operations. Then $`\stackrel{}{a}_\gamma ^{}`$ corresponds to the effective Hamiltonian, Eq. (13) (and may be determined using quantum process tomography, see M.S. Byrd and D.A. Lidar (2001) and references therein). The error vector $`\stackrel{}{e}`$ is given by their difference in the $`n^2`$-dimensional vector space where our geometric picture holds:

$$\stackrel{}{e}=\stackrel{}{a}_\gamma ^{}\stackrel{}{a}_\gamma ^t$$
(16)
The vector $`\stackrel{}{e}`$ gives us the magnitude and direction of the error (i.e., the basis elements $`\lambda _i`$ give the type of error, e.g., bit-flip and/or phase-flip). Now consider the magnitude of this error,

$$d(\stackrel{}{a}_\gamma ^{},\stackrel{}{a}_\gamma ^t)=(\stackrel{}{e}^{}\stackrel{}{e})^{1/2}$$
(17)
(in the case of Hamiltonian evolution there is no need for complex conjugation). This is the Euclidean distance between the two vectors in the adjoint representation space. For two two-state density matrices, it is proportional to the Euclidean distance between the two Bloch vectors, as is the trace distance. In general, computing $`d(\stackrel{}{a}_\gamma ^{},\stackrel{}{a}_\gamma ^t)`$ is more manageable than other measures of distance (e.g., fidelity M.A. Nielsen and I.L. Chuang (2000)), since it does not require diagonalization.

In the case of imperfect BB operations, the goal is to minimize the distance $`d`$. For the purposes of optimization, note that

$$d(\stackrel{}{a}_\gamma ^{},\stackrel{}{a}_\gamma ^t)=(\stackrel{}{e}\stackrel{}{e})^{1/2}=((a_\gamma ^{})^2+(a_\gamma ^t)^22M\stackrel{}{a}_\gamma \stackrel{}{a}_\gamma ^t)^{1/2},$$
(18)
whereas the ordinary trace distance (Hilbert-Schmidt norm) between $`(\stackrel{}{a}_\gamma ^{},\stackrel{}{a}_\gamma ^t)`$ gives \[using Eq. (6)\]

$$\text{Tr}[(\stackrel{}{a}_\gamma \stackrel{}{\lambda })(\stackrel{}{a}_\gamma ^t\stackrel{}{\lambda })]=M\stackrel{}{a}_\gamma \stackrel{}{a}_\gamma ^t.$$
(19)
So minimizing $`d`$ is equivalent to maximizing $`(\stackrel{}{a}_\gamma ^{}\stackrel{}{a}_\gamma ^t)`$. The advantage of using $`d`$ is an intuitive one since the error vector simply describes a Euclidean vector (in the adjoint representation space).

For obtaining a desired unitary evolution, note that trace-norm distance for matrices, $`U`$ and $`V`$ is defined by

$$d_u(U,V)=\sqrt{1(1/n)\text{Re}[\text{Tr}(U^{}V)]},$$
(20)
where $`U`$, and $`V`$ are $`n\times n`$ matrices. For BB controls a short-time approximation is relevant. For the case of unitary evolution, the two measures, Euclidean distance and trace-norm for matrices, are equivalent. Approximating $`U`$ (desired evolution) and $`V`$ (actual evolution) by $`1\mathrm{l}iHt=1\mathrm{l}it_ia_{i\gamma }^t\lambda _i`$ and $`1\mathrm{l}iH^{}t=1\mathrm{l}it_ia_{i\gamma _i}^{}\lambda _i`$ respectively,

$`d(U,V)`$ $``$ $`\sqrt{1(1/n)\text{Re}[\text{Tr}(1\mathrm{l}+HH^{}t^2)]}`$ (21)
$``$ $`\sqrt{\stackrel{}{a}_\gamma ^{}\stackrel{}{a}_\gamma ^t}`$
is an $`O(t)`$ approximation to the unitary evolution. This is equivalent to Eq. (17).

## V Examples

We now discuss the example of storing a single qubit. In this case our geometrical picture can be cast in the familiar Bloch sphere representation.

Consider the noisy evolution of a stored qubit. Suppose the evolution of the qubit is governed by the Hamiltonian

$$H=\underset{\gamma }{}\underset{i=1}{\overset{3}{}}a_{i\gamma }\sigma _iB_\gamma ,$$
(22)
where the $`a_{i\gamma }`$ are real coefficients, $`B_\gamma `$ are bath operators, and the $`\sigma _i`$ are the Pauli matrices with the usual identification $`\sigma _1=\sigma _x`$, $`\sigma _2=\sigma _y`$, $`\sigma _3=\sigma _z`$. (As above, the identity component is neglected.) For faithful storage, a set of BB operations $`\{U_k\}`$ should serve to eliminate this Hamiltonian. Under such controls, the evolution is described by

$$H_{eff}=\frac{1}{|𝒢|}\underset{ik\gamma }{}U_k^{}[(\stackrel{}{a}_\gamma \stackrel{}{\sigma })B_\gamma ]U_k=\underset{\gamma }{}(\stackrel{}{a}_\gamma ^{}\stackrel{}{\sigma })B_\gamma ,$$
(23)
where $`U_kSU(2)`$ and $`R^{(k)}SO(3)`$, and

$$\stackrel{}{a}_\gamma ^{}=\frac{1}{|𝒢|}\underset{k=0}{\overset{|𝒢|1}{}}R^{(k)}\stackrel{}{a}_\gamma .$$
(24)
In the subsections below we consider different choices of the subgroup $`𝒢`$. These equations then describe a sum of vectors on the Bloch sphere. The mapping from the unitary matrices $`U_k`$ to the rotation matrices $`R^{(k)}`$ is given by

$`U^{}\sigma _iU`$ $`=`$ $`e^{i\sigma _3\alpha /2}e^{i\sigma _2\beta /2}e^{i\sigma _3\gamma /2}\sigma _ie^{i\sigma _3\gamma /2}e^{i\sigma _2\beta /2}e^{i\sigma _3\alpha /2}`$ (25)
$`=`$ $`R_{ij}\sigma _j.`$
Explicitly, the rotation matrix is given by

$$R=\left[\begin{array}{ccc}\mathrm{cos}(\alpha )\mathrm{cos}(\beta )\mathrm{cos}(\gamma )& \mathrm{sin}(\alpha )\mathrm{cos}(\beta )\mathrm{cos}(\gamma )& \mathrm{sin}(\beta )\mathrm{cos}(\gamma )\\ \mathrm{sin}(\alpha )\mathrm{sin}(\gamma )& \mathrm{cos}(\alpha )\mathrm{sin}(\gamma )& \\ & & \\ \mathrm{cos}(\alpha )\mathrm{cos}(\beta )\mathrm{sin}(\gamma )& \mathrm{sin}(\alpha )\mathrm{cos}(\beta )\mathrm{sin}(\gamma )& \mathrm{sin}(\beta )\mathrm{sin}(\gamma )\\ +\mathrm{sin}(\alpha )\mathrm{cos}(\gamma )& +\mathrm{cos}(\alpha )\mathrm{cos}(\gamma )& \\ & & \\ \mathrm{cos}(\alpha )\mathrm{sin}(\beta )& \mathrm{sin}(\alpha )\mathrm{sin}(\beta )& \mathrm{cos}(\beta )\end{array}\right].$$
(26)
Alternatively, one may use

$$\mathrm{exp}(i(\theta /2)\widehat{n}\stackrel{}{\sigma })\stackrel{}{x}\stackrel{}{\sigma }\mathrm{exp}(i(\theta /2)\widehat{n}\stackrel{}{\sigma })=\stackrel{}{x}^{}\stackrel{}{\sigma },$$
(27)
where

$$\stackrel{}{x}^{}=(\widehat{n}\stackrel{}{x})\widehat{n}+[(\widehat{n}\times \stackrel{}{x})\times \widehat{n}]\mathrm{cos}(\theta )+[\widehat{n}\times \stackrel{}{x}]\mathrm{sin}(\theta ).$$
(28)
The correspondence between the unitary and orthogonal groups is made by

$`{\displaystyle \underset{j}{}}x_j^{}\sigma _j`$ $`=`$ $`{\displaystyle \underset{i,j}{}}x_i[R_{\widehat{n}}(\theta )]_{ij}\sigma _j`$ (29)
$`=`$ $`U(\widehat{n},\theta /2)\left({\displaystyle \underset{i}{}}x_i\sigma _i\right)U^{}(\widehat{n},\theta /2),`$
where $`R_{\widehat{n}}(\theta )`$ is a rotation by $`\theta `$ about the axis $`\widehat{n}`$. Although this notation is more compact, the Euler angle parameterization of $`SU(3)`$ and $`SU(4)`$ have been given M. S. Byrd (1998, 2000); Todd Tilma, et al. .

### V.1 Storing One Qubit

To be specific, consider an unwanted pure dephasing interaction described by the Hamiltonian

$$H=g\sigma _3B.$$
(30)
Using Eq. (8) we find that the coordinate (adjoint) representation of this Hamiltonian is

$$a_i=\frac{g}{2}\text{Tr}(\sigma _i\sigma _3)=g\delta _{i3},$$
(31)
i.e., $`\stackrel{}{a}=(0,0,g)`$. This dephasing could be corrected through the use of a single BB operation $`U_1=\mathrm{exp}(i\sigma _1\pi /2)=i\sigma _1`$, so that

$$U_1^{}\sigma _3U_1=\sigma _3.$$
(32)
For the geometric picture, using $`\beta =\pi ,\alpha =\gamma =0`$ in Eq. (26):

$$R^{(1)}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),$$
(33)
which inverts the adjoint vector $`\stackrel{}{a}`$. This is shown schematically in Fig. 1. It is simple to check that, as required, $`H_{eff}=0`$. This example uses the lowest dimensional finite order group $`C_2`$.

Now, let us use our geometric picture to derive another class of BB operations for pure dephasing on a single qubit. Clearly, the point is to find a set of rotations of $`\stackrel{}{a}`$ which when added sum up to zero. The next example is the group $`C_3`$, which consists of rotating $`\stackrel{}{a}`$ by $`2\pi /3`$ and $`2\pi /3`$ about a fixed axis, i.e., uses two non-trivial BB operations. This is depicted in Fig. 2 where we have chosen $`\sigma _1`$ as the fixed axis. The set of rotation matrices that accomplish this are $`R^{(1)}=R_{\widehat{x}}(2\pi /3)`$ and $`R^{(2)}=R_{\widehat{x}}(2\pi /3)=(R^{(1)})^2`$, where

$`R_{\widehat{x}}(\theta )`$ $`=`$ $`\left(\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}(\theta )& \mathrm{sin}(\theta )\\ 0& \mathrm{sin}(\theta )& \mathrm{cos}(\theta )\end{array}\right)`$ (37)
and the corresponding unitary transformations are $`U_k=\mathrm{exp}(i\sigma _1(\pm \pi /3))`$, $`k=1,2`$.

A large number of BB operations is undesirable due to time constraints so that large sets become successively more difficult to implement effectively. However, depending on the symmetries of the Hamiltonian and the experimentally available BB operations, using a larger number may be advantageous.

Let us consider the next higher order set of BB operations. This will be a set with 4 elements. The subgroup condition requires a set forming either the cyclic group of order 4 or the so-called Vierergruppe since these are the only two groups of order 4 Michael Tinkham (1964). An example of the cyclic group of order 4 would be the four-fold rotations about a single axis, (e.g., $`\pi /2`$ around the $`\sigma _1`$ axis, as in Fig. 3). An example of the the other fourth order group is the set of rotations by $`\pi `$ about three orthogonal symmetry axes.

Note that the set of vectors pointing to the vertices of a tetrahedron also will sum to zero and thus form a set of adjoint vectors, representing BB modified Hamiltonians, that will produce the desired decoupling effect, the elimination of the interaction Hamiltonian Eq. (30). This set is determined by

$$\underset{k=1}{\overset{4}{}}\stackrel{}{a}^{(k)}=0,\text{and}\stackrel{}{a}^{(k)}\stackrel{}{a}^{(k^{})}=\mathrm{const}\mathrm{cos}(\theta ).$$
(38)
This implies $`\theta =\mathrm{cos}^1(1/3)`$ so that for Eq. (30) the set of rotations can be $`1\mathrm{l},R_{\widehat{y}}(\theta ),R_{\widehat{a}_2}(2\pi /3),R_{\widehat{a}_2}(2\pi /3)`$ acting on the initial vector $`\stackrel{}{a}_1=(0,0,g)`$, and $`\widehat{a}_2`$ is the direction of $`\stackrel{}{a}_2=R_{\widehat{y}}(\theta )\stackrel{}{a}_1`$. These rotations will take $`\stackrel{}{a}`$ to different positions which correspond to the vertices of a tetrahedron (see Fig. 4). The corresponding $`U_k`$ are found using Eq. (29).

Note that this last example uses a set of rotations that does not satisfy the subgroup condition. This shows that the subgroup condition is sufficient, but not necessary. Though it is not a necessary condition, it is important due to its convenience. The conditions both necessary and sufficient for first order decoupling are that the sum of the modified Hamiltonians, defined by the modified adjoint vectors, sum to the desired effective Hamiltonian.

### V.2 Two Qubits

Now let us suppose a computation is to be performed on two-qubits using the Heisenberg exchange coupling in the presence of a collective dephasing mechanism. From the study of computation on decoherence free subspaces (DFS), we know that this is possible since the Heisenberg interaction commutes with the group elements that form the stabilizer of the DFS J. Kempe, D. Bacon, D.A. Lidar, and K.B. Whaley (2001). Our goal here will be to interpret this condition geometrically.

Let a basis for the Lie algebra be given by

$$\{\lambda _i\}_{i=0}^{15}=\{\sigma _i\sigma _j\}_{i,j=0}^3,$$
(39)
and labels correspond to

$$\begin{array}{ccc}\lambda _j,j=0,1,2,3\hfill & & \sigma _i1\mathrm{l},i=0,1,2,3,\hfill \\ \lambda _j,j=4,5,6\hfill & & 1\mathrm{l}\sigma _i,i=1,2,3,\hfill \\ \lambda _j,j=7,8,9\hfill & & \sigma _1\sigma _i,i=1,2,3,\hfill \\ \lambda _j,j=10,11,12\hfill & & \sigma _2\sigma _i,i=1,2,3,\hfill \\ \lambda _j,j=13,14,15\hfill & & \sigma _3\sigma _i,i=1,2,3.\hfill \end{array}$$
(40)
This forms an orthogonal basis with respect to the trace and has normalization given by

$$\text{Tr}(\lambda _i\lambda _j)=4\delta _{ij}.$$
(41)
The Heisenberg interaction can be written as:

$$H_{ex}=J\stackrel{}{\sigma }_1\stackrel{}{\sigma }_2\stackrel{}{v}_1\stackrel{}{\lambda },$$
(42)
where $`\stackrel{}{v}_1=(0,0,0,0,0,0,J,0,0,0,J,0,0,0,J)`$, so that

$`H_{ex}=J(\lambda _7+\lambda _{11}+\lambda _{15})`$ (43)
The collective dephasing is given by:

$$H_I=g(\sigma _31\mathrm{l}+1\mathrm{l}\sigma _3)B\stackrel{}{v}_2\stackrel{}{\lambda }B.$$

where $`\stackrel{}{v}_2=(0,0,g,0,0,g,0,0,0,0,0,0,0,0,0)`$. So

$$H_I=g(\lambda _3+\lambda _6)B.$$
(44)
Note that $`\mathrm{Tr}(H_{ex}H_I)=4\stackrel{}{v}_1\stackrel{}{v}_2=0`$ so $`\stackrel{}{v}_1\stackrel{}{v}_2`$. A method for achieving the effective decoupling without the loss of the desired Heisenberg interaction is to consider the little group of $`\stackrel{}{v}_1`$. (The set of rotations that leaves this vector fixed.) From that set of rotations, a subset of rotations exists which will rotate the interaction Hamiltonian since the two vectors lie in orthogonal subspaces. These rotations clearly must be about the axis defined by $`\stackrel{}{v}_1`$. Thus we may express this as $`R_{\widehat{v}_1}(\theta )`$. To limit the number of pulses in the sequence of BB operations, a parity-kick operation is desired. This further limits our choices to those operations that rotate $`\stackrel{}{v}_2`$ by an angle $`\pi `$. More specifically, we seek a rotation that inverts the components in the directions $`\lambda _3`$ and $`\lambda _6`$, since $`\stackrel{}{v}_2=g(\lambda _3+\lambda _6)`$. The directions $`\lambda _3`$ and $`\lambda _6`$ define a plane perpendicular to $`\widehat{v}_1`$, so that the desired rotation matrix is effectively an $`SO(3)`$ rotation matrix with a non-trivial component in this plane, i.e., $`R_{\widehat{v}_1}(\pi )`$. It is then simple to check that the corresponding unitary operation satisfying the parity kick condition $`U^{}H_IU=R_{\widehat{v}_1}(\pi )H_I=H_I`$, is (see Fig. 5)

$$UU_1U_2=\mathrm{exp}(i(\sigma _1^{(1)}+\sigma _1^{(2)})\pi /2)=\sigma _1^{(1)}\sigma _1^{(2)},$$
(45)
where the superscript indicates the qubit on which the operator acts. This interaction leaves $`H_{ex}`$ unaffected and provides decoupling equivalent to Eq. (23). Note that this is a useful means for achieving the desired decoupling, because exchange interactions can be turned on during a gate operation in a solid state device and the decoupling can be achieved during the process without interruption of the desired interaction. The geometric picture shows that the above $`U`$ is by no means unique: any discrete $`SO(3)`$ subgroup acting in the $`(\lambda _3,\lambda _6)`$ plane, and whose elements add up to zero, will do.

## VI Conclusion

A geometric treatment of bang-bang (BB) operations has been provided. This perspective provides an intuitive picture for BB operations and their imperfections. The group averaging is made explicit through the corresponding average over a set of coordinate vectors representing rotations of the Hamiltonian; the resultant vector is the sum over all the BB modified Hamiltonians. These quantities are useful for computations, complementing the somewhat more abstract approaches of previous treatments L. Viola, E. Knill, and S. Lloyd (2000); Zanardi (2001). Since after the application of an imperfect set of decoupling operations, one is concerned with remaining error(s), such tools are useful for visualization. The often promoted Bloch sphere representation used in some of the examples treated here provides a means for extending intuition beyond the the low dimensional cases.

The usual group-theoretic symmetrization description of BB operations assumes that the set of pulses forms a discrete subgroup L. Viola, E. Knill, and S. Lloyd (2000); Zanardi (2001). We showed here that this is not a necessary condition, through the example of symmetrizing by the vertices of a tetrahedron. This fact has been well appreciated in the context of recoupling schemes in NMR (see for example R.R. Ernst, G. Bodenhausen and A. Wokaun (1987)). Here we wish to emphasize it in the context of quantum information processing.

The two-qubit example in Section V.2 provides a way in which this geometric analysis aids in the problem of finding decoupling interactions. The similarity between this example and recoupling techniques in NMR and other systems is no coincidence. The BB operations were, after all, related to NMR techniques in the earliest papers describing such interactions for quantum error correction Viola and Lloyd (1998). The geometric viewpoint is quite general and provides an instructive way in which to decompose such problems. They may be particularly useful for the types of recoupling techniques one requires for reducing constraints on quantum computer proposals D.A. Lidar and L.-A. Wu (2001).

In the subgroup framework, our geometric picture uses a homomorphic mapping between the fundamental representation and the corresponding adjoint representation. The problem of inverting this map from the adjoint to the fundamental representation may well be difficult for groups of higher dimension than $`SU(4)`$. However, for universal quantum computation, one and two qubit gates are sufficient and fortunately the discrete subgroups of unitary groups have been classified up to $`SU(4)`$. (See Amihay Hanany and Yang-Hui He (2001); W.M. Fairbanks, T. Fulton and W.H. Klink (1964); D. Anselmi, M. Bill , P. Fr , L. Girardello and A. Zaffaroni (1994) and references therein). The determination of the appropriate subgroup could consist of searching a discrete solution space. This appears feasible since the lower order subgroups are more relevant given the strict time constraints of the BB assumptions.

Given the scarcity of qubits in current quantum computing systems, we believe that the BB method is an important tool. We hope that the work presented here will be helpful in constructing sequences of BB pulses and analyzing their imperfections.

###### Acknowledgements.

DAL acknowledges support from PREA, PRO, the Connaught Fund, and AFOSR (F49620-01-1-0468). MSB acknowledges S. Schneider and K. Khodjasteh for help creating the graphics. We would also like to thank Lorenza Viola for several useful comments.