Как удалить разрывы из комплексного угла компонентов собственного вектора NumPy?

Я использую linalg.eig NumPy для квадратных матриц. Мои квадратные матрицы являются функцией 2D-области, и я смотрю на комплексные углы ее собственных векторов вдоль параметризованной окружности в этой области. Пока рассматриваемый мной путь является гладким, я ожидаю, что комплексные углы компонентов каждого собственного вектора будут гладкими. Однако в некоторых случаях это не относится к Python (хотя это относится к другим языкам программирования). Для параметра M=0 (некоторый аргумент в моей матрице, который появляется на его диагонали) у меня есть компоненты, которые выглядят так:

когда они должны в идеале выглядеть так (M=0.1):

Что я пробовал:

• Я проверил, что в обоих случаях матрицы эрмитовы.
• Когда я использую linalg.eigh, M=0.1 становится прерывистым, а M=0 иногда становится непрерывным.
• Использование np.unwrap ничего не дало.
• Разница между фазами компонентов (т.е. np.angle(v1-v2) для собственного вектора v=[[v1],[v2]]) плавная / непрерывная, но это не то, что я хочу.
• Исправление семени NumPy перед решением ничего не дало для разных значений семени. Например: np.random.seed(1).

Что еще я могу сделать? Я пытаюсь использовать Sympy eigenvects только потому, что у меня заканчиваются варианты, и я задал еще один вопрос, касающийся другого потенциального подхода здесь: Как заставить первый компонент собственных векторов NumPy быть действительным?. Но я не знаю, что еще я могу попробовать.

Вот минимальный рабочий пример, который отлично работает в записной книжке Jupyter:

import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

M = 0.01; # nonzero M is okay
M = 0.0; # M=0 causes problems

def matrix_generator(kx,ky,M):
    a = 2.46;    t = 1;    k = np.array((kx,ky));
    d1 = (a/2)*np.array((1,np.sqrt(3)));d2 = (a/2)*np.array((1,-np.sqrt(3)));d3 = -a*np.array((1,0));
    sx = np.matrix([[0,1],[1,0]]);sy = np.matrix([[0,-1j],[1j,0]]);sz = np.matrix([[1,0],[0,-1]]);
    hx = np.cos(k@d1)+np.cos(k@d2)+np.cos(k@d3);hy = np.sin(k@d1)+np.sin(k@d2)+np.sin(k@d3);
    return -t*(hx*sx - hy*sy + M*sz)

n_segs = 200; #number of segments in (kx,ky) loop
evecs_along_loop = np.zeros((n_segs,2,2),dtype=float)
# parameterize circular loop
kx0 = 0.5; ky0 = 1; r1=0.2; r2=0.2; 
a = np.linspace(0.0, 2*np.pi, num=n_segs+2)
kloop=np.zeros((n_segs+2,2))
for i in range(n_segs+2):
    kloop[i,:]=np.array([kx0 + r1*np.cos(a[i]), ky0 + r2*np.sin(a[i])]) 
    
# assign eigenvector complex angles
for j in np.arange(n_segs):
    np.random.seed(2)
    H = matrix_generator(kloop[j][0],kloop[j][1],M)
    eval0, psi0 = LA.eig(H)
    evecs_along_loop[j,:,:] = np.angle(psi0)
    
# plot eigenvector complex angles
for p in np.arange(2):
    for q in np.arange(2):
        print(f"Phase for eigenvector element {p},{q}:")
        fig = plt.figure()
        ax = plt.axes()
        ax.plot((evecs_along_loop[:,p,q]))
        plt.show()

Пояснение к комментарию anon01: Для M=0 примерная матрица с некоторым значением (kx,ky) будет выглядеть так:

a = np.matrix([[0.+0.j, 0.99286437+1.03026667j],
 [0.99286437-1.03026667j, 0.+0.j]])

Для M =/= 0 диагональ будет отличной от нуля (но реальной).