Phase transition patterns for coupled complex scalar fields at finite temperature and density

The thermodynamical potential density is defined as

$$V_{\rm eff}(\mu,T)=-\frac{T}{{\cal V}}\ln Z(\beta,\mu),$$

(4)

where ${\cal V}$ is the space volume and $Z(\beta,\mu)$ is the grand partition function,

$$Z(\beta,\mu)={\rm tr}\left[e^{-\beta({\hat{H}}-\mu_{a}{\hat{Q}}_{a})}\right],$$

(5)

where $\beta=1/T$, ${\hat{H}}$ is the Hamiltonian operator, ${\hat{Q}}_{a}$ denotes the conserved charge operators and $\mu_{a}$ are the corresponding chemical potentials. The grand partition function can be written as a functional integral over the fields as usual in quantum field theory Kapusta:2006pm . In the integral functional field form, the grand partition function for the model Eq. (1) then becomes

$\displaystyle Z(\beta,\mu)=\int_{\rm periodic}D\phi_{1}D\phi_{2}D\psi_{1}D\psi_{2}\,e^{-S_{\rm Eucl}},$

(6)

where the functional integrals over the fields are performed under the periodic boundary conditions in imaginary Euclidean time, $\phi_{i}({\bf x},0)=\phi_{i}({\bf x},\beta)$ and $\psi_{i}({\bf x},0)=\psi_{i}({\bf x},\beta)$, and the Euclidean action $S_{\rm Eucl}$ is given by

	$\displaystyle S_{\rm Eucl}$	$\displaystyle=$	$\displaystyle\int_{0}^{\beta}d\tau\int d^{3}x\left[\frac{1}{2}\left(\frac{\partial\phi_{1}}{\partial\tau}-i\mu_{\phi}\phi_{2}\right)^{2}\right.$		(7)
		$\displaystyle+$	$\displaystyle\left.\frac{1}{2}\left(\frac{\partial\phi_{2}}{\partial\tau}+i\mu_{\phi}\phi_{1}\right)^{2}\right.$
		$\displaystyle+$	$\displaystyle\left.\frac{1}{2}(\nabla\phi_{1})^{2}+\frac{1}{2}(\nabla\phi_{2})^{2}+\frac{m_{\phi}^{2}}{2}(\phi_{1}^{2}+\phi_{2}^{2})\right.$
		$\displaystyle+$	$\displaystyle\left.\frac{\lambda_{\phi}}{24}(\phi_{1}^{2}+\phi_{2}^{2})^{2}\right.$
		$\displaystyle+$	$\displaystyle\left.\frac{1}{2}\left(\frac{\partial\psi_{1}}{\partial\tau}-i\mu_{\psi}\psi_{2}\right)^{2}+\frac{1}{2}\left(\frac{\partial\psi_{2}}{\partial\tau}+i\mu_{\psi}\psi_{1}\right)^{2}\right.$
		$\displaystyle+$	$\displaystyle\left.\frac{1}{2}(\nabla\psi_{1})^{2}+\frac{1}{2}(\nabla\psi_{2})^{2}+\frac{m_{\psi}^{2}}{2}(\psi_{1}^{2}+\psi_{2}^{2})\right.$
		$\displaystyle+$	$\displaystyle\left.\frac{\lambda_{\psi}}{24}(\psi_{1}^{2}+\psi_{2}^{2})^{2}\right.$
		$\displaystyle+$	$\displaystyle\left.\frac{\lambda}{4}\left(\phi_{1}^{2}+\phi_{2}^{2}\right)\left(\psi_{1}^{2}+\psi_{2}^{2}\right)\right].$

The expression (6) with Eq. (7) generalizes for the present case of a system with two complex coupled scalar fields the one obtained in the case of one only complex scalar field case, as given, e.g., in Ref. Kapusta:2006pm (see also Ref. Benson:1991nj ).

The general implementation of the OPT in the Lagrangian density happens through an interpolation defined as (see, e.g., Refs. Rosa:2016czs ; Farias:2021ult and references there in)

$$\mathcal{L}\rightarrow\mathcal{L}^{\delta}=(1-\delta)\mathcal{L}_{0}(\eta)+\delta\mathcal{L},$$

(8)

where $\mathcal{L}_{0}$ is the Lagrangian density of the free theory, which is modified by an arbitrary mass parameter $\eta$.

The standard interpolation procedure given by Eq. (8) gives for our model the following Lagrangian density

	$\displaystyle\mathcal{L}^{\delta}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left[\left(\partial_{\mu}\phi_{1}\right)^{2}+\left(\partial_{\mu}\phi_{2}\right)^{2}+\left(\partial_{\mu}\psi_{1}\right)^{2}+\left(\partial_{\mu}\psi_{2}\right)^{2}\right]$		(9)
		$\displaystyle-$	$\displaystyle\frac{m_{\phi}^{2}}{2}\left(\phi_{1}^{2}+\phi_{2}^{2}\right)-\frac{m_{\psi}^{2}}{2}\left(\psi_{1}^{2}+\psi_{2}^{2}\right)$
		$\displaystyle-$	$\displaystyle\frac{\delta}{4!}\lambda_{\phi}\left(\phi_{1}^{2}+\phi_{2}^{2}\right)^{2}-\frac{\delta}{4!}\lambda_{\psi}\left(\psi_{1}^{2}+\psi_{2}^{2}\right)^{2}$
		$\displaystyle-$	$\displaystyle\frac{\delta}{4}\lambda\left(\phi_{1}^{2}+\phi_{2}^{2}\right)\left(\psi_{1}^{2}+\psi_{2}^{2}\right)$
		$\displaystyle-$	$\displaystyle\left(1-\delta\right)\frac{\eta_{\phi}^{2}}{2}\left(\phi_{1}^{2}+\phi_{2}^{2}\right)-\left(1-\delta\right)\frac{\eta_{\psi}^{2}}{2}\left(\psi_{1}^{2}+\psi_{2}^{2}\right)$
		$\displaystyle=$	$\displaystyle\frac{1}{2}\left[\left(\partial_{\mu}\phi_{1}\right)^{2}+\left(\partial_{\mu}\phi_{2}\right)^{2}+\left(\partial_{\mu}\psi_{1}\right)^{2}+\left(\partial_{\mu}\psi_{2}\right)^{2}\right]$
		$\displaystyle-$	$\displaystyle\frac{\Omega_{\phi}^{2}}{2}\left(\phi_{1}^{2}+\phi_{2}^{2}\right)-\frac{\Omega_{\psi}^{2}}{2}\left(\psi_{1}^{2}+\psi_{2}^{2}\right)$
		$\displaystyle-$	$\displaystyle\frac{1}{4!}\delta\lambda_{\phi}\left(\phi_{1}^{2}+\phi_{2}^{2}\right)^{2}-\frac{1}{4!}\delta\lambda_{\psi}\left(\psi_{1}^{2}+\psi_{2}^{2}\right)^{2}$
		$\displaystyle-$	$\displaystyle\frac{1}{4}\delta\lambda\left(\phi_{1}^{2}+\phi_{2}^{2}\right)\left(\psi_{1}^{2}+\psi_{2}^{2}\right)$
		$\displaystyle+$	$\displaystyle\frac{1}{2}\delta\eta_{\phi}^{2}\left(\phi_{1}^{2}+\phi_{2}^{2}\right)+\frac{1}{2}\delta\eta_{\psi}^{2}\left(\psi_{1}^{2}+\psi_{2}^{2}\right),$

with $\Omega_{\phi}^{2}=m_{\phi}^{2}+\eta_{\phi}^{2}$, $\Omega_{\psi}^{2}=m_{\psi}^{2}+\eta_{\psi}^{2}$, where $\eta_{\phi,\psi}$ are mass parameters determined through a variational procedure (see below) and $\delta$ is a dimensionless parameter used as a bookkeeping parameter only to keep track of the order that the OPT is implemented and $\delta$ is set to one at the end. Note that the OPT interpolation changes the usual Feynman rules of the theory. The quartic vertices are changed to

	$\displaystyle-i\lambda_{\phi}\to-i\delta\lambda_{\phi},$
	$\displaystyle-i\lambda_{\psi}\to-i\delta\lambda_{\psi},$
	$\displaystyle-i\lambda\to-i\delta\lambda,$		(10)

and the OPT also leads to the additional vertices that come from the interpolation procedure and that are quadratic in the fields, which comes from the last two terms appearing in Eq. (9). Their Feynman rules are simply

$\displaystyle i\delta\eta_{\phi}^{2}\;\;\;{\rm and}\;\;i\delta\eta_{\psi}^{2},$

(11)

while the bare propagators for the fields have masses replaced in the OPT procedure by $m_{\phi}^{2}\to\Omega_{\phi}^{2}$ and $m_{\psi}^{2}\to\Omega_{\psi}^{2}$.

All calculations are carried out similarly as done in perturbation theory and can be evaluated at any order in $\delta$. Hence, up to this stage the results remain strictly perturbative and very similar to the ones obtained via an ordinary perturbative calculation. Since all quantities evaluated at any finite order $\delta^{k}$ in the OPT depends explicitly on $\eta_{\phi,\psi}$, these parameters need to be fixed appropriately. It is through the freedom in fixing $\eta_{\phi,\psi}$ that nonperturbative results can be generated in the OPT. Since $\eta_{\phi,\psi}$ do not belong to the original theory, these parameters need to be fixed considering some appropriate procedure. For instance, one may fix them by requiring that a given physical quantity, which is been calculated perturbatively to order-$\delta^{k}$, to be evaluated at the value where it is less sensitive to this parameter. This criterion is known as the principle of minimal sensitivity (PMS) Stevenson:1981vj . In this work we consider the effective thermodynamic potential (ETP), which is evaluated to order-$\delta^{k}$, $V_{\rm eff}^{\left(\delta^{k}\right)}$, as the appropriate quantity to be optimized and as considered in most of the OPT works in general¹¹1See, for example, Refs. Farias:2008fs ; Rosa:2016czs ; Yukalov:2019nhu , for examples of other different quantities that can be optimized, different optimization methods and a comparison between the results.. In this case, the PMS criterion translates into the variational relation

$$\frac{\partial V_{\rm eff}^{(\delta^{k})}}{\partial\eta_{\phi,\psi}}\Bigr{|}_{\eta_{\phi,\psi}=\bar{\eta}_{\phi,\psi},\delta=1}=0.$$

(12)

There can also be other optimization procedures that can be applied to fix the OPT mass parameters, but they have shown to be equivalent to the PMS one, including the convergence properties (see, see for instance, Ref. Rosa:2016czs for a discussion on these issues). The optimum value for $\bar{\eta}_{\phi}$ and $\bar{\eta}_{\psi}$ derived from Eq. (12) are now nontrivial functions of the original parameters of the theory. In particular, $\bar{\eta}_{\phi}$ and $\bar{\eta}_{\psi}$ become explicit functions of the couplings and, as a consequence of this, nonperturbative results are generated.

By shifting the fields around their respective background expectations values in Eq. (9), which can be taken along $\phi_{1}$ and $\psi_{1}$ without loss of generality, then $\phi_{1}\to\phi_{1}^{\prime}=\phi_{1}+\phi_{0}$ and $\psi_{1}\to\psi_{1}^{\prime}=\psi_{1}+\psi_{0}$, with $\langle\phi_{1}\rangle=\langle\phi_{2}\rangle=\langle\psi_{1}\rangle=\langle\psi_{2}\rangle=0$ and $\langle\phi_{1}^{\prime}\rangle=\phi_{0}$ and $\langle\psi_{1}^{\prime}\rangle=\psi_{0}$, we can, thus, derive the effective potential at first order in the OPT, i.e., at first order in $\delta$. All terms contributing to the ETP to first order in the OPT are given in Fig. 1. They are all explicitly derived in the Appendix A and the renormalized ETP at first order in the OPT is then given by

	$\displaystyle V_{\rm eff,R}^{(\delta)}$	$\displaystyle=$	$\displaystyle\frac{m_{\phi}^{2}-\mu_{\phi}^{2}}{2}\phi_{0}^{2}+\frac{m_{\psi}^{2}-\mu_{\psi}^{2}}{2}\psi_{0}^{2}$		(13)
		$\displaystyle+$	$\displaystyle\frac{\lambda_{\phi}}{4!}\phi_{0}^{4}+\frac{\lambda_{\psi}}{4!}\psi_{0}^{4}+\frac{\lambda}{4}\phi_{0}^{2}\psi_{0}^{2}$
		$\displaystyle+$	$\displaystyle Y(\Omega_{\phi},T,\mu_{\phi})+Y(\Omega_{\psi},T,\mu_{\psi})$
		$\displaystyle+$	$\displaystyle\left(\frac{\lambda_{\phi}}{3}\phi_{0}^{2}+\frac{\lambda}{2}\psi_{0}^{2}-\eta_{\phi}^{2}\right)X(\Omega_{\phi},T,\mu_{\phi})$
		$\displaystyle+$	$\displaystyle\left(\frac{\lambda_{\psi}}{3}\psi_{0}^{2}+\frac{\lambda}{2}\phi_{0}^{2}-\eta_{\psi}^{2}\right)X(\Omega_{\psi},T,\mu_{\psi})$
		$\displaystyle+$	$\displaystyle\frac{\lambda_{\phi}}{3}X^{2}(\Omega_{\phi},T,\mu_{\phi})+\frac{\lambda_{\psi}}{3}X^{2}(\Omega_{\psi},T,\mu_{\psi})$
		$\displaystyle+$	$\displaystyle\lambda X(\Omega_{\phi},T,\mu_{\phi})X(\Omega_{\psi},T,\mu_{\psi}),$

where the functions $Y(\Omega_{i},T,\mu_{i})$ and $X(\Omega_{i},T,\mu_{i})$ haven been defined in the Appendix A and given by Eqs.(45) and (46), respectively.

Applying the PMS procedure Eq. (12) to the renormalized ETP Eq. (13), we obtain that $\bar{\eta}_{\phi}$ and $\bar{\eta}_{\psi}$ are obtained from the coupled equations,

	$\displaystyle\bar{\eta}_{\phi}^{2}$	$\displaystyle=$	$\displaystyle\frac{\lambda_{\phi}}{3}\tilde{\phi}^{2}+\frac{\lambda}{2}\tilde{\psi}^{2}$
		$\displaystyle+$	$\displaystyle\frac{2\lambda_{\phi}}{3}X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}}+\lambda X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}},$
	$\displaystyle\bar{\eta}_{\psi}^{2}$	$\displaystyle=$	$\displaystyle\frac{\lambda_{\psi}}{3}\tilde{\psi}^{2}+\frac{\lambda}{2}\tilde{\phi}^{2}$
		$\displaystyle+$	$\displaystyle\frac{2\lambda_{\psi}}{3}X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}}+\lambda X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}},$

which are to be solved together with the ones defining the background field values $\tilde{\phi}$ and $\tilde{\psi}$, obtained from,

$$\frac{\partial V_{\rm eff,R}}{\partial\phi_{0}}\Bigr{|}_{\phi_{0}=\tilde{\phi},\psi_{0}=\tilde{\psi}}=0,\;\;\;\;\frac{\partial V_{\rm eff,R}}{\partial\psi_{0}}\Bigr{|}_{\phi_{0}=\tilde{\phi},\psi_{0}=\tilde{\psi}}=0,$$

(16)

which give, respectively, the expressions,

	$\displaystyle\tilde{\phi}\left[m_{\phi}^{2}-\mu_{\phi}^{2}+\frac{\lambda_{\phi}}{6}\tilde{\phi}^{2}+\frac{\lambda}{2}\tilde{\psi}^{2}\right.$
	$\displaystyle\left.+\frac{2\lambda_{\phi}}{3}X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}}+\lambda X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}}\right]=0,$
			(17)
	$\displaystyle\tilde{\psi}\left[m_{\psi}^{2}-\mu_{\psi}^{2}+\frac{\lambda_{\psi}}{6}\tilde{\psi}^{2}+\frac{\lambda}{2}\tilde{\phi}^{2}\right.$
	$\displaystyle\left.+\frac{2\lambda_{\psi}}{3}X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}}+\lambda X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}}\right]=0.$
			(18)

Equations (17) and (18) have the trivial solutions, $\tilde{\phi}=\tilde{\psi}=0$, while the coupled gap equations valid when $\tilde{\phi}\neq 0$ and $\tilde{\psi}\neq 0$ are given, respectively, by

	$\displaystyle m_{\phi}^{2}-\mu_{\phi}^{2}+\frac{\lambda_{\phi}}{6}\tilde{\phi}^{2}+\frac{\lambda}{2}\tilde{\psi}^{2}+\frac{2\lambda_{\phi}}{3}X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}}$
	$\displaystyle+\lambda X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}}=0,$		(19)
	$\displaystyle m_{\psi}^{2}-\mu_{\psi}^{2}+\frac{\lambda_{\psi}}{6}\tilde{\psi}^{2}+\frac{\lambda}{2}\tilde{\phi}^{2}+\frac{2\lambda_{\psi}}{3}X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}}$
	$\displaystyle+\lambda X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}}=0.$		(20)

Note that by combining Eqs. (LABEL:pmsetaphi) and (LABEL:pmsetapsi) with Eqs. (19) and (20), we obtain

$$\tilde{\phi}^{2}=\frac{3}{\lambda_{\phi}}\left(m_{\phi}^{2}-\mu_{\phi}^{2}+\bar{\eta}_{\phi}^{2}\right),$$

(21)

and

$$\tilde{\psi}^{2}=\frac{3}{\lambda_{\psi}}\left(m_{\psi}^{2}-\mu_{\psi}^{2}+\bar{\eta}_{\psi}^{2}\right),$$

(22)

which gives that at the critical point for, e.g., in the $\phi$-direction, $\tilde{\phi}(T_{c,\phi},\mu_{c,\phi})=0$, the critical chemical potential gets uniquely fixed by

$$\mu_{\phi,c}^{2}=m_{\phi}^{2}+\bar{\eta}_{\phi}^{2}(T_{c,\phi},\mu_{c,\phi}),$$

(23)

and, equivalently, at the critical point in the $\psi$-direction, $\tilde{\psi}(T_{c,\psi},\mu_{c,\psi})=0$, it leads to

$$\mu_{\psi,c}^{2}=m_{\psi}^{2}+\bar{\eta}_{\psi}^{2}(T_{c,\psi},\mu_{c,\psi}).$$

(24)

Equations (23) and (24) generalize to the present problem the usual condition for Bose-Einstein condensation (BEC) for an ideal Bose gas Kapusta:2006pm , $|\mu|=|m|$ at the critical temperature for BEC. Note that the contribution from the interactions to the BEC transition point enters implicitly in the OPT functions, $\bar{\eta}_{\phi}^{2}$ and $\bar{\eta}_{\psi}^{2}$, in Eqs. (23) and (24).

From the effective thermodynamical potential, we can compute the pressure,

$$P(\mu_{\phi},\mu_{\phi},T)=-V_{\rm eff,R}^{(\delta)}\Bigr{|}_{\bar{\eta}_{\phi},\bar{\eta}_{\psi},\tilde{\phi},\tilde{\psi}},$$

(25)

which is evaluated at the PMS values $\eta_{\phi}=\bar{\eta}_{\phi}$ and $\eta_{\psi}=\bar{\eta}_{\psi}$ and at the VEV values $\phi_{0}=\tilde{\phi}$ and $\psi_{0}=\tilde{\psi}$. Given the pressure, the densities are evaluated as

	$\displaystyle n_{\phi}=\frac{\partial P(\mu_{\phi},\mu_{\phi},T)}{\partial\mu_{\phi}},$		(26)
	$\displaystyle n_{\psi}=\frac{\partial P(\mu_{\phi},\mu_{\phi},T)}{\partial\mu_{\psi}}.$		(27)

Then, making use of Eq. (13) in combination with the PMS equation (12), we find

	$\displaystyle n_{\phi}$	$\displaystyle=$	$\displaystyle\mu_{\phi}\tilde{\phi}^{2}-\frac{\partial Y(\Omega_{\phi},\mu_{\phi},T)}{\partial\mu_{\phi}}\Bigr{|}_{\bar{\eta}_{\phi}},$		(28)
	$\displaystyle n_{\psi}$	$\displaystyle=$	$\displaystyle\mu_{\psi}\tilde{\psi}^{2}-\frac{\partial Y(\Omega_{\psi},\mu_{\psi},T)}{\partial\mu_{\psi}}\Bigr{|}_{\bar{\eta}_{\psi}}.$		(29)

Note that the above equations for $n_{\phi}$ and $n_{\psi}$ are to be solved simultaneously with those for the PMS, Eqs. (LABEL:pmsetaphi)and (LABEL:pmsetapsi), together with those for the background fields, Eqs. (17) and (18).

For our numerical results we will consider for illustration purposes the base parameters values for the couplings: $\lambda_{\phi}=0.018$, $\lambda_{\psi}=0.6$ and $\lambda=-0.03$. Note that though $\lambda<0$, we have that $\lambda_{\phi}\lambda_{\psi}=0.0108>9\lambda^{2}=0.0081$, thus, this choice of couplings satisfies the boundness condition for the potential. In other words, the tree potential is safely inside the stable region. Other choices can be made but the results are qualitatively similar under the conditions considered here. Thus, by considering $m_{\phi}^{2}=m_{\psi}^{2}=m^{2}>0$ and the set of coupling parameters given above, ISB can be shown to happen in the direction of $\phi$ at high temperatures, while $\psi$ remains in the symmetric phase. This is what simple perturbation theory at high temperatures would predict for the two-field coupled complex scalar model in the absence of chemical potentials Farias:2021ult . Note that we can recover the PT results from the OPT interpolated ETP Eq. (13) by setting $\eta_{\phi}=\eta_{\psi}=0$, which then gives the ETP at first order in the coupling constants in the first order PT approximation. In the OPT case, we do obtain though, quantitative differences when compared with the results from PT as we are going to illustrate.

The situation illustrated in the Fig. 3 demonstrates ISB in the direction of the field $\phi$. Starting from a symmetry restored (SR) state at $T=0,\,\mu_{\phi}=\mu_{\psi}=0$, the $\phi$ field eventually acquires a nonvanishing VEV at a critical temperature. Both OPT and PT are compared. One notices that the OPT tends to produce a higher critical temperature than in the PT approximation case. Given the parameters considered, the field $\psi$ remains in the SR phase. The effect of a finite chemical potential on the behavior for the VEV of the $\phi$ field is illustrated in the Fig. 4. Note that for our choice of parameters and for $\mu$ held constant at the values considered, the field $\psi$ remains always with a vanishing background expectation value, $\tilde{\psi}=0$, hence, we do not show it in Figs. 3 and 4. In Fig. 4, we restrict only to show the OPT case, since for PT the results are similar, with the same trend as seen in Fig. 3, leading to smaller critical temperatures when compared to the OPT. From the results shown in Fig. 4, we see that the chemical potential tends to decrease the critical temperature for ISB. Thus, the larger is the chemical potential, the easier is to reach symmetry inversion in the direction of $\phi$, i.e., the $\phi$ field acquires a VEV $\langle\phi\rangle\neq 0$ at lower temperatures as the chemical potential increases. As the field changes smoothly through the critical point, the phase transition associated with ISB here is of second order. The overall behavior for the critical temperature for ISB in the direction of $\phi$, $T_{c,\phi}$, as a function of the chemical potential, is shown in Fig. 5.

The Fig. 5 shows that the higher the chemical potential, the lower the value of the critical temperature. The presence of charge is seen here to favor symmetry inversion (ISB), causing it to happen at a lower critical temperature.

From the mass eigenvalues expressions given in Appendix B, we can define the corresponding Higgs and Goldstone effective modes for $\phi$ and $\psi$ in the context of the OPT. This can be done by introducing in the definitions Eq. (81), the thermal and chemical potential contributions such that Duarte:2011ph ; Farias:2021ult , $M_{H,\phi}^{2}\to M_{H,\phi}^{2}(T,\mu_{\phi},\mu_{\psi})$, $M_{G,i}^{2}\to M_{G,\phi}^{2}T,\mu_{\phi},\mu_{\psi})$ and similarly for the $\psi$ field, where, at first order in the OPT,

	$\displaystyle M_{H,\phi}^{2}(T,\mu_{\phi},\mu_{\psi})$	$\displaystyle=$	$\displaystyle m_{\phi}^{2}-\mu_{\phi}^{2}+\frac{\lambda_{\phi}}{2}\tilde{\phi}^{2}+\frac{\lambda}{2}\tilde{\psi}^{2}$		(30)
		$\displaystyle+$	$\displaystyle\frac{2\lambda_{\phi}}{3}X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}}$
		$\displaystyle+$	$\displaystyle\lambda X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}},$
	$\displaystyle M_{G,\phi}^{2}(T,\mu_{\phi},\mu_{\psi})$	$\displaystyle=$	$\displaystyle m_{\phi}^{2}+\frac{\lambda_{\phi}}{6}\tilde{\phi}^{2}+\frac{\lambda}{2}\tilde{\psi}^{2}$		(31)
		$\displaystyle+$	$\displaystyle\frac{2\lambda_{\phi}}{3}X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}}$
		$\displaystyle+$	$\displaystyle\lambda X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}},$
	$\displaystyle M_{H,\psi}^{2}(T,\mu_{\phi},\mu_{\psi})$	$\displaystyle=$	$\displaystyle m_{\psi}^{2}+\frac{\lambda_{\psi}}{2}\tilde{\psi}^{2}+\frac{\lambda}{2}\tilde{\phi}^{2}$		(32)
		$\displaystyle+$	$\displaystyle\frac{2\lambda_{\psi}}{3}X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}}$
		$\displaystyle+$	$\displaystyle\lambda X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}},$
	$\displaystyle M_{G,\psi}^{2}(T,\mu_{\phi},\mu_{\psi})$	$\displaystyle=$	$\displaystyle m_{\psi}^{2}+\frac{\lambda_{\psi}}{6}\tilde{\psi}^{2}+\frac{\lambda}{2}\tilde{\phi}^{2}$		(33)
		$\displaystyle+$	$\displaystyle\frac{2\lambda_{\psi}}{3}X(\Omega_{\psi},T,\mu_{\psi})\Bigr{|}_{\eta_{\psi}=\bar{\eta}_{\psi}}$
		$\displaystyle+$	$\displaystyle\lambda X(\Omega_{\phi},T,\mu_{\phi})\Bigr{|}_{\eta_{\phi}=\bar{\eta}_{\phi}},$

where $\tilde{\phi}$ and $\tilde{\psi}$ are obtained from the solution of Eq. (16), while $\bar{\eta}_{\phi}$ and $\bar{\eta}_{\psi}$ from the PMS equations (LABEL:pmsetaphi) and (LABEL:pmsetapsi). Note that the Higgs modes must remain positive in both symmetric and broken phases, while vanishing at the critical point for phase transition. The Goldstone modes, on the other hand, must remain massless in the broken phases, according to the Goldstone theorem concerning symmetry breaking of a continuous symmetry, while in the symmetric phase follows the Higgs modes.

Substituting Eqs. (30)-(33) in the mass eigenvalue equation (82)-(85), we obtain the corresponding ones at finite temperature and chemical potential, $\mathcal{M}_{1}^{2}\to m_{H,\phi}^{2},\;\mathcal{M}_{2}^{2}\to m_{G,\phi}^{2},\;\mathcal{M}_{3}^{2}\to m_{H,\psi}^{2},\;\mathcal{M}_{4}^{2}\to m_{G,\psi}^{2}$. These Higgs and Goldstone modes for each of the fields are plotted in Fig. 6 for the cases of vanishing chemical potentials (panel a) and also for nonvanishing chemical potential (panel b). In both cases we see that the Goldstone theorem is correctly reproduced. The fact that the OPT satisfies the Goldstone theorem has also been seen in previous applications Duarte:2011ph ; Farias:2021ult .

In the previous subsection, we have investigated the two coupled complex scalar field at finite temperature and chemical potential with respect to symmetry inversion (ISB). Let us now study the model for symmetry persistence at high temperatures, i.e., we will study the case for SNR in one of the field directions. We continue to use the same set of bare coupling constants as before for clarity, but now we consider that the symmetries for both the fields, in the vacuum, are both broken, $m_{\phi}^{2}<0$ and $m_{\psi}^{2}<0$, such that both fields have a nonvanishing VEV at $T=0$ and $\mu_{\phi}=\mu_{\psi}=0$ initially. In this case, we expect SNR in the direction of $\phi$, while $\psi$ should suffer the usual symmetry restoration (SR) at high temperatures. This is illustrated in Fig. 7.

In Fig. 8 we show the Higgs and Goldstone modes for the fields when $\mu_{\phi}=\mu_{\psi}=0$ (panel a) and when $\mu_{\phi}=\mu_{\psi}=0.5M$ (panel b) for the present case of SNR in the direction of $\phi$ and SR in the direction of $\psi$. As in the case studied in the previous subsection, we also see here that the Goldstone theorem is correctly reproduced.

Let us now turn our attention of how a finite density affects the results. In the previous two subsections we were interested in the effect of a finite chemical potential in the transition patterns displayed by the two coupled complex scalar field system. Here, we are interested in investigating the same effects but now at finite densities. This is mostly motivated by the seminal work performed in Refs. Bernstein:1990kf ; Benson:1991nj . In particular, in Ref. Benson:1991nj it was shown that finite density effects were already able to not only delay symmetry restoration at finite temperature, but also to promote symmetry nonrestoration for a sufficiently large charge (density). The novelty result was that both ISB and SNR could be possible already for a one-field model case. Here, we are interested in studying how the finite density effects will further affect the phase transition pattern when considering the case of more than one coupled complex scalar field. In order to facilitate the comparison with Ref. Benson:1991nj , we will also assume densities for the fields such that the ratio of number density to entropy density is kept fixed. The motivation here, as also in Ref. Benson:1991nj , stems from the fact that in an adiabatic expansion, with the entropy remaining constant, the ratio of charge density over entropy density remains constant. This is like the expected situation in the case of the expanding Universe in the radiation dominated phase when in the absence of entropy production processes. Thus, we assume from now on that the charge density $n_{i}$ of the fields over entropy energy density remains constant, $n_{i}/s={\rm constant}$. Since in the ultrarelativistic case, $s\propto T^{3}$ and $n_{i}\propto T^{3}$, we take Benson:1991nj

$$n_{i}=\tau_{i}T^{3},$$

(34)

By using Eq. (34) in conjunction with the equations defining the densities, Eqs. (28) and (29), together with the PMS and gap equations, we study the behavior of the fields expectations values $\tilde{\phi}$ and $\tilde{\psi}$ as a function of the temperature and fixed ratio $\tau$. In Fig. 9 it is shown the VEV of the fields, $\tilde{\phi}$ (panel a) and $\tilde{\psi}$ (panel b), as a function of temperature. Two representative values of $\tau$ have been used for illustration. Note that for the parameters considered, in the absence of finite charge effects ISB is expected to happen in the direction of the $\phi$ field (see, e.g., Fig. 3). In the presence of finite charge densities, we see that two effects appear. First, the ISB transition can happen earlier, i.e., at a lower critical temperature. Second, there is now also an ISB transition in the $\psi$ field direction, i.e., $\tilde{\psi}$ can now acquire a nonvanishing value at finite temperatures and also here we see that the larger is $\tau$, the lower is the critical temperature for ISB in the direction of $\psi$. Finite charges are then realizing the transition pattern (a) $\to$ (b) $\to$ (d) shown in Fig. 1 in the present case.

In Fig. 10 we now consider the case where both fields are initially in their broken symmetry states, i.e., $\tilde{\phi}\neq 0$ and $\tilde{\psi}\neq 0$ at $T=0$ and $\mu_{\phi}=\mu_{\psi}=0$. For the parameters considered, in the absence of finite charge densities, the transition pattern expected is the one shown in Fig. 7, i.e., SNR in the direction of $\phi$ and usual SR transition in the direction of $\psi$. From Fig. 10 we now see that besides the $\phi$ field remaining in a SNR state, now the finite charge density also induces a SNR in the direction of the $\psi$ field. For the parameters considered in Fig. 10, for a charge density over entropy density ratio $\tau\gtrsim 0.037$, we find that it ceases to exist a critical temperature for SR in the direction of $\psi$ and the field remains in a SNR state.

As a note to be remarked here, even though we have in this part of the work continued to work with a negative value for the intercoupling $\lambda$, in the presence of finite charge densities both situations shown in Figs. 9 and 10 are still realized. ISB and SNR happens independently of the sign of the coupling between $\phi$ and $\psi$. This happens exclusively because of the effect of considering large enough densities. Our results then generalize to the two-field case the situation found in Ref. Benson:1991nj , where it was studied for the one-field case.