matrix – How to make first component of eigenvectors on grid real and normalize using Map and Normalize?

I am trying to use the following approach: How get eigenvectors without phase jump?
to apply the same conditions to an array of complex-valued Mathematica eigenvectors.

The matrix solved for is a 3 x 3 matrix. I get my eigenvectors at points on a two-dimensional N x N grid, and so my array of eigenvectors has dimensions N x N x 3 x 3. The documentation for Map allows one to specify levelspec, but so far, I have not been able to successfully impose the conditions of normalization and real-first-component to my eigenvectors. I am still having trouble understanding multi-dimensional arrays and indexing. Ideally, I would have an N x N x 3 x 3 array of eigenvectors whose first component is real, and all eigenvectors normalized. With my limited understanding, I cannot seem to find a way to do this without flattening my array and losing information on the underlying N x N grid.

Almost every time, I end up with an array that is flattened (so, 100x100x3 becomes 100, etc). I tried working with just one eigenvector as in the referenced question, but I really want to preserve the underlying N x N structure, as I want to be able to associate each eigenvector to its respective 2D grid’s coordinates without losing that information by Flattening.

I have attached my code.

(* First generate matrix and get eigenvector array. *)
a = 3.19; 
e1 = 1.046; e2 = 2.104; t0 = -0.184; t1 = 0.401; t2 = 0.507; t11 = 
0.218; t12 = 0.338; t22 = 0.057;

h0 = 2*t0*(Cos(2*((1/2)*kx*a)) + 
      2*Cos(((1/2)*kx*a))*Cos((ky*a*Sqrt(3)/2))) + e1;
h1 = -2*Sqrt(3)*t2*Sin(((1/2)*kx*a))*Sin((ky*a*Sqrt(3)/2)) + 
   2*1 I*t1*(Sin(2*((1/2)*kx*a)) + 
      Sin(((1/2)*kx*a))*Cos((ky*a*Sqrt(3)/2)));
h1star = -2*Sqrt(3)*t2*Sin(((1/2)*kx*a))*Sin((ky*a*Sqrt(3)/2)) - 
   2*1 I*t1*(Sin(2*((1/2)*kx*a)) + 
      Sin(((1/2)*kx*a))*Cos((ky*a*Sqrt(3)/2)));
h2 = 2*t2*(Cos(2*((1/2)*kx*a)) - 
      Cos(((1/2)*kx*a))*Cos((ky*a*Sqrt(3)/2))) + 
   2*Sqrt(3)*1 I*t1*Cos(((1/2)*kx*a))*Sin((ky*a*Sqrt(3)/2));
h2star = 2*
    t2*(Cos(2*((1/2)*kx*a)) - 
      Cos(((1/2)*kx*a))*Cos((ky*a*Sqrt(3)/2))) - 
   2*Sqrt(3)*1 I*t1*Cos(((1/2)*kx*a))*Sin((ky*a*Sqrt(3)/2));
h11 = 2*t11*Cos(2*((1/2)*kx*a)) + (t11 + 3*t22)*Cos(((1/2)*kx*a))*
    Cos((ky*a*Sqrt(3)/2)) + e2;
h22 = 2*t22*Cos(2*((1/2)*kx*a)) + (3*t11 + t22)*Cos(((1/2)*kx*a))*
    Cos((ky*a*Sqrt(3)/2)) + e2;
h12 = Sqrt(3)*(t22 - t11)*Sin(((1/2)*kx*a))*Sin((ky*a*Sqrt(3)/2)) + 
   4*1 I*t12*
    Sin(((1/2)*kx*a))*(Cos(((1/2)*kx*a)) - Cos((ky*a*Sqrt(3)/2)));
h12star = 
  Sqrt(3)*(t22 - t11)*Sin(((1/2)*kx*a))*Sin((ky*a*Sqrt(3)/2)) - 
   4*1 I*t12*
    Sin(((1/2)*kx*a))*(Cos(((1/2)*kx*a)) - Cos((ky*a*Sqrt(3)/2)));
H(kx_, ky_) = {{h0, h1, h2}, {h1star, h11, h12}, {h2star, h12star, 
    h22}};

numsquares = 100;
kxs = N(Subdivide(-0.1, 0.1, numsquares)); kys = 
 N(Subdivide(-0.1, 0.1, numsquares));
evecs = Table(
   Eigenvectors(H(kxs((i)), kys((j)))), {i, 1, numsquares}, {j, 1, 
    numsquares});
Dimensions(wfcs)

(* Try to impose condition in several ways. *)
evecs2 = Map(Normalize(#/#((1))) &, evecs);
Dimensions(evecs2)

evecs2 = Map(Normalize(#/#((1))) &, evecs((All,All,1)));
Dimensions(evecs2)

Any thoughts? Thanks.