text:
# Recent Results on Particle Production from OPAL

## 1 Introduction

The large number of hadronic $`\mathrm{Z}^0`$ decays collected by the LEP experiments at $`\sqrt{s}=M_{\mathrm{Z}^0}`$ allow for detailed studies of the dynamics of the strong interaction. In three recent OPAL studies presented in this paper, production rates of various hadron species in the fragmentation of quark and gluon jets and of identified up, down or strange flavour jets have been studied. The results constitute a detailed test of models describing the fragmentation of partons in to hadrons in the final state.

## 2 $`𝝅^\mathrm{𝟎}`$, $`𝜼`$, $`𝐊^\mathrm{𝟎}`$ and charged multiplicities in quark and gluon jets

Because of the colour enhancement of the gluon-gluon coupling with respect to the quark-gluon coupling, the particle multiplicity in gluon jets is higher than in quark jets. As this multiplicity enhancement appears at the level of the couplings, it is expected to be largely independent of the particle species observed in the final state, with small corrections only due to e.g. the decay properties of heavy hadrons. In most fragmentation models, both for quark and gluon jet fragmentation, the neutralisation of colour fields occurs via the creation of quark-antiquark pairs. It has however been suggested  that in the fragmentation of a gluon jet, colour neutralisation could also occur via the creation of gluon pairs. In this case, an enhanced production of isoscalar mesons such as the $`\eta `$ would be expected in the fragmentation of gluon jets and, if they exist, also of glueballs.

OPAL has measured $`\pi ^0`$, $`\eta `$, $`\mathrm{K}^0`$ and charged particle multiplicities in quark and gluon jets as a function of the hardness scale of the jet defined as $`Q_{\mathrm{jet}}E_{\mathrm{jet}}\mathrm{sin}\frac{\theta }{2}`$. From hadronic $`\mathrm{Z}^0`$ decays with a 3-jet topology two samples of jets are selected consisting of the second and the third highest energy jet from each event, respectively. In the region $`8`$ GeV $`<Q_{\mathrm{jet}}<26`$ GeV, these samples contain different fractions of quark and gluon jets and the mean $`\pi ^0`$, $`\eta `$, $`\mathrm{K}^0`$ and charged particle multiplicities measured for the jets in these samples can be used to unfold the multiplicities for pure quark and gluon jets.

In Fig. 1 we show the mean charged particle, $`\pi ^0`$, $`\eta `$ and $`\mathrm{K}^0`$ multiplicities in quark and gluon jets as a function of $`Q_{\mathrm{jet}}`$. Also shown is the result of a fit to the charged particle multiplicities, which has been scaled by an appropriate normalisation factor in figures b,c and d.

Fig. 2 shows the ratio of the mean multiplicity in gluon and in quark jets for charged particles, $`\pi ^0`$, $`\eta `$ and $`\mathrm{K}^0`$. All cases can be described by the fit result obtained for charged particles. The multiplicity enhancement in gluon jets is thus found to be independent of the studied particle species. These results do not confirm an enhanced $`\eta `$ production in three-jet events reported by the L3 collaboration.

## 3 Leading particle production in light flavour jets

OPAL has measured leading $`\pi ^\pm `$, $`\mathrm{K}^\pm `$, $`\mathrm{K}_\mathrm{S}^0`$, $`\mathrm{p}(\overline{\mathrm{p}})`$ and $`\mathrm{\Lambda }(\overline{\mathrm{\Lambda }})`$ rates in up, down and strange flavour jets. As suggested in , these highest energy particles in a jet often carry the flavour of the quark from which the jet originated. So far few experimental results have confirmed this leading particle effect in the fragmentation of light flavour jets. The method used by OPAL to determine leading particle rates in light flavour jets was proposed in  and is based on counting all events in which a leading hadron is tagged in one of the event hemispheres and all events in which both hemispheres have a leading hadron tag. Using this information, our knowledge of the flavour composition of hadronic $`\mathrm{Z}^0`$ decays and imposing some constraints based on isospin symmetry and flavour independence, the production rates for the different hadron species in up, down and strange flavour jets can be unfolded.

In Fig. 3, the leading $`\pi ^\pm `$ and $`\mathrm{K}^\pm `$ rates in up, down and strange flavour jets are shown as a function of $`x_{\mathrm{cut}}`$, the lower threshold applied on the scaled momentum of the particle. The rates are compared to the corresponding JETSET and HERWIG predictions. In accordance with the expected leading particle pattern, the highest rates are those of pions in up and down flavour jets and kaons in strange flavour jets. The kaon rate in up flavour jets is lower than would naïvely be expected from the leading particle effect. This is due to the suppressed production of strange quark-antiquark pairs from the hadronic sea. The strangeness suppression factor, associated with this effect, can be determined directly from the ratio of leading $`\mathrm{K}^\pm `$ production in up and strange flavour jets and the ratio of leading $`\mathrm{K}^0`$ production in down and strange flavour jets, shown in Fig. 4a and 4b. The strangeness suppression factor obtained is

$$\gamma _\mathrm{s}=0.422\pm 0.049(\mathrm{stat}.)\pm 0.059(\mathrm{syst}.),$$

which is about one standard deviation high in comparison with other measurements, in most of which $`\gamma _\mathrm{s}`$ was determined rather indirectly, however.

We have also studied baryon production in the final state. In the Monte Carlo models to which we compare our data, baryon production is governed by the creation of diquark pairs when breaking up strings (in JETSET) or when decaying colourless clusters (in HERWIG). In Fig. 5, the ratio of proton over pion production in up flavour jets and of lambda over kaon production in strange flavour jets are compared to JETSET and HERWIG. The HERWIG Monte Carlo significantly overestimates the fraction of baryons produced.

## 4 Charged particle multiplicity in up, down and strange flavoured $`𝐙^\mathrm{𝟎}`$ decays

OPAL has studied mean charged particle multiplicities in up, down and strange flavoured $`\mathrm{Z}^0`$ decays. Because of the flavour independence of the strong interaction, these multiplicities are expected to be identical, with small corrections only due to the decay properties of heavy hadrons.

For the purpose of the measurement, the charged particle multiplicity has been determined for three event samples, defined by the observation of a leading $`\mathrm{K}^\pm `$ with fractional energy $`x_E>0.5`$, of a leading $`\mathrm{K}_\mathrm{S}^0`$ with $`x_E>0.4`$, or of a leading charged hadron with $`x_E>0.7`$, respectively. These samples differ in the relative fractions of up, down and strange flavoured events they contain and have been used to unfold the mean charged particle multiplicity for purely up, down or strange flavoured $`\mathrm{Z}^0`$ decays. These are found to be

$$n_\mathrm{u}=17.77\pm 0.51(\mathrm{stat}.){}_{1.20}{}^{+0.86}(\mathrm{syst}.),$$

$$n_\mathrm{d}=21.44\pm 0.63(\mathrm{stat}.){}_{1.17}{}^{+1.46}(\mathrm{syst}.),$$

$$n_\mathrm{s}=20.02\pm 0.13(\mathrm{stat}.){}_{0.37}{}^{+0.39}(\mathrm{syst}.),$$

where we point out that in particular the results for up and down flavoured events are strongly anti-correlated. Within the uncertainties, the measured mean multiplicities are consistent with being identical, as expected from the flavour independence of the strong interaction.

## 5 Conclusions

Mean $`\pi ^0`$, $`\eta `$, $`\mathrm{K}^0`$ and charged particle multiplicities have been determined for quark and gluon jets. The multiplicity enhancement in gluon jets is found to be independent of the studied particle species.

We have measured leading hadron production rates for various hadron species in up, down and strange flavour jets. The results confirm the leading particle effect in the fragmentation of light flavour jets and have been used for a direct determination of the strangeness suppression factor, yielding $`\gamma _\mathrm{s}=0.422\pm 0.049(\mathrm{stat}.)\pm 0.059(\mathrm{syst}.)`$. In addition the results constitute a detailed test of fragmentation models. In particular the HERWIG Monte Carlo model is found to strongly overestimate baryon production in the final state.

The mean charged particle multiplicity has been determined for up, down and strange flavoured $`\mathrm{Z}^0`$ decays. The results are identical within the uncertainties of the measurement, as expected from the flavour independence of the strong interaction.


format: text

text:
# Fast Numerical simulations of 2D turbulence using a dynamic model for Subgrid Motions

## 1 Introduction

The dynamics of high Reynolds number turbulent flow couples a large range of scales from the characteristic size of the domain to the dissipative scales. This range is usually too large to be fully resolved by Direct Numerical Simulation of the Navier-Stokes equations (DNS). This is the case for many applications like in aeronautics, geophysics or astrophysics where the typical Reynolds numbers are of the order of $`10^6`$ to $`10^{12}`$. The largest Reynolds numbers which can be achieved by DNS are of the order of $`10^4`$ to $`10^5`$ . The difficulty of direct simulation of turbulent flows at high Reynolds number arises because of two scaling laws: first, the memory requirement to be able to deal with all the scale of motions growths like $`Re^{9/4}`$, where $`Re`$ is the Reynolds number; second, the time stepping to be used to advance the equations has to be computed using the smallest resolved scale. As a result, computational time usually growths like $`Re^3`$. Despite of the growth of the computer power , the direct simulation of all realistic problems will not be conceivable before a long time.

Part of these difficulties could be circumvented if one could find an efficient approximate scheme to model the small scales of the turbulent flow, which usually monopolize 90 per cent of the computations and, for most applications, do not need to be known with the same precision as the large scales. From an algorithmic point of view, several methods have been proposed to try to describe the small turbulent scales with less degree of freedom than in a DNS like for example, sparse Fourier transform , rarefied modes , wavelets , self-adaptative grid mesh . From a theoretical point of view, the effort has been put on the model of the action of the small scales onto the larger scales of motion, or onto the small scales themselves, to be able to deal with simpler and cheaper systems to compute. In the extreme case of large eddy simulations (LES) or random averaged numerical simulations (RANS), the effect of the small scales at large scales is directly modeled as a function of the resolved or mean scales of motions, resulting in both memory and computational gain. This is often done at the price of introduction of arbitrary parameters, which require calibration and which lower the predictive power of the simulation. Also, some models are not theoretically satisfying, since they break for example the basic symmetries of the original equations .

In this paper, we continue a study initiated in  and develop a new dynamic model for the turbulent subgrid quantities which is directly derived from the constitutive equations of motions. This model respects all the basic symmetries of the original equations, as well as their conservation properties. The main assumption underlying the model is that the smaller turbulent scales are much less energetic than the resolved scales, and so, that their main interactions are with the resolved flow. It then appears reasonable to try to model these main interactions with as few approximations as possible, while more freedom is allowed regarding their mutual interactions. In the model we consider, the interactions between the resolved and subgrid scales are taken into account without any approximation, while the mutual interactions between subgrid scales are replaced by a turbulent viscosity. In 2D geometry, where energy condensation at large scale guarantees that the dynamics of the small scales is mainly non-local in scales this turbulent viscosity is so small that it can be set to zero in some simple cases . This allows a parameter free model of the subgrid scales, which will be described in the first part of this paper. In situations where dynamic processes at small scales, like vortex stretching, are important to regulate any small-scale instability, or even in 2D when ”large” and ”small” scales are not well enough separated, the turbulent viscosity cannot be set to zero . This option introduces a free parameter into the model, but, since it appears at a sub-dominant level, we may expect that its choice is not so critical to the success of the model.

After parameterization of the non-linear interactions between the subgrid scales, their resulting equation of motion becomes linear. This essential feature has two advantages: first, it enables, in certain simple flow geometries, analytical solutions for the subgrid-scale dynamics as a function of the resolved quantities, thereby closing the equations of motions at the resolved-scale level. This property was used to obtain analytical solution for mean profiles in channel flows  or for the Planetary Surface Layer . Second advantage of the linear description is that it allows for more efficient algorithms of integrations, using Lagrangian methods where the time stepping is done via a criterion based on the resolved scales. This allows numerical computations with a larger time step than DNS, and thus, at a lower computational cost.

An essential part of our model is the averaged Reynolds stresses describing the feedback of the subgrid scales onto the resolved component. In , we considered the scale separation parameter to be much less than the nonlinearity of the small scales and derived a feedback term in which the quadratic in the small-scale amplitude terms gave the main contribution into the Reynolds stress. Then, the model that consisted of two coupled “fluids” (the resolved and the subgrid one) was used to solve several problems with () and without scale separation (), the later being the typical problems such as the forced turbulence, the vortex merger and a turbulence decay. It was noticed that in the problems without a natural scale separation, such a model described very well the small-scale dynamics whereas improvements in modeling of the large scales (which is the main aim of the LES) was much more modest. The situation was clarified in  where direct numerical evaluations of the different contributions to the Reynolds stress were made. It was shown that in the problems without the scale separation the dominant contribution to the Reynolds stress comes from the linear rather than quadratic in the small-scale amplitude term.

In the present contribution, we focus on advancing the numerical approach introduced in  by introducing a better model for the turbulent Reynolds stress based on the results of the a priori tests performed in . Also, a more efficient procedure is used in this paper to generate the subgrid vorticity and velocity fields. Our goal here is to test our model both as a numerical method for an improved DNS, in the sense that we can compute the whole range of scales at a lower computational cost, and as pseudo LES method, in the sense that we compute only the larger scales, at an even lower computational cost. For the sake of simplicity, we shall consider here only the two-dimensional case, where both the hypothesis of the turbulent model  and its consequences  have been studied in detail. The generalization to the 3D case is the subject of an ongoing research.

## 2 The turbulence model

### 2.1 Scale decomposition

We consider a two-dimensional incompressible inviscid flow obeying the equations:

$$\{\begin{array}{ccc}_t\omega +\text{div}(𝐮\omega )\hfill & =\hfill & \nu \mathrm{\Delta }\omega ,\hfill \\ \text{div}𝐮\hfill & =\hfill & 0,\hfill \end{array}$$
(1)
where $`𝐮`$ is the velocity and $`\nu `$ is the viscosity. In 2D geometry, the vorticity has only one non-zero component that we denote by $`\omega `$. The resolved- and the subgrid-scale part of the velocity and the vorticity fields are defined via a filtering procedure:

$`𝐔(𝐱,t)`$ $`=`$ $`\overline{𝐮(𝐱,t)}={\displaystyle G(𝐱𝐱^{})𝐮(𝐱^{},t)𝑑𝐱^{}}`$ (2)
$`\mathrm{\Omega }(𝐱,t)`$ $`=`$ $`\overline{\omega (𝐱,t)}={\displaystyle G(𝐱𝐱^{})\omega (𝐱^{},t)𝑑𝐱^{}}`$ (3)
Here, $`G`$ is a filter, such that the resolved scales contain the main part of the total energy. In Section III.C, we shall propose a special shape of $`G`$ obeying this condition, and which is well adapted to our numerical method.

Each field is then decomposed as follows

$`𝐮(𝐱,t)`$ $`=`$ $`𝐔(𝐱,t)+𝐮^{}(𝐱,t),`$ (4)
$`\omega (𝐱,t)`$ $`=`$ $`\mathrm{\Omega }(𝐱,t)+\omega ^{}(𝐱,t),`$ (5)
where the primes denote subgrid-scale quantities. Inserting this decomposition into the Navier-Stokes equations and separating the resolved and the subgrid-scale parts, we get a set of coupled equations:

$`_t\mathrm{\Omega }`$ $`+`$ $`\text{div}(\overline{𝐔\mathrm{\Omega }})+\text{div}(\overline{𝐮^{}\mathrm{\Omega }})`$ (6)
$`+`$ $`\text{div}(\overline{𝐔\omega ^{}})+\text{div}(\overline{𝐮^{}\omega ^{}})=\nu \mathrm{\Delta }\mathrm{\Omega },`$
and

$`_t\omega ^{}`$ $`+`$ $`\text{div}(𝐔\mathrm{\Omega })\text{div}(\overline{𝐔\mathrm{\Omega }})`$ (7)
$`+`$ $`\text{div}(𝐔\omega ^{})\text{div}(\overline{𝐔\omega ^{}})`$
$`+`$ $`\text{div}(u^{}\mathrm{\Omega })\text{div}(\overline{u^{}\mathrm{\Omega }})`$
$`+`$ $`\text{div}(u^{}\omega ^{})\text{div}(\overline{u^{}\omega ^{}})=\nu \mathrm{\Delta }\omega ^{}.`$
The second up to the sixth terms of the l.h.s. of (7) are the contributions due to non-local interactions between the resolved and the subgrid scales. The last two terms in the l.h.s. of (7) are the contributions due to non-linear (local in scale space) interactions among the subgrid scales. In the favorable case where most of the energy is at the resolved scales, these last two terms can be expected to negligible with respect to, e.g., the second up to the fourth terms of the l.h.s. of (7). Indeed, in a recent detailed numerical analysis of the system (6) and (7), Laval et al  showed that for a very steep filter (when the filter $`G`$ is a cut-off in the spectral space), the small-scale dynamics (and thus the large-scale dynamics) is independent of the local interactions: when the later are neglected, the small-scale velocity and vorticity field are not significantly changed, in both forced and decaying turbulence, even after several eddy turn-over times . The analysis has also shown that in 2D turbulence, the leading order contribution at both large and small-scale comes from the correlations involving the large-scale velocity field, i.e. $`𝐔\mathrm{\Omega }`$ and $`𝐔\omega ^{}`$. The next leading order contribution comes from the correlations between small-scale velocity and vorticity $`u^{}\omega ^{}`$, while the correlation between small-scale velocity and large-scale vorticity $`u^{}\mathrm{\Omega }`$ gives the lowest contribution.

### 2.2 A priori numerical estimates

In the present case, the filter $`G`$ is smoother, and subgrid scales include both large and small scales. It is thus important to conduct additional numerical evaluation of the various terms of (6) and (7) using a smooth filter, to check the validity of the non-locality assumption for the subgrid scales. We performed an apriori test with the filter which will be used in our model. The choice of the filter will be discussed in section D. The test was done using a vorticity field from a DNS on a $`1024^2`$ grid. The Figures 1 and 2 give the distribution of each non-linear terms involved in the resolved scales and the subfilter-scales equations. The result of the filtering process one the initial field is shown fig. 3. Even with a smooth filter, the non-linear term involving only subgrid scales are still small compared to the higher order term.

### 2.3 The model

Keeping only the leading order contributions, the coupled equations (6) and (7) become:

$$_t\mathrm{\Omega }+\text{div}\overline{𝐔\mathrm{\Omega }}+\text{div}\overline{𝐔\omega ^{}}=\nu \mathrm{\Delta }\mathrm{\Omega },$$
(8)
$`_t\omega ^{}+\text{div}(𝐔\omega ^{})`$ $`=`$ $`F+\nu _t\mathrm{\Delta }\omega ^{},`$
$`F(𝐱,t)`$ $`=`$ $`(\text{div}(𝐔\mathrm{\Omega })\text{div}\overline{𝐔\mathrm{\Omega }})+\text{div}\overline{𝐔\omega ^{}}.`$ (9)
Here, $`\nu _t`$ is a turbulent viscosity, which will be introduced to damp the small-scale noise arising in our numerical scheme and $`F`$ is a force which describes the subgrid-scale generation via the enstrophy cascade (energy cascade in 3D).

## 3 Numerical implementation

### 3.1 Numerical strategy and performances

The resolved-scale equation has the form of an Euler equation with an additional forcing coming from the interaction with the subgrid-scale motions. One is therefore led to standard strategies to solve this equation, depending on the type of the flow (spectral methods for flows in simple periodic geometry, or finite difference or finite elements for more complicated geometries). Here, we shall consider the periodic case, thereby using a spectral code for solving the resolved-scale equation. This code is described in .

The situation is markedly different at the subgrid level, where the basic equation is linear in the subgrid motions, with an inhomogeneous part provided by the subgrid scale generation via the enstrophy cascade. This linearity suggests a solution strategy based on projection of the subgrid scales onto appropriate modes. Since the subgrid scales are usually very inhomogeneous (e.g., vorticity filaments), it seems logical to use a decomposition of the subgrid-scale velocity field into localized modes, thereby optimizing the memory requirements to store the subgrid-scale fields. A popular local mode decomposition uses wavelets (see e.g. ). Wavelets are however sometimes difficult to implement, and are not very handy to use in analytical computations. Here, we choose to use a Gabor decomposition, which provides a localized description while allowing theoretical manipulations similar to that obtained with Fourier modes. The linearity of the subgrid motions also provides room for a further computational time reduction via semi-Lagrangian methods of integration, using a time step related to the resolved scale. The numerical expected performances of our model, integrated using the semi-Lagrangian, Gabor method, at a given resolution and total integration time, is given in Table 1, and compared with both standard spectral methods based on FFT, and a LES approach based on FFT.

Computational time for several particular examples is given in table 2 which shows an obvious gain in speed for our method with respect to the DNS.

If we now compare with a classical LES, our method is somewhat slower because it is necessary to keep a sufficient number of localized modes for accurate description of the Reynolds stresses. However, a substantial gain in accuracy is achieved in comparison with traditional LES because of the better description of the nonlocal interaction of scales.

### 3.2 Description of the algorithm

Our algorithm of resolution of the turbulent model is based on five steps:

* Compute the force $`F`$ (9) using the resolved and subgrid fields at time $`t`$;
* Project this force onto a set of Gabor modes;
* Advance the subgrid field to time $`t+dt`$ in the Gabor space using this projection and the Lagrangian advection algorithm;
* Compute the resolved Reynolds stresses at time $`t+dt`$ using a procedure similar to an inverse Gabor transform;
* Advance the resolved velocity field at time $`t+dt`$ using this Reynolds stress.

Step 1 and 5 involve standard procedures linked with our resolved scale code (in the present case a Fourier spectral code with an Adams-Bashforth time stepping where a full desaliasing was introduced by keeping only the first 2/3 smaller wavenumbers in each direction, see ). Step 2, 3 and 4 involve new original procedures, based on interesting properties of the Gabor transform. For the sake of clarity , we present here only main results, leaving detailed computations in appendices.

### 3.3 The continuous Gabor transform

The Gabor transform (GT) is defined as

$`\widehat{𝐮^{}}(𝐱,𝐤,t)`$ $`=`$ $`{\displaystyle f(ϵ^{}(𝐱𝐱^{}))e^{i𝐤(𝐱𝐱^{})}𝐮^{}(𝐱^{})𝑑𝐱^{}},`$ (10)
where $`f`$ is a rapidly decreasing function at infinity. Note that $`1/ϵ^{}`$ has a meaning of the scale separation length. It has to be chosen to lie in between of the integral scale and the Kolmogorov scale, the later being far in the subgrid range in many applications. For LES purposes, $`1/ϵ^{}`$ has to be close to (but not less than) the minimal resolved scale (grid scale $`dh`$ in our case). Thus, to derive our model we need to assume the subfilter-scale wavenumber $`k`$ to be greater than $`ϵ^{}`$ and use $`ϵ^{}/k`$ as a small parameter. Technically, however, is convinient to keep $`k`$ fixed and perform expansions in small $`ϵ^{}`$.

In the special case where $`f=\sqrt{G}`$, where $`G`$ is the filter (see eq. (3)), we have the following simple reconstruction formulae:

$$𝐮^{}(𝐱,t)=\frac{1}{(2\pi )^2f(0)}\widehat{𝐮^{}}(𝐱,𝐤,t)𝑑𝐤,$$
(11)
$$\overline{𝐔\omega ^{}}(𝐱,t)=𝐔\frac{1}{(2\pi )^2}\mathrm{}\left[\widehat{\omega ^{}}(𝐱,𝐤,t)f(𝐤,t)\right]𝑑𝐤+O(ϵ^{}),$$
(12)
where $`\mathrm{}`$ means the real part, and finally

$$\overline{u_i^{}u_j^{}}(𝐱,t)=\frac{1}{(2\pi )^2}\mathrm{}\left[\widehat{u_i^{}}(𝐱,𝐤,t)\widehat{u_j^{}}(𝐱,𝐤,t)\right]𝑑𝐤+O(ϵ^{}).$$
(13)
One may also note the connection between the Gabor velocity and the Gabor vorticity (similar to the Fourier quantities):

$$\widehat{𝐮}(𝐤,t)=\frac{i}{k^2}\widehat{\omega }(𝐤,t)(e_z\times 𝐤)+O(ϵ^{}),$$
(14)
where $`e_z`$ is a unit vector in $`z`$-direction.

Finally, one may show that the GT of the subgrid-scale equation (2.3) is (see appendix A and  for more details)

$$\left(_t+𝐔_x_x\left(𝐔𝐤\right)_k\right)\widehat{\omega ^{}}(𝐱,𝐤,t)=\widehat{F}(𝐱,t)\nu _tk^2\omega ^{}(𝐱,𝐤,t).$$
(15)
### 3.4 The discretized optimum Gabor transform

The continuous Gabor transform is not suitable for numerical implementation. Moreover, because it involves both position and wavenumbers, it theoretically describes a given field via a much larger mode numbers than in traditional spectral methods (typically $`N^{2D}`$ versus $`N^D`$ in D dimensions). However, as discussed previously, inhomogeneous fields may be represented with a high precision by a much smaller set of Gabor modes than the theoretical number, a number even potentially smaller than $`N^D`$ (see e.g. ). The algorithm we devised uses $`N^D`$ Gabor modes and, therefore, there is a room for further code optimization.

#### 3.4.1 The discretization

The discretization of the subgrid field in the Gabor space is done via a Particle In Cell (PIC) method. This method is traditionally used in plasma physics, but was used previously in hydrodynamics by Nazarenko et al  to study the interaction of sound wave-packets with turbulence. Details about the method can also be found in . In this method, the Gabor modes of the subgrid-scale vorticity field are replaced with discrete wave-packets (particles). For example, the vorticity field $`\widehat{\omega ^{}}(𝐱,𝐤,t)`$ is discretized with $`N_p`$ wave-packets $`\alpha `$ in the Gabor space:

$$\widehat{\omega ^{}}(𝐱,𝐤,t)=\underset{\alpha =1}{\overset{N_p}{}}\widehat{\sigma }_\alpha (t)S_𝐱(𝐱𝐱_\alpha (t))S_𝐤(𝐤𝐤_\alpha (t)).$$
(16)
In (16), $`\alpha `$ labels the wave-packet. Note that our representation allows several wave-packets with different wavenumbers to be located at the same position $`𝐱`$, like in the Gabor transform. Here $`S_𝐱`$ and $`S_𝐤`$ are some interpolating functions (particle “size”). Since the computation of non-linear terms in the resolved-scale equation involves a second order spatial derivative, we are led to choose a linear interpolation in the $`𝐱`$ direction. In the $`𝐤`$ direction, a zero particle “thickness” turns out to be sufficient for our purpose. We thus adopt the following representation,

$`S_𝐱(𝐱)`$ $`=`$ $`S(x)S(y),`$ (17)
$`S_𝐤(𝐤)`$ $`=`$ $`\delta (p)\delta (q),`$ (18)
where $`𝐱=(x,y)`$, $`𝐤=(p,q)`$, $`\delta `$ is the Dirac function and the function $`S(\eta )`$ is defined by:

$`S(\eta )`$ $`=`$ $`\{\begin{array}{cc}(dh|\eta |)/dh\hfill & \mathrm{if}|\eta |<dh,\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}`$ (21)
Here, $`dh`$ is a length scale governing the accuracy of the discretization. Its choice will be discussed later on.

The PIC algorithm of reconstruction of the continuous fields from the discretized coefficients provides a natural filter $`G`$ and the function $`f`$ to be used in the Gabor transform through (see appendix D):

$$G(𝐱)=f^2(𝐱)=f^2(0)\times S_𝐱^2(𝐱),$$
(23)
with

$$f(0)=\frac{1}{\sqrt{S_𝐱^2(𝐱^{})𝑑𝐱^{}}}.$$
(24)
In the k space, the filtering operation is achieved by simply multiplication of the Fourier coefficients by the Fourier transform of $`G`$, which is

$$\widehat{G}(p,q)=\frac{36}{(pdh)^2(qdh)^2}\left\{\left(1\frac{sin(pdh)}{pdh}\right)\left(1\frac{sin(qdh)}{qdh}\right)\right\}$$
(25)
#### 3.4.2 Reconstruction of the subgrid scale correlations

With our discretization, the formulae of reconstruction of the subgrid correlations (12) can be obtained using (16). Taking into account the condition that $`\omega ^{}`$ is real, we have (see appendix C for details):

$$\overline{𝐔\omega ^{}}(𝐱,t)2𝐔\underset{\alpha _+=1}{\overset{N_p}{}}\mathrm{}\left[\widehat{\sigma _{\alpha _+}}f^{}(𝐤_{\alpha _+})\right]S_𝐱(𝐱𝐱_{\alpha _+}),$$
(26)
where the sum is only over the wave-packets with positive wavenumber $`p_\alpha `$<sup>1</sup><sup>1</sup>1. Since the vorticity $`\omega ^{}(𝐱,t)`$ is real, each wave-packet at position $`𝐱_\alpha `$ and with wavenumber $`𝐤_\alpha `$ will have a partner wave-packet with same amplitude, at the same location and with an opposite wavenumber $`𝐤_\alpha `$.. Similarly, one can use (47) and (16) to re-write (13) as (65) (see appendix C)

$`\overline{u^{}v^{}}(𝐱,t)=2{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}{\displaystyle \frac{q_{\alpha _+}p_{\alpha _+}}{(p_{\alpha _+}^{}{}_{}{}^{2}+q_{\alpha _+}^{}{}_{}{}^{2})^2}}|\widehat{\sigma }_{\alpha _+}|^2S_𝐱^2(𝐱𝐱_{\alpha _+})`$
$`\overline{u_{}^{}{}_{}{}^{2}}(𝐱,t)=2{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}{\displaystyle \frac{+q_{\alpha _+}^2}{(p_{\alpha _+}^{}{}_{}{}^{2}+q_{\alpha _+}^{}{}_{}{}^{2})^2}}|\widehat{\sigma }_{\alpha _+}|^2S_𝐱^2(𝐱𝐱_{\alpha _+})`$ (27)
$`\overline{v_{}^{}{}_{}{}^{2}}(𝐱,t)=2{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}{\displaystyle \frac{+p_{\alpha _+}^2}{(p_{\alpha _+}^{}{}_{}{}^{2}+q_{\alpha _+}^{}{}_{}{}^{2})^2}}|\widehat{\sigma }_{\alpha _+}|^2S_𝐱^2(𝐱𝐱_{\alpha _+})`$
#### 3.4.3 Optimum discrete Gabor transform

For a given field $`\omega ^{}(𝐱,t)`$, associated with a grid of size $`2\pi /N`$, our ”optimum” discrete Gabor transform only retains one wavenumber per position (i.e. $`N^2`$ Gabor modes, for a $`N^2`$ initial discrete field). This wavenumber, and the amplitude of the corresponding Gabor mode are chosen as follows: we use the exact relation between the field and its $`N^2`$ wave-packet components

$$\omega ^{}(𝐱,t)=\frac{2}{f(0)}\underset{\alpha _+=1}{\overset{N_p}{}}\mathrm{}\left[\sigma _{\alpha _+}(t)\right]S_𝐱(𝐱𝐱_\alpha ).$$
(28)
An easy way to satisfy (28) exactly is to create one wave-packet at each grid point, so that $`N_p=N^2`$. Applying equation (28) at each grid point $`𝐱_i`$ we then obtain:

$$\omega ^{}(𝐱_i,t)=\frac{2}{f(0)}\mathrm{}\left[\sigma _{\alpha _i}(t)\right].$$
(29)
This equation fixes the real part and the physical coordinates of the wave-packet. To find its wavenumber, we use the velocity at the grid point $`u^{}(𝐱_i)`$ to find two conditions:

$`u^{}(𝐱_i,t)`$ $`=`$ $`{\displaystyle \frac{2}{f(0)}}{\displaystyle \frac{q_{\alpha _i}}{(p_{\alpha _i}^2+q_{\alpha _i}^2)}}\mathrm{}\left[\sigma _{\alpha _i}(t)\right],`$ (30)
$`v^{}(𝐱_i,t)`$ $`=`$ $`{\displaystyle \frac{2}{f(0)}}{\displaystyle \frac{+p_{\alpha _i}}{(p_{\alpha _i}^2+q_{\alpha _i}^2)}}\mathrm{}\left[\sigma _{\alpha _i}(t)\right].`$ (31)
These conditions link the two components of the wavenumber with the imaginary part of the amplitude of the wave-packet, which is still a free parameter at this stage. We then select the phase of the GT by requiring that the imaginary part of the wave packet equals its real part. Other choices could have been made, but this one turns out to simplify the computations. The characteristics of each wave-packets $`\alpha _i`$ can finally be summarized as:

$$\{\begin{array}{ccc}\mathrm{}\left[\sigma _{\alpha _i}\right]\hfill & =\hfill & f(0)\omega _f^{}(𝐱_i)/2,\hfill \\ \mathrm{}\left[\sigma _{\alpha _i}\right]\hfill & =\hfill & \mathrm{}\left[\sigma _{\alpha _i}\right],\hfill \\ q_{\alpha _i}/p_{\alpha _i}\hfill & =\hfill & u_f^{}(𝐱_i)/v_f^{}(𝐱_i),\hfill \\ p_i\hfill & =\hfill & \frac{\omega _f^{}(𝐱_i)\left(q_{\alpha _i}/p_{\alpha _i}\right)^2}{v_f^{}(𝐱_i)(1+\left(q_{\alpha _i}/p_{\alpha _i}\right)^2)},\hfill \\ 𝐱_{\alpha _i}\hfill & =\hfill & 𝐱_i.\hfill \end{array}$$
(32)
This procedure creates $`N^2`$ wave-packets, from any vorticity field on a grid of size $`2\pi /N`$ and can be used to initialize the subgrid scale vorticity field from any given initial condition or to GT the force $`F`$ by transforming it into an equivalent vorticity field $`\omega _f^{}(𝐱,t)=F(𝐱,t)/dt`$.

### 3.5 The Lagrangian scheme

A Lagrangian interpretation of the subgrid scale equation (15) shows that its time integration is equivalent to the evolution of each wave-packet in the (x,k) space. Each wave-packet carries a complex amount of Gabor vorticity $`(\sigma _\alpha )`$ and is advected at the resolved-scale velocity $`𝐔(𝐱)`$, while its wavenumber and amplitude evolve according to the local resolved strain. For the trajectory of the wave-packet and its amplitude we have

$`\dot{𝐱}_\alpha `$ $`=`$ $`𝐔(𝐱_\alpha (t)),`$ (33)
$`\dot{𝐤}_\alpha `$ $`=`$ $`_x\left(𝐤_\alpha (t)𝐔(𝐱_\alpha (t))\right),`$ (34)
$`\dot{\sigma }_\alpha `$ $`=`$ $`\widehat{F}(𝐱_\alpha ,𝐤_\alpha ,t)\nu _tk_\alpha ^2\sigma _\alpha .`$ (35)
By these equations, the wave-packets evolve continuously in the physical space in between the grid points associated with the finite, resolved-scale resolution. To find the resolved-scale quantities at the position of each particle (i.e. possibly in between grid mesh points), an interpolation procedure is used via the functions $`S_𝐱`$ and $`S_𝐤`$ associated with the PIC method.

### 3.6 Initial conditions

All the simulations used the same initial random vorticity field with all the energy concentrated at very large scales ( the initial energy spectrum is given by $`E(k)=ke^{(kk_o)^2}`$ with $`k_o=\mathit{1}`$). In practise, the noise associated with the PIC method is a function of the number of wave-packets used for computation. If the initial condition is such that the initial grid on which the vorticity is defined is too coarse ($`M`$ too small), one can use a two step procedure to reach a reasonable number of wave-packets: in a first stage of the simulation, we create at each time step $`M^2`$ wave-packets with the vorticity $`\omega _f^{}(𝐱,t)=F(𝐱,t)/dt`$ created by the forcing $`F`$. During this stage, the wave-packets which were created at earlier times are moved into the phase space using the ray equations (33,34), but their amplitude is kept constant. This procedure is used until the total number of wave-packet reaches a desired number. From this point on, the wave-packet creation is shut down, and the wave-packets are evolved according to (33,34,35).

### 3.7 Noise reduction and effective turbulent viscosity

For inviscid simulations (when the turbulent viscosity is taken equal to $`0`$), the numerical procedure develops a noise at very subgrid scales, which is a function both of the integration time and the total number of particles. We have found that an efficient way to reduce this subgrid-scale noise could be achieved by periodically recreating a new set of $`N_p`$ wave-packets via first a rebuilding of the vorticity field in the physical space, on a grid of size $`\sqrt{N_p}`$ using (11), followed by a re-creation of wave-packets using (32). This procedure keeps the correlations associated with the subgrid-scale field unchanged, and therefore does not affect directly the resolved scale. In effect, this procedure acts as a filter for the subgrid scales which are smaller than the size of the reconstruction grid, and thus, can be seen as an effective turbulent viscosity, which cannot be estimated a priori but which adjusts itself to the noise level.

### 3.8 Test of the accuracy of the method

To test the accuracy of our numerical method, we have performed a priori tests using data from a DNS at high resolution. The total vorticity field was divided into a subgrid-and a resolved-scale field using the same filter than in our model. The subgrid-scale field is further discretized into wave-packets according to eq. (32). A comparison was then made between the various components of the stresses in the DNS, and in our discretization scheme. The result of is shown fig. 4, which shows the square modulus of each of the three Reynolds stress components with respect to their wavenumbers. The accuracy of the discretized reconstruction of each terms is remarkable for the large scales. At very small scales ($`k>42`$), the discretization tends to produce a noise in the components of the Reynolds stress involving the subgrid-scale velocity. In our case, we do not consider these terms so that we do not have to bother about this noise. In situations where one, or two of these troublesome terms are considered, the noise can be removed via a filtering at small scale. In spectral simulations, this truncature is provided naturally via the procedure to remove aliasing, which filters out modes larger than $`2k_{max}/3`$ where $`k_{max}`$ is the maximum wavenumber of the resolved-scale field.

The inset of figure 4 shows the comparison between the total real Reynolds stresses versus the total modeled Reynolds stress after desaliasing. The agreement is nearly perfect.

## 4 Numerical results

Our method can be used to perform two kind of simulations, depending on whether one is interested in the small-scale behavior or not. In the first case, one needs to model the subgrid scales with a large number of wave-packets, so as to be able to reconstruct the small-scale field with a good accuracy. Typically, one needs about $`N_p=N^2`$ wave-packets to be able to reconstruct faithfully details at scale $`2\pi /N`$ and to produce a result which may be compared with a Direct Numerical Simulation at resolution $`N^2`$. In the case when only large scales matter, one needs to keep the smallest number of wave-packet necessary to compute accurately the Reynolds stresses at the wavenumber cut-off. Since the resolved- scale field at the cut-off $`k_c=M/2`$ produces (by non-linear beating) information up to scale $`2\pi /M`$, we used a minimum of $`N_p=M^2`$ wave-packets in our ”Large Eddy Simulations”. Finally, note that our method allows nearly inviscid computations, since it does not require the existence of a viscosity at large scale and since it uses a minimal ”effective viscosity” which starts acting only at scales $`2\pi /M`$. Very large Reynolds number can then be achieved via an adequate number of wave-packets.

We present results illustrating these points in two classical situations, - decaying and forced turbulence. In each case, the performance of our model are discussed and compared with results from the DNS and other popular subgrid-scale parameterizations used in 2D turbulence.

### 4.1 Decaying turbulence

For the decaying case, the initial condition was chosen so that the energy is concentrated at very large scales. The reference DNS was performed at a resolution $`N^2=1024^2`$ with a viscosity $`\nu =\mathrm{1.8\hspace{0.33em}10}^4`$, leading to a Reynolds number $`Re10^4`$. The simulation was stopped after approximatively 50 turnover times, by which time the initial condition has evolved into a robust dipole structure. The separation between resolved and subgrid scale used for our model is taken at $`k_c=21`$, corresponding to a computation over a grid $`64^2`$.

Let us first compare the DNS to our model using a large number of wave-packets ($`N_p=512^2`$). We ran two different simulations: one in which the viscosity at resolved scale was set to $`\nu =\mathrm{1.8\hspace{0.33em}10}^4`$, like in the DNS (”viscous” simulation); another one in which $`\nu `$ was set to zero at resolved scale (”inviscid” simulation). The total simulation time of each simulation is about 1h30 on a Sun workstation, roughly 150 times less than the DNS.

The total vorticity field after 50 turn over times, reconstructed by adding the contribution from the resolved and subgrid scales, is shown in Fig. 5 and compared with the vorticity field of the DNS. A global comparison for small and large scales can also be done via energy spectra. This is done in Fig. 6.

Clearly, the two pictures display good overall similarities at large scales, and marked differences at smaller scales. The largest structures are very similar and they are well localized even after 50 turnover times. The spectra in the two simulations overlap. The viscous simulation gives large-scale structures which are in closer agreement with the (viscous) DNS (see Table 2). Together, these results confirm that at large scale, the dynamically important coupling term between resolved and subgrid scale is the term $`div(\overline{U}\omega ^{})`$, in agreement with the dynamical analysis of Laval et al , performed using a cut-off filter. At smaller scales, our model seems to produce thinner filaments and smaller structures. This effect is clearly visible on the spectra of the two simulations: in the DNS, the $`k^3`$ inertial law starts to level off towards $`k=60`$ (due to viscous effects), while in our model, the power law extends over a wider range of scales, up to approximately $`k=100`$. To check whether this difference comes from our subgrid scale scheme, we performed a direct comparison of the smallest scales of the simulation $`k>21`$. Fig. 7 shows such a comparison at some earlier time (after about 15 turnover times), when the small scales have not yet been washed out by viscosity.

At this time, differences are negligible, proving that our model also captures the dominant coupling mechanism at subgrid scales. This finding is also in agreement with the dynamical analysis of Laval et al . The further differences arising over longer time scales can therefore be due to two effects: one is the error accumulation due to sub-dominant neglected terms in our model (like the terms involving coupling of the resolved vorticity $`\mathrm{\Omega }`$ with the subgrid velocity $`𝐮^{}`$); the second is viscous effects. We believe that the first possibility is ruled out by the dynamical analysis of Laval et al , in which numerical simulations of our model equations were performed using the same viscosity than in the direct numerical simulations <sup>2</sup><sup>2</sup>2These simulations were not fast: they were based on spectral methods and were even more slowly than the direct method.. In that case, the differences at small scales appear to be negligible. We therefore interpret the differences between the two models as a viscous effect, and, actually, as an indication of the lower effective viscosity in our model. This would explain both the later bending of the spectra, and the finer structure of the filaments.

The second series of comparison were performed using the minimum number of wave-packet, equal in the present case to $`N_p=128^2`$. In such a case, the number of wave-packets does not allow accurate reconstruction for the scales less than $`2\pi /128`$. The simulation is however much faster, and takes only 8 minutes of computational time. The total large-scale vorticity field after 50 turn over times is compared with the corresponding large-scale vorticity field of the DNS in Fig. 8.

The agreement is still very good and to quantify this agreement, we computed the correlation coefficient between the two simulations, defined by:

$$C_{1,2}=\frac{_{1<i,j<64}\omega _1(i,j)\omega _2(i,j)}{\sqrt{_{1<i,j<64}\omega _1(i,j)^2_{1<i,j<64}\omega _1(i,j)^2}}$$
(36)
where $`\omega _1(i,j)`$ and $`\omega _2(i,j)`$ are the large-scale vorticity field. This correlation coefficient was computed at two different times, corresponding to 15 and 50 turnover time, and reported in table 2. At earlier time, our three models (viscous or not) are all characterized by a very good correlation coefficient (about 99 percent). At later time, a slight difference appear. In fact, the best correlations are achieved by the viscous model, and the model with the least number of subgrid-scale modes. This is not surprising, since our noise removing procedure produces an effective viscosity which is larger as the number of subgrid-scale modes decreases. The high resolution model M1 is then the less viscous model of all three, and therefore, can be expected to produces the largest difference with respect to the viscous DNS.

Two other popular 2D turbulent model were also tested along the same line. In the first one (HDNS), the viscous term $`\nu \mathrm{\Delta }\mathrm{\Omega }`$ of the Navier-Stokes equation is replaced by a ”hyper-viscous” term $`\nu _p\mathrm{\Delta }^p\mathrm{\Omega }`$. In our simulations, we took $`p=8`$ and $`\nu _p=10^{18}`$. The second model is the Anticipated Potential Vorticity Model (APVM) developed by Sadourny and Basdevant . Both models are very cheap, taking only about 2 minutes of computational time. They are also less accurate, as can be seen from both the Fig. 8, and the table. In the HDNS $`64^2`$, the correlation coefficient is only about 50 percent in the end of the simulation. For the APVM, it is higher (about 75 percent), but still lower than our ”minimal model”. The energy spectra of the 4 simulations can also be compared. This is done in Fig. 6. Our model develops an energy spectrum very close to the APVM one, and slightly less steep than the DNS at scales close to the cut-off. This is because at this scale, viscous effects start being felt and tend to bend the spectrum. The HDNS spectrum is much steeper near the cut-off than both the reference DNS and the two other models.

### 4.2 Forced turbulence

Similar simulations were performed in the case of forced turbulence. In such a case, our model was ran with $`M^2=64^2`$ resolved-scale Fourier modes, and $`N_p=512^2`$ subgrid-scale Gabor modes, and no additional viscosity. Run were also performed using a HDNS or the APVM model over $`64^2`$ Fourier modes. The initial condition is a vorticity field with an energy spectrum concentrated at the forced wavenumber (k=15) and with a small amount of total energy. The simulation was forced by keeping constant the energy of the mode k=(15,0). In this situation, the vorticity field is progressively built via the stochastic forcing, which is itself strongly dependent on the exact structure of the vorticity field (since it must aim at keeping one mode constant). Due to the developing spectral cascades, the scale interactions become less nonlocal and we may then expect sub-dominant terms to play an enhanced role (with respect to the decaying case). Indeed, we have observed that the vorticity fields in the DNS and in the model do not exactly correspond: both simulations display similar small scale intense vortices, but they are not located at the same places. This effect can also be seen more clearly in the spectra. They are shown in fig. 9.

Due to the chaoticity, they all differ at the mode $`k=1`$ which is the most sensitive to the exact position of each individual vortex. At smaller scales, marked differences appear between the models. Clearly, the HDNS gives the worst result, with a large deficit of energy over all scales. This is because in this case, the forcing scale is very close to the scale at which the dissipation takes place, and most of the energy is dissipated before the inverse cascade to larger scales can occur. The two other models, which do not introduce an explicit dissipation at the cut-off scale, perform much better. The APVM model tends to overestimate the rate of energy at the largest scale and to underestimate it in the inertial range. Our model slightly underestimates the amount of energy at the largest scale, but gives results very close to the DNS in the inertial range of scale. At scales smaller than the cut-off, both the DNS and our model produce an energy spectra steeper than the $`k^3`$ law. At the largest wavenumber, the beginning of the viscous range is clearly visible in the DNS but not in our model. This is again an indication of the large Reynolds number achieved by our model. As for computational times, the performances are similar to the decaying case.

## 5 Discussion

We developed a new dynamic model of subgrid-scale turbulence based on a simple hypothesis about the subgrid-scale evolution. This approach is different from traditional turbulent models since our model provides expressions of the turbulent Reynolds sub-grid stresses via estimates of the sub-grid velocities rather than velocities correlations. The subgrid-scale dynamic is given by a linear equation, describing the advection of subgrid-scale wave-packets with the mean flow. This feature allows a reduction of the time step used in the simulation, via the use of a pseudo-Lagrangian method. We thereby achieved a reduction of the computational time by a factor $`150`$ in the typical cases we considered. Our method can also be used for Large Eddy Simulation strategies, in which only the large scales are computed. This allows a large memory savings, and an additional reduction of the computational time (typically a factor 10 in the case we considered). The resulting simulation is more costly than a traditional LES simulation, based on Fast Fourier Transform algorithms (like hyperviscosity or vorticity dissipative schemes). It could however become more competitive for more complicated geometries, where finite difference schemes become more appropriate than FFT. This would make LES models based on hyperviscous schemes much more costly, while our method would proportionally keep at the same computational performance.

In the present paper, only 2D turbulence has been considered. As is well known, 2D turbulence is very special, because there is no vortex stretching. One may therefore wonder whether our model could be applied to 3D flows. To answer this question, Laval and Dubrulle  have recently performed a numerical analysis of the hypotheses pertaining our turbulent model. The main hypothesis is that the subgrid dynamics can be very well approximated by an equation in which non-local terms are retained, while local terms (non-linear in the subgrid velocities) are replaced by a turbulent viscosity. The numerical analysis, performed via a numerical simulation of the coupled resolved and approximate subgrid equations indeed showed that such hypothesis is valid, insofar as both the energy spectra, the structures and the statistical properties of such a model were very close to that of a DNS. Even so this analysis can only be performed at a rather moderate Reynolds number (80 to 200, based on the Taylor scale), we find this analysis encouraging.

## Appendix A Derivation of the subgrid scale equation in (x,k) space

Some properties of the GT will be useful in the sequel. They are (see  for details):

$`\widehat{_i𝐮^{}}(𝐱,𝐤,t)`$ $`=`$ $`_i\widehat{𝐮^{}}(𝐱,𝐤,t)`$ (37)
$`=`$ $`ik_i\widehat{𝐮^{}}(𝐱,𝐤,t)+O(ϵ^{})`$ (38)
$$\widehat{𝐔\omega ^{}}(𝐱,𝐤,t)𝐔(𝐱,t)\widehat{\omega ^{}}(𝐱,𝐤,t)+i\left(_x_k\right)𝐔\omega ^{}.$$
(39)
Let us derive the GT of the equation:

$$_t\omega ^{}+div(U_j\omega ^{})=F(𝐱,t)+\nu _t\mathrm{\Delta }\omega ^{}.$$
(40)
Using (38), we find that the GT of the viscous term is $`\nu _tk^2\widehat{\omega }^{}`$. The GT and the time derivative commute, so that the GT of the first term of the lhs gives :

$$\widehat{_t\omega ^{}}=_t\widehat{\omega ^{}}$$
(41)
Using the space derivative property (38), the GT of the second terms can be developed into:

$`\widehat{_jU_j\omega ^{}}`$ $`=`$ $`\widehat{U_j_j\omega ^{}}`$
$``$ $`U_j\widehat{_j\omega ^{}}+i_lU_j{\displaystyle \frac{}{k_l}}\widehat{_j\omega ^{}}`$
$``$ $`U_j_j\widehat{\omega ^{}}+i_lU_j{\displaystyle \frac{}{k_l}}\left(ik_j\widehat{\omega ^{}}\right)`$
$``$ $`U_j_j\widehat{\omega ^{}}+i_lU_jk_j{\displaystyle \frac{}{k_l}}\widehat{\omega ^{}}\widehat{\omega ^{}}_jU_j`$
Using the incompressibility ($`_jU_j=0`$), we finally obtain :

$$D_t\widehat{\omega ^{}}(𝐱,𝐤,t)=\widehat{F}(𝐱,t)\nu _tk^2\widehat{\omega }^{}(𝐱,𝐤,t)$$
(43)
with :

$`D_t`$ $`=`$ $`_t+𝐔_x\left(𝐔𝐤\right)_k,`$ (44)
## Appendix B Reconstruction formulae for the correlations

Using the formula 38 for the GT of space derivative, one can derive the subgrid scale velocity field in the (x,k) space with respect to the subgrid scale vorticity:

$`\widehat{\omega ^{}}(x,y,p,q,t)`$ $`=`$ $`\widehat{_xv^{}}(x,y,p,q,t)\widehat{_yu^{}}(x,y,p,q,t)`$
$`=`$ $`_x\widehat{v^{}}(x,y,p,q,t)_y\widehat{u^{}}(x,y,p,q,t)`$
$`=`$ $`ip\widehat{v^{}}(x,y,p,q,t)iq\widehat{u^{}}(x,y,p,q,t)+O(ϵ^{}),`$
where $`𝐤=(p,q)`$ and $`𝐱=(x,y)`$. Reversing the formula, we obtain the expression of the Gabor Transform of the velocity:

$`\widehat{u^{}}(x,y,p,q,t)`$ $`=`$ $`{\displaystyle \frac{iq}{p^2+q^2}}\widehat{\omega ^{}}(x,y,p,q,t)+O(ϵ^{}),`$ (47)
$`\widehat{v^{}}(x,y,p,q,t)`$ $`=`$ $`{\displaystyle \frac{ip}{p^2+q^2}}\widehat{\omega ^{}}(x,y,p,q,t)+O(ϵ^{}).`$ (48)
Consider now the following expression:

$`{\displaystyle \frac{1}{2}\left[\widehat{u^{}}(𝐱,𝐤,t)\widehat{\omega ^{}}(𝐱,𝐤,t)+\widehat{u^{}}(𝐱,𝐤,t)\widehat{\omega ^{}}(𝐱,𝐤,t)\right]𝑑𝐤}`$ (49)
$`=`$ $`{\displaystyle \mathrm{}\left[\widehat{u^{}}(𝐱,𝐤,t)\widehat{\omega ^{}}(𝐱,𝐤,t)\right]𝑑𝐤}`$
$`=`$ $`{\displaystyle \mathrm{}\left[f(ϵ^{}(𝐱𝐱^{}))e^{i𝐤(𝐱𝐱^{})}u^{}(𝐱^{},t)𝑑𝐱^{}f(ϵ^{}(𝐱𝐱^{\prime \prime }))e^{i𝐤(𝐱𝐱^{\prime \prime })}\omega ^{}(𝐱^{\prime \prime },t)𝑑𝐱^{\prime \prime }\right]𝑑𝐤}`$
$`=`$ $`{\displaystyle f(ϵ^{}(𝐱𝐱^{}))f(ϵ^{}(𝐱𝐱^{\prime \prime }))u^{}(𝐱^{},t)\omega ^{}(𝐱^{\prime \prime },t)\left(e^{i𝐤(𝐱^{\prime \prime }𝐱^{})}𝑑𝐤\right)𝑑𝐱^{}𝑑𝐱^{\prime \prime }}.`$
Using the definition of the Dirac function,

$$\frac{1}{(2\pi )^2}e^{i𝐤(𝐱^{\prime \prime }𝐱^{})}𝑑𝐤=\delta (𝐱^{}𝐱^{\prime \prime }),$$
(50)
and the fact that $`f^2=G`$, one simply gets

$$\frac{1}{(2\pi )^2}\mathrm{}\left[\widehat{u^{}}(𝐱,𝐤,t)\widehat{\omega ^{}}(𝐱,𝐤,t)\right]𝑑𝐤=\overline{u^{}\omega ^{}}(𝐱,t).$$
(51)
We may proceed along the same line for the average of the non-linear product of a large scale field with a subgrid scale field. Using the definition of the average $`\overline{𝐔\omega ^{}}(𝐱,t)`$ can be written as follows,

$$\overline{𝐔\omega ^{}}(𝐱,t)=f^2(ϵ^{}(𝐱𝐱^{}))𝐔(𝐱^{},t)\omega ^{}(𝐱^{},t)𝑑𝐱^{}.$$
(52)
Using the Taylor development of $`U`$ with respect to $`𝐱`$ at the first order, we obtain the first order approximation of this term:

$`\overline{𝐔\omega ^{}}(𝐱,t)`$ $`=`$ $`(𝐔\overline{\omega ^{}})(𝐱)+{\displaystyle f^2(ϵ^{}(𝐱𝐱^{}))(𝐱𝐱^{})𝐔(𝐱^{},t)\omega ^{}(𝐱^{},t)𝑑𝐱^{}}`$ (53)
$`=`$ $`(𝐔\overline{\omega ^{}})(𝐱)+O(ϵ^{}).`$
If we now apply the definition of the average for a product of two subgrid scales field (51) with the quantities $`1`$ et $`\omega ^{}`$, we finally obtain:

$`\overline{𝐔\omega ^{}}(𝐱,t)`$ $`=`$ $`𝐔{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle \mathrm{}\left[\widehat{\omega ^{}}(𝐱,𝐤,t)\widehat{1}(𝐱,𝐤,t)\right]𝑑𝐤}+O(ϵ^{})`$ (54)
$`=`$ $`𝐔{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle \mathrm{}\left[\widehat{\omega ^{}}(𝐱,𝐤,t)\widehat{f}(k)\right]𝑑𝐤}+O(ϵ^{}).`$ (55)
## Appendix C Numerical computation of the Reynolds Stress components

The subgrid scale field in physical space can be obtained from their discrete formula in the (x,k) by an integration with respect to $`𝐤`$:

$`\omega ^{}(𝐱,t)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^2f(0)}}{\displaystyle \widehat{\omega ^{}}(𝐱,𝐤,t)𝑑𝐤}`$
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^2f(0)}}{\displaystyle \underset{\alpha =1}{\overset{2N_p}{}}\widehat{\sigma }_\alpha (t)S_𝐱(𝐱𝐱_\alpha (t))\delta (𝐤𝐤_\alpha (t))d𝐤}`$
$`=`$ $`{\displaystyle \frac{1}{f(0)}}{\displaystyle \underset{\alpha =1}{\overset{2N_p}{}}}\widehat{\sigma }_\alpha (t)S_𝐱(𝐱𝐱_\alpha (t))`$
$`=`$ $`{\displaystyle \frac{1}{f(0)}}\left\{{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}\widehat{\sigma }_{\alpha _+}(t)S_𝐱(𝐱𝐱_{\alpha _+}(t))+{\displaystyle \underset{\alpha _{}=1}{\overset{N_p}{}}}\widehat{\sigma }_{\alpha _1}(t)S_𝐱(𝐱𝐱_\alpha _{}(t))\right\},`$
where

$$\underset{\alpha =1}{\overset{2N_p}{}}=\underset{\alpha _+=1}{\overset{N_p}{}}+\underset{\alpha _{}=1}{\overset{N_p}{}},$$
(57)
$`_{\alpha _+=1}^{N_p}`$ means sum over half the particles with the wavenumber $`𝐤_\alpha `$ and $`_{\alpha _{}=1}^{N_p}`$ is the sum over particles with an opposite wavenumber $`𝐤_\alpha `$. We choose the “positive” particles $`\alpha _+`$ with $`p_{\alpha _+}>0`$ and $`𝐤_\alpha =(p_{\alpha _+},q_{\alpha _+})`$. Because $`𝐱_{\alpha _+}=𝐱_\alpha _{}`$, $`\omega ^{}(𝐱,t)`$ can be written:

$$\omega ^{}(𝐱,t)=\frac{1}{f(0)}\underset{\alpha _+=1}{\overset{N_p}{}}\left(\widehat{\sigma }_{\alpha _+}(t)+\widehat{\sigma }_\alpha _{}(t)\right)S_𝐱(𝐱𝐱_{\alpha _+}(t)).$$
(58)
Since $`\omega ^{}(𝐱,t)`$ is real, the following property holds:

$$\widehat{\sigma }_{\alpha _+}(t)=\widehat{\sigma }_\alpha _{}^{}(t)$$
(59)
and the formula used to rebuild the vorticity field in physical space is :

$$\omega ^{}(𝐱,t)=\frac{2}{f(0)}\underset{\alpha _+=1}{\overset{N_p}{}}\mathrm{}\left[\widehat{\sigma }_{\alpha _+}(t)\right]S_𝐱(𝐱𝐱_{\alpha _+}(t))$$
(60)
where $`\mathrm{}\left[\widehat{\sigma }_{\alpha _+}(t)\right]`$ is the real part of $`\widehat{\sigma }_{\alpha _+}(t)`$. The same developments can be made for the two velocity component, using (47) and (48):

$`u^{}(𝐱,t)`$ $`=`$ $`{\displaystyle \frac{2}{f(0)}}{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}{\displaystyle \frac{q_{\alpha _+}}{p_{\alpha _+}^{}{}_{}{}^{2}+q_{\alpha _+}^{}{}_{}{}^{2}}}\mathrm{}\left[\widehat{\sigma }_{\alpha _+}(t)\right]S_𝐱(𝐱𝐱_{\alpha _+}(t))`$ (61)
$`v^{}(𝐱,t)`$ $`=`$ $`{\displaystyle \frac{2}{f(0)}}{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}{\displaystyle \frac{p_{\alpha _+}}{p_{\alpha _+}^{}{}_{}{}^{2}+q_{\alpha _+}^{}{}_{}{}^{2}}}\mathrm{}\left[\widehat{\sigma }_{\alpha _+}(t)\right]S_𝐱(𝐱𝐱_{\alpha _+}(t))`$ (62)
where $`\mathrm{}\left[\widehat{\sigma }_{\alpha _+}(t)\right]`$ is now the Imaginary part of $`\widehat{\sigma }_{\alpha _+}(t)`$.

The previous expressions can be used to compute velocity correlations. Using the analytical definition of these terms with respect the the subgrid scales field in (x,k) space (51,55), and using the definition of the subgrid scale discretization (16), we get:

$`\overline{u^{}v^{}}(𝐱,t)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle \mathrm{}\left[\widehat{u^{}}(𝐱,𝐤,t)\widehat{v^{}}(𝐱,𝐤,t)\right]𝑑𝐤}`$
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle }\mathrm{}[\left({\displaystyle \underset{\alpha =1}{\overset{2N_p}{}}}{\displaystyle \frac{iq}{p^2+q^2}}\widehat{\sigma }_\alpha S_𝐱(𝐱𝐱_\alpha )\delta (𝐤𝐤_\alpha )\right)`$
$`\left({\displaystyle \underset{\beta =1}{\overset{2N_p}{}}}{\displaystyle \frac{+ip}{p^2+q^2}}\widehat{\sigma }_\beta S_𝐱(𝐱𝐱_\beta )\delta (𝐤𝐤_\beta )\right)]d𝐤.`$
If all the wave-packets have a different wavenumber (at least on a domain equal to the “support” of $`S_x`$), eq C can be written:

$$\overline{u^{}v^{}}(𝐱,t)=2\underset{\alpha _+=1}{\overset{N_p}{}}\mathrm{}\left[\frac{q_{\alpha _+}p_{\alpha _+}}{(p_{\alpha _+}^{}{}_{}{}^{2}+q_{\alpha _+}^{}{}_{}{}^{2})^2}\widehat{\sigma }_{\alpha _+}\widehat{\sigma }_\alpha _{}S_𝐱^2(𝐱𝐱_{\alpha _+})\right].$$
(64)
Using the fact that $`\widehat{\sigma }_\alpha _{}`$ and $`\widehat{\sigma }_{\alpha _+}`$ are complex conjugate, we finally get:

$`\overline{u^{}v^{}}(𝐱,t)=2{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}{\displaystyle \frac{q_{\alpha _+}p_{\alpha _+}}{(p_{\alpha _+}^{}{}_{}{}^{2}+q_{\alpha _+}^{}{}_{}{}^{2})^2}}|\widehat{\sigma }_{\alpha _+}|^2S_𝐱^2(𝐱𝐱_{\alpha _+}),`$
$`\overline{u_{}^{}{}_{}{}^{2}}(𝐱,t)=2{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}{\displaystyle \frac{+q_{\alpha _+}^2}{(p_{\alpha _+}^{}{}_{}{}^{2}+q_{\alpha _+}^{}{}_{}{}^{2})^2}}|\widehat{\sigma }_{\alpha _+}|^2S_𝐱^2(𝐱𝐱_{\alpha _+}),`$ (65)
$`\overline{v_{}^{}{}_{}{}^{2}}(𝐱,t)=2{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}{\displaystyle \frac{+p_{\alpha _+}^2}{(p_{\alpha _+}^{}{}_{}{}^{2}+q_{\alpha _+}^{}{}_{}{}^{2})^2}}|\widehat{\sigma }_{\alpha _+}|^2S_𝐱^2(𝐱𝐱_{\alpha _+}).`$
The non-linear term involving the resolved scale velocity field can also be written in terms of wave-packets coordinates. Using eq. (55) and the definition of discretization (16) , one can write:

$`\overline{𝐔\omega ^{}}(𝐱,t)`$ $``$ $`{\displaystyle \frac{𝐔}{(2\pi )^2}}{\displaystyle \mathrm{}\left[\widehat{\omega ^{}}(𝐱,𝐤,t)\widehat{1}(𝐱,𝐤,t)\right]𝑑𝐤}`$ (66)
$``$ $`{\displaystyle \frac{𝐔}{(2\pi )^2}}{\displaystyle \mathrm{}\left[\widehat{\omega ^{}}(𝐱,𝐤,t)f(𝐤)\right]𝑑𝐤}`$
$``$ $`{\displaystyle \frac{𝐔}{(2\pi )^2}}{\displaystyle \mathrm{}\left[\underset{\alpha =1}{\overset{2N_p}{}}\widehat{\sigma _\alpha }f^{}(𝐤)S_𝐱(𝐱𝐱_\alpha )\delta (𝐤𝐤_\alpha )\right]𝑑𝐤}`$
$``$ $`2𝐔{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}\mathrm{}\left[\widehat{\sigma _{\alpha _+}}f^{}(𝐤_{\alpha _+})\right]S_𝐱(𝐱𝐱_{\alpha _+}).`$
## Appendix D Choice of the filter

The choice of the interpolating function $`S_x(𝐱)`$ (eq. 17) dictates the shape of the filter and, therefore, of the unknown function $`f(x)`$ in the definition of the Gabor Transform. The proof proceeds via the quantity $`\overline{\omega _{}^{}{}_{}{}^{2}}(𝐱,t)`$, by equating the formula of the filter discretized on a regular grid with a cell size $`\mathrm{\Delta }𝐱`$ and the formula 60:

$`\overline{\omega _{}^{}{}_{}{}^{2}}(𝐱,t)`$ $`=`$ $`{\displaystyle \underset{i}{}}\omega _{}^{}{}_{}{}^{2}(𝐱_i)f^2(𝐱𝐱_i)(\mathrm{\Delta }𝐱)^2`$ (67)
$`=`$ $`2{\displaystyle \underset{\alpha _+=1}{\overset{N_p}{}}}|\sigma _{\alpha _+}|^2S_𝐱^2(𝐱𝐱_{\alpha _+})`$ (68)
For this equality to be valid at all grid point $`X_i`$, the following link between $`f`$ and $`S`$ must hold:

$$f^2(𝐱)=C\times S_𝐱^2(𝐱).$$
(69)
Using now the normalization:

$`S_𝐱(0)`$ $`=`$ $`1.`$ (70)
$`{\displaystyle f^2(𝐱^{})𝑑𝐱^{}}`$ $`=`$ $`1.`$ (71)
the eq. 69 becomes:

$$f(𝐱)=f(0)S_𝐱(𝐱)$$
(72)
with

$$f(0)=\frac{1}{\sqrt{{\displaystyle S_𝐱^2(𝐱^{})𝑑𝐱^{}}}}.$$
(73)


text:
# Spin Environment Engineering for a Persistent Current Qubit

## I Introduction

The idea to perform computation with the help of quantum mechanics dates back to Feynman . Research in this direction got a considerable acceleration after the realisation  that quantum computers could solve in a polynomial time certain problems which require nonpolynomial time on classical computers. Two main obstacles to a universal quantum computer are decoherence and scalability. Decoherence, if sufficiently weak, can be dealt with by quantum error correcting codes, see Ref. for an experimental implementation. Scalability requires that a given technology can be upscaled to a quantum coherent circuit with hundreds or thousands of quantum gates, see Ref. for reviews on various proposals of scalable implementations. At the present stage solid state implementations have most to offer from the point of view of scalability but at the same time they are most likely to suffer from decoherence. The aim of this note is to reduce decoherence in the persistent current qubit (PCQ) proposal .

Two remarkable experiments  took place last year where a superconducting loop broken by Josephson junctions was forced into states which are quantum superpositions of two states with opposite macroscopic persistent currents: $`|`$ and $`|`$. The two considered combinations $`|\pm =|\pm |`$ were eigenstates separated by a gap of $`0.1K`$ so it is not surprising that each of them does not suffer from decoherence at the temperature of $`50mK`$. A more general superposition $`\alpha |++\beta |`$ would quickly decohere into a mixed state $`|\alpha |^2|++|+|\beta |^2||`$. To use a persistent current qubit for quantum computation any superposition must be free from decoherence. A design of a more microscopic scalable PCQ was proposed in Ref.. These authors went on to estimate decoherence times from different sources relevant for their design  and found out that the most dangerous are nuclear spins which couple to the magnetic field induced by persistent currents. Similar conclusion is reached in Ref.. In the calculations of Ref. the spins are treated as random static background fields as may be justified by the long relaxation time of nuclear spins which is of the order of minutes. Such a static background can be measured at the beginning of the quantum computation and compensated for by an adjustable counterterm in the qubit Hamiltonian. In this note we re-address this problem and find out that spins cannot be treated as a mere random background because the qubit state quickly gets entangled with the spins state. This observation is in agreement with Ref.. What is more, because of the large number of involved spins, even for the long relaxation time the spins background performs a random walk which leads to dephasing in the qubit state. The former more acute entanglement problem can be treated by exposing the nuclear spins to a strong external magnetic field $`B`$ parallel to the plane of the superconducting loop and penetrating through its thin layer, see Fig.1.

$`B`$ has to be less than the critical field in a type I superconducting film. If the $`N`$ nuclear spins are exposed to a magnetic field $`b`$ from the persistent current which is much weaker than $`B`$, as can be characterized by a quality factor

$$Q=\frac{B^2}{Nb^2},$$
(1)
then the entanglement between the qubit state and the spins state is negligible and leads to negligible partial decoherence. The means which treat entanglement incidentally also reduce the latter random walk problem, which can be furter eliminated by increasing the spin relaxation time $`T_r`$ which depends exponentially on temperature.

## II  Single spin environment

We begin with an elementary example where the qubit is coupled to just one environmental $`1/2`$ spin. For the purpose of this paper the spin-1/2 is a sufficient representation of the nuclear spin-9/2 of niobium. The spin is coupled to the external magnetic field $`B`$; its Hamiltonian is

$$H_\mathrm{E}=B\sigma _z^\mathrm{E}.$$
(2)
The current of the qubit induces a magnetic field which at the location of the environmental spin has strength $`b`$. For the time being we assume that the field is perpendicular to the external field in, say, $`x`$-direction. The Hamiltonian of interaction between the qubit and the spin is

$$H_{\mathrm{Q}\mathrm{E}}=b\sigma _z\sigma _x^\mathrm{E}.$$
(3)
The external magnetic field $`B`$, which is in plane of the superconducting loop, does not couple to the qubit state. For the sake of simplicity we set the qubit Hamiltonian $`H_\mathrm{Q}=0`$, what does not change the general conclusion of this note. The initial state of $`Q`$ and $`E`$ is a product

$$|\psi (0)\left(\alpha |+\beta |\right)|_\mathrm{E}.$$
(4)
The states $`|_{(\mathrm{E})},|_{(\mathrm{E})}`$ are $`\pm 1`$ eigenstates of $`\sigma _z^{(\mathrm{E})}`$. Under the total Hamiltonian $`H_{\mathrm{Q}\mathrm{E}}+H_\mathrm{E}`$ the initial $`|\psi (0)`$ evolves into

$`|\psi (t)`$ $`=`$ $`\alpha |e^{it(+b\sigma _x^\mathrm{E}+B\sigma _z^\mathrm{E})}|_\mathrm{E}+`$ (6)
$`\beta |e^{it(b\sigma _x^\mathrm{E}+B\sigma _z^\mathrm{E})}|_\mathrm{E}.`$
At the same time the initial pure reduced density matrix of the qubit

$$\rho (0)=\left(\begin{array}{cc}\alpha ^{}\alpha & \alpha \beta ^{}\\ \alpha ^{}\beta & \beta ^{}\beta \end{array}\right)$$
(7)
evolves into

$$\rho (t)=\left(\begin{array}{cc}\alpha ^{}\alpha & O(t)\alpha \beta ^{}\\ O(t)\alpha ^{}\beta & \beta ^{}\beta \end{array}\right),$$
(8)
where the overlap $`O(t)`$ is given by

$`O(t)`$ $``$ $`_\mathrm{E}|e^{+it(+b\sigma _x^\mathrm{E}+B\sigma _z^\mathrm{E})}e^{it(b\sigma _x^\mathrm{E}+B\sigma _z^\mathrm{E})}|_\mathrm{E}=`$ (10)
$`1{\displaystyle \frac{2b^2}{\mathrm{\Omega }^2}}\mathrm{sin}^2\sqrt{\mathrm{\Omega }}t.`$
Here $`\mathrm{\Omega }=\sqrt{B^2+b^2}`$. $`\rho (t)`$ is mixed when $`|O(t)|<1`$. The entanglement with a single spin is reversible because there are periodic revivals at the times when $`O(t)=1`$. We note that for $`b^2B^2`$ the overlap $`O(t)`$ remains $`1`$ and the qubit state stays pure.

Even just one environmental spin can lead to irreversible full decoherence when the spin interacts with its thermal bath. In the absence of the qubit, thermalisation of the spin density matrix $`\rho _\mathrm{E}`$ towards $`\frac{1}{2}I`$ can be described by a master equation

$$T_r\frac{d}{dt}\rho _\mathrm{E}=\sigma _{}\rho _\mathrm{E}\sigma _++\sigma _+\rho _\mathrm{E}\sigma _{}\rho _\mathrm{E},$$
(11)
where $`\sigma _{}=(\sigma _x^\mathrm{E}i\sigma _y^\mathrm{E})/2`$ are spin-$`1/2`$ lowering/raising operators and $`T_r`$ is the spin relaxation time of the order of minutes. This master equation can be formally unravelled by a stochastic evolution of the spin state. The spin state makes a down(up) flip within time $`dt`$ with a probability given by $`dt/T_r`$ times the probability that the spin is in the $`|_\mathrm{E}(|_\mathrm{E})`$ state. The master equation describes evolution of an average over such stochastic trajectories of the spin state.

Let us apply this stochastic unravelling to the state of the spin and the qubit. The $`|\psi (t)`$ in Eq.(6) can be rewritten as

$`|\psi (t)=\left({\displaystyle \frac{ib}{\mathrm{\Omega }}}\mathrm{sin}\mathrm{\Omega }t\right)\left(\alpha |\beta |\right)|_\mathrm{E}+`$ (12)
$`\left(\mathrm{cos}\sqrt{\mathrm{\Omega }}t+{\displaystyle \frac{iB}{\mathrm{\Omega }}}\mathrm{sin}\mathrm{\Omega }t\right)\left(\alpha |+\beta |\right)|_\mathrm{E}`$ (13)
The up transition at the time $`t`$ can be described as a projection of $`|\psi (t)`$ on $`_\mathrm{E}|`$ followed by a spin flip
$`|_\mathrm{E}|_\mathrm{E}`$. As a result $`|\psi (t)`$ jumps to

$$|\psi (t^+)=\left(\alpha |+\beta |\right)|_\mathrm{E}.$$
(14)
The resulting state of the system is the same as its initial state: this up transition causes no decoherence.

The probability of the down transition in the interval $`(t,t+dt)`$ is $`dt/T_r`$ times the probability that $`|\psi (t)`$ is in the $`|_\mathrm{E}`$ state which is $`|_\mathrm{E}|\psi (t)|^2(b^2/B^2)\mathrm{sin}^2Btb^2/2B^2`$, where we assume $`b^2B^2`$. The probability is $`dtb^2/2B^2T_r`$. The down transition is a projection on $`_\mathrm{E}|`$ followed by a spin flip $`|_\mathrm{E}|_\mathrm{E}`$. $`|\psi (t)`$ jumps to

$$|\psi (t^+)=\left(\alpha |_\mathrm{S}\beta |_\mathrm{S}\right)|_\mathrm{E}.$$
(15)
The resulting system state has a flipped sign of $`\beta `$ as compared to the initial system state (4). This sign flip also flips the sign of the off-diagonal elements in the reduced qubit density matrix, compare Eqs.(7,8). The evolution of $`\rho (t)`$ including spin relaxation is obtained after an average over stochastic trajectories. The off-diagonal elements of $`\rho (t)`$ average to zero on a timescale of $`\tau _B^{(2)}=2T_rB^2/b^2`$. This is the time of full decoherence due to spin relaxation. We note that for $`B^2b^2`$ the time is much longer than $`T_r`$.

## III  N-spin environment

Decoherence becomes more efficient when there are many copies of the environmental spin numbered by $`i=1,\mathrm{},N`$. Interaction Hamiltonians are $`H_{\mathrm{Q}\mathrm{E}_\mathrm{i}}=b_i\sigma _z\sigma _x^{\mathrm{E}_\mathrm{i}}`$. The initial product state

$$|\psi (0)=\left(\alpha |+\beta |\right)\underset{i=1}{\overset{N}{}}|_{\mathrm{E}_\mathrm{i}}$$
(16)
evolves into a $`|\psi (t)`$ such that the reduced state of the qubit has the form of Eq.(8) but with

$$O(t)=\underset{i=1}{\overset{N}{}}\left(1\frac{2b_i^2}{\mathrm{\Omega }_i^2}\mathrm{sin}^2\mathrm{\Omega }_it\right)$$
(17)
where $`\mathrm{\Omega }_i=\sqrt{B^2+b_i^2}`$. The early time overlap is $`O(t)=12Nb^2t^2`$, where $`b^2=_ib_i^2/N`$. $`O(t)`$ decays from $`1`$ on a timescale of $`\tau ^{(1)}=1/\sqrt{N}b`$. For $`b_i^2B^2`$ different factors in Eq.(17) oscillate with different frequencies dispersed in the range $`B\pm b^2/B`$. For $`t2\pi B/b^2`$ the phases of oscillating terms in different factors become effectively random and the overlap averages to

$$\overline{O}=1\frac{Nb^2}{B^2}=1\frac{1}{Q}.$$
(18)
For $`Q=B^2/Nb^21`$ the system state remains pure; there is only a negligible partial decoherence. Fig.2 shows $`O(t)`$ for a sample of $`b_i`$’s at two different values of $`B`$.

Let us include spin relaxation. To get full decoherence it is enough that only one out of $`N`$ spins makes the down transition so the decoherence time is $`N`$ times shorter than for a single spin, $`\tau ^{(2)}=(2T_rB^2/b^2)/N=2QT_r`$. $`Q1`$ makes this decoherence time is much longer than $`T_r`$.

It is worthwhile to compare the large $`Q`$ limit with the case of $`B=0`$ when Eq.(17) gives

$$O(t)=\underset{i=1}{\overset{N}{}}\mathrm{cos}2b_it.$$
(19)
The different factors, which oscillate with frequencies in a range proportional to $`b`$, go out of phase after a time of $`2\pi /b`$ when the overlap averages to zero. For $`B=0`$ the time $`\tau ^{(1)}=1/\sqrt{N}b`$ when $`O(t)`$ decays away from $`1`$ is a time of full decoherence.

## IV  Thermalized environment

So far we considered the single spin environment and the multispin environment with a fully polarized initial state. A thermal initial state is a weighted average over different partially polarized states, like e.g. $`|_{\mathrm{E}_1}|_{\mathrm{E}_2}\mathrm{}|_{\mathrm{E}_\mathrm{N}}`$. The overlap (17) is the same for any such partially polarized state as for the fully polarized state so the estimates of decoherence time for the polarized state go through unchanged for the thermal density matrix.

The mixedness of the initial environmental state does matter when we finally take into account the $`z`$-component of the qubit magnetic field at the location of spins which gives an extra term in the interaction Hamiltonian

$$H_{\mathrm{Q}\mathrm{E}}^z=\sigma _z\underset{i=1}{\overset{N}{}}b_i^z\sigma _z^{\mathrm{E}_\mathrm{i}}.$$
(20)
This Hamiltonian, when applied to a partially polarized spins state, gives $`H_{\mathrm{Q}\mathrm{E}}^z=\sigma _z_{i=1}^N|b_i^z|p_i`$, where $`p_i=\pm 1`$ is a polarization of the $`i`$-th spin. If $`p_i`$’s were static, then this $`H_{\mathrm{Q}\mathrm{E}}^z`$ could be measured at the beginning of the quantum computation and balanced by a $`\sigma _z`$-counterterm in the qubit Hamiltonian, as suggested in Ref.. The problem is that $`p_i`$’s are not static because they flip at the rate of one spin per $`T_r/N`$. The total polarization $`p(t)=_{i=1}^Np_i(t)`$ performs a random walk with one $`\pm 2`$ step every $`T_r/N`$ so that $`\overline{[p(t)p(0)]^2}=4tN/T_r`$, where the overline means here an average over random walks. The Hamiltonian $`H_{\mathrm{Q}\mathrm{E}}^z\sigma _zb[p(t)p(0)]`$ is random so the qubit states $`|`$ and $`|`$ accumulate opposite random phases $`\pm \varphi (t)`$, which grow as $`\overline{\varphi (t)^2}=4Nb^2t^3/3T_r`$. They lead to full dephasing after a dephasing time $`\tau ^{(3)}=(3\pi ^2T_r/Nb^2)^{1/3}`$. In this discussion we assume $`Q1`$ so that entanglement is negligible.

## V  Numerical estimates

The external magnetic field $`B`$ is bounded from above by the critical field of the type I superconducting material $`H_c`$. This parameter is the highest ($`0.206T`$) for niobium. In this respect niobium is much better than alluminium which has much lower $`H_c=0.01T`$. We choose niobium and set $`B=0.1T`$. There are materials like lead or tin that have abundant isotopes with zero nuclear spin. We do not consider them here because niobium and alluminium microtechnology is more advanced and, what is even more important, even lead or tin or their substrates would have spin impurities requiring environment engineering.

In the qubit design of Ref. the $`1\mu m\times 1\mu m`$ superconducting loop is made of an alluminium $`0.5\mu m\times 0.5\mu m`$ wire with persistent current of $`100nA`$. We suggest niobium instead of alluminium and the loop to be just $`0.1\mu m`$ thick. $`0.1\mu m`$ is roughly twice the penetration depth so that $`B`$ can penetrate through the superconductor. Note that the critical field is higher in such thin films than in bulk superconductor . Our superconducting loop will produce a flux of $`10^4\mathrm{\Phi }_0`$ which is an order of magnitude less than in the design of Ref.. The loop is also thinner so it contains fewer spins. The superconductor volume of $`10^{19}m^3`$ contains $`10^{10}`$ of niobium nuclear spins. They are subject to an average $`b^2=\mathrm{\hspace{0.33em}10}^{15}T^2`$. The quality factor $`Q=B^2/Nb^2=10^3`$ gives the overlap $`\overline{O}=0.999`$.

If the external field were $`B=0`$, then full decoherence would take place after

$$\tau ^{(1)}=\frac{1}{\sqrt{N}b}10^5s,$$
(21)
where we assume the nuclear spin magnetic moment of $`10^7Hz/T`$. For our $`B=0.1T`$ after $`\tau ^{(1)}`$ we get only negligible partial decoherence with $`\overline{O}=0.999`$. The full decoherence due to spin relaxation could happen after

$$\tau ^{(2)}=2QT_r10^5s,$$
(22)
where we assume $`T_r`$ to be of the order of minutes, but the dephasing time

$$\tau ^{(3)}=\left(\frac{3\pi ^2T_r}{Nb^2}\right)^{1/3}10^3s,$$
(23)
is much shorter. Even this $`1ms`$ time can be significantly increased by decreasing temperature because $`T_r\mathrm{exp}(\mathrm{\Delta }/k_\mathrm{B}T)`$, where $`\mathrm{\Delta }`$ is the superconducting gap. Reduction of temperature from $`50mK`$ down to $`5mK`$ increases $`\tau ^{(3)}`$ 25 times.

## VI  Conclusion

Maximalization of the quality factor $`Q=B^2/Nb^2`$ so that $`Q1`$ makes entanglement with spin environment negligible. The same technique can be used for any qubit design which involves superposition of states with different magnetic fields. According to our estimates, the present day persistent current qubit proposal  is not far from achieving a $`1ms`$ dephasing time due to the finite spin relaxation time $`T_r`$. $`1ms`$ is sufficient for the planned quantum computation at a $`GHz`$ rate. To achieve this desired decoherence time it is important to replace alluminium by niobium and make the circuit more microscopic than in Ref.. It is a common wisdom that more microscopic systems are better from the point of view of decoherence, compare e.g. our formula for $`\tau ^{(1)}`$. We showed that in addition to this much expected behaviour, there is a threshold below which a sufficiently microscopic qubit can be made free from entanglement with the spin bath.

Acknowledgements.  I would like to thank Diego Dalvit and Nikolay Prokof’ev for their comments on the manuscript. This work was supported in part by NSA.


text:
# Chandra X-ray Observations of Cygnus A and Pictor A

## 1. Introduction

The talk described observations of Cygnus A (z = 0.0562) and Pictor A (z = 0.035) with the Chandra X-ray Observatory. Since the early results from these observations are published (Wilson, Young & Shopbell 2000, 2001), we here confine ourselves to a summary of the main conclusions and a few remarks about other aspects of the results on Cygnus A, which will be amplified in a future paper.

## 2. Cygnus A

### 2.1. Cavity, Hot Spots, Jets and Nucleus

The Chandra image of Cygnus A is shown in Fig. 1. This figure shows only the region of the radio source and omits the large-scale, X-ray emitting gas which comprises the intracluster medium. The dominant feature of Fig. 1 is a limb-brightened, elliptical (presumably prolate-spheroidal in three dimensions) structure. The major axis coincides with that of the radio source and the radio hot spots are seen as compact X-ray sources at the ends of the major axis. This structure appears to be the observational manifestation of the cavity which is expected to be inflated in the intracluster medium by relativistic material which has passed through the hot spots at the ends of the jets. Its existence and dynamics were first inferred and explored in the classic paper by Scheuer (1974). The ratio of major to minor axis of the ellipse is $``$ 2 – 2.5, in agreement with the expectations of simple models (cf. Begelman & Cioffi 1989). Bright, curved bands of emission project along the minor axis of the ellipse. It is tempting to interpret these features as “belts” of gas extending around the equator of the prolate spheroidal structure.

A curved structure runs along the major axis from the nucleus to each pair of hot spots. These two features are apparently the X-ray manifestations of the jets. They are, however, broad, diffuse and of low contrast, quite unlike the sharply defined radio structure of the large scale NW (the near side of the radio source) jet (cf. Carilli et al. 1996). The X-ray jet almost disappears near the NW hot spots (A and B in the nomenclature of Hargrave & Ryle 1974). There is also an faint elliptical semi-ring in the X-ray image; this semi-ring extends $``$ 10′′ to the NE of hot spots A and B. The combination of this structure and the two hot spots resembles a set of ‘headphones’ (the hot spots mark the ‘earpieces’). Definitive statements cannot be made in view of the low contrast of the X-ray jets, but we speculate that either:

a) the NW jet is continuous and steady but much broader than would be inferred from the radio maps. In this interpretation, the hot spots and the semi-ring reflect the interaction of the outer annulus of the jet with the intracluster gas. The observed radio jets and hot spots might then be effects of limb-brightening at the edges of the broad jet (cf. Blandford  who suggests that the hot spots might be a ring vortex seen in projection). However, the great brightness of the hot spots and the close agreement between their radio synchrotron and X-ray synchrotron self-Compton emissions (see next subsection) challenge a limb-brightening interpretation. Alternatively, an absence of radio emission from the center of the jet may mean that the jet is actually hollow or that the jet’s center comprises material unconducive to radio emission.

or: b) the NW jet is intrinsically narrow, as inferred from the radio maps, but the direction in which it is launched is continuously changing (precession of the nozzle?), so that its terminus traces out a circle where it meets the intracluster gas. The elliptical semi-ring is then the projection of this terminal circle, and the hot spots are two locations along the circle where the radio emission is unusually bright, as a result of projection or other effects.

The X-ray nucleus is extended in the SE - NW direction. Its spectrum is shown in Fig. 2. The high energy continuum may be modeled as an absorbed power law, as known from previous work (e.g. Arnaud 1996). Notable features include: a) a K$`\alpha `$ line from neutral or lowly ionized iron (seen in the observed spectrum near 6.0 keV), b) strong K edge absorption from neutral or lowly ionized iron (observed near 7.0 keV), and c) emission at energies below 2 keV in excess of that expected from the absorbed power law. Emission lines, which are confirmed in an independent spectrum with higher signal to noise at low energies, are apparent and tentatively identified with K$`\alpha `$ transitions of silicon, neon and possibly magnesium. The soft X-ray emission apparently originates from an extended region viewed through a relatively low absorbing column.

### 2.2. Synchrotron Self-Compton Emission from the Hot Spots

Fig. 3 (left) shows X-ray contours on a grey scale of a 6 cm image with resolution 0$`.^{\prime \prime }`$35 (Perley, Dreher & Cowan 1984) in the vicinities of the western hot spots (A and B). The extents and morphologies of the X-ray and radio hot spots are very similar, with the directions of elongation agreeing to within a few degrees. Similar results are obtained for hot spots D and E. The X-ray spectra of hot spots A and D have been modeled with an absorbed power law. In both cases, the absorbing column is NH = 3.3 $`\times `$ 1021 cm-2, in excellent agreement with the Galactic column in the direction of Cygnus A. The photon indices are also similar at $`\mathrm{\Gamma }`$ = 1.8 $`\pm `$ 0.2.

Thermal models for the X-ray emission may be ruled out as they require too high gas densities. Fig. 3 (right) shows the results of calculation of a synchrotron self-Compton (SSC) model for hot spot A. As may be seen, the predicted SSC radiation is in excellent agreement with the Chandra-observed spectrum for a magnetic field of 1.5 $`\times `$ 10-4 gauss. This value may be compared with the equipartition value of 2.8 $`\times `$ 10-4 gauss in hot spot A, calculated assuming no relativistic protons, the broken power-law spectra of Carilli et al. (1991), low frequency cut-offs at 10 MHz and high frequency cut-offs at 400 GHz. The uncertainty in the magnetic field strength estimated from the SSC model is estimated to be a few tens of percent.

The most straightforward interpretation is that the relativistic gas is an electron-positron plasma and is close to equipartition with the magnetic field. If, on the other hand, the energy in relativistic protons dominates that in relativistic electrons and the X-rays are still SSC radiation, the energy density in relativistic protons must exceed that in the magnetic field. It is notable that the magnetic field cannot be less than 1.5 $`\times `$ 10-4 gauss since the SSC radiation would then exceed the observed X-radiation.

The alternative is that B $`>`$ 1.5 $`\times `$ 10-4 gauss in which case the predicted SSC emission would be too weak to account for the observed X-ray emission. The X-rays would then have to be synchrotron radiation. However, in view of the excellent agreement between the SSC model and the observations, we consider a synchrotron model implausible.

## 3. Pictor A

Fig. 4 shows the Chandra image of the nucleus, western jet and western hot spot of Pictor A (much fainter X-ray emission is found to the east of the nucleus - see Wilson, Young & Shopbell 2001). The jet extends at least 1$`.^{}`$9 (110 kpc) westward from the nucleus and is spatially coincident with the faint radio jet (Perley, Röser & Meisenheimer 1997). The X-ray jet “points” at the western hot spot, which is 4$`.^{}`$2 (240 kpc) from the nucleus. When the nuclear X-ray source is shifted slightly to coincide with the nuclear radio source, the peak of the X-ray emission from the hot spot is within 1′′ of its radio peak.

Images of the western hot spot at radio, optical (Röser 1989; Perley, Röser & Meisenheimer 1997) and X-ray wavelengths with similar resolutions are shown in Fig. 5. As may be seen, the overall X-ray morphology is remarkably similar to the radio and optical. The high linear polarization of the hot spot at both radio and optical wavelengths indicates that the emission is synchrotron radiation in these wavebands.

Fig. 6 shows the broad band spectrum of the western hot spot. The radio spectrum of the hot spot is well described by a power law with $`\alpha `$ = 0.740 $`\pm `$ 0.015 (Meisenheimer, Yates & Röser 1997), but there must be a break or turnover in the spectrum at 1013-14 Hz to accommodate the near infrared and optical measurements. It is also apparent that the X-ray spectrum is not a simple extension of the radio and optical measurements to higher frequencies.

Wilson, Young & Shopbell (2001) reached the following conclusions about the X-ray emission of the hot spot, assuming that its bulk outward motion is non-relativistic. (i) Inverse Compton scattering of the synchrotron radio photons by the relativistic electrons responsible for the radio emission (i.e. a synchrotron self-Compton model) may be ruled out for the hot spot’s X-ray emission, as the predicted spectrum differs from that observed. (ii) A successful inverse Compton model may be constructed for the X-ray emission by invoking a low energy population of relativistic electrons with an energy index appropriate to the X-ray spectral index. However, the magnetic field must be reduced to $``$ 1% of equipartition, and there is no evidence for such an electron population. (iii) The emission may be reproduced in a composite synchrotron plus synchrotron self-Compton model (solid line in Fig. 6). However, the model is contrived, requiring similar fluxes from the two components in the Chandra band. (iv) Synchrotron radiation is possible, but the electron population must be distinct from that which produces the radio and optical emission. However, the X-ray spectrum is in excellent agreement with the expected form if the electrons are accelerated by strong shocks. Relativistic electrons may also result from a ‘proton induced cascade’ (e. g. Mannheim, Krülls & Biermann 1991).

If the jet is non-relativistic, inverse Compton scattering is an implausible model for its X-ray emission, for it requires a magnetic field a factor of 30 below equipartition. Further, it is hard to understand why the jet is brighter than the lobe (the opposite of the situation in the radio) in such a model. If the jet is relativistic, these difficulties are eased, and we consider inverse Compton scattering by such a jet off the microwave background a viable mechanism. However, the magnetic field must still be well below equipartition for plausible angles between the jet and the line of sight in this lobe-dominated, FRII radio galaxy. Synchrotron radiation is another possibility.

## 4. The difference between the hot spots of Cygnus A and Pictor A

It is striking that the X-ray emissions from the hot spots of Cygnus A, 3C 123 (Hardcastle, Birkinshaw & Worrall 2000) and 3C 295 (Harris et al. 2000) conform to an SSC model with a magnetic field close to equipartition (calculated assuming that only electrons are present). In three other galaxies - Pictor A, 3C 120 (Harris et al. 1999) and 3C 390.3 (Harris, Leighly & Leahy 1998) - the X-ray emission is orders of magnitude too strong to be SSC emission with an equipartition magnetic field. While this is a small sample, we (Wilson, Young & Shopbell 2000) noted that 3C 295 and Cygnus A are in clusters with prominent cooling flows and 3C 123 is in a cluster with strong, extended X-ray emission and thus may be within a cooling flow. In contrast, Pictor A, 3C 120 and 3C 390.3 are not in cooling flows. Such an environmental difference could affect the X-ray power from the hot spots in two ways: (a) the presence of high density surrounding gas may inhibit production of X-ray synchrotron emitting electrons in hot spots, or (b) the environments of the “overluminous” X-ray hot spots might have such low gas density that the bulk outward motion of the hot spots is relativistic. Luminous X-ray emission could then originate through inverse Compton scattering of the microwave background radiation, as we have discussed for the jet of Pictor A. In the second case, objects with high X-ray to radio flux ratios would tend to be viewed more pole-on than objects with low X-ray to radio flux ratios. Chandra observations of a larger sample of hot spots should cast light on these issues.

This research was supported by NASA grant NAG 81027 and by the Graduate School of the University of Maryland. We are grateful to Rick Perley for providing his radio images.

## References

Arnaud, K. 1996, ‘X-ray evidence for a buried quasar in Cygnus A’, in Cygnus A - Study of a Radio Galaxy, ed. C. L. Carilli and D. E. Harris (Cambridge University Press), 51–55

Begelman, M. C. & Cioffi, D. F. 1989, ‘Overpressured Cocoons in Extragalactic Radio Sources’ ApJ, 345, L21–L24

Blandford, R. D. 1996, ‘The Untangling of Cygnus A’, in Cygnus A - Study of a Radio Galaxy, ed. C. L. Carilli and D. E. Harris (Cambridge University Press), 264–274

Carilli, C. L., Perley, R. A., Dreher, J. W. & Leahy, J. P. 1991, ‘Multifrequency Radio Observations of Cygnus A: Spectral Aging in Powerful Radio Galaxies’, ApJ, 383, 554–573

Carilli, C. L., Perley, R. A., Bartel, N. & Dreher, J. W. 1996, ‘The Jets in Cygnus A: from pc- to kpc-scales’ in Cygnus A \- Study of a Radio Galaxy, ed. C. L. Carilli and D. E. Harris (Cambridge University Press), 76–85

Hardcastle, M. J., Birkinshaw, M. & Worrall, D. M. 2000, ‘A Chandra detection of the radio hotspot of 3C123’, MNRAS, submitted (2000)

Hargrave, P. J. & Ryle, M. 1974, ‘Observations of Cygnus A with the 5-km Radio Telescope’ MNRAS, 166, 305–327

Harris, D. E., Leighly, K. M. & Leahy, J. P. 1998, ‘X-ray Emission from a Radio Hot Spot in 3C390.3: Evidence for the Deflection of a Radio Jet by a Neighboring Galaxy’ ApJ, 499, L149–L152

Harris, D. E., Hjorth, J., Sadun, A. C., Silverman, J. D. & Vestergaard, M. 1999, ‘X-ray Emission from the Radio Jet in 3C120’, ApJ, 518, 213–218

Harris, D. E. et al. 2000, ‘Chandra X-ray Detection of the Radio Hot Spots of 3C295’, ApJ, 530, L81–L84

Mannheim, K., Krülls, W. M. & Biermann, P. L. 1991, ‘A novel mechanism for non-thermal X-ray emission’, A&A, 251, 723–731

Meisenheimer, K., Yates, M. G. & Röser, H-J. 1997, ‘The synchrotron spectra of radio hot spots. II. Infrared imaging’, A&A, 325, 57–73

Perley, R. A., Dreher, J. W. & Cowan, J. J. 1984, ‘The Jet and Filaments in Cygnus A’ ApJ, 285, L35–L38

Perley, R. A., Röser, H-J. & Meisenheimer, K. 1997, ‘The Radio Galaxy Pictor A - a Study with the VLA’ A&A, 328, 12–32

Röser, H-J. 1989, ‘Continuum Observations of Hot Spots at Wavelengths $`<`$ 1 cm’ in Lecture Notes in Physics vol. 327, Hot Spots in Extragalactic Radio Sources, eds. Meisenheimer, K. & Röser, H.-J. (Springer-Verlag) 91–114

Scheuer, P. A. G. 1974, ‘Models of Extragalactic Radio Sources with a Continuous Energy Supply from a Central Object” MNRAS, 166, 513–528

Wilson, A. S., Young, A. J. & Shopbell, P. L. 2000, ‘Chandra Observations of Cygnus A: Magnetic Field Strengths in the Hot Spots of a Radio Galaxy’ ApJ, 544, L27–L30

Wilson, A. S., Young, A. J. & Shopbell, P. L. 2001, ‘Chandra X-ray Observations of Pictor A: High Energy Cosmic Rays in a Radio Galaxy?’ ApJ(in press for Feb 1 2001 issue; also available as astro-ph/0008467)

Figure Captions

Figure 1. A grey scale representation of the Chandra X-ray image of Cygnus A. The shading is proportional to the square root of the intensity. Coordinates are for epoch J2000.0 throughout this paper.

Figure 2. Chandra X-ray spectrum of the nucleus of Cygnus A. The x axis represents observed photon energy. The upper panel shows the data \[points with error bars\] and a model of an absorbed power law spectrum folded through the instrumental response (solid line). The lower panel shows the deviations of the observed spectrum from the model. Our study of the nuclear spectrum is in collaboration with Keith Arnaud.

Figure 3. a: X-ray emission (contours) superposed on a 6 cm VLA radio map (grey scale) of the region of the western hot spots (A, the brighter, and B, $``$ 6′′ SE of A) in Cygnus A. Contours are plotted at 2, 4, 8, 12, 16, 24 and 32 counts per pixel (0$`.^{\prime \prime }`$5 $`\times `$ 0$`.^{\prime \prime }`$5). The grey scale is proportional to the square root of the radio brightness. b: Spectrum of hot spot A. The points show the radio fluxes and the line through them the model of the synchrotron radiation. The “bow tie” is the Chandra measured boundary of the X-ray spectrum (these error lines are 90% confidence after freezing NH at its best fit value, which coincides with the Galactic column). The solid line is the predicted SSC spectrum for $`\gamma _{\mathrm{min}}`$ = 1 and the dashed line for $`\gamma _{\mathrm{min}}`$ = 100.

Figure 4. Grey scale representation of the full resolution Chandra image of the nucleus, jet and western hot spot of Pictor A.

Figure 5. The morphology of the western hotspot of Pictor A, at 3.6 cm radio (left panels; PRM), R band optical (center panels; Röser (1989) and PRM) and X-ray (right panels; this paper) wavelengths. The upper panels are grey scale representations and the lower panels are contour plots in which the contours are logarithmically spaced and separated by a factor of 2. The radio and optical images have a FWHM resolution $``$ 1$`.^{\prime \prime }`$5 (shown as the filled circles at the bottom right) and the X-ray image has a FWHM of 1$`.^{\prime \prime }`$2 $`\times `$ 1$`.^{\prime \prime }`$9 (the PSF is shown at the bottom right; no smoothing has been applied to the X-ray image).

Figure 6. The broad band spectrum, plotted as $`\nu S_\nu `$, of the western hotspot of Pictor A from radio through X-ray wavelengths. Note that the scale along the y axis is expanded compared to that along the x axis. The circular and triangular symbols are measured data points, while the “bow tie” represents the Chandra spectrum. The lines represent various models discussed in Wilson, Young & Shopbell (2001).