Particular Elliptic pde of divergence form with indicator boundary data and Feynman-Kac formula Consider divergence form elliptic pde in smooth boundary domain D
$$  Au:=\sum_{i,j}\partial_{i}(a_{ij}(x)\partial_{i}u(x)), $$
with boundary data $u|_{\partial D}:=1_{A}$ for $A\subset \partial D$. 
By Riesz representation there is a measure called the L-measure s.t.
$$u(x)=\int_{\partial D} \phi(y)d\omega_{x}^{L}(y)$$
if $u|_{\partial D}=\phi\in C(\partial D)$. For the harmonic case (Dirichlet problem) we have that the harmonic measure is harmonic and satisfies
$$\omega_{x}(A,\partial D)=P_{x}[B_{T_{\partial D}}\in A],$$
where B is Brownian motion. 

Q1: Do we have any such formula for nice enough sets A and pde Au that give a Feynman-Kac formula for some Ito process X:
  $$u(x)=\omega^{L}(x,A)=P_{x}[X_{T_{\partial D}}\in A].$$

In modern literature they mainly deal with continuous boundary data. For a related question in parabolic pdes see this question. 
The particular equation I have is in the upper half plane $\{(x,y):y>0\}$
$$\frac{1}{\beta}y \Delta u+\partial_{y}u=0$$
with $\beta>0$ and boundary data $u|_{R}=1_{R^{-}}$ the indicator on the half line. Indeed as mentioned below as well the divergence form I had in mind was $\nabla(y^{\beta}\nabla u)=0$ (which is the same in our case because $y>0$). 

Q2: If we can prove that we can apply Feynman-Kac then we obtain:
  $$u(x)=P[(X_{1,T_{\mathbb{H}}},X_{2,T_{\mathbb{H}}  })\in \mathbb{R}^{-}],$$
  where $dX_{1}=dB_{1}, dX_{2}=\frac{\beta}{X_{2}}dt+dB_{2}.$ 

Attempts
1)mainly going through the corresponding proofs for the harmonic measure. In Doob's potential theory book 2.IX.7/13 he proves the harmonic measure statement above. 
2)As with dirichlet we have the same representation for continuous boundary data:
$$E_{x}[\phi(X_{T_{\partial D}})]=u(x)=\int_{\partial D} \phi(y)d\omega_{x}^{L}(y).$$
So at least heuristically by approximating the indicator function $1_{R^{-}}$ by continuous functions we obtain something close to the desired formula above. 
 A: There is some ambiguity in your question. The operator $L$ (which you call $A$, but $A$ is also used for a set) is first given in a divergence form. Then it appears that you are interested in $L_1 = \beta^{-1} y \Delta + \partial_y$, which is not given in divergence form. Finally, in Q2, you consider the diffusion generated by $L_2 = \Delta + \beta y^{-1} \partial_y$. Both $L_1$ and $L_2$ have the same harmonic functions as the divergence form operator is $L_3 = \nabla \cdot (y^\beta \nabla)$.
All three operators clearly correspond to certain diffusion processes in the upper half-plane, which differ one from another by a random time change.
The process generated by $L_3$ has Brownian motion as the $x$ coordinate and Bessel diffusion as the $y$ coordinate (both coordinates are independent). If $\beta \ge 1$, the second coordinate will never hit $0$ and it will converge to $\infty$ as $t \to \infty$. Therefore, the Dirichlet problem in the half-space is ill-posed: any non-constant solution of $L_3 u = 0$ in the upper half-plane is unbounded.
If $\beta \in (-\infty, 1)$, then the problem is well-posed, and the harmonic measure for $L_3$ is indeed the distribution of the diffusion at its life-time. In fact, the solution can be given explicitly: $$u(x, y) = \int_{-\infty}^\infty u_0(x) \, \frac{c_\beta y^{1 - \beta}}{(x^2 + y^2)^{1 - \beta/2}} \, dx,$$ where $c_\beta$ is an appropriate normalisation constant. This can be verified directly, or by means of Fourier transform. For $\beta \in (-1, 1)$ this is precisely the Caffarelli–Silverstre extension technique for the fractional Laplacian, studied in the probabilistic context by Molchanov and Ostrovskii in 1969.
Time change does not essentially change the answer, except that the life-time may turn out to be infinite. If I am not mistaken, it will be finite for $L_2$ and $L_3$ and infinite for $L_1$.
The answer to your general question depends on your assumptions on the coefficients of $L$. If they are regular enough (uniformly elliptic is enough), then the answer is essentially "yes": there is a diffusion behind your PDE.
Final remark: strictly speaking, your probabilistic expression for the harmonic measure is not the Feynman–Kac formula (unless you take zero potential in F–K). In the classical case it is, I believe, due to Kakutani.
