text:
# 1 Introduction

## 1 Introduction

In the previous paper, one of the authors shows that the SUSY flavor problem can be solved in $`E_6`$ unification by using non-abelian horizontal gauge symmetry and anomalous $`U(1)_A`$ gauge symmetry, whose anomaly is cancelled by the Green-Schwarz mechanism, even if large neutrino mixing angles are obtained. It is essential that the fundamental representation $`\mathrm{𝟐𝟕}`$ of $`E_6`$ has two $`\overline{\mathrm{𝟓}}`$ fields of $`SU(5)`$. Actually $`\mathrm{𝟐𝟕}`$ is decomposed as

$$\mathrm{𝟐𝟕}\underset{\mathrm{𝟏𝟔}_1}{\underset{}{[\mathrm{𝟏𝟎}_{(1,1)}+\overline{\mathrm{𝟓}}_{(1,3)}+\mathrm{𝟏}_{(1,5)}]}}+\underset{\mathrm{𝟏𝟎}_2}{\underset{}{[\overline{\mathrm{𝟓}}_{(2,2)}+\mathrm{𝟓}_{(2,2)}]}}+\underset{\mathrm{𝟏}_4}{\underset{}{[\mathrm{𝟏}_{(4,0)}]}}$$
(1.1)
under $`E_6SO(10)\times U(1)_V^{}SU(5)\times U(1)_V^{}\times U(1)_V`$, where the representation of $`SO(10)\times U(1)_V^{},SU(5)\times U(1)_V^{}\times U(1)_V`$ are explicitly denoted in the above. If we introduce three $`\mathrm{𝟐𝟕}`$ fields $`\mathrm{\Psi }_i`$ $`(i=1,2,3)`$ for three generation quarks and leptons, three of six $`\overline{\mathrm{𝟓}}`$ fields become massive with three $`\mathrm{𝟓}`$ fields after breaking $`E_6`$ into $`SU(5)`$, and the remaining fields ($`3\times \overline{\mathrm{𝟓}}`$) remain massless. In the $`6\times 3`$ mass matrix for $`\overline{\mathrm{𝟓}}`$ and $`\mathrm{𝟓}`$ fields, it is natural to expect that the elements for the third generation field $`\mathrm{\Psi }_3`$ become larger to realize larger Yukawa couplings than the first and second generation fields $`\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$. Therefore, all the three massless modes of $`\overline{\mathrm{𝟓}}`$ come mainly from the first two generation fields $`\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$. This structure is interesting because it can explain larger mixing angles of lepton sector than of quark sector as discussed in Ref.. Moreover, if we introduce non-abelian horizontal symmetry $`SU(2)_H`$ and take the first two generation fields as doublet, then all three generation $`\overline{\mathrm{𝟓}}`$ fields have degenerate sfermion masses, which are very important to suppress flavor changing neutral current (FCNC) processes with large neutrino mixing angles as discussed in Ref..

In the above arguments, $`E_6`$ gauge group plays an important role. Actually, it is essential that a single field includes two $`\overline{\mathrm{𝟓}}`$ fields to realize large neutrino mixing angles with suppressing FCNC processes. However, in order to break $`E_6`$ into the standard model (SM) gauge group $`SU(3)_C\times SU(2)_L\times U(1)_Y`$, adjoint Higgs fields $`\mathrm{𝟕𝟖}`$ are required, which may not be so easily realized in the framework of superstring models. To avoid the adjoint Higgs, it is a simple way to adopt non-simple group as a unification group. Which kinds of non-simple group do not spoil the above interesting features? The answer is simple. In order to satisfy the essential point that two fields with the same quantum number under the SM gauge group are included in a single multiplet, $`SU(2)_E`$, which is a subgroup of $`E_6`$ group and rotates $`(\overline{\mathrm{𝟓}}_{(1,3)},\overline{\mathrm{𝟓}}_{(2,2)})`$ and $`(\mathrm{𝟏}_{(1,5)},\mathrm{𝟏}_{(4,0)})`$ as doublets, is sufficient. Therefore, it is interesting to consider the unification group which include $`SU(2)_E`$. The $`SU(3)^3E_6`$ is in the case and we know that realistic $`SU(3)^3`$ model can be straightforwardly constructed, in which doublet-triplet splitting problem is solved and realistic quark and lepton mass matrices are obtained including large neutrino mixing angles. Therefore, if we introduce non-abelian horizontal symmetry in addition to $`SU(3)^3`$, FCNC processes can be naturally suppressed with large neutrino mixing angles. In this paper, we consider another non-simple gauge group, $`SO(10)_F\times U(1)_{V_F^{}}`$, which can include $`SU(2)_E`$ because of the unusual embedding of the SM gauge group. We show that in this model, doublet-triplet splitting is realized by missing partner mechanism. The original missing partner mechanism was introduced in $`SU(5)`$ unification group, but it requires several large dimensional representation Higgs fields. To avoid the large dimensional Higgs fields, flipped $`SU(5)`$ has been considered. It is known that the gauge group $`SU(5)_F\times U(1)_X`$ cannot be unified into $`SO(10)`$ without spoiling the missing partner mechanism, but we show that $`SO(10)_F\times U(1)_{V_F^{}}E_6`$ can embed the flipped $`SU(5)`$ without spoiling the missing partner mechanism. As noted in the above, the flipped $`SO(10)`$ gauge group includes $`SU(2)_E`$, that is important to solve the SUSY flavor problem by introducing non-abelian horizontal gauge symmetry and anomalous $`U(1)_A`$ gauge symmetry.

## 2 Review of flipped $`SU(5)`$ model

We briefly review the flipped $`SU(5)`$ model and the reason why the flipped $`SU(5)`$ model cannot be embedded in $`SO(10)`$ GUT.

It is well-known that one family standard model fermions $`Q(\mathrm{𝟑},\mathrm{𝟐})_{\frac{1}{6}}`$, $`U^c(\overline{\mathrm{𝟑}},\mathrm{𝟏})_{\frac{2}{3}}`$, $`D^c(\overline{\mathrm{𝟑}},\mathrm{𝟏})_{\frac{1}{3}}`$, $`L(\mathrm{𝟏},\mathrm{𝟐})_{\frac{1}{2}}`$, and $`E^c(\mathrm{𝟏},\mathrm{𝟏})_1`$ plus the right-handed neutrino $`N^c(\mathrm{𝟏},\mathrm{𝟏})_0`$ under the SM gauge group $`SU(3)_C\times SU(2)_L\times U(1)_Y`$ are unified into an $`SO(10)`$-spinorial $`\mathrm{𝟏𝟔}`$ superfield:

$$\mathrm{\Psi }(\mathrm{𝟏𝟔})\mathrm{𝟏𝟎}_\mathrm{\Psi }(\mathrm{𝟏𝟎}_1)+\overline{\mathrm{𝟓}}_\mathrm{\Psi }(\overline{\mathrm{𝟓}}_3)+\mathrm{𝟏}_\mathrm{\Psi }(\mathrm{𝟏}_5),$$
(2.1)
where the decomposition is specified into $`SU(5)\times U(1)_V`$. The matter content of the flipped $`SU(5)`$ models can be obtained from the corresponding assignment of the standard $`SU(5)`$ GUT model by means of the “flipping” $`U^cD^c`$, $`N^cE^c`$:

$`\mathrm{𝟏𝟎}_\mathrm{\Psi }`$ $`=`$ $`(Q,D^c,N^c)`$
$`\overline{\mathrm{𝟓}}_\mathrm{\Psi }`$ $`=`$ $`(U^c,L)`$ (2.2)
$`\mathrm{𝟏}_\mathrm{\Psi }`$ $`=`$ $`E^c.`$
It is important that if $`\mathrm{𝟏𝟎}_1`$ representation Higgs $`\mathrm{𝟏𝟎}_C`$ is introduced, $`SU(5)\times U(1)_X`$ can be broken into the standard model gauge group $`SU(3)_C\times SU(2)_L\times U(1)_Y`$ by the vacuum expectation value (VEV) of the component of $`N^c`$. Here, the hypercharge operator is written

$$Y=\frac{1}{5}(XY^{}),$$
(2.3)
where $`Y^{}`$ is the generator of $`SU(5)_F`$ which commutes with $`SU(3)_C\times SU(2)_L`$. Then the $`SO(10)`$-vectorial $`\mathrm{𝟏𝟎}`$ superfield decomposed as

$$H(\mathrm{𝟏𝟎})\mathrm{𝟓}_H(\mathrm{𝟓}_2)+\overline{\mathrm{𝟓}}_H(\overline{\mathrm{𝟓}}_2)$$
(2.4)
includes the SM doublet Higgs $`H_d=L^{}`$ and $`H_u=\overline{L}^{}`$ as

$`\mathrm{𝟓}_H`$ $`=`$ $`(\overline{D}^c,L^{})`$
$`\overline{\mathrm{𝟓}}_H`$ $`=`$ $`(D^c,\overline{L}^{}),`$ (2.5)
where $`D^c`$ and $`L^{}`$ have the same quantum number of SM gauge group as $`D^c`$ and $`L`$, respectively. If we introduce interactions in the superpotential as

$$W_{MP}=\mathrm{𝟏𝟎}_C\mathrm{𝟏𝟎}_C\mathrm{𝟓}_H+\overline{\mathrm{𝟏𝟎}}_{\overline{C}}\overline{\mathrm{𝟏𝟎}}_{\overline{C}}\overline{\mathrm{𝟓}}_H,$$
(2.6)
only the triplet Higgs $`\overline{D}^c`$ and $`D^c`$ can be superheavy with $`D^c`$ in $`\mathrm{𝟏𝟎}_C`$ and $`\overline{D}^c`$ in $`\overline{\mathrm{𝟏𝟎}}_{\overline{C}}`$, respectively, by developing the VEVs of $`\mathrm{𝟏𝟎}_C`$ and $`\overline{\mathrm{𝟏𝟎}}_{\overline{C}}`$, but the doublet Higgs $`L^{}`$ and $`\overline{L}^{}`$ have no partner and remain massless. This is essential of the missing partner mechanism in the flipped $`SU(5)`$ model.

Unfortunately, this missing partner mechanism in the flipped $`SU(5)`$ model cannot be extended to $`SO(10)`$ unification. In $`SO(10)`$ unification the interactions (2.6) are included in the $`SO(10)`$ symmetric interactions $`C(\mathrm{𝟏𝟔})C(\mathrm{𝟏𝟔})H(\mathrm{𝟏𝟎})`$ and $`\overline{C}(\overline{\mathrm{𝟏𝟔}})\overline{C}(\overline{\mathrm{𝟏𝟔}})H(\mathrm{𝟏𝟎})`$, which include also

$$\mathrm{𝟏𝟎}_C\overline{\mathrm{𝟓}}_C\overline{\mathrm{𝟓}}_H+\overline{\mathrm{𝟏𝟎}}_{\overline{C}}\mathrm{𝟓}_{\overline{C}}\mathrm{𝟓}_H.$$
(2.7)
Through these interactions, the doublet Higgs $`(\overline{L}^{})_H`$ and $`(L^{})_H`$ become superheavy with $`L_C`$ and $`(L^{})_{\overline{C}}`$, respectively, by developing the VEVs of $`\mathrm{𝟏𝟎}_{\overline{C}}`$ and $`\overline{\mathrm{𝟏𝟎}}_{\overline{C}}`$. (In this paper, $`X^{}`$ is a component of $`\overline{\mathrm{𝟏𝟔}}`$ of $`SO(10)`$ and denotes the complex conjugate representation of $`X`$ which is a component of $`\mathrm{𝟏𝟔}`$ of $`SO(10)`$. ) Therefore, doublet-triplet splitting is spoiled by this extension.111Of course, if we neglect the component fields $`\overline{\mathrm{𝟓}}_C`$ and $`\mathrm{𝟓}_{\overline{C}}`$ by hand, such extension becomes possible .

In the next section, we show that the missing partner mechanism of the flipped $`SU(5)`$ model can be embedded in $`SO(10)_F\times U(1)_{V_F^{}}`$ unification group.

## 3 Flipped $`SO(10)`$ model

As noted in the introduction, $`\mathrm{𝟐𝟕}`$ of $`E_6`$ is decomposed as

$$\mathrm{𝟐𝟕}\underset{\mathrm{𝟏𝟔}_1}{\underset{}{[\mathrm{𝟏𝟎}_{(1,1)}+\overline{\mathrm{𝟓}}_{(1,3)}+\mathrm{𝟏}_{(1,5)}]}}+\underset{\mathrm{𝟏𝟎}_2}{\underset{}{[\overline{\mathrm{𝟓}}_{(2,2)}+\mathrm{𝟓}_{(2,2)}]}}+\underset{\mathrm{𝟏}_4}{\underset{}{[\mathrm{𝟏}_{(4,0)}]}}$$
(3.1)
under $`E_6SO(10)\times U(1)_V^{}SU(5)\times U(1)_V^{}\times U(1)_V`$. There are two ways to embed the flipped $`SU(5)`$ matters $`\mathrm{𝟏𝟎}_\mathrm{\Psi }=(Q,D^c,N^c)`$, $`\overline{\mathrm{𝟓}}_\mathrm{\Psi }=(U^c,L)`$ and $`\mathrm{𝟏}_\mathrm{\Psi }=E^c`$ in the above decomposition of $`\mathrm{𝟐𝟕}`$ of $`E_6`$ into $`SO(10)\times U(1)_V^{}`$. As discussed in the previous section, the usual embedding $`SU(5)_F\times U(1)_X`$ in $`SO(10)`$,

$$\underset{\mathrm{𝟏𝟔}_1}{\underset{}{[\mathrm{𝟏𝟎}_\mathrm{\Psi }+\overline{\mathrm{𝟓}}_\mathrm{\Psi }+\mathrm{𝟏}_\mathrm{\Psi }]}}+\underset{\mathrm{𝟏𝟎}_2}{\underset{}{[\overline{\mathrm{𝟓}}_H+\mathrm{𝟓}_H]}}+\underset{\mathrm{𝟏}_4}{\underset{}{[\mathrm{𝟏}_S]}},$$
(3.2)
where $`\mathrm{𝟓}_H=(\overline{D}^c,L^{})`$, $`\overline{\mathrm{𝟓}}_H=(D^c,\overline{L}^{})`$ and $`\mathrm{𝟏}_S`$ is singlet under $`SU(5)_F\times U(1)_X`$, spoils the missing partner mechanism. The other embedding can be obtained by means of the “flipping” $`\overline{\mathrm{𝟓}}_\mathrm{\Psi }\overline{\mathrm{𝟓}}_H`$ and $`\mathrm{𝟏}_\mathrm{\Psi }\mathrm{𝟏}_S`$:

$$\underset{\mathrm{𝟏𝟔}_1}{\underset{}{[\mathrm{𝟏𝟎}_\mathrm{\Psi }+\overline{\mathrm{𝟓}}_H+\mathrm{𝟏}_S]}}+\underset{\mathrm{𝟏𝟎}_2}{\underset{}{[\overline{\mathrm{𝟓}}_\mathrm{\Psi }+\mathrm{𝟓}_H]}}+\underset{\mathrm{𝟏}_4}{\underset{}{[\mathrm{𝟏}_\mathrm{\Psi }]}}.$$
(3.3)
In this embedding, if $`\mathrm{𝟏}_S`$ component of $`\mathrm{𝟏𝟔}_1`$ field have non-vanishing VEV, $`SO(10)_F\times U(1)_{V_F^{}}`$ is broken into $`SU(5)_F\times U(1)_X`$. Here, the operator $`X`$ is obtained as

$$X=\frac{1}{4}(5V_F^{}V_F),$$
(3.4)
where $`V_F`$ is the generator of $`SO(10)_F`$ which commute with $`SU(5)_F`$. The hypercharge operator is

$$Y=\frac{1}{5}(XY^{})=\frac{1}{20}(5V_F^{}V_F4Y^{}).$$
(3.5)
Note that the each $`SU(2)_E`$ doublet $`(D^c,D^c)`$, $`(L^{},L)`$ and $`(N^c,S)`$, which has the same quantum number of SM gauge group, is included into a single multiplet $`\mathrm{𝟏𝟔}_1`$, $`\mathrm{𝟏𝟎}_2`$ and $`\mathrm{𝟏𝟔}_1`$, respectively. This means that $`SU(2)_E`$ is embedded in $`SO(10)_F`$.

We introduce two pairs of Higgs fields $`[\mathrm{\Phi }(\mathrm{𝟏𝟔}_1),\overline{\mathrm{\Phi }}(\overline{\mathrm{𝟏𝟔}}_1)]`$ and $`[C(\mathrm{𝟏𝟔}_1),\overline{C}(\overline{\mathrm{𝟏𝟔}}_1)]`$ to break $`SO(10)_F\times U(1)_{V_F^{}}`$ into the SM gauge group. Supposing that the VEVs $`|\mathrm{\Phi }|=|\overline{\mathrm{\Phi }}|`$ breaks $`SO(10)_F\times U(1)_{V_F^{}}`$ into $`SU(5)_F\times U(1)_X`$, the components $`\mathrm{𝟏𝟎}_\mathrm{\Phi }`$ and $`\overline{\mathrm{𝟏𝟎}}_{\overline{\mathrm{\Phi }}}`$ are absorbed by the Higgs mechanism. The VEVs $`|C|=|\overline{C}|`$ break $`SU(5)_F\times U(1)_X`$ into the SM gauge group, and the components $`Q`$ and $`N^c`$ are absorbed by the Higgs mechanism. All the remaining components $`\overline{\mathrm{𝟓}}_\mathrm{\Phi }`$, $`\mathrm{𝟓}_{\overline{\mathrm{\Phi }}}`$, $`\overline{\mathrm{𝟓}}_C`$, $`\mathrm{𝟓}_{\overline{C}}`$, $`(D^c)_C`$ and $`(D^c)_{\overline{C}}`$ must be massive except a pair of doublets. For example, through the interactions in the superpotential,

$$W_{SO(10)}=\overline{\mathrm{\Phi }}\overline{\mathrm{\Phi }}CC+\overline{C}\overline{C}\mathrm{\Phi }\mathrm{\Phi },$$
(3.6)
which include the interactions (2.6) after developing the VEVs $`|\mathrm{\Phi }|=|\overline{\mathrm{\Phi }}|`$, pairs $`[(D^c)_{\overline{\mathrm{\Phi }}},(D^c)_C]`$ and $`[(D^c)_\mathrm{\Phi },(D^c)_{\overline{C}}]`$ become massive. If we introduce the mass term for $`C`$ and $`\overline{C}`$, then only $`(\overline{L}^{})_\mathrm{\Phi }`$ and $`(\overline{L}^{})_{\overline{\mathrm{\Phi }}}`$ remain massless, namely, doublet-triplet splitting is realized. There are several interactions which unstabilize the doublet-triplet splitting. For example, the terms $`\overline{\mathrm{\Phi }}\mathrm{\Phi }F(\overline{C}C,\overline{\mathrm{\Phi }}\mathrm{\Phi })`$ give directly the doublet Higgs mass, so they must be forbidden. (We will return to this subject lator in a concrete model.)

We assume that three generation matter fields $`\mathrm{\Psi }_i(\mathrm{𝟐𝟕})=\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}+\mathrm{𝟏𝟎}_{\mathrm{\Psi }_i}+\mathrm{𝟏}_{\mathrm{\Psi }_i}`$ $`(i=1,2,3)`$ respect $`E_6`$ symmetry. It is an easy way to guarantee the cancellation of gauge anomaly. Among the three generation matter fields $`\mathrm{\Psi }_i`$, there are six fields which have the same quantum number under the SM gauge group as $`(D^c,L)`$. Only three linear combinations of these fields become quarks and leptons, and other modes become superheavy with the three $`(\overline{D}^c,\overline{L}^{})`$ fields through the interactions $`\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟎}_{\mathrm{\Psi }_j}\mathrm{\Phi }`$ and $`\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟎}_{\mathrm{\Psi }_j}C`$ by developing the VEVs of $`\mathrm{\Phi }`$ and $`C`$. It is interesting that up-type Yukawa coupling can be obtained from the renormalizable interactions $`\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟎}_{\mathrm{\Psi }_j}\mathrm{\Phi }`$, because $`O(1)`$ top Yukawa coupling can be naturally realized. But Yukawa couplings of down quark sector and of charged lepton sector are obtained from the higher dimensional interactions $`\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟔}_{\mathrm{\Psi }_j}\overline{C}\overline{\mathrm{\Phi }}`$ and $`\mathrm{𝟏𝟎}_{\mathrm{\Psi }_i}\mathrm{𝟏}_{\mathrm{\Psi }_j}\overline{C}\overline{\mathrm{\Phi }}`$, respectively. Because we have six singlets $`N_i^c`$ and $`S_i`$ in the matter sector, the mass matrix for right-handed neutrino becomes $`6\times 6`$ matrix which are obtained from the interactions $`\mathrm{\Psi }_i\mathrm{\Psi }_j\overline{\mathrm{\Phi }}\overline{\mathrm{\Phi }}`$, $`\mathrm{\Psi }_i\mathrm{\Psi }_j\overline{\mathrm{\Phi }}\overline{C}`$ and $`\mathrm{\Psi }_i\mathrm{\Psi }_j\overline{C}\overline{C}`$. Yukawa couplings of Dirac neutrino sector are obtained from the interactions $`\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟎}_{\mathrm{\Psi }_j}\mathrm{\Phi }`$. Therefore, the mass terms of all quarks and leptons can be obtained in this scenario.

Unfortunately, as in the flipped $`SU(5)`$ model, this missing partner mechanism in the flipped $`SO(10)`$ model cannot be extended to $`E_6`$ unification. In $`E_6`$ unification the interactions (3.6) are included in the $`E_6`$ symmetric interactions $`\mathrm{\Phi }(\mathrm{𝟐𝟕})\mathrm{\Phi }(\mathrm{𝟐𝟕})\overline{C}(\overline{\mathrm{𝟐𝟕}})\overline{C}(\overline{\mathrm{𝟐𝟕}})`$ and $`\overline{\mathrm{\Phi }}(\overline{\mathrm{𝟐𝟕}})\overline{\mathrm{\Phi }}(\overline{\mathrm{𝟐𝟕}})C(\mathrm{𝟐𝟕})C(\mathrm{𝟐𝟕})`$, which include also $`\mathrm{𝟏𝟔}_\mathrm{\Phi }\mathrm{𝟏𝟎}_\mathrm{\Phi }\mathrm{𝟏𝟎}_{\overline{C}}\overline{\mathrm{𝟏𝟔}}_{\overline{C}}`$ and $`\mathrm{𝟏𝟔}_C\mathrm{𝟏𝟎}_C\mathrm{𝟏𝟎}_{\overline{\mathrm{\Phi }}}\overline{\mathrm{𝟏𝟔}}_{\overline{\mathrm{\Phi }}}`$ of $`SO(10)_F`$. After developing the VEVs $`|\mathrm{\Phi }|=|\overline{\mathrm{\Phi }}|`$, these interactions give $`\mathrm{𝟓}_\mathrm{\Phi }\mathrm{𝟓}_{\overline{C}}\overline{\mathrm{𝟏𝟎}}_{\overline{C}}`$ and $`\overline{\mathrm{𝟓}}_{\overline{\mathrm{\Phi }}}\overline{\mathrm{𝟓}}_C\mathrm{𝟏𝟎}_C`$ of $`SU(5)_F`$, which give mass terms to doubet Higgs by taking non-vanishing VEVs $`|C|=|\overline{C}|`$. Therefore, doublet-triplet splitting is spoiled in this extension.

## 4 Flipped $`SO(10)`$ model with anomalous $`U(1)_A`$

It is important to find a concrete flipped $`SO(10)`$ model in which doublet-triplet splitting is realized with generic interactions and to examine whether the realistic quark and lepton mass matrices are realized or not. In a series of papers, we have pointed out that anomalous $`U(1)_A`$ symmetry plays an important role in solving various problems in SUSY grand unified theory (GUT) with generic interactions. This is mainly because the SUSY zero mechanism (holomorphic zero)222Note that if the total charge of an operator is negative, the $`U(1)_A`$ invariance and analytic property of the superpotential forbids the existence of the operator in the superpotential, since the Froggatt-Nielsen  field $`\mathrm{\Theta }`$ with negative charge cannot compensate for the negative total charge of the operator (the SUSY zero mechanism). can control various terms which must be forbidden.

In this section, we present a concrete flipped $`SO(10)`$ model with generic interaction by introducing anomalous $`U(1)_A`$ symmetry.

### 4.1 Higgs sector

The Higgs contents are listed in Table I.

Table I. The typical values of anomalous $`U(1)_A`$ charges are listed. $`\pm `$ is $`Z_2`$-parity and $`i=1,2`$.

|  | non-vanishing VEV | vanishing VEV |
| --- | --- | --- |
| $`\mathrm{𝟏𝟔}_1`$ | $`\mathrm{\Phi }(\varphi =0,)`$ $`C(c=2,+)`$ | $`\mathrm{\Phi }_i^{}(\varphi _i^{}=5,)`$ |
| $`\overline{\mathrm{𝟏𝟔}}_1`$ | $`\overline{\mathrm{\Phi }}(\overline{\varphi }=1,)`$ $`\overline{C}(\overline{c}=2,+)`$ | $`\overline{\mathrm{\Phi }}_i^{}(\overline{\varphi }_i^{}=4,)`$ |
| 1 | $`\mathrm{\Theta }(\theta =1,+)`$ $`\overline{Z}_i(\overline{z}_i=1,+)`$ $`Z(z=4,)`$ | $`S^{}(s^{}=8,+)`$ |

Following the general discussion on the determination of VEVs of the models with anomalous $`U(1)_A`$ charges, only the negatively charged fields can have non-vanishing VEVs. The scale of these VEVs are determined by the anomalous $`U(1)_A`$ charges as

$$\overline{\mathrm{\Phi }}\mathrm{\Phi }\lambda ^{(\varphi +\overline{\varphi })},\overline{C}C\lambda ^{(c+\overline{c})},$$
(4.1)
where $`\lambda `$ is the ratio of the VEV of Froggatt-Nielsen field $`\mathrm{\Theta }`$, which is essentially determined by the Fayet-Illiopoulos $`D`$-term parameter, to the cutoff $`\mathrm{\Lambda }`$. In this paper, we take $`\lambda `$ as around the Cabbibo angle $`\mathrm{sin}\theta _W0.22`$. If the $`\mathrm{𝟏}_{(1,5)}`$ component of $`\mathrm{\Phi }`$ and the $`\mathrm{𝟏}_{(1,5)}`$ component of $`\overline{\mathrm{\Phi }}`$ have non-vanishing VEVs, $`SO(10)_F\times U(1)_{V_F^{}}`$ is broken into $`SU(5)_F\times U(1)_X`$. The $`\mathrm{𝟏𝟎}_{(1,1)}`$ of $`\mathrm{\Phi }`$ and $`\overline{\mathrm{𝟏𝟎}}_{(1,1)}`$ of $`\overline{\mathrm{\Phi }}`$ are absorbed by the Higgs mechanism at that time. Moreover, if the $`\mathrm{𝟏𝟎}_{(1,1)}`$ component of $`C`$ and the $`\overline{\mathrm{𝟏𝟎}}_{(1,1)}`$ component of $`\overline{C}`$ have non-vanishing VEVs, $`SU(5)_F\times U(1)_X`$ is broken into the SM gauge group. Then the $`Q`$ component of $`\mathrm{𝟏𝟎}_{(1,1)}`$ of $`C`$ and the $`\overline{Q}`$ component of $`\overline{\mathrm{𝟏𝟎}}_{(1,1)}`$ of $`\overline{C}`$ are absorbed by the Higgs mechanism. Therefore, the remaining negatively charged fields except singlets under the SM gauge group are the $`\overline{\mathrm{𝟓}}_{(1,3)}`$ components of $`\mathrm{\Phi }`$ and $`C`$, the $`D^c`$ component of $`C`$, and the mirror components of $`\overline{\mathrm{\Phi }}`$ and $`\overline{C}`$. Among these negatively charged fields, no mass term appears because of the SUSY zero (holomorphic zero) mechanism. In order to make them massive, we have to take account of the positively charged fields $`\mathrm{\Phi }_i^{}`$ and $`\overline{\mathrm{\Phi }}_i^{}`$. Note that in a $`\mathrm{𝟏𝟔}_1`$ field, there are two colored Higgs $`D^c`$ and $`D^c^{}`$ because of $`SU(2)_E`$ symmetry, but only one doublet $`\overline{L}^{}`$. Therefore, the colored Higgs mass matrix becomes $`7\times 7`$ matrix $`M_T`$ which is given by

$$\begin{array}{cccccccccc}\overline{D}^c\backslash D^c& & \mathrm{𝟏𝟎}_C& \overline{\mathrm{𝟓}}_C& \overline{\mathrm{𝟓}}_\mathrm{\Phi }& \mathrm{𝟏𝟎}_{\mathrm{\Phi }_1^{}}& \mathrm{𝟏𝟎}_{\mathrm{\Phi }_2^{}}& \overline{\mathrm{𝟓}}_{\mathrm{\Phi }_1^{}}& \overline{\mathrm{𝟓}}_{\mathrm{\Phi }_2^{}}& \\ \overline{\mathrm{𝟏𝟎}}_{\overline{C}}& (\mathrm{}& 0& 0& 0& 0& 0& \lambda ^{\overline{c}+\varphi _1^{}\mathrm{\Delta }}& \lambda ^{\overline{c}+\varphi _2^{}\mathrm{\Delta }}& )\mathrm{}\\ \mathrm{𝟓}_{\overline{C}}& 0& 0& 0& \lambda ^{\overline{c}+\varphi _1^{}+\mathrm{\Delta }}& \lambda ^{\overline{c}+\varphi _2^{}+\mathrm{\Delta }}& 0& 0\\ \mathrm{𝟓}_{\overline{\mathrm{\Phi }}}& 0& 0& 0& 0& 0& \lambda ^{\overline{\varphi }+\varphi _1^{}}& \lambda ^{\overline{\varphi }+\varphi _2^{}}\\ \overline{\mathrm{𝟏𝟎}}_{\overline{\mathrm{\Phi }}_1^{}}& 0& \lambda ^{\overline{\varphi }_1^{}+c\mathrm{\Delta }}& 0& \lambda ^{\overline{\varphi }_1^{}+\varphi _1^{}}& \lambda ^{\overline{\varphi }_1^{}+\varphi _2^{}}& \lambda ^{\overline{\varphi }_1^{}+\varphi _1^{}\mathrm{\Delta }}& \lambda ^{\overline{\varphi }_1^{}+\varphi _1^{}\mathrm{\Delta }}\\ \overline{\mathrm{𝟏𝟎}}_{\overline{\mathrm{\Phi }}_2^{}}& 0& \lambda ^{\overline{\varphi }_2^{}+c\mathrm{\Delta }}& 0& \lambda ^{\overline{\varphi }_2^{}+\varphi _1^{}}& \lambda ^{\overline{\varphi }_2^{}+\varphi _2^{}}& \lambda ^{\overline{\varphi }_2^{}+\varphi _1^{}\mathrm{\Delta }}& \lambda ^{\overline{\varphi }_2^{}+\varphi _1^{}\mathrm{\Delta }}\\ \mathrm{𝟓}_{\overline{\mathrm{\Phi }}_1^{}}& \lambda ^{\overline{\varphi }_1^{}+c+\mathrm{\Delta }}& 0& \lambda ^{\overline{\varphi }_1^{}+\varphi }& \lambda ^{\overline{\varphi }_1^{}+\varphi _1^{}+\mathrm{\Delta }}& \lambda ^{\overline{\varphi }_1^{}+\varphi _2^{}+\mathrm{\Delta }}& \lambda ^{\overline{\varphi }_1^{}+\varphi _1^{}}& \lambda ^{\overline{\varphi }_1^{}+\varphi _2^{}}\\ \mathrm{𝟓}_{\overline{\mathrm{\Phi }}_2^{}}& \lambda ^{\overline{\varphi }_2^{}+c+\mathrm{\Delta }}& 0& \lambda ^{\overline{\varphi }_2^{}+\varphi }& \lambda ^{\overline{\varphi }_2^{}+\varphi _1^{}+\mathrm{\Delta }}& \lambda ^{\overline{\varphi }_2^{}+\varphi _2^{}+\mathrm{\Delta }}& \lambda ^{\overline{\varphi }_2^{}+\varphi _1^{}}& \lambda ^{\overline{\varphi }_2^{}+\varphi _2^{}}\end{array},$$
(4.2)
where $`\mathrm{\Delta }\frac{1}{2}(\overline{\varphi }\varphi \overline{c}+c)`$. The rank becomes seven for the charge assignment in Table I. On the other hand, the mass matrix for doublet Higgs becomes $`4\times 4`$ matrix. The charges in Table I lead to

$$M_D=\begin{array}{ccccccc}(\overline{L}^{})\backslash \overline{L}^{}& & C& \mathrm{\Phi }& \mathrm{\Phi }_1^{}& \mathrm{\Phi }_2^{}& \\ \overline{C}& (\mathrm{}& 0& 0& 0& 0& )\mathrm{}\\ \overline{\mathrm{\Phi }}& 0& 0& \lambda ^{\varphi _1^{}+\overline{\varphi }}& \lambda ^{\varphi _2^{}+\overline{\varphi }}\\ \overline{\mathrm{\Phi }}_1^{}& 0& \lambda ^{\varphi +\overline{\varphi }_1^{}}& \lambda ^{\varphi _1^{}+\overline{\varphi }_1^{}}& \lambda ^{\varphi _2^{}+\overline{\varphi }_1^{}}\\ \overline{\mathrm{\Phi }}_2^{}& 0& \lambda ^{\varphi +\overline{\varphi }_2^{}}& \lambda ^{\varphi _1^{}+\overline{\varphi }_2^{}}& \lambda ^{\varphi _2^{}+\overline{\varphi }_2^{}}\end{array}.$$
(4.3)
It is obvious that the rank is reduced to three, and therefore one pair of doublet Higgs appears in this model. The massless modes are written

$`H_u`$ $`=`$ $`(\overline{L}^{})_C,`$ (4.4)
$`H_d`$ $`=`$ $`(\overline{L}^{})_{\overline{C}}`$ (4.5)
where $`H_u`$ and $`H_d`$ are the doublet Higgs for up-quark sector and for down-quark sector, respectively.

### 4.2 Quark and lepton sector

In this subsection, we use the standard definition of $`\overline{\mathrm{𝟓}}(D^c,L)`$ field. If we introduce three generation matter fields $`\mathrm{\Psi }_i(\mathrm{𝟐𝟕})=\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}+\mathrm{𝟏𝟎}_{\mathrm{\Psi }_i}+\mathrm{𝟏}_{\mathrm{\Psi }_i}`$ $`(i=1,2,3)`$ with their charges $`(\psi _1,\psi _2,\psi _3)=(4,3,1)`$ in addition to the Higgs sector in Table I, the massless modes of $`\overline{\mathrm{𝟓}}`$ fields, where we have used the usual definition for $`\overline{\mathrm{𝟓}}`$, become

$`\overline{\mathrm{𝟓}}_1`$ $`=`$ $`\overline{\mathrm{𝟓}}_{}^{}{}_{\mathrm{\Psi }_1}{}^{}+\lambda ^3\overline{\mathrm{𝟓}}_{}^{}{}_{\mathrm{\Psi }_3}{}^{}+\lambda ^{1.5}\overline{\mathrm{𝟓}}_{\mathrm{\Psi }_2}+\lambda ^{3.5}\overline{\mathrm{𝟓}}_{\mathrm{\Psi }_3}`$
$`\overline{\mathrm{𝟓}}_2`$ $`=`$ $`\overline{\mathrm{𝟓}}_{\mathrm{\Psi }_1}+\lambda ^{2.5}\overline{\mathrm{𝟓}}_{}^{}{}_{\mathrm{\Psi }_3}{}^{}+\lambda ^1\overline{\mathrm{𝟓}}_{\mathrm{\Psi }_2}+\lambda ^3\overline{\mathrm{𝟓}}_{\mathrm{\Psi }_3}`$ (4.6)
$`\overline{\mathrm{𝟓}}_3`$ $`=`$ $`\overline{\mathrm{𝟓}}_{}^{}{}_{\mathrm{\Psi }_2}{}^{}+\lambda ^2\overline{\mathrm{𝟓}}_{}^{}{}_{\mathrm{\Psi }_3}{}^{}+\lambda ^{0.5}\overline{\mathrm{𝟓}}_{\mathrm{\Psi }_2}+\lambda ^{2.5}\overline{\mathrm{𝟓}}_{\mathrm{\Psi }_3},`$
where $`\overline{\mathrm{𝟓}}^{}(D^c,L^{})`$ and we fix the three bases of the massless modes $`(\overline{\mathrm{𝟓}}_1,\overline{\mathrm{𝟓}}_2,\overline{\mathrm{𝟓}}_3)`$ to $`(\overline{\mathrm{𝟓}}_{}^{}{}_{\mathrm{\Psi }_1}{}^{},\overline{\mathrm{𝟓}}_{\mathrm{\Psi }_1},\overline{\mathrm{𝟓}}_{}^{}{}_{\mathrm{\Psi }_2}{}^{})`$. These are obtained from the mass matrix of three $`\mathrm{𝟓}`$ fields and six $`\overline{\mathrm{𝟓}}`$ fields which are given from the interactions $`\mathrm{\Psi }_i\mathrm{\Psi }_j\mathrm{\Phi }Z`$ and $`\mathrm{\Psi }_i\mathrm{\Psi }_jC`$ by developing the VEVs of $`\mathrm{\Phi }`$, $`C`$ and $`Z`$. Then we can estimate the Yukawa couplings of quarks and leptons.

The Yukawa couplings of up quark sector are obtained as

$$Y_u=\begin{array}{cccccc}& & U_{\mathrm{\Psi }_1}^c& U_{\mathrm{\Psi }_2}^c& U_{\mathrm{\Psi }_3}^c& \\ Q_{\mathrm{\Psi }_1}& (\mathrm{}& \lambda ^6& \lambda ^5& \lambda ^3& )\mathrm{}\\ Q_{\mathrm{\Psi }_2}& \lambda ^5& \lambda ^4& \lambda ^2\\ Q_{\mathrm{\Psi }_3}& \lambda ^3& \lambda ^2& 1\end{array}$$
(4.7)
from the interactions $`\lambda ^{\psi _i+\psi _j+c}\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟎}_{\mathrm{\Psi }_j}C`$. The Yukawa couplings of down quark sector and of charged lepton sector are given as

$$Y_d^T(Y_e)=\begin{array}{cccccc}& & Q_{\mathrm{\Psi }_1}(E_{\mathrm{\Psi }_1}^c)& Q_{\mathrm{\Psi }_2}(E_{\mathrm{\Psi }_2}^c)& Q_{\mathrm{\Psi }_3}(E_{\mathrm{\Psi }_3}^c)& \\ \overline{\mathrm{𝟓}}_1& (\mathrm{}& \lambda ^6& \lambda ^5& \lambda ^3& )\mathrm{}\\ \overline{\mathrm{𝟓}}_2& \lambda ^{5.5}& \lambda ^{4.5}& 0\\ \overline{\mathrm{𝟓}}_3& \lambda ^5& \lambda ^4& \lambda ^2\end{array}$$
(4.8)
from the higher dimensional interactions $`\lambda ^{\psi _i+\psi _j+2\overline{c}}\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟔}_{\mathrm{\Psi }_j}\overline{C}\overline{C}`$ and $`\lambda ^{\psi _i+\psi _j+2\overline{c}}\mathrm{𝟏𝟎}_{\mathrm{\Psi }_i}\mathrm{𝟏}_{\mathrm{\Psi }_j}\overline{C}\overline{C}`$, respectively. Note that only $`\overline{\mathrm{𝟓}}^{}`$ fields can have non-vanishing Yukawa couplings through the interactions. This is because the interactions $`\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟔}_{\mathrm{\Psi }_j}\overline{C}\overline{\mathrm{\Phi }}`$ and $`\mathrm{𝟏𝟎}_{\mathrm{\Psi }_i}\mathrm{𝟏}_{\mathrm{\Psi }_j}\overline{C}\overline{\mathrm{\Phi }}`$ are forbidden by $`Z_2`$-parity. The above mass matrices give almost good values for masses and mixings for quark sector and charged lepton sector.

The Yukawa couplings for the Dirac neutrino are given as

$$Y_{n_D}=\begin{array}{ccccccccc}& & N_{\mathrm{\Psi }_1}^c& N_{\mathrm{\Psi }_2}^c& N_{\mathrm{\Psi }_3}^c& S_{\mathrm{\Psi }_1}& S_{\mathrm{\Psi }_2}& S_{\mathrm{\Psi }_3}& \\ \overline{\mathrm{𝟓}}_1& (\mathrm{}& \lambda ^{6.5}& \lambda ^{5.5}& \lambda ^{3.5}& \lambda ^6& \lambda ^5& \lambda ^3& )\mathrm{}\\ \overline{\mathrm{𝟓}}_2& \lambda ^6& \lambda ^5& \lambda ^3& \lambda ^{5.5}& \lambda ^{4.5}& \lambda ^{2.5}\\ \overline{\mathrm{𝟓}}_3& \lambda ^{5.5}& \lambda ^{4.5}& \lambda ^{2.5}& \lambda ^5& \lambda ^4& \lambda ^2\end{array}$$
(4.9)
through the interactions $`\lambda ^{\psi _i+\psi _j+c}\mathrm{𝟏𝟎}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟔}_{\mathrm{\Psi }_j}C`$. The vanishing component is caused by SUSY zero (holomorphic zero). The right-handed neutrino mass matrix becomes

$`M_{nR}`$ $`=`$ $`\begin{array}{cccccccccc}& & N_{\mathrm{\Psi }_1}^c& N_{\mathrm{\Psi }_2}^c& N_{\mathrm{\Psi }_3}^c& & S_{\mathrm{\Psi }_1}& S_{\mathrm{\Psi }_2}& S_{\mathrm{\Psi }_3}& \\ N_{\mathrm{\Psi }_1}^c& (\mathrm{}& \lambda ^8& \lambda ^7& \lambda ^5& & \lambda ^{7.5}& \lambda ^{6.5}& 0& )\mathrm{}\\ N_{\mathrm{\Psi }_2}^c& \lambda ^7& \lambda ^6& \lambda ^4& & \lambda ^{6.5}& 0& 0\\ N_{\mathrm{\Psi }_3}^c& \lambda ^5& \lambda ^4& 0& & 0& 0& 0\\ & \\ S_{\mathrm{\Psi }_1}& \lambda ^{7.5}& \lambda ^{6.5}& 0& & \lambda ^7& \lambda ^6& \lambda ^4\\ S_{\mathrm{\Psi }_2}& \lambda ^{6.5}& 0& 0& & \lambda ^6& \lambda ^5& \lambda ^3\\ S_{\mathrm{\Psi }_3}& 0& 0& 0& & \lambda ^4& \lambda ^3& \lambda \end{array}\mathrm{\Lambda }`$ (4.10)
through the interactions $`\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟔}_{\mathrm{\Psi }_j}\overline{C}\overline{C}`$, $`\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟔}_{\mathrm{\Psi }_j}\overline{C}\overline{\mathrm{\Phi }}Z`$ $`\mathrm{𝟏𝟔}_{\mathrm{\Psi }_i}\mathrm{𝟏𝟔}_{\mathrm{\Psi }_j}\overline{\mathrm{\Phi }}\overline{\mathrm{\Phi }}`$. Here vanishing components are caused by the SUSY zero (holomorphic zero) mechanism. Then the neutrino mass matrix is given by

$$M_\nu =Y_{n_D}M_{nR}^1Y_{n_D}^TH_u^2\eta ^2\lambda ^3\left(\begin{array}{ccc}\lambda ^2& \lambda ^{1.5}& \lambda \\ \lambda ^{1.5}& \lambda & \lambda ^{0.5}\\ \lambda & \lambda ^{0.5}& 1\end{array}\right)\frac{H_u^2\eta ^2}{\mathrm{\Lambda }},$$
(4.11)
where $`\eta `$ is a renormalization factor. This gives bi-large neutrino mixings but to realize the mass scale for the neutrino, we have to take the cutoff $`\mathrm{\Lambda }10^{13}`$ GeV if we take $`H_u\eta 200`$ GeV. Such a small cutoff scale leads to too short neucleon life-time via dimension six operators. Therefore, the charge assignment in Table I looks unrealistic.

However, because the neutrino scale is determined by the anomalous $`U(1)_A`$ charges as

$`M_\nu `$ $``$ $`\lambda ^{5l}\left(\begin{array}{ccc}\lambda ^2& \lambda ^{1.5}& \lambda \\ \lambda ^{1.5}& \lambda & \lambda ^{0.5}\\ \lambda & \lambda ^{0.5}& 1\end{array}\right){\displaystyle \frac{H_u^2\eta ^2}{\mathrm{\Lambda }}},`$ (4.12)
$`l`$ $`=`$ $`2c+\overline{c}10,`$ (4.13)
there may be other realistic models with other charge assignments. To obtain larger value of $`l`$, smaller $`c`$ and/or larger $`\overline{c}`$ is needed. Because $`C`$ includes $`H_u`$, the charge $`c`$ is determined as $`c=2\psi _3=2n`$ so that the top Yukawa coupling becomes $`O(1)`$. Here we take $`\psi _i=\delta _i+n`$ \[$`(\delta _1,\delta _2,\delta _3)=(3,2,0)`$\] to obtain realistic Cabbibo-Kobayashi-Maskawa matrix. To realize bi-large neutrino mixings (i.e., $`\overline{\mathrm{𝟓}}`$ fields in Eq. (4.6)), we must take

$$r\frac{1}{2}[(c\overline{c})(\varphi \overline{\varphi })]\frac{1}{2}.$$
(4.14)
Once we fix the mixing structure of $`\overline{\mathrm{𝟓}}`$ fields, the Yukawa couplings for down quarks are proportional to $`\lambda ^{\psi _i+\psi _j+2\overline{c}}\overline{C}\lambda ^{\delta _i+\delta _j+\frac{3}{2}(\overline{c}c)}`$. Therefore, roughly speaking, $`\mathrm{tan}\beta \frac{H_u}{H_d}`$ is proportional to $`\lambda ^{\frac{3}{2}(\overline{c}c)}`$. Then, smaller $`c`$ and/or larger $`\overline{c}`$ lead to smaller $`\mathrm{tan}\beta `$. For fixed $`\mathrm{tan}\beta `$, smaller $`c`$ and $`\overline{c}`$ lead to larger $`l`$. However, unless the condition

$$c2\overline{c}2$$
(4.15)
is satisfied, $`(Y_d)_{33}`$ vanishes by the SUSY zero mechanism. A compromised charge assignment is $`(\varphi ,\overline{\varphi },c,\overline{c},\varphi _i^{},\overline{\varphi }_i^{},\overline{z}_i,z,s^{})=(1,1,4,3,8,8,1,6,12)`$. Then $`l`$ becomes $`5`$, so the cutoff scale can be larger than the $`10^{15}`$ GeV. Actually, the running gauge couplings of $`SU(3)_C`$ and $`SU(2)_L`$, which should meet at the cutoff scale in this flipped $`SO(10)`$ scenario, meet around the scale in this charge assignment. And the Yukawa coupling of bottom quark becomes $`\lambda ^{3.5}`$ which can be realistic although the large ambiguity of $`O(1)`$ coefficients is required.

## 5 Summary

In this paper, we have shown that the missing partner mechanism in flipped $`SU(5)`$ model can be embedded in flipped $`SO(10)`$ model whose gauge group is $`SO(10)_F\times U(1)_{V_F^{}}E_6`$. It is interesting that the gauge group includes $`SU(2)_E`$, that plays an important role in solving SUSY flavor problem by the horizontal gauge symmetry and anomalous $`U(1)_A`$ gauge symmetry. As a proof of existence of a concrete model, we build a flipped $`SO(10)`$ model by introducing anomalous $`U(1)_A`$ gauge symmetry.

## Acknowledgement

N.M. is supported in part by Grants-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan.


format: text

text:
# A First-Principles Investigation of Phonon Softenings and Lattice Instabilities in the Shape-Memory System Ni2MnGa

## Abstract

Ferromagnetic Ni2MnGa has unique magnetoelastic properties. These are investigated by detailed computational studies of the phonon dispersion curves for the non-modulated cubic $`\mathrm{L2}_1`$ and tetragonal structures. For the $`\mathrm{L2}_1`$ structure, a complete softening of the transverse acoustic mode has been found around the wave vector $`𝐪=[1/3,1/3,0](2\pi /a)`$. The softening of this $`\mathrm{TA}_2`$ phonon mode leads to the premartensitic modulated super-structure observed experimentally. Further phonon anomalies, related to other structural transformations in Ni2MnGa, have also been found and examined. These anomalies appear to be due to the coupling of particular acoustic phonon modes and optical modes derived from Ni.

Displacive, diffusionless structural transformations of the martensitic type are known to occur in many metallic alloys. These transformations involve cooperative rather than diffusive displacements of atoms and are often associated with phonon anomalies in the parent phase and related precursor phenomena. It is a general challenge of the fundamental physics to explain driving forces of the martensitic transformations.

Since the discovery of a martensitic transformation in ferromagnetic Ni2MnGa Webster , this material has attracted strong interest Nature2 . This Heusler alloy is one of the very rare magnetic materials which undergo a martensitic transformation below the Curie temperature, whereby the combination of magnetic and structural features is responsible for its unique magnetomechanical properties. Shear deformations of more than 5 % have been obtained in magnetic fields Ullakko2 ; Murray ; Sozinov . These features together make Ni2MnGa very efficient for magnetic-shape-memory (MSM) technology Ullakko . The MSM technology is based on the magnetic field induced redistribution of martensitic domains in the sample. From the technological point of view, Ni2MnGa is much more promising than other materials being presently in commercial use, for example the well-known material Tb-Dy-Fe (Terfenol-D) which exhibits magnetostrictive strains of about 0.1 %. Design of new efficient MSM magneto-mechanical actuator devices is in progress Aaltio . Also the search for new materials with magnetic shape-memory effects is underway, even in antiferromagnets Nature1 .

Despite the experimental and technological success, a microscopic theory is missing, which would be able to show how far the premartensitic phase transition will result from the coupling of soft phonon modes and homogeneous strains associated with the shear constant $`C^{}`$ Gooding ; Brown . A complete understanding requires the evaluation of the phonon spectrum. The aim of this letter is to present supercell phonon calculations for the $`\mathrm{L2}_1`$ and $`\mathrm{T}`$ structures (see Fig. 1 and Table 1) using a state-of-the-art first-principles method based on density-functional theory.

Ni2MnGa is ferromagnetic at room temperature ($`T_\mathrm{C}380`$ K) and undergoes a two-step martensitic phase transformation. As a matter of fact, a number of different thermal and stress-induced martensitic structures have been observed in Ni2MnGa Mart-92 . Table 1 gives a summary of the structures found for Ni2MnGa Webster ; Brown ; Mart-92 . In the parent phase, several phonon anomalies can be observed with decreasing temperature. At the temperature $`T260`$ K, a nearly complete softening of the $`[1/3,1/3,0](2\pi /a)`$ transverse acoustic $`\mathrm{TA}_2`$ mode with polarization along \[11̄0\] leads to a premartensitic phase transition, which has been evidenced by neutron scattering Zhel95 ; Zhel ; Stuhr ; Brown , X-ray Fritsch , electron microscopy Cesari , and ultrasonic measurements Worgull ; Manosa97  or a combination of the previous methods Manosa . This structural transformation involves a commensurate periodic distortion of the parent phase with a propagation vector equal to that of the soft mode $`𝐪=[1/3,1/3,0](2\pi /a)`$ corresponding to six atomic planes or three lattice spacings. The precursor phenomena involves the magneto-elastic coupling as it has been described using Landau-type models Planes . These phenomenological models take as input very important lattice-dynamical properties. Thus, explicit calculation of the lattice dynamics is required since in transition metals, the d-electrons allow for many types of coupling, and several modes can be involved.

For the sake of simplicity, we present phonon calculations for the direction , which is the most interesting one as seen from the measurements for the acoustic modes Zhel . We have used the direct method for the evaluation of the phonon dispersion curves Parlinski ; Phonon , whereby the Hellmann-Feynman forces are calculated with the Vienna Ab-initio Simulation Package (VASP) Kresse-96 ; Kresse-99PAW  and the implemented projected augmented wave formalism (PAW) PAW . Within density-functional theory, the electronic exchange and correlation are treated by using the generalized gradient approximation (GGA). The 3d electrons of Ga have been included as valence states. The importance of using PAW for Ni2MnGa instead of pseudopotentials has been pointed out earlier Ayuela-priv . An orthorhombic supercell of 40 atoms formed by five tetragonal crystallographic unit cells, as seen in Fig. 1(b), has been used to calculate the phonon spectrum for the  direction detail1 . The direct method implies that we calculate forces induced on all atoms of the supercell when a single atom is displaced from its equilibrium position. Displacement of the atoms in only one tetragonal unit cell along all directions allows one to derive the force constant matrix and the dynamical matrix. Diagonalization of the dynamical matrix then leads to a set of eigenvalues for the phonon frequencies and corresponding eigenvectors. The calculations have been done by using an amplitude for the displacements of $`u=0.03`$ Å. This value is sufficient to calculate the forces with required precision and to fulfill the conditions of the harmonic approximation.

Phonon dispersion curves calculated for Ni2MnGa in the cubic direction are shown in Fig. 2(a). The acoustic modes are qualitatively and quantitatively in very good agreement with experimental results obtained from inelastic neutron scattering Zhel . For instance, following the acoustic branch, the values of the $`\mathrm{TA}_2`$ branch at $`\zeta =1`$, of the crossing between $`\mathrm{LA}`$ and $`\mathrm{TA}_1`$ branch, and of the maximum of the $`\mathrm{LA}`$ branch are 2.33, 4.43 and 5.13 THz, in good agreement with the experimental ones, 2.67, 4.35 and 5.07 THz Zhel . The initial slopes of the curves in the longitudinal ($`v_\mathrm{L}=5.23\times 10^5\mathrm{m}/\mathrm{s}`$) and transverse ($`v_{\mathrm{TA}_1}=3.35\times 10^5\mathrm{m}/\mathrm{s}`$; $`v_{\mathrm{TA}_2}=1.02\times 10^5\mathrm{m}/\mathrm{s}`$) modes agree well with the sound velocities measured via the neutron scattering dispersion curves Zhel . The slope of the phonon dispersion curves around the $`\mathrm{\Gamma }`$ point is positive, which is consistent with the stability around this phase.

We find a complete softening of the transverse acoustic $`\mathrm{TA}_2`$ phonon mode between $`\zeta =0.2`$ and $`\zeta =0.55`$. This softening occurs around the wave vector coordinate $`\zeta =1/3`$, which corresponds to the soft-mode phonon anomaly observed in the experimental studies Zhel . We point out that the experimentally observed premartensitic transition is related to this soft mode.

The softening of the $`\mathrm{TA}_2`$ mode means that the cubic $`\mathrm{L2}_1`$ high-temperature structure is unstable at zero temperature with respect to a particular atomic displacement leading to the formation of a modulated premartensitic structure. The $`3\mathrm{L}`$ structure can be modeled by a large supercell composed of three basic bct supercells in the  direction (see Fig. 4) detail2 . The atoms of the input structure have a modulated amplitude with a maximum around 0.05 Å where the Ni atoms are in opposite phase to the Mn and Ga atoms. Full relaxation of the three lattice parameters, $`a`$, $`b`$ and $`c`$ has been allowed. Also, lower magnetizations of 60-80 % of the full magnetization value for the whole supercell has been taken into account. In addition, two periods commensurate with the lattice were tested, but the structure with one period, as shown in Fig. 4, gives the minimum energy. The static displacements of the $`3\mathrm{L}`$ modulated structure are shown in Fig. 4(b). The modulation of the atoms stays in phase which differs from the input data, but it is typical for an acoustic branch such as $`\mathrm{TA}_2`$. The amplitudes for all the atoms are nearly the same, although they differ slightly but not according to the atom masses. Rather, the order is according to the sequence of the optical phonons energies, where Ni has the lowest energy, optical modes of Ga are in the next energy windows, and Mn optical phonons are of the highest energy (see Fig. 2).

We have considered so far only the acoustic mode $`\mathrm{TA}_2`$, because this mode is responsible for the phase transformation $`\mathrm{L2}_1`$ $``$ $`3\mathrm{L}`$. However, (see Fig. 2) the other two acoustic modes, $`\mathrm{TA}_1`$ and $`\mathrm{LA}`$, are coupled to the low - lying optical modes of Ga. Fig. 3 shows how this coupling affects the behaviour of the polarization vectors of Ga for its acoustic modes. The Ga atom type have been chosen for clarity, while acoustic modes of the other atoms show similar trends. Figure 3 illustrates that the acoustic-optical hybridization (compare with Fig. 2) leads to a number of anomalies, which occur at the corresponding wave vector coordinates $`\zeta =0.3`$, $`\zeta =0.43`$ and $`\zeta =0.715`$. A longitudinal-phonon optical mode can be associated to a charge density wave, which gives a sharp peak in the magnetic susceptibility $`\chi (q)`$. The crossing at $`\zeta =0.43`$ compares well with those for which anomalies are observed in the susceptibility calculations, which have been associated with the lattice instabilities in the $`5\mathrm{L}`$ and $`7\mathrm{L}`$ structures Velikokhatni ; Harmon . These crossings correlate in a natural way with the change of the $`C_{44}`$ elastic constant, which is, as it turns out, not related to the $`\mathrm{TA}_2`$ branch through the premartensitic transformation. It comes as a novelty that the involvement of the optical modes is required in order to understand the phase transformations in these alloys.

It is interesting to note that the anomaly in the $`\mathrm{TA}_2`$ branch has been found around the soft mode $`[1/3,1/3,0](2\pi /a)`$ in our calculations, which is the wave vector associated with the phonon anomaly of the premartensitic transition Zhel . Also the phonon branches on the  direction for the $`\mathrm{T}`$ structure are given in Fig. 2(b). In this latter case, the branches $`\mathrm{TA}_2`$ and $`\mathrm{TA}_1`$ become similar and the softening of the $`\mathrm{TA}_2`$ branch, present in the cubic phase, disappears. This finding is consistent with the stability of the $`\mathrm{T}`$ structure Ayuela  and adds insight into interpreting the experiments Manosa . The sequence of the optical modes is different in this case and they are not degenerate at the $`\mathrm{\Gamma }`$ and $`X`$ points. The instability of the cubic structure at zero temperature is the clue to understand how the martensite phase may nucleate. Instead of twinnings and stacking faults, the modulation of the martensitic tetragonal structure of Ni2MnGa corresponds to a smaller structural change easier to accommodate or to derive from the parent cubic structure. With respect to the electronic structure, the lowering of the symmetry allows for a splitting of the density of states peak at the Fermi level Brown ; Ayuela , thereby decreasing the valence band energy contribution to the total energy. Thus, structural modulation and lowering of the electronic energy explains the 5L structure of Ni2MnGa Za .

On the other hand, the recent susceptibility calculations Velikokhatni ; Harmon  have shown that the peak in the generalized susceptibility at $`𝐪=[\zeta ,\zeta ,0](2\pi /a)`$ for $`\zeta =0.33`$ corresponds to a magnetization of about 60%. This suceptibility peak moves toward higher values of $`\zeta `$ as the magnetization increases: full magnetization leads to a peak at $`\zeta =0.4`$. However our calculations show no need for the extra magnetization hypothesis, since all the interesting $`𝐪`$ vectors are already present. This casts serious doubts on the role of magnetization in these intermediate or premartensitic transitions. Secondly, this finding reinforces the arguments that the elastic entropy term drives the transition apl  and the spin spirals do not show anomalies involved in the shuffling of the $`\mathrm{L2}_1`$ structure prb . In view of these results, the intermediate and the martensitic transition seem to be really independent, although, of course, the $`\mathrm{L2}_1`$ structure is necessary in order to have the intermediate transformation at $`\zeta =0.33`$.

In summary, we have calculated the phonon dispersion for two structures of Ni2MnGa: the cubic $`\mathrm{L2}_1`$ and the tetragonal $`\mathrm{T}`$ at zero temperature. The cubic structure is found to be unstable with respect to a specific re-arrangement (shuffling) of the atoms which leads to a modulated super-structure $`3\mathrm{L}`$. We have calculated stability of the $`3\mathrm{L}`$ structure with high accuracy and found that it reflects perfectly our phonon calculations. Furthermore, we have shown that the subtle competition of low-lying optical modes of Ni and acoustic modes account for the experimentally observed wave vectors, for which anomalies have been observed ($`3\mathrm{L}`$, $`5\mathrm{L}`$ and $`7\mathrm{L}`$ structures). These calculations have been done by using first-principles methods, and are in very good agreement with experimental findings.

###### Acknowledgements.

This work has been supported by the Graduate School “Structure and Dynamics of Heterogeneous Systems” of the Deutsche Forschungsgemeinschaft (DFG) and by the Academy of Finland through its Centers of Excellence Program (2000-2005). Computer facilities of the Research Centre Jülich (Project ”Lattice dynamics”) are acknowledged. We thank Prof. M. Acet, Prof. K. Parlinski, Dr. A. Postnikov and Prof. A. Planes for valuable discussions.


text:
# Nonequilibrium Processes in Non-Hamiltonian Statistical Ensembles

(November 10, 2003)

## Abstract

A nonequilibrium statistical operator method is developed for ensembles of particles obeying non-Hamiltonian equations of motion in classical phase space. The main consequences of non-zero compressibility of phase space are examined in terms of time-dependent dynamic quantities. The generalized transport equations involve the phase-space compressibility in a non-trivial way. Our results are useful in molecular dynamics simulation studies as well as nonequilibrium or quasiclassical approximations of quantum-classical dynamics.

preprint: Submitted to Journal of Mathematical Physics

Theoretical interest in the non-Hamiltonian motion appears in the early works Ashurst and Hoover (1973); Andersen (1980); Evans (1983), where it was realized that the constant-temperature simulations require the introduction of additional non-Hamiltonian acceleration into dynamical equations. This has opened new research directions associated with extended systems. For the past few years the physics of non-Hamiltonian systems was partially developed for molecular dynamics simulation purposes Nosé (1984); Hoover (1985); Martyna et al. (1992); Cao and Voth (1996); Yang and Cao (2002). Particular applications of non-Hamiltonian dynamical systems are overviewed in a recent paper Tuckerman et al. (2001), where the authors exploit their formulation of non-Hamiltonian statistical mechanics Tuckerman et al. (1999). In this regard it is important to mention also the interesting studies of mixed quantum-classical systems Tully (1998); Burant and Tully (2000); Kapral and Ciccotti (1999); Nielsen et al. (2000a, b), where the non-Hamiltonian behavior plays an essential role.

One of the most popular non-Hamiltonian ensembles is the so-called Nosé-Hoover chain Martyna et al. (1992). In case with two thermostat variables, the corresponding Hamiltonian is

$$_{NH}=\frac{p^2}{2m}+\frac{p_{\eta _1}^2}{2M_{\eta _1}}+\frac{p_{\eta _2}^2}{2M_{\eta _2}}+V(q)+gk_bT(\eta _1+\eta _2),$$
(1)
where $`\eta _i`$ and $`p_{\eta _i}`$ are the conjugate thermostat variables, $`g`$ is the number of degrees of freedom, $`T`$ is the thermostat temperature. It can be shown that the Hamiltonian (1) is conserved by the following (non-Hamiltonian) equations of motion:

$$\dot{q}=\frac{p}{m},$$
(2)
$$\dot{\eta }_i=\frac{p_{\eta _i}}{M_{\eta _i}},$$
(3)
$$\dot{p}=\frac{p_{\eta _1}}{M_{\eta _1}}p\frac{V(q)}{q},$$
(4)
$$\dot{p}_{\eta _1}=\frac{p^2}{m}gk_bTp_{\eta _1}\frac{p_{\eta _2}}{M_{\eta _2}},$$
(5)
$$\dot{p}_{\eta _2}=\frac{p_{\eta _1}^2}{M_{\eta _1}}gk_bT.$$
(6)
In real MD simulations the particular choice of non-Hamiltonian ensemble can vary depending on the physical problem under study. But all of them can be described in general formalism we present here. Recently using a generalization of the symplectic form of the Hamilton equations of motion the authors of Ref. Sergi and Ferrario, 2001 showed that there is a unique general structure that underlies the majority of the equations of motion for extended systems. Furthermore, Sergi Sergi (2003, 2004) proposed a consistent theory of statistical mechanics for non-Hamiltonian systems under equilibrium conditions. The tempting feature of the theory Sergi and Ferrario (2001); Sergi (2003, 2004) is that it is based on the introduction of a suitable Poisson bracket, and has the general structure of classical formalism Morrison (1998) with a natural appearance of the nonzero phase-space compressibility. The main feature of the systems of interest is energy conservation, which is a reason of non-Hamiltonian behavior. It is interesting therefore to consider the situation, when the energy is still conserved, but the system itself is out of equilibrium, particularly when the dynamical quantities depend on time explicitly.

The starting point in describing the non-Hamiltonian statistical ensembles is the introduction of the proper non-Hamiltonian Poisson bracket:

$$\{A(t),C(t)\}=\underset{ik}{\overset{2N}{}}\frac{A}{x_i}_{ik}\frac{C}{x_k},$$
(7)
where $`𝐱=(q,p)=(q_1,q_2\mathrm{}q_N;p_1,p_2\mathrm{}p_N)`$ is the $`2N`$-dimensional vector in phase space, $`A(t)`$ and $`C(t)`$ are arbitrary dynamical quantities (we allow them to depend on time explicitly), $`_{ik}=_{ki}`$ is an antisymmetric matrix, whose elements are functions of $`𝐱`$ and, possibly, time. Following Ref. Sergi, 2003 we write the equations of motion in standard form $`\dot{x}_j=\{x_j,\}`$, where $``$ is the Hamiltonian, which is a constant of motion by construction due to of antisymmetry of the matrix $`_{ik}`$.

The main feature of the non-Hamiltonian system is the failure of the Jacobi relation, which implies time-irreversibility of non-Hamiltonian mechanics. This is formally a consequence of our requirement for eq. (7) and the standard equation of motion to be satisfied simultaneously. By defining the Liouville operator

$$i\widehat{L}A=_{ij}\frac{A}{x_i}\frac{}{x_j}=\{A,\},$$

we obtain the generalized Liouville equation (see Ref. Sergi and Ferrario, 2001) in the form

$$\frac{}{t}\rho (q,p;t)=\left(i\widehat{L}+\varkappa \right)\rho (q,p;t),$$
(8)
where

$$\varkappa =\underset{ij}{\overset{2N}{}}\frac{_{ij}}{x_i}\frac{}{x_j}$$

is the phase space compressibility. For the purpose of generality we also allow it to depend on time explicitly (through $`_{ij}`$). For particular determination of the matrix $`_{ij}`$ (and therefore $`\varkappa `$) it is useful to write the equation of motion in the form

$$\dot{x}_i=\underset{j=1}{\overset{2N}{}}_{ij}\frac{}{x_j}.$$

For the particular ensemble, described by Eqs. (1) - (6), we obtain

$$=\left(\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 1& 0& 0& 0& p& 0\\ 0& 1& 0& p& 0& p_{\eta _1}\\ 0& 0& 1& 0& p_{\eta _1}& 0\end{array}\right),$$

so that the compressibility can be calculated as

$$\varkappa =\frac{p_{\eta _1}}{M_{\eta _1}}\frac{p_{\eta _2}}{M_{\eta _2}}.$$

In general for a system of $`N`$ identical particles (more exactly for $`N`$ degrees of freedom) we define the statistical time-dependent average as

$$A^t=\underset{N}{}A_N(q,p;t)\rho _N(q,p;t)𝑑\mathrm{\Gamma }_N\mathrm{Tr}\{A(q,p;t)\rho (q,p;t)\},$$
(9)
where $`A_N(q,p;t)`$ denotes the dynamic variable $`A`$ for a system of $`N`$ particles, i.e., $`A(q_1,\mathrm{},q_N;p_1,\mathrm{},p_N;t)`$, $`d\mathrm{\Gamma }_N`$ is the element of phase volume.

The time evolution of an arbitrary dynamical variable obeys the usual equation

$$\frac{dA}{dt}=\frac{A}{t}+\{A,\}=\frac{A}{t}+i\widehat{L}A.$$
(10)
One interesting consequence of the time-irreversibility is that the average derivative of a dynamical variable is not necessarily equal to the derivative of its average. In fact, taking into account eq. (8) we obtain the relation

$$\frac{dA^t}{dt}=\left(\frac{A}{t}+\{A,\}\varkappa A\right)\rho 𝑑\mathrm{\Gamma },$$
(11)
or, alternatively

$$\frac{dA}{dt}^t=\frac{dA^t}{t}+\varkappa A^t.$$
(12)
Equations (11) and (12) are particularly important in identifying the real measured mean quantity. Another important peculiarity of non-Hamiltonian systems follows if one is interested in the time evolution operator for the initial statistical distribution $`\rho (q,p;t)=U_\pm ^{(\varkappa )}(t,t^{})\rho (q,p;0)`$. Following the standard procedure Zubarev (1971) we obtain the time-ordered evolution operator $`U_\pm ^{(\varkappa )}(t,t^{})`$ in the following form

$$U_\pm ^{(\varkappa )}(t,t^{})=\mathrm{exp}_\pm \left\{\underset{t^{}}{\overset{t}{}}\left(i\widehat{L}(\tau )+\varkappa \right)𝑑\tau \right\}$$
(13)
where $`\mathrm{exp}_\pm `$ denote the correspondingly time-ordered exponents. We emphasize that the evolution operator in Eq. (13) does not describe the evolution of dynamic variables. Eq. (10) gives $`A(q,p;t)=U_\pm ^{}(t,t^{})A(q,p;t^{})`$, where

$$U_\pm ^{}(t,t^{})=U_\pm ^{(0)}(t,t^{})=\mathrm{exp}_\pm \left\{i\underset{t^{}}{\overset{t}{}}\widehat{L}(\tau )𝑑\tau \right\}.$$
(14)
Let us suppose that the given nonequilibrium non-Hamiltonian ensemble can be characterized by the average values $`\zeta _n^t`$ of a finite set of dynamical variables $`\{\zeta _n(t)\}`$. Then the time evolution of the nonequilibrium macroscopic state is described by the generalized transport equations (see Eqs. (11) and (12))

$$\frac{\zeta _n^t}{t}=\mathrm{Tr}\left\{\left(\dot{\zeta }_n\varkappa \zeta _n\right)\rho (t)\right\}=\mathrm{Tr}\left\{\left(\{\zeta _n,\}\varkappa \zeta _n\right)\rho (t)\right\},$$
(15)
where the trace operation is defined in Eq. (9). Let $`\stackrel{~}{\rho }(t)`$ be a trial statistical operator, such that $`\zeta _n^t=\mathrm{Tr}\left(\zeta _n\stackrel{~}{\rho }(t)\right)`$, and $`\mathrm{Tr}\left(\stackrel{~}{\rho }(t)\right)=1`$. We choose the entropy functional in the standard form Zubarev (1971),

$$\stackrel{~}{S}[\stackrel{~}{\rho }(t)]=\mathrm{Tr}\left(\stackrel{~}{\rho }(t)\mathrm{ln}\stackrel{~}{\rho }(t)\right)\underset{n}{}\phi _n(t)\mathrm{Tr}\left(\zeta _n\stackrel{~}{\rho }(t)\right)\lambda (t)\mathrm{Tr}\stackrel{~}{\rho }(t),$$
(16)
where $`\phi _n(t)`$ and $`\lambda (t)`$ are the corresponding Lagrange multipliers. Varying (16) with respect to the trial distribution $`\stackrel{~}{\rho }(t)`$ we obtain the relevant distribution,

$$\rho _r=\mathrm{exp}\left\{\mathrm{\Phi }(t)\underset{n}{}\phi _n(t)\zeta _n\right\},$$
(17)
where $`\mathrm{\Phi }(t)`$ is the Massieu-Plank function. Here we use the notion “relevant” following Zubarev Zubarev (1971); Zubarev et al. (1996); Zubarev (1996) for the distribution that maximizes the entropy functional (16) at a particular time. The relevant distribution is not yet the required nonequilibrium distribution because, in general, it does not satisfy the Liouville equation. Nevertheless, as shown below the relevant statistical distribution will serve as an auxiliary distribution to determine special solutions of the Liouville equation that describe irreversible macroscopic processes.

The Lagrange multipliers $`\phi _n`$ are determined from the self-consistency conditions,

$$\zeta _n^t=\zeta _n_r^t\mathrm{Tr}\left(\zeta _n\rho _r(t)\right),$$
(18)
so that $`S(t)=\mathrm{\Phi }(t)+_n\phi _n(t)\zeta _n^t`$. The meaning of the nonequilibrium entropy is the same as the corresponding entropy in the Hamiltonian systems Zubarev (1971); Zubarev et al. (1996). In contrast with the Gibbs entropy, $`S(t)`$ can change in time. Let us demonstrate the specific feature of this entropy in the non-Hamiltonian system by calculating the derivative

$$\frac{dS(t)}{dt}=\underset{n}{}\frac{\delta S(t)}{\delta \zeta _n^t}\frac{\zeta _n^t}{t}=\underset{n}{}\phi _n\frac{\zeta _n^t}{t}=\underset{n}{}\phi _n\mathrm{Tr}\left\{\left(\dot{\zeta }_n\varkappa \zeta _n\right)\rho (t)\right\}.$$
(19)
It seems confusing that $`dS(t)/dt0`$ even for constant $`\zeta _n`$ due to the compressibility of the phase space. This of course is misunderstanding, because the entropy and its evolution are essentially determined by the mean values of $`\zeta _n`$ and by the derivatives of their mean values. This is particularly important in applications of the theory since the interpretations of the desired quantities depend on the definition of the averaging procedure.

To find the formal solution of the generalized Liouville equation (8), we assume that at some initial instant of time the distribution coincides with the relevant one (17), satisfying the self-consistency relations (18). Following the method proposed by Zubarev Zubarev (1996), the retarded solutions of the Liouville equation in our case is found by solving the following perturbed equation

$$\frac{\rho (t)}{t}+\left(i\widehat{L}+\varkappa \right)\rho (t)=ϵ\left(\rho (t)\rho _r(t)\right)$$
(20)
with $`ϵ0`$ in the final expressions. The formal solution of the initial Liouville equation (8) reads

$$\rho (t)=\underset{ϵ+0}{lim}ϵ\underset{\mathrm{}}{\overset{t}{}}e^{(ϵ+\varkappa )(tt^{})}e^{i(tt^{})\widehat{L}}\rho _r(t^{})𝑑t^{}.$$
(21)
Integrating eq. (21) by parts we arrive at

$$\rho (t)=\rho _r(t)\underset{ϵ+0}{lim}ϵ\underset{\mathrm{}}{\overset{t}{}}e^{(ϵ+\varkappa )(tt^{})}e^{i(tt^{})\widehat{L}}\left(\frac{}{t^{}}+i\widehat{L}+\varkappa \right)\rho _r(t^{})𝑑t^{}.$$
(22)
In principle equation (22) answers the question how to construct nonequilibrium distribution, which is a functional of dynamic variables and the phase space compressibility, if the system develops from a relevant state.

Let us further write the statistical operator as a sum of relevant distribution and local-time deviation $`\rho (t)=\rho _r(t)+\mathrm{\Delta }\rho (t)`$, so that equation (20) becomes

$$\left(\frac{}{t}+i\widehat{L}+\varkappa +ϵ\right)\mathrm{\Delta }\rho (t)=\left(\frac{}{t}+i\widehat{L}+\varkappa \right)\rho _r(t).$$
(23)
Then, keeping in mind that $`\rho _r(t)`$ depends on time only through the relevant variables $`\zeta _n`$, we can calculate its time derivative

$`{\displaystyle \frac{\rho _r(t)}{t}}={\displaystyle \underset{n}{}}{\displaystyle \frac{\delta \rho _r(t)}{\delta \zeta _n^t}}{\displaystyle \frac{\zeta _n^t}{t}}={\displaystyle \underset{n}{}}{\displaystyle \frac{\delta \rho _r(t)}{\delta \zeta _n^t}}\mathrm{Tr}\left\{\left(\dot{\zeta }_n\varkappa \zeta _n\right)\rho (t)\right\}`$
$`={\displaystyle \underset{n}{}}{\displaystyle \frac{\delta \rho _r(t)}{\delta \zeta _n^t}}\mathrm{Tr}\left\{\zeta _n{\displaystyle \frac{\rho (t)}{t}}\right\}={\displaystyle \underset{n}{}}{\displaystyle \frac{\delta \rho _r(t)}{\delta \zeta _n^t}}\mathrm{Tr}\left\{\zeta _n\left(i\widehat{L}+\varkappa \right)\rho (t)\right\}.`$ (24)
We define the Kawasaki-Gunton operator Kawasaki and Gunton (1973) in the standard manner

$$𝒫_{KG}A=\rho _r(t)\mathrm{Tr}A+\underset{n}{}\frac{\delta \rho _r(t)}{\delta \zeta _n^t}\left\{\mathrm{Tr}(A\zeta _n)\zeta _n^t\mathrm{Tr}A\right\}.$$
(25)
Using the relations $`\mathrm{Tr}\{i\widehat{L}\rho (t)\}=0`$, $`\mathrm{Tr}\left\{\dot{\rho }(t)\right\}=\varkappa ^t`$, one can prove that the projection operator (21) conserves its features in the non-Hamiltonian case, namely $`𝒫_{KG}(t)\rho (t)=\rho _r(t)`$, $`𝒫_{KG}(t)\dot{\rho }(t)=\dot{\rho }_r(t)`$. Equation (24) can now be rewritten in the projected form

$$\frac{\rho _r(t)}{t}=𝒫_{KG}(t)\left(i\widehat{L}+\varkappa \right)\rho (t).$$
(26)
Inserting equation (26) into Eq. (23) and integrating over time, we obtain a solution for $`\mathrm{\Delta }\rho (t)`$ as follows

$$\mathrm{\Delta }\rho (t)=\underset{\mathrm{}}{\overset{t}{}}e^{ϵ(tt^{})}U_r^{(\varkappa )}(t,t^{})𝒬_{KG}(t^{})\left(i\widehat{L}+\varkappa \right)\rho _r(t^{})𝑑t^{},$$
(27)
where $`𝒬_{KG}=1𝒫_{KG}`$, and (compare with Eq. (13))

$$U_r^{(\varkappa )}(t,t^{})=\mathrm{exp}\left\{\underset{t^{}}{\overset{t}{}}𝒬_{KG}(t)\left(i\widehat{L}+\varkappa \right)𝑑\tau \right\}.$$
(28)
Equation (27) allows us to express the right-hand side of Eq. (15) in terms of relevant distribution, i.e. in terms of variables $`\zeta _n^t`$. So we get the closed set of generalized transport equations for observables.

To obtain more transparent equations we define the entropy operator

$$\widehat{S}(t)=\mathrm{\Phi }(t)+\underset{n}{}\phi _n(t)\zeta _n(t),$$

so that $`\rho _r(t)=\mathrm{exp}\{\widehat{S}(t)\}`$. Introducing the Mori projection operator Mori (1965)

$$𝒫_M(t)A=A_r^t+\underset{n}{}\frac{\delta A_r^t}{\delta \zeta _n^t}\left(\zeta _n\zeta _n^t\right),$$
(29)
which satisfy the identities $`𝒫_M^2(t)=𝒫_M(t)`$, $`𝒫_M(t)\zeta _n=\zeta _n`$, one can prove the relation

$$𝒬_{KG}(t)\left\{(i\widehat{L}+\varkappa )\rho _r\right\}=𝒬_M(t)\left\{\varkappa i\widehat{L}\widehat{S}\right\}\rho _r,$$
(30)
where $`𝒬_M(t)=1𝒫_M(t)`$ is complementary to the Mori projector (29). Now the generalized transport equation (15) reads

$$\frac{}{t}\zeta _n^t=\frac{}{t}\zeta _n_r^t+\mathrm{Tr}\left\{𝒬_M(t)\left(\dot{\zeta }_n\varkappa \zeta _n\right)\mathrm{\Delta }\rho \right\},$$
(31)
where $`\mathrm{\Delta }\rho `$ must be taken from equation (27). It is useful to define new dynamical variables (generalized fluxes)

$$_n(t)=𝒬_M(t)\left(\dot{\zeta }_n\varkappa \zeta _n\right),$$

which are orthogonal to the space of relevant variables in a sense that $`𝒫_M_n(t)=0`$. The variables $`_n(t)`$ determine the effects of macroscopic degrees of freedom of the evolution on the nonequilibrium state of the non-Hamiltonian system described by a set of relevant variables and a given phase space compressibility $`\varkappa `$.

Incorporating equations (27), (30), and (31) we obtain the final result for the generalized transport equations

$$\frac{}{t}\zeta _n^t=\frac{}{t}\zeta _n_r^t+\underset{n}{}\underset{\mathrm{}}{\overset{t}{}}e^{ϵ(tt^{})}𝒟_{mn}(t,t^{})\phi _n(t^{})𝑑t^{}+\underset{\mathrm{}}{\overset{t}{}}e^{ϵ(tt^{})}_{mn}(t,t^{})𝑑t^{},$$
(32)
where

$$𝒟_{mn}(t,t^{})=\mathrm{Tr}\left\{_m(t)U_r^{(\varkappa )}(t,t^{})_m(t^{})\rho _r(t^{})\right\}$$
(33)
are the normal generalized kinetic coefficients, and

$$_{mn}(t,t^{})=\mathrm{Tr}\left\{_m(t)U_r^{(\varkappa )}(t,t^{})𝒬_M(t^{})\varkappa \left(\underset{n}{}\phi _n\zeta _n1\right)\rho _r(t^{})\right\}$$
(34)
are the effective forces determining the anomalous contribution into the system evolution associated with nonlinear phase space compressibility.

The last term in Eq. (32) is associated with the initial irreversibility of the process due to non-Hamiltonian behavior, while the second term is due to the nonequilibrium process itself. It does not vanish even when $`\varkappa =0`$, and coincides with the well-known result for Hamiltonian ensembles. These terms demonstrate the general structure of the solution and clearly show the contributions from two physically different processes. In practice equations (33) and (34) can be used with standard simplifications for particular systems.

In conclusion, we have proposed a possible realization of the nonequilibrium statistical mechanics for ensembles of particles obeying non-Hamiltonian dynamics. Starting with the non-Hamiltonian Poisson bracket (7) and the equations of motion in standard form, we derived the basic relations for time-dependent dynamical variables with conserved Hamiltonian. Introducing the reduced description method in analogy with Hamiltonian dynamics we found the retarded solution (22) of the generalized Liouville equation (8). We emphasize however, that Eq. (22) is the specific solution, not general, and it describes the evolution from a specific auxiliary relevant distribution (this is reflected in the introduction of parameter $`ϵ`$). Using the projection operator technique we obtained the generalized transport equations (32), which contain two physically different contributions due to the statistical irreversibility (33) and the time-irreversibility (34) associated with non-Hamiltonian features of the system.

###### Acknowledgements.

The research is supported by the NSF career award and the ACS petroleum research fund.


text:
# Equation of state and phonon frequency calculations of diamond at high pressures.

## I Introduction

Diamond is the archetype of the covalently bonded, tetrahedrally coordinated insulators. It’s extreme hardness is highly valued in technology and is also exploited in high pressure research when using the diamond anvil cell. The elastic properties of diamond near ambient conditions are well characterized through ultrasonic and Brillouin techniques. However, as a consequence of the small compressibility, the changes in elastic properties under high hydrostatic pressure are not well confined experimentally. This, for instance, applies to the variation of the bulk modulus $`B`$ with pressure $`P`$, $`B_0^{}=(dB/dP)_{P=0}`$, a basic parameter in the equation-of-state (EOS) modelling. A property closely related to the compression behavior of diamond is the pressure-induced frequency shift of the threefold degenerate F2g zone-center optical phonon mode; its shift with pressure provides an approximate measure of the change in relative density, because the mode Grüneisen parameter is close to oneMitra69 ; Whalley76 ; Grims78 ; Hanfl85b .

High-pressure x-ray diffraction experimentsAleks87 ; Gille99 ; Occelli02  yield the ambient-pressure bulk modulus $`B_0`$ in good agreement with acoustic measurementsMcSki72 ; Grims75 ; Ramda93 ; Vogel96 ; Zoubo98 . The analysis of diffraction data is usually based on adopting $`B_0^{}=4.0`$ obtained from ultrasonic measurements up to 0.2 GPaMcSki72 . The only exception is the recent diffraction study of diamond to 140 GPa by Occelli et al.Occelli02 , who report $`B_0^{}=3.0`$, at variance with the ultrasonic measurement. It was subsequently arguedHolza02  that the ruby pressure calibrationMao86  employed in the diffraction work of Ref. \[Occelli02, \] may need a revision. If so, the pressure shift of the F2g phonon frequency, also studied by Occelli et al. up to 140 GPa, should be affected in a similar manner. The F2g phonon mode behavior at high pressures was frequently studied by Raman spectroscopyMitra69 ; Whalley76 ; Grims78 ; Hanfl85b ; Boppa85 ; Aleks87 ; Gonch85 ; Aleks86 ; Tardi90 ; Muino94 ; Schif97  and its possible role in pressure calibration was addressed early onHanfl85b ; Aleks87 ; Sherm85 .

We report the calculation of the EOS and optical phonon frequency of diamond at high pressures within Density Functional Theory (DFT). Extensive theoretical work on diamond under pressure by ab initio methods has addressed changes in bonding, elasticity, lattice dynamics, thermodynamical properties, phase stability, and electronic excitations, see, e.g., Refs. \[Zunger77, ; Yin83, ; Cheli84, ; Cohen85, ; Cheli87, ; Hanfl85b, ; Niels86, ; Cardo86, ; Fahy87, ; Fahy87a, ; Fahy87b, ; Fahy88, ; Mailh91, ; Vanca92, ; Pavon93, ; Willa94, ; Clark95, ; Scand96, ; Serra98, ; Xie99, ; Wu99, ; Zhao99, \]. Here, we are interested in a specific question: What are, based on different implementations and approximations of DFT, the constraints on $`B_0^{}`$ of diamond and on the nonlinear pressure shift of the optical phonon frequency. The calculated values of $`B_0^{}`$ reported in the literature scatter by about 25%, spanning a range similar to that of the experimental results. Calculated pressure effects on the optical phonon frequency were reported in Refs. \[Hanfl85b, \] and \[Niels86, \]. Results of a more recent calculationWu99  show some disagreement with Refs. \[Hanfl85b, \] and \[Niels86, \]. In view of new experimental resultsOccelli02  and the – at least apparent – uncertainties in the previous theoretical predictions we considered it worthwhile to revisit the calculation of the EOS and optical phonon frequency of diamond under pressure, combined with accurate procedures to extract the parameters of interest.

## II Details of the calculations

The total energy calculations performed in this work are based on DFTHohenberg64  within, on one hand, the Local Density ApproximationKohn65  (LDA) and, on the other hand, the Generalized Gradient ApproximationPW92 ; PBE96  (GGA). In order to grasp the uncertainties consequent to the choice of the computational method and of its inherent assumptions, we are using simultaneously three different approaches, viz. the pseudopotential method applied in the basis of plane wavesGonze02 ; Goe97 ; PTAA92 ; Gon96  (PW), an all-electron method relying on the Augmented Plane Waves plus Local Orbitalssoft:WIEN2k ; SNS00 ; MBSS01  (APW+lo), and an intermediate approach implemented in the basis of Projector Augmented Wavessoft:VASP-short ; Kress99  (PAW); the latter approachBloch94  treats the valence states as part of an all-electron problem and describes them by all-electron wave-functions.

The pseudopotential employed is the “dual-space” separable pseudopotential of Hartwigsen, Goedecker, and HutterHartw98  and its GGA counterpart constructed by X. GonzeGonze02 ; Gonze ; the PAW potentials were constructed by G. Kresse and denoted as C\_h in Ref. \[soft:VASP-short, \]. All methods take into account scalar-relativistic corrections (though no meaningful contribution expected in diamond), either explicitly (APW+lo) or through the construction of the (pseudo-)potentials.

The numerical convergence of all three methods, with respect to the size of the basis set and k-space sampling was thoroughly tested. The plane-wave cutoffs of 75 Hartree and 1500 eV were applied with, respectively, the PW and the PAW-basis. These rather large values are required for getting reliable results for the pressure which is calculated analytically, in both approaches, with the aid of the stress theoremNielsen83 ; Nielsen85 . In the APW+lo method a muffin-tin radius of $`R_{MT}=1.2`$ a.u. was used. The plane-wave expansion was limited by $`R_{MT}\times K_{max}=9`$ and the charge density Fourier expansion by $`G_{max}=20`$ Ry1/2. In the APW+lo method the $`P(V)`$ results were generated from $`E^{\text{tot}}(V)`$ fitted as described below.

The integration over the Brillouin zone is performed using the same set of 28 k-points in the irreducible wedge in both plane-wave based approaches; the set was generated by the “special points” approach Monkh76  using a $`6\times 6\times 6`$ mesh with 4 different fractional shifts. In the APW+lo approach the tetrahedron integration was employedJA71 ; BJA94 , based on a uniform, $`\mathrm{\Gamma }`$-centered $`10\times 10\times 10`$ mesh of 47 k-points.

The frozen-phonon approach, which consists in the evaluation of the total energy $`E^{\mathrm{tot}}`$ of the crystal with frozen-in atomic displacements, was applied (in the two plane-wave based methods only) with the same cutoffs and with the same k-point mesh ($`6\times 6\times 6`$ \+ 4 shifts) which, due to the lowered symmetry, results in a 91 k-points set. Small displacements $`\stackrel{}{u}(1)=\pm (u,u,u)`$ and $`\stackrel{}{u}(2)=(u,u,u)`$ with $`u/a`$ = 0.002 were applied to the two atoms of the basis – thus either compressing or stretching the C(1) – C(2) bond, like in the zone-center F2g mode – and the two values obtained for $`\mathrm{\Delta }E^{\mathrm{tot}}`$ \[viz. the $`\mathrm{\Delta }E^{\mathrm{tot}}`$(outward) and $`\mathrm{\Delta }E^{\mathrm{tot}}`$(inward)\] were averaged. The eigenfrequency is then found from the expression for the energy of a harmonic oscillator

$$\mathrm{\Delta }E^{\mathrm{tot}}=\frac{1}{2}M\omega ^2\left[|\stackrel{}{u}(1)|^2+|\stackrel{}{u}(2)|^2\right],$$
(1)
where $`M`$ is the atomic mass. The above-mentioned averaging procedure of the $`\mathrm{\Delta }E^{\mathrm{tot}}`$(outward) and $`\mathrm{\Delta }E^{\mathrm{tot}}`$(inward) eliminates the cubic contribution to $`\mathrm{\Delta }E^{\mathrm{tot}}`$. It turns out that the remaining quartic anharmonicity is small at $`u/a=0.002`$ and contributes to the uncertainty of the resulting eigenfrequencies by less than 0.1 cm-1 (checked by repeating the same calculations with displacements $`u/a`$=0.0015 and 0.003).

Total energies, pressures, and phonon frequencies were calculated at 12 different volumes ranging from 6.35 Å3/atom down to 3.35 Å3/atom, i.e., for about 10% volume expansion to 40% compression relative to the experimental equilibrium volume of 5.6725 Å3/atom at 300 K Sato02 . The lower volume limit corresponds to roughly 600 GPa maximum pressure.

## III Pressure-volume relationship

### III.1 Analytical form

Figure 1 shows calculated pressures as a function of volume. An excellent consistency of the data points obtained by the three different methods is immediately apparent. Also, the present results substantiate the claims made in Ref. \[Kress99, \] concerning the equivalence of the PAW approach with all-electron methods such as that of Ref. \[soft:WIEN2k, \].

The proper analysis of the $`P(V)`$ results yields the equilibrium properties $`V_0`$, $`B_0=V_0(dP/dV)_0`$, and $`B_0^{}=(dB/dP)_0=(d\mathrm{ln}B/d\mathrm{ln}V)_0`$ (throughout this manuscript the subscript ’zero’ refers to zero (ambient) pressure or equilibrium volume, either calculated or experimental depending on context). One also needs to identify an analytical form of the $`P(V)`$ behavior which, using the above parameters, fits the calculated pressures over the full volume range.

We have tried several of the common equation-of-state formsMurna44 ; Birch78 ; Vinet86 ; Holza96 . The best analytical expression, in the least-squares sense, was identified as

$$P(V)=3B_0X^n(1X)\mathrm{exp}[\eta (1X)]$$
(2)
where $`X=(V/V_0)^{1/3}`$, $`\eta =3B_0^{}/2+1/2n`$, and $`n=7/2`$. With $`n=7/2`$ Eq. 2 is a blend of the VinetVinet86  (or RydbergStace81 ; Gaura02 ) form and the HolzapfelHolza96  expression, for which $`n=2`$ and $`n=5`$, respectively. Actually, the parameter $`n`$ varies between 3.4 and 3.6 for individual sets of calculated $`P(V)`$ data of diamond. Without loss of significant digits in the fitted parameters of interest, we fix $`n`$ at a value of 7/2; in this case the energy versus volume relation, obtained by integration of Eq. 2, can be written without invoking special functions other than the error function:

$`E(V)`$ $`=`$ $`E_0+9B_0V_0\left[f(V)f(V_0)\right]\mathrm{exp}(\eta )/\sqrt{\eta }`$ (3)
$`f(V)`$ $`=`$ $`\sqrt{\pi }(2\eta +1)\mathrm{erf}(\sqrt{\eta }X^2)+2\sqrt{\eta }\mathrm{exp}(\eta X)/X^2.`$
Here, $`E_0`$ is the energy at $`V=V_0`$. Equation 2 was used to fit the $`P(V)`$ results obtained via the stress theorem (PW and PAW methods), while Eq. 3 served to determine the EOS parameters from total energies (APW+lo method). For the PW and PAW methods it was carefully checked that the directly calculated $`P(V)`$ data are consistent with those obtained by differentiation of $`E(V)`$. For all sets of calculated $`P(V)`$ or $`E(V)`$ data the rms deviations were less than 0.04 GPa or 0.05 meV, respectively.

Aleksandrov et al.Aleks87  pointed out that a quadratic dependence of pressure on change in relative density $`\mathrm{\Delta }\rho /\rho _0=V_0/V1`$, written as

$$P(\rho )=B_0\frac{\mathrm{\Delta }\rho }{\rho _0}\left[1+\frac{B_0^{}1}{2}\frac{\mathrm{\Delta }\rho }{\rho _0}\right],$$
(4)
applies in the case of diamond for $`P100`$ GPa. This expression is found to fit the calculated results up to 600 GPa quite well. It is only in the statistical sense that Eq. 4 is slightly inferior to Eq. 2, i.e., the rms deviations are larger by a factor of two, but this is not relevant for the present discussion. The bulk modulus, obtained by differentiation of Eq. 4 with respect to normalized volume, is given by

$$B(\rho )=B_0\frac{\rho }{\rho _0}\left[1+(B_0^{}1)\frac{\mathrm{\Delta }\rho }{\rho _0}\right].$$
(5)
Équation 5 corresponds to a quadratic dependence of the bulk modulus on relative change in density.

### III.2 Calculated EOS parameters

The obtained pairs of ($`V_0`$,$`B_0`$) values are displayed in Fig. 2. Within a given approximation for the exchange-correlation functional (GGA or LDA), the three methods employed here (PW, PAW, APW+lo) yield results for $`V_0`$ and $`B_0`$ which are in very good agreement with each other. The choice of the exchange-correlation functional, however, leads to differences of 4% in $`V_0`$ and 5% in $`B_0`$, reflecting in part the well-known overbinding of the LDA. All $`B_0^{}`$ values (cf. inset to Fig. 2) fall into a narrow range between 3.6 and 3.7 , the difference between averaged LDA and GGA results being close to 1%.

Since the calculated EOS parameters clearly split into two groups, i.e. the GGA and LDA results, with nearly identical parameter values within each group, we have combined the $`P(V)`$ points for GGA and LDA, respectively, to obtain the GGA and LDA parameters $`V_0`$, $`B_0`$, and $`B_0^{}`$ listed in Table 1. The corresponding $`P(V)`$ relations are shown by solid and broken lines in Fig. 1, where the inset illustrates deviations of individual calculated results from the average curves. Similarly, the solid and broken lines in Fig. 2 and its inset show the corresponding volume dependences of the bulk modulus and its pressure derivative.

The larger $`B_0`$ value for LDA compared to GGA correlates with the smaller equilibrium volume. However, the $`B(V)`$ curve for LDA falls below that for GGA (see Fig. 2). This also applies to $`B_0`$ values obtained in other self-consistent LDA calculations for diamond Niels86 ; Fahy87 ; Mailh91 ; Pavon93  (represented by open circles in Fig. 2).

### III.3 Zero-point motion effects

The calculated results for $`V_0`$ and $`B_0`$ should be compared to experimental properties after ’correction’ for vibrational effects. Thermal effects are almost negligible at 300 K because of the high Debye temperature of diamond ($``$2000 K), but zero-point motion needs to be considered. According to Monte Carlo simulationsHerre00  the static-lattice equilibrium volume for diamond is about 5.59 Å3. All the data for $`V_0`$ shown in Fig. 2 fall into a range of $`\pm 2`$% around this value. The static-lattice bulk modulus at $`V=5.59`$ Å3 is estimated to be 469 GPa (this follows from the experimental values of $`V_0`$ and $`B_0`$ and the calculated $`B_0^{}`$). Thus, relative to static-lattice properties, the GGA results cannot be considered superior to the LDA ones; they appear better when compared to experimental data because the GGA errors happen to mimic the zero-point effects.

The zero-point vibrational pressure $`P^{ZP}`$, i.e., the isochoric change in pressure when zero-point motion is switched on, is about 6 GPa at the static-lattice value of $`V_0`$Herre00 . An approximate relation for the volume dependence of $`P^{ZP}`$ is

$$P^{ZP}(V)\frac{dP(V)}{dV}\mathrm{\Delta }V^{ZP}(V),$$

where $`\mathrm{\Delta }V^{ZP}`$ is the isobaric volume expansion due to zero-point motion. For getting $`\mathrm{\Delta }V^{ZP}(V)`$ within the quasi-harmonic thermodynamics, we simply refer to the calculated pressure dependence of the zero-point expansion presented in Fig. 10 of Ref. \[Herre00, \]. The corresponding volume dependence of $`P^{ZP}`$ is displayed in the inset of Fig. 1. $`P^{ZP}`$ increases by about 6 GPa at 600 GPa. Adding such small changes in $`P^{ZP}`$ to the calculated $`P(V)`$ results leads to an increase of the $`B_0`$ values by about 1%, but would not affect the calculated $`B_0^{}`$ values (within the estimated uncertainty of the optimum value given below).

### III.4 The essence of the EOS calculations

In Fig. 3, the calculated $`P(V)`$ relations are plotted in reduced coordinates, i.e. pressure normalized by $`B_0(\mathrm{\Delta }\rho /\rho _0)`$ as a function of $`\mathrm{\Delta }\rho /\rho _0`$. In this representation small differences in ambient-pressure volumes and bulk moduli are suppressed and a linear slope corresponds to a quadratic dependence of pressure on change in relative density (Eq. 4). The LDA and GGA results hardly differ in slope. In this way, Fig. 3 illustrates the main information we extract from our EOS calculations for diamond:

(1) The value of $`B_0^{}`$ does not depend much on crystal potential and basis set issues and on approximations for the exchange-correlation functional. Our average value is

$$B_0^{}=3.65\pm 0.05,$$

where the small uncertainty reflects the scatter for the different methods of calculation.

(2) Within the volume range considered here, the most appropriate 3-parameter analytical form for the $`P(V)`$ relation of diamond (Eq. 2 or 4) is transferable between LDA and GGA solutions; only $`V_0`$ and $`B_0`$ need to be adjusted.

(3) Inserting the calculated $`B_0^{}`$ value and the experimental data for $`V_0`$ and $`B_0`$ (see Table 1) into Eq. 2 (Eq. 4) is considered to yield, on the basis of this work, the optimum representation of the EOS of diamond at 300 K. Actually, at the experimental $`V_0`$ the GGA bulk modulus (calculated value of 440 GPa plus 1% correction for zero-point motion effects) happens to be very close to the experimental $`B_0`$, and the calculated pressure at the experimental $`V_0`$ is only 2 GPa. Therefore, the optimum EOS corresponds to the $`P(V)`$ relation obtained from GGA calculations, with only a small correction applied.

### III.5 Comparison with other results for $`B_0^{}`$

Calculated values of $`B_0^{}`$ reported in the literature are 3.54 \[Cheli84, \], 3.6, 4.5 \[Niels86, \], 3.24 \[Fahy87, \], 3.5 \[Pavon93, \] and 3.97 \[Herre00, \]. In some cases, the differences to our result are small. Larger deviations may in part result from the procedure used to extract the EOS parameters (e.g., fit of a Murnaghan equation which is inadequate).

The experimental value $`B_0^{}`$ =4 $`\pm `$ 0.5 (Ref. \[McSki72, \]) derived from sound speed measurements at low pressures (0.2 GPa) is not sufficiently accurate to test the calculated result.

The only other experimental value of $`B_0^{}`$ stems from the recent x-ray diffraction experiments up to 140 GPaOccelli02 . The EOS data appear to be of high quality. They were measured using helium as a pressure medium which is considered to provide almost hydrostatic conditions in DAC experiments. The obtained value of $`B_0^{}=3.0(1)`$ is inconsistent with the ultrasonic result and also significantly smaller than the calculated value. The latter difference is illustrated in Fig. 3. At nominally 140 GPa, about the maximum pressure reached in the experiments, the difference between our calculated $`B_0^{}`$ value and that Ref. \[Occelli02, \] translates to a pressure difference of about +11 GPa. This means that either repulsion is slightly overestimated in the calculations (independent of the exchange-correlation functional) or that the experimental data suffer from systematic errors, a combination of these effects not being ruled out. We will return to this issue below.

## IV Optical phonon frequency

### IV.1 Calculated phonon frequency versus volume

Figure 4 shows calculated F2g phonon frequencies as a function of volume. At a given volume, the frequencies of all calculations are quite consistent with each other. The frequencies are too low by about 3% if compared to the experimental phonon frequency at ambient pressure. This small discrepancy does not originate from anharmonicity connected to the frozen-phonon displacement.

The volume dependence of phonon frequencies is usually characterized by the simple scaling law

$$\frac{\omega (V)}{\omega _0}=\left(\frac{V}{V_0}\right)^\gamma =\left(\frac{\rho }{\rho _0}\right)^\gamma ,$$
(6)
which assumes that the mode Grüneisen parameter $`\gamma =d\mathrm{ln}\omega /d\mathrm{ln}V`$ is independent of volume. The frozen-phonon results indicate some volume dependence of $`\gamma `$, an observation which was also noted by NielsenNiels86 . The expression

$$\frac{\mathrm{\Delta }\omega (V)}{\omega _0}=\frac{\gamma _0}{\gamma _0^{}}\left[\left(\frac{V}{V_0}\right)^{\gamma _0^{}}1\right].$$
(7)
with $`\mathrm{\Delta }\omega (V)=\omega (V)\omega _0`$ yields a slightly better match of our frozen-phonon results (3 cm-1 rms deviation or better for individual sets). Within standard deviations, the parameters $`\gamma _0`$ and $`\gamma _0^{}`$ for the combined GGA and LDA results are identical (see Table 1). The upper inset to Fig. 4 illustrates the small volume dependence of the mode Grüneisen parameter resulting from $`\gamma _0^{}<1`$.

With $`\rho /\rho _0\omega /\omega _0`$ ($`\gamma 1`$) it follows from Eq. 5 that the ratio of normalized bulk modulus to normalized phonon frequency is approximately linear in relative density and the slope is $`B_0^{}1`$:

$$\frac{B}{B_0}\left(\frac{\omega }{\omega _0}\right)^11+(B_0^{}1)\frac{\mathrm{\Delta }\rho }{\rho _0}.$$
(8)
This scaling between bulk modulus and optical phonon frequency could come handy as a simple relation between an elastic and a dynamical property of diamond.

### IV.2 Calculated vs. experimental mode Grüneisen parameter

To compare the calculated phonon results with experimental data, we first consider the measured linear pressure coefficient of the phonon frequency. According to Schiferl et al.Schif97 , the ’best’ experimental value in the literature, taken as the average of Refs. \[Hanfl85b, ; Boppa85, ; Muino94, \], should be $`(d\omega /dP)_0=2.90\pm 0.05`$ cm-1/GPa ($`\gamma _00.96`$). This value originates from experiments which employed the ruby pressure calibration with a linear coefficient $`A=(dP/d\mathrm{ln}\lambda )_0=1905(10)`$ GPa for the R1 line wavelength ($`\lambda `$) shift near zero pressure as determined by Piermarini et al.Pierm75  Recent high-precision measurements of the ruby line shift up to 1 GPaGrass01  and a reinterpretationHolza02  of R1 line shift data measured up to $``$20 GPaPierm75 ; Nakan00  indicate that the parameter $`A`$ is smaller ($`A=1820\pm 30`$ GPa) compared to the previously accepted value. This revision leads to a corresponding increase of the phonon pressure coefficient. In this context it is helpful that the ratio of the diamond phonon frequency shift to ruby wavelength shift is explicitly given in Refs. \[Hanfl85b, \] and \[Aleks87, \]. The results, 7.94(10) cm-1/nm and 7.85(6) cm-1/nm, agree quite well with each other. Taking the average value (7.9(1) cm-1/nm) in combination with $`A`$ = 1820(30) GPa, the corrected value from DAC experiments should be $`(d\omega /dP)_0=3.0(1)`$ cm-1/GPa. Within experimental uncertainties, this linear pressure coefficient agrees with results obtained by methods which do not involve any ruby calibrationMitra69 ; Whalley76 ; Grims78 , in particular that of Whalley et al.Whalley76  With $`B_0`$=444.8 GPa the corresponding mode Grüneisen parameter becomes $`\gamma _0=1.00(3)`$, the error being mainly due to the uncertainty of the ruby coefficient $`A`$. We note the excellent agreement of the corrected experimental value for $`\gamma _0`$ with the results obtained within the LDA and GGA approximations (see Table 1).

Occelli et al. report an average mode Grüneisen parameter of $`\overline{\gamma }=0.97`$ for the range 0–140 GPa. It should be pointed out that their value of $`\overline{\gamma }`$ is independent of the pressure scale used in the experiments because the volume was measured directly. Taking into account the predicted volume dependence $`\gamma (V)`$ (see inset to Fig. 4), the calculations are fully consistent with the experimental volume dependence of the phonon frequency. This is no longer the case if we turn to the pressure dependence.

### IV.3 The pressure – phonon-frequency relationship

One can solve Eq. 7 for the relative volume, insert the corresponding expression into Eq. 2 (Eq. 4), and use the ’GGA’ and ’LDA’ parameters from Table 1 to obtain an analytical form for pressure as a function of phonon frequency. On the other hand, in the analysis of earlier phonon calculationsHanfl85b ; Niels86 , the $`P(\omega )`$ behavior was parametrized using the analog of a Birch expressionBirch78

$$P(X)=\frac{3}{2}b_0(X^7X^5)[1+\eta (1X^2)]$$
(9)
with $`X=(\omega /\omega _0)^{1/3}`$. Equation 9 (more or less an ad hoc choice in the earlier work) happens to yield excellent representations of the present calculated $`P(\omega )`$ results. To facilitate direct comparison with Refs. \[Hanfl85b, \] and \[Niels86, \], the parameters $`b_0`$ and $`\eta `$ are summarized in Table 2. Note that for our calculated data $`\eta `$ is the only adjustable parameter, because $`b_0=B_0/\gamma _0`$.

Applicability of Eq. 4 in combination with $`\gamma 1`$ implies a nearly quadratic dependence of pressure on change in phonon frequency. Thus, we compare calculated and experimental results in terms of reduced coordinates $`P/\mathrm{\Delta }\omega `$ and $`\mathrm{\Delta }\omega `$, cf. Fig. 5. Most of the calculated results represented in Fig. 5, i.e., the present ones within LDA and GGA and those of Refs. \[Hanfl85b, \] and \[Niels86, \], agree very well with respect to the slope. This is of course reflected in the small scatter of the $`\eta `$ values given in Table 2. The average value from our calculations is

$$\eta =0.068\pm 0.03.$$

Inserting this value and the experimental data for $`\omega _0`$ and $`b_0=\omega _0(dP/d\omega )_0`$ (Table 2) into Eq. 9 yields the $`P(\omega )`$ dependence consistent with experimental data near ambient pressure and with the nonlinear behavior predicted by the present calculations.

In the reduced coordinates of Fig. 5 the calculated $`P(\omega )`$ results of Wu and XuWu99  do not agree well with those of other calculations. They use the quadratic function

$$\omega (P)=\omega _0+a_1P+a_2P^2$$
(10)
to fit their data. The choice of pressure as the independent variable in a quadratic expression is not appropriate in view of Eq. 4 and $`\gamma 1`$. It leads to parameter correlation which is the possible reason that their results exhibit curvature in Fig. 5.

Occelli et al.Occelli02  have measured the optical phonon frequency for pressures up to 140 GPa. They also give frequency as a quadratic function of pressure, i.e, $`a_1=2.83`$ cm-1/GPa and $`a_2=3.65\times 10^3`$ cm-1/GPa2 in Eq. 10. The curvature of their data in Fig. 5 could again be related to the particular choice of the fitting function. The main observation, however, is that the average slope for the experimental data differs from the one predicted by theory. At the frequency corresponding to 140 GPa experimental pressure the pressure difference between calculations and experiment again amounts to about 11 GPa.

The sign and magnitude of this difference is very similar to that encountered when comparing experimental and calculated EOS results. Thus, in terms of pressure we find differences between calculated and experimental results which are about the same for different physical quantities considered. This hints to an explanation where the discrepancies are caused by the same systematic error(s).

## V Remarks on the ruby calibration

We consider the possibility that the discrepancy in the high pressure regime between experimental and calculated EOS and phonon frequency results for diamond are caused by some error in the experimental pressure scale, i.e. the calibration of the ruby R1 line shift according to Ref. \[Mao86, \].

HolzapfelHolza02  recently proposed a revised ruby pressure scale, based on an analysis of published EOS data of selected elemental solids, including those of Ref. \[Occelli02, \] for diamond. The revised ruby calibration was cast into a three-parameter analytical expression for pressure as a function of the R1 line wavelength. With the recommended parameter valuesHolza02  and when restricted to about 200 GPa pressure (a range beyond the upper limit of most DAC studies), the revised calibration can be represented by a simple second-order polynomial in frequency shift $`\mathrm{\Delta }\nu `$ of the R1 line. This is evident from Fig. 6 which shows plots of $`(|\mathrm{\Delta }\nu |/\nu _0)^1\times P`$ versus $`|\mathrm{\Delta }\nu |/\nu _0`$ for Holzapfel’s revised scale and for the 1978/1986 scales of Mao et al.Mao78 ; Mao86  We write the quadratic polynomial aslambda

$$P(\mathrm{\Delta }\nu )=A\frac{|\mathrm{\Delta }\nu |}{\nu _0}\left[1+\frac{B}{2}\frac{|\mathrm{\Delta }\nu |}{\nu _0}\right].$$
(11)
A linear regression for the ’Revised’ line in Fig. 6 with $`A=1820`$ GPa (fixed value as recommended in Ref. \[Holza02, \]) gives a slope parameter of $`B`$=15.8.

From the experimental $`P(V)`$ and $`P(\omega )`$ data of diamondOccelli02  one can recover the ruby line shift $`\mathrm{\Delta }\nu `$ according to the ruby calibrationMao86  used in the experimental work. Combining the experimental $`V(\mathrm{\Delta }\nu )`$ and $`\omega (\mathrm{\Delta }\nu )`$, respectively, with the calculated $`P(V)`$ and $`P(\omega )`$ relations reported here yields $`P(\mathrm{\Delta }\nu )`$ as shown by the two solid lines in Fig. 6. These curves run more or less parallel to Holzapfel’s revised calibration. Obviously, both the calculated $`P(V)`$ and $`P(\omega )`$ relations for diamond in combination with the experimental data of Ref. \[Occelli02, \] support the proposed revision of the ruby pressure scale. It should be noted that the ruby calibration discussed by Aleksandrov et al.Aleks87  would be a little higher in pressure, but in the coordinates of Fig. 6 it exhibits a slope similar to that of the ’revised’ line.

## VI Conclusions

We offer the following conclusions:

(1) The results of the first-principles EOS and phonon frequency calculations for diamond reported here do not depend on the computational method, i.e., the choice of the crystal potential and basis sets in the PW, PAW, and APW+lo methods.

(2) The calculated equilibrium values for volume ($`V_0`$) and bulk modulus ($`B_0`$) of diamond do depend on the exchange-correlation functional (LDA or GGA), a well-known fact in density functional theory. However, the value of the pressure derivative of the bulk modulus \[$`B_0^{}`$=3.65(5)\] and the analytical form of the $`P(V)`$ relation are not affected within the volume range (40% compression) covered here. In other words, the calculated nonlinear component in the $`P(V)`$ behavior is independent of the method of calculation and the exchange-correlation approximation.

(3) The different theoretical methods and approximations also yield very similar results for the nonlinearity in the pressure versus phonon frequency relation $`P(\omega )`$.

(4) Both the calculated $`P(V)`$ and $`P(\omega )`$ relations exhibit a nonlinear behavior which is more pronounced than that observed in recent x-ray diffraction and Raman measurements of diamond up to 140 GPaOccelli02 . In terms of absolute pressure, the differences between theoretical and experimental results are essentially the same for the $`P(V)`$ and $`P(\omega )`$ behavior, of the order of 10 GPa at an experimental pressure of 140 GPa.

(5) Our calculated nonlinearity of $`P(\omega )`$ agrees well with that obtained in earlier ab initio phonon frequency calculations for diamond under pressureHanfl85b ; Niels86 . The work of Hanfland et al.Hanfl85b  indicated the need to reconsider the calibration of the ruby pressure scale used at that timeMao78 . The ’quasi-hydrostatic’ (1986) ruby calibrationMao86  reduced the discrepancy between theoretical predictions and experimental results. The remaining difference discussed here at least calls for some caution when experimental data measured above 50 GPa and based on the 1986 ruby pressure scale are compared to ab initio calculations or reduced shock wave data.

###### Acknowledgements.

The authors thank P. Loubeyre and W. B. Holzapfel for making available preprints of Refs. \[Occelli02, \] and \[Holza02, \], respectively. Part of the computer resources used in this work were provided by the Scientific Committee of IDRIS (Institut du Développement et des Ressources en Informatique Scientifique), Orsay (France).