What does `jacobian=TRUE` do in model$optimize(data = data_list, jacobian = TRUE)

No — this I learned from Aki now!

jacobian=FALSE gives the optimum (MAP) in the constrained space. jacobian=TRUE gives the optimum (MAP) in the unconstrained space. (This is what the case study says anyway.)

Both calculate in the unconstrained space as this is what Stan does. And what is prior and what is likelihood in your model, if such identities of the terms of the log-posterior are even defined, doesn’t matter.

Theoretically the target probability measure (posterior) is the same in both spaces but practicalities like gradients are with densities, not measures. And the base measure (implicit Lebesque measure, constant density), against which you define your densities in Stan, changes due to the parameter transformation, and you have to compensate for that with the Jacobian. Then one can somehow conceptualize the polarity of the Jacobian binary option in two ways, of which the one used in Stan options seems to be the unintuitive one for many (incl. me, and you).