  Stellar limb darkening affects a wide range of astronomical measurements and
is frequently modelled with a parametric model using polynomials in the cosine
of the angle between the line of sight and the emergent intensity.
Two-parameter laws are particularly popular for cases where one wishes to fit
freely for the limb darkening coefficients (i.e. an uninformative prior) due to
the compact prior volume and the fact that more complex models rarely obtain
unique solutions with present data. In such cases, we show that the two limb
darkening coefficients are constrained by three physical boundary conditions,
describing a triangular region in the two-dimensional parameter space. We show
that uniformly distributed samples may be drawn from this region with optimal
efficiency by a technique developed by computer graphical programming:
triangular sampling. Alternatively, one can make draws using a uniform,
bivariate Dirichlet distribution. We provide simple expressions for these
parametrizations for both techniques applied to the case of quadratic,
square-root and logarithmic limb darkening laws. For example, in the case of
the popular quadratic law, we advocate fitting for q_1 = (u_1+u_2)^2 and q_2 =
0.5u_1(u_1+u_2)^{-1} with uniform priors in the interval [0,1] to implement
triangular sampling easily. Employing these parametrizations allows one to
derive model parameters which fully account for our ignorance about the
intensity profile, yet never explore unphysical solutions, yielding robust and
realistic uncertainty estimates. Furthermore, in the case of triangular
sampling with the quadratic law, our parametrization leads to significantly
reduced mutual correlations and provides an alternative geometric explanation
as to why naively fitting the quadratic limb darkening coefficients
precipitates strong correlations in the first place.
