# System of linear differential equations

Just for fun I would like to enter the Dirac equation into Mathematica to see what shape of solution it will return. So I define the Dirac matrices:

id = IdentityMatrix[4];
Subscript[\[Gamma], 0] = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}};
Subscript[\[Gamma], 1] = {{0, 0, 0, 1}, {0, 0, 1, 0}, {0, -1, 0, 0}, {-1, 0, 0, 0}};
Subscript[\[Gamma], 2] = {{0, 0, 0, -I}, {0, 0, I, 0}, {0, I, 0, 0}, {-I, 0, 0, 0}};
Subscript[\[Gamma], 3] = {{0, 0, 1, 0}, {0, 0, 0, -1}, {-1, 0, 0, 0}, {0, 1, 0, 0}};


Dirac differential operator:

Subscript[Dop, m_][Y_] := Hold[(I Subscript[\[Gamma], 0].D[in, t] - I Subscript[\[Gamma], 1].D[in, x] - I Subscript[\[Gamma], 2].D[in, y] - I Subscript[\[Gamma], 3].D[in, z] + m id.in)] /. in :> Y // ReleaseHold


And a Dirac spinor to solve for:

\[Psi] = {u1[t, x, y, z], u2[t, x, y, z], v1[t, x, y, z], v2[t, x, y, z]};


Finally, I ask for the solution:

DSolve[Table[Subscript[Dop,m][\[Psi]][[it]] == 0, {it, 1, 4}], \[Psi], {t, x, y, z}]


But, unfortunately, Mathematica just instantly returns back the input without even attempting to solve it. Does this mean Mathematica cannot solve this? Or did I do something wrong?

Maybe it might help to say how the solutions are usually obtained: First one transforms the operator to a second order one (by acting on it with its complex conjugate). The second order operator then has plane wave solutions, and therefore so does the original. Then it is solved in Fourier space (like anything involving waves).

EDIT:

The second order differential operator is obtained by first defining the complex conjugate of the original:

Subscript[cDop, m_][Y_] := Hold[(-I Subscript[\[Gamma], 0].D[in, t] + I Subscript[\[Gamma], 1].D[in, x] + I Subscript[\[Gamma], 2].D[in, y] + I Subscript[\[Gamma], 3].D[in, z] + m id.in)] /. in :> Y // ReleaseHold


And then applying them consecutively:

Subscript[cDop,m][Subscript[Dop,m][\[Psi]]]


where this should turn out to be the Klein-Gordon (plane wave) operator

KG:=(D[#,{t,2}]-D[#,{x,2}]-D[#,{y,2}]-D[#,{z,2}]+m^2 #)&


in each of the four components.

I am basically wondering whether the application of the second complex conjugate operator projects the solution space onto plane waves. And whether the original first order operator might have additional non-plane wave solutions as well?

• Could you add the second-order equations to the question? Thanks. Nov 7, 2016 at 22:18
• Added the second order operator. Nov 7, 2016 at 22:31
• I didn't see any need to use your edit containing the operator product because the Fourier approach already works for each factor separately.
– Jens
Nov 8, 2016 at 4:48

This is a system of (linear) partial differential equations which Mathematica indeed seems to be unable to handle with DSolve. However, the linearity of the problem guarantees that the Fourier approach will provide a basis for all possible solutions.

To follow that approach, you just have to replace $\psi$ by a plane-wave ansatz, and then the solution method would go like this (I removed subscripts to improve readability):

id = IdentityMatrix[4];
γ0 = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}};
γ1 = {{0, 0, 0, 1}, {0, 0, 1, 0}, {0, -1, 0, 0}, {-1, 0, 0, 0}};
γ2 = {{0, 0, 0, -I}, {0, 0, I, 0}, {0, I, 0, 0}, {-I, 0, 0, 0}};
γ3 = {{0, 0, 1, 0}, {0, 0, 0, -1}, {-1, 0, 0, 0}, {0, 1, 0, 0}};

Dop[Y_] :=
Hold[(I γ0.D[in, t] - I γ1.D[in, x] -
I γ2.D[in, y] - I γ3.D[in, z] + m id.in)] /. in :> Y // ReleaseHold

ψ = {u1, u2, u3, u4} Exp[ I (px x + py y + pz z) - I ℰ t];

secular =
Subtract @@@
Assuming[{px, py, pz, ℰ, t, x, y, z} ∈
Reals, FullSimplify[Table[Dop[ψ][[it]] == 0, {it, 1, 4}]]];

TableForm[Collect[secular, {u1, u2, u3, u4}]]


$$\begin{array}{c} \text{u1} (m+\mathcal{E})+\text{u4} (\text{px}-i \text{py})+\text{pz}\, \text{u3} \\ \text{u2} (m+\mathcal{E})+\text{u3} (\text{px}+i \text{py})-\text{pz}\, \text{u4} \\ \text{u3} (\mathcal{E}-m)+\text{u2} (\text{px}-i \text{py})+\text{pz}\,\text{u1} \\ \text{u4} (\mathcal{E}-m)+\text{u1} (\text{px}+i \text{py})-\text{pz}\, \text{u2} \\ \end{array}$$

This is the set of simultaneous quantities that have to be equal to zero in order for the plane wave to provide a solution. From this condition we can get the dispersion relation by setting the determinant of the coefficient matrix of the plane-wave amplitudes {u1, u2, u3, u4} to zero:

j = D[secular, {{u1, u2, u3, u4}}];

{ℰMinus, ℰPlus} = ℰ /. {ToRules[Reduce[Det[j] == 0, ℰ]]}


$\left\{-\sqrt {m^2+\text{px}^2+\text{py}^2+\text{pz}^2},\sqrt{m^2+\text{px}^2+\text{py}^2+\text{pz}^2}\right\}$

Here, Det[j] is the secular determinant (obtained by differentiating the vector secular with respect to the vector {u1, u2, u3, u4}).

The result gives the negative and positive energy eigenvalues of the operator in the original question, as a function of the plane-wave momentum components and the parameter m.