The statement about the convergence of $\sqrt{n}(g(\overline{X}) - g(\mu))$ can also be found under the name "Cramér's Theorem", see for example "A Course in Large Sample Theory" (Thomas Ferguson, 1996), Theorem 7 (Cramer). It needs the distribution of $\sqrt{n}(\overline{X}-\mu)$ which will often be given by the CTL. Strictly speaking the CTL is only a result on the asymptotic distribution of the sample mean, NOT of functions of sample means.

while $\theta_0$ theoretically should be close to the true parameters $\theta$ you have to take an educated guess in practice, since $\theta$ is unknown. as a naive approach if you don't have any prior informations: why don't you just try a range of different starting values and see if the sequence of estimators from the optimization algorithm always converge to the "same" final estimator? (keep in mind that if you know that the log likelihood function is strictly concave, then the maximum likelihood estimator will be unique.)

i would be happy about a discussion topic on this matter.

as a first try I would also suggest to model the situation using a linear regression. with $p=50$ variables a best subset selection is nowadays feasible (see for example arxiv.org/abs/1803.01454 and its introduction for more state of the art algorithms). however, in any way you should be clear about your goals: are you interested in variable selection or in prediction?

The OLS estimator does not need to be the only BLUE estimator. For example, the maximum likelihood estimator in a regression setup with normal distributed errors is BLUE too, since the closed form of the estimator is identical to the OLS (but as a method, ML-estimation is clearly different from OLS.). The Gauss–Markov Theorem however tells you that in the class of linear unbiased estimators you don't have too look further than OLS, since every other estimator in this class can not do better under the assumptions.

there are indeed different forms of modelling the influence of the factor variable and different corresponding $X$-matrices. however you want to measure two effects (the effect of eacht of the levels of the gender-factor) and hence you need two columns. Try for example: model <- lm(salary ~ gender, data = df,contrasts = list(gender = "contr.sum")) and take a look at the model.matrix(model).