(* basis states {upA-upB,upA-downB,downA-upB,downA-downB} *)

(* spin operator ma along the vector a for qubit a and mb along b for qubit b *)

ma={{az,ax-I*ay},{ax+I*ay,-az}}; mb={{bz,bx-I*by},{bx+I*by,-bz}};

m=ArrayFlatten[KroneckerProduct[ma,mb]];

MatrixForm[m]

(* Bell state *)

psi={0,1,-1,0}/Sqrt[2];

(* expectation value of the correlator of the spins along a and b *)

psi.m.psi//FullSimplify

(* this is the angle between the unit vectors a and b *)