Simulation of a measurement

Quantum Mechanics
                       K. Konishi, G. Paffuti

Solve Problem 2.3 using Mathematica. Write a simulation for the proposed set of measurements. In a finite number of measurements the outcomes will be distributed with a binomial distribution. Compare the results of the simulation with the expected distribution.

We report here the text of the problem 2.3:

Text

A two-level system is described by the Hamiltonian,

H = (Subscript[E, 0])(({{1, 0}, {0, -1}}))

We consider in this system the measurement of two variables described by the operators A and B,

A =(({{0, -i}, {i, 0}}));  B = (({{2, -(Sqrt[2])i}, {(Sqrt[2])i, 1}})).

Solution with Mathematica

ClearAll["Global`*"];
Off[General::"spell",General::"spell1"]

Let us define the 3 operators H, A, B as two by two matrices in Mathematica language:

H = E0{{1,0},{-1,0}};
A = {{0,-I},{I,0}};
B = {{2,-(Sqrt[2])I},{(Sqrt[2])I,1}};

Question 1)

Eigenvalues and eigenvectors of A are obtained diagonilizing the matrix with the Mathematica function Eigensystem

{λA,ψA}=Eigensystem[A]

Normalized eigenvectors are obtained with Norm

ψA = ψA/(Norm/@ψA)

An analogous computation gives eigenvalues and eigenvectors of B

{λB,ψB}=Eigensystem[B]

ψB=ψB/(Norm/@ψB)

The largest eigenvalue of B  (b1) is 3.

Question 2)

A and B are not compatible, i.e. they do not commute:

A.B == B.A

nor they are conserved, i.e. they do not commute with H:

Simplify[{A.H == H.A, B.H== H.B},Assumptions-> E0>0]

We used the function Simplify as otherwise Mathematica consider E0 as a free variable, and for E0 = 0 clearly the operators commute.

Question 3)

The measurement of B can give only one of the eigenvalues of this variable. For an arbitrary state ψ the probability of obtaining the k-eigenvalue eig[[k]] is given by the expression

prob =|⟨eigvec[[k]] | ψ ⟩((|)^2)

For operators A and B respectively:

probA[k_,psi_]:= Abs[Conjugate[ψA[[k]]].psi]^2
probB[k_,psi_]:= Abs[Conjugate[ψB[[k]]].psi]^2

If a measurement of A gives 1 as a result the state is projected into the corresponding eigenstate (eigenvecA[[2]]). The probabilities for the two possible values of B are

Simplify[probB[#1,ψA[[2]]]&/@{1,2}]

The sum is 1, as expected:

Simplify[Total[%]]

Question 4)

After the first measurement of A, the system is projected in the state ψA[[2]]. Let us call ψ this state

ψ = ψA[[2]];

a)
The measurement of A immediately after gives again A=1

b)
At time t the state evolves. Vectors (1,0) and (0,1) are eygenstates of H and we can write the evolution as

evol[psi_,t_]:= {psi[[1]] Exp[-I E0 t] , psi[[2]] Exp[+I E0 t]}

The probability of getting again A=1 is

FullSimplify[Abs[Conjugate[ψ].evol[ψ,t]]^2,Assumptions->{ t > 0,E0>0}]

c)
The probability of getting b1 measuring B on the state ψ is

|(ψ | ψB[[1]])((|)^2).

After the measurement the system is in the state ψB[[1]] and the (relative) probability of measuring now A=1 is |(ψB[[1]] | ψ)|^2. In this way the probability of the whole sequence is:

|(ψB[[1]] | ψ)((|)^2)

probAequal1ForBequal1=Simplify[Abs[Conjugate[ψB[[1]]].ψ]^2*Abs[Conjugate[ψ].ψB[[1]]]^2]

or

N[%]

d)
The measurement of B acts as a filter, if we accept both outcomes of the measurement the total probability of measuring A =1 after measuring B is:

probAequal1ForAllB=Simplify[Abs[Conjugate[ψB[[1]]].ψ]^2*Abs[Conjugate[ψ].ψB[[1]]]^2 + 
Abs[Conjugate[ψB[[2]]].ψ]^2*Abs[Conjugate[ψ].ψB[[2]]]^2]

Note that this probability is not 1, i.e. it is different from the result of measurement  a). This is a crucial point. If the measurement of B is not performed the state ψ can always be analyzed in terms of eigenstates of B but for the final measurement of A=1 the probability is given by

probA =|⟨ψ | ψ ⟩((|)^2)=1= |⟨ψ |b1⟩⟨b1 ψ ⟩ + ⟨ψ |b2⟩⟨b2 ψ ⟩((|)^2)

i.e. we have an interference of amplitudes. 
In the present case (d) the measurement of B is done and the corrisponding result is b1 or b2 (exclusive) with relative probability

probb1 =|⟨b1 | ψ ⟩((|)^2);  probb2 =|⟨b2 | ψ ⟩((|)^2);

The state after measurement is either |b1⟩ or |b2⟩, with no inteference, and the (relative) probability of measuring again A=1 in the two distinct situations is

probA[b1] =|⟨ψ | b1 ⟩((|)^2);  probA[b2] =|⟨ψ | b2 ⟩((|)^2);

If we do not distinguish between the two cases the total probability is

probA =probb1*probA[b1]+probb2 * probA[b2] ;

this is the expression computed above.

A simulation of mesurements

A simple way of generating an event with probability ξ is to extract a random number r, with uniform distribution in [0,1] and accept the event if r < ξ. Our variables, A and B, have only two possible outcomes and the effect of a mesurement on a given state ψ is to project on an eigenstate vec[[k]] with probabilities given by Q.M.. A simple way to do this is given by the functions measA and measB below:

measA[ψ_]:=If[RandomReal[]≤((Abs[Conjugate[ψA[[1]]].ψ])^2),ψA[[1]],ψA[[2]]];
measB[ψ_]:=If[RandomReal[]≤((Abs[Conjugate[ψB[[1]]].ψ])^2),ψB[[1]],ψB[[2]]];

To fix the sequence of random numbers (for reproducibility) let us use SeedRandom

SeedRandom[100];

We simulate the experiment preparing a beam of  identical systems (1000 for example) in the state ψA[[2]]:

numSystems = 1000;
PsiBeam = Table[ψA[[2]],{numSystems}];

We can measure B on each system:

beamAfterB =  measB[#]&/@ PsiBeam;

To analyse the results we write a function which select a given beam using the eigenvectors eigvec of a given variable

selectBeam[beam_,eigvec_]:= {Select[beam,#== eigvec[[1]]&],Select[beam,#== eigvec[[2]]&]};

In particular for the system obtained after a measurement of B we have

beamAfterBSelected = selectBeam[beamAfterB,ψB];

The length of the two selected list of vectors gives the number ofoutcomes for each value of B:

Length/@beamAfterBSelected

In this particular case we obtained 973 times B=b1, and 27 times B = b2.

Measurements

We can now simulate the 3 different measurements. The measurement of A on ψ gives evidently always the value 1, let us check this point

beamAfterA = measA[#]&/@PsiBeam;
Length/@selectBeam[beamAfterA,ψA]

We already have produced the result of the measurement of B, i.e. the list beamAfterB. Suppose now to measure again A only on the positive outcomes for B=b1.

The corresponding systems are contained in the selected sublist beamAfterBSelected[[1]] and the measurement of A gives

(* This measures A on the particles on which B = b1 has been measured *)
beamAfterA = measA[#]&/@beamAfterBSelected[[1]]; 
(* Here we select the outcomes *)
beamAfterASelected = selectBeam[beamAfterA,ψA];
foundA=Length/@beamAfterASelected

foundAforBequal1 = foundA[[2]];

The sum of the numbers above is less than the total number of systems beacause the outcomes B=b2 have been excluded.

The second type of experiment accepts all the results of the first measurement of B, i.e. measA acts on the whole list beamAfterB, the result is

beamAfterA = measA[#]&/@beamAfterB;
beamAfterASelected = selectBeam[beamAfterA,ψA];
foundA=Length/@beamAfterASelected

foundAforAllB= foundA[[2]];

The probability of expected  A=1 measurements in the two cases has been computed above and the expected outocomes are, in the two cases:

numSystems*N[{probAequal1ForBequal1,probAequal1ForAllB}]

to be compared with

{foundAforBequal1,foundAforAllB}

Note that in general these numbers are different from the "experimental" ones, only in the limit of infinitely many systems the Q.M. prediciton is reached, for finitely many systems there are statistical fluctuations. These fluctuations have nothing to do with Q.M., they reflect the fact that the frequency definition of probability  (num favorable cases) / (num possible cases) tends to the true probability only asimptotically.

Distribution of the results

To clarify the last point let us repeat many times ( 500) the experiment of measuring B with a beam of 250 systems in the state ψ.

nsample= 500;nbeam = 250;

PsiBeam = Table[ψA[[2]],{nbeam}];

We make a table with the numbers of outocomes for B=b1 and B=b2 and extract only the first result, with [[All,1]], i.e. the numbers of favorable outcomes for each experiment.

numberb1Results=Table[Length/@selectBeam[ measB[#]&/@ PsiBeam,ψB],{nsample}][[All,1]];

The statistical analysis of the data (mean and variance) is:

N[{Mean[numberb1Results],Variance[numberb1Results]}]

The probablity of obtaining b1 has been computed above and is

(Subscript[p, 1])=|⟨ψB[[1]]|ψA[[2]]⟩((|)^2) = (3+2 (Sqrt[2]))/6,

p1 = Simplify[probB[1,ψA[[2]]]]

The probability of obtaining k success on n trials is given by the binomial distribution

P[k,n] = (({{n}, {k}}))(p1^k) (((1-p1))^(n-k))

it is well known, and easy to verify, that

Mean[k] = n p1;  Mean[(k^2)] - (((Mean[k]))^2)= n p1(1-p1)

In our case

p1= (1/6) (3+2 (Sqrt[2]));
N[{p1*nbeam,p1*(1-p1)*nbeam}]

We can compare the distribution of obtained positive outcomes (i.e. B=b1) with the expected binomial distribution

Show[
Histogram[numberb1Results,{nbeam-14.5`+Range[1,15]}],ListPlot[N[Table[{j,nsample Binomial[nbeam,j] (p1^j) (((1-p1))^(nbeam-j))},{j,nbeam-15,nbeam}]],PlotStyle->PointSize[0.02`]],AxesLabel-> {"Pos. outcomes","distribution"}]