Fundamental solution of an elliptic PDE in divergence form with non-symmetric matrix I am looking for the fundamental solution of the following PDE
$$\partial_i (a^{ij}\partial_j u)=f$$
where $a^{ij}(x)$ is a non-symmetric matrix with possibly non-constant coefficients.
I could find a paper Clements - A fundamental solution for linear second-order elliptic systems with variable coefficients for the symmetric case, I am wondering does anyone know a paper for non-symmetric case? Or it is much harder than the symmetric case to find its fundamental solution? I searched a lot but I am asking because I thought maybe I have missed something.
Any help would be very much appreciated.
 A: If you assume that the coefficients $a^{ij}$ are smooth functions and let $$b^{ij} = \frac{1}{2}(a^{ij} + a^{ji}),$$ then the PDE can be written as
$$
b^{ij}\partial^2_{ij}u + \partial_ia^{ij}\partial_ju = f.
$$
If the symmetric matrix $[b^{ij}]$ is positive definite, then the PDE is elliptic. In that case, there exists a fundamental solution on $\mathbf{R}^n$. One (but probably not the best way) to prove this is to construct a right inverse as a pseudodifferential operator. This operator has a kernel, which therefore is the fundamental solution.
ADDED: You can find the details of how to construct a right inverse in a book that presents pseudodifferential operators and how to use them to prove regularity and existence of solutions to elliptic PDEs. Roughly, it goes like this: First, you use pseudodifferential calculus (too long to explain here but Bazin explains the key ideas) to solve for, on an open domain, a parametrix, which is a pseudodifferential operator of order $-2$ such that
$$
PQ = I + R,
$$
where $P$ is your differential operator (of order $2$), $I$ is the identity operator, and $R$ is a pseudodifferential operator of order $-k$, for some sufficiently large $k$ (usually called a smoothing operator). The exact symbols for $Q$ and $R$ can in principle be computed. If $\Omega$ is chosen small enough, then the operator norm of $R$ in the appropriate Sobolev space will be less than $1$. Therefore, $I+R$ can be inverted and $Q(I+R)^{-1}$ is the right inverse. By the Schwarz kernel theorem, there is a fundamental solution corresponding to this. The catch here is that the only explicit formulas for $(I+R)^{-1}$ involve infinite sums, including the obvious one using a geometric series.
A: Let me start with a constant coefficient operator
$$
P(D)=\sum_{1\le j, k\le n} a_{jk} D_jD_k,\quad D_j=\frac{\partial }{i\partial x_j}.
$$
Note that in two dimensions, you have elliptic operators with are squares of linear forms, such as $(D_1+iD_2)^2$. To avoid this particular case, we shall
assume that $P$ is elliptic and $n\ge 3$, then defining $A=(a_{jk} )$, you may assume that $\Re A=(A+A^*)/2$ is positive-definite: in fact in dimension $\ge 3$ the range of $A$ is a cone with aperture $<π$ in the complex plane, so that after multiplication by a complex number, you may assume that this cone is symmetric with respect to the real line in $\mathbb C$. We have thus
$$
P(D)=\langle A_1 D, D\rangle+i\langle A_2 D, D\rangle,\quad A_1=(A+A^*)/2 \quad A_2=(A-A^*)/(2i).
$$
As a result, $P(D)$ is the Fourier multiplier
$$
\langle A_1 \xi, \xi\rangle+i\langle A_2 \xi, \xi\rangle,\quad
\langle A_1 \xi, \xi\rangle\ge c_0\vert \xi\vert^2,\quad c_0
>0.$$
A parametrix for the operator $P(D)$ is thus the Fourier multiplier
$$
E(\xi)=\bigl(\langle A_1 \xi, \xi\rangle+i\langle A_2 \xi, \xi\rangle\bigr)^{-1}
=\langle A \xi, \xi\rangle^{-1}.
$$
There is no reason to be worried by the singularity of $E$, since in the first place $E$ is vanishing only at $\xi=0$ so that the operator $E(D)$ defined by 
$$
(E(D) u)(x)=\int e^{ix\cdot \xi} \langle A \xi, \xi\rangle^{-1}\hat u(\xi) d\xi(2π)^{-n},
$$
is well-defined for $u$ in the Schwartz space and even for a tempered distribution such that $\hat u$ is locally bounded since the singularity of $E(\xi)$ is of type $\vert \xi\vert^{-2}$ which is locally integrable since $n>2$ (you will need as well some decay at infinity for $\hat u$, but much less than fast decay). There is even a nicer version, using the homogeneity of $E(\xi)$ (homogeneous of degree $-2$), so that its Fourier transform is homogeneous with degree $2-n$: The operator $E(D)$ appears as a convolution operator with $\hat E(-x)$.
Of course you can generalize this for an elliptic operator with $C^\infty$ coefficients, using pseudo-differential operators: then the principal symbol in the asymptotic expansion will be
$
 \langle A(x) \xi, \xi\rangle^{-1}.
$
If you stick with general constant coefficient operator, you can define a parametrix using essentially the same idea of inverting $ \langle A \xi, \xi\rangle$, but then you will have to use a theorem of division of distributions by a polynomial (something that can always be performed, you may even divide by an analytic function).
On the other hand, if you take a general operator with variable coefficients and  complex-valued symbol, you may run into trouble, because some of them are not even locally solvable, so that no decent parametrix could be defined.
