(*Hamiltonian *)

H =((-ℏ)/(2m))(∇+A)·(∇+A)+ μ (Subscript[B, z])

S = Simplify[#,TimeConstraint->1]&;

(* quantities we're trying to find *)

F =((-(ℏ^2))/(2m)) ∇∘A + A∘A;
G=-((ℏ^2)/m) A.∇

Clear[B1,z]
B1= {β x - (γ/2)x z, -β y -(γ/2)y z,B - (γ/4)((x^2)+(y^2))+(γ/2)(z^2)};
b= (B1/(Sqrt[B1.B1]));
zh= {0,0,1};
δ [i_]:= IdentityMatrix[i];
ϵ= LeviCivitaTensor[3]//Normal;
σ=Table[PauliMatrix[k],{k,3}];
U = (1/(Sqrt[2]))((Sqrt[1+b.zh])δ[2]-I((Cross[b,zh].σ )/(Sqrt[1+ b.zh]))  )//FullSimplify;
Udagger=  (1/(Sqrt[2]))((Sqrt[1+b.zh])δ[2]+I((Cross[b,zh].σ )/(Sqrt[1+ b.zh]))  )//FullSimplify;

(U.Udagger)//FullSimplify

A=Table[U.( Grad[Udagger,{x,y,z}][[All,All,i]]//Simplify),{i,3}];

F=-((ℏ^2)/(2m))(  (∑_(i=1)^3)A[[i]].A[[i]] +(∑_(i=1)^3)D[A[[i]],{x,y,z}[[i]]] );

(* A trick for calculating series expansion of three parameters, rescaling all of them and doing series expansion for t, credit to Nikhil) *)

Simplify[Series[F/.{x-> t x,y-> t y,z-> t z},{t,0,1}]//Normal,Assumptions-> B>0&& β>0&& (β^2)/B1 > γ/2&& γ>0]/.{t-> 1}//MatrixForm

(* For evaluating G, we can expand A first *)

Aapprox= Simplify[Series[A/.{x-> t x,y->t y,z->t z},{t,0,1}]//Normal,Assumptions-> B>0&& β>0&& (β^2)/B1 > γ/2&& γ>0]/.{t-> 1};

Aapprox[[1]]//MatrixForm//S
Aapprox[[2]]//MatrixForm//S
Aapprox[[3]]//MatrixForm//S

G:=(-((ℏ^2)/(2m)))2(∑_(i=1)^3) Aapprox[[i]](∂_({x,y,z}[[i]]))#&