Assume the initial state is all B0 at time t=0.  Then:

ampB1at0 = 1/Sqrt[2]

ampB2at0 = 1/Sqrt[2]

Now propogate the two states in time:

ampb1[t_] = Simplify[ampB1at0 Exp[-I(m1- I gamma)t]]

ampb2[t_] = Simplify[ampB2at0 Exp[-I(m2- I gamma)t]]

After a time t, the B0 amplitude is the sum of B1 and B2:

amp1 =Simplify[1/Sqrt[2] (ampb1[t]+ampb2[t])/.
{m1->m+dm,m2->m-dm}]

amp2=Simplify[Expand[amp1/(Exp[-gamma t - I m t]/2)]] *
	(Exp[-gamma t - I m t]/2)

ampB0[t] = amp2[[2]] Cos[dm t]

Get the probability after a time t:

probB0 = ampB0[t]*(ampB0[t]/.{I->-I,-I->I})

Plug in the actual parameters:

probactual=(probB0/.{t->tau/gamma})/.{dm->3.1*10^(-10),
	gamma->6.58*10^(-22) /(1.56*10^(-12))}

Get just a pure decay spectrum:

probzero=(probB0/.{t->tau/gamma})/.{dm->0,
	gamma->6.58*10^(-22) /(1.56*10^(-12))}

Plot the probability of having a B0:

plot1=Plot[probactual,{tau,0,2}]

plot2=Plot[probzero,{tau,0,2},
	PlotRange->All]

And finally the B0bar distribution:

plot3=Plot[probzero-probactual,{tau,0,2},
	PlotRange->All]

Show all three:

Show[plot1,plot2,plot3]