Theory and experiment for resource-efficient joint weak-measurement

Incompatible observables underlie pillars of quantum physics such as contextuality and entanglement. The Heisenberg uncertainty principle is a fundamental limitation on the measurement of the product of incompatible observables, a `joint' measurement. However, recently a method using weak measurement has experimentally demonstrated joint measurement. This method [Lundeen, J. S., and Bamber, C. Phys. Rev. Lett. 108, 070402, 2012] delivers the standard expectation value of the product of observables, even if they are incompatible. A drawback of this method is that it requires coupling each observable to a distinct degree of freedom (DOF), i.e., a disjoint Hilbert space. Typically, this `read-out' system is an unused internal DOF of the measured particle. Unfortunately, one quickly runs out of internal DOFs, which limits the number of observables and types of measurements one can make. To address this limitation, we propose and experimentally demonstrate a technique to perform a joint weak-measurement of two incompatible observables using only one DOF as a read-out system. We apply our scheme to directly measure the density matrix of photon polarization states.


Introduction
Modern quantum measurement techniques have pushed forward our understanding and ability to manipulate quantum particles. Often, fundamental and practical measurements involve the product of two or more observables of a quantum system. In particular, correlations of incompatible or non-commuting observables A and B, defined by [A, B] ≡ AB − BA = 0, are central to our understanding of entanglement [1,2] and the Heisenberg uncertainty principle. A 'joint' measurement of A and B refers to the measurement process which outputs the expectation value of the product of the two observables BA = tr (BAρ), where ρ is the density matrix of the system. The standard procedure to perform a joint measurement would be to measure observable A, then measure B. This fails for incompatible observables since the first measurement collapses the state of a particle into an eigenstate of A, erasing the information about B and randomizing its value.
In contrast, weak measurement decreases the disturbance caused by the measurement process and thereby mostly preserves the quantum state of the system, thus allowing one to obtain correlations between any chosen set of general observables, including incompatible ones [3][4][5][6][7][8][9][10]. To perform such a measurement, the observable is weakly coupled to a separate read-out system (the 'pointer') that indicates the average result of the measurement. Even though this approach refers to an individual system, weak measurement requires repeating the measurement on identically prepared systems, and averaging. This compensates for the little information that is extracted in a single trial. Weak measurement is a type of non-destructive quantum measurement that minimizes disturbance of the measured system [11]. As shown in Ref. [12], as one decreases the disturbance caused by the measurement process, one also decreases the 'predictability' of the measurement. Weak measurement has a broad range of applications from amplifying tiny signals [13][14][15] to fundamental studies on the meaning of a quantum state [16,17]. Particularly relevant to this paper are Refs. [4,5,18,19], which showed that if two observables are weakly measured, the average measurement outcome is simply the expectation value of the product of those two observables, BA . Remarkably, this holds even if A and B are incompatible, which would make BA non-Hermitian and, nominally, unobservable.
More recently the weak measurement formalism was expanded to deal with composite systems and performing a measurement of the product of two or more observables, known as a joint weak-measurement.
Joint weak-measurement has proven to be useful, for example, in experimental realizations of the Cheshire cat [20,21], the Hardy's paradox [22,23], the study of quantum dynamics, and to give insight into the role of time ordering in the quantum domain [24,25]. The ability to jointly measure incompatible observables has also shown to have many applications in the field of quantum metrology [26][27][28][29][30][31]. Joint weak-measurement of multiple observables enables sequentially probing a quantum system for characterizing its quantum evolution [9,32]. Another example is the test of the Leggett-Garg inequalities for sequential measurements of multiple observables in a single system [33].
Known methods for the realization of a joint weak-measurement are resource-intensive. Specifically, they require either interactions that involve three or more particles or a separate readout system for each observable. With a few exceptions [22], due to the absence of two-particle interactions, even single-observable weak measurement resorts to a strategy of using internal degrees of freedom (DOF) as the read-out systems [4,5,18,19]. For example, one can measure the polarization of a photon by using its posi-tion DOF as a read-out [34]. For a joint measurement, this strategy is particularly limiting given that quantum particles have a limited number of DOF. For instance, for a photon there are just four DOF: polarization, and a threedimensional wavevector (which, in turn, incorporates frequency-time and transverse positionmomentum). Due to this limitation, joint weakmeasurement experiments have never progressed beyond the product of two observables [7,12]. To overcome this constraint, the present work theoretically introduces and experimentally demonstrates a technique to perform a joint weakmeasurement of multiple observables using a single DOF as the read-out system.
We implement our technique to directly measure quantum states. This is a type of quantum state estimation where the state is fully determined by the shift of the pointer. Quantum state estimation has become an invaluable tool in the subject of quantum information, which requires verification of the quality (i.e., 'fidelity') of resource quantum states. The experimental demonstration of the direct measurement of the wave function opened up new research lines in quantum state estimation. The directness of the method means that one can obtain the complex amplitudes of a quantum state, in any chosen basis [16]. An important aspect of direct state estimation is that no optimization or complicated inversion is involved. Solving such problems is one of the key goals of current research in quantum state estimation [35][36][37]. Further work demonstrated how to estimate a general quantum state by directly measuring the density matrix [6,19], or by directly measuring phase-space quasiprobability distributions of states, such as the Dirac distribution [38,39]. Similarly, we apply our singlepointer joint weak-measurement method to directly determine any chosen element of the density matrix. Specifically, we obtain the density matrix of photon polarization states using a single pointer for the two requisite observables.
The rest of the paper is organized as follows. We start by describing weak measurement in terms of raising and lowering operators. Then, we outline the theory of our technique to perform a joint weak-measurement and introduce an important ingredient, the fractional Fourier transform. Next, we present the experimental demonstration of our technique and an application to quantum state estimation. Finally, we summarize our work and point out some future possible directions. An overview of the Fractional Fourier Transform (FrFT) is given in Appendix A.

Weak measurement
In this section, we introduce a theoretical model for quantum measurement, von Neumann's model, that is typically used to describe weak measurement. The model involves a measured quantum system S and a pointer system P [40]. The latter indicates the measured value, the read-out, on a meter. A key aspect of the model is that the pointer is also quantum mechanical. Before the measurement, S and P are in an initial product state, |I S ⊗ |φ P , here ⊗ indicates a tensor product between different Hilbert spaces and the subscript is a label of the system. Both of these symbols will be omitted in the rest of the paper. As usual, we assume that the pointer's initial spatial wave function φ(x) is a Gaussian centered at zero [40]: where σ x is the standard deviation of the position probability-distribution. The pointer's initial state happens to be same as the ground state of a harmonic oscillator. Thus, following Ref. [5], we define a lowering operator a as the operator that annihilates this pointer state, a |φ = 0. By this logic, from here on we label the pointer's initial state as |0 = |φ . As a standard lowering operator, a can be written in terms of the position x and momentum p of the pointer as follows a = x/(2σ x ) + ipσ x /h. (Note, we use the natural length-scale σ x of the system in place of the mass m and angular frequency ω that usually appear in the harmonic oscillator: σ x = h/(2mω).) Associated with a, there is a raising operator a † that fulfills [a, a † ] = 1. Similarly we can define number states |n = (a † ) n √ n! |0 . Formulating the model in terms of lowering, raising operators and number states has proven fruitful in the past [5,18,19] and will be important for what follows.
Suppose we want to measure observable A of S. Then, in the von Neumann model, one couples S to P by the following Hamiltonian, here g is a real parameter that indicates the interaction strength and we have used the usual decomposition of p in terms of a and a † . We stress that there is no physical harmonic potential in the system and thus no quantum harmonic oscillator. We are following Ref. [5] and simply using the formalism of raising and lowering operators to analyze the effect of the interaction on the pointer state. We now consider the state of the total system after the unitary evolution induced by H: 2σx is a unitless parameter that quantifies the measurement strength. In general, the evolved system is in an entangled state between S and P. For a strong interaction (γ 1), in each trial, a measurement of the position of the pointer will unambiguously indicate the value of A (though it is not particularly obvious in this harmonic oscillator formulation). So far the model is general and independent of the measurement strength. Now we consider the weak measurement regime. A weak measurement is characterized by γ 1, which allows one to approximate the evolved state as U A |I |0 = |I |0 + γA |I |1 to first order in γ. In the weak regime, the entanglement between the pointer and measured system is reduced, and the initial state |I of the particle is largely preserved. Following the work in [11], a post-selection on a final sytem state |F is performed. Mathematically, this amounts to projecting onto F | and renormalizing, after which the pointer's final state is |φ = |0 + γ F |A|I F |I |1 . Thus the pointer's final state is largely left unchanged. That is, it is mostly left in |0 , but a small component proportional to γ, is transferred to |1 due to the interaction with S.
Our goal is to identify in what manner the pointer is shifted by the interaction. To this end, we find the expectations of the position and momentum of the final pointer. These respectively appear as the real and imaginary parts of a ≡ φ |a|φ = 1 2σx x + i σx h p . Thus, using |φ from just above, one finds Consequently, the pointer is shifted from having x = p = 0 to indicating an average outcome A w = 1 2σxγ x + i σx hγ p . This average pointer shift was introduced by Aharonov, Albert, and Vaidman in Ref. [11] and is called the 'weak value'. Unlike in the standard expectation value, |F = |I and, thus, the weak value is a potentially complex quantity. In summary, the real and imaginary parts of the weak value are the average shifts of the position and momentum of the pointer, which, in turn, are given by the expectation value of the lowering operator.

Joint weak-measurement
For composite systems, one is interested in the average value of the product of observables such as BA . Universally, this involves correlations between two measurement outcomes (e.g., as in Bell's inequalities). In the von Neumann model, this corresponds to correlations between pointer distributions. This is true for both strong and weak measurements. In the latter case, the average outcome should be the joint weak-value, A number of techniques have been proposed and demonstrated to observe the joint weak-value in pointer correlations. We now briefly review these techniques. First, we review the case of compatible operators A and B. These could be two different observables of a single particle or observables acting on two different particles. Ref. [4] proposed using a separate von Neumann interaction (i.e., Eq. 2) and pointer for each observable (pointers 1 and 2). This was simplified in [5], which found that a 1 a 2 = BA w /γ 2 . This strategy of performing two separate weak measurements was experimentally demonstrated in [22].
A more challenging case, and the subject of this work, is the one in which A and B act on the same particle, but are incompatible e.g., two complementary observables such as position and momentum. Furthermore, the product BA is not Hermitian, thus it is not considered a valid observable in standard quantum mechanics. For example, naively replacing A with BA in the von Neumann Hamiltonian, Eq. 2, results in nonunitary time evolution. However, in the weak regime, a measurement of A largely preserves the quantum state of the particle allowing a subsequent measurement of B. The correlations between the outcomes of the two measurements give BA w . A technique along these lines was proposed in [18]. As with the compatible observable case above, it used a separate von Neumann interaction and pointer for each observable (pointer 1 for B and pointer 2 for A). In [19], the required correlation between the pointers was shown to be a 1 a 2 = BA w /γ 2 and experimentally demonstrated in [6,7]. In summary, for both compatible and incompatible observables, the same technique works. The drawback of the technique is that it requires one pointer for each observable.
In particular, this requirement of one pointer per observable is resource-intensive. In most implementations of weak measurement, pointers are internal DOF of the measured particle. For example, in [34] a photon's polarization is measured by coupling it to the same photon's transverse spatial DOF. In absence of inter-particle interactions, this facilitates the use of weak measurement, but quickly uses up all available internal DOF. In turn, this limits the number of observables in the product and the number of DOF that can be used in the measured system for other quantum information tasks. It is natural to ask: can we perform a joint weak-measurement with a single pointer?
The main contribution of the present work is to introduce and experimentally demonstrate such a technique. Our technique uses a sequence of two standard von Neumann interactions, each given by Eq. 2. Unlike the previous techniques, the two interactions couple the system to the same pointer. As in section 2, the total initial state is |I |0 . The first interaction U A couples the pointer to A, while the second U B couples the same pointer to B. The action of two von Neumann unitaries with equal interaction strength γ . This is the final state of the total system after the two interactions. which used correlations between two different lowering operators a 1 a 2 , we will aim to find the expectation of the product of two identical lowering operators, a 2 . Thus, we must expand the pointer state after the interaction to second order in the interaction strength γ. There are three second-order terms: m = n = 1; m = 0, n = 2; and m = 2, n = 0. Along with the zero and first order terms, this gives Now we post-select the system on a final state |F . To second order in γ, the renormalized pointer's final state is As per our aim, we now calculate the expectation value a 2 for |φ : This equation contains the weak value of the product observable BA w but also other nontrivial weak values, A 2 w and B 2 w . However, if we limit the two observables to be projectors, then A 2 = A and B 2 = B. This turns the nontrivial weak values into single-observable weak values, which we can replace with A + B w = a /γ. Using this and rearranging Eq. 7 to solve for BA w we arrive at In this way, we have expressed the joint weakvalue solely in terms of expectation values on the pointer's final state. However, an additional step is still necessary. While the expectation value of a single lowering operator is easily measured in an experiment by measuring x and p in separate trials, powers of lowering operators cannot be measured as easily. To solve this, we express a 2 using a = x 2σx + i p 2σp , where we have used σ x σ p =h 2 (which is valid since the pointer is in the minimum uncertainty state |0 ). Doing so, leads to the appearance of cross terms such as xp + px, which do not correspond to a straightforwardly physical read-out system observable.
To overcome this problem, we can use the Hermitian observable d which is an equally weighted combination of x and p: d = σ d √ 2 x σx + p σp . Here, σ x , σ p and σ d are the standard deviations of the pointer in x, p and d spaces, respectively. The d observable naturally appears in a variety of quantum systems. In the Heisenberg picture in quantum optics, the x field quadrature rotates to d after an eighth of a period of oscillation; this is equivalent to an x-p phase-space rotation of Rπ/2, with R = 1/2 where R is the rotation order. Similarly, x rotates to p after a quarter period (R = 1). Just as the Fourier Transform links x and p, the fractional Fourier Transform (FrFT) was introduced to calculate the effect of a rotation order R on a state in the Schrödinger picture [41]. In summary, there are established practical methods to physically implement FrFTs and measure d.
The reason we have introduced this new observable is that the square of d will contain the desired cross terms. Calculating d 2 and solving for the cross terms we find Upon substituting Eq. 9 in Eq. 8, we obtain an expression for the real and imaginary parts of BA w : and (11) Note that every term in Eqs. 10 -11 is a ratio of two variables with the same units, therefore each term is unitless. For the same reason, experimental scaling factors e.g., a magnification in the x domain, cancel out. Hence, characterization of experimental scaling factors is not required for the use of our technique.  In summary, Eqs. 10 -11 express the full complex joint weak-value for product observable BA in terms of Hermitian observables on the pointer's final state. As expected, the joint weak-value appears in second order powers of x and p and our new observable d. This comprises our proposed technique to weakly measure the product of incompatible observables using only a single pointer.

Experiment
In this section, we present the experimental demonstration of our proposed technique using photons. Specifically we perform a joint weak-measurement of incompatible polarization projectors. The experimental setup is shown in Fig. 1. The measured observable will be in the photon's polarization DOF. The pointer is the photon's transverse x position with probabilitydistribution given by the absolute square of the wave function in Eq. 1 with σ x = 403 µm. The photon source is a He:Ne laser at 633 nm with a power of 1.19 mW. The setup can be divided into state preparation, weak measurements, strong measurement stages, and a read-out apparatus section. In order to test our technique, we prepare a range of polarization states |I = α |H + β |V , where |H (|V ) is the horizontal (vertical) polarization. For state preparation, we use a polariz-ing beam splitter (PBS) followed by a half-wave plate (HWP), set at an angle of θ/2 with respect to the |H polarization, and a quarter-wave plate (QWP) (see the caption in Fig. 1 for setting details).
A von Neumann measurement of polarization can be performed with a birefringent crystal (e.g., a BBO crystal) acting as a beam displacer. This optical component transversely shifts the photon by ∆x = gt = 150 µm if the photon is in the |H polarization state and leaves it unshifted if it is in |V . In this way, the crystal couples the polarization observable A = |H H| to the photon's transverse spatial position x that plays the role of the pointer. The strong measurement regime is characterized by ∆x greater than σ x , in which the eigenstates of A are fully separated. Our experiment is performed in the weak measurement regime where ∆x is less than σ x .
In order to measure the product of two observables with our technique, the setup performs two weak measurements in a row. Each von Neumann interaction (i.e., Eq. 2) is achieved with a separate BBO crystal. Both crystals are aligned such that they shift the transverse profile of horizontally polarized photons in the horizontal direction x, leaving the transverse profile in the y direction unchanged. Thus, they couple to the same pointer, the x DOF. The first BBO implements a measurement of A = π H = |H H|. Before the second BBO, there is a HWP oriented at 22.5 • . This effectively rotates the second measured observable to B = π 45 • = |45 • 45 • |, with |45 • = |H +|V √ 2 . These two measurements and their read-out constitute an experimental application of our joint weak-measurement technique that uses a single pointer. Lastly, a strong measurement of polarization observable π j (j = H or V ) is performed.
In our experiment, we need the ability to measure three incompatible observables of the pointer, x, p, and d. This is the read-out of the result of the weak measurement. As we will explain, lens transformations will allow us to switch between these spatial observables, transforming them to a final transverse position x on a camera. We measure the probability distribution of the observables in Eqs. 10 -11 on a monochrome 8 bit CMOS camera with a pixel width in x of 2.2 µm. To make room for the optical lengths required for the lens transformations we add a 4f lens-pair to the imaging system (f 1 = 100 cm and f 2 = 120 cm). The 4f is positioned such that f 1 is 100 cm after the crystals. This ensures that the spatial wave function at the exit surface of the second crystal is recreated 120 cm after the f 2 lens. Our goal is to leave the camera fixed in place while different lenses are inserted in order to measure x, p, and d.
To measure p, and d we use an optical FrFT of the spatial DOF. The special case of rotation order R = 1 (a standard Fourier Transform), is already widely used; the transverse position x at one focal length after lens f p = 100 cm is proportional to p at any distance before the lens.
Hence, lens f p can be placed at any distance after the 4f lens pair as long as it is f p distance from the camera. Less common is the optical spatial FrFT, which was introduced in [42-44] (more details in Appendix A). At a distance z after lens f d = 100 cm, x will be proportional to the d observable a distance z before the lens. Here, For d, the phase-space rotation parameter R equals 1/2 making z = 29 cm. This d lens transformation fixes the distance of the camera from the 4f lens pair. Lastly, a single lens (f x = 12.5 cm) placed 160 cm after f 2 relays the image from the 4f lens pair. This ensures that 18 cm after f x , x on the camera is proportional to x at the crystals. The values of f x , f p , and f d were chosen so that each measured x distribution spans many pixels. By switching in one lens at a time, f x , f p , or f d , the camera effectively measures the corresponding observable.
A key experimental simplification is that we do not need to experimentally or theoretically determine the proportionality constants between x at the camera and x, p, and d. The imaging magnification between x and x is an example of such a proportionality constant. Since they depend on the focal lengths and lens-camera distances, these constants are difficult to experimentally determine precisely. Instead, Eqs. 10 -11 show that each observable is divided by the width of the pointer's initial distribution in that observable, e.g. p/σ p . Consequently, the units cancel and all calculations can be conducted directly in terms of x , i.e. camera pixel index.
The data acquisition consisted of taking five camera images per pointer observable (i.e., per lens configuration). A background image, taken with the laser blocked, was subtracted from each.  The resulting image was integrated along the vertical direction y , and normalized to the brightest image obtained in that configuration. The resulting one dimensional probability distribution P (x ), corresponds to the probability of detecting a photon in position x with final polarization |H or |V . With the f x lens in place, this is effectively an x read-out. The expectation values required for the joint weak-value (Eqs. 10 -11) can be obtained as x /σ x = x /σ x , where x = P (x )x dx . For p and d, a similar procedure is followed.
As our first demonstration of the technique, we weakly measure the non-Hermitian product observable π 45 • π V for a range of input states, |I . Specifically, we set the state preparation HWP at an angle of θ/2 and the QWP at 45 • in order to produce the state: For each input state and each image, the expectation values for the joint weak-value (Eqs. 10 -11) were evaluated. The uncertainties were estimated by the standard deviation in the joint weak-value across the five recorded images.
Curves for the real and imaginary parts of the joint weak-value are shown in Fig. 2. The experimental values closely follow the expected curves calculated from the nominal input state |I . However, they do not agree within error. These deviations are likely due to imperfections in the waveplates. These imperfections will also propagate to the alignment of the displacement axes of the BBO crystals since the waveplates are used in the alignment process. Such imperfections have been shown to be the dominant source of systematic error in similar past experiments [6,45]. Nonetheless, the results demonstrate the validity of our proposed technique using a single pointer.

Direct Measurement of the Quantum State
We now move to a more sophisticated demonstration of our technique, the direct measurement of each element of the density matrix of polarization states. Such a direct measurement was introduced in [6,19]. An important advantage of this direct state estimation approach is the number of measurement bases it requires. To obtain a given element of the density matrix, this method requires joint weak-measurements in two complementary bases independently of the dimension of the quantum system. This contrasts with standard quantum state tomography, which requires O(m) bases for an m-dimensional system. The direct estimation approach determines the density matrix of a quantum system element-byelement. One can envision a scenario where the off-diagonal elements of a system's density matrix (known as coherences) are monitored (via direct estimation) as a way to detect decoherence [46]. Unlike in Refs. [6,19], which used two pointers for the joint weak-measurements, here we use only a single pointer (measured in three different bases) for the same task.
A joint weak-measurement of the product π i π 45 • π j , with i, j = H or V (with no postselection) gives the element ρ(i, j) of the density matrix. As shown in [19], the average outcome of a weak measurement without post-selection is the 'weak average' (rather than the weak value), which is equal to the expectation value of the measured observable C = tr [Cρ]. Thus, the direct measurement procedure results in a joint weak-average, ρ(i, j) = 2 tr [π i π 45 • π j ρ]. Therefore, by varying the first and last projectors, the density matrix can be directly determined element-by-element.
To measure the density matrix experimentally, we changed the HWPs settings to scan over the projectors π i and π j . As shown in [19], the final observable in the product can be measured . In frames (b) to (d), we show experimental elements of the density matrix ρ of the polarization states |ψ 1 = cos θ |H + sin θ |V , |ψ 2 = cos θ |H + i sin θ |V and |ψ 3 = 1 √ 2 |H − ie 2iθ |V , respectively. Solid lines correspond to the theory. States following path 1 corresponds to linear polarization (no QWP in the setup), states along paths 2 and 3 were obtained by setting the QWP at 0 • and 45 • , respectively. Error bars are calculated using the standard deviation in the joint weak-value across the five recorded images, as previously employed in Fig. 2. either weakly or strongly. The last PBS implements a strong measurement. For each pair of projectors, we measure the required expectation values in Eqs. 10 -11. We test our direct measurement with three sets of pure polarization states. A general polarization state |ψ = α |H + β |V has density matrix ρ ≡ |ψ ψ| = |α| 2 |H H| + αβ * |H V | + βα * |V H| + |β| 2 |V V | . The states sets are |ψ 1 = cos θ |H + sin θ |V , |ψ 2 = cos θ |H + i sin θ |V and |ψ 3 = 1 √ 2 |H − ie 2iθ |V . These are visualized in the Poincaré sphere in Fig. 3a. For all cases, the first HWP varies the parameter θ scanning the interval [0 • , 180 • ]. The QWP is removed for |ψ 1 , and it is set at 0 • and 45 • for |ψ 2 and |ψ 3 , respectively.
Our results are shown in Fig. 3b -d. The solid lines correspond to the real and imaginary parts of the elements of the theoretical density matrix and the points are the corresponding experimental joint weak-values. The latter should be equal to the real and imaginary parts of the density matrix elements and indeed follow the expected curve for each element. As before, deviations are thought to be the result of systematic errors in the polarization optics. This direct determination of the density matrix demonstrates the utility of our technique for weak measurement applications in quantum information.

Discussion and Conclusion
Before summarizing, we discuss some special cases and extensions of our technique. First, the special case where B = A and the observable is general, i.e., not necessarily a projector. In this case, the unitary evolution (Eq. 5) of the two von Neumann interactions (i.e., Eq. 2) is equivalent to a single unitary given by U A U A = e γA(a † +a) e γA(a † +a) = e 2γA(a † +a) . This corresponds to a unitary of a single von Neumann interaction of the A observable with a doubled interaction strength 2γ. By measuring x, p, and d on the pointer and using Eqs. 10 -11, we can then find the weak value A 2 w . This behaviour can be generalized so that a single von Neumann interaction can be used to measure A N w . To do so, one will need to measure corresponding powers of observables on the pointer, e.g., x N as well as the N -th power of hybrid observables such as d.
In the derivation of our technique, we focused on the case that A and B are projectors. However, for general A and B, measuring the product of projectors is enough to obtain the joint weakvalue BA w . This can be seen if we express A in its spectral decomposition A = α απ α . Here α is an eigenvalue corresponding to the eigenstate |α , and π α = |α α| [47]. Analogously, for B we have B = β βπ β . Therefore, BA ≡ α,β αβπ β π α . Thus, by measuring each of the products π β π α and adding the results, the joint weak-value BA w can be obtained. In summary, our method can be used to weakly measure the product of general incompatible observables A and B.
Our technique is also applicable to observables on separate quantum systems, e.g., A is measured on a first particle and B is measured on a second particle. The standard procedure for weak measurement would couple BA to a single pointer using H = gBAp, which is Hermitian now. This, however, requires a three-particle interaction which is challenging. Our method uses a two-particle interaction on each system while still only using a single pointer. One would first use the standard von Neumann interaction Eq. 2, to couple the pointer to A, then couple the same pointer to B. This double coupling can be challenging to implement. Particularly, in photons one would need an optical nonlinear effect at the single-photon level. As our technique largely preserves the initial quantum state, potential applications include demonstrating contextuality [48][49][50][51], and tests of Leggett-Garg inequalities [33,[52][53][54][55][56][57]. Motivated by the direct measurement of entangled systems [58,59] and entanglement witnesses [60,61], a potential future research direction is the use of this technique for characterizing entanglement.
A possible extension of our technique is measuring the product of m observables of a quantum system, using a single pointer. One would perform subsequent couplings between each of the m observables, and the pointer. The product of the m observables will appear in the expectation value of the m power of the lowering operator of the pointer a m . Performing such a read-out potentially requires full-tomography on the pointer. The advantage of this approach is that only twoparticle interactions are employed, and a single pointer is required for measuring the product of m observables.
As an outline of future work, the technique introduced in this paper can also be extended to general types of pointers such as spin pointers. Indeed, previous work showed that the lowering operator formalism can be extended to spin pointers (e.g., polarization) and spin lowering operators [5]. To measure the product of N observables, one would need to sequentially couple the N observables to the spin pointer. The product of the N observables will appear as a coefficient of the N th excited state of the pointer (just as the product BA appears as a coefficient of the second excited state in Eq. 5). Thus the key requirement will be that one needs N spin levels (i.e., N = 2S + 1, where S is the spin) to measure the product of N observables. This shows that the method trades the resource of N dimensions in the pointer for an N -particle interaction (Eq. 2) or, alternately, N separate pointer systems as in Ref. [5]. Instead of encoding the measurement information in many separate pointers as in Ref. [5,6,22], we encode the information from the measurement in a single highdimensional pointer.
As the main contribution of this paper, we theoretically derived and experimentally demonstrated a method to perform a joint weakmeasurement of two incompatible observables using a single pointer. We then employed this method to directly and individually measure each element of a system's density matrix. Since product observables are ubiquitous in quantum information processing, our technique may be useful for probing and characterizing such processors in situ without substantially disturbing them. Our work optimizes the use of resources needed to perform a joint-weak measurement freeing degrees of freedom of a quantum particle for quantum information tasks. We hope our technique facilitates the use of weak measurement in complex dynam-ics and new studies in the quantum realm.

A Appendix
In this appendix we give an overview of the Fractional Fourier Transform (FrFT). Since its inception, the FrFT has been framed in different contexts: Condon described the FrFT mathematically in [62], while Namias framed it in quantum mechanical terms [41]. This versatility has allowed a broad range of applications ranging from differential equations to quantum optics [63], passing by signal processing, [64,65] and telecommunications [66]. Further details of the FrFT, and code for its numerical implementation can be found in Ref. [64].
The FrFT of order R of a function f (x) is denoted by f R (u) ≡ FrFT(f, R), and it is defined as the following transform: where for all real values of R, except for even integers. For R = 4j, with j an integer, K R (u, u ) = δ(u − u ), and for R = 4j + 2, K R (u, u ) = δ(u + u ).
Particularly important are the R = 0 and R = 1 cases of the FrFT. These correspond, respectively, to the identity and the standard Fourier Transform (FT) operators: f 0 = f (x), and f 1 = FT(f )(p). A useful property of the FrFT is that it is additive in R, i.e., FrFT(FrFT(f, R 2 ), R 1 ) = FrFT(f, R 1 + R 2 ). For example, FrFT(FrFT(f, 1), 1) = FrFT(f, 2) = f (−x) is the parity operator. Therefore, the FrFT is a transformation that interpolates the identity and a normal FT operators. a b Figure 4: The Fractional Fourier Transform (FrFT) allows a continuous transformation of a spatial function to its Fourier transform, passing through other FrFT domains. The spatial function f (x) in this figure is a rectangular function. In frame (a), we are plotting |FrFT(f, R)| 2 , for different FrFT orders R. Some of them are highlighted in orange, each plot is renormalized for better visualization. In frame (b), we highlight the cases of R = 0, 0.5 and 1, which correspond to a transformation to the x, d and p domains respectively. In the x-axis, we have a square function, while on the p-axis we observe a sinc function as we expect to be the standard Fourier Transform of the square function, in our case FrFT for R = 1 corresponding to a π/2 rotation. The d domain is characterized by an equal weight of x and p, and it shows an intensity profile between a square and a sinc functions.
The FrFT is strongly rooted in position-momentum phase-space. Associated with every FrFT order R, there is a different quadrature u that corresponds to a different superposition of position x and momentum p. For dimensionless variables, such domain can be written as u = x cos (Rπ/2)+p sin (Rπ/2). The connection to phase-space is made through the Wigner distribution of a function. The action of a FrFT of order R is a clockwise rotation in phase-space of the Wigner distribution by an angle Rπ/2. A relevant result states that a projection of the Wigner distribution (a marginal of the Wigner distribution) on a domain at an angle Rπ/2, corresponds to the absolute squared of the FrFT of order R. In other words, performing a FrFT on a function allows one to obtain the function representation in a different FrFT domain. Fig. 4 illustrates an example. In the x-axis we have a square function, while on the p-axis we observe a sinc function as we expect to be the standard Fourier Transform of the square function, in our case FrFT for R = 1 corresponding to a π/2 rotation.
It is widely known that a FT can be performed optically by allowing free-space propagation to the far field. The FrFT completes this picture of wave propagation. Diffraction is a continuous process of FrFTs as demonstrated in [67]. Optical implementations of the FrFT have been proposed utilizing lenses, graded-index media, and waveguide arrays [42,44,68,69]. Now we describe the setup to perform a FrFT utilized in our implementation.
Lohman [42] proposed a setup based on a lens that performs a FrFT of order R. Fig. 5 reproduces such setup. The transverse initial position x of an optical mode is fractional Fourier transformed by a lens of focal lenght f , and free-space propagation by a distance z = f tan Rπ 4 sin Rπ 2 before and after the lens f . We utilize this in our experiment for obtaining d.