I have a cost function \chi^2(p) \sim \sum_i^{N_\text{im}} [\beta_i(p) - \langle \beta(p)\rangle]^2 and associated likelihood \mathscr{L}(p) \sim e^{-\chi^2/2}.

Given a parameter p, the set of \beta are calculated from a blackbox function, which involves a series of only mathematical functions.

I would like to use Julia’s Turing to maximize the likelihood function and estimate p, assuming it is drawn from a normal distribution.

I have gone through multiple tutorials and guides about Turing, yet I am pretty lost how to do this. Could anyone give me some pointers how this is done? I appreciate any help.

As a probabilistic programming language, Turing is very useful to turn a model into a likelihood function and perform inference of unobserved variables. Here you already have the likelihood function and there are no unobserved variables, so I think your problem is a pure optimization problem. And more precisely a least-squares problem.
See this post for an overview of techniques: