The objective of this paper is to quantify uncertainties in thermoacoustic stability analysis with a Helmholtz solver and its adjoint. Thermoacoustic combustion instability may be described by the Helmholtz equation combined with a model for the flame dynamics. Typically, such a formulation leads to an eigenvalue problem in which the eigenvalue appears under nonlinear terms, such as exponentials related to time delays that result from the flame model. Consequently, the standard adjoint sensitivity formulation should be augmented by first- and second-order correction terms that account for the nonlinearities. Such a formulation is developed in the present paper, and applied to the model of a combustion test rig with a premix swirl burner. The uncertainties considered concern plenum geometry, outlet acoustic reflection coefficient, as well as gain and phase of the flame response. The nonlinear eigenvalue problem and its adjoint are solved by an in-house adjoint Helmholtz solver, based on an axisymmetric finite volume approach. In addition to first-order correction terms of the adjoint formulation, which are often used in literature, second-order terms are also taken into account. It is found that one particular second-order term has significant impact on the accuracy of results. Finally, the Probability Density Function of the growth rate in the presence of uncertainties in input paramters is calculated with Monte Carlo simulations. It is found that the second-order adjoint method, while giving quantitative agreement, requires far less compute resources than Monte Carlo sampling for the full nonlinear eigenvalue problem.