  We present a novel, iterative method using an empirical Bayesian approach for
modeling the limb darkened WASP-121b transit from the TESS light curve. Our
method is motivated by the need to improve $R_{p}/R_{\ast}$ estimates for
exoplanet atmosphere modeling, and is particularly effective with the limb
darkening (LD) quadratic law requiring no prior central value from stellar
atmospheric models. With the non-linear LD law, the method has all the
advantages of not needing atmospheric models but does not converge. The
iterative method gives a different $R_{p}/R_{\ast}$ for WASP-121b at a
significance level of 1$\sigma$ when compared with existing non-iterative
methods. To assess the origins and implications of this difference, we generate
and analyze light curves with known values of the limb darkening coefficients
(LDCs). We find that non-iterative modeling with LDC priors from stellar
atmospheric models results in an inconsistent $R_{p}/R_{\ast}$ at 1.5$\sigma$
level when the known LDC values are as those previously found when modeling
real data by the iterative method. In contrast, the LDC values from the
iterative modeling yields the correct value of $R_{p}/R_{\ast}$ to within
0.25$\sigma$. For more general cases with different known inputs, Monte Carlo
simulations show that the iterative method obtains unbiased LDCs and correct
$R_{p}/R_{\ast}$ to within a significance level of 0.3$\sigma$. Biased LDC
priors can cause biased LDC posteriors and lead to bias in the $R_{p}/R_{\ast}$
of up to 0.82$\%$, 2.5$\sigma$ for the quadratic law and 0.32$\%$, 1.0$\sigma$
for the non-linear law. Our improvement in $R_{p}/R_{\ast}$ estimation is
important when analyzing exoplanet atmospheres.
