Why do stellar boundaries exist? What theorems or insights are available to decide wether a star with some given interdependence between its density, pressure and temperature distributions ought to have a boundary at a finite distance from its center? I know that in very ideal situations, we have full clarity. But of course, we should do better than that, no?
As a motivation to this question (Skip this if you don't mind) let me read Hawking & Ellis' (1973) slick attempt to a proof of the upper mass limit for $n=3$-polytrope (spherically symmetric, static) white dwarves, p.304:
*The hydrostatic equilibrium equation reads
\begin{equation}
\frac{dp}{dr}(r)=-\rho(r) GM(r) r^{-2}.
\end{equation}
where $M(r)$ is the mass witin a shell of radius $r$ around the origin.
*Multiply both sides by $r^4$ and integrate over $r$. Do integration by parts on the LHS:
\begin{equation}
p(R)R^4 - 4 \int_0^R p(r)r^3 \text{d}r = - \frac{GM(R)^2}{8\pi}
\end{equation}
If $R$ is the stellar boundary, the first term vanishes. If I'm correct, the argument may also proceed if we can find a sequence of radii $R_n \to \infty$ such that $P(R_n)R_n^4 \to 0$. For the rest of the argument we require that the first term is negligible
*On the other hand
\begin{equation}
\frac{d}{dr}\left(\int_0^r pr'^3\text{d}r\right)^{\frac{3}{4}} = \frac{3}{4}\left(\int_0^r pr'^3\right)^{-\frac{1}{4}}pr^3 = \frac{3}{4}\left(\frac{1}{4}pr^4-\frac{1}{4}\int_0^r \frac{dp}{dr}r'^4\text{d}r'\right)^{-\frac{1}{4}}pr^3 < \frac{3\sqrt{2}}{4}p^{\frac{3}{4}}r^2
\end{equation}
where in the last line we used the negativity of $\frac{dp}{dr}$ (follows immediately from the equilibrium equation).
*Since $p\leq C\rho^{\frac{4}{3}}$, we then have
\begin{equation}
\int_0^R pr^3\text{d}r \leq C \left(\int_0^R\rho r^2\text{d}r\right)^{\frac{4}{3}}=C\left(\frac{M(R)}{4\pi}\right)^{\frac{4}{3}}
\end{equation}
*Together with the first line we deduce that
\begin{equation}
M(R) < \frac{(8C)^{\frac{3}{2}}}{(4\pi)^{\frac{1}{2}}}.
\end{equation}
I've already checked the usual references (Chandrasekhar (1939), Horedt), but not really found anything I need or like. Again, these references seem to discuss only very ideal situations.
EDIT: Many of the comments and answers below invoke too much to the many details of the surrounding physics (radiation/chemical effects) of the problem and therefore lack generality (Imagine gaseous bodies with negligible radiation: e.g. a white dwarf in a universe a few billion years from now. Maybe I shouldn't have said "stellar" boundary in the title./ maybe some gas configurations are ruled out in Nature by formation processes instead of equilibrium constraints). My question can really be interpreted in a more well defined way than that: something like "suppose $p(r)$ and $\rho(r)$ solve the hydrostatic equation throughout the star and suppose $f(\rho(r)) < p(r) < g(\rho(r))$ (where $f$ and $g$ are some specified functions), then $\lim_{r\to R-} \rho(r)=\lim_{r\to R-} p(r)=0$ for some $R>0$?"
2nd EDIT: I just realized that the argument of the mass upper bound easily bypasses the ``quick-pressure-decay" requirement. The point is that from
\begin{equation}
p(R)R^4 - 4 \int_0^R p(r)r^3 \text{d}r = - \frac{GM(R)^2}{8\pi}
\end{equation}
(the first line)
we have that
\begin{equation}
4 \int_0^R p(r)r^3 \text{d}r \geq \frac{GM(R)^2}{8\pi}
\end{equation}
which is an inequality in the right direction. Together with the other computations, we therefore have that for any $R>0$
\begin{equation}
M(R) < \frac{(8C)^{\frac{3}{2}}}{(4\pi)^{\frac{1}{2}}}.
\end{equation}
which implies that even if the star has no boundary, then $M=\lim_{R \to \infty} M(R)$ exists and is smaller than or equal to $\frac{(8C)^{\frac{3}{2}}}{(4\pi)^{\frac{1}{2}}}$.
 A: 
Why do stellar boundaries exist?

Stars do not have a hard boundary or surface in the sense that your question seems to suggest you are thinking.  The sun looks like a nice, discrete sphere due to optical effects.  What we call the surface of the sun is the photosphere.
We call this the surface, because it is what we effectively see when we look at the sun.  However, this surface is just the point where the optical depth approaches unity.  Meaning, it is the region where the ionized gas becomes opaque to visible light photons.
Technically, the sun's atmosphere can be said to encompass what is known as the heliosphere, so we, at Earth, are technically within the sun's atmosphere.  Thus, the upper boundary is more like the termination shock than the photosphere, but this depends upon the question you wish to address (more on that below).
More massive stars have an even more ambiguous surface due to various effects.  For instance, Wolf-Rayet and O-Type stars often have very extended coronas making it difficult to identify a surface.

What theorems or insights are available to decide wether a star with some given interdependence between its density, pressure and temperature distributions ought to have a boundary at a finite distance from its center?

I am going to read into your question because as several comments have already stated, a star cannot be infinite.
Think of the atmosphere of a planet like Earth or Venus.  We generally describe these using a model because in reality, they are not homogenous nor always continuous (i.e., I am thinking about sharp density gradients that can occasionally arise).  The model is often of the mass density and it follows an exponential with the following form:
$$
\rho\left( h \right) = \rho_{o} e^{-h/h_{o}} \tag{1}
$$
where $\rho\left( h \right)$ is the mass density at altitude $h$, $\rho_{o}$ is some reference point or known mass density (e.g., average mass density at sea level might be a nice choice), and $h_{o}$ is the scale height or e-folding distance (i.e., it is a similar quantity to the half-life of radioactive materials).  The $h_{o}$ parameter is often related to the soft upper boundary of an atmosphere.  It can be determined rather easily but its physical significance is the important factor here.  You can see there is no hard boundary anywhere within Equation 1 (i.e., the magnitude $\rho\left( h \right)$ asymptotically approaches zero but never reaches it) but that model does a very good job of describing most planetary atmospheres.
Answer
To answer your question, the upper boundary is often determined from a model (like those you present in your question or Equation 1 in my answer) and the physical interpretation is that at altitudes above $h_{o}$, the density (or whatever other parameter is relevant) is so much lower than everything below that we can approximate it as being small, negligible, or zero depending on the level of accuracy needed for the given problem.
This is not satisfying, I know, but physics is about finding ways to approximate nature without ignoring relevant effects.  So the upper boundary is really determined by the problem you hope to address and your choice of model, since in reality there is no hard upper boundary to a gaseous body like a star (Note:  I am ignoring stellar cores and exotic cases like neutron stars.).
Example
In the case of the sun, note that the density is also described by an exponential model, but a slightly more complex one due to the ionization of the particles.  Regardless, in the low chromosphere, the total hydrogen number density can be on the order of $n_{H} \sim 10^{14} \ cm^{-3}$ or $\sim 10^{20} \ m^{-3}$, which is already four orders of magnitude more tenuous than the photosphere.  Above roughly one solar radii, $R_{\odot}$, the number density falls even more down to $\sim 10^{4} \ cm^{-3}$($\sim 10^{10} \ m^{-3}$) at only an altitude of $\sim 5 \ R_{\odot}$, or roughly ten orders of magnitude.
Thus, as you can see the amount of matter per unit volume starts to become negligibly small above some scale height but does not go to zero.  We cannot model everything perfectly, so the trick is to approximate where the atmosphere no longer matters within the limits of the question you are trying to address.
Side Note:  At the altitude of 1 AU (i.e., roughly the location of Earth's orbit), the number density of the sun's atmosphere has dropped to $\sim 1-10 \ cm^{-3}$($\sim 10^{6}-10^{7} \ m^{-3}$).
A: Let me give an example I just computed:
Preliminaries:
1) Through Poisson's equation the hydrostatic equation can be rewritten as
\begin{equation}
\frac{1}{r^2}\frac{d}{dr}\left(\frac{r^2}{\rho}\frac{dp}{dr}\right)=-4\pi G \rho
\end{equation}
2) We can simplify this by introducing an enthalpy $h$ by integrating
\begin{equation}
\frac{dh}{dr}=\rho \frac{dp}{dr}
\end{equation}
which implies that $h$, like $p$, is nonincreasing in $r$. $h$ is defined up to a constant term which we set by the requirement that $h(R) = 0$ (which is not always possible, but okay in the case which I will consider. BTW, $R=\infty$ is of course possible). The hydrostatic equation now reads
\begin{equation}
\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dh}{dr}\right)=-4\pi G \rho
\end{equation}
3) To facilitate matters, we introduce dimensionless quantities
\begin{equation}
\begin{cases}
& r(x)=Rx \\
& \rho(r(x))=\rho(0) u(x) \\
& h(r(x))=h(0) v(x)
\end{cases}
\end{equation}
so that $u(0)=v(1)=1$ and where we have chosen the scale $R$ so that the equilibrium equation nicely reads
\begin{equation}
\frac{1}{x^2}(x^2v'(x))' = - u(x)
\end{equation}
(primes denote differentiation with respect to $x$).
Theorem
Suppose that $u$ and $v$ are a solution of the above differential equation with the stated boundary conditions and $\forall x > 0$ we have that $v(x)^2=:f(v(x))\leq u(x)\leq g(v(x)):=1-\epsilon + \epsilon v(x)$ ($\epsilon>0$ arbitrarily small), then $v$ has a root for some finite $x^*$. (Note: for a star made of nondegenerate cold matter, we would have $u=v^{\frac{3}{2}}$, so this case is treated presently)
proof:
The proof proceeds by assuming that $v>0$ on the entire positive real axis (so that $v$ also satisfies the hydrostatic equation on that whole range) and then deriving a contradiction.
I leave it to the reader to check that $\lim_{x\to 0+} u(x)=1=\lim_{x\to 0+} v(0)$ implies that $\lim_{x\to 0+} v'(x)=0$. On the interval $I:=(0,x^*)$ where $u$ is positive, we have for sure that $\forall x \in I$
\begin{equation}
x^2v'(x)=-\int_0^x s^2u(s)\text{d}s < 0
\end{equation}
so $v$ is surely decreasing on the same interval. Also note that
\begin{equation}
-1=\lim_{x \to 0+}-u(x)=\lim_{x \to 0+} \left(v''(x)+2\frac{v'(x)}{x}\right)=3\lim_{x \to 0+}v''(x)
\end{equation}
So we can definitely find a $\delta>0$ such that $v(x)\leq 1-\frac{1}{12}x^2$ as long as $x < \delta$. Together with the fact that $v$ is decreasing, this implies $v < v_1:=1-e$ where 
\begin{equation}
e(x)=\begin{cases} &\frac{1}{12}x^2 \text{ when } x < \delta \\ &\frac{1}{12}\delta \text{  otherwise }\end{cases}.
\end{equation}
This implies that $u \leq g(v) < g(v_1)=1-\epsilon e$. So integrating
\begin{equation}
x^{-2}(x^2v')'=-u > 1 - \epsilon e
\end{equation}
gives
\begin{equation}
v(x)>max(0,1-\frac{x^2}{2}+d(x))\geq \begin{cases}& 1-\frac{x^2}{2}+d(x) \text{ when } x < \sqrt{6} \\ & 0 \text{ otherwise.}\end{cases}=:v_2(x)
\end{equation}
where I defined the positive continuous function $d(x):=\int_0^x \frac{\text{d}s}{s^2}\int_0^s \text{d}t(t^2 e(t))$ whose details are not important. Now in turn we have that
\begin{equation}
x^{-2}(x^2v')'=-u \leq -f(v) \leq -f(v_2)=-v_2^2.
\end{equation}
Again integrating, we get that
\begin{equation}
v(x) \leq \begin{cases}
& 1-\int_0^x \frac{\text{d}s}{s^2}\int_0^s (t^2(1-\frac{t^2}{6}+d(t))^2)\text{d}t \text{ when } x < \sqrt{6} \\
& 1-\int_0^\sqrt{6} \frac{\text{d}s}{s^2}\int_0^s (t^2(1-\frac{t^2}{6}+d(t))^2)\text{d}t+(\frac{1}{x}-\frac{1}{\sqrt{6}})\int_0^{\sqrt{6}} (t^2(1-\frac{t^2}{6}+d(t))^2)\text{d}t \text{ otherwise }
\end{cases}
< \begin{cases}
& 1-\int_0^x \frac{\text{d}s}{s^2}\int_0^s (t^2(1-\frac{t^2}{6})^2)\text{d}t \text{ when } x < \sqrt{6} \\
& 1-\underbrace{\int_0^\sqrt{6} \frac{\text{d}s}{s^2}\int_0^s (t^2(1-\frac{t^2}{6})^2)\text{d}t}_{=:a}+(\frac{1}{x}-\frac{1}{\sqrt{6}})\underbrace{\int_0^{\sqrt{6}} (t^2(1-\frac{t^2}{6}+d(t))^2)\text{d}t}_{=:b} \text{ otherwise }
\end{cases}
\end{equation}
where it can be calculated that $a=\frac{19}{35}$ and $b>\frac{16}{35}\sqrt{6}$ with the inequality due to the extra positive contribution of the positive function $d$. So the second half of the upper bound reads
\begin{equation}
v(x)<\frac{16}{35}+(\frac{1}{x}-\frac{1}{\sqrt{6}})b=\frac{16}{35}+(\frac{\sqrt{6}}{x}-1)(\frac{16}{35}+\kappa)
\end{equation}
(when $x > \sqrt{6}$)
where $\kappa>0$ is a small constant. This upper bound has a zero at finite $\tilde{x}$. Hence $v$ has a zero at a finite $x'\leq\tilde{x}$, which is a contradiction. Hence $v$ has a zero at a finite $x^*$ (This logic sounds strange, but is correct)
I have similarly used the tricks which got me this theorem to determine, numerically, that likewise taking $f(v)=0.95 v^3$ (recall that $u=v^3$ is the polytrope of degenerate matter) and $g(v)=v^{2.5}$ also implies a finite radius. Combinations like $(g,f)=(v^{1.5},v^3)$ and $(g,f)=(v^2,v^3)$ do not seem to yield any conclusion (using this method numerically).
