[proofplan]
The $l = 0$ case follows from connectedness: a smooth function with vanishing differential is locally constant, and path-connectedness forces it to be globally constant. For $l \ge 1$, we construct an explicit chain homotopy $\iota : \Omega^l(M) \to \Omega^{l-1}(M)$ satisfying the homotopy identity $d \circ \iota + \iota \circ d = \mathrm{id}_{\Omega^l(M)}$. This identity immediately implies that every closed $l$-form is exact: if $d\omega = 0$ then $\omega = d(\iota\omega)$. The formula for $\iota$ integrates a $k$-form along the radial ray from the origin to each point, exploiting the star-shaped geometry of $M$.
[/proofplan]
[step:Compute $H^0_{\mathrm{dR}}(M)$ from connectedness]
A $0$-form is a smooth function $f \in C^\infty(M)$. The coboundary condition is $df = 0$, which means
\begin{align*}
\frac{\partial f}{\partial x_i}(x) = 0 \quad \text{for all } x \in M,\ i = 1, \dots, k.
\end{align*}
Since $M$ is an open subset of $\mathbb{R}^k$ it is locally path-connected, and since $M$ is convex (star-shaped implies connected) it is path-connected. A function on a path-connected open set whose total derivative vanishes identically is constant: for any two points $p, q \in M$, join them by a $C^1$ path $\gamma : [0,1] \to M$ and compute
\begin{align*}
f(q) - f(p) = \int_0^1 \frac{d}{dt}(f \circ \gamma)(t) \, d\mathcal{L}^1(t) = \int_0^1 df_{\gamma(t)}(\dot\gamma(t)) \, d\mathcal{L}^1(t) = 0.
\end{align*}
Thus $\ker(d : \Omega^0(M) \to \Omega^1(M)) = \{\text{constant functions}\} \cong \mathbb{R}$. Since there are no $(-1)$-forms, $H^0_{\mathrm{dR}}(M) = \ker d / \{0\} \cong \mathbb{R}$.
[guided]
The de Rham cohomology in degree zero is $H^0_{\mathrm{dR}}(M) = \ker(d : \Omega^0 \to \Omega^1) / \mathrm{im}(d : \Omega^{-1} \to \Omega^0)$. Since there are no forms of degree $-1$, the denominator is $\{0\}$ and we simply need to identify the kernel of $d$ on smooth functions.
A smooth function $f : M \to \mathbb{R}$ satisfies $df = 0$ if and only if $\partial_{x_i} f(x) = 0$ for all $x \in M$ and all $i = 1, \dots, k$. Why does this force $f$ to be constant? Because $M$ is path-connected (any star-shaped set is convex, and convex sets are path-connected), and for any $C^1$ path $\gamma : [0,1] \to M$ from $p$ to $q$, the fundamental theorem of calculus gives
\begin{align*}
f(q) - f(p) = \int_0^1 \frac{d}{dt}(f \circ \gamma)(t) \, d\mathcal{L}^1(t) = \int_0^1 df_{\gamma(t)}(\dot\gamma(t)) \, d\mathcal{L}^1(t) = 0,
\end{align*}
where the last equality uses $df \equiv 0$. So $f(p) = f(q)$ for all $p, q \in M$, i.e., $f$ is constant.
The map $f \mapsto f(x_0)$ (for any fixed $x_0 \in M$) defines a linear isomorphism from $\ker d$ to $\mathbb{R}$. Hence $H^0_{\mathrm{dR}}(M) \cong \mathbb{R}$.
[/guided]
[/step]
[step:Define the contracting homotopy $\iota : \Omega^l(M) \to \Omega^{l-1}(M)$]
Fix a basepoint $0 \in M$ (which lies in $M$ since $M$ is star-shaped with respect to $0$). For $l \ge 1$ and $\omega \in \Omega^l(M)$, define the operator
\begin{align*}
\iota : \Omega^l(M) &\to \Omega^{l-1}(M)
\end{align*}
as follows. Write $\omega$ in the global coordinates $(x_1, \dots, x_k)$ on $M \subset \mathbb{R}^k$ as
\begin{align*}
\omega = \sum_{I} \omega_I \, dx_I,
\end{align*}
where the sum is over strictly increasing multi-indices $I = (i_1 < \cdots < i_l)$ and $dx_I = dx_{i_1} \wedge \cdots \wedge dx_{i_l}$. For each such $I$ and each $j = 1, \dots, l$, define
\begin{align*}
(\iota\omega)_x := \sum_{I} \sum_{j=1}^{l} (-1)^{j-1} \left( \int_0^1 t^{l-1} \omega_I(tx) \, d\mathcal{L}^1(t) \right) x_{i_j} \, dx_{I \setminus \{i_j\}},
\end{align*}
where $dx_{I \setminus \{i_j\}}$ denotes the $(l-1)$-form obtained by deleting $dx_{i_j}$ from $dx_I$. This integral is well-defined: since $M$ is star-shaped with basepoint $0$, for each $x \in M$ the ray $\{tx : t \in [0,1]\}$ lies entirely in $M$, so $\omega_I(tx)$ is a smooth function of $t \in [0,1]$. The operator $\iota$ is $\mathbb{R}$-linear and maps smooth forms to smooth forms.
[guided]
We want to prove that every closed $l$-form on $M$ is exact, i.e., $H^l_{\mathrm{dR}}(M) = 0$ for $l \ge 1$. The standard strategy is to construct an explicit chain homotopy — a linear map $\iota : \Omega^l(M) \to \Omega^{l-1}(M)$ satisfying $d \circ \iota + \iota \circ d = \mathrm{id}$. If such $\iota$ exists, then for any closed $\omega$ (with $d\omega = 0$) we get $\omega = d(\iota\omega) + \iota(d\omega) = d(\iota\omega)$, so $\omega$ is exact.
Why does star-shapedness allow us to build $\iota$? A star-shaped domain $M$ can be "contracted" to a point along straight-line paths. The operator $\iota$ integrates $\omega$ along these radial rays, accumulating the $(l-1)$-dimensional content of $\omega$ as $t$ ranges from $0$ to $1$.
The formula looks complicated but has a clean geometric interpretation. For a point $x \in M$, the value $\iota\omega$ at $x$ is obtained by integrating the contraction of $\omega$ with the radial vector field $R_x = \sum_i x_i \partial_{x_i}$ (which points from $0$ toward $x$) along the path $t \mapsto tx$, weighted by $t^{l-1}$ to account for the scaling of $\omega$ under $x \mapsto tx$. More precisely,
\begin{align*}
(\iota\omega)_x = \int_0^1 t^{l-1} (\iota_R \omega)_{tx} \, d\mathcal{L}^1(t),
\end{align*}
where $\iota_R \omega$ denotes interior multiplication of $\omega$ by $R$. Expanding $(\iota_R \omega)_{tx}$ in coordinates using $R_{tx} = \sum_i (tx_i) \partial_{x_i}$ and the formula for interior multiplication of a wedge product yields exactly the written expression.
The integrals converge because $\omega_I$ is smooth and $M$ is star-shaped (so $tx \in M$ for all $t \in [0,1]$), making the integrand smooth in $t$ on a compact interval.
[/guided]
[/step]
[step:Verify the homotopy identity $d \circ \iota + \iota \circ d = \mathrm{id}$]
We show that for every $\omega = \sum_I \omega_I \, dx_I \in \Omega^l(M)$,
\begin{align*}
(d \circ \iota + \iota \circ d)(\omega) = \omega.
\end{align*}
By linearity it suffices to treat a single term $\omega = \omega_I \, dx_I$ for a fixed strictly increasing multi-index $I = (i_1 < \cdots < i_l)$.
We compute $(d \circ \iota + \iota \circ d)(\omega_I \, dx_I)$ by tracking what happens when we apply $\iota$ and $d$ to a monomial form, using the homotopy formula
\begin{align*}
\frac{d}{dt}[t^l \omega_I(tx)] = l \, t^{l-1} \omega_I(tx) + t^l \sum_{m=1}^k x_m (\partial_{x_m} \omega_I)(tx).
\end{align*}
Integrating from $t = 0$ to $t = 1$ gives, by the fundamental theorem of calculus,
\begin{align*}
\omega_I(x) - 0 = \int_0^1 \frac{d}{dt}[t^l \omega_I(tx)] \, d\mathcal{L}^1(t) = l \int_0^1 t^{l-1} \omega_I(tx) \, d\mathcal{L}^1(t) + \sum_{m=1}^k x_m \int_0^1 t^l (\partial_{x_m} \omega_I)(tx) \, d\mathcal{L}^1(t).
\end{align*}
The first integral is precisely the coefficient of $dx_I$ in $d(\iota(\omega_I \, dx_I))$ (obtained by differentiating the $\iota$ expression and collecting $dx_I$ terms). The second integral is precisely the coefficient of $dx_I$ in $\iota(d(\omega_I \, dx_I))$ (obtained by applying $\iota$ to $d\omega_I \wedge dx_I = \sum_m (\partial_{x_m}\omega_I) dx_m \wedge dx_I$ and collecting $dx_I$ terms). Since the sum of the two expressions equals $\omega_I(x)$, we conclude
\begin{align*}
(d \circ \iota + \iota \circ d)(\omega_I \, dx_I) = \omega_I \, dx_I,
\end{align*}
and by linearity $(d \circ \iota + \iota \circ d)(\omega) = \omega$ for all $\omega \in \Omega^l(M)$.
[guided]
We want to verify $d(\iota\omega) + \iota(d\omega) = \omega$ directly. By linearity, fix $\omega = f \, dx_I$ with $f = \omega_I \in C^\infty(M)$ and $I = (i_1 < \cdots < i_l)$. We compute each piece.
**Computing $\iota(\omega)$.** By the definition of $\iota$,
\begin{align*}
\iota(f \, dx_I) = \sum_{j=1}^l (-1)^{j-1} \left(\int_0^1 t^{l-1} f(tx) \, d\mathcal{L}^1(t)\right) x_{i_j} \, dx_{I \setminus \{i_j\}}.
\end{align*}
**Computing $d(\iota\omega)$.** Let $A_j(x) = \int_0^1 t^{l-1} f(tx) \, d\mathcal{L}^1(t)$. Then
\begin{align*}
d(\iota\omega) = \sum_{j=1}^l (-1)^{j-1} d(A_j(x) \, x_{i_j} \, dx_{I \setminus \{i_j\}}).
\end{align*}
Applying the Leibniz rule and using $d(x_{i_j}) = dx_{i_j}$:
\begin{align*}
d(A_j \cdot x_{i_j} \, dx_{I \setminus \{i_j\}}) = dA_j \wedge x_{i_j} \, dx_{I\setminus\{i_j\}} + A_j \, dx_{i_j} \wedge dx_{I\setminus\{i_j\}}.
\end{align*}
The term $A_j \, dx_{i_j} \wedge dx_{I\setminus\{i_j\}} = A_j \cdot (-1)^{j-1} \, dx_I$ (reinserting $dx_{i_j}$ at position $j$). The term $dA_j \wedge x_{i_j} \, dx_{I\setminus\{i_j\}}$ involves $\partial_{x_m} A_j = \int_0^1 t^l (\partial_{x_m} f)(tx) \, d\mathcal{L}^1(t)$ (differentiating under the integral). Summing over $j$ and collecting the $dx_I$ component gives
\begin{align*}
\text{coefficient of } dx_I \text{ in } d(\iota\omega) = \sum_{j=1}^l (-1)^{j-1} \cdot (-1)^{j-1} A_j(x) = l \int_0^1 t^{l-1} f(tx) \, d\mathcal{L}^1(t).
\end{align*}
(The cross terms $dA_j \wedge x_{i_j} dx_{I\setminus\{i_j\}}$ contribute to basis elements other than $dx_I$ and will cancel with corresponding terms from $\iota(d\omega)$.)
**Computing $\iota(d\omega)$.** We have $d\omega = \sum_{m \notin I} (\partial_{x_m} f) \, dx_m \wedge dx_I$. Applying $\iota$ to each term $(\partial_{x_m} f) \, dx_m \wedge dx_I$ (an $(l+1)$-form with index set $\{m\} \cup I$) and retaining only the $dx_I$ component (which comes from the $j = 1$ term in $\iota$, i.e., removing $dx_m$) gives coefficient
\begin{align*}
\text{coefficient of } dx_I \text{ in } \iota(d\omega) = \sum_{m=1}^k x_m \int_0^1 t^l (\partial_{x_m} f)(tx) \, d\mathcal{L}^1(t).
\end{align*}
**Summing.** The coefficient of $dx_I$ in $d(\iota\omega) + \iota(d\omega)$ is
\begin{align*}
l \int_0^1 t^{l-1} f(tx) \, d\mathcal{L}^1(t) + \sum_{m=1}^k x_m \int_0^1 t^l (\partial_{x_m} f)(tx) \, d\mathcal{L}^1(t) = \int_0^1 \frac{d}{dt}[t^l f(tx)] \, d\mathcal{L}^1(t),
\end{align*}
where we used the product rule: $\frac{d}{dt}[t^l f(tx)] = l t^{l-1} f(tx) + t^l \sum_m x_m (\partial_{x_m} f)(tx)$. By the fundamental theorem of calculus,
\begin{align*}
\int_0^1 \frac{d}{dt}[t^l f(tx)] \, d\mathcal{L}^1(t) = \bigl[t^l f(tx)\bigr]_{t=0}^{t=1} = f(x) - 0 = f(x).
\end{align*}
So $(d \circ \iota + \iota \circ d)(f \, dx_I) = f \, dx_I$. By linearity, $(d \circ \iota + \iota \circ d)(\omega) = \omega$ for all $\omega \in \Omega^l(M)$.
[/guided]
[/step]
[step:Conclude that every closed $l$-form is exact for $l \ge 1$]
Let $l \ge 1$ and let $\omega \in \Omega^l(M)$ satisfy $d\omega = 0$. Applying the homotopy identity established in the previous step,
\begin{align*}
\omega = (d \circ \iota + \iota \circ d)(\omega) = d(\iota\omega) + \iota(d\omega) = d(\iota\omega) + \iota(0) = d(\iota\omega).
\end{align*}
Thus $\omega = d(\iota\omega)$ is exact, with the explicit primitive $\iota\omega \in \Omega^{l-1}(M)$. Since every closed $l$-form is exact, $H^l_{\mathrm{dR}}(M) = \ker d / \mathrm{im}\, d = 0$ for all $l \ge 1$. Together with the computation $H^0_{\mathrm{dR}}(M) \cong \mathbb{R}$ from the first step, this completes the proof.
[/step]