fundamental solution of an elliptic PDE in divergence form with non-symmetric matrix

This question requires an articulated answer, since the topic dealt is complex and ramified. A fundamental solution for a not necessarily divergence form $2$nd order elliptic system with $C^{2,h}$ coefficients was first constructed by Georges Giraud in 1932 ([5]) by using the theory of multidimensional singular integrals he developed at the same time and independently of Solomon Mikhlin. Here, following Gaetano Fichera ([4], pp. 136-137 or pp. 672-674 of the English translation) we will offer a sketch of the construction of the local fundamental solution and then some hints on where to find the construction of the global one.

The problem
The systems of elliptic operators considered, assuming known the standard multiindex notation, is the following one $$ \sum_{0\le |\alpha| \le 2 } A_\alpha(x) D^\alpha u=\sum_{j=0}^2 A_j(x,D)u=A(x,D)u=f(x)\label{1}\tag{1} $$ where $\alpha\in\mathbb{N}^n$ is a multiindex, $A_\alpha(x)$, $0\le |\alpha|\leq 2$ and $A_j(x,D)$, $0\le j \le 2$ are $n\times n$ matrices, and $f,u$ are $n$ dimensional vector functions. The system \eqref{1} is assumed to be only elliptic in the sense that $$ \det\sum_{|\alpha| = 2} A_\alpha\xi^\alpha\neq 0\quad \forall \xi\in\mathbb{R}^n,\;\xi\neq 0\label{2}\tag{2} $$ We will also define the adjugate matrix differential operator $B_2$ as $$ B_2(x,D)=\mathrm{adj}\,A_2(x,D)= \mathrm{cof}\big[\,{^T\!A_2(x,D)}\big] $$

Construction of the parametrix for \eqref{1}
In this section and in the following one, we assume that the matrices in \eqref{1} and the non-homogeneous datum have a lower regularity, precisely in $C^{0,h}(E)$.
Let $L(x,D)$ be the scalar differential operator of second order whose characteristic polynomial is the first side of \eqref{2} and let $\psi(z,x-y)$ the parametrix of $L(x,D)$ constructed for $x=z$: then $$ S(z,x-y)=B_2(z,D_x)\psi(z,x-y)\label{3}\tag{3} $$ satisfies the following relationship $$ D_x^\alpha S(y,x-y)=O(|x-y|^{2-n-|\alpha|}\log|x-y|) $$ and therefore $$ A(x,D)_x S(y,x-y)=O(|x-y|^{h-n}\log|x-y|). $$ If $f\in [C^{0,h}(E)]^n$, then $$ u(x)=\int\limits_E S(y,x-y)f(y)\mathrm{d}y\in C^2(E) $$ and putting $K(x,y)=A(x,D)_x S(y,x-y)$ we have $$ A(x,D)u(x) =f(x)+\int\limits_E K(x-y)f(y)\mathrm{d}y $$ thus $S$ is a parametrix for the matrix differential operator $A(x,D)$: note that $S$ is globally defined on the whole $E$.

Construction of a local fundamental solution
Now fixing an arbitrary $x^o\in E$ and considering a ball $B(x^o,R)$ such that $$ \int\limits_{B(x^o,R)}\Vert K(x,y)\Vert \mathrm{d}y \le p <1 $$ we can construct the functional series $$ H(x,y)=\sum_{s=0}^\infty (-1)^s K_s(x,y)\:\text{ where }\: K_s(x,y)=\!\!\!\!\int\limits_{B(x^o,R)} K(x,t) K(t,y)_{s-1}\mathrm{d}t $$ which solves the following Fredholm integral equation of the second kind $$ H(x,y)+K(x,y)+\!\!\!\!\!\int\limits_{B(x^o,R)} K(x,t)H(t,y)\mathrm{d}t=0 $$ This implies that the function $$ F(x,y)=S(y,x-y)+\!\!\!\!\!\int\limits_{B(x^o,R)} S(t,x-t)H(t,y)\mathrm{d}y \label{4}\tag{4} $$ is the sought for local fundamental solution for the system \eqref{1}.

Construction of the global fundamental solution
In this section we shift to a still more colloquial style of exposition: the only technical detail we remark is that for the global construction we must assume that the coefficients of $A(x,D)$ belong to $C^{2,h}(E)$ since we should be able to define its transposed ${^T\!A}(x,D)$.
Basically there are two ways of constructing a global fundamental solution in the spirit of the development sketched above. The older one was developed by Giraud [5], by modifying the approach leading to \eqref{4}: by suitably defining a $K(x,y)$ with a nice behavior at infinity, he is able to define again a Fredholm integral equation, whose solution defined on all $E$. A survey of his work can be found in the book by Miranda ([8] §20 pp. 67-73) which, however address to the original work of Giraud.
The more recent method is due to Fichera (see [2] and [3]): he simplifies the approach by considering only bounded domains $E$ enclosed in a suitable ball $B(x,R)$ and by imposing suitable boundary conditions on $\partial B(x,R)$. Despite not defining the fundamental solutions for unbounded $E$, Fichera's approach was developed by him in order to work also for higher order elliptic systems.

Final notes
The case of a $2$nd order divergence form elliptic system with a symmetric matrix of coefficients of very low regularity (i.e. bounded measurable) was studied by Littman, Weinberger and Stampacchia in [7] (see also [9]). These results were later generalized to the case of a $2$nd order uniformly elliptic system with a non necessarily symmetric matrix of bounded measurable coefficients by Grüter and Widman (see [6]) and to higher order systems by Ariel Barton ([1]) who nicely surveys also the recent progresses of the studies related to the construction of Green's functions and fundamental solutions of elliptic systems with measurable coefficients. It seems that all of the works cited in this last section use the divergence structure (and thus implicitly the variational structure) of the elliptic system analyzed: I do not know if there are studies (for example using the very weak solution concept) where the fundamental solution/Green function is constructed in analogy with the work of Giraud and Fichera, i.e. without assuming any particular structure, and with bonded measurable coefficients.

Bibliography

[1] Ariel Barton (2016), "Gradient estimates and the fundamental solution for higher-order elliptic systems with rough coefficients", Manuscripta Mathematica, 151, 3-4, pp. 375-418, MR3556825, Zbl 1358.35035.

[2] Gaetano Fichera (1961), "Linear elliptic equations of higher order in two independent variables and singular integral equations, with applications to anisotropic inhomogeneous elasticity", in Langer, Rudolph E., PARTIAL DIFFERENTIAL EQUATIONS AND CONTINUUM MECHANICS, Proceedings of an International Conference Conducted by the Mathematics Research Center at the University of Wisconsin, Madison, June 7-15, 1960, Publication of the Mathematics Research Center, United States Army, the University of Wisconsin., no. 5, Madison: The University of Wisconsin Press, pp. 55–80, MR0156084, Zbl 0111.29602.

[3] Gaetano Fichera (1962), "La soluzione fondamentale principale per una equazione differenziale ellittica di ordine superiore". Bulletin Mathématique de La Société des Sciences Mathématiques et Physiques de la République Populaire Roumaine, 6 (54)(3/4), nouvelle série, pp. 139-149. Retrieved from JSTOR. MR0185267, Zbl 0196.40702.

[4] Gaetano Fichera (1964), "Problemi elastostatici con vincoli unilaterali: il problema di Signorini con ambigue condizioni al contorno", (Italian), Atti della Accademia Nazionale dei Lincei. Memorie. Classe di Scienze Fisiche, Matematiche e Naturali, Serie VIII, 7 (2): 91–140, Zbl 0146.21204. Translated in English as Gaetano Fichera(1964), "Elastostatic problems with unilateral constraints: the Signorini problem with ambiguous boundary conditions", Seminari dell'istituto Nazionale di Alta Matematica 1962–1963, Rome: Edizioni Cremonese, pp. 613–679 (this last version is the one available also in his "Opere scelte").

[5] Georges Giraud (1932, "Généralisation des problèmes sur les opérations du type elliptique", (French) Bulletin des Sciences Mathématiques, II Série, vol. 56, part 1, pp. 248-272, pp. 281-312, pp. 316-352 and 384, JFM 58.0494.02

[6] Michael Grüter and Kjell-Ove Widman (1982), "The Green function for uniformly elliptic equations", Manuscripta Mathematica 37, pp. 303-342 Zbl 0485.35031.

[7] Walter Littman, Hans Weinberger and Guido Stampacchia (1962), "Regular points for elliptic equations with discontinuous coefficients", Annali della Scuola Normale Superiore di Pisa, Classe di Scienze, serie III, Vol. 17, n° 1-2, pp. 43-77, MR161019, Zbl0116.30302.

[8] Carlo Miranda (1970) [1955], Partial Differential Equations of Elliptic Type, Ergebnisse der Mathematik und ihrer Grenzgebiete – 2 Folge, Band 2, translated by Motteler, Zane C. (2nd Revised ed.), Berlin – Heidelberg – New York: Springer Verlag, pp. XII+370, doi:10.1007/978-3-642-87773-5, ISBN 978-3-540-04804-6, MR 0284700, Zbl 0198.14101

[9] Guido Stampacchia (1966), "Équations elliptiques du second ordre à coefficients discontinus" (notes du cours donné à la 4me session du Séminaire de mathématiques supérieures de l'Université de Montréal, tenue l'été 1965), (in French), Séminaire de mathématiques supérieures 16, Montréal: Les Presses de l'Université de Montréal, pp. 326, ISBN 0-8405-0052-1, MR0251373, Zbl 0151.15501.

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.

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.