Let $\Phi_t(x)$ be the flow of your vector field $\nabla d_{\mathcal{A}}$, starting at $x \in \mathcal{A}$. It will exist until the time $t_0(x)$, when it hits the set where $\nabla d_{\mathcal{A}}$ is not smooth.
Consider the ODE
$$2 \frac{\partial}{\partial t} z(t, x) + z(t, x)\cdot \Delta d_{\mathcal{A}}(\Phi_t(x)) = 0, ~~~~~~~ z(0, x) = W(x)$$
depending on the parameter $x \in \partial \mathcal{A}$. It is a linear ODE, hence it has a unique solution that exists for all time, or rather until $t_0(x)$ (because afterwards the equation does not make sense anymore). It is also well-known that the solution depends smoothly on the parameter $x$, as all the data is smooth.
Now let $\mathcal{B} = \{\Phi_t(x) \mid x \in \partial \mathcal{A}, t < t_0(x)\} \subseteq \mathcal{A}$ be the set of points $y$ that lie on a unique flowline, i.e. there exist a unique $t$ and a unique $x \in \partial \mathcal{A}$ such that $y = \Phi_t(x)$ (that this set has the claimed description is because flowlines cannot cross). This is a neighborhood of $\partial \mathcal{A}$ in $\mathcal{A}$, and its complement in $\mathcal{A}$ has measure zero.
For points in $\mathcal{B}$, set
$$Z(\Phi_t(x)) := z(t, x).$$
This gives you a solution that has the desired boundary values and is smooth on $\mathcal{B}$.
However, you know nothing about the smoothness in $\mathcal{A} \setminus \mathcal{B}$. For example, if you take a disc in $\mathbb{R}^n$, your function $Z$ will be continuous if and only if $W(x) \equiv \mathrm{const}$. Otherwise it will not be continuous at the origin.
It is true that for any initial datum $u_0\in C^\infty(\mathbb{R})$ there exists a solution $u\in C^\infty(\mathbb{R}^+\times\mathbb{R})$ to the heat equation with initial condition $u(0,x)=u_0(x)$. As you point out, this will not be unique.
I can give a method of constructing such solutions now. The idea is to show that we can write $u=\sum_{n=1}^\infty f_n$ where $f_n$ are carefully constructed solutions chosen such that the partial sums $\sum_{n=0}^mf_n(0,x)$ eventually agree with $u_0(x)$ any bounded subset of the reals, and for which $f_n$ tends to zero arbitrarily quickly in the compact-open topology. First a bit of notation. I use $\mathbb{R}^+=[0,\infty)$ for the nonnegative reals. For a space $X$ then $C_0(X)$, $C^\infty(X)$, $C^\infty_0(X)$, and $C^\infty_K(X)$, represent the continuous real-valued functions on $X$ which are respectively vanishing at infinity, smooth, smooth and vanishing at infinity, and smooth with compact support. Let $(P_t)_{t\geq 0}$ be the kernels
$$
\begin{align}
&P_t\colon C_0(\mathbb{R})\to C_0(\mathbb{R}),\\\\
&P_tu(x)=\frac{1}{\sqrt{4\pi t}}\int_\mathbb{R}e^{-(x-y)^2/4t}u(y)\,dy
\end{align}
$$
for $t > 0$, and $P_0u=u$. This is the Markov transition function for Brownian motion (more precisely, for standard Brownian motion scaled by $\sqrt{2}$, because of the normalization used here). For $u\in C_K^\infty(\mathbb{R})$, then $f(t,x)=P_tu(x)$ is in $C_0^\infty(\mathbb{R}^+\times\mathbb{R})$, and is a solution to the heat equation with initial condition $f(0,x)=u(x)$, agreeing with the classical solution stated in the question. I'll also consider initial conditions $u\in C^\infty_K((a,\infty))$ (for $a\in\mathbb{R}$) by setting $u(x)\equiv0$ for all $x\le a$. The first step in the construction is to find initial conditions supported in $(a,\infty)$ so that $P_tu(0)$ approximates any given continuous function of time that we like.
1) For any $a > 0$ and $h\in C_0((0,\infty])$, there exists a sequence $u_1,u_2,\ldots\in C^\infty_K((a,\infty))$ such that $\sqrt{t}P_tu_n(0)$ converges uniformly to $h(t)$ (over $t > 0$) as $n\to\infty$.
Consider the closure, $V$, in $C_0((0,\infty])$ (under the uniform norm) of the space of functions $t\mapsto\sqrt{4\pi t}P_tu(0)$ for $u\in C^\infty_K((a,\infty))$. Note that the limit as $t\to\infty$ does always exist and is just the integral of $u$. Then $V$ is a closed linear subspace of $C_0((0,\infty])$. Consider a sequence $u_n\in C_0((a,\infty))$ tending to the delta function $\delta_b$ at a point $b > a$, in the sense that $u_n$ all have support in the same compact set, and converge in distribution to $\delta_b$. Then, from the expression defining $P_t$, $\sqrt{4\pi t}P_tu_n(0)$ converges uniformly over $t > 0$ to $\exp(-b^2/4t)$. So, the function $t\mapsto\exp(-b^2/4t)$ is in $V$. As the set of functions of the form $t\mapsto\exp(-b^2/4t)$ for $b > a$ is closed under multiplication and separates points, the locally compact version of the Stone-Weierstrass theorem says that $V=C_0((0,\infty])$. Then, (1) follows.
2) For any $u\in C^\infty_K((0,\infty))$ and $a,T > 0$, there exists a sequence $u_1,u_2,\ldots\in C^\infty_0((a,\infty))$ such that $f_n(t,x)\equiv P_t(u+u_n)(x)$ converges uniformly to zero (along with all its partial derivatives to all orders) over $t\in[0,T]$ and $x\le0$.
Choosing $0 < \epsilon < a$ so that the support of $u$ is contained in $(\epsilon,\infty)$, (1) implies that we can choose $u_n\in C^\infty_K((a,\infty)$ so that $\sqrt{t}P_tu_n(\epsilon)$ converges uniformly to $-\sqrt{t}P_tu(\epsilon)$ over $t\ge0$ as $n\to\infty$. Then, $f_n(t,x)\equiv P_t(u+u_n)(x)$ is a bounded solution to the heat equation with boundary conditions $f_n(0,x)=0$ for $x\le\epsilon$ and $f_n(t,\epsilon)=P_t(u+u_n)(\epsilon)$. It is then standard that the solution is given by an integral over the boundary,
$$
f_n(t,x)=\int_0^t\frac{\epsilon-x}{\sqrt{4\pi (t-s)^3}}e^{-(\epsilon-x)^2/4(t-s)}\sqrt{s}P_s(u_n+u)(\epsilon)\frac{ds}{\sqrt{s}}
$$
for $x\le0$. Differentiating this wrt $x$ and $t$ an arbitrary number of times, and using dominated convergence as $n\to\infty$, it follows that $f_n(t,x)$ (and all its partial derivatives) converge uniformly to zero over $x\le0$ and $t\in[0,T]$ as $n\to\infty$.
3) Suppose that $T > 0$, $0 < a < b$ and $u\in C^\infty_K(\mathbb{R})$ has support contained in $(-\infty,-a)\cup(a,\infty)$. Then, there exists a sequence $u_n\in C^\infty_K(\mathbb{R})$ with supports in $(-\infty,-b)\cup(b,\infty)$ such that $P_t(u+u_n)(x)$ (and all its partial derivatives) tends to 0 uniformly over $\vert x\vert\le a$ and $t\in[0,T]$.
Set $u^+(x)=1_{\{x > 0\}}u(x)$ and $u^-(x)=1_{\{x < 0\}}u(x)$. Applying (2) to $u^+(x+a)$, there exists a sequence $u^+_n\in C^\infty_K(\mathbb{R})$ with supports in $(b,\infty)$ such that $P_t(u^++u^+_n)(x)$ tends to 0 uniformly over $x\le a$ and $t\in[0,T]$. Applying the same argument to $u^-(-x-a)$, there exists $u^-_n\in C^\infty_K(\mathbb{R})$ with support in $(-\infty,-b)$ such that $P_t(u^-+u^-_n)(x)$ tends to zero uniformly over $x\ge -a$ and $t\in[0,T]$. A sequence satisfying the statement of (3) is $u_n=u^+_n+u^-_n$.
4) If $u\in C^\infty(\mathbb{R})$ then there exists $f\in C^\infty(\mathbb{R}^+\times\mathbb{R})$ solving the heat equation with initial condition $f(0,x)=u(x)$.
We can inductively choose a sequence $u_n\in C^\infty_K(\mathbb{R})$ such that $\sum_{m=1}^nu_m(x)=u(x)$ for $\vert x\vert\le n+1$ and $n\ge1$. Let $u_1=u$ on $[-2,2]$ and then, for each $n\ge2$, apply the following step.
- As $\tilde u=u-\sum_{m=1}^{n-1}u_m$ is zero on $[-n,n]$, it has support in $(-\infty,-a)\cup(a,\infty)$ for $a=n-1/2$. Choosing $b=n+1$, by (3), there exists $v\in C^\infty_K(\mathbb{R})$ with support in $(-\infty,-b)\cup(b,\infty)$ such that $P_t(\tilde u +v)(x)$ along with all its partial derivatives up to order $n$ are bounded by $2^{-n}$ over $\vert x\vert\le n-1/2$ and $t\le n$. Take $u_n=\tilde u + v$.
Setting $f_n(t,x)=P_tu_n(x)$, then $f_n$ are smooth functions satisfying the heat equation and the initial conditions $\sum_{m=1}^nf_m(0,x)=u(x)$ for $\vert x\vert\le n+1$. Also, by the choice of $u_n$, for $n > 1$ then $f_n$ together with all its derivatives up to order $n$ is bounded by $2^{-n}$ on $[0,n]\times[1-n,n-1]$. As $\sum_n2^{-n} < \infty$, the limit $f=\sum_nf_n$ exists and is smooth with all partial derivatives commuting with the summation. Then $f$ has the required properties.
Best Answer
I've played with such things 15 years ago. Here's what I remember...
Sure you need some semi-boundedness condition on K. Uniqueness should then follow via semigroup domination (aka Kato's inequality) from scalar case.
I. Shigekawa, L^p Contraction Semigroups for Vector Valued Functions, J. Funct. Analysis 147 (1997), 69-108
Shigekawa had rediscovered Barry Simon's semigroup domination criterion (i.e. Kato's inequality) and generalized it to the vector valued things.
See also:
R.S. Strichartz, Analysis of the Laplacian on the Complete Riemannian Manifold, J. Funct. Analysis 52 (1983), 48-79. Scalar case: Theorem 3.5 - but has uniqueness only for 1< p<∞.
(Some details in Strichartz' proof(s) for vector valued Laplacians not optimal (classic problems with doing tensor calculus and Stokes...) not even in the follow-up paper "L^p contractive projections and the heat semigroup for differential forms" J. Funct. Analysis, 65 (1985), 348-357)
Perhaps look also for El Maati Ouhabaz ca. 1999. He's written a book, "Analysis of Heat Equations on Domains" (2004) but I never got hold of it.