There's a relation between twodimensional Brownian motion and conformal maps, see e.g. Thurston's answer to this question. Given two nonempty simplyconnected domains $U$ and $V$ in the complex plane which are not equal to the complex plane, the Riemann mapping theorem guarantees the existence of a conformal bijection $U\to V$. Is there a way to use Brownian motion to construct such a map?
Let $U$ be a Jordan domain and $b_0,b_1,b_2$ be points appearing in order on $\partial U$. Let $f : U \to \mathbb{H}$ be the Riemann map normalized so that $f(b_0)=0$, $f(b_1)=1$ and $f(b_2)=\infty$. Chris Bishop taught me a poor man's way to construct $f$ using Brownian motion.
Let $A_k$ be the arc of $\partial U$ between $b_k$ and $b_{k+1}$, where the index is taken modulo 3. Similarly, let $B_0$ denote the interval $(0,1)$, $B_1$ the interval $(1,+\infty)$ and $B_2$ the interval $(\infty,0)$, so that $f$ maps $A_k$ to $B_k$ for each $k$. For $z\in U$, let $p_k$ denote the probability that Brownian motion started at $z$ first hits $\partial U$ along $A_k$. The probability $p_k=\omega(z,U)(A_k)$ is known as the harmonic measure of $A_k$ with respect to $(z,U)$ and is invariant under conformal maps. One can show there is a unique point $w \in \mathbb{H}$ such that $\omega(w,\mathbb{H})(B_k)=p_k$ for each $k$. By conformal invariance, we have $w=f(z)$.
Here is how to find $w$ given $(p_0,p_1,p_2)$. If $C_0,C_1,C_2$ are consecutive arcs on the unit circle such that $C_k$ measures $2\pi p_k$ radians, then the harmonic measure $\omega(0,\mathbb{D})(C_k)$ is obviously $p_k$. Now, let $c_k$ be the point at the junction between $C_{k1}$ and $C_k$ and let $M$ be the Moebius transformation sending $c_0,c_1,c_2$ to $0,1,\infty$ respectively. You can write $M$ explicitly as $$ M(\zeta) = \frac{(c_1c_2)(\zetac_0)}{(c_1c_0)(\zetac_2)}. $$ It is easy to see that $M$ maps the unit disk to the upper halfplane. By conformal invariance, we have $$\omega(M(0),\mathbb{H})(B_k)=\omega(0,\mathbb{D})(C_k)=p_k$$ for each $k$ and hence $M(0)=w=f(z)$.
I implemented this in Matlab some years ago. What you do is for each $z$ where you want to evaluate the Riemann map, throw a bunch of Brownian motions to estimate the probabilities $p_0$,$p_1$,$p_2$. It was not very efficient, but was quite easy to program.
I believe that a way of doing this is discussed in Chapter 5 of Richard Bass' "Probabilistic Techniques in Analysis." If I remember correctly, the main place Brownian motion is used is to construct the Green function for the domain.
This chapter of this book discusses how the proofs of many theorems in complex analysis can be proven with Brownian motion... Picards theorem and Carleson's corona theorem in particular. Most of the proofs stem from Thurston's observation that Brownian motion is conformally invariant, but it is a good read for those interested.

$\begingroup$ Is there any discussion in the literature about the relation between BM and Conformal maps in higher dimensions? Or this is strictly planar phenomenon? $\endgroup$ Nov 3 '19 at 11:07
I think the problem in higher dimensions is not BM but rather conformal maps, or lack thereof. You can see the details in this post.
Conformal maps in higher dimensions
I don't think BM has anything much to do with these maps, and it seems to me that BM wouldn't even be invariant under inversions in dimension 3 and above, since it is transient there, so $B_t$ eventually wanders off to infinity, but then the inversion of it would wander off to the origin, which makes no sense.