Eigenvectors in the limit $ \mu\rightarrow 0 $ are not the same as eigenvectors when setting $ \mu=0 $ from the beginning

Question

I would like to find the eigenvectors of a matrix and see what the eigenvectors look like in the limit of $ \mu\rightarrow 0 $:

M = {{0, 1, 0, 0},
     {-ω0^2 - μ^2 - 4*g*x, 0, 0, μ},
     {-μ, 0, 0, 1},
     {0, 0, -ωk^2 + μ^2, 0}};

{s, v} = Eigensystem[-I*M];

One can see that most of the components of the eigenvectors have a $ 1/\mu^c\, (c=1, 2, 3) $ dependence which makes the vector components diverge for $ \mu\rightarrow 0 $. What I do not understand is that if I set $ \mu=0 $ before calculating the eigensystem, the eigenvectors do not diverge at all.

So the question that comes up is, for small $ \mu $ I would expect the eigenvectors to converge to the case in which $ \mu=0 $ from the beginning. Why is this not the case?

Answer 1

As I see it, it is the case. Let me show you how. The "trick" is just like what @SjoerdSmit has mentioned in the comment: normalization.

In addition to M, s, and v, I add two new

{sp, vp} = Eigensystem[-I M /. μ -> 0];

Let me take the first eigenvalue and its corresponding vector as an example:

So the problem now is when $ \mu\rightarrow0 $, are s[[1]] and v[[1]] approaching them? My answer is, yes, they are.

Eigenvalue:

One sees that Out[10] is the same as Out[35].

Eigenvector:

One also sees that Out[43] together with Out[44] expresses the same thing as Out[36].

Answer 2

I'm going to try and show that indeed the question title is false!

• First we'll need to define the original M matrix in a modular form:

  M[μ_] := {
   {0, 1, 0, 0},
   {-ω0^2 - μ^2 - 4 g x, 0, 0, μ},
   {-μ, 0, 0, 1},
   {0, 0, -ωk^2 + μ^2, 0}
  }
  

  This will prove useful later on.

• Next, we are going to define a function that operates on the M matrix and returns its eigenvalues. We are not going to use Eigenalues or Eigensystem, this time.

  evals, the function we are going to define, makes use of the modular definition of M (see above).

  (* evals returns a list of eigenvalues for M *)
  
  evals[μ_] := evals[μ] = Solve[0 == Det[M[μ] - λ IdentityMatrix[4]], λ]
  

Discussion

The eigenvalues of M when $μ\rightarrow0$ can be obtained by first solving for the roots of the characteristic equation in the general case ($μ\neq0$) and then substitute for the value of $μ=0$. On the other hand, solving for the eigenvalues in the specific case where $μ=0$ is straightforward. These two operations are performed by the code below:

(* operators applied after the outcome is derived *)
postOps = MapAt[fAux, \[Bullet], {All, -1, 1, -1}] /* FullSimplify /* PowerExpand

(* first solve for the eigenvalues and then set μ = 0 *)
λ /. evals[μ] /. μ -> 0 // Simplify /* postOps 

{-Sqrt[-4 g x - ω0^2], Sqrt[-4 g x - ω0^2], -I ωk, I ωk}

Similarly,

(* solve for the eigenvalues of M[0] *)
λ /. evals[0]

{-Sqrt[-4 g x - ω0^2], Sqrt[-4 g x - ω0^2], -I ωk, I ωk}

Please note how the two evaluations produce the same output. Hence, it seems that the eigenvalues of M[0] are indeed the same with the respective eigenvalues of M[μ]/.μ->0.

"OK, so, what happens with the eigenvectors, bud?"

First we'll consider the case for M[0] and then we'll discuss the general case M[μ].

In the case of M[0], the eigenvectors can be obtained in the following way:

(* just a vector of auxiliary variables *)
xs = {x1, x2, x3, x4};

(* the eigensystem in the general case for μ *)
eqs[μ_] := eqs[μ] = M[μ].xs - λ xs == 0 // Thread

(* the eigenvalues @ μ = 0 *)
v = λ /. evals[0];

(* post-ops *)
(postOps = {\[Bullet], 1/v} /* 
   MapThread[xs /. Quiet@Solve[#1, xs] /. Thread[xs -> #2] &, \[Bullet]] /* 
   Flatten[\[Bullet], 1]; 

(* the eigenvectors @ μ = 0 *)
eqs[0] /. Transpose[{Thread[λ -> v]}] // Simplify /* postOps

{ 
     {-(1/Sqrt[-4 g x - ω0^2]), 1, 0, 0}, 
     {1/Sqrt[-4 g x - ω0^2], 1, 0, 0}, 
     {0, 0, I/ωk, 1},
     {0, 0, -(I/ωk), 1}
    }

It so happens that these are the eigenvectors one obtains when using the built-ins eg Eigensystem[M[0]].

Now, it's not clear to me if this is actually a part of the original question but I feel I should address what happens when the calculations are performed for the general case ($μ\neq0$) and then we impose the restriction.

In that case the eigenvectors are different from the ones we previously calculated (I won't present those calculations here since they are essentially the same with the ones for the special case presented above; the only difference is that the eigenvalues used should be the ones corresponding to the general case).

My guess as to why this happens has to do with the matrix M[μ] itself. Notice how it is possible to write the general-case matrix M[μ] as a sum of two matrices:

M[μ] == M[0] + {{0, 0, 0, 0}, 
                {-μ^2, 0, 0, μ}, 
                {-μ, 0, 0, 0}, 
                {0, 0, μ^2, 0}}

In the special case where $μ=0$ the second matrix disappears and we obtain the familiar results as in above.

In the general case, where $μ\neq0$ the simple-case matrix M[0] is perturbed by the non-zero second matrix (see above) in such ways that essentially it is transformed into a completely different matrix with completely distinct characteristics from the former one.

That's as far as I can go. I would be really interested to know more.

---

Auxiliary code

• The Bullet operator (obtained from here)

  \[Bullet] /: f_[pre___, \[Bullet], post___] := With[{n=Length[List@pre], m=Length[List@post]}, 
    Curry[f, Join[Range[n], {n + m + 1}, Range[m] + n]][pre, post]
   ]
  

• The fAux auxiliary function

  (* just an error message *)
  
  fAux::herr = "Something went wrong. Please, look into it!";
  

  and

  (* relevant for input that match -Sqrt[x__] *)
  
  fAux[pattn : -Sqrt[x__]] := With[{sim = Simplify[x]},
    Which[
      MatchQ[sim, Power[p__, 2]], -sim /. Power[p__, 2] :> p,
      True, (Message[h::herr]; pattn)
     ]
   ]
  

  and

  (* relevant for input that match Sqrt[x__] *)
  
  fAux[pattn : Sqrt[x__]] := With[{sim = Simplify[x]},
    Which[
      MatchQ[sim, Power[p__, 2]], sim /. Power[p__, 2] :> p,
      True, (Message[h::herr]; pattn)
     ]
    ]
  

Answer 3

The premise of the question is wrong. The eigenvectors for generic $\mu$ do become those for $\mu=0$ in the $\mu\to0$ limit, but with a qualification. To illustrate my point, take the $1\times 1$ matrix $$ M=1+\mu $$

The eigenvector of this matrix is just any non-zero complex number. I will take, just because I want to, $v=i/\mu$. It is clear that in the $\mu\to0$ limit the eigenvector blows up. But the eigenvectors are only defined up to a global multiplicative coefficient, so I could also take $v=i$. This no longer blows up. So the reason the original eigenvector was singular at $\mu\to0$ was because I chose a bad normalisation.

Now consider the same matrix, but taking $\mu=0$, $$ M=1 $$

Again, any non-zero complex number is an eigenvector. I will now choose $v=3$. This is different from the $\mu\to0$ limit in the previous case. Is this a problem? Well, again, no: eigenvectors are only defined up to a global coefficient, so I could redefine $v\to iv/3$, and now this would match the previous case.

More generally, if you have a matrix $M(\mu)$, its eigenvectors at $\mu\to\mu_\star$ will only match those at $\mu=\mu_\star$ up to a global coefficient. (The exception is when an eigenspace becomes degenerate at $\mu_\star$, but this does not happen in the example in the OP, so we will not worry about that).

Now this is precisely what happens in the case of OP. To compare the eigenvectors at $\mu\to0$ with those at $\mu=0$, one must specify some convention for the global coefficient; otherwise the comparison is just meaningless. For example, one may normalise the vectors to unit norm (but then the global phase is still ambiguous). In order to speed up the computations I will take $g=\pi$, $x=e$, $\omega_0=\gamma$ and $\omega_k=\zeta(3)$. OP is invited to repeat the computations above with symbolic parameters, but the calculations will take forever. With this, let

M[μ_] := {{0, 1, 0, 0}, {-ω0^2 - μ^2 - 4 g x, 0, 0, μ}, {-μ, 0, 0, 1}, {0, 0, -ωk^2 + μ^2, 0}} /. {g -> π, x -> E, ω0 -> Zeta[3], ωk -> EulerGamma, μ -> 0}

Eigenvectors[M[0]] // Map[Normalize] // N // Chop
Eigenvectors[M[1/1000]] // Map[Normalize] // N // Chop

$$ \left( \begin{array}{cccc} -0.16 i & 0.98 & 0 & 0 \\ +0.16 i & 0.98 & 0 & 0 \\ 0 & 0 & 0.\, -0.86 i & 0.49 \\ 0 & 0 & 0.\, +0.86 i & 0.49 \\ \end{array} \right)\qquad \left( \begin{array}{cccc} -0.16 & 0.98 i & -0.00002 i& 1.56\cdot10^{-6} \\ -0.16 & -0.98 i & 0.00002 i& 1.56\cdot10^{-6} \\ 0.00001 & 8.1\cdot10^{-6} i & -0.86 i & 0.49 \\ 0.00001 & -8.1\cdot10^{-6} & -0.86 i & 0.49 \\ \end{array} \right) $$

OP may appreciate that these two matrices of eigenvectors are basically the same (up to factors of $i$ which, again, are arbitrary). The agreement becomes better for smaller values of $\mu$, and the matrices become identical at $\mu\to0$.