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.