Obtaining the effective Lagrangian for axions from the spontaneous breakdown of Peccei-Quinn symmetry I am reading Section $3$ of this review titled The Strong CP Problem and Axions by R. D. Peccei and also the post here.
A famous solution of the Strong CP problem in QCD is offered by the proposal of the Peccei-Quinn symmetry and its spontaneous breakdown. In a nutshell, the idea is that a global anomalous chiral ${\rm U(1)}$ symmetry, called the Peccei-Quinn symmetry, exists at very high energies. It is also broken at very high energies to produce a low-energy effective Lagrangian $$\mathscr{L}_{\rm eff}=\mathscr{L}_{{\rm SM}}+\Big(\bar{\theta}+\xi\frac{a}{f_a}\Big)\frac{g^2}{32\pi^2}G^a_{\mu\nu}\tilde{G}^{\mu\nu}_a-\frac{1}{2}\partial_\mu a\partial^\mu a+\mathscr{L}_{\rm int}(\partial_\mu a,\psi)\tag{1}$$ where $\bar{\theta}$ is called the effective theta parameter of QCD, $G^a_{\mu\nu}$ are the gluon field strength tensor, $a$ is a Goldstone boson field (called axion), $f_a$ is the scale at which the symmetry is broken and $\xi$ is a dimensionless coupling.

I wonder how the Lagrangian of Eq.$(1)$ can be derived. I tried to sketch a possible way to motivate a derivation but it remains incomplete. Any suggestions about how to proceed further?

My attempt We start by defining the ${\rm U(1)}_{\rm PQ}$ transformation on the Standard Model (SM) fields and also extend the SM with an additional complex scalar field $\phi$ (with a potential $V(\phi)$ and some interaction terms with the SM fields).
Assuming $V(\phi)$ is minimized at $|\phi|=f_a\big/\sqrt{2}$, we can write $\phi$ as $$\phi=\frac{1}{\sqrt{2}}(f_a+\rho)\exp\Big[i\frac{a}{f_a}\Big].\tag{2}$$ When expanded, the gradient term for $\phi$ gives rise to $$(\partial_\mu\phi^*)(\partial^\mu\phi)=\frac{1}{2}(\partial_\mu\rho)(\partial^\mu\rho)+\frac{1}{2}(\partial_\mu a)(\partial^\mu a)+\frac{1}{2f_a^2}(\rho^2+2f_a\rho)(\partial_\mu a)(\partial^\mu a).\tag{3}$$ A typical potential becomes $$V(\phi)=\frac{\lambda}{4}\Big(|\phi|^2-\frac{f_a^2}{2}\Big)^2=\frac{\lambda}{16}\Big(\rho^2+2f_a\rho\Big)^2\tag{4}$$
Question Since $a$ be the axion field, we get the desired kinetic term $-\frac{1}{2}\partial_\mu a\partial^\mu a$ (apart from a sign). But we should also obtain another physical field $\rho$. Why is the field $\rho$ absent in the Lagrangian $(1)$?
Note I do not know of any reference that starts with the high energy Lagrangian and from that systematically derives the effective Lagrangian $(1)$. If there is any reference containing a detailed derivation of $(1)$, please suggest.
 A: There are many ways to generate the topological term
$$(\bar\theta + \theta_a)G\tilde G\,.$$
Let's consider the simplest example, known as the KSVZ mechanism. Let's introduce a new vector-like quark $q$ to the standard model Lagrangian. "Vector-like" means that the quark does not have chiral symmetry, and so it is an electroweak singlet. It is a dirac fermion, so it may have a bare mass, allowing us to evade collider constraints on new strongly coupled particles.
Let $\Phi$ be the axion parent field, so it will obey a (global) $U(1)$ axial symmetry, known as Peccei Quinn symmetry. This new quark will similarly be charged under the global $U(1)$, leading to the following Lagrangian for $\Phi$ and $q$
$$\mathscr{L} = \partial_\mu\Phi\partial^\mu\Phi + \lambda\left(|\Phi|^2 - \frac12f_a\right)^2 + {\rm i}\bar q\gamma^\mu\partial_\mu q - m \bar q q  + (y\Phi \bar q q + \text{h.c.})\,.$$
The scale $f_a$ must be extremely large to avoid constraints, so $\Phi$ may safely be assumed to have acquired its VEV
$$\Phi\to \frac{1}{\sqrt{2}}f_a e^{{\rm i}\theta}\,.$$
In this expression, I've neglected the fact that the radial mode still is a degree of freedom. This choice amounts to saying its excitations will be suppressed in our effective theory of the axion, which is accurate as long as $f_a$ is much bigger than the energy scales we wish to deal with.
Thus, the angular degree of freedom couples to the quark with the effective Lagrangian
$$\mathscr{L} = \frac{1}{2}f_a^2\partial_\mu\theta\partial^\mu\theta + {\rm i}\bar q\gamma^\mu\partial_\mu q - m \bar q q + \frac{f_a}{\sqrt{2}}(y e^{{\rm i} \theta}\bar q q + \text{h.c.})\,.$$
Now, we recall the following non-trivial fact about QCD. A chiral rotation of the quarks $q\to e^{{\rm i}\gamma^5\alpha}q$ introduces the following CP violating term into the QCD Lagrangian
$$-\alpha \frac{g_s^2}{32\pi^2}G\tilde G\,.$$
This originally caused a lot of confusion because this term is a total derivative, and therefore it cannot be derived at any order in perturbation theory. However, in the path integral formulation, it is not hard to check that the chiral rotation of quarks introduces this term through the path integral measure (see this reference).
Thus, if we want the quark mass matrix to be real, and thus conserve CP, we must rotate away the imaginary part of the quark Yukawas. The term $\bar\theta$ is the CP violating phase of the Standard Model, and it can be written
$$\bar\theta = \arg\det M$$
where $M$ is the quark mass matrix.
Earlier, however, we introduced a new source of CP violation: the axion parent field phase $\theta,$ coupled to QCD through a new vector-like quark. Introducing a diagonal pure-imaginary entry to the quark mass matrix leads to
$$\arg\det M\to \bar\theta + \theta\,,$$
and thus we arrive at the final axion Lagrangian
$$\mathscr{L} = \frac{1}{2}f_a^2\partial_\mu\theta\partial^\mu\theta +(\bar\theta + \theta)G\tilde G + ...\,.$$
