I recently studied similar questions. In what follows, I try to give a summery of the information I gathered. I hope this is somewhat helpful.
2.) I try to answer your second question first. I also find it hard to find Information about this topic and the reason for this is in my opinion, that there are no deep theorems to prove about Duhamel's Principle. To the best of my knowledge, Duhamel's Principle is more an idea and it has to be applied to every case specifically using other already well-known tools. I hope the following makes it more clear.
1.) The thing to keep in mind concerning your first question is that Duhamel's Principle is used to construct solutions for non-homogeneous linear evolution equations. In your case the term $f(u)$ makes the equation in general non-linear, thus, you can't use Duhamel's principle to construct a solution. However, assume you already have a solution $u$ satisfying your equation then $u$ satisfies the representation formula that you describe in 3. In the following I will derive this representation formula. First, we have to know what is meant by $u(\cdot,0) = u_0(\cdot)$ in the case where $u_0$ is a measure. $u_0$ is a measure, so we can't understand this equation in a point-wise sense. This is why one usually interprets this equation in a distribution sense, i.e., we require that for all $\phi \in C_0^\infty(\mathbb{R})$, i.e. for all $\phi$ that are smooth and vanish at infinity, the following $$\lim_{t\rightarrow 0} \int u(x,t) \phi(x) \, \mathrm{d}x = \int \phi(x) \, \mathrm{d}u_0(x). $$
Note that in particular this is true for $\phi(x) = K(x,s)$ for any $s>0$. Now assume $u$ is a $C^{2,1} \cap L^\infty$ solution for your problem and it satisfies the boundary condition in a distributional sense. First, we look at the following:
$$\mathrm{d}/\mathrm{d}s \int K(x-y,t-s) u(y,s) \, \mathrm{d}y = \int K(x-y,t-s)(u_t- u'')(y,s) \, \mathrm{d}y,$$
where we interchanged the limit and derivative (by dominated convergence), applied integration by parts, and used that $K_t - K'' = 0$ for $t>0$. Next, we multiply $u_t -u'' = f(u)$ with $K(x,t-s)$ and we integrate. We get the following
\begin{align}
&\int_\epsilon^{t-\epsilon} \int_{\mathbb{R}} K(x-y,t-s) (u_t(y,s)-u''(y,s)) \, \mathrm{d}y\mathrm{d}s \\ = &\int_\epsilon^{t-\epsilon}\int_{\mathbb{R}} K(x-y,t-s) f(u(y,s)) \, \mathrm{d}y\mathrm{d}s \\ = &\int_{\mathbb{R}} K(x-y,\epsilon)u(y,t -\epsilon) \, \mathrm{d}y - \int_{\mathbb{R}} K(x-y,t-\epsilon)u(y,\epsilon) \, \mathrm{d}y,
\end{align}
where we used the equation derived above and the fundamental theorem of calculus. Letting $\epsilon \rightarrow 0$ in the last equation, yields $$ u(x,t) = \int_\mathbb{R} K(x-y,t) \, \mathrm{d}u_0(y) + \int_0^t \int_\mathbb{R} K(x-y,t-s) f(u(y,s)) \, \mathrm{d}y \mathrm{d}s,$$ where we used the initial value condition in a distributional sense and that $\lim_{t\rightarrow 0} \int K(x-y,t) g(y) \, \mathrm{d}y = g(x)$ for a continuous bounded function $g$.
3.) Let's illustrate this point with an example. Assume you have a non-negative measurable function $U$ satisfying the integral equation above, i.e.,
$$ U(x,t) = \int_\mathbb{R} K(x-y,t) \, \mathrm{d}u_0(y) + \int_0^t \int_\mathbb{R} K(x-y,t-s) f(U(y,s)) \, \mathrm{d}y \mathrm{d}s,$$
but we don't suppose any regularity on $U$ in the sense of differentiability. If you now can additionally show that $U \in L^\infty_{loc}(0,T;L^\infty(\mathbb{R}))$, then you can for instance deduce that $U$ is hölder continuous and therefore, using parabolic Schauder estimates you get that $U\in C^{2,1}$ is actually a regular solution. So basically if you can find a positive measurable function $U$ that satisfies the integral equation and you can show that the integral converges, then you can bootstrap to get more regularity, thus a classical solution.
Best Answer
Let $u$ and $v$ be two solutions with the stated properties. Let $\eta$ be a standard mollifier on $\mathbb{R}^d$ and let $u_\varepsilon :=u \ast \eta_\epsilon$ and $v_\varepsilon := v \ast \eta_\epsilon$, where we only convolve with respect to the $x$-variable.
Claim: We have $u_\varepsilon, v_\varepsilon \in C ^\infty((0, \infty) \times \mathbb{R} ^d) \cap C ( [0, \infty) \times \mathbb{R} ^d)$, $\partial_t u_\varepsilon - \Delta u_\varepsilon = \partial_t v_\varepsilon - \Delta v_\varepsilon = 0$ and $u_\varepsilon(0) = v_\varepsilon(0) = u_0 \ast \eta_\epsilon \in C_b ( \mathbb{R} ^d)$.
Proof of the Claim: As $u$ and $v$ are bounded, we can pass the differentiation into the integrals and get \begin{equation*} \partial_t u_\varepsilon - \Delta u_\varepsilon = (u_t - \Delta u) \ast \eta_\epsilon = 0 \end{equation*} and similarly $\partial_t v_\varepsilon - \Delta v_\varepsilon = 0$. Furthermore, $u_\varepsilon, v_\varepsilon \in C ^\infty((0, \infty) \times \mathbb{R} ^d)$ and $u_0 \ast \eta_\epsilon \in C_b ( \mathbb{R} ^d)$. Thus, we are left to show continuity in $t = 0$. This follows by the weak$^*$ convergence: Let $x \in \mathbb{R} ^d$. Then \begin{align*} u_\varepsilon (t,x) - u_\varepsilon (0,x) & = \int_{ \mathbb{R} ^d }\! ( u(t,y) - u_0 (y) ) \eta_\epsilon (x - y) \, dy \\ & \to 0 \end{align*} as $\eta_\epsilon (x - \cdot) \in L^1(\mathbb{R} ^d)$ and $u(t,\cdot) \to u_0$ weak$^*$ in $L^\infty ( \mathbb{R}^d )$ as $t \to 0$. In the same fashion, we see that $v_\varepsilon$ is continuous in $t = 0$, which proofs the claim.
The heat equation for regular initial data in $C_b(\mathbb{R}^d)$ has a unique bounded classical solution, see Thm. 6 in Evans, Chapter 2.3. Hence, $u_\varepsilon = v_\varepsilon$ for all $\varepsilon > 0$. As $u_\varepsilon \to u$ and $v_\varepsilon \to v$ locally uniformly in $(0,\infty) \times \mathbb{R} ^d$ as $\varepsilon \downarrow 0$, we get $u = v$ in $(0,\infty) \times \mathbb{R} ^d$.