text:
# Analysis of 𝑆𝑈(3)_𝑐⊗𝑆𝑈(3)_𝐿⊗𝑈(1)_𝑋 local gauge theory

## 1 Introduction

The remarkable experimental success of the standard model (SM) local gauge group $`G_{SM}SU\left(3\right)_cSU\left(2\right)_LU\left(1\right)_Y`$ with the flavor sector $`SU\left(2\right)_LU\left(1\right)_Y`$ hidden and $`SU\left(3\right)_c`$ confined, lies in its accurate predictions at energies below a few hundreds GeV. However, the SM is not the only model for which this is true and many physicists believe that it does not represent the final theory, but serves merely as an effective theory originating from a more fundamental one. So, extensions of the SM are always welcome.

One can extend the SM either by adding new fermion fields (adding a right-handed neutrino field constitute its simplest extension), by augmenting the scalar sector to more than one higgs representation, or by enlarging the local gauge group. In this last direction, $`SU\left(3\right)_LU\left(1\right)_X`$ as a flavor group has been studied previously by many authors who have explored many possible fermion and higgs-boson representation assignments, either as identical replicas of one family structures as in the SM  or as a multi-family structure which points to a natural explanation of the total number of families in nature.

With regard to the different models in Ref., most of them are plagued with physical inconsistencies such as gauge anomalies, right-handed currents at low energies, unwanted flavor changing neutral currents, violation of universality, etc.. The model in Refs. for three families of quarks and leptons is consistent with the low energy phenomenology and it is anomaly free thanks to the introduction of quarks with exotic electric charges $`4/3`$ and $`5/3`$. On the other hand, the model in Refs., also for three families, is consistent with low energy phenomenology and does not include particles with exotic electric charges.

In this paper we present an analysis of the local gauge group $`SU\left(3\right)_cSU\left(3\right)_LU\left(1\right)_X`$. We find six models which are anomaly free, do not include fermions (quarks and leptons) with exotic electric charges and are consistent with the low energy phenomenology. Two of the models are one family models and are natural extensions of the SM (one of them is an $`E_6`$ subgroup), while the other four are models for three families of quarks and leptons; two of them, up to our knowledge, new in the literature. The models under consideration get their symmetries broken via the most economical set of higgs fields. We analyze also the limit in which the neutral currents reproduce the SM phenomenology.

Our paper is organized in the following way: In section two we introduce the characteristics of the gauge group and present the six different models mentioned above; in section three we describe the scalar sector needed to break the symmetry; in section four we analyze the gauge boson sector paying special attention to the two neutral currents and their mixing, in section five we analyze the fermion masses for one particular model and in the last section we give our conclusions.

## 2 The model

In what follows we assume that the electroweak gauge group is $`SU\left(3\right)_LU\left(1\right)_XSU\left(2\right)_LU\left(1\right)_Y`$. We also assume that the left handed quarks (color triplets) and left-handed leptons (color singlets) transform under the two fundamental representations of $`SU\left(3\right)_L`$ (the 3 and $`3^{}`$). Two classes of models will be discussed: one family models where the anomalies cancel in each family as in the SM, and family models where the anomalies cancel by an interplay between the families. As in the SM, $`SU\left(3\right)_c`$ is vectorlike.

All the models analyzed have the same gauge boson sector, but they differ in their fermion content and may differ in the scalar sector too.

### 2.1 One family models

The most general expression for the electric charge generator in $`SU\left(3\right)_LU\left(1\right)_X`$ is a linear combination of the three diagonal generators of the gauge group

$$Q=aT_{3L}+\frac{2}{\sqrt{3}}bT_{8L}+XI_3,$$
(1)
where $`T_{iL}=\lambda _{iL}/2;\lambda _{iL}`$ being the Gell-Mann matrices for $`SU\left(3\right)_L`$ normalized as Tr$`(\lambda _i\lambda _j)=2\delta _{ij}`$, $`I_3=Dg(1,1,1)`$ is the diagonal $`3\times 3`$ unit matrix, and $`a`$ and $`b`$ are arbitrary parameters to be calculated ahead. Notice that we have absorbed an eventual coefficient for $`X`$ in its definition.

Now, having in mind the canonical iso-doublets for $`SU\left(2\right)_L`$ in one family, we start by defining two $`SU\left(3\right)_L`$ triplets

$$\chi _L=\left(\begin{array}{c}u\\ d\\ q\end{array}\right)_L,\psi _L=\left(\begin{array}{c}e^{}\\ \nu _e\\ l\end{array}\right)_L,$$

where $`q_L`$ and $`l_L`$ are $`SU\left(2\right)_L`$ singlet exotic quark and lepton fields respectively of electric charges to be fixed ahead. This structure implies that $`a=1`$ in Eq.(1) and one gets a one-parameter set of models. Now if the {$`SU\left(3\right)_L,U\left(1\right)_X`$} quantum numbers for $`\chi _L`$ and $`\psi _L`$ are {$`3,X_\chi `$} and {$`3^{},X_\psi `$} respectively, then by using Eq.(1) we have the relationship:

$$X_\chi +X_\psi =Q_q+Q_l=1/3,$$
(2)
where $`Q_q`$ and $`Q_l`$ are the electric charge values of the $`SU\left(2\right)_L`$ singlets $`q`$ and $`l`$ respectively, in units of the absolute value of the electron electric charge.

Now in order to cancel the $`\left[SU\left(3\right)_L\right]^3`$ anomaly, two more $`SU\left(3\right)_L`$ lepton anti-triplets with quantum numbers {$`3^{},X_i\},i=1,2`$, must be introduced (together with their corresponding right-handed charged components). Each one of those multiplets must include one $`SU\left(2\right)_L`$ doublet and one singlet of new leptons. The quarks fields $`u_L^c,d_L^c`$ and $`q_L^c`$ color anti-triplets and $`SU\left(3\right)_L`$ singlets, with $`U\left(1\right)_X`$ quantum numbers $`X_u,X_d`$ and $`X_q`$ respectively, must also be introduced in order to cancel the $`\left[SU\left(3\right)_c\right]^3`$ anomaly. Then the hypercharges $`X_\alpha `$ with $`\alpha =\chi ,\psi ,1,2,u,d,q,\mathrm{}`$ are fixed using Eqs. (1), (2) and the anomaly constraint equations coming from the vertices $`\left[SU\left(3\right)_c\right]^2U\left(1\right)_X,\left[SU\left(3\right)_L\right]^2U\left(1\right)_X,\left[grav\right]^2U\left(1\right)_X`$ and $`\left[U\left(1\right)_X\right]^3`$, which are:

$`\left[SU\left(3\right)_c\right]^2U\left(1\right)_X`$ $`:`$ $`3X_\chi +X_u+X_d+X_q=0`$
$`\left[SU\left(3\right)_L\right]^2U\left(1\right)_X`$ $`:`$ $`3X_\chi +X_\psi +X_1+X_2=0`$
$`\left[grav\right]^2U\left(1\right)_X`$ $`:`$ $`9X_\chi +3X_u+3X_d+3X_q+3X_\psi +3X_1+3X_2+{\displaystyle \underset{singl}{}}X_{ls}=0`$
$`\left[U\left(1\right)_X\right]^3`$ $`:`$ $`9X_\chi ^3+3X_u^3+3X_d^3+3X_q^3+3X_\psi ^3+3X_1^3+3X_2^3+{\displaystyle \underset{singl}{}}X_{ls}^3=0,`$
where $`X_{ls}`$ are the hypercharges of the right-handed charged lepton singlets needed in order to have a consistent field theory.

What we have so far is an infinite number of possible models each one characterized by the parameter $`b`$ in Eq.(1); the value of $`b`$ is the key factor in determining the electric charge of the extra particles in the several models to be presented. We are going to drastically limit this number of possible models by imposing the constraint of excluding models with particles with exotic electric charges; that is, we are going to allow only models with quarks of electric charges $`\pm 2/3`$ and $`\pm 1/3`$ and leptons of electric charges $`\pm 1`$ and 0. We will see that this requirement render us with only two different sets of models $`\left(b=\pm 1/2\right)`$ with equivalent gauge sectors.

#### 2.1.1 Model A

Let us start with a model with an extra down type quark $`q=D`$ of electric charge $`Q_q=Q_D=1/3`$ ($`b=1/2`$) which in turn implies $`Q_l=0`$, that is, $`l_L`$ is a new neutral lepton $`N_{1L}^0`$. Eq.(1) then implies $`X_q=X_d=1/3`$, $`X_u=2/3`$, which combined with the anomaly constraint equations and Eq.(2) gives $`X_\chi =0,X_\psi =1/3,_{singl}X_{ls}=0`$ and $`X_1+X_2=1/3`$. By demanding for leptons of electric charges $`\pm 1`$ and 0 only, we have for the simplest solution that $`X_1=1/3,X_2=2/3`$ and $`X_{ls}=0`$, with this last constraint meaning that we do not need right-handed charged leptons in our simplest anomaly-free model.

Putting all this together we end up with the following multiplet structure for this model:

$$\begin{array}{cccc}& & & \\ & & & \\ \chi _L=\left(\begin{array}{c}u\\ d\\ D\end{array}\right)_L& u_L^c& d_L^c& D_L^c\\ & & & \\ (3,3,0)& (3^{},1,\frac{2}{3})& (3^{},1,\frac{1}{3})& (3^{},1,\frac{1}{3})\\ & & & \end{array}$$

$$\begin{array}{ccc}& & \\ & & \\ \psi _L=\left(\begin{array}{c}e^{}\\ \nu _e\\ N_1^0\end{array}\right)_L& \psi _{1L}=\left(\begin{array}{c}E^{}\\ N_2^0\\ N_3^0\end{array}\right)_L& \psi _{2L}=\left(\begin{array}{c}N_4^0\\ E^+\\ e^+\end{array}\right)_L\\ & & \\ (1,3^{},\frac{1}{3})& (1,3^{},\frac{1}{3})& (1,3^{},\frac{2}{3})\\ & & \end{array},$$

where the numbers inside the parenthesis refer to $`(SU\left(3\right)_c,SU\left(3\right)_L,U\left(1\right)_X)`$ quantum numbers. This anomaly-free structure is the simplest one we can construct for a single family in $`SU\left(3\right)_LU\left(1\right)_X`$. As a matter of fact, the 27 states above are just the 27 states in the fundamental representation of the electroweak-strong unification group $`E_6`$, so this gauge and fermion structure is such that $`SU\left(3\right)_cSU\left(3\right)_LU\left(1\right)_XE_6`$. A phenomenological analysis of this model has been published already in Ref..

#### 2.1.2 Model B

For this model we start with an extra up type quark $`q=U`$ of electric charge $`Q_q=Q_U=2/3\left(b=1/2\right)`$ which in turn implies $`Q_l=1`$, that is, $`l_L`$ is now an exotic electron $`E^{}`$. Following the same steps as for model A we end up with the following multiplet structure:

$$\begin{array}{cccc}& & & \\ & & & \\ \chi _L=\left(\begin{array}{c}u\\ d\\ U\end{array}\right)_L& d_L^c& u_L^c& U_L^c\\ & & & \\ (3,3,\frac{1}{3})& (3^{},1,\frac{1}{3})& (3^{},1,\frac{2}{3})& (3^{},1,\frac{2}{3})\\ & & & \end{array}$$

$$\begin{array}{cccccc}& & & & & \\ & & & & & \\ \psi _L=\left(\begin{array}{c}e^{}\\ \nu _e\\ E_1^{}\end{array}\right)_L& \psi _{1L}=\left(\begin{array}{c}N_1^0\\ E_2^+\\ \nu _e^c\end{array}\right)_L& \psi _{2L}=\left(\begin{array}{c}E_2^{}\\ N_2^0\\ E_3^{}\end{array}\right)_L& e_L^+& E_{1L}^+& E_{3L}^+\\ & & & & & \\ (1,3^{},\frac{2}{3})& (1,3^{},\frac{1}{3})& (1,3^{},\frac{2}{3})& (1,1,1)& (1,1,1)& (1,1,1)\\ & & & & & \end{array}$$

A simple check shows that this multiplet structure is also free of anomalies. A phenomenological analysis of this model has been started already in Ref., where it is shown that model B as presented here is a subgroup of $`SU\left(6\right)U\left(1\right)`$, an electroweak-strong unification group which has not been considered in the literature so far.

The gauge boson content of models A and B are equivalent, and they become the same just by replacing $`33^{}`$ in the irreducible representations of the fermion fields ($`bb`$ when the complex conjugate of the covariant derivative is taken).

#### 2.1.3 Other one-family Models

Following the same steps as for the two previous cases, we attempt to construct models where $`q_L`$ has electric charges $`2/3`$ or 1/3. Eq.(2) then implies that $`Q_l=1/3`$ and $`2/3`$ respectively which correspond to leptons with exotic electric charges. Not only fractionally charged free particles has not been detected at low energies, but the phenomenology of those models could become tremendously confusing with leptons with electric charges equal to the antiquarks.

In a similar way by asking for a model with $`Q_l=1`$ we will have, according to Eq.(2), that $`Q_q=4/3`$, a model with a quark with an exotic electric charge which we have excluded from the models discussed here (a model with quarks with exotic electric charges is presented in Ref. for example).

### 2.2 Family models

For these models each individual family possesses non-vanishing anomalies and the anomaly cancellation takes place between families and, for some models, only with a matching of the number of families with the number of quark colors, does the overall anomaly vanish. It is also a feature of this type of models that the third family is treated different to the other two, or either that all the three families are treated independently.

An algebraic manipulation of Eqs.(1) and (2) and the anomaly constraint equations, allows us to combine the fermion multiplets of the two models A and B to produce the following models (with the replacement $`33^{}`$ in model B in order to assure a unique covariant derivative):

#### 2.2.1 Model C

All the left-handed lepton generations belong to the representation $`(1,3,2/3)`$ of $`(SU\left(3\right)_c,SU\left(3\right)_L,U\left(1\right)_X)`$, that is:

$$\begin{array}{ccc}& & \\ & & \\ \psi _L^\alpha =\left(\begin{array}{c}\nu _\alpha \\ \alpha ^{}\\ E_\alpha ^{}\end{array}\right)_L& \alpha _L^+& E_{\alpha L}^+\\ & & \\ (1,3,2/3)& (1,1,1)& (1,1,1)\\ & & \end{array}$$

for $`\alpha =e,\mu ,\tau `$; while quarks transform as follows:

$$\begin{array}{cccc}& & & \\ & & & \\ \chi _L^a=\left(\begin{array}{c}d^a\\ u^a\\ U^a\end{array}\right)_L& u_L^{ac}& d_L^{ac}& U_L^{ac}\\ & & & \\ (3,3^{},1/3)& (3^{},1,\frac{2}{3})& (3^{},1,\frac{1}{3})& (3^{},1,\frac{2}{3})\\ & & & \end{array}$$

for $`a=1,2`$ two of the families. For the other family we have:

$$\begin{array}{cccc}& & & \\ & & & \\ \chi _{3L}=\left(\begin{array}{c}u_3\\ d_3\\ D\end{array}\right)_L& u_{3L}^c& d_{3L}^c& D_L^c\\ & & & \\ (3,3,0)& (3^{},1,\frac{2}{3})& (3^{},1,\frac{1}{3})& (3^{},1,\frac{1}{3})\\ & & & \end{array}.$$

The arithmetic shows that all the anomalies vanish for this fermion content. As far as we know the study of this model is relatively new in the literature; it was introduced for the first time in Ref..

At first glance this structure does not allow for neutrino masses; even though, a variant of this model, with the capability to explain the main features of the atmospheric and solar neutrino experimental results, has been presented in Ref..

#### 2.2.2 Model D

In a similar way we get the following multiplet structure:

$$\begin{array}{cccc}& & & \\ & & & \\ \chi _L^a=\left(\begin{array}{c}u_a\\ d_a\\ D_a\end{array}\right)_L& u_{aL}^c& d_{aL}^c& D_{aL}^c\\ & & & \\ (3,3,0)& (3^{},1,\frac{2}{3})& (3^{},1,\frac{1}{3})& (3^{},1,\frac{1}{3})\\ & & & \end{array}$$

for $`a=1,2`$, the quarks in two of the three families. For the quarks in the other family we have:

$$\begin{array}{cccc}& & & \\ & & & \\ \chi _L^3=\left(\begin{array}{c}d_3\\ u_3\\ U\end{array}\right)_L& u_{3L}^c& d_{3L}^c& U_L^c\\ & & & \\ (3,3^{},\frac{1}{3})& (3^{},1,\frac{2}{3})& (3^{},1,\frac{1}{3})& (3^{},1,\frac{2}{3})\\ & & & \end{array}.$$

The three lepton generations transform now as anti-triplets of $`SU\left(3\right)_L`$ as follows:

$$\begin{array}{cc}& \\ & \\ \psi _L^\alpha =\left(\begin{array}{c}\alpha ^{}\\ \nu _\alpha \\ N_\alpha ^0\end{array}\right)_L& \alpha _L^+\\ & \\ (1,3^{},\frac{1}{3})& (1,1,1)\\ & \end{array},$$

for $`\alpha =e,\mu ,\tau `$ the three families. This model has been largely studied in the literature (see Refs.). Again, this model has been used recently in connection with neutrino oscillations.

#### 2.2.3 Other models

Contrary to the one family models, we can now play the game of canceling the anomalies in several different ways.

We start by defining the following closed set of fermions (closed in the sense that they include the antiparticles of the charged particles):
$`S_1=[(\nu _\alpha ,\alpha ^{},E_\alpha ^{});\alpha ^+;E_\alpha ^+]`$ with quantum numbers $`(1,3,2/3);(1,1,1)`$ and $`(1,1,1)`$ respectively.
$`S_2=[(\alpha ^{},\nu _\alpha ,N_\alpha ^0);\alpha ^+]`$ with quantum numbers $`(1,3^{},1/3)`$ and $`(1,1,1)`$ respectively.
$`S_3=[(d,u,U);u^c;d^c;U^c]`$ with quantum numbers $`(3,3^{},1/3);(3^{},1,2/3);(3^{},1,1/3)`$ and $`(3^{},1,2/3)`$ respectively.
$`S_4=[(u,d,D);d^c;u^c;D^c]`$ with quantum numbers $`(3,3,0);(3^{},1,1/3);(3^{},1,2/3)`$ and $`(3^{},1,1/3)`$ respectively.
$`S_5=[(e^{},\nu _e,N_1^0);(E^{},N_2^0,N_3^0)`$; $`(N_4^0,E^+,e^+)]`$ with quantum numbers $`(1,3^{},1/3)`$; $`(1,3^{},1/3)`$ and $`(1,3^{},2/3)`$ respectively.
$`S_6=[(\nu _e,e^{},E^{});(E_2^+,N_1^0,N_2^0);(N_3^0,E_2^{},E_3^{});e^+,E_1^+;E_3^+]`$ with quantum numbers $`(1,3,2/3)`$; $`(1,3,1/3)`$; $`(1,3,2/3);(1,1,1);(1,1,1)`$ and $`(1,1,1)`$ respectively.

Now we calculate the four anomalies for each set of particles. The results are presented in Table I.

TABLE I. Anomalies for $`S_i`$.

| Anomalies | $`S_1`$ | $`S_2`$ | $`S_3`$ | $`S_4`$ | $`S_5`$ | $`S_6`$ |
| --- | --- | --- | --- | --- | --- | --- |
| $`\left[SU\left(3\right)_c\right]^2U\left(1\right)_X`$ | 0 | 0 | 0 | 0 | 0 | 0 |
| $`\left[SU\left(3\right)_L\right]^2U\left(1\right)_X`$ | $`2/3`$ | $`1/3`$ | 1 | 0 | 0 | $`1`$ |
| $`\left[grav\right]^2U\left(1\right)_X`$ | 0 | 0 | 0 | 0 | 0 | 0 |
| $`\left[U\left(1\right)_X\right]^3`$ | 10/9 | 8/9 | $`12/9`$ | $`6/9`$ | 6/9 | 12/9 |

Notice from Table I that model A is just $`\left(S_4+S_5\right)`$ and model B is $`\left(S_3+S_6\right)`$. Model C is represented by $`\left(3S_1+2S_3+S_4\right)`$ and model D by $`\left(3S_2+S_3+2S_4\right)`$, but what is most remarkable is that we can now construct new anomaly-free models for two, three, four and more families. For example two new three family models are:
Model E: $`S_1+S_2+S_3+S_4`$ plus Model A = $`\left(S_1+S_2+S_3+2S_4+S_5\right)`$
Model F $`S_1+S_2+S_3+S_4`$ plus Model B=$`\left(S_1+S_2+2S_3+S_4+S_6\right)`$.
A model for four families will be given for example by: $`2\left(S_1+S_2+S_3+S_4\right)`$, etc..

The main feature of models E and F above is that, contrary to models C and D, each one of the three families is treated in a different way. As far as we know, these two models have not been studied in the literature so far.

## 3 The scalar sector

Even though the representation content for the fermions may vary significantly from model to model, all $`SU\left(3\right)_LU\left(1\right)_X`$ models presented so far have the same gauge boson sector as it will be discussed in the following section. Now, to achieve an spontaneous breaking of the symmetry in the most economic way, using the chain

$$SU\left(3\right)_cSU\left(3\right)_LU\left(1\right)_XSU\left(3\right)_cSU\left(2\right)_LU\left(1\right)_YSU\left(3\right)_cU\left(1\right)_Q,$$

we need two complex higgs scalars $`\varphi _i(1,3^{},1/3)=(\varphi _i^{},\varphi _i^0,\varphi _i^{}_{}{}^{}0),i=1,2`$, with Vacuum Expectation Values (VEV) $`\varphi _1=(0,0,V)^T`$ and $`\varphi _2=(0,v/\sqrt{2},0)^T`$, with the hierarchy $`V>>v250`$ GeV the electroweak mass scale. Now, to break the symmetry and at the same time give masses to all the fermion fields is a model dependent analysis. So, let us outline in the following three sections the analysis for model A for example, for which a third higgs field $`\varphi _3(1,3^{},2/3)=(\varphi _3^0,\varphi _3^+,\varphi _3^{}_{}{}^{}+)`$ with VEV $`\varphi _3=(v^{}/\sqrt{2},0,0)^T`$ is also needed, where $`v^{}v`$.

The analysis for model B is done in Ref. and for model D in Ref..

## 4 The gauge boson sector

There are a total of 17 gauge bosons in the gauge group under consideration; they are: one gauge field $`B^\mu `$ associated with $`U\left(1\right)_X`$, the 8 gluon fields associated with $`SU\left(3\right)_c`$ which remain massless after breaking the symmetry, and another 8 associated with $`SU\left(3\right)_L`$ and that we write for convenience in the following way:

$$\frac{1}{2}\lambda _\alpha A_\alpha ^\mu =\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}D_1^\mu & W^{+\mu }& K^{+\mu }\\ W^\mu & D_2^\mu & K^{0\mu }\\ K^\mu & \overline{K}^{0\mu }& D_3^\mu \end{array}\right),$$

where $`D_1^\mu =A_3^\mu /\sqrt{2}+A_8^\mu /\sqrt{6},D_2^\mu =A_3^\mu /\sqrt{2}+A_8^\mu /\sqrt{6}`$, and $`D_3^\mu =2A_8^\mu /\sqrt{6}`$. $`\lambda _i,i=1,2,\mathrm{},8`$ are the eight Gell-Mann matrices normalized as mentioned in Section 2.1 This allows us to write now the charge operator as

$$Q=\frac{\lambda _3}{2}+\frac{\lambda _8}{2\sqrt{3}}+XI_3,$$

where $`I_3`$ is the $`3\times 3`$ unit matrix.

In model A, after breaking the symmetry with $`\varphi _i,i=1,2,3`$, and using for the covariant derivative for triplets $`D^\mu =^\mu i\frac{g}{2}\lambda _\alpha A_\alpha ^\mu ig^{}XB^\mu `$, we get the following mass terms for the charged gauge bosons in the electroweak sector: $`M_{W^\pm }^2=\frac{g^2}{4}\left(v^2+v^2\right),M_{K^\pm }^2=\frac{g^2}{4}\left(2V^2+v^2\right),M_{K^0\left(\overline{K}^0\right)}^2=\frac{g^2}{4}\left(2V^2+v^2\right)`$. For the neutral gauge bosons we get a mass term of the form:

$$M=V^2\left(\frac{g^{}B^\mu }{3}\frac{gA_8^\mu }{\sqrt{3}}\right)^2+\frac{v^2}{8}\left(\frac{2g^{}B^\mu }{3}gA_3^\mu +\frac{gA_8^\mu }{\sqrt{3}}\right)^2+\frac{v^2}{8}\left(gA_3^\mu \frac{4g^{}B^\mu }{3}+\frac{gA_8^\mu }{\sqrt{3}}\right)^2$$

By diagonalizing $`M`$ we get the physical neutral gauge bosons which are defined through the mixing angle $`\theta `$ and $`Z_\mu ,Z_\mu ^{}`$ by:

$`Z_1^\mu `$ $`=`$ $`Z_\mu \mathrm{cos}\theta +Z_\mu ^{}\mathrm{sin}\theta `$
$`Z_2^\mu `$ $`=`$ $`Z_\mu \mathrm{sin}\theta +Z_\mu ^{}\mathrm{cos}\theta `$
$`\mathrm{tan}\left(2\theta \right)`$ $`=`$ $`{\displaystyle \frac{\sqrt{12}C_W\left(1T_W^2/3\right)^{1/2}\left[v^2\left(1+T_W^2\right)v^2\left(1T_W^2\right)\right]}{3\left(1T_W^2/3\right)\left(v^2+v^2\right)C_W^2\left[8V^2+v^2\left(1T_W^2\right)^2+v^2\left(1+T_W^2\right)^2\right]}},`$
where the photon field $`A^\mu `$ and the fields $`Z_\mu `$ and $`Z_\mu ^{}`$ are given by

$`A^\mu `$ $`=`$ $`S_WA_3^\mu +C_W\left[{\displaystyle \frac{T_W}{\sqrt{3}}}A_8^\mu +\left(1T_W^2/3\right)^{1/2}B^\mu \right],`$
$`Z^\mu `$ $`=`$ $`C_WA_3^\mu S_W\left[{\displaystyle \frac{T_W}{\sqrt{3}}}A_8^\mu +\left(1T_W^2/3\right)^{1/2}B^\mu \right],`$
$`Z^\mu `$ $`=`$ $`\left(1T_W^2/3\right)^{1/2}A_8^\mu +{\displaystyle \frac{T_W}{\sqrt{3}}}B^\mu .`$ (4)
$`S_W`$ and $`C_W`$ are the sine and cosine of the electroweak mixing angle respectively ($`T_W=S_W/C_W`$) defined by $`S_W=\sqrt{3}g^{}/\sqrt{3g^2+4g^2}`$. Also we can identify the $`Y`$ hypercharge associated with the SM gauge boson as:

$$Y^\mu =\left[\frac{T_W}{\sqrt{3}}A_8^\mu +\left(1T_W^2/3\right)^{1/2}B^\mu \right].$$

In the limit $`\theta 0`$, $`M_Z=M_{W^\pm }/C_W`$, and $`Z_1^\mu =Z^\mu `$ is the gauge boson of the SM. This limit is obtained either by demanding $`V\mathrm{}`$ or $`v^2=v^2\left(C_W^2S_W^2\right)`$. In general $`\theta `$ may be different from zero although it takes a very small value, determined from phenomenology for each particular model.

The former definitions for $`A^\mu ,Z^\mu ,Z^{}_{}{}^{}\mu ,Y^\mu `$ and $`S_W`$ are the same for all the six models in Section 2. The value for $`\mathrm{tan}\left(2\theta \right)`$ and the expressions for the masses of the gauge bosons are model dependent.

### 4.1 Currents

The currents for fermions are different for each model and also they are different from those of the SM. As an example let us present the analysis for model A; a similar analysis for model B is presented in Ref. and for model D in Ref..

#### 4.1.1 Charged currents

The interactions among the charged vector fields with leptons for model A are

$`H^{CC}`$ $`=`$ $`{\displaystyle \frac{g}{\sqrt{2}}}[W_\mu ^+(\overline{u}_L\gamma ^\mu d_L\overline{\nu }_{eL}\gamma ^\mu e_L^{}\overline{N}_{2L}^0\gamma ^\mu E_L^{}\overline{E}_L^+\gamma ^\mu N_{4L}^0)`$ (5)
$`+K_\mu ^+\left(\overline{u}_L\gamma ^\mu D_L\overline{N}_{1L}^0\gamma ^\mu e_L^{}\overline{N}_{3L}^0\gamma ^\mu E_L^{}\overline{e}_L^+\gamma ^\mu N_{4L}^0\right)`$
$`+K_\mu ^0(\overline{d}_L\gamma ^\mu D_L\overline{N}_{1L}^0\gamma ^\mu \nu _{eL}\overline{N}_{3L}^0\gamma ^\mu N_{2L}^0\overline{e}_L^+\gamma ^\mu E_L^+)]+H.c.,`$
which implies that the interactions with $`K^\pm `$ and $`K^0\left(\overline{K}^0\right)`$ bosons violate the lepton number and the weak isospin. Notice also that the first two terms in the previous expression constitute the charged weak current of the SM as far as we identify $`W^\pm `$ as the $`SU\left(2\right)_L`$ charged left-handed weak bosons.

#### 4.1.2 Neutral currents

The neutral currents $`J_\mu \left(EM\right),J_\mu \left(Z\right)`$ and $`J_\mu \left(Z^{}\right)`$, associated with the Hamiltonian $`H^0=eA^\mu J_\mu \left(EM\right)+\frac{g}{C_W}Z^\mu J_\mu \left(Z\right)+\frac{g^{}}{\sqrt{3}}Z^\mu J_\mu \left(Z^{}\right)`$ are:

$`J_\mu \left(EM\right)`$ $`=`$ $`{\displaystyle \frac{2}{3}}\overline{u}\gamma _\mu u{\displaystyle \frac{1}{3}}\overline{d}\gamma _\mu d{\displaystyle \frac{1}{3}}\overline{D}\gamma _\mu D\overline{e}^{}\gamma _\mu e^{}\overline{E}^{}\gamma _\mu E^{}={\displaystyle \underset{f}{}}q_f\overline{f}\gamma _\mu f`$
$`J_\mu \left(Z\right)`$ $`=`$ $`J_{\mu ,L}\left(Z\right)S_W^2J_\mu \left(EM\right)`$
$`J_\mu \left(Z^{}\right)`$ $`=`$ $`T_WJ_\mu \left(EM\right)J_{\mu ,L}\left(Z^{}\right),`$ (6)
where $`e=gS_W=g^{}C_W\sqrt{1T_W^2/3}>0`$ is the electric charge, $`q_f`$ is the electric charge of the fermion $`f`$ in units of $`e`$, $`J_\mu \left(EM\right)`$ is the electromagnetic current (vectorlike as it should be), and the left-handed currents are

$`J_{\mu ,L}\left(Z\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\overline{u}_L\gamma _\mu u_L\overline{d}_L\gamma _\mu d_L+\overline{\nu }_{eL}\gamma _\mu \nu _{eL}\overline{e}_L^{}\gamma _\mu e_L^{}+\overline{N}_2^0\gamma _\mu N_2^0\overline{E}^{}\gamma _\mu E^{}\right)`$
$`=`$ $`{\displaystyle \underset{f}{}}T_{3f}\overline{f}_L\gamma _\mu f_L`$
$`J_{\mu ,L}\left(Z^{}\right)`$ $`=`$ $`S_{2W}^1\left(\overline{u}_L\gamma _\mu u_L\overline{e}_L^{}\gamma _\mu e_L^{}\overline{E}_L^{}\gamma _\mu E_L^{}\overline{N}_{4L}^0\gamma _\mu N_{4L}^0\right)`$ (7)
$`T_{2W}^1\left(\overline{d}_L\gamma _\mu d_L\overline{E}_L^+\gamma _\mu E_L^+\overline{\nu }_{eL}\gamma _\mu \nu _{eL}\overline{N}_{2L}^0\gamma _\mu N_{2L}^0\right)`$
$`T_W^1\left(\overline{D}_L\gamma _\mu D_L\overline{e}_L^+\gamma _\mu e_L^+\overline{N}_{1L}^0\gamma _\mu N_{1L}^0\overline{N}_{3L}^0\gamma _\mu N_{3L}^0\right)`$
$`=`$ $`{\displaystyle \underset{f}{}}T_{9f}\overline{f}_L\gamma _\mu f_L,`$
where $`S_{2W}=2S_WC_W,T_{2W}=S_{2W}/C_{2W},C_{2W}=C_W^2S_W^2,\overline{N}_2^0\gamma _\mu N_2^0=\overline{N}_{2L}^0\gamma _\mu N_{2L}^0+\overline{N}_{2R}^0\gamma _\mu N_{2R}^0=\overline{N}_{2L}^0\gamma _\mu N_{2L}^0\overline{N}_{2L}^{0c}\gamma _\mu N_{2L}^{0c}=\overline{N}_{2L}^0\gamma _\mu N_{2L}^0\overline{N}_{4L}^0\gamma _\mu N_{4L}^0`$, and similarly $`\overline{E}\gamma _\mu E=\overline{E}_L^{}\gamma _\mu E_L^{}\overline{E}_L^+\gamma _\mu E_L^+`$. In this way $`T_{3f}=Dg.(1/2,1/2,0)`$ is the third component of the weak isospin acting on the representation 3 of $`SU\left(3\right)_L`$ (the negative when acting on $`3^{}`$), and $`T_{9f}=Dg.(S_{2W}^1,T_{2W}^1,T_W^1)`$ is a convenient $`3\times 3`$ diagonal matrix acting on the representation 3 of $`SU\left(3\right)_L`$ (the negative when acting on $`3^{}`$). Notice that $`J_\mu \left(Z\right)`$ is just the generalization of the neutral current present in the SM, which allows us to identify $`Z_\mu `$ as the neutral gauge boson of the SM.

The couplings of the physical states $`Z_1^\mu `$ and $`Z_2^\mu `$ are then given by:

$`H^{NC}`$ $`=`$ $`{\displaystyle \frac{g}{2C_W}}{\displaystyle \underset{i=1}{\overset{2}{}}}Z_i^\mu {\displaystyle \underset{f}{}}\left\{\overline{f}\gamma _\mu \left[a_{iL}\left(f\right)\left(1\gamma _5\right)+a_{iR}\left(f\right)\left(1+\gamma _5\right)\right]f\right\}`$ (8)
$`=`$ $`{\displaystyle \frac{g}{2C_W}}{\displaystyle \underset{i=1}{\overset{2}{}}}Z_i^\mu {\displaystyle \underset{f}{}}\left\{\overline{f}\gamma _\mu \left[g\left(f\right)_{iV}g\left(f\right)_{iA}\gamma _5\right]f\right\},`$
where

$`a_{1L}\left(f\right)`$ $`=`$ $`\mathrm{cos}\theta \left(T_{3f}q_fS_W^2\right){\displaystyle \frac{g^{}\mathrm{sin}\theta C_W}{g\sqrt{3}}}\left(T_{9f}q_fT_W\right)`$
$`a_{1R}\left(f\right)`$ $`=`$ $`q_fS_W\left(\mathrm{cos}\theta S_W{\displaystyle \frac{g^{}\mathrm{sin}\theta }{g\sqrt{3}}}\right)`$
$`a_{2L}\left(f\right)`$ $`=`$ $`\mathrm{sin}\theta \left(T_{3f}q_fS_W^2\right){\displaystyle \frac{g^{}\mathrm{cos}\theta C_W}{g\sqrt{3}}}\left(T_{9f}q_fT_W\right)`$
$`a_{2R}\left(f\right)`$ $`=`$ $`q_fS_W\left(\mathrm{sin}\theta S_W+{\displaystyle \frac{g^{}\mathrm{cos}\theta }{g\sqrt{3}}}\right),`$ (9)
and

$`g\left(f\right)_{1V}`$ $`=`$ $`\mathrm{cos}\theta \left(T_{3f}2S_W^2q_f\right){\displaystyle \frac{g^{}\mathrm{sin}\theta }{g\sqrt{3}}}\left(T_{9f}C_W2q_fS_W\right)`$
$`g\left(f\right)_{2V}`$ $`=`$ $`\mathrm{sin}\theta \left(T_{3f}2S_W^2q_f\right){\displaystyle \frac{g^{}\mathrm{cos}\theta }{g\sqrt{3}}}\left(T_{9f}C_W2q_fS_W\right)`$
$`g\left(f\right)_{1A}`$ $`=`$ $`\mathrm{cos}\theta T_{3f}{\displaystyle \frac{g^{}\mathrm{sin}\theta }{g\sqrt{3}}}T_{9f}C_W`$
$`g\left(f\right)_{2A}`$ $`=`$ $`\mathrm{sin}\theta T_{3f}{\displaystyle \frac{g^{}\mathrm{cos}\theta }{g\sqrt{3}}}T_{9f}C_W,`$ (10)
to be compared with the SM values $`g\left(f\right)_{1V}^{SM}=T_{3f}2S_Wq_f`$ and $`g\left(f\right)_{1A}^{SM}=T_{3f}`$. The values of $`g_{iV},g_{iA};i=1,2`$ are listed in Tables II and III. As we can see, in the limit $`\theta =0`$ the couplings of $`Z_1^\mu `$ to the ordinary leptons and quarks are the same as in the SM. Because of this, we can test the new phenomenology beyond the SM.

TABLE II. The $`Z_1^\mu \overline{f}f`$ couplings.     f $`g_{1V}`$ $`g_{1A}`$    u $`\left(\frac{1}{2}\frac{4S_W^2}{3}\right)\left[\mathrm{cos}\theta \mathrm{sin}\theta /\left(4C_W^21\right)^{1/2}\right]`$ $`\frac{\mathrm{cos}\theta }{2}\mathrm{sin}\theta /\left[2\left(4C_w^21\right)^{1/2}\right]`$  d $`\mathrm{cos}\theta \left(\frac{1}{2}+\frac{2S_W^2}{3}\right)\frac{\mathrm{sin}\theta }{\left(4C_W^21\right)^{1/2}}\left(\frac{1}{2}\frac{S_W^2}{3}\right)`$ $`\frac{1}{2}\left\{\mathrm{cos}\theta +\mathrm{sin}\theta C_{2W}/\left[2\left(4C_W^21\right)^{1/2}\right]\right\}`$  D $`\frac{2S_W^2\mathrm{cos}\theta }{3}+\mathrm{sin}\theta \left(1\frac{5}{3}S_W^2\right)/\left(4C_W^21\right)^{1/2}`$ $`C_W^2\mathrm{sin}\theta /\left(4C_W^21\right)^{1/2}`$  $`e^{}`$ $`\mathrm{cos}\theta \left(\frac{1}{2}+2S_W^2\right)+\frac{3\mathrm{sin}\theta }{\left(4C_W^21\right)^{1/2}}\left(\frac{1}{2}S_W^2\right)`$ $`\frac{\mathrm{cos}\theta }{2}+\frac{\mathrm{sin}\theta }{\left(4C_W^21\right)^{1/2}}\left(\frac{1}{2}C_W^2\right)`$  $`E^{}`$ $`\mathrm{cos}\theta \left(1+2S_W^2\right)\frac{S_W^2\mathrm{sin}\theta }{\left(4C_W^21\right)^{1/2}}`$ $`C_W^2\mathrm{sin}\theta /\left(4C_W^21\right)^{1/2}`$  $`\nu _e,N_2^0`$ $`\frac{1}{2}\left[\mathrm{cos}\theta +\mathrm{sin}\theta \left(12S_W^2\right)/\left(4C_W^21\right)^{1/2}\right]`$ $`\frac{1}{2}(\mathrm{cos}\theta +\mathrm{sin}\theta (12S_W^2)/(4C_W^21)^{1/2}`$  $`N_1^0,N_3^0`$ $`C_W^2\mathrm{sin}\theta /\left(4C_W^21\right)^{1/2}`$ $`C_W^2\mathrm{sin}\theta /\left(4C_W^21\right)^{1/2}`$  $`N_4^0`$ $`\frac{1}{2}\left[\mathrm{cos}\theta \mathrm{sin}\theta /\left(4C_W^21\right)^{1/2}\right]`$ $`\frac{1}{2}\left[\mathrm{cos}\theta \mathrm{sin}\theta /\left(4C_W^21\right)^{1/2}\right]`$

TABLE III. The $`Z_2^\mu \overline{f}f`$ couplings.     f $`g_{2V}`$ $`g_{2A}`$    u $`\left(\frac{1}{2}\frac{4S_W^2}{3}\right)\left[\mathrm{sin}\theta \mathrm{cos}\theta /\left(4C_W^21\right)^{1/2}\right]`$ $`\frac{\mathrm{sin}\theta }{2}\mathrm{cos}\theta /\left[2\left(4C_w^21\right)^{1/2}\right]`$  d $`\mathrm{sin}\theta \left(\frac{1}{2}+\frac{2S_W^2}{3}\right)\frac{\mathrm{cos}\theta }{\left(4C_W^21\right)^{1/2}}\left(\frac{1}{2}\frac{S_W^2}{3}\right)`$ $`\frac{1}{2}\left\{\mathrm{sin}\theta +\mathrm{cos}\theta C_{2W}/\left[2\left(4C_W^21\right)^{1/2}\right]\right\}`$  D $`\frac{2S_W^2\mathrm{sin}\theta }{3}+\mathrm{cos}\theta \left(1\frac{5}{3}S_W^2\right)/\left(4C_W^21\right)^{1/2}`$ $`C_W^2\mathrm{cos}\theta /\left(4C_W^21\right)^{1/2}`$  $`e^{}`$ $`\mathrm{sin}\theta \left(\frac{1}{2}+2S_W^2\right)+\frac{3\mathrm{cos}\theta }{\left(4C_W^21\right)^{1/2}}\left(\frac{1}{2}S_W^2\right)`$ $`\frac{\mathrm{sin}\theta }{2}+\frac{\mathrm{cos}\theta }{\left(4C_W^21\right)^{1/2}}\left(\frac{1}{2}C_W^2\right)`$  $`E^{}`$ $`\mathrm{sin}\theta \left(1+2S_W^2\right)\frac{S_W^2\mathrm{cos}\theta }{\left(4C_W^21\right)^{1/2}}`$ $`C_W^2\mathrm{cos}\theta /\left(4C_W^21\right)^{1/2}`$  $`\nu _e,N_2^0`$ $`\frac{1}{2}\left[\mathrm{sin}\theta +\mathrm{cos}\theta \left(12S_W^2\right)/\left(4C_W^21\right)^{1/2}\right]`$ $`\frac{1}{2}(\mathrm{sin}\theta +\mathrm{cos}\theta (12S_W^2)/(4C_W^21)^{1/2}`$  $`N_1^0,N_3^0`$ $`C_W^2\mathrm{cos}\theta /\left(4C_W^21\right)^{1/2}`$ $`C_W^2\mathrm{cos}\theta /\left(4C_W^21\right)^{1/2}`$  $`N_4^0`$ $`\frac{1}{2}\left[\mathrm{sin}\theta +\mathrm{cos}\theta /\left(4C_W^21\right)^{1/2}\right]`$ $`\frac{1}{2}\left[\mathrm{sin}\theta +\mathrm{cos}\theta /\left(4C_W^21\right)^{1/2}\right]`$

## 5 Masses for fermions

Again this subject is model dependent. Just for the sake of completeness let us write the Yukawa lagrangian that the three higgs scalars in Section 3 produce for the fermion fields in model A:

$`_Y`$ $`=`$ $`_Y^Q+_Y^l`$
$`_Y^Q`$ $`=`$ $`\chi _L^TC\left(h_u\varphi _3u_L^c+h_D\varphi _1D_L^c+h_d\varphi _2d_L^c+h_{dD}\varphi _2D_L^c+h_{Dd}\varphi _1d_L^c\right)+h.c.`$ (11)
$`_Y^l`$ $`=`$ $`ϵ_{abc}\left[\psi _L^aC\left(h_1\psi _{1L}^b\varphi _3^c+h_2\psi _{2L}^b\varphi _1^c+h_3\psi _{2L}^b\varphi _2^c\right)+\psi _{1L}^aC\left(h_4\psi _{2L}^b\varphi _1^c+h_5\psi _{2L}^b\varphi _2^c\right)\right]`$ (12)
$`+H.c.,`$
where $`h_\eta ,\eta =u,d,D,dD,Dd,1,2,3,4,5`$ are Yukawa couplings of order one, $`a,b,c`$ are $`SU\left(3\right)_L`$ tensor indices and $`C`$ is the charge conjugation operator.

Using the VEV as in section 3 and assuming that we are referring to the third family, we see that $`m_t=h_uv^{}/\sqrt{2}`$, $`m_Dh_DV`$ but it mixes with the $`b`$ quark producing a kind of see-saw mechanism that implies $`m_b<<m_t`$. Also for leptons we have $`m_Eh_4V`$ but again it mixes with the $`\tau `$ lepton producing also a kind of see saw mechanism which implies that $`m_\tau m_b<<m_t`$. The neutral sector is more complicated; the analysis of the $`5\times 5`$ mass matrix gives: first two eigenvalues $`\pm h_1v^{}/\sqrt{2}`$ which correspond to a Dirac neutrino with a mass of the order of the electroweak mass scale; other two are $`\pm V+\eta `$, where $`\eta `$ is a small see-saw quotient, which correspond to a very massive pseudo-Dirac neutrino, and finally a tiny mass Majorana neutrino.

So the higgs fields and VEV used break the symmetry in the appropriate way, and produce a realistic pattern of masses for the fermion fields (at least for one of the families).

## 6 Conclusions

In this paper we have studied the theory of $`SU\left(3\right)cSU\left(3\right)_LU\left(1\right)_X`$ in detail. By restricting the fermion field representations to particles without exotic electric charges we end up with six different models, two one family models and four models for three families. The two one family models are sketched in the papers by K.T. Mahanthappa and P.K. Mohapatra in Ref., but enough attention was not paid to the anomaly cancellation constraints in their analysis. The four three family models are relatively new in the literature, with two of them (models E and F) introduced here for the first time, as far as we know.

If we allow for particles with exotic electric charges in our analysis, we end up with an infinite number of models, where the model in Refs. is just one of them (probably the most elegant one!).

The low energy predictions of the six models discussed here are not the same. All of them have in common a new neutral current which mixes with the SM neutral current which is also included as part of each model. When the mixing angle between the two neutral currents is zero ($`\mathrm{sin}\theta =0)`$, exact agreement with the SM predictions is achieved, but the use of experimental results from LEP, SLAC and atomic parity violation bound the mixing angle to values which are model dependent. For partial analysis see for example Refs. and .

Detailed analysis in each model of flavor changing neutral currents, GIM mechanism, mass scales of the new gauge bosons, mass spectrum for the neutral spin 1/2 particles etc., are model dependent and they will be presented elsewhere.

Finally let us mention that the most remarkable result of our analysis is the existence of models E and F, where the three families are treated different. In these models it should be simple to implement the horizontal survival hypothesis, that is, to provide masses at tree level only for the particles in the third family, as done for example in the previous section, with the known particles in the first and second families getting masses as radiative corrections.

## 7 Acknowledgments

This work was partially supported by BID and Colciencias in Colombia. We thank C. García Canal for a critical reading of the original manuscript. WAP acknowledges warm hospitality from the Theoretical Physics Laboratory at the Universidad de la Plata in Argentina, where this work was completed.


format: text

text:
# Thermodynamics of hot dense H-plasmas: Path integral Monte Carlo simulations and analytical approximations

## I Introduction

The thermodynamics of strongly correlated Fermi systems at high pressures are of growing importance in many fields, including shock and laser plasmas, astrophysics, solids and nuclear matter, see Refs.  for an overview. In particular the thermodynamical properties of hot dense plasmas under high pressure are of importace for the description of plasmas relevant for laser fusion . Further among the phenomena of current interest are Fermi liquids, metallic hydrogen , plasma phase transition, e.g.  and references therein, bound states etc., which occur in situations where both Coulomb and quantum effects are relevant. There has been significant progress in recent years to study these systems analytically and numerically, see e.g. . Due to the enormeous difficulties to develop analytical descriptions for hydrogen plasmas with strong coupling, e.g. , there is still an urgent need to test the analytical theory by an independent numerical approach.

An approach which is particularly well suited to describe thermodynamic properties in the region of high pressure, characterized by strong coupling and strong degeneracy, is the path integral quantum Monte Carlo (PIMC) method. There has been remarkable recent progress in applying these techniques to Fermi systems, for an overview see e.g. Refs. . However, these simulations are essentially hampered by the fermion sign problem. To overcome this difficulty, several strategies have been developed to simulate macroscopic Coulomb systems : the first is the restricted PIMC concept where additional assumptions on the density operator $`\widehat{\rho }`$ are introduced which reduce the sum over permutations to even (positive) contributions only. This requires knowledge of the nodes of the density matrix which is available only in a few special cases, e.g. . However, for interacting macroscopic systems, these nodes are known only approximately, e.g. , and the accuracy of the results is difficult to assess from within this scheme. An alternative are direct fermionic PIMC simulations which have occasionally been attempted by various groups  but which were not sufficiently precise and efficient for practical purposes. Recently, three of us have proposed a new path integral representation for the N-particle density operator  which allows for direct fermionic path integral Monte Carlo simulations of dense plasmas in a broad range of densities and temperatures. Using this concept we computed the pressure and energy of a degenerate strongly coupled hydrogen plasma  and the pair distribution functions in the region of partial ionization and dissociation . This scheme is rather efficient when the number of time slices (beads) in the path integral is less or equal 50 and was found to work well for temperatures $`k_BT0.1Ry`$. In this paper we derive further improved formulas for the pressure and energy and give, for the first time, a detailed derivation of all main results and rigorously justify the use of the effective quantum pair potential (Kelbg potential) in direct PIMC simulations. Further, in the present work this method will be applied to high-pressure plasmas ($`p10^110^2`$Mbar) in such temperature regions were considerable deviations from the ideal behavior are observed.

One difficulty of PIMC simulations is that reliable error estimates are often not available, in particular for strongly coupled degenerate systems. Moreover, in this region no reliable data from other theories such as density functional theory or quantum statistics, e.g. , are available which would allow for an unambiguous test. Furthermore, results from classical molecular dynamics simulations exist, but they apply only to fully ionized and weakly degenerate plasmas, e.g. , which is outside the range of interest for this work. Also, new quantum molecular dynamics approaches are being developed, e.g. , but they are only beginning to produce accurate results.

Therefore, it is of high interest to perform quantitative comparisons of independent simulations, such as restricted and direct fermionic PIMC, and to develop improved analytical approximations, which is the aim of this paper. We compare recent results of Militzer et al.  for pressure and energy isochors ($`n2.510^{23}\mathrm{cm}^3`$) of dense hydrogen to our own direct PIMC results. This is a non-trivial comparison since the two approaches employ independent sets of approximations. Nevertheless, we find very good agreement for temperatures ranging from $`10^6K`$ to as low as $`50,000K`$. This is remarkable since there the coupling and degeneracy parameters reach rather large values, $`\mathrm{\Gamma }3`$ and $`n_e\mathrm{\Lambda }^310`$, and the plasma contains a substantial fraction of bound states.

Further, we use the new data to make a comparison with analytical estimates which are based on Padé approximations for strongly ionized plasmas. These formulae were constructed on the basis of the known analytical results for the limiting cases of low density  and high density . These Padé approximations are exact up to quadratic terms in the density and interpolate between the virial expansions and the high-density asymptotics . We find that the results for the internal energy and for the pressure agree well with the PIMC results in the region of the density temperature plane, where $`\mathrm{\Gamma }1.6`$ and $`n\mathrm{\Lambda }^35`$.

## II Path integral representation of thermodynamic quantities

We now our direct PIMC scheme. All thermodynamic properties of a two-component plasma are defined by the partition function $`Z`$ which, for the case of $`N_e`$ electrons and $`N_p`$ protons, is given by

$`Z(N_e,N_p,V,\beta )`$ $`=`$ $`{\displaystyle \frac{Q(N_e,N_p,\beta )}{N_e!N_p!}},`$ (1)
$`\text{with}Q(N_e,N_p,\beta )`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{V}{}}𝑑q𝑑r\rho (q,r,\sigma ;\beta ),`$ (2)
where $`\beta =1/k_BT`$. The exact density matrix is, for a quantum system, in general, not known but can be constructed using a path integral representation ,

$`{\displaystyle \underset{V}{}}𝑑R^{(0)}{\displaystyle \underset{\sigma }{}}\rho (R^{(0)},\sigma ;\beta )`$ $`=`$ $`{\displaystyle \underset{V}{}}𝑑R^{(0)}\mathrm{}𝑑R^{(n)}\rho ^{(1)}\rho ^{(2)}\mathrm{}\rho ^{(n)}`$ (3)
$`\times `$ $`{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{P}{}}(\pm 1)^{\kappa _P}𝒮(\sigma ,\widehat{P}\sigma ^{})\widehat{P}\rho ^{(n+1)},`$ (4)
where $`\rho ^{(i)}\rho (R^{(i1)},R^{(i)};\mathrm{\Delta }\beta )R^{(i1)}|e^{\mathrm{\Delta }\beta \widehat{H}}|R^{(i)}`$, whereas $`\mathrm{\Delta }\beta \beta /(n+1)`$ and $`\mathrm{\Delta }\lambda _a^2=2\pi \mathrm{}^2\mathrm{\Delta }\beta /m_a`$, $`a=p,e`$. $`\widehat{H}`$ is the Hamilton operator, $`\widehat{H}=\widehat{K}+\widehat{U}_c`$, containing kinetic and potential energy contributions, $`\widehat{K}`$ and $`\widehat{U}_c`$, respectively, with $`\widehat{U}_c=\widehat{U}_c^p+\widehat{U}_c^e+\widehat{U}_c^{ep}`$ being the sum of the Coulomb potentials between protons (p), electrons (e) and electrons and protons (ep)\]. Further, $`R^{(i)}=(q^{(i)},r^{(i)})(R_p^{(i)},R_e^{(i)})`$, for $`i=1,\mathrm{}n+1`$, $`R^{(0)}(q,r)(R_p^{(0)},R_e^{(0)})`$, and $`R^{(n+1)}R^{(0)}`$ and $`\sigma ^{}=\sigma `$. This means, the particles are represented by fermionic loops with the coordinates (beads) $`[R][R^{(0)};R^{(1)};\mathrm{};R^{(n)};R^{(n+1)}]`$, where $`q`$ and $`r`$ denote the electron and proton coordinates, respectively. The spin gives rise to the spin part of the density matrix $`𝒮`$, whereas exchange effects are accounted for by the permutation operator $`\widehat{P}`$, which acts on the electron coordinates and spin projections, and the sum over the permutations with parity $`\kappa _P`$. In the fermionic case (minus sign), the sum contains $`N_e!/2`$ positive and negative terms leading to the notorious sign problem. Due to the large mass difference of electrons and ions, the exchange of the latter is not included. The matrix elements $`\rho ^{(i)}`$ can be rewritten identically as

$`R^{(i1)}|e^{\mathrm{\Delta }\beta \widehat{H}}|R^{(i)}=`$ (5)
$`{\displaystyle 𝑑\stackrel{~}{p}^{(i)}𝑑\overline{p}^{(i)}R^{(i1)}|e^{\mathrm{\Delta }\beta \widehat{U}_c}|\stackrel{~}{p}^{(i)}\stackrel{~}{p}^{(i)}|e^{\mathrm{\Delta }\beta \widehat{K}}|\overline{p}^{(i)}\overline{p}^{(i)}|e^{\frac{\mathrm{\Delta }\beta ^2}{2}[\widehat{K},\widehat{U}_c]}\mathrm{}|R^{(i)}}.`$ (6)
To compute thermodynamic functions, the logarithm of the partition function has to be differentiated with respect to thermodynamic variables. In particular, for the equation of state $`p`$ and internal energy $`E`$ follows,

$`\beta p`$ $`=`$ $`\mathrm{ln}Q/V=[\alpha /3V\mathrm{ln}Q/\alpha ]_{\alpha =1},`$ (7)
$`\beta E`$ $`=`$ $`\beta \mathrm{ln}Q/\beta ,`$ (8)
where $`\alpha `$ is a length scaling parameter $`\alpha =L/L_0`$. This means, in the path integral representation (4), each high-temperature density matrix has to be differentiated in turn. For example, the result for the energy will have the form

$`\beta E`$ $`=`$ $`{\displaystyle \frac{1}{Q}}{\displaystyle \underset{V}{}}𝑑R^{(0)}\mathrm{}𝑑R^{(n)}`$ (9)
$`\times `$ $`{\displaystyle \underset{k=1}{\overset{n+1}{}}}\rho ^{(1)}\mathrm{}\rho ^{(k1)}\left[\beta {\displaystyle \frac{\rho ^{(k)}}{\beta }}\right]\rho ^{(k+1)}\mathrm{}\rho ^{(n)}{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{P}{}}(\pm 1)^{\kappa _P}𝒮(\sigma ,\widehat{P}\sigma ^{})\widehat{P}\rho ^{(n+1)},`$ (10)
and, analogously for other thermodynamic functions.

There are two different approaches to evaluate this expression. One is to first choose an approximation for the high-temperature density matrices $`\rho ^{(i)}`$ and then to perform the differentiation. The other way is to first differentiate the operator expression for $`\rho ^{(k)}`$ and use an approximation for the matrix elements only in the final result. As we checked, the second method is more accurate and will be used in the following.

To evaluate the derivatives in Eq. (10), it is convenient to indroduce dimensionless integration variables $`\eta ^{(k)}=(\eta _p^{(k)},\eta _e^{(k)})`$, where $`\eta _a^{(k)}=\kappa _a(R_a^{(k)}R_a^{(k1)})`$ for $`k=1,\mathrm{},n`$ and $`a=p,e`$, and $`\kappa _a^2m_ak_BT/(2\pi \mathrm{}^2)=1/\lambda _a^2`$, . This has the advantage that now the differentiation of the density matrix affects only the interaction terms. Indeed, one can show that

$`\beta {\displaystyle \frac{\rho ^{(k)}}{\beta }}`$ $`=`$ $`\beta {\displaystyle \frac{\mathrm{\Delta }\beta U_c(X^{(k1)})}{\beta }}\rho ^{(k)}+\beta \stackrel{~}{\rho }_\beta ^{(k)},`$ (11)
where $`X^{(0)}(\kappa _pR_p^{(0)},\kappa _eR_e^{(0)})`$, $`X^{(k)}(X_p^{(k)},X_e^{(k)})`$, with $`X_a^{(k)}=\kappa _aR_a^{(0)}+_{l=1}^k\eta _a^{(l)}`$, and $`k`$ runs from $`1`$ to $`n`$. Further, $`X^{(n+1)}(\kappa _pR_p^{(n+1)},\kappa _eR_e^{(n+1)})=X^{(0)}`$, and we denoted

$`\stackrel{~}{\rho }_\beta ^{(k)}={\displaystyle 𝑑p^{(k)}X^{(k1)}|e^{\mathrm{\Delta }\beta \widehat{U}_c}|p^{(k)}e^{\frac{p^{(k)}|p^{(k)}}{4\pi (n+1)}}p^{(k)}|\frac{}{\beta }e^{\frac{(\mathrm{\Delta }\beta )^2}{2}[\widehat{K},\widehat{U}_c]}\mathrm{}|X^{(k)}},`$ (12)
where $`p_a^{(k)}=\stackrel{~}{p}_a^{(k)}/(\kappa _a\mathrm{})`$, $`p(k)(p_p^{(k)},p_e^{(k)})`$ and use has been made of Eq. (6). For $`k=n+1`$, we have

$`\beta {\displaystyle \frac{}{\beta }}{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{P}{}}(\pm 1)^{\kappa _P}𝒮(\sigma ,\widehat{P}\sigma ^{})\widehat{P}\rho ^{(n+1)}={\displaystyle \underset{\sigma }{}}{\displaystyle \underset{P}{}}(\pm 1)^{\kappa _P}𝒮(\sigma ,\widehat{P}\sigma ^{})\times `$ (13)
$`\times `$ $`\left\{\beta {\displaystyle \frac{\mathrm{\Delta }\beta U_c(X^{(n)})}{\beta }}\widehat{P}\rho ^{(n+1)}+\widehat{P}\left[\beta \stackrel{~}{\rho }_\beta ^{(n+1)}\right]\right\}.`$ (14)
Further, $`U_c(X^{(k1)})U_c^{(1)}(X^{(k1)})+U_c^{(2)}(X^{(k1)})`$, with $`U_c^{(1)}`$ and $`U_c^{(2)}`$ denoting the interaction between identical and different particle species, respectively, $`U_c^{(1)}(X)=U_c^e(X)+U_c^p(X)`$ and $`U_c^{(2)}(X)=U_c^{ep}(X)`$.

Using these results and Eq. (10), we obtain for the energy

$`\beta E`$ $`=`$ $`{\displaystyle \frac{3}{2}}(N_e+N_p){\displaystyle \frac{1}{Q}}{\displaystyle \frac{1}{\lambda _p^{3N_p}\lambda _e^{3N_e}}}{\displaystyle \underset{V}{}}𝑑R^{(0)}𝑑\eta ^{(1)}\mathrm{}𝑑\eta ^{(n)}{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{P}{}}(\pm 1)^{\kappa _P}𝒮(\sigma ,\widehat{P}\sigma ^{})`$ (15)
$`\times `$ $`\{{\displaystyle \underset{k=1}{\overset{n+1}{}}}\rho ^{(1)}\mathrm{}\rho ^{(k1)}[\beta {\displaystyle \frac{\mathrm{\Delta }\beta U_c^{(1)}(X^{(k1)})}{\beta }}\beta {\displaystyle \frac{\mathrm{\Delta }\beta U_c^{(2)}(X^{(k1)})}{\beta }}+\beta \stackrel{~}{\rho }_\beta ^{(n+1)}]`$ (16)
$`\times `$ $`\rho ^{(k)}\mathrm{}\rho ^{(n)}\widehat{P}\rho ^{(n+1)}\left\}\right|_{X^{(n+1)}=X^{(0)},\sigma ^{}=\sigma }.`$ (17)
This way, the derivative of the density matrix has been calculated, and we turn to the next point - to find approximations for the high-temperature density matrix.

## III High-temperature asymptotics of the density matrix in the Path integral approach. Kelbg potential

In this section we derive an approximation for the high-tempature density matrix which is suitable for direct PIMC simulations. Further, we demonstrate that the proper choice of the effective quantum pair potential is given by the Kelbg potential. Following Refs. , we derive a modified representation for the density matrix. The mains steps are:

1. The N-particle density matrix is expanded in terms of 2-particle, 3-particle etc. contributions from which only the first, $`\rho _{ab}`$, is retained ;
2. In the high-temperature limit, $`\rho _{ab}`$ factorizes into a kinetic ($`\rho _0`$) and an interaction term ($`\rho _U^{ab}`$), $`\rho _{ab}\rho _0\rho _U^{ab}`$, because it can be shown that

$`e^{\frac{(\mathrm{\Delta }\beta )^2}{2}[\widehat{K},\widehat{U}_c]}=\widehat{I}+O\left({\displaystyle \frac{1}{(n+1)^2}}\right),`$ (18)
where $`\widehat{I}`$ is the unity operator. In this way we get the following representation for the two-particle density matrix

$$\rho _{ab}=\left(\frac{(m_am_b)^{3/2}}{(2\pi \mathrm{}\beta )^3}\right)\mathrm{exp}[\frac{m_a}{2\mathrm{}^2\beta }(𝐫_a𝐫_a^{})^2]\mathrm{exp}[\frac{m_b}{2\mathrm{}^2\beta }(𝐫_b𝐫_b^{})^2]\mathrm{exp}[\beta \mathrm{\Phi }_{ab}]$$
(19)
where $`\mathrm{\Phi }_{ab}(𝐫_a,𝐫_a^{},𝐫_b,𝐫_b^{})`$ is the off-diagonal two-particle effective potential.
3. In the following, the off-diagonal matrix elements of the effective binary potentials will be approximated by the diagonal ones by taking the Kelbg potential at the center coordinate, $`\mathrm{\Phi }^{ab}(r,r^{};\mathrm{\Delta }\beta )\mathrm{\Phi }^{ab}(\frac{r+r^{}}{2};\mathrm{\Delta }\beta )`$.;
4. For the plasma parameter region of interest, the protons can be treated classically, and $`\mathrm{\Phi }^{ii}`$ may be approximated by the Coulomb potential.

We will now comment on these steps in some more detail. We calculated the effective potential by solving a Bloch equation by first order perturbation theory. The method has been described in detail in . This procedure defines an effective off-diagonal quantum pair potential for Coulomb systems, which depends on the inter-particle distances $`𝐫_{ab},𝐫_{ab}^{}`$. As a result of first-order perturbation theory we get explicitely

$`\mathrm{\Phi }^{ab}(𝐫_{ab},𝐫_{ab}^{},\mathrm{\Delta }\beta )e_ae_b{\displaystyle _0^1}{\displaystyle \frac{d\alpha }{d_{ab}(\alpha )}}\mathrm{erf}\left({\displaystyle \frac{d_{ab}(\alpha )}{2\lambda _{ab}\sqrt{\alpha (1\alpha )}}}\right),`$ (20)
where $`d_{ab}(\alpha )=|\alpha 𝐫_{ab}+(1\alpha )𝐫_{ab}^{}|`$, $`\mathrm{erf}(x)`$ is the error function $`\mathrm{erf}(x)=\frac{2}{\sqrt{\pi }}_0^x𝑑te^{t^2}`$, and $`\lambda _{ab}^2=\frac{\mathrm{}^2\mathrm{\Delta }\beta }{2\mu _{ab}}`$ with $`\mu _{ab}^1=m_a^1+m_b^1`$. It is interesting to note, that a simple approximation of the complicated integral over $`\alpha `$ by the length of the interval multiplied with the integrand in the center (Mittelwertsatz) leads us to the so-called KTR-potential due to Klakow, Toepffer and Reinhard which (in the diagonal approximation) is often used in quasi-classical MD simulations

$`\mathrm{\Phi }^{ab}(𝐫_{ab},𝐫_{ab}^{},\mathrm{\Delta }\beta ){\displaystyle \frac{e_ae_b}{d_{ab}(1/2)}}\mathrm{erf}\left({\displaystyle \frac{d_{ab}(1/2)}{\lambda _{ab}}}\right),`$ (21)
In our direct PIMC calculations we used the full expression for the interaction potential, keeping the $`\alpha `$-integration but, in order to save computer time, we approximated the two-particle interaction potential by its diagonal elements. The diagonal element ($`𝐫_{ab}^{}=𝐫_{ab}`$) of $`\mathrm{\Phi }^{ab}`$ is just the familiar Kelbg potential, given by (we will use the same notation for the potential)

$`\mathrm{\Phi }^{ab}(|𝐫_{ab}|,\mathrm{\Delta }\beta )\mathrm{\Phi }^{ab}(𝐫_{ab},𝐫_{ab},\mathrm{\Delta }\beta )={\displaystyle \frac{e_ae_b}{\lambda _{ab}x_{ab}}}[1e^{x_{ab}^2}+\sqrt{\{}\pi \}x_{ab}(1\mathrm{erf}(x_{ab}))],`$ (22)
where $`x_{ab}=|𝐫_{ab}|/\lambda _{ab}`$, and we underline that the Kelbg potential is finite at zero distance. The error of the above approximations, for each of the high-temperature factors on the r.h.s. of Eq. (4), is of the order $`1/(n+1)^2`$.

With these approximations, we obtain the result $`\rho ^{(k)}=\rho _0^{(k)}\rho _U^{(k)}+O[(1/n+1)^2]`$, where $`\rho _0^{(k)}`$ is the kinetic density matrix, while $`\rho _U^{(k)}=e^{\mathrm{\Delta }\beta U(X^{(k1)})}\delta (X^{(k1)}X^{(k)})`$, where $`U`$ denotes the following sum of Coulomb and Kelbg potentials, $`U(X^{(k)})=U_c^p(X_p^{(k)})+U^e(X_e^{(k)})+U^{ep}(X_p^{(k)},X_e^{(k)})`$. Notice that special care has to be taken in performing the derivatives with respect to $`\beta `$ of the Coulomb potentials which appear in Eq. (17). Indeed, products $`\rho ^{(1)}\mathrm{}\rho ^{(n)}\widehat{P}\rho ^{(n+1)}\beta \frac{\mathrm{\Delta }\beta U_c(X^{(k1)})}{\beta }`$ have a singularity at zero interparticle distance which is integrable but leads to difficulties in the simulations. Due to the Kelbg potential, for the e-e and p-p interaction, this singularity is weakenend, but it is enhanced for the e-p interaction. In order to assure efficient simulations we, therefore, further transform the e-p contribution in the following way:

$`{\displaystyle _0^1}𝑑\alpha {\displaystyle 𝑑R^{(k1)}R^{(k2)}|e^{\mathrm{\Delta }\beta \alpha \widehat{K}}|R^{(k1)}\left[\beta \frac{}{\beta }\left(\mathrm{\Delta }\beta U_c^{(2)}(R^{(k1)})\right)\right]}`$ (23)
$`\times R^{(k1)}|e^{\mathrm{\Delta }\beta (1\alpha )\widehat{K}}|R^{(k)}`$ (24)
$``$ $`R^{(k1)}|e^{\mathrm{\Delta }\beta \widehat{K}}|R^{(k)}\left[\beta {\displaystyle \frac{}{\beta }}\left(\mathrm{\Delta }\beta U^{(2)}(R^{(k1)})\right)\right]+O\left[(1/n+1)^2\right].`$ (25)
This means, within the standard error of our approximation $`O\left[(\mathrm{\Delta }\beta )^2\right]`$, we have replaced the e-p Coulomb potential $`U_c^{(2)}`$ by the corresponding Kelbg potential $`U^{(2)}`$, which is much better suited for MC simulations.

Thus, using $`\lambda _p\lambda _e`$, we finally obtain for the energy:

$`\beta E={\displaystyle \frac{3}{2}}(N_e+N_p)+{\displaystyle \frac{1}{Q}}{\displaystyle \frac{1}{\lambda _p^{3N_p}\mathrm{\Delta }\lambda _e^{3N_e}}}{\displaystyle \underset{s=0}{\overset{N_e}{}}}{\displaystyle }dqdrd\xi \rho _s(q,[r],\beta )\times `$ (26)
$`\{{\displaystyle \underset{p<t}{\overset{N_p}{}}}{\displaystyle \frac{\beta e^2}{|q_{pt}|}}+{\displaystyle \underset{l=0}{\overset{n}{}}}[{\displaystyle \underset{p<t}{\overset{N_e}{}}}{\displaystyle \frac{\mathrm{\Delta }\beta e^2}{|r_{pt}^l|}}+{\displaystyle \underset{p=1}{\overset{N_p}{}}}{\displaystyle \underset{t=1}{\overset{N_e}{}}}\mathrm{\Psi }_l^{ep}]`$ (27)
$`+{\displaystyle \underset{l=1}{\overset{n}{}}}\left[{\displaystyle \underset{p<t}{\overset{N_e}{}}}C_{pt}^l{\displaystyle \frac{\mathrm{\Delta }\beta e^2}{|r_{pt}^l|^2}}+{\displaystyle \underset{p=1}{\overset{N_p}{}}}{\displaystyle \underset{t=1}{\overset{N_e}{}}}D_{pt}^l{\displaystyle \frac{\mathrm{\Delta }\beta \mathrm{\Phi }^{ep}}{|x_{pt}^l|}}\right]`$ (28)
$`{\displaystyle \frac{1}{\mathrm{det}|\psi _{ab}^{n,1}|_s}}{\displaystyle \frac{\mathrm{det}|\psi _{ab}^{n,1}|_s}{\beta }}\},`$ (29)
$`\mathrm{with}C_{pt}^l={\displaystyle \frac{r_{pt}^l|y_{pt}^l}{2|r_{pt}^l|}},D_{pt}^l={\displaystyle \frac{x_{pt}^l|y_p^l}{2|x_{pt}^l|}},`$ (30)
and $`\mathrm{\Psi }_l^{ep}\mathrm{\Delta }\beta [\beta ^{}\mathrm{\Phi }^{ep}(|x_{pt}^l|,\beta ^{})]/\beta ^{}|_{\beta ^{}=\mathrm{\Delta }\beta }`$ contains the electron-proton Kelbg potential $`\mathrm{\Phi }^{ep}`$. Here, $`\mathrm{}|\mathrm{}`$ denotes the scalar product, and $`q_{pt}`$, $`r_{pt}`$ and $`x_{pt}`$ are differences of two coordinate vectors: $`q_{pt}q_pq_t`$, $`r_{pt}r_pr_t`$, $`x_{pt}r_pq_t`$, $`r_{pt}^l=r_{pt}+y_{pt}^l`$, $`x_{pt}^lx_{pt}+y_p^l`$ and $`y_{pt}^ly_p^ly_t^l`$, with $`y_a^n=\mathrm{\Delta }\lambda _e_{k=1}^n\xi _a^{(k)}`$. Here we introduced dimensionless distances between neighboring vertices on the loop, $`\xi ^{(1)},\mathrm{}\xi ^{(n)}`$, thus, explicitly, $`[r][r;y_e^{(1)};y_e^{(2)};\mathrm{}].`$ Further, the density matrix $`\rho _s`$ in Eq. (30) is given by

$`\rho _s(q,[r],\beta )=C_{N_e}^se^{\beta U(q,[r],\beta )}{\displaystyle \underset{l=1}{\overset{n}{}}}{\displaystyle \underset{p=1}{\overset{N_e}{}}}\varphi _{pp}^l\mathrm{det}|\psi _{ab}^{n,1}|_s,`$ (31)
where $`U(q,[r],\beta )=U_c^p(q)+\{U^e([r],\mathrm{\Delta }\beta )+U^{ep}(q,[r],\mathrm{\Delta }\beta )\}/(n+1)`$ and $`\varphi _{pp}^l\mathrm{exp}[\pi |\xi _p^{(l)}|^2]`$. We underline that the density matrix (31) does not contain an explicit sum over the permutations and thus no sum of terms with alternating sign. Instead, the whole exchange problem is contained in a single exchange matrix given by

$`\psi _{ab}^{n,1}_se^{\frac{\pi }{\mathrm{\Delta }\lambda _e^2}\left|(r_ar_b)+y_a^n\right|^2}_s.`$ (32)
As a result of the spin summation, the matrix carries a subscript $`s`$ denoting the number of electrons having the same spin projection. For more details, we refer to Refs. .

In similar way, we obtain the result for the equation of state,

$`{\displaystyle \frac{\beta pV}{N_e+N_p}}=1+{\displaystyle \frac{1}{N_e+N_p}}{\displaystyle \frac{(3Q)^1}{\lambda _p^{3N_p}\mathrm{\Delta }\lambda _e^{3N_e}}}{\displaystyle \underset{s=0}{\overset{N_e}{}}}{\displaystyle }dqdrd\xi \rho _s(q,[r],\beta )\times `$ (33)
$`\{{\displaystyle \underset{p<t}{\overset{N_p}{}}}{\displaystyle \frac{\beta e^2}{|q_{pt}|}}+{\displaystyle \underset{p<t}{\overset{N_e}{}}}{\displaystyle \frac{\mathrm{\Delta }\beta e^2}{|r_{pt}|}}{\displaystyle \underset{p=1}{\overset{N_p}{}}}{\displaystyle \underset{t=1}{\overset{N_e}{}}}|x_{pt}|{\displaystyle \frac{\mathrm{\Delta }\beta \mathrm{\Phi }^{ep}}{|x_{pt}|}}`$ (34)
$`+{\displaystyle \underset{l=1}{\overset{n}{}}}\left[{\displaystyle \underset{p<t}{\overset{N_e}{}}}A_{pt}^l{\displaystyle \frac{\mathrm{\Delta }\beta e^2}{|r_{pt}^l|^2}}{\displaystyle \underset{p=1}{\overset{N_p}{}}}{\displaystyle \underset{t=1}{\overset{N_e}{}}}B_{pt}^l{\displaystyle \frac{\mathrm{\Delta }\beta \mathrm{\Phi }^{ep}}{|x_{pt}^l|}}\right]`$ (35)
$`+{\displaystyle \frac{\alpha }{\mathrm{det}|\psi _{ab}^{n,1}|_s}}{\displaystyle \frac{\mathrm{det}|\psi _{ab}^{n,1}|_s}{\alpha }}\},`$ (36)
$`\mathrm{with}A_{pt}^l={\displaystyle \frac{r_{pt}^l|r_{pt}}{|r_{pt}^l|}},B_{pt}^l={\displaystyle \frac{x_{pt}^l|x_{pt}}{|x_{pt}^l|}}.`$ (37)
The structure of Eqs. (30, 37) is obvious: we have separated the classical ideal gas part (first term). The ideal quantum part in excess of the classical one and the correlation contributions are contained in the integral term, where the second line results from the ionic correlations (first term) and the e-e and e-i interaction at the first vertex (second and third terms respectively). The third and fourth lines are due to the further electronic vertices and the explicit temperature dependence \[in Eq. (30) and volume dependence in Eq. (37)\] of the exchange matrix, respectively. The main advantage of Eqs. (30, 37) is that the explicit sum over permutations has been converted into the spin determinant which can be computed very efficiently using standard linear algebra methods. Furthermore, each of the sums in curly brackets in Eqs. (30, 37) is bounded as the number of vertices increases, $`n\mathrm{}`$, and is thus well suited for efficient Monte Carlo simulations. Notice also that Eqs. (30, 37) contain the important limit of an ideal quantum plasma in a natural way .

## IV Comparison of direct and restricted PIMC simulations

Expressions (30, 37) are well suited for numerical evaluation using Monte Carlo techniques, e.g. . In our Monte Carlo scheme we used three types of steps, where either electron or proton coordinates, $`r_i`$ or $`q_i`$ or inidividual electronic beads $`\xi _i^{(k)}`$ were moved until convergence of the calculated values was reached. Our procedure has been extensively tested. In particular, we found from comparison with the known analytical expressions for pressure and energy of an ideal Fermi gas that the Fermi statistics is very well reproduced . Further, we performed extensive tests for few–electron systems in a harmonic trap where, again, the analytically known limiting behavior (e.g. energies) is well reproduced . For the present simulations of dense hydrogen, we varied both the particle number and the number of time slices (beads). As a result of these tests, we found that to obtain convergent results for the thermodynamic properties of dense hydrogen, particle numbers $`N_e=N_p=50`$ and beads numbers in the range of $`n=6\mathrm{}20`$ are adequate .

We will now compare our results with some available results obtained by the Monte Carlo technique developed by the Urbana group . We may first state that both Monte Carlo techniques differ in several fundamental points, so that they are essentially independent approaches. Let us briefly outline the main differences between the technique developed in Urbana, known as the restricted PIMC scheme  and references therein, and the approach described here. These authors performed simulations with 32 electrons and protons; their restricted PIMC scheme required to use a rather small time step assuring $`1/\mathrm{\Delta }\beta 210^6K`$. Also, the treatment of the interactions differs from our scheme: the authors of Ref. perform a numerical solution of the Bloch equation for the two-particle density matrix whereas we use an analytical approximation for the effective pair interaction (based on the Kelbg potential, see above). Finally, Ref. approximately computes the nodal surface of the density matrix using a variational ansatz which is then used to restrict the integrations to the region of positive density matrix. For more details regarding the restricted PIMC simulations, see Refs. .

Let us now turn to a comparison of the numerical results. The restricted PIMC simulation data for dense hydrogen are taken from Ref. . A comparison of results for the pressure and the internal energy for a fixed value of the density ($`r_s=1.86`$) is shown in Fig. 1 and TABLE I. At high temperatures, above 50,000 K, where only a small fraction of atoms is expected, the agreement is rather good. This is remarkable since the nonideality and the degeneracy reach values of 3 and 10, respectively. This result demonstrates that, at least for $`r_s1`$ and for $`T50,000K`$ both methods yield results which are more or less equivalent. At $`T<50,000K`$, where partial ionization is expected, we still observe a reasonable agreement of both approaches, however, we see also that the differences start to grow. The reasons for that are manyfold. From our results we conclude that the main problem is not the bound state formation - atoms and molecules are well described by the two PIMC simulations which use a physical picture which does not involve any artificial distinction between free and bound electrons, e.g. . On the other hand, with growing degeneracy $`n\mathrm{\Lambda }^3`$, both PIMC methods become less reliable, and a detailed analysis, although being very desirable, will have to be based on more extensive calculations in the future.

Further we present in TAB. I also Padé results for the weakly nonideal region. We find good agreement with the PIMC results for $`T>10^5`$ K. Details on the method of these analytical calculations will be discussed in the next section.

## V Comparison with analytical approximations for the thermodynamic functions of strongly ionized dense plasmas

In this section we give a comparison of the available data points from direct PIMC calculations with analytical estimates based on Padé approximations for strongly ionized plasmas . The comparison concentrates on H-plasmas in a region in the density-temperature plane with the following borders

$`0.2`$ $``$ $`\mathrm{\Gamma }1.6,`$ (38)
$`0.2`$ $``$ $`n_e\mathrm{\Lambda }_e^35,`$ (39)
which will be called “rhomb of moderate nonideality and moderate degeneracy“ (see bold rhomb in Fig. 1). With respect to analytical treatment, this rhombic region is of particular difficulty since none of the known analytic limiting expressions is valid. Further we calculated several points for $`r_s=1.86`$ and $`\mathrm{\Gamma }2`$ which correspond to the PIMC data discussed in the previous section and also an isotherm at $`T=50,000K`$ including some data at higher density, outside the rhomb, cf. Fig. 1 for an overviev.

We demonstrate below that the Padé approximations which interpolate between the limits where theoretical results are available are a useful tool for the description of the available data points, at least for the case of moderate nonideality $`\mathrm{\Gamma }1.6`$ and moderate degeneracy $`n_e\mathrm{\Lambda }^35`$. The Padé approximations which we use here were constructed in earlier work, , from the known analytical results for limiting cases of low density  and high density . The structure of the Padé approximations was devised in such a way that they are analytically exact up to quadratic terms in the density (up to the second virial coefficient) and interpolate between the virial expansions and the high-density asymptotic expressions . The formation of bound states was taken into account by using a chemical picture. This means the plasma is considered as a mixture of free electrons, free ions, atoms and molecules which are in chemical equilibrium, being described by mass action laws or minimization of the free energy .

We follow in large here this cited work, only the contribution of the ion-ion interaction which is, in most cases, the largest one, was substantially improved following recent work of Kahlbaum, who succeeded in describing the available classical Monte Carlo data for the ions by accurate Padé approximations . By using Kahlbaum’s formulas we achieve a rather accurate description of the thermodynamics in the region where the plasma behaves like a classical one-component ion plasma imbedded into a sea of nearly ideal electrons. This is the region where the electrons are strongly degenerate

$$n_e\mathrm{\Lambda }_e^31\text{and}r_s1,$$
(40)
and the ions are still classical but nonideal

$`\mathrm{\Gamma }1\text{and}n_i\mathrm{\Lambda }_i^31.`$ (41)
This region lies in the upper left corner of Fig. 1.

With respect to the chemical picture we restrict ourselves to the region of strong ionization where the number of atoms is still relatively low and where the fraction of molecules is small as well, see below. We will discuss here only the general structure of the Padé formulae. For example, the internal energy density of the plasma is given by

$$u=u_{id}+u_{int}.$$
(42)
Here $`u_{id}`$ is the internal energy of an ideal plasma consisting of Fermi electrons, classical protons and classical atoms, and $`u_{int}`$ is the interaction energy

$$u_{int}=u_{ii}+u_{ee}+u_{ie}+u_{\mathrm{vdW}}.$$
(43)
The interaction contribution to the internal energy consists of four terms:

* Ion-ion interaction contribution: this term which, in general, yields the largest contribution is generated by the OCP subsystem of the protons. For the OCP energy of protons many expressions are available, e.g. . We have used here the most precise formula due to Kahlbaum  which interpolates between the Debye region, $`u_{ii}\mathrm{\Gamma }^{3/2}`$, and the high density fluid, $`u_{ii}\mathrm{\Gamma }`$.
* Electron-electron interaction: This term corresponds to the OCP energy of the electron subsystem. We used the rather simple expressions used in earlier work .
* Electron-proton interaction: This term corresponds to the interaction between the two OCP subsystems which is mostly due to polarization effects. Again, we used the rather simple expressions proposed in earlier work .
* Van der Waals contribution: In the region of densities and temperatures defined above this contribution gives only a small correction. Therefore, this term was approximated here in the simplest way by a second virial contribution. The neutral particles were treated as hard spheres.

In the region of densities which are studied here, molecules do not play a role, therefore, the formation of molecules was taken into account only in a very rough approximation according to Ref. . The number density of the neutrals was calculated on the basis of a nonideal Saha equation. We restricted this comparison to a region where the number density of neutrals is relatively small, the degree of ionization being larger than $`75\%`$.

The contributions to the pressure were calculated, in part, from scaling relations e.g. we used $`p_{ii}=u_{ii}/3`$, and, for the other (smaller) contributions, by numerical differentiation of the free energy given earlier . In a similar way, the chemical potential which appears in the nonideal Saha equation was obtained. For the partition function in the Saha equation we used the Brillouin-Planck-Larkin expression . The solution of the nonideal Saha equation which determines the degree of ionization (the density of the atoms) was solved by up to 100 iterations starting from the ideal Saha equation.

Since all the expression described so far are given in analytic form, the calculation of about 1000 data points for energy and pressure takes less than a minute on a PC. The result of our calculations for density-temperature points in the “rhomb of moderate nonideality and moderate degeneracy” are given in Figs. 2,3. Further, we give in TAB. I several data points obtained from the Padé formulas. Since the Padé formulas used here do not apply to low temperatures, we included in TAB. I only Padé data for $`T>10^5K`$.

Summarizing the results for the internal energy and for the pressure we find that the Padé results, with a few exceptions, agree well with the PIMC data in the region of the density temperature plane, where $`\mathrm{\Gamma }1.6`$ and $`n\mathrm{\Lambda }^35`$. The agreement is particularly good for the energies. \[The larger deviations for the pressure may be due to the numerical differentiation.\] In fact, the Padé formulas used here in combination with the chemical picture works only in the case that the plasma is strongly ionized, i.e. the degree of ionization is larger than $`75\%`$. The description of the region where a higher percentage of atoms and, due to this, also molecules is present needs a more refined chemical picture .

Finally, we compare the Padé and PIMC data along the isotherm $`T=50,000K`$ which is given in Fig. 5. This figure shows the transition from a classical ideal gas (low density) to a nearly ideal quantum gas (limit of high density). In the central part, $`n10^{19}\mathrm{cm}^310^{25}\mathrm{cm}^3`$, Coulomb interaction leads to strong deviations from the behavior of an ideal plasma. The strong increase of the energy at high density is due to the Mott effect and to the increase of the ideal quantum contribution to the electron energy. Comparing the Padé and PIMC results, we find good agreement up to electron densities $`n=10^{22}\mathrm{cm}^3`$. For higher densities, the deviations are growing. For intermediate densities, $`n10^{22}\mathrm{cm}^310^{24}\mathrm{cm}^3`$ the PIMC data are more reliable. On the other hand, in the limit of very high density, $`r_s1`$, the Padé results are known to correctly approach the ideal quantum plasma limit whereas the PIMC data should be regarded as preliminary due to the extremely high electron degeneracy. Interestingly, we find that at high density the Padé data approaches the ideal curves earlier than the PIMC data which is important for further improvement of the presented Monte Carlo approach.

## VI Discussion

This work is devoted to the investigation of the thermodynamic properties of hot dense partially ionized plasmas in the pressure range between 0.1 and 100 Mbar. Most of the new results are based on a Quantum Monte Carlo study of a correlated proton-electron system with degenerate electrons and classical protons. In this paper, we gave a detailed derivation of improved estimators for the internal energy and the equation of state for use in direct fermionic path integral simulations. Also, we gave a rigorous justification for the use of an effective quantum pair potential (Kelbg potential) in PIMC simulations.

Further, we compared our direct PIMC results with independent restricted PIMC data of Militzer and Ceperley for one isochor corresponding to $`r_s=1.86`$, Fig. 2. We found very good quantitative agreement between the two PIMC methods for temperatures in the range of $`50,000KT10^6K`$, where $`\mathrm{\Gamma }3`$ and $`n_e\mathrm{\Lambda }_e^310`$. This region is particularly complicated as here pressure and temperature ionization occur and, therefore, an accurate and consistent treatment of scattering and bound states is crucial. This agreement is remarkable because the two simulation methods are completely independent and use essentially different approximations. We, therefore, expect that the results for the thermodynamic properties of high pressure hydrogen plasmas in this temperature-density range are reliable within the limits of the simulation accuracy. This is the main result of the present paper.

In future work, it will be important to extend the range of agreement. To analyze the deviations between the simulation methods, we also included some data for $`r_s=1.86`$ and lower temperatures, $`10,000KT50,000K`$, Fig. 2. At this point, no conclusive answer about the reasons of the deviations can be given. For these parameters, the electron degeneracy is growing rapidly and, therefore, each of the simulation methods is becoming less reliable. So these data should be regarded as preliminary results which will be useful for future improvements of the simulations.

Furthermore, the Monte Carlo results allowed us to develop and test analytical approximations of Padé-type which are improvements of earlier approximations  in a region in the density-temperature plane bounded by $`\mathrm{\Gamma }1.6`$ and $`n_e\mathrm{\Lambda }_e^35`$. This is a region of moderate nonideality and degeneracy and high degree of ionization. We have shown that for these parameters, the Padé approximations which interpolate between the limits where theoretical results are available agree well with the Monte Carlo data, cf. Figs. 2-4 and Table I. Thus, these approximations provide a useful tool for the description of these plasmas which include hydrogen at a pressur between 0.1 and 100 Mbar. At lower temperature, deviations from the Monte Carlo data are growing, cf. Fig. 2. This is mostly due to the growing role of bound states. Whether the Padé approximations, in combination with an improved chemical picture (mass action law), continue to work at lower temperatures, has still to be explored, first steps are under way .

Also, we showed some data for $`T=50,000K`$ and higher pressure, up to $`p10^6`$ Mbar, Fig. 5. Here the Monte Carlo simulations are particularly difficult due to the high electron degeneracy, and they can benefit from the Pade simulations, as the latter correctly reproduce the high-density limit, $`r_s1`$.

## VII Acknowledgements

We acknowledge stimulating discussions with W.D. Kraeft, D. Kremp and M. Schlanges. We thank D.M. Ceperley and B. Militzer for discussions on PIMC concepts and for providing us with the data of Ref.  prior to publication. Further we thank J. Ortner for informing us about an alternative derivation for the off-diagonal elements of the interaction potential. This work was made possible by generous support from the Deutsche Forschungsgemeinschaft (Mercator-Programm) for VSF and by a grant for CPU time at the NIC Jülich.