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

Question

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 = 2t0(Cos[2((1/2)kx*a)] + 
      2Cos[((1/2)kxa)]Cos[(kyaSqrt[3]/2)]) + e1;
h1 = -2Sqrt[3]t2Sin[((1/2)kxa)]Sin[(kyaSqrt[3]/2)] + 
   21 It1(Sin[2((1/2)kxa)] + 
      Sin[((1/2)kxa)]Cos[(kyaSqrt[3]/2)]);
h1star = -2Sqrt[3]t2Sin[((1/2)kxa)]Sin[(kyaSqrt[3]/2)] - 
   21 It1(Sin[2((1/2)kx*a)] + 
      Sin[((1/2)kxa)]Cos[(kyaSqrt[3]/2)]);
h2 = 2t2(Cos[2((1/2)kxa)] - 
      Cos[((1/2)kxa)]Cos[(kyaSqrt[3]/2)]) + 
   2Sqrt[3]1 It1Cos[((1/2)kxa)]Sin[(kyaSqrt[3]/2)];
h2star = 2*
    t2(Cos[2((1/2)kxa)] - 
      Cos[((1/2)kxa)]Cos[(kyaSqrt[3]/2)]) - 
   2Sqrt[3]1 It1Cos[((1/2)kxa)]Sin[(kyaSqrt[3]/2)];
h11 = 2t11Cos[2((1/2)kx*a)] + (t11 + 3t22)Cos[((1/2)kxa)]*
    Cos[(kyaSqrt[3]/2)] + e2;
h22 = 2t22Cos[2((1/2)kx*a)] + (3t11 + t22)Cos[((1/2)kxa)]*
    Cos[(kyaSqrt[3]/2)] + e2;
h12 = Sqrt[3](t22 - t11)Sin[((1/2)kxa)]Sin[(kyaSqrt[3]/2)] + 
   41 It12
    Sin[((1/2)kxa)](Cos[((1/2)kxa)] - Cos[(kyaSqrt[3]/2)]);
h12star = 
  Sqrt[3](t22 - t11)Sin[((1/2)kxa)]Sin[(kyaSqrt[3]/2)] - 
   41 It12*
    Sin[((1/2)kxa)](Cos[((1/2)kxa)] - Cos[(kya*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.

Answer 1

Quoting @LouisB's comment:

Tell Map to operate on all elements that can be accessed with 3 subscripts, like this: Map[Normalize[#/SelectFirst[#, # != 0 &]] &, evecs, {3}] The SelectFirst is used to avoid divide-by-zero.