text:
# Low Frequency VLA Observations of Abell 754: Evidence for a Cluster Radio Halo and Possible Radio Relics

## 1 INTRODUCTION

Clusters of galaxies are the most massive gravitationally bound objects in the universe which are in a state of quasi-equilibrium. Our current understanding of structure formation suggests that clusters form via merging of smaller entities, such as galaxy groups and small clusters. During a merging process clusters are able to release considerable energy into particle acceleration via shock waves and turbulence (De Young, 1992). The accelerated relativistic particles have a short radiative lifetime but can be reaccelerated by recurrent merger events. For several galaxy clusters, radio observations of the relativistic plasma reveal the presence of extended ($`>`$ 500 kpc) diffuse emission which has a steep spectral index and no observable optical counterpart (Feretti & Giovannini, 1996, and references therein). The synchrotron nature of the emission allows us to trace large regions of relativistic charged electrons and magnetic fields within the X-ray emitting intracluster medium (ICM). In general, this diffuse emission appears to fall within two categories: halos which are centrally located, regular in shape and unpolarized, and relics which are located in the peripheral regions of the cluster, irregular in shape and generally highly polarized (Feretti & Giovannini, 1996).

The origin of these regions of diffuse emission is still a matter of great debate. The large physical scales involved combined with the relatively short radiative lifetimes suggest that the electrons must be re-accelerated within the ICM (Jaffe, 1977). Theoretical models propose that this particle acceleration may occur in the turbulence and shocks associated with major cluster merger events (De Young, 1992; Tribble, 1993; Feretti, 1999; Brunetti et al., 2001), although alternative models include secondary particle production from proton-proton interactions (Dennison, 1980), and particle diffusion from nearby head-tail galaxies (Giovannini et al., 1993). Recent work suggests that there is a positive correlation between the presence of radio halos and relics, and current (or recent) merger activity associated with clusters (Feretti & Giovannini, 1996; Feretti, 1999).

Clusters of galaxies contain a considerable fraction of hot plasma, the ICM, which is heated up to several $`10^710^8`$ K due to virialization. This hot thermal plasma traces the gravitational potential of clusters and is thus an ideal tool to determine the dynamical state of a cluster. Because of its high temperature, which corresponds to several keV in energy, the ICM plasma is only directly observable with X-ray telescopes. The X-ray emission allows us not only to trace the potential but also to trace the presence of non-uniform heating, which is caused by adiabatic compression or shock waves from merger events. This kind of heating, which was difficult to detect with past X-ray telescopes is now observable in detail with the latest X-ray satellites, such as Chandra (Vikhlinin, Markevitch, & Murray, 2000a, b) or XMM-Newton (Arnaud et al., 2000).

In this paper we present new radio observations of the galaxy cluster Abell 754 (Abell, Corwin & Olowin, 1989), which was observed at 74 and 330 MHz with the NRAO VLA observatory111The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.. X-ray observations of this cluster (Henry & Briel, 1995; Henriksen & Markevitch, 1996), which trace the thermal ICM, indicate that it is undergoing a violent merger event. Recent 3D MHD/N-body simulations of Roettiger et al. (1999a, b) indicate that shocks and turbulence associated with the cluster merger could provide the magnetic field amplification and particle re-acceleration necessary to generate diffuse radio relic emission from regions of the ICM. Further, the energy input from the merger may accelerate relativistic particles in the cluster center to produce a radio halo. Therefore A754 seems to be the ideal host for diffuse radio emission. However, up to now radio observations have not been able to provide an unambiguous detection of such emission. Its detection is very important since it strongly supports the hypothesis of a connection between radio halos/relics and merger activity.

### 1.1 Past Observations of A754

At a redshift of z = 0.054 (Bird, 1994), A754 is a rich galaxy cluster, $`\sigma _v`$= 900 km $`\mathrm{s}^1`$ (Zabludoff & Zaritsky, 1995), which has been extensively studied in the optical, radio and X-ray bands. Observations in the 0.5-10 keV ASCA band (Henriksen & Markevitch, 1996) find an X-ray luminosity of $`6.4\times 10^{44}h_{75}^2\mathrm{ergs}\mathrm{s}^1`$ and an emission-weighted temperature of 8.5$`\pm `$0.5 keV for this cluster. Maps of the X-ray surface brightness and X-ray temperature from ROSAT (Henry & Briel, 1995), as well as X-ray temperature maps from ASCA (Henriksen & Markevitch, 1996) show temperature variations that indicate that this cluster is going through or has recently gone through a violent merger event.

A754 has been the object of low frequency radio observations for over 20 years. In no reported observation has a relic been detected and evidence for a radio halo has been inconclusive. Extended emission was reported at 2.7 GHz (Wielebinski et al., 1977) and at 408 MHz (Mills, Hunstead & Skellern, 1978). In both cases the resolution of 3 to 4 made it difficult to distinguish apparent extended flux from the blended images of point sources associated with galaxies in or behind the cluster. A higher resolution observation at 610 MHz (Harris et al., 1980) showed no extended emission, yet the authors conclude that a halo likely exists because, with just the point sources they detected, they were unable to account for all emission previously reported at lower frequencies. This issue was later re-investigated with the VLA in a B and C array observation (40 minutes combined) at 330 MHz (Giovannini, Tordi & Feretti, 1999; Giovannini & Feretti, 2000). This observation detected no extended flux in the central regions of this cluster. We note however that the short integration times obtained at 330 MHz would provide very poor $`uv`$-coverage which could easily lead to missing large scale structure in the maps.

## 2 Results

### 2.1 The Data and Reduction Methods

VLA C-array observations of A754 were conducted on March 21, 2000 simultaneously at 330 MHz and 74 MHz. The observational parameters are summarized in Table 1. At both frequencies, observations were conducted in spectral line mode to reduce the effects of bandwidth smearing and allow for more accurate and less costly radio frequency interference (RFI) excision. Ultimately 8% and 23% of the data were flagged at 330 MHz and 74 MHz, respectively. Cygnus A was used as a flux density, bandpass, and initial phase calibrator, and we estimate flux densities we quote as accurate to $`15\%`$. Successive loops of self-calibration (AIPS task CALIB) and non-coplanar, wide-field image deconvolution (AIPS task IMAGR) were used to mitigate confusion and achieve maximum sensitivity in our images (Perley, 1999). The realized rms sensitivities of $`6.5`$ mJy/beam and $`200`$ mJy/beam at 330 and 74 MHz, respectively, are both approximately ten times higher than the expected thermal noise levels. We attribute the difference to a combination of confusion and poorly understood, broad-band, mainly VLA-generated RFI which especially affects low-frequency observations in the compact C and D configurations.

### 2.2 The Radio Halo at 74 MHz and 330 MHz.

Our 74 MHz and 330 MHz images are presented as contours in Figures 1 and 2, respectively, while a 1.4 GHz image from the NRAO VLA sky survey (NVSS:Condon et al. (1998)) is superimposed on Figure 1 as grey-scale. At both frequencies we detect extended emission in the central region of the cluster which is consistent with a radio halo, although the 5’ resolution at 74 MHz is clearly too poor on its own to distinguish between genuine halo-emission and a blend of point source emission. However, analysis of the 74 MHz image in conjunction with information from the 330 MHz and 1.4 GHz data indicates that there is clearly excess 74 MHz emission beyond just a blending of the four main NVSS point sources in the halo region. To estimate their 74 MHz contribution we determined their spectral indices between 330 MHz and 1.4 GHz and extrapolated these to 74 MHz. The prediction for the combined flux density of these four point sources at 74 MHz is 2.9 Jy, while we detect 7 Jy of emission at the cluster center. Therefore we detect an excess of about 4 Jy of 74 MHz flux in addition to the point sources.

The case for extended emission is even stronger at 330 MHz (Figure 2) where the resolution is sufficient to clearly separate the brightest point sources from the diffuse emission in the region. Apart from the four main point sources, we measure a total of 750 mJy of flux density in the region. Of course it is possible that even this emission is made up of a blending of even more faint point sources in the center of the cluster which are undetected at 1.4 GHz, and close inspection of the NVSS image does suggest more faint point sources in the region (Figure 3). However, there are regions of the halo which appear strongly at 330 MHz but not at all at 1.4 GHz. If we examine region A which we define in Figure 2, the fact that it has a peak flux density of 49 mJy/beam yet contains no observable flux at 1.4 GHz requires a spectral index of at least 1.5, $`\mathrm{S}_\nu \nu ^\alpha `$. Taking the halo as a whole, our estimates of 4 Jy at 74 MHz and 750 mJy at 330 MHz indicate a spectral index of $`1.1`$. This is unusual for point sources associated with radio galaxies. The four brighter point sources all have spectral indices ranging from 0.5 to 1.0. We therefore assert that while some of this emission could be due to point sources, a significant fraction corresponds to a diffuse cluster radio halo with 3$`\sigma `$ major and minor axes of $`450^{\prime \prime }\times 400^{\prime \prime }`$, respectively. At the redshift of the cluster this corresponds to a scale of $`430\times 380h_{75}^1`$ kpc, a factor of 2-3 smaller than the largest known giant Mpc halos but still not the smallest halo known (Feretti, 1999). From our 330 MHz image we estimate the halo center at $`09^h08^m43^s`$,$`09^{}36^{}57^{\prime \prime }`$, indicating it is shifted westward of the X-ray maximum and also of the cluster center, which we place at $`09^h09^m03^s`$,$`09^{}39^{}46^{\prime \prime }`$ with reference to an elliptical beta-model fit. We note that the $`\alpha `$=1.1 spectrum of the diffuse emission we infer is consistent with previous estimates (Harris et al., 1980) though the total flux levels are still somewhat ($`40\%`$) lower than predicted.

Initially puzzling, 330 MHz observations by Giovannini & Feretti (2000) at comparable rms noise failed to detect the halo emission. However, a plausible explanation emerges in light of the fundamentally different spatial characteristics of the two data sets. Our observations were 3 hours in the VLA C array, while theirs was a combination of B and C array snapshots (40 min total) with consequently less thorough $`uv`$ coverage. A possible explanation is therefore that their superior resolution provided by B-array baselines sets a lower basic confusion level and thus provided a nominal sensitivity comparable to our longer, lower resolution integration. However, their snapshot C array $`uv`$ coverage would at the same time have been much less sensitive to the large-scale extended emission we detect. Indeed imaging of a snapshot sub-section of our own C array data failed to detect the halo.

### 2.3 Discovery of Radio Relics at 74 MHz.

Our 74 MHz image (Figure 1) reveals two, newly detected objects without counterparts at either 330 MHz or 1.4 GHz which we label as the “East Relic” and “West Relic” (see Table 2). While the $`6\sigma `$ (integrated flux $`1.45`$ Jy) and $`4\sigma `$ ($`850`$ mJy, unresolved) detections are weak, it is nevertheless highly unusual to detect sources at 74 MHz which are not also seen by the VLA at either 330 MHz or 1.4 GHz. In fact all other objects in a 4$`{}_{}{}^{}\times 4^{}`$ 74 MHz image that appear at the brightness of the fainter West Relic or higher have clear counterparts in the 1.4 GHz NVSS images, and 330 MHz counterparts within the $`2^{}`$ radius region in which the 74 and 330 MHz fields of view overlap. This suggests an unusually high spectral index for these objects. In order to explain the lack of a detection at above three times the rms sensitivity in the 330 MHz or 1.4 GHz images, the “East Relic” would need a spectral index of at least 1.8 or 1.4 respectively. Such high spectral indices are typical of cluster radio relics.

There is further evidence that these objects detected at 74 MHz are indeed radio relics. First the morphology of the “East Relic”, with an elongation perpendicular to its separation from the cluster X-ray maximum is highly suggestive of a relic. Its major and minor axes are $`600^{\prime \prime }\times 230^{\prime \prime }`$, respectively, which corresponds to scale sizes of $`580\times 220h_{75}^1`$ kpc at the redshift of the cluster. The reality of the mainly unresolved “West Relic” is clearly much more speculative, though it is noteworthy that its location with respect to the cluster X-ray maximum is opposite to that of the “East Relic”. An overlay of the positions of these relics with ROSAT hard X-ray data (Figure 4) reveals that the East relic falls exactly at the location of flattened, compressed X-ray contours. Such a rapid gradient in X-ray surface brightness is expected at the location of merger shock waves. Further, Roettiger, Stone & Mushotzky (1998) have predicted locations of shock waves based on their model of the cluster merger in A754, and it is especially noteworthy that the location of the East relic coincides with a jump in the X-ray temperature predicted for the shock waves in this cluster. We note however that the signal-to-noise ratio on both of these new objects is still very low and thus the objects (the West relic in particular) should be regarded with caution. We plan to undertake deep follow-up observations of these sources with the VLA at higher resolution to confirm these detections and provide tighter constraints on their spectra. Additional follow-up observations will be undertaken to identify the compact X-ray source located at the position of the East relic and determine if the two objects are connected.

## 3 Discussion

In the numerical model of Roettiger, Stone & Mushotzky (1998) a large subcluster (cluster to subcluster mass ratio 2.5:1) has passed from East to West slightly off-axis from the core of the A754 system. This merger event results in the formation of two shock waves which are visible as temperature jumps in Fig. 2a of their work. The Eastern shock heats the infalling, stripped subcluster gas in a relatively small area. The location of this shock at the ram pressure flattened edge of the X-ray peak is coincident with the observed location of the East relic source. The western shock heats the outer regions of the main cluster ICM and is much more extended. In principle, this shock could have caused the West relic, but we note that the West relic is further away from A754’s X-ray peak ($`1.8h_{75}^1\mathrm{Mpc}`$) than the shock wave in the numerical model ($`1.2h_{75}^1\mathrm{Mpc}`$). The measured (Henry & Briel, 1995; Henriksen & Markevitch, 1996) and simulated (Roettiger, Stone & Mushotzky, 1998) temperature maps show that mainly the outer western region is heated by this shock wave, but not the innermost cluster core region.

The temperature variation suggests that the impact of the infalling gas combined with the change in the gravitational potential of the primary cluster due to the core passage of the dark matter subcluster was not violent enough to produce a shock wave in the dense main cluster core. The core was only compressed adiabatically, and relaxed nearly to its original temperature in a later merger stage. When the impact moved into outer, less dense regions, a shock wave was created and heated the gas non-adiabatically. Such emergence of outgoing shock waves is observed in numerical simulations of cluster merger events (Schindler, private communication).

A fraction of the energy dissipated in shock waves is often transfered into relativistic particle populations. Accelerated relativistic electrons have short radiative lifetimes ($`0.1`$ Gyr) and should therefore produce the radio relic emission close to the current location of the shock waves. The electron acceleration required to produce the relic emission could result from Fermi-I diffusive shock acceleration of ICM electrons (Enßlin et al., 1998; Roettiger et al., 1999b), or adiabatic energization of relativistic electrons confined in fossil radio plasma, released by a former active radio galaxy (Enßlin & Gopal-Krishna, 2001).

Observationally, the radio halo of A754 appears to be displaced from the main cluster core in the direction of the West relic. This suggests that the halo emission results only from regions which were heated by the re-emerging western shock wave after core passage. Since the radiative cooling time of the radio halo electrons observed at 330 MHz is smaller than 0.3 Gyr, but the shock passage was more than 0.3 Gyrs ago (at least in the numerical model), the radio electrons must have been re-populated after shock passage. In situ Fermi II re-acceleration (powered by some residual plasma turbulence after the shock passage) of shock accelerated electrons could be able to re-accelerate the electrons (Brunetti et al., 2001), if the turbulence is able to remain for a few 100 Myr. Alternatively, relativistic protons might be efficiently accelerated at shock waves and would have lifetimes comparable to the cluster age. Since these protons can seed relativistic electrons and positrons into the ICM by hadronic collisions with the ICM ions, they might be responsible for the radio halo emission which glows for a long time after the shock passage (Dennison, 1980; Vestrand, 1982; Blasi & Colafrancesco, 1999; Dolag & Enßlin, 2000). This process would only be efficient in the dense cluster core regions and thus is likely not important to the formation of peripheral radio relics.

The above scenarios are speculative but make clear predictions which further observations can confirm. First, further observations can reveal that either or both of the relic structures are real. Second, we expect a high degree of source-intrinsic radio polarization in the East relic, due to the nearly perpendicular line of sight to the shock compression direction. The compressed magnetic fields should be roughly aligned with the shock plane and relic major axis, thus the synchrotron emission should be highly polarized. For the West relic, the geometry is less constrained, but polarization is also very likely.

Further, the steep spectral index of the East relic can be understood, but would likely have different spectral curvature, in the two relic formation models. In the Fermi I shock acceleration model a low Mach number shock would produce a steep power-law radio spectrum (Mach number $`<3.5`$ for a spectral index steeper than 1.8). In the revived fossil radio plasma model the steep spectrum would be a synchrotron/inverse Compton cooling cut off, which was shifted by the shock compression close to the observed frequency. Additional information on the spectral curvature of the relics might therefore help to discriminate these two scenarios.

Finally, if radio halos occur only in the shocked ICM – as A754 seems to indicate –, then one would expect a spatial correlation of shock-heated gas and diffuse synchrotron emission in other clusters which are in the early stages of major merger events. Such a correlation could be observable with high-resolution observations from X-ray satellites such as Chandra. At later merger stages, the turbulent gas motion, which is stirred by the violently relaxing dark matter cores, should have erased many of the temperature and radio halo substructures.

The hypothesis that radio halos occur only in shocked ICM could explain why the post merger cluster A3667 has no radio halo, despite the fact that two giant cluster radio relics (Röttgering et al., 1997) indicate the presence of peripheral merger shock waves (Enßlin et al., 1998). In this case the main cluster core seems to be unshocked (Roettiger et al., 1999b, see temperature map in their Fig. 3) because the merger was likely less violent. There only a mass ratio of the merging clusters of 5:1 is needed in order to reproduce its X-ray morphology, compared to 2.5:1 for A754 (Roettiger, Stone & Mushotzky, 1998). In the re-acceleration scenario (Brunetti et al., 2001), the shock waves would have either injected energetic electrons only in peripheral, low magnetic field regions, so that no radio halo occurs, or the necessary re-acceleration turbulence has decayed in the late merger state of A3667 (1 Gyr after core passage). The re-acceleration model can be tested due to the high energy cutoff expected in the radio electron spectrum at energies where the cooling time-scale is shorter than the acceleration time-scale. A stronger halo suppression can be expected in the scenario of radio halo electrons being produced by shock accelerated protons. In this scenario the shock waves in A3367 would have emerged at radii which were too peripheral to host a sufficiently dense thermal proton population required as targets for the hadronic secondary electron production to operate efficiently – in addition to the weaker radio emissivities there due to weaker magnetic fields . The hadronic halo scenario will be tested with the next generation Gamma ray telescopes for the unavoidable Gamma radiation from neutral pion decays following hadronic proton proton collisions (Vestrand, 1982; Enßlin et al., 1997; Colafrancesco & Blasi, 1998; Blasi, 1999; Dolag & Enßlin, 2000).

## 4 Conclusions

A754 has been successfully imaged at both 330 MHz and 74 MHz. We conclude these images strongly suggest the existence of both a radio halo and at least one and possibly two radio relics in this galaxy cluster. This finding is especially relevant in light of recent observational and theoretical evidence in favor of a merger-relic connection. Deeper 74 and 330 MHz imaging at higher resolution are clearly warranted to confirm the extended nature of the radio halo and the relic detections. These observations suggest that further study of this system and similar ones at low frequencies will open a new window of understanding on cluster dynamics and cluster formation.

## 5 Acknowledgments

The authors would like to thank Kurt Roettiger for extensive help with the observations and at early stages of the project. We also thank Rick Perley and Frazer Owen for useful discussions, as well as our anonymous referee whose comments and suggestions improved our paper. NEK also thanks D. Harris and L. Feretti for the Ringberg Castle invitation which inspired the project. We have made use of the ROSAT Data Archive of the Max-Planck-Institut für extraterrestrische Physik (MPE) at Garching, Germany. Basic research in radio astronomy at the Naval Research Laboratory is funded by the Office of Naval Research.


format: text

text:
# I Introduction

## I Introduction

The ability to create, manipulate and characterize quantum states is becoming an increasingly important area of physical research, with implications for areas of technology such as quantum computing, quantum cryptography and communications. With a series of measurements on a large enough number of identically prepared copies of an quantum system, one can infer,to a reasonable approximation, the quantum state of the system. Arguably, the first such experimental technique for determining the state of quantum system was devised by George Stokes in 1852 . His famous four parameters allow an experimenter to determine uniquely the polarization state of a light beam. With the insight provided by nearly 150 years of progress in optical physics, we can consider coherent light beams to be an ensemble of two-level quantum mechanical systems, the two levels being the two polarization degrees of freedom of the photons; the Stokes parameters allow one to determine the density matrix describing this ensemble. More recently, experimental techniques for the measurement of the more subtle quantum properties of light have been the subject of intensive investigation (see ref. for a comprehensive and erudite exposition of this subject). In various experimental circumstances is has been found reasonably straightforward to devise a simple linear tomographic technique in which the density matrix (or Wigner function) of a quantum state is found from a linear transformation of experimental data. However, there is one important drawback to this method, in that the recovered state might not correspond to a physical state because of experimental noise. For example, density matrices for any quantum state must be Hermitian, positive semi-definite matrices with unit trace. The tomographically measured matrices often fail to be positive semi-definite, especially when measuring low-entropy states. To avoid this problem the “maximum likelihood” tomographic approach to the estimation of quantum states has been developed . In this approach the density matrix that is “mostly likely” to have produced a measured data set is determined by numerical optimization.

In the past decade several groups have successfully employed tomographic techniques for the measurement of quantum mechanical systems. In 1990 Risley et al. at North Carolina State University reported the measurement of the density matrix for the nine sublevels of the $`n=3`$ level of hydrogen atoms formed following collision between $`\mathrm{H}^+`$ ions and He atoms, in conditions of high symmetry which simplified the tomographic problem . Since then, in 1993 Smithey et al. (University of Oregon), made a homodyne measurement of the Wigner function of a single mode of light . Other explorations of the quantum states of single mode light fields have been made by Mlynek et al. (University of Konstanz, Germany)  and Bachor et al. (Australian National University) . Other quantum systems whose density matrices have been investigated experimentally include the vibrations of molecules , the motion ions and atoms  and the internal angular momentum quantum state of the $`F=4`$ ground state of a cesium atom . The quantum states of multiple spin-half nuclei have been measured in the high-temperature regime using NMR techniques , albeit in systems of such high entropy that the creation of entangled states is necessarily precluded . The measurement of the quantum state of entangled qubit pairs, realized using the polarization degrees of freedom of a pair of photons created in a parametric down-conversion experiment was reported by us recently .

In this paper we will examine techniques in detail for quantum state measurement as it applies to multiple correlated two-level quantum mechanical systems (or “qubits” in the terminology of quantum information). Our particular emphasis is qubits realized via the two polarization degrees of freedom of photons, data from which we use to illustrate our results. However, these techniques are readily applicable to other technologies proposed for creating entangled states of pairs of two-level systems. Because of the central importance of qubit systems to the emergent discipline of quantum computation, a thorough explanation of the techniques needed to characterize the qubit states will be of relevance to workers in the various diverse experimental fields currently under consideration for quantum computation technology . This paper is organized as follows: In Section II we explore the analogy with the Stokes parameters, and how they lead naturally to a scheme for measurement of an arbitrary number of two-level systems. In Section III, we discuss the measurement of a pair of qubits in more detail, presenting the validity condition for an arbitrary measurement scheme and introducing the set of 16 measurements employed in our experiments. Section IV deals with our method for maximum likelihood reconstruction and in Section V we demonstrate how to calculate the errors in such measurements, and how these errors propagate to quantities calculated from the density matrix.

## II The Stokes parameters and quantum state tomography

As mentioned above, there is a direct analogy between the measurement of the polarization state of a light beam and the measurement of the density matrix of an ensemble of two-level quantum mechanical systems. Here we explore this analogy in more detail.

### A Single qubit tomography

The Stokes parameters are defined from a set of four intensity measurements : (i) with a filter that transmits 50$`\%`$ of the incident radiation, regardless of its polarization; (ii) with a polarizer that transmits only horizontally polarized light; (iii) with a polarizer that transmits only light polarized at $`45^o`$ to the horizontal; and (iv) with a polarizer that transmits only right circularly polarized light. The number of photons counted by a detector, which is proportional to the classical intensity, in these four experiments are as follows:

$`n_0`$ $`=`$ $`{\displaystyle \frac{𝒩}{2}}\left(H|\widehat{\rho }|H+V|\widehat{\rho }|V\right)={\displaystyle \frac{𝒩}{2}}\left(R|\widehat{\rho }|R+L|\widehat{\rho }|L\right)`$ (2.1)
$`n_1`$ $`=`$ $`𝒩\left(H|\widehat{\rho }|H\right)`$ (2.2)
$`=`$ $`{\displaystyle \frac{𝒩}{2}}\left(R|\widehat{\rho }|R+R|\widehat{\rho }|L+L|\widehat{\rho }|R+L|\widehat{\rho }|L\right)`$ (2.3)
$`n_2`$ $`=`$ $`𝒩\left(\overline{D}|\widehat{\rho }|\overline{D}\right)`$ (2.4)
$`=`$ $`{\displaystyle \frac{𝒩}{2}}\left(R|\widehat{\rho }|R+L|\widehat{\rho }|LiL|\widehat{\rho }|R+iR|\widehat{\rho }|L\right)`$ (2.5)
$`n_3`$ $`=`$ $`𝒩\left(R|\widehat{\rho }|R\right)`$ (2.6)
Here $`|H`$, $`|V`$, $`|\overline{D}=\left(|H|V\right)/\sqrt{2}=\mathrm{exp}(i\pi /4)\left(|R+i|L\right)/\sqrt{2}`$ and $`|R=\left(|Hi|V\right)/\sqrt{2}`$ are the kets representing qubits polarized in the linear horizontal, linear vertical, linear diagonal ($`45^o`$) and right-circular senses respectively, $`\widehat{\rho }`$ is the $`(2\times 2)`$ density matrix for the polarization degrees of the light (or, for a two-level quantum system) and $`𝒩`$ is constant dependent on the detector efficiency and light intensity. The Stokes parameters, which fully characterize the polarization state of the light, are then defined by

$`𝒮_0`$ $``$ $`2n_0=𝒩\left(R|\widehat{\rho }|R+L|\widehat{\rho }|L\right)`$ (2.8)
$`𝒮_1`$ $``$ $`2\left(n_1n_0\right)=𝒩\left(R|\widehat{\rho }|L+L|\widehat{\rho }|R\right)`$ (2.9)
$`𝒮_2`$ $``$ $`2\left(n_2n_0\right)=𝒩i\left(R|\widehat{\rho }|LL|\widehat{\rho }|R\right)`$ (2.10)
$`𝒮_3`$ $``$ $`2\left(n_3n_0\right)=𝒩\left(R|\widehat{\rho }|RL|\widehat{\rho }|L\right).`$ (2.11)
We can now relate the Stokes parameters to the density matrix $`\widehat{\rho }`$ by the formula

$$\widehat{\rho }=\frac{1}{2}\underset{i=0}{\overset{3}{}}\frac{𝒮_i}{𝒮_0}\widehat{\sigma }_i,$$
(2.13)
where $`\widehat{\sigma }_0=|RR|+|LL|`$ is the single qubit identity operator and $`\widehat{\sigma }_1=|RL|+|LR|`$, $`\widehat{\sigma }_2=i(|LR||RL|`$ and $`\widehat{\sigma }_3=|RR||LL|`$ are the Pauli spin operators. Thus the measurement of the Stokes parameters can be considered equivalent to a tomographic measurement of the density matrix of an ensemble of single qubits.

### B Multiple beam Stokes parameters: multiple qubit tomography

The generalization of the Stokes scheme to measure the state of multiple photon beams (or multiple qubits) is reasonably straightforward. One should, however, be aware that importance differences exist between the one-photon and the multiple photon cases. Single photons, at least in the current context, can be described in a purely classical manner, and the density matrix can be related to the purely classical concept of the coherency matrix . For multiple photon one has the possibility of non-classical correlations occurring, with quintessentially quantum-mechanical phenomena such as entanglement being present. We will return to the concept of entanglement and how it may be measured later in this paper.

An n-qubit state is characterized by a density matrix which may be written as follows:

$$\widehat{\rho }=\frac{1}{2^n}\underset{i_1,i_2,\mathrm{}i_n=0}{\overset{3}{}}r_{i_1,i_2,\mathrm{}i_n}\widehat{\sigma }_{i_1}\widehat{\sigma }_{i_2}\mathrm{}\widehat{\sigma }_{i_n},$$
(2.14)
where the $`4^n`$ parameters $`r_{i_1,i_2,\mathrm{}i_n}`$ are real numbers. The normalization property of the density matrices requires that $`r_{0,0,\mathrm{}0}=1`$, and so the density matrix is specified by $`4^n1`$ real parameters. The symbol $``$ represents the tensor product between operators acting on the Hilbert spaces associated with the separate qubits.

As Stokes showed, the state of a single qubit can be determined by taking a set of four projection measurements which are represented by the four operators $`\widehat{\mu }_0=|HH|+|VV|`$, $`\widehat{\mu }_1=|HH|`$, $`\widehat{\mu }_2=|\overline{D}\overline{D}|`$, $`\widehat{\mu }_3=|RR|`$. Similarly, the state of two qubits can be determined by the set of 16 measurements represented by the operators $`\widehat{\mu }_i\widehat{\mu }_j`$ ($`i,j=0,1,2,3`$). More generally the state of an n-qubit system can be determined by $`4^n`$ measurements given by the operators $`\widehat{\mu }_{i_1}\widehat{\mu }_{i_2}\mathrm{}\widehat{\mu }_{i_n}`$ ($`i_k=0,1,2,3`$) and ($`k=1,2,\mathrm{}n`$). This ‘tree’ structure for multi-qubit measurement is illustrated in fig.1.

The proof of this conjecture is reasonably straightforward. The outcome of a measurement is given by the formula

$$n=𝒩\mathrm{Tr}\left\{\widehat{\rho }\widehat{\mu }\right\},$$
(2.15)
where $`\widehat{\rho }`$ is the density matrix, $`\widehat{\mu }`$ is the measurement operator and $`𝒩`$ is a constant of proportionality which can be determined from the data. Thus in our n-qubit case, the outcomes of the various measurement are

$$n_{i_1,i_2,\mathrm{}i_n}=𝒩\mathrm{Tr}\left\{\widehat{\rho }\left(\widehat{\mu }_{i_1}\widehat{\mu }_{i_2}\mathrm{}\widehat{\mu }_{i_n}\right)\right\}.$$
(2.16)
Substituting from eq.(2.14) we obtain

$$n_{i_1,i_2,\mathrm{}i_n}=\frac{𝒩}{2^n}\underset{j_1,j_2,\mathrm{}j_n=0}{\overset{3}{}}\mathrm{Tr}\left\{\widehat{\mu }_{i_1}\widehat{\sigma }_{j_1}\right\}\mathrm{Tr}\left\{\widehat{\mu }_{i_2}\widehat{\sigma }_{j_2}\right\}\mathrm{}\mathrm{Tr}\left\{\widehat{\mu }_{i_n}\widehat{\sigma }_{j_n}\right\}r_{i_1,i_2,\mathrm{}i_n}.$$
(2.17)
As can be easily verified, the single qubit measurement operators $`\widehat{\mu }_i`$ are linear combinations of the Pauli operators $`\widehat{\sigma }_j`$, i.e. $`\widehat{\mu }_i=_{j=0}^3\mathrm{{\rm Y}}_{ij}\widehat{\sigma }_j`$, where $`\mathrm{{\rm Y}}_{ij}`$ are the elements of the matrix

$$\mathrm{{\rm Y}}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 1/2& 1/2& 0& 0\\ 1/2& 0& 1/2& 0\\ 1/2& 0& 0& 1/2\end{array}\right).$$
(2.18)
Further, we have the relation $`\mathrm{Tr}\left\{\widehat{\sigma }_i\widehat{\sigma }_j\right\}=2\delta _{ij}`$ (where $`\delta _{ij}`$ is the Kronecker delta). Hence eq.(2.17) becomes

$$n_{i_1,i_2,\mathrm{}i_n}=𝒩\underset{j_1,j_2,\mathrm{}j_n=0}{\overset{3}{}}\mathrm{{\rm Y}}_{i_1j_1}\mathrm{{\rm Y}}_{i_2j_2}\mathrm{}\mathrm{{\rm Y}}_{i_nj_n}r_{i_1,i_2,\mathrm{}i_n}.$$
(2.19)
Introducing the left-inverse of the matrix $`\mathrm{{\rm Y}}`$, defined so that $`_{k=0}^3(\mathrm{{\rm Y}}^1)_{ik}\mathrm{{\rm Y}}_{kj}=\delta _{ij}`$ and whose elements are

$$\mathrm{{\rm Y}}^1=\left(\begin{array}{cccc}1& 0& 0& 0\\ 1& 2& 0& 0\\ 1& 0& 2& 0\\ 1& 0& 0& 2\end{array}\right),$$
(2.20)
we can find a formula for the parameters $`r_{i_1,i_2,\mathrm{}i_n}`$ in terms of the measured quantities $`n_{i_1,i_2,\mathrm{}i_n}`$, viz.:

$`𝒩r_{i_1,i_2,\mathrm{}i_n}`$ $`=`$ $`{\displaystyle \underset{j_1,j_2,\mathrm{}j_n=0}{\overset{3}{}}}\left(\mathrm{{\rm Y}}^1\right)_{i_1j_1}\left(\mathrm{{\rm Y}}^1\right)_{i_2j_2}\mathrm{}\left(\mathrm{{\rm Y}}^1\right)_{i_nj_n}n_{i_1,i_2,\mathrm{}i_n}`$ (2.21)
$``$ $`𝒮_{i_1,i_2,\mathrm{}i_n}.`$ (2.22)
In eq.(2.22) we have introduced the n-photon Stokes parameter $`𝒮_{i_1,i_2,\mathrm{}i_n}`$, defined in an analogous manner to the single photon Stokes parameters give in eq.(LABEL:repulse).

Since, as already noted, $`r_{0,0,\mathrm{}0}=1`$, we can make the identification $`𝒮_{0,0,\mathrm{}0}=𝒩`$, and so the density matrix for the n-qubit system can be written in terms of the Stokes parameters as follows:

$$\widehat{\rho }=\frac{1}{2^n}\underset{i_1,i_2,\mathrm{}i_n=0}{\overset{3}{}}\frac{𝒮_{i_1,i_2,\mathrm{}i_n}}{𝒮_{0,0,\mathrm{}0}}\widehat{\sigma }_{i_1}\widehat{\sigma }_{i_2}\mathrm{}\widehat{\sigma }_{i_n}.$$
(2.23)
This is a recipe for measurement of the density matrices which, assuming perfect experimental conditions and the complete absence of noise, will always work. It is important to realize that the set of four Stokes measurements $`\{\widehat{\mu }_0,\widehat{\mu }_1,\widehat{\mu }_2,\widehat{\mu }_3\}`$ are not unique: there may be circumstances in which it is more convenient to use some other set, which are equivalent. A more typical set, at least in optical experiments, is $`\widehat{\mu }_0^{}=|HH|`$, $`\widehat{\mu }_1^{}=|VV|`$, $`\widehat{\mu }_2^{}=|DD|`$, $`\widehat{\mu }_3^{}=|RR|`$.

In the following section we will explore more general schemes for the measurement of two qubits, starting with a discussion, in some detail, of how the measurements are actually performed.

## III Generalized Tomographic Reconstruction of the Polarization State of Two Photons

### A Experimental set-up

The experimental arrangement used in our experiments is shown schematically in Fig.1. An optical system consisting of lasers, polarization elements and non-linear optical crystals (and collectively characterized for the purposes of this paper as a “black-box”,) is used to generate pairs of qubits in an almost arbitrary quantum state of their polarization degrees of freedom. A full description of this optical system and how such quantum states can be prepared can be found in ref. It is important to realize that the entangled photon pairs are produced in a non-deterministic manner: one cannot specify with certainly when a photon pair will be emitted; indeed there is a small probability of generating four, or six or higher numbers of photons. Thus we can only post-selectively generate entangled photon pairs: i.e. one only knows that the state was created after if has been measured. The output of the black box consists of a pair of beams of light, whose quanta can be measured by means of photo-detectors. To project the light beams onto a polarization state of the experimenter’s choosing, three optical elements are placed in the beam in front of each detector: a polarizer (which transmits only vertically polarized light), a quarter-wave plate and a half-wave plate. The angles of the fast axes of both of the waveplates can be set arbitrarily, allowing the $`|V`$ projection state fixed by the polarizer to be rotated into any polarization state that the experimenter may wish.

Using the Jones calculus notation, with the following convention,

$$\left(\begin{array}{c}0\\ 1\end{array}\right)=|\mathrm{V},\left(\begin{array}{c}1\\ 0\end{array}\right)=|\mathrm{H},$$
(3.1)
where $`|\mathrm{V}`$ ($`|\mathrm{H}`$) is the ket for a vertically (horizontally) polarized beam, the effect of quarter and half wave plates whose fast axes are at angles $`q`$ and $`h`$ with respect to the vertical axis, respectively, are given by the 2 $`\times `$ 2 matrices

$`\widehat{U}_{QWP}(q)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cc}i\mathrm{cos}(2q)& \mathrm{sin}(2q)\\ \mathrm{sin}(2q)& i+\mathrm{cos}(2q)\end{array}\right),`$ (3.4)
$`\widehat{U}_{HWP}(q)`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{cos}(2h)& \mathrm{sin}(2h)\\ \mathrm{sin}(2h)& \mathrm{cos}(2h)\end{array}\right).`$ (3.7)
Thus the projection state for the measurement in one of the beams is given by

$`|\psi _{proj}^{(1)}(h,q)`$ $`=`$ $`\widehat{U}_{QWP}(q)\widehat{U}_{HWP}(h)\left(\begin{array}{c}0\\ 1\end{array}\right)`$ (3.10)
$`=`$ $`a(h,q)|\mathrm{H}+b(h,q)|\mathrm{V},`$ (3.11)
where, neglecting an overall phase, the functions $`a(h,q)`$ and $`b(h,q)`$ are given by

$`a(h,q)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\mathrm{sin}(2h)i\mathrm{sin}[2(hq)]\right),`$ (3.12)
$`b(h,q)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\mathrm{cos}(2h)+i\mathrm{cos}[2(hq)]\right).`$ (3.13)
The projection state for the two beams is given by

$`|\psi _{proj}^{(2)}(h_1,q_1,h_2,q_2)`$ $`=`$ $`|\psi _{proj}^{(1)}(h_1,q_1)|\psi _{proj}^{(1)}(h_2,q_2)`$ (3.14)
$`=`$ $`a(h_1,q_1)a(h_2,q_2)|\mathrm{HH}+a(h_1,q_1)b(h_2,q_2)|\mathrm{HV}+`$ (3.16)
$`b(h_1,q_1)a(h_2,q_2)|\mathrm{VH}+b(h_1,q_1)b(h_2,q_2)|\mathrm{VV}.`$
We shall denote the projection state corresponding to one particular set of waveplate angles $`\{h_{1,\nu },q_{1,\nu },h_{2,\nu },q_{2,\nu }\}`$ Here the first subscript on the waveplate angle refers one of the two photon beams; the second subscript distinguishes which of the sixteen different experimental states is under consideration. by the ket $`|\psi _\nu `$; thus the projection measurement is represented by the operator $`\widehat{\mu }_\nu =|\psi _\nu \psi _\nu |`$. Consequently, the average number of coincidence counts that will be observed in a given experimental run is

$$n_\nu =𝒩\psi _\nu |\widehat{\rho }|\psi _\nu $$
(3.17)
where $`\widehat{\rho }`$ is the density matrix describing the ensemble of qubits, and $`𝒩`$ is a constant dependent on the photon flux and detector efficiencies. In what follows, it will be convenient to consider the quantities $`s_\nu `$ defined by

$$s_\nu =\psi _\nu |\widehat{\rho }|\psi _\nu .$$
(3.18)
### B Tomographically Complete set of Measurements

In Section II we have given one possible set of projection measurements $`\{|\psi _\nu \psi _\nu |\}`$ which which uniquely determine the density matrix $`\widehat{\rho }`$. However, one can conceive of situations in which these will not be the most convenient set of measurements to make. Here we address the problem of finding other sets of suitable measurements. The smallest number of states required for such measurements can be found by a simple argument: there are 15 real unknown parameters which determine a $`4\times 4`$ density matrix, plus there is the single unknown real parameter $`𝒩`$, making a total of 16.

In order to proceed it is helpful to convert the $`4\times 4`$ matrix $`\widehat{\rho }`$ into a 16-dimensional column vector. To do this we use a set of 16 linearly independent $`4\times 4`$ matrices $`\{\widehat{\mathrm{\Gamma }}_\nu \}`$ which have the following mathematical properties:

$`\mathrm{Tr}\left\{\widehat{\mathrm{\Gamma }}_\nu \widehat{\mathrm{\Gamma }}_\mu \right\}`$ $`=`$ $`\delta _{\nu ,\mu }`$ (3.19)
$`\widehat{A}`$ $`=`$ $`{\displaystyle \underset{\nu =1}{\overset{16}{}}}\widehat{\mathrm{\Gamma }}_\nu \mathrm{Tr}\left\{\widehat{\mathrm{\Gamma }}_\nu \widehat{A}\right\}\widehat{A},`$ (3.20)
where $`\widehat{A}`$ is an arbitrary $`4\times 4`$ matrix. Finding a set of $`\widehat{\mathrm{\Gamma }}_\nu `$ matrices is in fact reasonably straightforward: for example, the set of (appropriately normalized) generators of the Lie algebra $`SU(2)SU(2)`$ fulfill the required criteria (for reference, we list this set in Appendix A). These matrices are of course simply a re-labeling of the two-qubit Pauli matrices $`\widehat{\sigma }_i\widehat{\sigma }_j`$ ($`i,j=0,1,2,3`$) discussed above. Using these matrices the density matrix can be written as

$$\widehat{\rho }=\underset{\nu =1}{\overset{16}{}}\widehat{\mathrm{\Gamma }}_\nu r_\nu ,$$
(3.21)
where $`r_\nu `$ is the $`\nu `$-th element of a sixteen element column vector, given by the formula

$$r_\nu =\mathrm{Tr}\left\{\widehat{\mathrm{\Gamma }}_\nu \widehat{\rho }\right\}$$
(3.22)
Substituting from eq.(3.21) into eq.(3.17), we obtain the following linear relationship between the measured coincidence counts $`n_\nu `$ and the elements of the vector $`r_\mu `$:

$$n_\nu =𝒩\underset{\mu =1}{\overset{16}{}}B_{\nu ,\mu }r_\mu $$
(3.23)
where the $`16\times 16`$ matrix $`B_{\nu ,\mu }`$ is given by

$$B_{\nu ,\mu }=\psi _\nu |\widehat{\mathrm{\Gamma }}_\mu |\psi _\nu .$$
(3.24)
Immediately we find a necessary and sufficient condition for the completeness of the set of tomographic states $`\{|\psi _\nu \}`$: if the matrix $`B_{\nu ,\mu }`$ is non- singular, then eq.(3.23) can be inverted to give

$$r_\nu =(𝒩)^1\underset{\mu =1}{\overset{16}{}}\left(B^1\right)_{\nu ,\mu }n_\mu .$$
(3.25)
The set of sixteen tomographic states which we employed are given in Table 1. They can be shown to satisfy the condition that $`B_{\nu ,\mu }`$ is non- singular. By no means are these states unique in this regard: these were the states chosen principally for experimental convenience.

These states can be realized by setting specific values of the half- and quarter-wave plate angles. The appropriate values of these angles (measured from the vertical) are given in Table 1. Note that overall phase factors do not affect the results of projection measurements.

Substituting eq.(3.25) into eq.(3.21), we find that

$`\widehat{\rho }`$ $`=`$ $`(𝒩)^1{\displaystyle \underset{\nu =1}{\overset{16}{}}}\widehat{M}_\nu n_\nu `$ (3.26)
$`=`$ $`{\displaystyle \underset{\nu =1}{\overset{16}{}}}\widehat{M}_\nu s_\nu ,`$ (3.27)
where the sixteen $`4\times 4`$ matrices $`\widehat{M}_\nu `$ are defined by

$$\widehat{M}_\nu =\underset{\nu =1}{\overset{16}{}}\left(B^1\right)_{\nu ,\mu }\widehat{\mathrm{\Gamma }}_\mu .$$
(3.28)
The introduction of the $`\widehat{M}_\nu `$ matrices allows a compact form of linear tomographic reconstruction, eq.(3.27), will be most useful in the error analysis that follows. These $`\widehat{M}_\nu `$ matrices, valid for our set of tomographic states, are listed in Appendix B, together with some of their important properties. We can use one of these properties, eq.(B.6), to obtain the value of the unknown quantity $`𝒩`$. That relationship implies

$$\underset{\nu }{}\mathrm{Tr}\left\{\widehat{M}_\nu \right\}|\psi _\nu \psi _\nu |\widehat{\rho }=\widehat{\rho }.$$
(3.29)
Taking the trace of this formula, and multiplying by $`𝒩`$ we obtain:

$$\underset{\nu }{}\mathrm{Tr}\left\{\widehat{M}_\nu \right\}n_\nu =𝒩.$$
(3.30)
For our set of tomographic states, it can be shown that

$$\underset{\nu }{}\mathrm{Tr}\left\{\widehat{M}_\nu \right\}=\{\begin{array}{cc}1\hfill & \text{if }\nu =1,2,3,4\hfill \\ 0\hfill & \text{if }\nu =5,\mathrm{}16,\hfill \end{array}$$
(3.31)
hence the value of the unknown parameter $`𝒩`$ in our experiments is given by:

$`𝒩`$ $`=`$ $`{\displaystyle \underset{\nu =1}{\overset{4}{}}}n_\nu `$ (3.32)
$`=`$ $`𝒩\left(HH|\widehat{\rho }|HH+HV|\widehat{\rho }|HV+VH|\widehat{\rho }|VH+VV|\widehat{\rho }|VV\right).`$ (3.33)
Thus we obtain the final formula for the tomographic reconstruction of the density matrices of our states:

$$\widehat{\rho }=\left(\underset{\nu =1}{\overset{16}{}}\widehat{M}_\nu n_\nu \right)/\left(\underset{\nu =1}{\overset{4}{}}n_\nu \right).$$
(3.34)
As an example, the following set of 16 counts were taken for the purpose of tomographically determining the density matrix for an ensemble of qubits all prepared in a specific quantum state: $`n_1=34749,n_2=324,n_3=35805,n_4=444,n_5=16324,n_6=17521,n_7=13441,n_8=16901,n_9=17932,n_{10}=32028,n_{11}=15132,n_{12}=17238,n_{13}=13171,n_{14}=17170,n_{15}=16722,n_{16}=33586`$ Applying eq.(3.34) we find

$$\widehat{\rho }=\left(\begin{array}{cccc}0.4872& 0.0042+i0.0114& 0.0098i0.0178& 0.5192+i0.0380\\ 0.0042i0.0114& 0.0045& 0.0271i0.0146& 0.0648i0.0076\\ 0.0098+i0.0178& 0.0271+i0.0146& 0.0062& 0.0695+i0.0134\\ 0.5192i0.0380& 0.0648+i0.0076& 0.0695i0.0134& 0.5020\end{array}\right)$$
(3.35)
This matrix is shown graphically in fig.3a.

Note that, by construction, the density matrix is normalized, i.e. $`\mathrm{Tr}\{\widehat{\rho }\}=1`$ and Hermitian, i.e. $`\widehat{\rho }^{}=\widehat{\rho }`$ . However, when one calculates the eigenvalues of this measured density matrix, one finds the values 1.02155, 0.0681238, -0.065274 and -0.024396; and also, $`\mathrm{Tr}\{\widehat{\rho }^2\}=1.053`$ . Density matrices for all physical states must have the property of positive semi-definiteness, which (in conjunction with the normalization and Hermiticity properties) imply that all of the eigenvalues must lie in the interval $`[0,1]`$, their sum being 1; this in turn implies that $`0\mathrm{Tr}\{\widehat{\rho }^2\}1`$. Clearly the density matrix reconstructed above by linear tomography violates these condition. From our experience of tomographic measurements of various mixed and entangled states prepared experimentally, this seems to happen roughly 75$`\%`$ of the time for low entropy, highly entangled states; it seems to have a higher probability of producing the correct result for states of higher entropy, but the cautious experimenter should check every time. The obvious culprit for this problem is experimental inaccuracies and statistical fluctuations of coincidence counts, which mean that the actual numbers of counts recorded in a real experiment differ from those that can be calculated by eq.(3.17). Thus the linear reconstruction is of limited value for states of low entropy (which are of most experimental interest because of their application to quantum information technology); however as we shall see, the linear approach does provide a useful starting point for the numerical optimization approach to density matrix estimation which we will discuss in the next section.

## IV Maximum Likelihood Estimation

As mentioned in Section III, the tomographic measurement of density matrices can produce results which violate important basic properties such as positivity. To avoid this problem, the maximum likelihood estimation of density matrices may be employed. Here we describe a simple realization of this technique.

### A Basic approach

Our approach to the maximum likelihood estimation of the density matrix is as follows:

(i) Generate a formula for an explicitly “physical” density matrix, i.e. a matrix which has the three important properties of normalization, Hermiticity and positivity. This matrix will be a function of 16 real variables (denoted $`\{t_1,t_2,\mathrm{}t_{16}\}`$). We will denote the matrix as $`\widehat{\rho }_p(t_1,t_2,\mathrm{}t_{16})`$.

(ii) Introduce a “likelihood function” which quantifies how good the density matrix $`\widehat{\rho }_p(t_1,t_2,\mathrm{}t_{16})`$ is in relation to the experimental data. This likelihood function is a function of the 16 real parameters $`t_\nu `$ and of the 16 experimental data $`n_\nu `$. We will denote this function as $`(t_1,t_2,\mathrm{}t_{16};n_1,n_2,\mathrm{}n_{16})`$.

(iii) Using standard numerical optimization techniques, find the optimum set of variables $`\{t_1^{(opt)},t_2^{(opt)},\mathrm{}t_{16}^{(opt)}\}`$ for which the function $`(t_1,t_2,\mathrm{}t_{16};n_1,n_2,\mathrm{}n_{16})`$ has its maximum value. The best estimate for the density matrix is then $`\widehat{\rho }(t_1^{(opt)},t_2^{(opt)},\mathrm{}t_{16}^{(opt)})`$.

The details of how these three steps can be carried out are described in the next three sub-sections.

### B Physical Density Matrices

The property of non-negative definiteness for any matrix $`\widehat{𝒢}`$ is written mathematically as

$$\psi |\widehat{𝒢}|\psi 0|\psi .$$
(4.1)
Any matrix that can be written in the form $`\widehat{𝒢}=\widehat{T}^{}\widehat{T}`$ must be non-negative definite. To see that this is the case, substitute into eq.(4.1):

$$\psi |\widehat{T}^{}\widehat{T}|\psi =\psi ^{}|\psi ^{}0,$$
(4.2)
where we have defined $`|\psi ^{}=\widehat{T}|\psi `$. Furthermore $`(\widehat{T}^{}\widehat{T})^{}=\widehat{T}^{}(\widehat{T}^{})^{}=\widehat{T}^{}\widehat{T}`$, i.e. $`\widehat{𝒢}=\widehat{T}^{}\widehat{T}`$ must be Hermitian. To ensure normalization, one can simply divide by the trace: thus the matrix $`\widehat{g}`$ given by the formula

$$\widehat{g}=\widehat{T}^{}\widehat{T}/\mathrm{Tr}\{\widehat{T}^{}\widehat{T}\}$$
(4.3)
has all three of the mathematical properties which we require for density matrices.

For the two qubit system, we have a $`4\times 4`$ density matrix with 15 independent real parameters. Since it will be useful to be able to invert relation (4.3), it is convenient to choose a tri-diagonal form for $`\widehat{T}`$:

$$\widehat{T}(t)=\left(\begin{array}{cccc}t_1& 0& 0& 0\\ t_5+it_6& t_2& 0& 0\\ t_{11}+it_{12}& t_7+it_8& t_3& 0\\ t_{15}+it_{16}& t_{13}+it_{14}& t_9+it_{10}& t_4\end{array}\right),$$
(4.4)
Thus the explicitly “physical” density matrix $`\widehat{\rho }_p`$ is given by the formula

$$\widehat{\rho }_p(t)=\widehat{T}^{}(t)\widehat{T}(t)/\mathrm{Tr}\{\widehat{T}^{}(t)\widehat{T}(t)\}.$$
(4.5)
For future reference, the inverse relationship, by which the elements of $`\widehat{T}`$ can be expressed in terms of the elements of $`\widehat{\rho }`$, is as follows:

$$\widehat{T}=\left(\begin{array}{cccc}\sqrt{\frac{\mathrm{\Delta }}{_{11}^{(1)}}}& 0& 0& 0\\ & & & \\ \frac{_{12}^{(1)}}{\sqrt{_{11}^{(1)}_{11,22}^{(2)}}}& \sqrt{\frac{_{11}^{(1)}}{_{11,22}^{(2)}}}& 0& 0\\ & & & \\ \frac{_{12,23}^{(2)}}{\sqrt{\rho _{44}}\sqrt{_{11,22}^{(2)}}}& \frac{_{11,23}^{(2)}}{\sqrt{\rho _{44}}\sqrt{_{11,22}^{(2)}}}& \sqrt{\frac{_{11,22}^{(2)}}{\rho _{44}}}& 0\\ & & & \\ \frac{\rho _{41}}{\sqrt{\rho _{44}}}& \frac{\rho _{42}}{\sqrt{\rho _{44}}}& \frac{\rho _{43}}{\sqrt{\rho _{44}}}& \sqrt{\rho _{44}}\end{array}\right).$$
(4.6)
Here we have used the notation $`\mathrm{\Delta }=\mathrm{Det}(\widehat{\rho })`$; $`_{ij}^{(1)}`$ is the first minor of $`\widehat{\rho }`$, i.e. the determinant of the $`3\times 3`$ matrix formed by deleting the $`i`$-th row and $`j`$-th column of $`\widehat{\rho }`$; $`_{ij,kl}^{(2)}`$ is the second minor of $`\widehat{\rho }`$, i.e. the determinant of the $`2\times 2`$ matrix formed by deleting the $`i`$-th and $`k`$-th rows and $`j`$-th and $`l`$-th columns of $`\widehat{\rho }`$ ($`ik`$ and $`jl`$).

### C The Likelihood Function

The measurement data consists of a set of 16 coincidence counts $`n_\nu (\nu =1,2,\mathrm{}16)`$ whose expected value is $`\overline{n}_\nu =𝒩\psi _\nu |\widehat{\rho }|\psi _\nu `$. Let us assume that the noise on these coincidence measurements has a Gaussian probability distribution. Thus the probability of obtaining a set of 16 counts $`\{n_1,n_2,\mathrm{}N_{16}\}`$ is

$$P(n_1,n_2,\mathrm{}n_{16})=\frac{1}{\mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}}\underset{\nu =1}{\overset{16}{}}\mathrm{exp}\left[\frac{(n_\nu \overline{n}_\nu )^2}{2\sigma _\nu ^2}\right],$$
(4.7)
where $`\sigma _\nu `$ is the standard deviation for $`\nu `$-th coincidence measurement (given approximately by $`\sqrt{\overline{n}_\nu }`$) and $`\mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}`$ is the normalization constant. For our physical density matrix $`\widehat{\rho }_p`$ the number of counts expected for the $`\nu `$-th measurement is

$$\overline{n}_\nu (t_1,t_2,\mathrm{}t_{16})=𝒩\psi _\nu |\widehat{\rho }_p(t_1,t_2,\mathrm{}t_{16})|\psi _\nu .$$
(4.8)
Thus the likelihood that the matrix $`\widehat{\rho }_p(t_1,t_2,\mathrm{}t_{16})`$ could produce the measured data $`\{n_1,n_2,\mathrm{}n_{16}\}`$ is

$$P(t_1,t_2,\mathrm{}t_{16})=\frac{1}{\mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}}\underset{\nu =1}{\overset{16}{}}\mathrm{exp}\left[\frac{(𝒩\psi _\nu |\widehat{\rho }_p(t_1,t_2,\mathrm{}t_{16})|\psi _\nu n_\nu )^2}{2𝒩\psi _\nu |\widehat{\rho }_p(t_1,t_2,\mathrm{}t_{16})|\psi _\nu n_\nu )}\right],$$
(4.9)
where $`𝒩=_{\nu =1}^4N_\nu `$.

Rather than find maximum value of $`P(t_1,t_2,\mathrm{}t_{16})`$ it simplifies things somewhat to find the maximum of its logarithm (which is mathematically equivalent) <sup>§</sup><sup>§</sup>§Note that here we neglect the dependence of the normalization constant on $`t_1,t_2,\mathrm{}t_{16}`$, which only weakly effects solution for the most likely state. Thus the optimization problem reduces to finding the minimum of the following function:

$$(t_1,t_2,\mathrm{}t_{16})=\underset{\nu =1}{\overset{16}{}}\frac{(𝒩\psi _\nu |\widehat{\rho }_p(t_1,t_2,\mathrm{}t_{16})|\psi _\nu n_\nu )^2}{2𝒩\psi _\nu |\widehat{\rho }_p(t_1,t_2,\mathrm{}t_{16})|\psi _\nu n_\nu )}.$$
(4.10)
This is the “likelihood” function which we employed in our numerical optimization routine.

### D Numerical Optimization

We used the Mathematica 4.0 routine FindMinimum which executes a multidimensional Powell direction set algorithm (see ref. for a description of this algorithm). To execute this routine, one requires an initial estimate for the values of $`t_1,t_2,\mathrm{}t_{16}`$. For this, we used the tomographic estimate of the density matrix in the inverse relation (4.6), allowing us to determine a set of values for $`t_1,t_2,\mathrm{}t_{16}`$. Since the tomographic density matrix may not be non-negative definite, the values of the $`t_\nu `$’s deduced in this manner are not necessarily real. Thus for our initial guess we used the real parts of the $`t_\nu `$’s deduced from the tomographic density matrix.

For the example given in Section 2, the maximum likelihood estimate is

$$\widehat{\rho }=\left(\begin{array}{cccc}0.5069& 0.0239+i0.0106& 0.0412i0.0221& 0.4833+i0.0329\\ 0.0239i0.0106& 0.0048& 0.0023+i0.0019& 0.0296i0.0077\\ 0.0412+i0.0221& 0.0023i0.0019& 0.0045& 0.0425+i0.0192\\ 0.4833i0.0329& 0.0296+i0.0077& 0.0425i0.0192& 0.4839\end{array}\right).$$
(4.11)
This matrix is illustrated in Fig.2b. In this case, the matrix has eigenvalues 0.986022, 0.0139777, 0 and 0; and $`\mathrm{Tr}\{\widehat{\rho }^2\}=0.972435`$, indicating, while the linear reconstruction gave a non-physical density matrix, the maximum-likelihood reconstruction gives a legitimate density matrix.

## V Error Analysis

In this section we present an analysis of the errors inherent in the tomographic scheme described in Section III. Two sources of errors are found to be important: the shot noise error in the measured coincidence counts $`n_\nu `$ and the uncertainty in the settings of the angles of the waveplates used to make the tomographic projection states. We will analyze these two sources separately.

In addition to determining the density matrix of a pair of qubits, one is often also interested in quantities derived from the density matrix, such as the entropy or the entanglement of formation. For completeness, we will also derive the errors in some of these quantities.

### A Errors due to Count Statistics

From eq.(3.34) we see that the density matrix is specified by a set of sixteen parameters $`s_\nu `$ defined by

$$s_\nu =n_\nu /𝒩,$$
(5.1)
where $`n_\nu `$ are the measured coincidence counts and $`𝒩=_{\nu =1}^4n_\nu `$. We can determine the errors in $`s_\nu `$ using the following formula

$$\overline{\delta s_\nu \delta s_\mu }=\underset{\lambda ,\kappa =1}{\overset{16}{}}\left(\frac{s_\nu }{n_\lambda }\right)\left(\frac{s_\mu }{n_\kappa }\right)\overline{\delta n_\lambda \delta n_\kappa },$$
(5.2)
where the over-bar denotes the ensemble average of the random uncertainties $`\delta s_\nu `$ and $`\delta n_\lambda `$. The measured coincidence counts $`n_\lambda `$ are statistically independent Poissonian random variables, which implies the following relation:

$$\overline{\delta n_\lambda \delta n_\kappa }=n_\lambda \delta _{\lambda ,\kappa },$$
(5.3)
where $`\delta _{\lambda ,\kappa }`$ is the Kronecker delta.

Taking the derivative of eq.(5.1), we find that

$$\frac{s_\mu }{n_\nu }=\frac{1}{𝒩}\delta _{\mu \nu }\frac{n_\mu }{𝒩^2}D_\nu ,$$
(5.4)
where

$$D_\nu =\underset{\lambda =1}{\overset{4}{}}\delta _{\lambda ,\nu }=\{\begin{array}{cc}1\hfill & \text{if}\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1}\nu 4\hfill \\ & \\ 0\hfill & \text{if}\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}5}\nu 16.\hfill \end{array}$$
(5.5)
Substituting from eq.(5.4) into eq.(5.2) and using eq.(5.3), we obtain the result

$$\overline{\delta s_\nu \delta s_\mu }=\frac{n_\mu }{𝒩^2}\delta _{\nu ,\mu }+\frac{n_\nu n_\mu }{𝒩^3}(1D_\mu D_\nu ).$$
(5.6)
In most experimental circumstances $`𝒩1`$, and so the second term in eq.(5.6) is negligibly small in comparison to the first. We shall therefore ignore it, and use the approximate expression in the subsequent discussion;

$$\overline{\delta s_\nu \delta s_\mu }\frac{n_\mu }{𝒩^2}\delta _{\nu ,\mu }\frac{s_\mu }{𝒩}\delta _{\nu ,\mu }.$$
(5.7)
### B Errors due to Angular Settings Uncertainties

Using the formula (3.18) for the parameters $`s_\nu `$ we can find the dependence of the measured density matrix on errors in the tomographic states. The derivative of $`s_\nu `$ with respect to some generic wave-plate settings angle $`\theta `$ is

$$\frac{s_\nu }{\theta }=\left\{\frac{}{\theta }\psi _\nu |\right\}\widehat{\rho }|\psi _\nu +\psi _\nu |\widehat{\rho }\left\{\frac{}{\theta }|\psi _\nu \right\},$$
(5.8)
where $`|\psi _\nu `$ is the ket of the $`\nu `$-th projection state \[see eq.(3.16)\]. Substituting from eq.(3.27) we find

$$\frac{s_\nu }{\theta }=\underset{\mu =1}{\overset{16}{}}s_\mu \left[\left\{\frac{}{\theta }\psi _\nu |\right\}\widehat{M_\mu }|\psi _\nu +\psi _\nu |\widehat{M_\mu }\left\{\frac{}{\theta }|\psi _\nu \right\}\right].$$
(5.9)
For convenience, we shall label the four waveplate angles $`\{h_{1,\nu },q_{1,\nu },h_{2,\nu },q_{2,\nu }\}`$ which specify the $`\nu `$-th state by $`\{\theta _{\nu ,1},\theta _{\nu ,2},\theta _{\nu ,3},\theta _{\nu ,4}\}`$ respectively. Clearly the $`\mu `$-th state does not depend on any of the $`\nu `$-th set of angles. Thus we obtain the following expression for the derivatives of $`s_\nu `$ with respect to waveplate settings:

$$\frac{s_\nu }{\theta _{\lambda ,i}}=\delta _{\nu ,\lambda }\underset{\mu =1}{\overset{16}{}}s_\mu f_{\nu ,\mu }^{(i)},$$
(5.10)
where

$$f_{\nu ,\mu }^{(i)}=\left\{\frac{}{\theta _{\nu ,i}}\psi _\nu |\right\}\widehat{M_\mu }|\psi _\nu +\psi _\nu |\widehat{M_\mu }\left\{\frac{}{\theta _{\nu ,i}}|\psi _\nu \right\}.$$
(5.11)
The 1024 quantities $`f_{\nu ,\mu }^{(i)}`$ can be determined by taking the derivatives of the functional forms of the tomographic states given by eqs.(3.13) and (3.16), and evaluating those derivatives at the appropriate values of the arguments (see Table 1).

The errors in the angles are assumed to be uncorrelated, as would be the case if each wave-plate were adjusted for each of the 16 measurements. In reality, for qubit experiments, only one or two of the four waveplates are adjusted between one measurement and the next. However the assumption of uncorrelated angular errors greatly simplifies the calculation (which is, after all, only an estimate of the errors), and seems to produce reasonable figures for our error bars In other experimental circumstances, such as the measurement of the joint state of two spin $`1/2`$ particles, the tomography would be realized by performing unitary operations on the spins prior to measurement. In this case, an assumption analogous to ours will be wholly justified. Thus with the following assumption

$$\overline{\delta \theta _{\nu ,i}\delta \theta _{\mu ,j}}=\delta _{\nu ,\mu }\delta _{i,j}(\mathrm{\Delta }\theta )^2,$$
(5.12)
(where $`\mathrm{\Delta }\theta `$ is the RMS uncertainty in the setting of the waveplate, with an estimated value of $`0.25^o`$ for our apparatus) we obtain the following expression for the errors in $`s_\nu `$ due to angular settings:

$$\overline{\delta s_\nu \delta s_\mu }=\delta _{\nu ,\mu }\underset{i=1}{\overset{4}{}}\underset{ϵ,\lambda =1}{\overset{16}{}}f_{\nu ,ϵ}^{(i)}f_{\nu ,\lambda }^{(i)}s_ϵs_\lambda $$
(5.13)
Combining eqs.(5.13) and (5.7) we obtain the following formula for the total error in the quantities $`s_\nu `$:

$$\overline{\delta s_\nu \delta s_\mu }=\delta _{\nu ,\mu }\mathrm{\Lambda }_\nu $$
(5.14)
where

$$\mathrm{\Lambda }_\nu =\left[\frac{s_\nu }{𝒩}+\underset{i=1}{\overset{4}{}}\underset{ϵ,\lambda =1}{\overset{16}{}}f_{\nu ,ϵ}^{(i)}f_{\nu ,\lambda }^{(i)}s_ϵs_\lambda \right].$$
(5.15)
These sixteen quantities can be calculated using the parameters $`s_\nu `$ and the constants $`f_{\nu ,ϵ}^{(i)}`$. Note that the same result can be obtained by assuming a priori that the errors in the $`s_\nu `$ are all uncorrelated, with $`\mathrm{\Lambda }_\nu =\delta s_\nu ^2`$; the more rigorous treatment given here is however necessary to demonstrate this fact. For a typical number of counts, say $`𝒩=\mathcal{10000}`$ it is found that the contribution of errors from the two causes is roughly comparable; for larger numbers of counts, the angular settings will become the dominant source of error.

Based on these results, the errors in the values of the various elements of the density matrix estimated by the linear tomographic technique described in Section 3 are as follows:

$`\left(\mathrm{\Delta }\rho _{i,j}\right)^2`$ $`=`$ $`{\displaystyle \underset{\nu ,\mu =1}{\overset{16}{}}}{\displaystyle \frac{\rho _{i,j}}{s_\nu }}{\displaystyle \frac{\rho _{i,j}}{s_\mu }}\overline{\delta s_\nu \delta s_\mu }`$ (5.16)
$`=`$ $`{\displaystyle \underset{\nu =1}{\overset{16}{}}}\left(M_{\nu (i,j)}\right)^2\mathrm{\Lambda }_\nu `$ (5.17)
where $`M_{\nu (i,j)}`$ is the $`i,j`$ element of the matrix $`\widehat{M}_\nu `$.

A convenient way in which to estimate errors for a maximum likelihood tomographic technique (rather than a linear tomographic technique) is to employ the above formulae, with the slight modification that the parameter $`s_\nu `$ should be recalculated from eq.(3.18) using the estimated density matrix $`\widehat{\rho }_{est}`$. This does not take into account errors inherent in the maximum likelihood technique itself.

### C Errors in Quantities Derived from the Density Matrix

When calculating the propagation of errors, it is actually more convenient to use the errors in the $`s_\nu `$ parameters, (given by eq.(5.15), rather than the errors in the elements of density matrix itself (which have non-negligible correlations).

#### 1 von Neumann Entropy

The von Neumann entropy is an important measure of the purity of a quantum state $`\widehat{\rho }`$. It is defined by

$`𝒮`$ $`=`$ $`\mathrm{Tr}\left\{\widehat{\rho }\mathrm{log}_2\left(\widehat{\rho }\right)\right\}`$ (5.18)
$`=`$ $`{\displaystyle \underset{a=1}{\overset{4}{}}}p_a\mathrm{log}_2p_a,`$ (5.19)
where $`p_a`$ is an eigenvalue of $`\widehat{\rho }`$, i.e.

$$\widehat{\rho }|\varphi _a=p_a|\varphi _a,$$
(5.20)
$`|\varphi _a`$ being the $`a`$-th eigenstate ($`a`$= 1,…4). The error in this quantity is given by

$$\left(\mathrm{\Delta }𝒮\right)^2=\underset{\nu =1}{\overset{16}{}}\left(\frac{𝒮}{s_\nu }\right)^2\mathrm{\Lambda }_\nu .$$
(5.21)
Applying the chain rule, we find

$$\left(\frac{𝒮}{s_\nu }\right)=\underset{a=1}{\overset{4}{}}\left(\frac{p_a}{s_\nu }\right)\left(\frac{𝒮}{p_a}\right).$$
(5.22)
The partial differential of a eigenvalue can be easily found by perturbation theory. As is well known (e.g. ) the change in the eigenvalue $`\lambda _a`$ of a matrix $`\widehat{W}`$ due to a perturbation in the matrix $`\widehat{\delta W}`$ is

$$\delta \lambda _a=\varphi _a|\widehat{\delta W}|\varphi _a,$$
(5.23)
where $`|\varphi _a`$ is the eigenvector of $`\widehat{W}`$ corresponding to the eigenvalue $`\lambda _a`$. Thus the derivative of $`\lambda _a`$ with respect to some variable $`x`$ is given by

$$\frac{\lambda _a}{x}=\varphi _a|\frac{\widehat{W}}{x}|\varphi _a.$$
(5.24)
Since $`\widehat{\rho }=_{\nu =1}^{16}\widehat{M}_\nu s_\nu `$, we find that

$$\frac{p_a}{s_\nu }=\varphi _a|\widehat{M}_\nu |\varphi _a$$
(5.25)
and so, taking the derivative of eq.(5.19), eq.(5.22) becomes

$$\left(\frac{𝒮}{s_\nu }\right)=\underset{a=1}{\overset{4}{}}\varphi _a|\widehat{M}_\nu |\varphi _a\frac{\left[1+\mathrm{ln}p_a\right]}{\mathrm{ln}2}.$$
(5.26)
Hence

$$\left(\mathrm{\Delta }𝒮\right)^2=\underset{\nu =1}{\overset{16}{}}\left(\underset{a=1}{\overset{4}{}}\varphi _a|\widehat{M}_\nu |\varphi _a\frac{\left[1+\mathrm{ln}p_a\right]}{\mathrm{ln}2}\right)^2\mathrm{\Lambda }_\nu .$$
(5.27)
For the experimental example given above, $`𝒮=0.106\pm 0.049`$.

#### 2 Linear Entropy

The “linear entropy” is used to quantify the degree of mixture of a quantum state in an analytically convenient form, although unlike the von Neumann entropy it has no direct information theoretic implications. In a normalized form (defined so that its value lies between zero and one), the linear entropy for a two qubit system is defined by:

$`𝒫`$ $`=`$ $`{\displaystyle \frac{4}{3}}\left(1\mathrm{Tr}\left\{\widehat{\rho }^2\right\}\right)`$ (5.28)
$`=`$ $`{\displaystyle \frac{4}{3}}\left(1{\displaystyle \underset{a=1}{\overset{4}{}}}p_a^2\right).`$ (5.29)
To calculate the error in this quantity, we need the following partial derivative:

$`{\displaystyle \frac{𝒫}{s_\nu }}`$ $`=`$ $`{\displaystyle \frac{8}{3}}{\displaystyle \underset{a=1}{\overset{4}{}}}p_a{\displaystyle \frac{p_a}{s_\nu }}`$ (5.30)
$`=`$ $`{\displaystyle \frac{8}{3}}{\displaystyle \underset{a=1}{\overset{4}{}}}p_a\varphi _a|\widehat{M}_\nu |\varphi _a`$ (5.31)
$`=`$ $`{\displaystyle \frac{8}{3}}\mathrm{Tr}\left\{\widehat{\rho }\widehat{M}_\nu \right\}`$ (5.32)
$`=`$ $`{\displaystyle \frac{8}{3}}{\displaystyle \underset{\mu =1}{\overset{16}{}}}\mathrm{Tr}\left\{\widehat{M}_\mu \widehat{M}_\nu \right\}s_\mu `$ (5.33)
Hence the error in the linear entropy is

$`\left(\mathrm{\Delta }𝒫\right)^2`$ $`=`$ $`{\displaystyle \underset{\nu =1}{\overset{16}{}}}\left({\displaystyle \frac{𝒫}{s_\nu }}\right)^2\mathrm{\Lambda }_\nu ,`$ (5.34)
$`=`$ $`{\displaystyle \underset{\nu }{\overset{16}{}}}\left({\displaystyle \frac{8}{3}}{\displaystyle \underset{\mu =1}{\overset{16}{}}}\mathrm{Tr}\left\{\widehat{M}_\mu \widehat{M}_\nu \right\}s_\mu \right)^2\mathrm{\Lambda }_\nu .`$ (5.35)
For the example given in Sections III and IV, $`𝒫=0.037\pm 0.026`$.

#### 3 Concurrence, Entanglement of Formation and Tangle

The concurrence, entanglement of formation and tangle are measures of the quantum-coherence properties of a mixed quantum state . For two qubits The analysis in this subsection applies to the two qubit case only. Measures of entanglement for mixed n-qubit systems are a subject of on-going research: see, for example,  for a recent survey. , concurrence is defined as follows: consider the non-Hermitian matrix $`\widehat{R}=\widehat{\rho }\widehat{\mathrm{\Sigma }}\widehat{\rho }^\mathrm{T}\widehat{\mathrm{\Sigma }}`$ where the superscript $`\mathrm{T}`$ denotes transpose and the “spin flip matrix” $`\widehat{\mathrm{\Sigma }}`$ is defined by:

$$\widehat{\mathrm{\Sigma }}=\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\end{array}\right).$$
(5.36)
Note that the definition of $`\widehat{\mathrm{\Sigma }}`$ depends on the basis chosen; we have assumed here the “computational basis” $`\{|HH,|HV,|VH,|VV\}`$. In what follows, it will be convenient to write $`\widehat{R}`$ in the following form:

$$\widehat{R}=\frac{1}{2}\underset{\mu ,\nu =1}{\overset{16}{}}\widehat{q}_{\mu ,\nu }s_\mu s_\nu ,$$
(5.37)
where $`\widehat{q}_{\mu ,\nu }=\widehat{M}_\mu \widehat{\mathrm{\Sigma }}\widehat{M}_\nu ^\mathrm{T}\widehat{\mathrm{\Sigma }}+\widehat{M}_\nu \widehat{\mathrm{\Sigma }}\widehat{M}_\mu ^\mathrm{T}\widehat{\mathrm{\Sigma }}`$. The left and right eigenstates and eigenvalues of the matrix $`\widehat{R}`$ we shall denote by $`\xi _a|`$, $`|\zeta _a`$ and $`r_a`$, respectively, i.e.:

$`\xi _a|\widehat{R}`$ $`=`$ $`r_a\xi _a|`$ (5.38)
$`\widehat{R}|\zeta _a`$ $`=`$ $`r_a|\zeta _a.`$ (5.39)
We shall assume that these eigenstates are normalized in the usual manner for bi-orthogonal expansions, i.e. $`\xi _a|\zeta _b=\delta _{a,b}`$. Further we shall assume that the eigenvalues are numbered in decreasing order, so that $`r_1r_2r_3r_4`$. The concurrence is then defined by the formula

$`C`$ $`=`$ $`\mathrm{Max}\{0,\sqrt{r_1}\sqrt{r_2}\sqrt{r_3}\sqrt{r_4}\}`$ (5.40)
$`=`$ $`\mathrm{Max}\{0,{\displaystyle \underset{a=1}{\overset{4}{}}}\mathrm{sgn}\left({\displaystyle \frac{3}{2}}a\right)\sqrt{r_a}\},`$ (5.41)
where $`\mathrm{sgn}(x)=1`$ if $`x>0`$ and $`\mathrm{sgn}(x)=1`$ if $`x<0`$. The tangle is given by $`T=C^2`$ and the Entanglement of Formation by

$$E=h\left(\frac{1+\sqrt{1C^2}}{2}\right),$$
(5.42)
where $`h(x)=x\mathrm{log}_2x(1x)\mathrm{log}_2(1x)`$. Because $`h(x)`$ is a monotonically increasing function, these three quantities are to some extent equivalent measures of the entanglement of a mixed state.

To calculate the errors in these rather complicated functions, we must employ the perturbation theory for non-Hermitian matrices (see Appendix C for more details). We need to evaluate the following partial derivative,

$`{\displaystyle \frac{C}{s_\nu }}`$ $`=`$ $`{\displaystyle \underset{a=1}{\overset{4}{}}}\mathrm{sgn}\left({\displaystyle \frac{3}{2}}a\right){\displaystyle \frac{1}{2\sqrt{r_a}}}{\displaystyle \frac{r_a}{s_\nu }}`$ (5.43)
$`=`$ $`{\displaystyle \underset{a=1}{\overset{4}{}}}\mathrm{sgn}\left({\displaystyle \frac{3}{2}}a\right){\displaystyle \frac{1}{2\sqrt{r_a}}}\xi _a|{\displaystyle \frac{\widehat{R}}{s_\nu }}|\zeta _a`$ (5.44)
$`=`$ $`{\displaystyle \underset{a=1}{\overset{4}{}}}{\displaystyle \underset{\mu =1}{\overset{16}{}}}\mathrm{sgn}\left({\displaystyle \frac{3}{2}}a\right){\displaystyle \frac{1}{2\sqrt{r_a}}}\xi _a|\widehat{q}_{\mu ,\nu }s_\mu |\zeta _a,`$ (5.45)
where the function $`\mathrm{sgn}\left(x\right)`$ is the sign of the quantity $`x`$: it takes the value $`1`$ if $`x>0`$ and $`1`$ if $`x<0`$. Thus $`\mathrm{sgn}\left(3/2a\right)`$ is equal to $`+1`$ if $`a=1`$ and $`1`$ if $`a=2,3\mathrm{or}\mathrm{\hspace{0.17em}\hspace{0.17em}4}`$. Hence the error in the concurrence is

$`\left(\mathrm{\Delta }C\right)^2`$ $`=`$ $`{\displaystyle \underset{\nu =1}{\overset{16}{}}}\left({\displaystyle \frac{C}{s_\nu }}\right)^2\mathrm{\Lambda }_\nu `$ (5.46)
$`=`$ $`{\displaystyle \underset{\nu =1}{\overset{16}{}}}\left({\displaystyle \underset{a=1}{\overset{4}{}}}{\displaystyle \underset{\mu =1}{\overset{16}{}}}\mathrm{sgn}\left({\displaystyle \frac{3}{2}}a\right){\displaystyle \frac{1}{2\sqrt{r_a}}}\xi _a|\widehat{q}_{\mu ,\nu }s_\mu |\zeta _a\right)^2\mathrm{\Lambda }_\nu .`$ (5.47)
For our example the concurrence is $`0.963\pm 0.018`$.

Once we know the error in the concurrence, the errors in the tangle and the entanglement of formation can be found straightforwardly:

$`\mathrm{\Delta }T`$ $`=`$ $`2C\mathrm{\Delta }C`$ (5.48)
$`\mathrm{\Delta }E`$ $`=`$ $`{\displaystyle \frac{C}{\sqrt{1C^2}}}h^{}\left({\displaystyle \frac{1+\sqrt{1C^2}}{2}}\right)\mathrm{\Delta }C,`$ (5.49)
where $`h^{}(x)`$ is the derivative of $`h(x)`$. For our example the the tangle is $`0.928\pm 0.034`$ and the entanglement of formation is $`0.947\pm 0.025`$.

## VI Conclusions

In conclusion, we have presented a technique for reconstructing density matrices of qubit systems, including a full error analysis. We have extended the latter through to calculation of quantities of interest in quantum information, such as the entropy and concurrence. Without loss of generality, we have used the example of polarization qubits of entangled photons, but we stress that these techniques can be adapted to any physical realization of qubits.

## Acknowledgements

The authors would like to thank Joe Altepeter, Mauro d’Ariano, Zdenek Hradil, Kurt Jacobs, Poul Jessen, Michael Neilsen, Mike Raymer, Sze Tan, and Jaroslav Řeháček for useful discussions and correspondence. This work was supported in part by the U.S. National Security Agency, and Advanced Research and Development Activity (ARDA), by the Los Alamos National Laboratory LDRD program and by the Australian Research Council.

## Appendix A: The $`\widehat{\mathrm{\Gamma }}`$-matrices

One possible set of $`\widehat{\mathrm{\Gamma }}`$-matrices are generators of $`SU(2)SU(2)`$, normalized so that the conditions given in eq.3.20 are fulfilled. These matrices are:

$$\begin{array}{cccc}\widehat{\mathrm{\Gamma }}_1=\frac{1}{2}\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_2=\frac{1}{2}\left(\begin{array}{cccc}0& i& 0& 0\\ i& 0& 0& 0\\ 0& 0& 0& i\\ 0& 0& i& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_3=\frac{1}{2}\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right),& \widehat{\mathrm{\Gamma }}_4=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right),\\ & & & \\ \widehat{\mathrm{\Gamma }}_5=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_6=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 0& i\\ 0& 0& i& 0\\ 0& i& 0& 0\\ i& 0& 0& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_7=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_8=\frac{1}{2}\left(\begin{array}{cccc}0& 0& i& 0\\ 0& 0& 0& i\\ i& 0& 0& 0\\ 0& i& 0& 0\end{array}\right),\\ & & & \\ \widehat{\mathrm{\Gamma }}_9=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 0& i\\ 0& 0& i& 0\\ 0& i& 0& 0\\ i& 0& 0& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_{10}=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_{11}=\frac{1}{2}\left(\begin{array}{cccc}0& 0& i& 0\\ 0& 0& 0& i\\ i& 0& 0& 0\\ 0& i& 0& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_{12}=\frac{1}{2}\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right),\\ & & & \\ \widehat{\mathrm{\Gamma }}_{13}=\frac{1}{2}\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_{14}=\frac{1}{2}\left(\begin{array}{cccc}0& i& 0& 0\\ i& 0& 0& 0\\ 0& 0& 0& i\\ 0& 0& i& 0\end{array}\right),& \widehat{\mathrm{\Gamma }}_{15}=\frac{1}{2}\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right),& \widehat{\mathrm{\Gamma }}_{16}=\frac{1}{2}\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right).\end{array}$$
(A.1)
As noted in the text, this is only one possible choice for these matrices, and the final results are independent of the choice.

## Appendix B: The $`\widehat{M}`$-matrices and some of their properties

The $`\widehat{M}`$ matrices, defined by eq.(3.28), are as follows:

$$\begin{array}{cc}\widehat{M}_1=\frac{1}{2}\left(\begin{array}{cccc}2& (1i)& (1+i)& 1\\ (1+i)& 0& i& 0\\ (1i)& i& 0& 0\\ 1& 0& 0& 0\end{array}\right),& \widehat{M}_2=\frac{1}{2}\left(\begin{array}{cccc}0& (1i)& 0& 1\\ (1+i)& 2& i& (1+i)\\ 0& i& 0& 0\\ 1& (1+i)& 0& 0\end{array}\right),\\ & \\ \widehat{M}_3=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& i& (1+i)\\ 0& i& 0& (1i)\\ 1& (1i)& (1+i)& 2\end{array}\right),& \widehat{M}_4=\frac{1}{2}\left(\begin{array}{cccc}0& 0& (1+i)& 1\\ 0& 0& i& 0\\ (1i)& i& 2& (1i)\\ 1& 0& (1+i)& 0\end{array}\right),\\ & \\ \widehat{M}_5=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 2i& (1+i)\\ 0& 0& (1i)& 0\\ 2i& (1+i)& 0& 0\\ (1i)& 0& 0& 0\end{array}\right),& \widehat{M}_6=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 0& (1+i)\\ 0& 0& (1i)& 2i\\ 0& (1+i)& 0& 0\\ (1i)& 2i& 0& 0\end{array}\right),\\ & \\ \widehat{M}_7=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 0& (1+i)\\ 0& 0& (1i)& 2\\ 0& (1+i)& 0& 0\\ (1i)& 2& 0& 0\end{array}\right),& \widehat{M}_8=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 2& (1+i)\\ 0& 0& (1i)& 0\\ 2& (1+i)& 0& 0\\ (1i)& 0& 0& 0\end{array}\right),\\ & \\ \widehat{M}_9=\left(\begin{array}{cccc}0& 0& 0& i\\ 0& 0& i& 0\\ 0& i& 0& 0\\ i& 0& 0& 0\end{array}\right),& \widehat{M}_{10}=\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\end{array}\right),\\ & \\ \widehat{M}_{11}=\left(\begin{array}{cccc}0& 0& 0& i\\ 0& 0& i& 0\\ 0& i& 0& 0\\ i& 0& 0& 0\end{array}\right),& \widehat{M}_{12}=\frac{1}{2}\left(\begin{array}{cccc}0& 2& 0& (1+i)\\ 2& 0& (1+i)& 0\\ 0& (1i)& 0& 0\\ (1i)& 0& 0& 0\end{array}\right),\\ & \\ \widehat{M}_{13}=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 0& (1+i)\\ 0& 0& (1+i)& 0\\ 0& (1i)& 0& 2\\ (1i)& 0& 2& 0\end{array}\right),& \widehat{M}_{14}=\frac{1}{2}\left(\begin{array}{cccc}0& 0& 0& (1i)\\ 0& 0& (1i)& 0\\ 0& (1+i)& 0& 2i\\ (1+i)& 0& 2i& 0\end{array}\right),\\ & \\ \widehat{M}_{15}=\frac{1}{2}\left(\begin{array}{cccc}0& 2i& 0& (1i)\\ 2i& 0& (1i)& 0\\ 0& (1+i)& 0& 0\\ (1+i)& 0& 0& 0\end{array}\right),& \widehat{M}_{16}=\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\end{array}\right).\end{array}$$
(B.1)
The form of these matrices is independent of the chosen set of matrices $`\left\{\widehat{\mathrm{\Gamma }}_\nu \right\}`$ used to convert the density matrix into a column vector. However the $`\widehat{M}_\nu `$ matrices do depend on the set of tomographic states $`|\psi _\nu `$.

There are some useful properties of these matrices which we will now derive. From eq.(3.28), we have

$$\psi _\mu |\widehat{M}_\nu |\psi _\mu =\underset{\lambda }{}\psi _\mu |\widehat{\mathrm{\Gamma }}_\lambda |\psi _\mu \left(B^1\right)_{\lambda ,\nu }.$$
(B.2)
From eq.(3.24) we have $`\psi _\mu |\widehat{\mathrm{\Gamma }}_\lambda |\psi _\mu =B_{\mu ,\lambda }`$, thus we obtain the result

$$\psi _\mu |\widehat{M}_\nu |\psi _\mu =\delta _{\mu ,\nu }.$$
(B.3)
If we denote the basis set for the four-dimensional Hilbert space by $`\left\{|i(i=1,2,3,4)\right\}`$, then eq.(3.27) can be written as follows:

$$i|\widehat{\rho }|j=\underset{k,l}{}\underset{\nu }{}i|\widehat{M}_\nu |j\psi _\nu |kl|\psi _\nu k|\widehat{\rho }|l.$$
(B.4)
Since eq.(B.4) is valid for arbitrary states $`\widehat{\rho }`$, we obtain the following relationship:

$$\underset{\nu }{}i|\widehat{M}_\nu |j\psi _\nu |kl|\psi _\nu =\delta _{ik}\delta _{jl}.$$
(B.5)
Contracting eq.(B.5) over the indices $`(i,j)`$ we obtain:

$$\underset{\nu }{}\mathrm{Tr}\left\{\widehat{M}_\nu \right\}|\psi _\nu \psi _\nu |=\widehat{I},$$
(B.6)
where $`\widehat{I}`$ is the identity operator for our four dimensional Hilbert space.

A second relationship can be obtained by contracting eq.(B.5), viz:

$$\underset{\nu }{}i|\widehat{M}_\nu |j=\delta _{ij},$$
(B.7)
or, in operator notation,

$$\underset{\nu }{}\widehat{M}_\nu =\widehat{I}.$$
(B.8)
## Appendix C: Perturbation theory for non-Hermitian matrices

Whereas perturbation theory for Hermitian matrices is covered in most quantum mechanics text-books, the case of non-Hermitian matrices is less familiar, and so we will present it here. The problem is, given the eigenspectrum of a matrix $`\widehat{R}_0`$ , i.e.:

$`\xi _a|\widehat{R}_0`$ $`=`$ $`r_a\xi _a|`$ (C.1)
$`\widehat{R}_0|\zeta _a`$ $`=`$ $`r_a|\zeta _a,`$ (C.2)
where

$$\xi _a|\zeta _b=\delta _{a,b}$$
(C.3)
we wish to find expressions for the eigenvalues $`r_a^{}`$ and eigenstates $`\xi _a^{}|`$ and $`|\zeta _a^{}`$ of the perturbed matrix $`\widehat{R}^{}=\widehat{R}_0+\delta \widehat{R}`$.

We start with the standard assumption of perturbation theory, i.e. that the perturbed quantities $`r_a^{}`$, $`\xi _a^{}|`$ and $`|\zeta _a^{}`$ can be expressed as power series of some parameter $`\lambda `$:

$`r_a^{}`$ $`=`$ $`r_a^{(0)}+\lambda r_a^{(1)}+\lambda ^2r_a^{(2)}+\mathrm{}`$ (C.4)
$`|\zeta _a^{}`$ $`=`$ $`|\zeta _a^{(0)}+\lambda |\zeta _a^{(1)}+\lambda ^2|\zeta _a^{(2)}+\mathrm{}`$ (C.5)
$`\xi _a^{}|`$ $`=`$ $`\xi _a^{(0)}|+\lambda \xi _a^{(1)}|+\lambda ^2\xi _a^{(2)}|+\mathrm{}`$ (C.6)
Writing $`\widehat{R}^{}=\widehat{R}_0+\lambda \delta \widehat{R}`$, and comparing terms of equal powers of $`\lambda `$ in the eigen equations, one obtains the following formulae:

$`\widehat{R}_0|\zeta _a^{(0)}`$ $`=`$ $`r_a^{(0)}|\zeta _a^{(0)}`$ (C.7)
$`\xi _a^{(0)}|\widehat{R}_0`$ $`=`$ $`r_a^{(0)}\xi _a^{(0)}|`$ (C.8)
$`\left(\widehat{R}_0r_a^{(0)}\widehat{I}\right)|\zeta _a^{(1)}`$ $`=`$ $`\left(\delta \widehat{R}r_a^{(1)}\right)|\zeta _a^{(0)}`$ (C.9)
$`\xi _a^{(1)}|\left(\widehat{R}_0r_a^{(0)}\widehat{I}\right)`$ $`=`$ $`\xi _a^{(0)}|\left(\delta \widehat{R}r_a^{(1)}\right).`$ (C.10)
Equations (C.7) and (C.8) imply that, as might be expected,

$`|\zeta _a^{(0)}`$ $`=`$ $`|\zeta _a`$ (C.11)
$`\xi _a^{(0)}|`$ $`=`$ $`\xi _a|`$ (C.12)
$`r_a^{(0)}`$ $`=`$ $`r_a`$ (C.13)
Taking the inner product of eq.(C.9) with $`\xi _a|`$, and using the bi-orthogonal property eq.(C.3), we obtain

$$r_a^{(1)}=\xi _a|\delta \widehat{R}|\zeta _a.$$
(C.14)
This implies that

$`\delta r_a`$ $``$ $`r_a^{}r_a`$ (C.15)
$``$ $`\xi _a|\delta \widehat{R}|\zeta _a.`$ (C.16)
Thus, dividing both sides by some differential increment $`\delta x`$ and taking the limit $`\delta x0`$, we obtain

$$\frac{r_a}{x}=\xi _a|\frac{\widehat{R}}{x}|\zeta _a.$$
(C.17)
Using the completeness property of the eigenstates, $`_b|\zeta _b\xi _b|=\widehat{\mathrm{I}}`$ and the identity $`\widehat{R}_0=_br_b|\zeta _b\xi _b|`$ we obtain the following formula

$$\left(\widehat{R}_0r_a\widehat{\mathrm{I}}\right)^1=\underset{\stackrel{b}{ba}}{}\frac{1}{r_br_a}|\zeta _b\xi _b|.$$
(C.18)
Applying this to eq.(C.9) we obtain

$`|\delta \zeta _a^{(1)}`$ $``$ $`|\zeta _a^{}|\zeta _a`$ (C.19)
$``$ $`{\displaystyle \underset{\stackrel{b}{ba}}{}}\left({\displaystyle \frac{\xi _b|\delta \widehat{R}|\zeta _a}{r_br_a}}\right)|\zeta _b.`$ (C.20)
Similarly, eqs.(C.10) and (C.18) imply

$`\delta \xi _a|`$ $``$ $`\delta \xi _a^{}|\delta \xi _a|`$ (C.21)
$``$ $`{\displaystyle \underset{\stackrel{b}{ba}}{}}\left({\displaystyle \frac{\xi _a|\delta \widehat{R}|\zeta _b}{r_br_a}}\right)\xi _b|.`$ (C.22)
| $`\nu `$ | Mode 1 | Mode 2 | $`h_1`$ | $`q_1`$ | $`h_2`$ | $`q_2`$ |
| --- | --- | --- | --- | --- | --- | --- |
| 1 | $`|\mathrm{H}`$ | $`|\mathrm{H}`$ | $`45^o`$ | $`0`$ | $`45^o`$ | $`0`$ |
| 2 | $`|\mathrm{H}`$ | $`|\mathrm{V}`$ | $`45^o`$ | $`0`$ | $`0`$ | $`0`$ |
| 3 | $`|\mathrm{V}`$ | $`|\mathrm{V}`$ | $`0`$ | $`0`$ | $`0`$ | $`0`$ |
| 4 | $`|\mathrm{V}`$ | $`|\mathrm{H}`$ | $`0`$ | $`0`$ | $`45^o`$ | $`0`$ |
| 5 | $`|\mathrm{R}`$ | $`|\mathrm{H}`$ | $`22.5^o`$ | $`0`$ | $`45^o`$ | $`0`$ |
| 6 | $`|\mathrm{R}`$ | $`|\mathrm{V}`$ | $`22.5^o`$ | $`0`$ | $`0`$ | $`0`$ |
| 7 | $`|\mathrm{D}`$ | $`|\mathrm{V}`$ | $`22.5^o`$ | $`45^o`$ | $`0`$ | $`0`$ |
| 8 | $`|\mathrm{D}`$ | $`|\mathrm{H}`$ | $`22.5^o`$ | $`45^o`$ | $`45^o`$ | $`0`$ |
| 9 | $`|\mathrm{D}`$ | $`|\mathrm{R}`$ | $`22.5^o`$ | $`45^o`$ | $`22.5^o`$ | $`0`$ |
| 10 | $`|\mathrm{D}`$ | $`|\mathrm{D}`$ | $`22.5^o`$ | $`45^o`$ | $`22.5^o`$ | $`45^o`$ |
| 11 | $`|\mathrm{R}`$ | $`|\mathrm{D}`$ | $`22.5^o`$ | $`0`$ | $`22.5^o`$ | $`45^o`$ |
| 12 | $`|\mathrm{H}`$ | $`|\mathrm{D}`$ | $`45^o`$ | $`0`$ | $`22.5^o`$ | $`45^o`$ |
| 13 | $`|\mathrm{V}`$ | $`|\mathrm{D}`$ | $`0`$ | $`0`$ | $`22.5^o`$ | $`45^o`$ |
| 14 | $`|\mathrm{V}`$ | $`|\mathrm{L}`$ | $`0`$ | $`0`$ | $`22.5^o`$ | $`90^o`$ |
| 15 | $`|\mathrm{H}`$ | $`|\mathrm{L}`$ | $`45^o`$ | $`0`$ | $`22.5^o`$ | $`90^o`$ |
| 16 | $`|\mathrm{R}`$ | $`|\mathrm{L}`$ | $`22.5^o`$ | $`0`$ | $`22.5^o`$ | $`90^o`$ |
| Table 1 | | |  |  |  |  |

TABLE 1: The tomographic analysis states used in our experiments. The number of coincidence counts measured in projections measurements provide a set of 16 data that allow the density matrix of the state of the two modes to be estimated. We have used the notation $`|\mathrm{D}\left(|\mathrm{H}+|\mathrm{V}\right)/\sqrt{2}`$, $`|\mathrm{L}\left(|\mathrm{H}+i|\mathrm{V}\right)/\sqrt{2}`$ and $`|\mathrm{R}\left(|\mathrm{H}i|\mathrm{V}\right)/\sqrt{2}`$. Note that, when the measurement are taken in the order given by the table, only one waveplate angle had to be changed between each measurement.


text:
# Quantum Theory Cannot Forbid Superluminal Signaling

(April 7,2001)
This note analyzes the angular distributions of the probabilities of two-photon states come out of the single-photon’s stimulated emission amplification by means of a single-atom amplifier, to see that whether the quantum theory can forbid us exploiting EPR photon pairs combined with stimulated emission to realize superluminal signaling. Besides, we stressed that superluminal signaling will leads to a dilemma of causality in a system with two long superluminal channels and two short light channels.

It has been believed that the mathematical inseparability of the quantum theoretical representation is an essential part of nature, not a mere accident of the formalism. However, the attempts to realize superluminal signaling over the last twenty years by means of EPR pairs has not been successful because one cannot clone a single particle in an unknown state sufficiently well . Now let us see whether quantum theory can forbid superluminal signaling through a careful analysis of a concrete physical process.

Consider the stimulated emission of a single excited atom, or a single-atom light amplifier. Suppose the angle between the polarization direction of an incoming single-photon flow and the atom’s transition dipole moment $`\stackrel{}{\mu }`$ is $`\theta `$ , where $`\stackrel{}{\mu }`$ is perpendicular to the photon’s wave vector. We note the photon’s state by $`|\theta `$ . After scattering, the system’s two-photon state is of the form

$$|\mathrm{\Psi }_f^\theta =\alpha _\theta |2,0^\theta |g_\theta +\beta _\theta |1,1^\theta |g_{\theta +\frac{\pi }{2}},(1)$$

where $`|2,0^\theta `$ indicates that both photons are in the state $`|\theta `$ , $`|1,1^\theta `$indicates a photon in each of the states $`|\theta `$ and $`|\theta +\frac{\pi }{2}`$ , and $`|g_\theta `$ and $`|g_{\theta +\frac{\pi }{2}}`$ are the atom’s final states. The state $`|0,2^\theta `$ would indicate that both photons are in the state $`|\theta +\frac{\pi }{2}`$ . Let $`|\mathrm{\Omega }^\theta `$ be a vector with three components $`|2,0^\theta `$ ,$`|0,2^\theta `$ and $`|1,1^\theta `$ . We have $`|\mathrm{\Omega }^\theta =U(\theta )|\mathrm{\Omega }^0`$ , where $`U(\theta )`$ is a unitary transformation,

$$|2,0^\theta =\mathrm{cos}^2\theta |2,0+\mathrm{sin}^2\theta |0,2+\frac{1}{\sqrt{2}}\mathrm{sin}2\theta |1,1,(2a)$$

$$|0,2^\theta =\mathrm{sin}^2\theta |2,0+\mathrm{cos}^2\theta |0,2\frac{1}{\sqrt{2}}\mathrm{sin}2\theta |1,1,(2b)$$

$$|1,1^\theta =\frac{1}{\sqrt{2}}\mathrm{sin}2\theta |2,0+\frac{1}{\sqrt{2}}\mathrm{sin}2\theta |0,2+\mathrm{cos}2\theta |1,1.(2c)$$

Suppose $`|g_\theta |g_{\theta +\frac{\pi }{2}}`$ and the incoming photons flow are a mixture of photons in $`|\theta `$ state and in $`|\theta +\frac{\pi }{2}`$ state with the same probability. The four transition probabilities is

$$w_{2,0}^\theta =2\lambda ^2\mathrm{cos}^2\theta d\mathrm{\Omega },(3a)$$

$$w_{1,1}^\theta =\lambda ^2\mathrm{sin}^2\theta d\mathrm{\Omega },(3b)$$

$$w_{0,2}^{\theta +\frac{\pi }{2}}=2\lambda ^2\mathrm{sin}^2\theta d\mathrm{\Omega },(3c)$$

$$w_{1,1}^{\theta +\frac{\pi }{2}}=\lambda ^2\mathrm{cos}^2\theta d\mathrm{\Omega }.(3d)$$

where $`w_{2,0}^\theta `$ is the probability of generating a two-photon state $`|2,0^\theta `$ when incoming single-photon is in the state $`|\theta `$, $`w_{1,1}^\theta `$ is the probability of generating a two-photon state $`|1,1^\theta `$ when incoming single-photon is in the state $`|\theta `$, and so on; and

$$\lambda ^2=\frac{\omega ^3\mu ^2}{8\pi ^2\mathrm{}c^3}.(4)$$

Then by means of the unitary transformation $`U(\theta )`$ we can find the probability of the state $`|2,0`$ after scattering is

$$d\overline{\sigma _{2,0}^\theta }=\frac{1}{2}\lambda ^2(1+\mathrm{cos}^22\theta )d\mathrm{\Omega }.(5)$$

If we restrict our statistics to two-photon states, we have the probability of $`|2,0`$ state

$$\overline{P_{2,0}^\theta }=\frac{1}{3}(1+\mathrm{cos}^22\theta ),(6a)$$

the other two probabilities are

$$\overline{P_{1,1}^\theta }=\frac{1}{3},(6b)$$

$$\overline{P_{0,2}^\theta }=\frac{1}{3}\mathrm{sin}^22\theta .(6c)$$

\[If let $`|g_\theta =|g_{\theta +\frac{\pi }{2}}`$, we have

$$\overline{P_{2,0}^\theta }=\frac{2}{3}\mathrm{cos}^22\theta ,(7a)$$

$$\overline{P_{1,1}^\theta }=\frac{1}{3},(7b)$$

$$\overline{P_{0,2}^\theta }=\frac{2}{3}\mathrm{sin}^22\theta .(7c)]$$

We can find from the analysis that one origin of the phenomenon that $`\overline{P^\theta }`$ is dependent on parameter $`\theta `$ is the zero point energy of the light field, and the other one is likely to be the Bose-Einstein statistics of photons.

The fact that $`\overline{P^\theta }`$ is dependent on the parameter $`\theta `$ means that with the single-atom light amplifier one can distinguish the parameter $`\theta `$ of an incoming single-photon flow by measuring the probabilities of $`|2,0`$ state or $`|0,2`$ state. If Alice and Bob share a sufficiently large number of two-photon EPR pairs in the Bell states with rotation invariance, they do not need another channel to complete their communication, then there is no law of physics which will obviously stop the superluminal signaling between Alice and Bob.

It can be shown that the superluminal signaling between long distance leads to a dilemma of causality in a system with two pairs of Alice and Bob belonged to two different inertial frames separately, and with four channels: two long superluminal channels \[Alice(1),Bob(1)\] and \[Alice(2),Bob(2)\],and two short light channels \[Bob(1),Alice(2)\]and \[Bob(2),Alice\]. Then we know that if any phenomenon of superluminal signaling be found in experiment or theoretical analysis, it means directly something wrong with our understanding of nature.

Because of the success of Electrodynamics, General Relativity, and Quantum Field Theories we believe that the Lorentz covariance of our theories is right. Based on Lorentz covariance we can prove that any superluminal signaling will lead to the dilemma of causality. Then we believe that there is no superluminal signaling in nature. However, this conclusion does not imply that any non-relativistic theory cannot give out a result related with superluminal signaling at all. It is difficult to understand some people’s belief that the non-relativistic quantum theory itself can forbid superluminal signaling automatically. One may argue that he believes there is no conflict between the theory of relativity and the theory of quantum mechanics, but many people(for example, J. A. Wheeler) believe that the conflict is unavoidable. Considering the great difficulties we meet in the quantum field theory, especially in the quantum theory of gravity, we should realize that the conflict is evident and essential.