Why do Python NumPy/SciPy numerical eigenvectors act incorrect compared to MATLAB when there is at least one eigenvalue degenerate point?

Question

For a physics problem, I numerically solve an eigenvalue equation with a square matrix (that is usually complex-valued, but is guaranteed to have real eigenvalues by construction). I have a 2D parameter space (x,y), over which I solve this equation at every point. Then I take the dot product of an eigenvector and the discretized time derivative of another eigenvector along a trajectory in (x,y). For example, if my matrix is 4 x 4, and I get eigenvectors A, B, C, D with eigenvalues a, b, c, d, I am talking about (A).(dB/dt) (which is usually complex-valued). Then I calculate its complex argument/angle.

As long as my (x,y) points are sufficiently close enough (i.e. 'continuous'), and as long as I know that the eigenvalues of each eigenfunction considered are unique at every point along the trajectory (i.e. a =/= b at for all (x,y) making the trajectory), I can expect this angle to be a continuous quantity as well, without any sharp jumps (as in the first figure). However, when using Numpy/Scipy it seems as if this 'continuous' quality is ruined by a degenerate/repeating eigenvalue even if it is not on the trajectory. That is, if my trajectory comprises {(x1,y1),(x2,y2),(x3,y3),(x4,y4),(x5,y5)}, with unique eigenvalues at each point (i.e. a1 =/= b1, a2 =/= b2, ...), it seems as if that a9 = b9 at (x9, y9) causes trouble even if (x9, y9) is not on the trajectory considered. I deduced this by changing my matrix so that the eigenvalues of each eigenfunction are unique at each point (i.e. I forced a9 =/= b9). Doing so gave me a 'continuous' angle. On the contrary, both cases are 'continuous' in MATLAB, as expected in theory.

It seems as if the issue is with the complex phases assigned to the numerical eigenfunctions. I used np.unwrap to delve into this more (last figure below), but it wasn't too fruitful. I have spent weeks trying to understand why, to no avail. Any suggestions would be appreciated.

Expected behavior of complex argument along a trajectory without any repeating eigenvalues (both Python and MATLAB shows this when eigenvalues are unique everywhere):

Observed behavior when using numpy.linalg.eig or scipy.linalg.eig:

Observed behavior when using numpy.linalg.eig AND np.unwrap - there aren't supposed to be jumps: