[proofplan] We take an arbitrary direction $\phi \in H^1_0(\Omega)$ and examine the one-variable function obtained by restricting $I_\varepsilon$ to the affine line $u_\varepsilon+t\phi$. Since $u_\varepsilon$ is an unconstrained minimizer in the [Hilbert space](/page/Hilbert%20Space) $H^1_0(\Omega)$, this one-variable function has a minimum at $t=0$, and the two-sided difference quotients force its derivative there to vanish once the derivative exists. We compute the derivative term by term: the Dirichlet and load terms are elementary quadratic and linear differentiations, while the penalty term follows from the identity $B_\varepsilon'=\beta_\varepsilon$ and dominated convergence using the global Lipschitz bound on $\beta_\varepsilon$. The assumption $\varepsilon>0$ fixes the penalization parameter in the notation; the variational argument itself uses only the stated structural hypotheses on $B_\varepsilon$ and $\beta_\varepsilon$. Rearranging the resulting first-variation identity gives the weak penalized equation, and restricting to test functions gives the distributional formulation. [/proofplan] [step:Verify that the penalty term has the required integrability] We use the functional declared in the theorem statement: \begin{align*} I_\varepsilon: H^1_0(\Omega) &\to \mathbb{R} \end{align*} with \begin{align*} I_\varepsilon[v]=\frac{1}{2}\int_\Omega |\nabla v(x)|^2\,d\mathcal L^n(x)+\int_\Omega B_\varepsilon(v(x)-\psi(x))\,d\mathcal L^n(x)-\int_\Omega f(x)v(x)\,d\mathcal L^n(x) \end{align*} for $v\in H^1_0(\Omega)$. Let $\mathcal{L}^1$ denote one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb{R}$. Let $L_\varepsilon \geq 0$ denote a global Lipschitz constant for $\beta_\varepsilon$, so that \begin{align*} |\beta_\varepsilon(a)-\beta_\varepsilon(b)| \leq L_\varepsilon |a-b| \end{align*} for all $a,b \in \mathbb{R}$. Since $\beta_\varepsilon(0)=0$, this gives \begin{align*} |\beta_\varepsilon(s)| \leq L_\varepsilon |s| \end{align*} for every $s \in \mathbb{R}$. Let $w_\varepsilon: \Omega \to \mathbb{R}$ denote the function defined by \begin{align*} w_\varepsilon(x)=u_\varepsilon(x)-\psi(x). \end{align*} Because $u_\varepsilon \in H^1_0(\Omega)$ and $\psi \in H^1(\Omega)$, we have $w_\varepsilon \in H^1(\Omega) \subset L^2(\Omega)$. The same Lipschitz bound also makes the energy finite on each admissible line. Indeed, define $C_\varepsilon:=|B_\varepsilon(0)|$. For every $s\in\mathbb R$, the [Fundamental Theorem of Calculus](/theorems/632) applied to $r\mapsto B_\varepsilon(rs)$ gives \begin{align*} |B_\varepsilon(s)|\leq C_\varepsilon+\frac{L_\varepsilon}{2}|s|^2. \end{align*} Thus $B_\varepsilon(v-\psi)\in L^1(\Omega)$ for every $v\in H^1_0(\Omega)$, since $v-\psi\in L^2(\Omega)$ and $\Omega$ is bounded. Hence \begin{align*} |\beta_\varepsilon(w_\varepsilon(x))| \leq L_\varepsilon |w_\varepsilon(x)| \end{align*} for $\mathcal{L}^n$-a.e. $x \in \Omega$, so $\beta_\varepsilon(w_\varepsilon) \in L^2(\Omega)$. Therefore, for every $\phi \in H^1_0(\Omega) \subset L^2(\Omega)$, the product $\beta_\varepsilon(w_\varepsilon)\phi$ belongs to $L^1(\Omega)$ by the [Cauchy-Schwarz inequality](/theorems/432), and the penalty integral in the claimed weak equation is well-defined. [/step] [step:Differentiate the functional along an arbitrary admissible line] Fix $\phi \in H^1_0(\Omega)$. Define \begin{align*} g_\phi: \mathbb{R} &\to \mathbb{R} \end{align*} by \begin{align*} g_\phi(t)=I_\varepsilon[u_\varepsilon+t\phi]. \end{align*} Since $H^1_0(\Omega)$ is a [vector space](/page/Vector%20Space), $u_\varepsilon+t\phi \in H^1_0(\Omega)$ for every $t \in \mathbb{R}$. Since $u_\varepsilon$ minimizes $I_\varepsilon$ over $H^1_0(\Omega)$, the function $g_\phi$ has a minimum at $t=0$. Once we compute $g_\phi'(0)$, minimality gives the sign of the difference quotients: for $h>0$, $(g_\phi(h)-g_\phi(0))/h\geq 0$, while for $h<0$, $(g_\phi(h)-g_\phi(0))/h\leq 0$. If the two-sided derivative exists, these two inequalities force $g_\phi'(0)=0$. The Dirichlet part satisfies \begin{align*} \frac{d}{dt}\bigg|_{t=0} \frac{1}{2}\int_\Omega |\nabla u_\varepsilon+t\nabla\phi|^2\,d\mathcal{L}^n(x) = \int_\Omega \nabla u_\varepsilon\cdot\nabla\phi\,d\mathcal{L}^n(x). \end{align*} The load part satisfies \begin{align*} \frac{d}{dt}\bigg|_{t=0} \left( -\int_\Omega f(u_\varepsilon+t\phi)\,d\mathcal{L}^n(x) \right) = -\int_\Omega f\phi\,d\mathcal{L}^n(x). \end{align*} It remains to justify the penalty derivative. For $t \neq 0$, define the difference quotient \begin{align*} Q_t: \Omega &\to \mathbb{R} \end{align*} by \begin{align*} Q_t(x)=\frac{B_\varepsilon(w_\varepsilon(x)+t\phi(x))-B_\varepsilon(w_\varepsilon(x))}{t}. \end{align*} For $\mathcal{L}^n$-a.e. $x \in \Omega$, apply the [Fundamental Theorem of Calculus](/theorems/632) to the $C^1$ function $r\mapsto B_\varepsilon(w_\varepsilon(x)+r\phi(x))$ on the interval with endpoints $0$ and $t$, using $B_\varepsilon'=\beta_\varepsilon$. After the change of variables $r=\theta t$ in the one-dimensional integral, with $d\mathcal{L}^1(r)=|t|\,d\mathcal{L}^1(\theta)$ and orientation accounted for by division by $t$, this gives \begin{align*} Q_t(x)=\phi(x)\int_0^1 \beta_\varepsilon(w_\varepsilon(x)+\theta t\phi(x))\,d\mathcal{L}^1(\theta). \end{align*} Thus $Q_t(x) \to \beta_\varepsilon(w_\varepsilon(x))\phi(x)$ for $\mathcal{L}^n$-a.e. $x \in \Omega$ as $t \to 0$. For $|t|\leq 1$, using $\beta_\varepsilon(0)=0$ and the Lipschitz bound gives \begin{align*} |Q_t(x)| \leq |\phi(x)|\int_0^1 L_\varepsilon |w_\varepsilon(x)+\theta t\phi(x)|\,d\mathcal{L}^1(\theta). \end{align*} Hence \begin{align*} |Q_t(x)| \leq L_\varepsilon |\phi(x)|\,|w_\varepsilon(x)| + L_\varepsilon |\phi(x)|^2. \end{align*} The right-hand side belongs to $L^1(\Omega)$ because $w_\varepsilon,\phi \in L^2(\Omega)$. The [Dominated Convergence Theorem](/theorems/4) therefore yields \begin{align*} \frac{d}{dt}\bigg|_{t=0} \int_\Omega B_\varepsilon(w_\varepsilon+t\phi)\,d\mathcal{L}^n(x) = \int_\Omega \beta_\varepsilon(w_\varepsilon)\phi\,d\mathcal{L}^n(x). \end{align*} [guided] We fix a direction $\phi \in H^1_0(\Omega)$ and reduce the variational problem to a one-dimensional calculus problem. The admissible perturbations are \begin{align*} u_\varepsilon+t\phi \in H^1_0(\Omega) \end{align*} for $t \in \mathbb{R}$, because $H^1_0(\Omega)$ is a vector space. Define \begin{align*} g_\phi: \mathbb{R} &\to \mathbb{R} \end{align*} by \begin{align*} g_\phi(t)=I_\varepsilon[u_\varepsilon+t\phi]. \end{align*} Since $u_\varepsilon$ minimizes $I_\varepsilon$, the function $g_\phi$ has a minimum at $t=0$. We do not need to quote a separate one-variable theorem here: for $h>0$, minimality gives $(g_\phi(h)-g_\phi(0))/h\geq 0$, and for $h<0$, minimality gives $(g_\phi(h)-g_\phi(0))/h\leq 0$. Therefore, if $g_\phi$ is differentiable at $0$, the common two-sided limit of these quotients must be $g_\phi'(0)=0$. The Dirichlet term is quadratic in $t$. Expanding the integrand gives \begin{align*} |\nabla u_\varepsilon+t\nabla\phi|^2 = |\nabla u_\varepsilon|^2 + 2t\,\nabla u_\varepsilon\cdot\nabla\phi + t^2|\nabla\phi|^2. \end{align*} Since $\nabla u_\varepsilon,\nabla\phi \in L^2(\Omega;\mathbb{R}^n)$, the product $\nabla u_\varepsilon\cdot\nabla\phi$ belongs to $L^1(\Omega)$ by the [Cauchy-Schwarz inequality](/theorems/432). Hence \begin{align*} \frac{d}{dt}\bigg|_{t=0} \frac{1}{2}\int_\Omega |\nabla u_\varepsilon+t\nabla\phi|^2\,d\mathcal{L}^n(x) = \int_\Omega \nabla u_\varepsilon\cdot\nabla\phi\,d\mathcal{L}^n(x). \end{align*} The load term is linear in $t$. Since $f,\phi \in L^2(\Omega)$, the product $f\phi$ belongs to $L^1(\Omega)$, and therefore \begin{align*} \frac{d}{dt}\bigg|_{t=0} \left( -\int_\Omega f(u_\varepsilon+t\phi)\,d\mathcal{L}^n(x) \right) = -\int_\Omega f\phi\,d\mathcal{L}^n(x). \end{align*} The penalty term is the only point requiring a convergence argument. Let \begin{align*} w_\varepsilon: \Omega &\to \mathbb{R} \end{align*} be defined by \begin{align*} w_\varepsilon(x)=u_\varepsilon(x)-\psi(x). \end{align*} Then $w_\varepsilon \in H^1(\Omega)\subset L^2(\Omega)$. For $t\neq 0$, define \begin{align*} Q_t: \Omega &\to \mathbb{R} \end{align*} by \begin{align*} Q_t(x)=\frac{B_\varepsilon(w_\varepsilon(x)+t\phi(x))-B_\varepsilon(w_\varepsilon(x))}{t}. \end{align*} For each point $x$ where $w_\varepsilon(x)$ and $\phi(x)$ are finite, we apply the one-dimensional [Fundamental Theorem of Calculus](/theorems/632) to the $C^1$ function $r \mapsto B_\varepsilon(w_\varepsilon(x)+r\phi(x))$ on the interval between $0$ and $t$, using $B_\varepsilon'=\beta_\varepsilon$. This gives an integral over the variable $r$. We then set $r=\theta t$, so $\theta$ ranges over $[0,1]$ and one-dimensional Lebesgue measure transforms as $d\mathcal{L}^1(r)=|t|\,d\mathcal{L}^1(\theta)$, with the orientation encoded by the factor $1/t$ in the difference quotient. Hence \begin{align*} Q_t(x)=\phi(x)\int_0^1 \beta_\varepsilon(w_\varepsilon(x)+\theta t\phi(x))\,d\mathcal{L}^1(\theta). \end{align*} As $t\to 0$, continuity of $\beta_\varepsilon$ gives the pointwise limit \begin{align*} Q_t(x)\to \beta_\varepsilon(w_\varepsilon(x))\phi(x) \end{align*} for $\mathcal{L}^n$-a.e. $x \in \Omega$. To pass this pointwise limit through the integral, we need an $L^1$ majorant independent of small $t$. Let $L_\varepsilon$ be a global Lipschitz constant for $\beta_\varepsilon$. Since $\beta_\varepsilon(0)=0$, we have \begin{align*} |\beta_\varepsilon(s)|\leq L_\varepsilon |s| \end{align*} for every $s \in \mathbb{R}$. Therefore, for $|t|\leq 1$, \begin{align*} |Q_t(x)| \leq |\phi(x)|\int_0^1 L_\varepsilon |w_\varepsilon(x)+\theta t\phi(x)|\,d\mathcal{L}^1(\theta). \end{align*} Using $|w_\varepsilon(x)+\theta t\phi(x)|\leq |w_\varepsilon(x)|+|\phi(x)|$, we obtain \begin{align*} |Q_t(x)| \leq L_\varepsilon |\phi(x)|\,|w_\varepsilon(x)| + L_\varepsilon |\phi(x)|^2. \end{align*} This majorant belongs to $L^1(\Omega)$ because $w_\varepsilon,\phi \in L^2(\Omega)$ and products of $L^2$ functions are $L^1$ by the [Cauchy-Schwarz inequality](/theorems/432). The [Dominated Convergence Theorem](/theorems/4) applies and yields \begin{align*} \frac{d}{dt}\bigg|_{t=0} \int_\Omega B_\varepsilon(w_\varepsilon+t\phi)\,d\mathcal{L}^n(x) = \int_\Omega \beta_\varepsilon(w_\varepsilon)\phi\,d\mathcal{L}^n(x). \end{align*} [/guided] [/step] [step:Use the vanishing first variation to obtain the weak equation] Combining the three derivative computations gives \begin{align*} g_\phi'(0) = \int_\Omega \nabla u_\varepsilon\cdot\nabla\phi\,d\mathcal{L}^n(x) + \int_\Omega \beta_\varepsilon(u_\varepsilon-\psi)\phi\,d\mathcal{L}^n(x) - \int_\Omega f\phi\,d\mathcal{L}^n(x). \end{align*} Since $g_\phi$ has a minimum at $0$, we have $g_\phi'(0)=0$. Therefore \begin{align*} \int_\Omega \nabla u_\varepsilon\cdot\nabla\phi\,d\mathcal{L}^n(x) + \int_\Omega \beta_\varepsilon(u_\varepsilon-\psi)\phi\,d\mathcal{L}^n(x) = \int_\Omega f\phi\,d\mathcal{L}^n(x). \end{align*} Because $\phi \in H^1_0(\Omega)$ was arbitrary, this proves the weak formulation for every test direction in $H^1_0(\Omega)$. [/step] [step:Interpret the weak identity as a distributional equation] Let $C_c^\infty(\Omega)$ denote the space of smooth functions $\varphi:\Omega\to\mathbb R$ with compact support contained in $\Omega$, and let $\mathcal{D}'(\Omega)$ denote its [topological dual](/page/Topological%20Dual), the space of distributions on $\Omega$. Fix $\varphi \in C_c^\infty(\Omega)$. Since $C_c^\infty(\Omega)\subset H^1_0(\Omega)$, the weak identity applies with $\phi=\varphi$: \begin{align*} \int_\Omega \nabla u_\varepsilon\cdot\nabla\varphi\,d\mathcal{L}^n(x) + \int_\Omega \beta_\varepsilon(u_\varepsilon-\psi)\varphi\,d\mathcal{L}^n(x) = \int_\Omega f\varphi\,d\mathcal{L}^n(x). \end{align*} By the definition of the distributional Laplacian, the distribution $-\Delta u_\varepsilon$ acts on $\varphi$ by \begin{align*} (-\Delta u_\varepsilon)(\varphi) = \int_\Omega \nabla u_\varepsilon\cdot\nabla\varphi\,d\mathcal{L}^n(x). \end{align*} The functions $\beta_\varepsilon(u_\varepsilon-\psi)$ and $f$ both belong to $L^2(\Omega)\subset L^1_{\mathrm{loc}}(\Omega)$, so they define regular distributions on $\Omega$. Hence the preceding identity is exactly \begin{align*} -\Delta u_\varepsilon+\beta_\varepsilon(u_\varepsilon-\psi)=f \end{align*} in $\mathcal{D}'(\Omega)$. This completes the proof. [/step]