text:
# References

Double Symmetries in Field Theories

L.M.Slad 111E-mail: slad@theory.npi.msu.su

D.V.Skobeltsyn Institute of Nuclear Physics,
Moscow State University, Moscow 119899

Abstract

In the paper a concept of a double symmetry is introduced, and its qualitative characteristics and rigorous definitions are given. We describe two ways to construct the double-symmetric field theories and present an example demonstrating the high efficiency of one of them. In noting the existing double-symmetric theories we draw attention at a dual status of the group $`SU(2)_LSU(2)_R`$ as a secondary symmetry group, and in this connexion we briefly discuss logically possible aspects of the $`P`$-violation in weak interactions.

1. Qualitative characteristics of the double symmetry

Extension of symmetry approaches to the field theories construction is of an ordinary practice in physics. In the present paper we propose some generalization of already existing approaches which consist in constructing a new group $`𝒢_T`$, called a double symmetry group, on the basis of the given global group $`G`$ (the primary symmetry group) and some its representation $`T`$. This construction is carried out in the framework of group transformations in field-vectors space. Namely, the transformations of the double symmetry group consist of transformations of the primary and secondary symmetries, realized in the same field-vectors space of some representation $`S`$ of the group $`G`$. The secondary symmetry transformations, which can be global as well as local, have two features. Firstly, the parameters of these transformations belong to the space of the representation $`T`$ of the group $`G`$. Secondly, the secondary symmetry transformations do not violate the primary symmetry.

As a rule, the secondary symmetry transformations couple different irreducible representations, into which the representations $`S`$ of the group $`G`$ is decomposed, removing or reducing the arbitrariness of a field theory allowed by its invariance in respect to $`G`$. If one or another secondary symmetry corresponds to the real physical world, then this correspondence is almost inevitably concealed by its spontaneous breaking.

We restrict ourselves with formulations of the global symmetries. The knowledge of the secondary symmetry group allows to introduce the relevant local transformations and gauge fields in the standard way, if there is any need of this.

2. Rigorous definitions of the secondary and double symmetries

Definition 1. Assume that there are a symmetry group $`G`$ of some field theory and two its representations $`T`$ and $`S`$. Let $`\theta =\{\theta _a\}`$ be a vector in the representation space of $`T`$, $`\mathrm{\Psi }(x)`$ be any field vector in the representation space of $`S`$, and let $`D^a`$ be such operators that the field $`\mathrm{\Psi }^{}(x)`$, obtaned by the transformation

$$\mathrm{\Psi }^{}(x)=\mathrm{exp}(iD^a\theta _a)\mathrm{\Psi }(x),$$
(1)
belongs again to the representation space of $`S`$, i.e. for any $`gG`$

$$\mathrm{exp}(iD^b(T(g)\theta )_b)S(g)\mathrm{\Psi }(x)=S(g)\mathrm{\Psi }^{}(x).$$
(2)
Then the transformations (1) and their products will be called secondary symmetry transformations produced by the representation $`T`$ of the group $`G`$.

The relations (1) and (2) allows to give a more general definition.

Definition 2. Let the group $`𝒢_T`$ contain the subgroup $`G`$ and the invariant subgroup $`H_T`$, with $`G𝒢_T`$, $`H_T𝒢_T`$, and $`𝒢_T=H_TG`$. If any element $`h_TH_T`$ can be written in the form

$$h_T=h(\theta _1)h(\theta _2)\mathrm{}h(\theta _n),n\{1,2,\mathrm{}\},$$
(3)
where $`\theta _i`$ $`(i=1,2,\mathrm{},n)`$ is some vector of the representation space of $`T`$ of the group $`G`$, and if for any element $`gG`$

$$gh(\theta )g^1=h(T(g)\theta ),$$
(4)
then the group $`G`$ will be called the primary symmetry group, and the groups $`H_T`$ and $`𝒢_T`$ will be respectively called the secondary and double symmetry groups produced by the representation $`T`$ of the group $`G`$.

It follows from the relations (1) and (2) that the operators $`D^a`$ must satisfy the condition

$$D^a=S^1(g)D^bS(g)[T(g)]_b{}_{}{}^{a},$$
(5)
i.e., as one uses to say, such operators $`D^a`$ transform as the representation $`T`$ of the group $`G`$.

If $`G`$ is a Lie group, then within some vicinity $`U`$ of the unity element $`e`$ of the group $`G`$ the transformation operators $`S(g)`$ and $`T(g)`$ can be written as

$$S(g)=\mathrm{exp}(iL^jϵ_j),$$
(6)
$$T(g)=\mathrm{exp}(iM^jϵ_j),$$
(7)
where $`L^j`$ and $`M^j`$ are generators of the group $`G`$ of the representations $`S`$ and $`T`$, $`ϵ_j=ϵ_j(g)`$ are group parameters corresponding to the given element $`gU`$.

By using Eqs. (6) and (7), condition (5) takes the form

$$[L^j,D^a]=D^b(M^j)_b{}_{}{}^{a}.$$
(8)
Expanding the exponent of Eq. (1) in a power series, we see that the first term $`(\mathrm{\Psi })`$ transforms as $`S`$ and the second term $`(\theta \mathrm{\Psi })`$ does as a direct product $`TS`$. The following statement is getting obvious. Nontrivial ($`D^a0`$) secondary symmetry transformations (1) exist if and only if among irreducible representations of the direct product $`TS`$ there is at least one which belongs to the representation $`S`$. If the operators $`D^a`$ are matrix ones, then their elements are nothing but the unnormalized Clebsch-Gordon coefficients for relevant irreducible representations of the group $`G`$ to be found from the relation (8). Namely, let $`\omega _t`$, $`\omega _p`$, and $`\omega _q`$ denote three irreducible representations of the group $`G`$, such that $`\omega _tT`$, $`\omega _pS`$, and $`\omega _qTS`$, and let $`\alpha _t`$, $`\alpha _p`$ and $`\alpha _q`$ be a set of indices that characterize vectors in the spaces of relevant irreducible representations. Then we have the well-known relations

$$D_{\omega _q\alpha _q,\omega _p\alpha _p}^{\omega _t\alpha _t}=0,\mathrm{if}\omega _qS,$$
(9)
$$D_{\omega _q\alpha _q,\omega _p\alpha _p}^{\omega _t\alpha _t}=d_{\omega _q\omega _p}^{\omega _t}(\omega _t\alpha _t\omega _p\alpha _p|\omega _q\alpha _q),\mathrm{if}\omega _qS,$$
(10)
where $`(\omega _t\alpha _t\omega _p\alpha _p|\omega _q\alpha _q)`$ are the properly normalized Clebsch-Gordon coefficients, $`d_{\omega _q\omega _p}^{\omega _t}`$ are arbitrary quantities independent on the indices $`\alpha _t`$, $`\alpha _p`$ and $`\alpha _q`$.

3. Two ways to construct double-symmstric field theories

Structure of the groups $`H_T`$ and $`𝒢_𝒯`$, produced by the representation $`T`$ of the group $`G`$, depends on the representation $`S`$ of the group $`G`$ and on numerical values of the quantities determining the operators $`D^a`$.

If the secondary symmetry is produced by the adjoint representation of the group $`G`$ and the operators $`D^a`$ of Eq. (1) coinside with the group generators, then obviously the group $`H_T`$ is locally isomorphic to the group $`G`$. In such a case we shall say that the double symmetry is degenerated one.

In the general case there are at least two ways to complete the construction of the group $`𝒢_𝒯`$ and to construct the double-symmetric field theories, i.e. the theories whose Lagrangians are invariant under transformations of the group $`𝒢_𝒯`$.

The first way is to close the algebra of the operators $`D^a`$ starting from those or other considerations on the fields in question and their interactions. Possessing the Lie algebra of the group $`𝒢_𝒯`$ one can find then its representations and allowed double-symmetric theories. One proceeds so in the supersymmetry theory \[1-3\].

Desribing the second way for some fixed representation $`T`$, we consider, for simplicity, the matrix realization of the operators $`D^a`$. First we choose some class of representations of the group $`G`$ which seem to be admissible for the field $`\mathrm{\Psi }(x)`$. Some set of arbitrary quantities $`d_{\omega _q\omega _p}^{\omega _t}`$ determining the operators $`D^a`$ through Eq. (10) corresponds to each representation of this class. It is required the Lagrangian of the considered theory to be invariant under transformations of the group $`G`$ as well as under global transformations (1). Then, obviously, it will be invariant under any transformations of the groups $`H_T`$ and $`𝒢_𝒯`$. These requirements can lead to some selection of the representations of the considered class, and restrict the arbitrariness of quantities $`d_{\omega _q\omega _p}^{\omega _t}`$ and the arbitrariness of Lagrangian constants allowed by the invariance in respect to the group $`G`$. If the arbitrariness of quantities $`d_{\omega _q\omega _p}^{\omega _t}`$ is removed completely (up to a common constant) for each of selected representations of the group $`G`$, then the structure of the group $`H_T`$ becomes automatically fully defined, but, generally speaking, it is unlike for different representations of the group $`G`$. If it is not so, only then a question arises on closing the algebra of the operators $`D^a`$ taking into account the obtained restrictions. Thus, on this way the first place is taken by constraints established by the double symmetry produced by the representation $`T`$ of the group $`G`$, and the problem on the Lie algebra of the group $`𝒢_𝒯`$ is resolved automatically or moved onto the last place. It seems to be reasonable to demonstrate the efficiency of the described way with the help of a nontrivial example. Let us do this.

4. An example of constructing double-symmetric field theory: the efficiency of selecting representations and removing ambiguities

Consider the relativistically invariant 222Here, as well as in the monograph , the relativistic invariance means, in modern titles and notations , an invariance in respect to the orthochronous Lorentz group $`L^{}`$ generated by the proper Lorentz group $`L_+^{}`$ and spatial reflection $`P`$. Lagrangian of the general type  for some free fermionic field $`\mathrm{\Psi }(x)`$

$$_0=\frac{i}{2}[(_\mu \mathrm{\Psi },L^\mu \mathrm{\Psi })(\mathrm{\Psi },L^\mu _\mu \mathrm{\Psi })]\kappa (\mathrm{\Psi },\mathrm{\Psi }),$$
(11)
where $`(\mathrm{\Psi }_1,\mathrm{\Psi }_2)`$ is a relativistically invariant belinear form, $`L^\mu `$ are matrix operators, and $`\kappa `$ is a constant.

Require the Lagrangian (11) to be invariant under the global secondary symmetry transformations produced by a polar 4-vector of the group $`L^{}`$

$$\mathrm{\Psi }^{}(x)=\mathrm{exp}(iD^\mu \theta ^\mu )\mathrm{\Psi }(x),$$
(12)
where $`D^\mu `$ are matrix operators. This requirement will be fulfiled if

$$[L^\mu ,D^\nu ]=0,$$
(13)
$$(D^0\mathrm{\Psi }_1,\mathrm{\Psi }_2)=(\mathrm{\Psi }_1,D^0\mathrm{\Psi }_2).$$
(14)
The condition (5) for the operators $`D^\mu `$ and the condition, which the operators $`L^\mu `$ of the Lagrangian (11) obey, are equivalent. Therefore the matrix elements of both operators $`L^\mu `$ and $`D^\mu `$ are given by relations of the type (9-10). In the monograph  there are given their explicit forms. We shall use notations 333In these notations all irreducible representations of the group $`L_+^{}`$ are elegantly described: finite- and infinite-dimensional ones. Transition to the often used notations $`(j_1,j_2)`$, which are connected with the group $`SO(4)=SO(3)SO(3)`$ and describe only finite-dimensional irreducible representations, is the following: $`j_1=(l_1+l_01)/2`$, $`j_2=(l_1l_01)/2`$. $`\tau =(l_0,l_1)`$ of Ref.  for irreducible representations of the proper Lorentz group and denote the arbitrary quantities of Eq. (10), characterizing the operators $`L^\mu `$ and $`D^\mu `$, as $`c_{\tau ^{}\tau }`$ and $`d_{\tau ^{}\tau }`$, respectively.

Let the representation $`S`$ of the group $`L`$, as which the field $`\mathrm{\Psi }(x)`$ of the Lagrangian (11) transforms, be decomposible into a finite or infinite direct sum of irreducible finite-dimentional representations with half-integer spins, and let the multiplicity of any of these irreducible representations do not exeed 1. Then the Lagrangian (11) will be invariant under global transformations (12) if and only if $`S`$ is one of the infinite-dimentional representations $`S^{k_1}`$

$$S^{k_1}=\underset{n_1=0}{\overset{\mathrm{}}{}}\underset{k_0=k_1+1}{\overset{k_11}{}}(k_0,k_1+n_1),k_1=3/2,5/2,\mathrm{}$$
(15)
or is the representation $`S^F`$ containing all irreducible representations of the group $`L_+^{}`$ with half-integer spins.

For each of these representations all the quantities $`|c_{\tau ^{}\tau }|`$ and $`|d_{\tau ^{}\tau }|`$ have, up to their common constants, fully defined values. For any of the representations (15) we get in result

$$[D^\mu ,D^\nu ]=0,$$
(16)
i.e. the secondary symmetry group $`H_T`$ is the four-parameter Abelian group. Notice that, although in this case the Lie algebra of the double symmetry group $`𝒢_T`$ coinsides with the Lie algebra of the Poincarè group, the operators $`D^\mu `$ cannot be identified with the translation operators $`P^\mu `$ because $`D^\mu `$ act only on spin variables of a field but $`P^\mu `$ do not act on them.

The mass spectra corresponding to each of the obtained Lagrangians (11) are infinetely degenerated in spin and continious. In order to eliminate this degeneration it is necessary to break the secondary symmetry spontanously keeping the orthochronous Lorentz group symmetry.

5. On a perspective to use the double symmetry

The considered example of the double symmetry, playing here only the methodological role, is used by us as one of elements for constructing a theory with the fields which transform as orthochronous Lorentz group representations decomposible into an infinite direct sum of finite-dimentional irreducible representations 444Two our articles on the theory of infinite-component fields with the double symmetry, produced by the polar and axial 4-vectors of the group $`L^{}`$, are preparing for publication. Former the attemps to analize the theory with such class fields were not undertaken because of the infinite number of arbitrary constants and the absence of criterions for removing this arbitrariness. Carried out sometime ago investigations of the infinite-component field theories, which were assigned to an alternative describtion of hadrons as composite particles, were based on those representations of the Lorentz group that could be decomposed into a finite direct sum of infinite-dimentional irreducible representations or on a special representation of the group $`SO(4,1)`$ or $`SO(4,2)`$. It was established, however, that such theories possesses some properties (peculiarities of mass spectra, nonlocality, violation of $`CPT`$-invariance and the connection between spin and statistics) which are not admissible for particle physics. Problems, covered by these investigations as well as a sufficiently full list of relevant papers, can be found in the monograph .

6. Existing undegenerated double symmetry: supersymmetry

Notice now that sypersymmetry both in the $`x`$-space and in the superspace can be considered as a double symmetry produced by the bispinor representation $`T`$ of the proper Lorentz group. An element $`h(\theta )`$ of the corresponding secondary symmetry group $`H_T`$ has a form

$$h(\theta )=\mathrm{exp}(iQ_\alpha \theta ^\alpha +i\overline{Q}_{\dot{\alpha }}\overline{\theta }^{\dot{\alpha }}),$$
(17)
the parameters $`\theta ^\alpha `$ and $`\overline{\theta }^{\dot{\alpha }}`$ being anticommutating elements of the Grassmann algebra and belonging to the representation spaces of $`(1/2,0)`$ and $`(0,1/2)`$ of the proper Lorentz group, respectively (in the notations connected with the group $`SO(4)`$). A fulfilment of the relation (8), as a criterion of the secondary symmetry, was already required in respect to the bispinor operators $`Q_\alpha `$ and $`\overline{Q}_{\dot{\alpha }}`$ in the first paper on the supersymmetry algebra . In the supersymmetry theory the operators $`Q_\alpha `$ and $`\overline{Q}_{\dot{\alpha }}`$ are given in the form of sum of terms containing a first or zero power of the differential operator $`_\mu `$ (cf. \[2-3\]). Such realization of the operators $`Q_\alpha `$ and $`\overline{Q}_{\dot{\alpha }}`$ leads to that the secondary symmetry group $`H_T`$, generated by the elements (17), contains the space-time translation group, and the double symmetry group $`𝒢_𝒯`$ contains the Poincarè group. If one chose a matrix realization for the operators $`Q_\alpha `$ and $`\overline{Q}_{\dot{\alpha }}`$, then we would come out of the supersymmetry theory standards.

Let for a given representation $`S`$ of the group $`L_+^{}`$ the relation (8) for the spinor operators $`Q_\alpha `$ in some its realization keep arbitrary $`N`$ constants. Then there are $`N`$ linearly independent spinor operators $`Q_\alpha ^1`$, $`Q_\alpha ^2`$, $`\mathrm{}`$, $`Q_\alpha ^N`$, in this realization, satisfying the relation (8). In this case, on one hand, one can introduce a secondary symmetry produced by the $`N`$-multiple bispinor representation of the proper Lorentz group, and, on another hand, treat the numbers $`1,2,\mathrm{},N`$ of the spinor operators as an index related to some inner symmetry. If as well the operators $`Q_\alpha ^j`$ and $`\overline{Q}_{\dot{\alpha }}^j`$ ($`j=1,\mathrm{},N`$) keep a linear dependence on the operator $`_\mu `$, then the relevant double symmetry is an extended supersymmetry (cf. \[2-3\]).

7. Existing undegenerated double symmetry: the œ-model symmetry

Another example of an undegenerated double symmetry, used in particle physics, is the $`\sigma `$-model symmetry, the most accurately described by Gell-Mann and Levy . Infinitesimal transformations in the $`\sigma `$-model have the forms

$$N^{}=(1\frac{i}{2}\gamma ^5\text{ø}\text{`})N,$$
(18)
$$\text{ß}^{}=\text{ß}+i\text{`}\sigma ,$$
(19)
$$\sigma ^{}=\sigma i\text{`}\text{ß},$$
(20)
where the field $`N`$ is a nucleon isodoublet, the field ß is a pseudoscalar isotriplet, the field $`\sigma `$ is a scalar isosinglet, and the transformation parameter \` is a pseudoscalar isotriplet.

The listed transformation properties of the field and the parameter \` allow to state that the transformations (18) and (19-20) are transformations of the secondary symmetry produced by the representation $`T`$ = (isotriplet, pseudoscalar) or, that is equivalent, by the representation $`T`$ = (isotriplet, scalar) $``$ (isotriplet, pseudoscalar) of the group $`G=SU(2)L^{}`$. We wish especially emphasize that, firstly, transformations (18) and (19-20) do not violate the spatial reflection symmetry, and so the group $`G`$ contains the orthochronous Lorentz group, and, secondly, the parity of fields, involved in the transformations (19-20), are necessarily different due to a pseudoscalar character of the paremeter \`. $`G_T=SU(2)_LSU(2)_R`$ is a group of the secondary symmetry generated by the transformations (18) and (19-20), the parameters of one of the group $`SU(2)`$ being given by sum of the space scalar and pseudoscalar, and the other group parameters being given by their difference.

8. Dual status of the group $`𝑺𝑼\mathbf{(}\mathrm{𝟐}\mathbf{)}_𝑳\mathbf{}𝑺𝑼\mathbf{(}\mathrm{𝟐}\mathbf{)}_𝑹`$

In numerous works including the pioneer papers by Schwinger , Gürsey  and Touschek , which used the chiral symmetry group $`SU(N)_LSU(N)_R`$ or the group $`U(1)_A`$, it is difficult, if it’s really possible, to find any commentaries on transformation properties of these group parameters in respect to the spatial reflection. As a rule it follows from the contents of the papers that the transformations of the group $`SU(N)_LSU(N)_R`$ or the group $`U(1)_A`$ do not violate $`P`$-symmetry. This rule, however, is broken in the left-right symmetric model of electroweak interactions \[10-12\] and in the unifited models of strong and electroweak interactions including the first one as its element. Actually, in the papers \[10-12\] there are not at all indications that any component of one or another Higgs multiplet is the space pseudoscalar. We are forced to admit that all components of all Higgs multiplets, and, consequently, all parameters of the used group $`SU(2)_LSU(2)_R`$ are space scalars. Transformations of this group, generally speaking, violate the $`P`$-symmetry. If, for example, the initial state of a spinor field possesses a definite parity, then the state, obtained due to a transformation of the type (18), already does not possess it. The group in question itself does not establish any relation between coupling constants of two $`W`$-bosons to fermions. Therefore in the left-right symmetric model one introduces some discrete symmetry which transforms $`SU(2)_L`$ to $`SU(2)_R`$ but it is not identified with the spatial reflection. So one can state that the group $`SU(2)_LSU(2)_R`$ of this model is a secondary symmetry group produced by 2-multiple scalar isotriplet of the group $`G=SU(2)L_+^{}`$ (A discrete symmetry group), this discrete symmetry being not fully defined.

9. On $`P`$-properties of the physical vacuum and the gauge fields of electroweak interactions

Clarify now principal points related to the $`P`$-symmetry and its violation in electroweak interactions, if one corrects the left-right symmetric model so that it would be $`P`$-invariant before the spontaneous symmetry breaking. Note at once, that this correction does not affect values of the cross-sections and the decay probabilities.

Initial $`P`$-invariance and its observed violation is ensured by the local double symmetry produced by the representation $`T`$ = (isotriplet, scalar) $``$ (isotriplet, pseudoscalar) $``$ (isosinglet, scalar) $``$ $`(1,s)(1,p)(0,s)`$ of the group $`G=SU(2)L^{}`$. Corresponding parameters of the secondary symmetry transformations and the gauge fields will be denoted as $`\text{`}^{1s}=\{\theta _j^{1s}\}`$, $`\text{`}^{1p}=\{\theta _j^{1p}\}`$, $`\theta ^{0s}`$; $`𝐁_\mu ^{1V}=\{B_{j\mu }^{1V}\}`$, $`𝐁_\mu ^{1A}=\{B_{j\mu }^{1A}\}`$, $`B_\mu ^{0V}`$ ($`j=1,2,3`$; the indices $`V`$ and $`A`$ mean the polar and axial 4-vectors, respectively).

In the fermionic sector we restrict ourselves by isodoublet consisting of the fields of electronic neutrino $`\nu _e`$ and electron $`e`$. Its seconary symmetry transformations, as well as all anothers, are written in the form, explicitly satisfying the conditions of Definition 1

$$\psi ^{}=\mathrm{exp}\left(\frac{i}{2}\text{ø}\text{`}^{1s}\frac{i}{2}\gamma ^5\text{ø}\text{`}^{1p}+\frac{i}{2}\theta ^{0s}\right)\psi ,$$
(21)
with $`\psi ^T=(\nu _e,e)`$.

Checking that the transformations (21) constitute a group we introduce the above mentioned gauge fields with such phases that in the Lagrangian of their interaction with leptonic field $`\psi `$

$$_{int}=\frac{1}{2\sqrt{2}}\overline{\psi }\left(g_{1V}\gamma ^\mu \text{ø}𝐁_\mu ^{1V}+g_{1A}\gamma ^\mu \gamma ^5\text{ø}𝐁_\mu ^{1A}g_0\gamma ^\mu B_\mu ^{0V}\right)\psi $$
(22)
the coupling constants $`g_{1V}`$, $`g_{1A}`$ and $`g_0`$ are positive. The formulae (21) and (22) give also a knowledge of covariant derivatives corresponding those or others gauge transformations of the secondary symmetry group.

The most general form of the global secondary symmetry transformations of the fields $`𝐁_\mu ^{1V}`$ and $`𝐁_\mu ^{1A}`$ is

$$\left(\genfrac{}{}{0pt}{}{𝐁_\mu ^{1V}}{𝐁_\mu ^{1A}}\right)^{}=\mathrm{exp}\left[i\left(\begin{array}{cc}a_1𝐭& 0\\ 0& a_2𝐭\end{array}\right)\text{`}^{1s}i\left(\begin{array}{cc}0& b_1𝐭\\ b_2𝐭& 0\end{array}\right)\text{`}^{1p}\right]\left(\genfrac{}{}{0pt}{}{𝐁_\mu ^{1V}}{𝐁_\mu ^{1A}}\right),$$
(23)
where $`𝐭=\{t_j\}`$ ($`j=1,2,3`$) are the generators of adjoint representation of the group $`SU(2)`$, and $`a_1`$, $`a_2`$, $`b_1`$ and $`b_2`$ are arbitrary constants. A requirement the transformation (23) to be orthogonal and the Lagrangian (22) to be invariant under the global transformations (21) and (23) is filfiled if and only if $`a_1=a_2=b_1=b_2=1`$ and $`g_{1V}=g_{1A}g`$.

We introduce the Higgs field $`\mathrm{\Phi }`$ consisting of scalar $`\varphi ^{\frac{1}{2}s}`$ and pseudoscalar $`\varphi ^{\frac{1}{2}p}`$ isodoublets and taking some vacuum expectation values of its neutral components with isospin projection $`1/2`$: $`<\varphi _{1/2}^{\frac{1}{2}s}>=v_s`$, $`<\varphi _{1/2}^{\frac{1}{2}p}>=v_p`$. Relative phases of the fields $`\varphi ^{\frac{1}{2}s}`$ and $`\varphi ^{\frac{1}{2}p}`$ are fixed so that the secondary symmetry transformations have the form

$$\left(\genfrac{}{}{0pt}{}{\varphi ^{\frac{1}{2}s}}{\varphi ^{\frac{1}{2}p}}\right)^{}=\mathrm{exp}\left[\frac{i}{2}\left(\begin{array}{cc}\text{ø}& 0\\ 0& \text{ø}\end{array}\right)\text{`}^{1s}\frac{i}{2}\left(\begin{array}{cc}0& \text{ø}\\ \text{ø}& 0\end{array}\right)\text{`}^{1p}\frac{i}{2}\theta ^{0s}\right]\left(\genfrac{}{}{0pt}{}{\varphi ^{\frac{1}{2}s}}{\varphi ^{\frac{1}{2}p}}\right).$$
(24)
The Higgs fields $`\varphi ^{\frac{1}{2}s}`$ and $`\varphi ^{\frac{1}{2}p}`$ are sufficient to reproduce all results of the Weinberg-Salam model in the region of existing energies, except for the generation of fermionic masses. It’s just needed that the relation $`|v_sv_p|<<|v_s+v_p|`$ would be fulfiled and the vacuum expectation value of any other field of Higgs’s type would be much less then $`|v_sv_p|`$. Concerning fermionic masses one can suggest up to some moment that they appear as a result of interaction of fermions with Higgs fields constituting the multiplets $`\{\text{Œ}^{1s},\varphi ^{0p}\}`$ and $`\{\text{Œ}^{1p},\varphi ^{0s}\}`$, whose transformations are generated by the elements of the type (19-20).

In further formulae we neglect the vaccum expectation values af all fields apart from $`\varphi ^{\frac{1}{2}s}`$ and $`\varphi ^{\frac{1}{2}p}`$. From the Lagrangian $`_\mathrm{\Phi }=|𝒟_\mu \mathrm{\Phi }|^2`$, where $`𝒟_\mu `$ is a covariant derivative corresponding to the gauge transformations (24), we get that the electromagtetic field $`A_\mu `$, the fields of light $`W_\mu ^{(1)\pm }`$, $`Z_\mu ^{(1)}`$ and heavy $`W_\mu ^{(2)\pm }`$, $`Z_\mu ^{(2)}`$ intermediate bosons are described by the following relations

$$A_\mu =\frac{1}{\sqrt{g^2+g_0^2}}(g_0B_{3\mu }^{1V}+gB_\mu ^{0V}),$$
(25)
$$W_\mu ^{(i)\pm }=\frac{1}{2}[(B_{1\mu }^{1V}\eta _iB_{1\mu }^{1A})i(B_{2\mu }^{1V}\eta _iB_{2\mu }^{1A})],$$
(26)
$$Z_\mu ^{(i)}=\frac{\alpha _i}{\sqrt{g^2+g_0^2}}(gB_{3\mu }^{1V}g_0B_\mu ^{0V})\beta _iB_{3\mu }^{1A},$$
(27)
$$m_{W^{(i)}}^2=\frac{g^2}{4}|v_s\eta _iv_p|^2,$$
(28)
$$m_{Z^{(i)}}^2=\frac{|v_s|^2+|v_p|^2}{4}\left(g^2+\frac{g_0^2}{2}\eta _i\frac{g_0^2\sqrt{\gamma ^2+1}}{2\gamma }\right),$$
(29)
where

$$\eta _1=1,\eta _2=1,\gamma =\frac{g_0^2(|v_s|^2+|v_p|^2)}{2g\sqrt{g^2+g_0^2}(v_sv_p^{}+v_s^{}v_p)},$$

$$\alpha _i=\sqrt{\frac{1}{2}\left(1\eta _i\frac{\gamma }{\sqrt{\gamma ^2+1}}\right)},\beta _i=\eta _i\sqrt{\frac{1}{2}\left(1+\eta _i\frac{\gamma }{\sqrt{\gamma ^2+1}}\right)}.$$
(30)
From Eqs. (22) and (25-27) we get

$$_{int}=e_0\overline{e}\gamma ^\mu eA_\mu \frac{1}{2\sqrt{2}}\underset{i=1}{\overset{2}{}}[g\overline{\nu }_e\gamma ^\mu (1\eta _i\gamma ^5)eW_\mu ^{(i)+}+\mathrm{H}.\mathrm{c}.$$

$$+\overline{\nu }_e(\alpha _i\sqrt{g^2+g_0^2}\gamma ^\mu \beta _ig\gamma ^\mu \gamma ^5)\nu _eZ_\mu ^{(i)}+\overline{e}(\alpha _i\frac{g^2g_0^2}{\sqrt{g^2+g_0^2}}\gamma ^\mu +\beta _ig\gamma ^\mu \gamma ^5)eZ_\mu ^{(i)}].$$
(31)
It follows from Eqs. (28) and (31) that the observed domination of the left-hand week charged current is possible if and only if $`v_s0`$, $`v_p0`$ and $`\mathrm{arg}(v_p/v_s)\pm \pi /2`$. This means that the physical vacuum does not possess of a definite $`P`$-parity because $`Pv_s=v_s`$, $`Pv_p=v_p`$. This is the first principal point which has not been clear in the left-right symmetric model. It’s interesing, that at the time when the spontaneous breaking of the gauge symmetry and a related to it understanding of physical vacuum was not been formulated yet, Nambu and Jona-Lasinio noted in their paper , ”that the $`\gamma ^5`$ transformation changes the parity of the vacuum which will be in general a superposition of states of opposite parities”.

It follows from Eqs. (26-27) that the fields of all intermediate bosons constitute a superposition of polar and axial 4-vectors, these vectors having an equal weight in the fields of $`W`$-bosons. This is just the second principal point which has not been clear in the left-right symmetric model, though, on our opinion, it has the same powerful significance as the form of weak currents.

If we anywhere considered the intermediate boson masses and the coupling constants of $`Z`$-bosons to fermions as functions of the field values $`v_s`$ and $`v_p`$, then we would get from Eqs. (26)-(31) that under the spatial reflection the $`Z^{(1)}`$-boson coupling constants turn into the $`Z^{(2)}`$-boson coupling constants and vise versa, as well as

$$Pm_{W^{(1)}}^2=m_{W^{(2)}}^2,Pm_{Z^{(1)}}^2=m_{Z^{(2)}}^2,PW_\mu ^{(1)\pm }=(1)^{\delta _\mu }W_\mu ^{(2)\pm },PZ_\mu ^{(1)\pm }=(1)^{\delta _\mu }Z_\mu ^{(2)\pm },$$
(32)
where $`\delta _\mu =0`$ if $`\mu =0`$, and $`\delta _\mu =1`$ if $`\mu =1,2,3`$. This would cause an invariance of the Lagrangian (31).

Thus, the corrected version of the left-right symmetric model of electroweak interactions leads to the logically completed interpretation of the $`P`$-invariance violation, and reveals a full analogy between transformation properties in respect to the spatial reflection of week currents and corresponding gauge fields of the intermediate bosons.

Ackowledgements

I am very grateful to V.I. Savrin, I.P. Volobuev and N.P. Yudin for useful discussions and their support of my work.


format: text

text:
# Comparison of direct and Fourier space techniques in time-dependent density functional theory

## I Introduction

The time-dependent local density approximation has proven to be a useful tool to calculate the optical properties of finite systems such as atoms, molecules, and atomic clusters . The basic equation to be solved is conceptually very simple, little more than the time-dependent Schrödinger equation for a particle in a time-varying external field. Many numerical methods are in use to solve the equations. On the one side there are quantum chemistry methods based on atomic orbital representation for the wave function, and on another side there are methods based on mesh representations. We only consider the latter here, but even in this category there are a number of published techniques. Most fundamentally, the time evolution can be calculated directly or in Fourier space, i.e. in terms of frequencies. The former method is practically a necessity for dealing with very strong external fields and has been applied by two of us (K.Y. and G.B.) for the weak-field response as well. We shall call this approach the “nuclear physics”(NP) method, since the algorithms were originally developed in that field for describing nuclear reactions. The other methods we will consider solve equations in frequency space. The method described in ref.  had its origins in condensed matter theory and uses Fourier representation for both space and time; we shall call this the “condensed matter” (CMP) method. We also comment on ref.  which uses Fourier space for the time but a real space mesh for the spatial dependence. Here the problem is cast into a matrix diagonalization in the particle-hole representation; we shall call it the diagonalization method.

In this work we will compare the CMP code and the NP code for a specific system and present arguments for the scaling properties of the respective algorithms for larger systems. The system we choose to study is the atomic cluster Na8, and in particular the surface plasmon excitation which is seen as a strong peak at 2.5 eV excitation. The TDLDA is not an exact theory and it predicts a excitation energy at about 2.7 eV. We shall demand of both methods that they achieve within 0.1 eV of the converged value. It makes little sense to calculate to higher precision in view of the intrinsic limitations of the theory.

We shall now describe the various methods from a computational point of view. We shall use the symbol $`N`$ with subscripts for quantities that scale roughly as the size of the physical system under study, and $`M`$ for quantities that may be large but are independent of the size of the system. Important quantities common to the two codes are the number of electrons $`N_e`$ and the number of mesh points, $`N_G`$ and $`N_R`$ for real space and reciprocal space, respectively. Additional quantities that play a role are the number of frequencies to be calculated $`M_\omega `$, and the number of time steps to evolute the wave function in the real-time method, $`M_T`$. Also, in methods that rely on sparse matrix multiplication, we need the number of nonzero entries in a row of the Hamiltonian, $`M_H`$, and in iterative methods to solve large matrix equations we need the number of iterations for convergence, $`M_{it}`$. Finally, the response function method usually requires a sum over unoccupied states, $`N_c`$. This notation is summarized in Table I.

We will use same energy functional for all methods, so the choice of specific functional is not an issue in comparing the methods. As is commonly done, we calculate only the dynamics of the valence electrons. The core electrons are frozen and their presence is treated by using a pseudopotential to describe the ionic potential. We use the pseudopotential construction of Troullier and Martins , taking the nonlocal part by the method of Kleinman and Bylander  and including partial core corrections for the exchange-correlation energy. In this method, the local pseudopotential is fixed to the value in a particular angular momentum channel, and a nonlocal correction is made for other channels. Here we use the $`l=1`$ potential as the local potential, and apply the nonlocal correction to the $`l=0`$ and $`l=2`$ channels. The electron-electron interaction is taken in the simple local-density approximation (LDA) given by Perdew and Zunger . More complicated functionals have better predictive power for ground state properties , but give only small improvement to the optical response of neutral molecules. The proper description of the 1/r asymptotic behavior of the potential is going to be very important to describe charged systems, however for the aim of the present work this LDA deficiency is not relevant.

The geometry of the Na8 cluster was computed in ref., and the lowest energy structure found to be the bicapped octahedron (D2d symmetry). We use this structure in our comparison here. It has an average Na-Na bond length of 3.38 Å and a slight deviation from the spherical symmetry. This leads to a polarizability tensor with two different components and two close-lying peaks are obtained in the photoabsorption cross section.

## II Theoretical methods

Before describing in detail each of the two methods for representing the wave functions (direct and Fourier space), we need to comment on the choice of the spatial cell size and mesh size as well as the time/frequency parameters (all are summarized in Table I).

Since the wave functions are sensitive to boundaries, the calculations must be made in a volume several Angstroms larger than the size of the molecule or cluster. Using both methods, we determined how large a volume is needed to achieved 0.1 eV accuracy on the various excitation energies of interest in the system. We found that this is achieved in a spherical volume of radius $`R=8`$ Å using the NP code, and in a simple cubic supercell of side 12.7 Å using the CMP code. These have nearly the same volume, and thus the same average distance from the cluster to the boundary. We have checked the convergence of the results by increasing the volume to a sphere of 12 Å radius. The value of the plasma frequency is reduced by a maximum of 0.1 eV, that is, within the required accuracy.

We have used an uniform spatial grid with $`\mathrm{\Delta }x=0.5`$ Å spacing. This corresponds to a plane-wave cutoff energy of 6 Hartrees in the Fourier space method (see below). Within this parameters, a stable time-step to perform the time-evolution in the NP method is $`\mathrm{\Delta }t=0.003\mathrm{}/eV<<\mathrm{}(\mathrm{\Delta }x)^2/m`$. The required 0.1eV accuracy in energy is obtained for total simulation times of 10 $`\mathrm{}`$/eV. Similarly in the Fourier space method we have taken a uniform frequency grid of $`M_\omega =100`$ between 0 and 5 eV. Note that if the response is required for larger frequencies we need to increase the number of points. The whole response is obtained at once in the time evolution method (unless up to energies of the order of $`(\mathrm{\Delta }t)^1`$). This is a great advantage when the whole response is needed.

### A  NP method

This method uses a direct solution of the time-dependent single-electron Schrödinger equation,

$$i\mathrm{}\frac{\varphi _i(𝐫,t)}{t}=H_{KS}(t)\varphi _i(𝐫,t)(\mathrm{i}=1\mathrm{}\mathrm{occ}.)$$
(1)
where $`H_{KS}`$ is the Kohn-Sham Hamiltonian operator

$$H_{KS}(t)=\frac{\mathrm{}^2}{2m}^2+V_{ion}(𝐫)+e^2d^3r^{}\frac{n(𝐫^{},t)}{|𝐫𝐫^{}|}+V_{xc}(𝐫,t)$$
(2)
and $`n`$ is the time-dependent electron density $`n(𝐫,t)=_{i=1}^{occ}\varphi _i^{}(𝐫,t)\varphi _i(𝐫,t)`$. In the solution of this equation in the spatial and time variables following the algorithm of ref., there are two time-consuming operations. One is multiplying the single-electron Hamiltonian operator by the vector representing the wave function. The dimensionality of the vector is the number of mesh points $`N_R`$ times the numbers of electron orbitals $`N_e`$. The operator is a sparse matrix with $`M_H`$ nonzero elements per row. Thus the basic operation requires about $`N_eN_RM_H`$ complex floating point operations. The time evolution operator in the NP code is implements by a power series expansion of the exponential operator $`\mathrm{exp}(iH\mathrm{\Delta }t)`$ to fourth order. A predictor-corrector cycle requires two such operations. Thus the method requires 8 Hamiltonian multiplications per time step. Thus for $`M_T`$ time steps the total number of floating point operations is given by

$$NPFPO:10N_eN_RM_HM_T$$

The sparseness of the Hamiltonian matrix in a real space formulation is determined by the finite difference formula for kinetic energy (nine-point formula in our case); and by the nonlocal-projection parts of the potential. In total we have a number of non-zero elements of the each Hamiltonian row $`M_H100`$ for the grid parameters used for Na8.

The other time-consuming part of the NP algorithm is solving the Poisson equation, which must be done twice at each time step. The NP code uses a multipole expansion combined with a relaxation method to deal with the higher multipoles. It is hard to estimate the scaling properties of this part, but in the present study this part of the computation takes 1.5 times as many operations as the Hamiltonian multiplication operation. We shall assume the same factor for estimating the scaling properties of the algorithm. In principle, the Poisson equation can be solved by methods that are of order $`N_R`$ or $`N_R\mathrm{log}N_R`$, as multigrid or fast-Fourier transformation, so this part should not dominate for large system.

Storage requirements are small: the vector wave function plus $`V_{Hartree}`$ and $`V_{ion}`$ local potentials in Hamiltonian, charge densities and some intermediate arrays. $`V_{Hartree}`$ requires a slightly larger volume because of the way the Poisson equation is solved.

$$\mathrm{NP}\mathrm{storage}:N_R(N_e+4.5)$$

This NP method is ideal to be combined with molecular dynamics simulations for the ions because it uses only ground-state occupied information and would scale roughly linearly with the number of atoms in the system. There is not so much book-keeping as in the usual perturbative formalism (no need for storing the large set of unoccupied wave-functions and the large dielectric matrices).

### B CMP method

Here the basic object of the calculation is the linear response to an external field of some frequency $`\omega `$. The linear response matrix $`\chi `$ is constructed in momentum space with the following matrix inversion

$$\chi =(1\chi _0K)^1\chi _0$$
(3)
where the independent particle response $`\chi _0`$ and the interaction $`K`$ are matrices defined as follows. The $`\chi _0`$ has elements $`G,G^{}`$ given by

$$\chi _0(𝐆,𝐆^{},\omega )=\frac{1}{\mathrm{\Omega }}\underset{kj}{}(f_kf_j)\frac{k|e^{i𝐆𝐫}|ii|e^{i𝐆^{}𝐫}|k}{\omega ϵ_j+ϵ_i+i\eta }$$
(4)
where $`\mathrm{\Omega }`$ denotes the unit-cell volume, $`i,k`$ label Kohn-Sham eigenfunctions and $`ϵ_k`$ and $`f_k`$ are the corresponding eigenenergies and occupancy factors. The sum goes over $`N_e`$ occupied orbitals and $`N_c`$ empty orbitals. The interaction $`K`$ is the Fourier transform of the electron-electron interaction in the Kohn-Sham equation, which is given in coordinate space by

$$K(𝐫,𝐫^{})=\frac{e^2}{|𝐫𝐫^{}|}+\frac{\delta V_{xc}(𝐫)}{\delta n(𝐫^{})}$$
(5)
We now describe the computation starting with the Kohn-Sham wave functions and energies in a momentum space representation. To evaluate the independent particle response $`\chi _0`$ in eq.( 4), one first calculates the particle-hole matrix elements of the momentum operator and stores them in a table (or in disk). This computational effort is of the order of $`N_eN_cN_G^2`$ operations, and the table size to be stored is $`N_eN_cN_G`$ complex numbers. Then the evaluation of eq. (4) requires $`N_G^2`$ matrix elements to be calculated, each requiring particle-hole summation, to give $`2N_G^2N_eN_c`$ operations for each frequency. If one were to make full space calculation, the number of empty orbitals summed in eq. (4) would be of the same order as the dimensionality of the space. However, the number of empty orbitals can be severely truncated without effecting the long-wavelength dipole response. In the example, we find $`N_c=320`$ is adequate, which is more than an order of magnitude smaller than the size of the space and corresponds to include unoccupied states up to 20 eV above the highest-occupied orbital. This is a reasonable approximation as we are interested only in getting the optical spectra for excitation energies below 10 eV. This approximation is an important saving in building up the response matrix.

One also truncates the calculation of the response matrix in another way. We have also assumed that the off-diagonal elements of the response function are zero for G-vectors outside an sphere of 1.25 Å (that is to consider $``$ 3200 points in the G-space). This corresponds to reducing the number of matrix elements to be computed and stored to $`N_G(N_G+4)/18`$. Note that the necessity to store the $`N_G^2`$ matrix puts a higher demand on the computer memory than the NP method. The memory required to store the $`N_G^2`$ complex, double-precision numbers in the example problem is $`164`$Mb.

There are now three steps to evaluate eq. (3), two matrix multiplications and a matrix inversion. The matrices are not sparse, so the matrix multiplications each cost $`2(N_G/3)^3`$ arithmetical operations §§§A small technical point should be mentioned, associated with the divergence of the Coulomb interaction at $`𝐆=𝐆^{}`$. This is dealt with by taking a numerical limit as $`|𝐆𝐆^{}|0`$, and this adds about 10% to the number of operations for computing the matrix product. . The matrix inversion is of the same order, requiring $`(N_G/3)^3`$ operations. The total is $`5(N_G/3)^3.`$ These represent the most computationally demanding steps in the CMP method, given the truncation in the $`N_c`$. The computed $`\chi `$ is next transformed to the coordinate space representation. Using the fast Fourier transform, this takes $`N_G^2\mathrm{log}N_G`$ operations. The dynamical polarizability can be now computed from $`\alpha (\omega )=V_{ext}\chi V_{ext}`$ as a matrix times vector multiplication. From this one can easily extract the photoabsorption cross section $`\sigma (\omega )=\frac{4\pi \omega }{c}Im\alpha (\omega )`$.

Then the total computational effort in the CM method is:

$$\mathrm{CM}\mathrm{FPO}:M_\omega (N_cN_e(N_G/3)^2+5(N_G/3)^3)$$

with the last term dominant. The storage requirements for all the occupied and unoccupied wave functions plus the whole complex response matrix is

$$\mathrm{CM}\mathrm{storage}:(N_e+N_c)N_G+2(N_G+N_G^2/9)+N_cN_eN_G/3$$

To achieve the targeted energy convergence with this algorithm, the momentum space mesh was chosen to correspond to a simple cubic supercell of $`L=12.7`$ Å on a side. This implies that the mesh spacing in momentum space is $`\delta k=2\pi /L=0.137`$ Å. The momentum space representation takes all the points within a sphere of radius $`k_{max}=1.83`$ Å (that corresponds to a plane-wave cutoff energy of 12 Ry). The size of the vector in the momentum representation is thus $`N_G=4\pi (k_{max}/\mathrm{\Delta }k)^3/310,000`$. Note that this is slightly smaller than the number required for the coordinate space representation, however we need to stress that a larger number of G-vectors are needed to describe the action of the potential on a wave-function ($`V\psi `$ corresponds to a convolution in Fourier space). Finally, an additional numerical parameter is the imaginary part of the frequency $`\eta `$, which we have taken as $`\eta =0.05`$ eV to produce a resolution of 0.1 eV in the spectral features.

In the discussion below we have not include the computational requirements to perform the ground state calculations, occupied and unoccupied orbitals. This could be a major storage bottle-neck for very large systems as the calculation of a large set of unoccupied wave functions has a cubic scaling of the number of atoms in memory and computing time. In the present calculation this initialization process takes 10% of the total computational time.

### C Other methods

We mention here two other methods from a computational point of view. Since we have not carried out numerically computations on our test problem with these methods, the discussion will be brief.

#### 1 Modified Sternheimer method

The modified Sternheimer method was first applied to the time-dependent Kohn-Sham equation for atomic excitations, and has since been applied to the dielectric response of crystals using the momentum space representation  and to the finite system C60  using the coordinate space representation. Here one solves an inhomogeneous equation for the perturbed wave functions $`\varphi _i^\pm `$ using an iterative method. The perturbation is a sinusoidal potential field combining the external field $`V_{ext}`$ and the internal field from the time-varying electron density. The equations are

$$(ϵ_IH_{KS}^0\pm \omega +i\eta )\varphi _i^\pm =\widehat{P}V_i$$
(6)
where

$$V_i=(V_{ext}+K\delta n)\varphi _i$$

and

$$\delta n=\mathrm{Re}\underset{i}{}\varphi _i(\varphi _i^++\varphi _i^{}).$$
(7)
$`\widehat{P}`$ is a projection operator removing occupied orbitals. In ref. , the two equations are constructed in coordinate space and solved with a double iteration. One makes a guess for the density $`\delta n`$, and solves eq. (6) by the conjugate gradient method. $`\delta n`$ is refined from the resulting $`\varphi _i^\pm `$ again with the conjugate gradient method, and the process is repeated to convergence. The numerical cost will thus depend largely on the cost of the Hamiltonian operation which is $`M_HN_RN_e`$ in coordinate space, and the number of iterations $`M_{it}`$ required to get a converged solution. Remembering also that frequency space methods need $`M_\omega `$, the number of frequencies to be examined, the computational cost of this method is

$$\mathrm{Modified}\mathrm{Sternheimer}(\mathrm{real}\mathrm{space}):M_\omega M_{it}M_HN_RN_e$$
(8)
The method can be used in this form for nonresonant frequencies, but near the eigenfrequencies the nearby singularities in eq. (6) must be removed for the conjugate gradient method to converge. Thus this method would be similar to methods utilizing the particle-hole representation in needing a considerable number of the wave functions and eigenenergies of unoccupied states. The singularities are removed by projecting on the unoccupied wave function subspace the right hand side of eq. (6),

$$V_i^{}=V_i\underset{j}{}\varphi _i(\varphi _i,V_i).$$

The desired wave functions $`\varphi _i^\pm `$ are obtained from the projected solutions $`\varphi _i^\pm `$ by

$$\varphi _i^\pm =\varphi _i^\pm +\underset{j}{}\frac{\varphi _i(\varphi _i,V_i)}{ϵ_jϵ_i\omega i\eta }.$$

It is difficult to give an a priori estimate of $`M_{it}`$ or its size-scaling properties (although with our notation we have assumed that it does not grow with $`N`$). Unfortunately, our implementation of eq. (7) still left the convergence somewhat erratic. Typically it takes of the order of $`M_{it}1000`$ iterations of the double loop to get convergence. Thus it would require some improvement of the algorithm to make it attractive to apply to large systems.

The momentum space implementation of the modified Sternheimer method is similar. This method also needs the conditioning step for convergence of the CG iteration. The main difference is in the Hamiltonian multiplication, which here requires $`2(N_G/3)^3`$ operations as discussed in Sect. IIB. Thus the total is

$$\mathrm{Modified}\mathrm{Sternheimer}(\mathrm{momentum}\mathrm{space}):\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}M_\omega M_{it}(N_G/3)^3$$
(9)
Because the Hamiltonian operation is more costly in momentum space, this method is probably not competitive to the others, unless it were the case that the convergence of the iteration were intrinsically much more reliable.

#### 2 Diagonalization method

The frequency-space methods discussed so far have relied in some way on operator inversion. It is also possible to cast the problem as one of matrix diagonalization. This method was applied to cluster excitations in the TDLDA by Vasiliev et al. . The authors start from a basis in coordinate space and construct Kohn-Sham orbitals for both occupied and empty states as is done in the CMP method, but representing the orbitals in coordinate space mesh, as in the NP method. The storage requirement for the orbitals is $`(N_c+N_e)N_R`$, which is larger than in the NP method but smaller than in the CMP method.

The next step of the calculation is to construct the matrix to be diagonalized. The eigenvalue equation to be solved is

$$\mathrm{𝐑𝐅}_n=\omega _n^2𝐅_n$$
(10)
where $`𝐅_𝐧`$ are the eigenvectors and $`𝐑`$ is a matrix. Its elements are

$$R_{\alpha ,\alpha ^{}}=(ϵ_iϵ_j)^2\delta _{\alpha ,\alpha ^{}}+2\sqrt{(ϵ_iϵ_j)(ϵ_i^{}ϵ_j^{})}K_{\alpha ,\alpha ^{}}$$
(11)
where the indices $`\alpha =(ij),\alpha =(i^{},j^{})`$ label combinations of unoccupied orbitals $`i`$ and occupied orbitals $`j`$. The interaction matrix elements $`K_{\alpha ,\alpha ^{}}`$ are simply the particle-hole matrix elements of the residual interaction, eq.(5). There is a substantial computational cost in construct the interaction matrix $`K`$. A straightforward transformation from the coordinate space to the particle-hole representation requires $`N_R^2N_e^2N_c^2`$ operations for the Coulomb interaction. However, this is reduced considerably by using an efficient method to solve the Poisson equation. For example, using the fast Fourier transform one may find the Coulomb field for a given particle-hole state taking only $`N_R\mathrm{log}N_R`$ operations. Saving the Coulomb field in the coordinate representation, the matrix element to a given final state takes $`N_R`$ operations. The effort of solving the Poisson equation is thus distributed over the number of final states, and the operations to construct the full matrix has a leading dependence $`N_RN_e^2N_c^2`$, the scaling appropriate for the local part of the interactionHowever, in the implementation of ref. , the Poisson solver in fact is the most costly operation.. Once the matrix is constructed, the diagonalization requires $`(N_cN_e)^3`$ operations. However, taking the $`N`$ values from Table I, the matrix diagonalization effort is small compared to that needed to construct the matrix. We have therefore taken that step to assign this method’s size scaling in Table II.

## III Numerical results

We will discuss in detail the physical quantities computed in the NP and CMP methods and refer to  for the results using the diagonalization method. We want to stress that the three approaches must give the same values if the numerical parameters are chosen with fine enough grids and large enough cutoffs to get converged results.

With the parameter sets chosen for the two methods, the results are quite similar. In Table IV we show calculated Kohn-Sham energies and the surface plasmon energy. The first entry $`ϵ_1`$ is the Kohn-Sham eigenvalue of the most bound orbital. The absolute energies have no significance in the supercell method, because the absolute Coulomb potential is undefined. Therefore, for this entry we give the value from the NP code and set the scale of the CMP energies at that value. The next three rows correspond to the other bound orbitals use the $`G=0`$ point of the Brillioun zone for the CMP values. We can see that the methods agree to within less than 0.1 eV. The next entry is the lowest unoccupied orbital. This is significantly different for the two methods. This orbital has sufficient extension to have its energy sensitive to the boundary, which of course is different for the two methods. We confirm the boundary sensitivity in the CMP code by calculating the energies at other points in the Brillouin zone. Differences are less than 0.1 eV for occupied orbitals, but reach 0.2 eV for the lowest unoccupied orbital. This last point indicates the fact that the empty orbitals are more sensitive to the boundary conditions and in the periodic supercell they feel the potential from the other clusters.

We have also presented in Table IV the results of the NP method The plasmon frequency is sensitive to the core-exchange correction at the level of 0.1 eV. We have included that correction in $`H_{KS}`$ it improves the description of the structural properties of Na metal. We note that the result without core corrections (2.89 eV) it is very close to the jellium value (2.9 eV). . We have also checked the convergence of the plasmon frequency with respect to the cell size and found that this value is converged to less than 0.01eV for a sphere of R=12 Å. The fully converged value in the NP method is 2.65 eV. The difference with the experimental value of 2.53 eV can be attributed to deficiencies in the LDA approximation as well as for finite temperature effects in the experiments.

In Table V we summarize the results for the static averaged electrical polarizability of Na8 obtained by the different methods. The agreement among the different approaches is very good and the remaining difference with experiments can be again assigned to core polarization, exchange-correlation and temperature effects. These effects tends to increase the polarizability bringing the computed values close to the experiments .

## IV Conclusions

In the theory of electronic excitations of finite many-electron systems, the time-dependent Kohn-Sham equation with an adiabatic local density approximation for the interaction energy function offers an attractive compromise towards the goals of accuracy and computational practicality. But even within the TDLDA scheme there are several methods in use, and our purpose was to compare them on the same footing by applying them to the same physical problem, and demanding the same accuracy. The goal is to gain a general understanding of the numerical resources (total numbers of arithmetic operations and computer memory) required by the different methods. One can then extrapolate to large systems and make a judgment on which methods offer the best prospects.

We have only considered methods based on a grid representation of the electron wave functions, and have concentrated on two algorithms, the NP method in real time and real space, and the CMP method in Fourier transformed time and space.

We chose to study the response of the Na8 cluster around the surface plasmon excitation energy. The two methods turned out to have comparable requirement on arithmetic operations. However, it should also be noted that the computational work increases with the range of frequencies that one studies in the CMP method, but not in the NP method. With latter, the entire response is obtained from a single calculation.

In comparing the two methods to ascertain their scaling with the size of the system $`N`$, we have deliberately ignored the first task in either method, the construction of the eigenstates of the static Kohn-Sham operator. In the NP method only the occupied orbitals are needed, but in the CMP method one also needs a large number of unoccupied orbitals as well. Their calculation scales like $`N_e^3`$ in principle, but in practice this phase of the computation is short compared to the dynamic calculation and so we ignore it. Let us now compare the scalings by taking the expressions in Table II, dropping the subscripts on the $`N`$ quantities. The NP method thus scales as $`N^2`$. This behavior was also found studying the excitations of long carbon molecules. The CMP method has a poorer scaling behavior, namely $`N^3`$. We also considered two other methods without however examining them in as much detail. In principle, the modified Sternheimer method in coordinate space can achieve $`N^2`$ scaling without the cost of the large $`M_T`$ factor of the real-time method. However, we did not find a reliably converging iteration procedure to solve the basic inhomogeneous linear equation set. The final method we discussed, the diagonalization method using real space and Fourier time, seems to have a poorer $`N`$-scaling than the others, but may be advantageous in some circumstances (see below).

Besides arithmetic operations, storage can play a role in the practicality of the different algorithms for large systems. Here we find that the storage requirements are grossly different for the NP and CMP methods, favoring the NP approach. From Table II, it has a $`N^2`$ scaling while the CMP method has an $`N^3`$ behavior. This is already significant in the Na8 system we studied, as may be seen from Table III.

Thus our results favor the real-time and real-space methods, offering economy in both storage and arithmetic operations. However, there are a number of caveats. We have not considered the suitability of the different algorithms for parallel computing. In a parallel computing environment, the frequency-space methods gain favor because the $`M_\omega `$ factor can be trivially absorbed in the parallel processing. In addition, the diagonalization method can benefit from the parallel computation of different rows of the matrix. Also the sparseness of the Hamiltonian matrix is important for the real space method; this would be lost if for example the energy functional used the full Fock exchange interaction.

Finally, we mention two nonnumerical benefits of the real-time method: as was said earlier, it is nonperturbative and therefore allows effects of large fields to be calculated with the same effort. And it uses the same energy functional (permitting the program to call the same subroutine) for the dynamic calculation as for the static calculation to prepare the ground state.

## V Acknowledgment

We are grateful to J. Chelikowsky and I. Vasiliev for communications and providing us with their computer code. This work was supported by the Department of Energy under Grant FG06-90ER-40561, by the DGES (PB98-0345) and JCyL (VA28/99), and by he Grant-in-Aid for Scientific Research from the Ministry of Education, Science and Culture (Japan), No. 11640372. AR acknowledges the hospitality of the Institute for Nuclear Theory where this work was started and the computer time provided by the C4 (Centre de Computació i Comunicacions de Catalunya).


text:
# Hermitian quark mass matrices with four texture zeros

## I Introduction

The standard model is known to have a large number of parameters. However, most of them are contained in the pair of quark mass matrices, $`M_U`$ and $`M_D`$. Of the thirty-six parameters in these two complex matrices, only ten \[six quark masses, three Cabibbo-Kobayashi-Maskawa (CKM) angles and a CP phase\] are physical. This redundancy arises from unphysical, right-handed (RH) rotations in both the $`U`$\- and $`D`$-sectors, in addition to common left-handed (LH) rotations, which cancel out in the CKM matrix, $`V^{\text{CKM}}=V_U^{}V_D`$, where $`V_U^{}M_UU_U=M_U^{\text{diag}}`$ and $`V_D^{}M_DU_D=M_D^{\text{diag}}`$. To eliminate these unphysical degrees of freedom, it was pointed out earlier  that RH rotations can be used to reduce both $`M_U`$ and $`M_D`$ to the upper triangular form, which in the hierarchical basis exhibits most clearly the quark masses and left-handed (LH) rotation angles,

$$M_UU_U^R=M_U^t=\left(\begin{array}{ccc}X& X& X\\ 0& X& X\\ 0& 0& X\end{array}\right),$$
(1)
$$M_DU_D^R=M_D^t=\left(\begin{array}{ccc}X& X& X\\ 0& X& X\\ 0& 0& X\end{array}\right).$$
(2)
The matrices $`U_U^R`$ and $`U_D^R`$ can be constructed explicitly as $`R_{13}R_{23}R_{12}`$, such that $`R_{13}`$ eliminates the (3,1) element of $`M`$, etc. Another way to obtain the triangular matrix elements directly is by a geometric argument. Let us write a general (real) matrix $`M`$ in the form:

$$M=\left(\begin{array}{ccc}a_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right).$$
(3)
We may regard the three rows as components of the vectors $`\stackrel{}{a},\stackrel{}{b}`$ and $`\stackrel{}{c}`$. Then the transformation into the triangular form amounts to a rotation into a new coordinate system where the axes are aligned in the directions of the unit vectors: $`(\stackrel{}{b}\times \stackrel{}{c})/|\stackrel{}{b}\times \stackrel{}{c}|`$, $`\stackrel{}{c}\times (\stackrel{}{b}\times \stackrel{}{c})/(|\stackrel{}{c}||\stackrel{}{b}\times \stackrel{}{c}|)`$, and $`\stackrel{}{c}/|\stackrel{}{c}|`$, respectively. It follows that

$$M^t=\left(\begin{array}{ccc}\frac{\stackrel{}{a}(\stackrel{}{b}\times \stackrel{}{c})}{\stackrel{}{b}\times \stackrel{}{c}}& \frac{(\stackrel{}{a}\times \stackrel{}{c})(\stackrel{}{b}\times \stackrel{}{c})}{\stackrel{}{c}\stackrel{}{b}\times \stackrel{}{c}}& \frac{\stackrel{}{a}\stackrel{}{c}}{\stackrel{}{c}}\\ 0& \frac{\stackrel{}{b}\times \stackrel{}{c}}{\stackrel{}{c}}& \frac{\stackrel{}{b}\stackrel{}{c}}{\stackrel{}{c}}\\ 0& 0& \stackrel{}{c}\end{array}\right).$$
(4)
For complex vectors $`\stackrel{}{a},\stackrel{}{b}`$ and $`\stackrel{}{c}`$, each transforms as a 3 under RH SU(3) rotations. From $`[\mathrm{𝟑}\times \mathrm{𝟑}]_{anti}\overline{\mathrm{𝟑}}`$, $`[\mathrm{𝟑}\times \mathrm{𝟑}\times \mathrm{𝟑}]_{anti}\mathrm{𝟏}`$, etc., we readily find:

$$M^t=\left(\begin{array}{ccc}\frac{\stackrel{}{a}(\stackrel{}{b}\times \stackrel{}{c})}{\stackrel{}{b}\times \stackrel{}{c}}& \frac{(\stackrel{}{a}\times \stackrel{}{c})(\stackrel{}{b^{}}\times \stackrel{}{c^{}})}{\stackrel{}{c}\stackrel{}{b}\times \stackrel{}{c}}& \frac{\stackrel{}{a}\stackrel{}{c^{}}}{\stackrel{}{c}}\\ 0& \frac{\stackrel{}{b}\times \stackrel{}{c}}{\stackrel{}{c}}& \frac{\stackrel{}{b}\stackrel{}{c^{}}}{\stackrel{}{c}}\\ 0& 0& \stackrel{}{c}\end{array}\right).$$
(5)
Note that this construction is unique up to a diagonal phase matrix multiplying from the right: $`M^tM^tP`$, $`P=`$ diagonal phase matrix.

After we have transformed $`M_U`$ and $`M_D`$ to the upper triangular form, we may reduce them further by common LH rotations on both. Since there are three degrees of freedom in these rotations, we can use them to generate three additional zeros in the pair. Note that this result is only approximately correct. The LH rotations will generate small elements in the lower-left part of the matrices, which can be removed by even smaller RH rotations. Thus, the process actually consists of a sequence of rotations $`R_i`$, which can be schematically expressed as $`R_3((\frac{m_1}{m_2})^2V_{12},(\frac{m_2}{m_3})^2V_{23},(\frac{m_1}{m_3})^2V_{13})R_1(V_{12},V_{23},V_{13})MR_2(\frac{m_1}{m_2}V_{12},\frac{m_2}{m_3}V_{23},\frac{m_1}{m_3}V_{13})`$.

The above argument provides an algorithm to reduce any pair of $`M_U`$ and $`M_D`$ into upper triangular forms, with the minimal nine non-vanishing elements between the two matrices. As was shown earlier, the most interesting property of this reduction is that to a good approximation, all of the remaining nine non-vanishing matrix elements in $`M^U`$ and $`M^D`$ are physical and are simple products of a quark mass and some CKM matrix elements. This method can thus provide us with a powerful tool to assess the viability of any proposed mass matrices. It was applied to hermitian mass matrices with five texture zeros, and a unique pair is identified as most favorable to present data .

In this paper, we will apply the same technique to the study of hermitian mass matrices with four texture zeros, which generally entail one testable relation for each pair. The study of 4-texture-zero hermitian mass matrices has been a subject of interest recently \[3-8\]. Attentions have been focused on specific models or structures, which can lead to predictive values for the CKM matrix elements. However, due to the complexity involved in the analysis, no complete investigation has been performed so far. Using the triangular basis can greatly reduce the amount of work involved in analyzing the viability of 4-zero textures. With this efficient tool, we are able to provide the first complete study of all viable four-zero texture pairs that exhibit hierarchical structure. The results are presented in Tables I through IV.

## II Hermitian Mass Matrices and Texture Zeros

We start with quark mass matrices in the triangular form. Using the mass relations $`m_u:m_c:m_t\lambda ^8:\lambda ^4:1`$ and $`m_d:m_s:m_b\lambda ^4:\lambda ^2:1`$, and the CKM elements $`V_{us}=\lambda `$, $`V_{cb}\lambda ^2`$, $`V_{ub}\lambda ^4`$, and $`V_{td}\lambda ^3`$, the properly normalized Yukawa matrices for $`U`$ and $`D`$ can be put into the most general triangular form,

$$T^U=\left(\begin{array}{ccc}a_U\lambda ^8& b_U\lambda ^6& c_U\lambda ^4\\ 0& d_U\lambda ^4& e_U\lambda ^2\\ 0& 0& 1\end{array}\right),T^D=\left(\begin{array}{ccc}a_D\lambda ^4& b_D\lambda ^3& c_D\lambda ^3\\ 0& d_D\lambda ^2& e_D\lambda ^2\\ 0& 0& 1\end{array}\right).$$
(6)
Here, all of the coefficients ($`a,b,\mathrm{}`$) are assumed to be of order unity or less. Without loss of generality, we also take the diagonal parameters (i.e. $`a_U`$, $`a_D`$, $`d_U`$, $`d_D`$) to be real throughout this paper. We offer the following remarks.

1. The hierarchical structures of $`T^U`$ and $`T^D`$ are manifest. As was shown earlier, we can diagonalize them approximately with only a LH rotation, $`M_{U,D}^{\text{diag}}V_{U,D}^{}T^{U,D}`$, with the matrix elements $`(V_{U,D})_{ij}T_{ij}^{U,D}/T_{jj}^{U,D}`$ ($`i<j`$).
2. The diagonal elements of $`T^U`$ and $`T^D`$ are essentially the quark masses. The CKM matrix elements are simply given by $`V_{ij}^{\text{CKM}}(V_D)_{ij}(V_U)_{ij}`$,   ($`i<j`$). These simple relations are lost when we go to other bases, including the hermitian basis. Triangular matrices thus stand out as the unique basis whose matrix elements correspond directly to quark masses and LH rotation angles. This feature makes the triangular form especially useful for analyzing the texture of quark mass matrices.
3. To avoid fine tuning in generating the CKM mixings, naturalness criteria  has been imposed when writing the above triangular matrices. This simply implies that, for the LH rotation angles, $`|V_{U,D}|_{ij}|V_{ij}^{\text{CKM}}|\times 𝒪(1)`$. Note that this condition can always be implemented by applying a common LH rotation on $`T^U`$ and $`T^D`$.
4. We can put $`T^U`$ and $`T^D`$ in the minimal parameter basis (m.p.b.), with only three nonzero off-diagonal elements between them. This is achieved from Eq. (6) by common LH rotations to generate three zeros. For example, the (1,3) element of $`T^D`$ can be set to zero by a common LH rotation $`R_{13}(c_D\lambda ^3)`$. This method produces ten pairs of triangular matrices, as listed in Table I of Ref.. The elements of each pair consist of a simple product of quark mass and CKM elements. Nothing is unphysical in the m.p.b.. As a result, the viability of a given texture can be readily tested by turning it into one of the ten.

Turning now to hermitian matrices, the consequences of their texture zeros can be easily understood in terms of the triangular parameters of Eq. (6). For this purpose, one can start from Eq. (6) and, by a RH rotation, generate the corresponding hermitian form,

$`Y^U`$ $`=`$ $`\left(\begin{array}{ccc}\left(a_U+c_Uc_U^{}+\frac{b_Ub_U^{}}{d_U}\right)\lambda ^8& (b_U+c_Ue_U^{})\lambda ^6& c_U\lambda ^4\\ (b_U^{}+c_U^{}e_U)\lambda ^6& (d_U+e_Ue_U^{})\lambda ^4& e_U\lambda ^2\\ c_U^{}\lambda ^4& e_U^{}\lambda ^2& 1\end{array}\right)\times (1+𝒪(\lambda ^4)),`$ (10)
$`Y^D`$ $`=`$ $`\left(\begin{array}{ccc}\left(a_D+\frac{b_Db_D^{}}{d_D}\right)\lambda ^4& b_D\lambda ^3& c_D\lambda ^3\\ b_D^{}\lambda ^3& d_D\lambda ^2& e_D\lambda ^2\\ c_D^{}\lambda ^3& e_D^{}\lambda ^2& 1\end{array}\right)\times (1+𝒪(\lambda ^2)).`$ (14)
We immediately see that whereas off-diagonal hermitian zeros (e.g. $`Y_{13}^U`$, $`Y_{23}^U`$ or $`Y_{ij}^D`$ ($`ij`$)) have a one-to-one correspondence to zeros in the triangular form, the diagonal hermitian zeros (i.e. $`Y_{11}^U`$, $`Y_{22}^U`$, and $`Y_{11}^D`$) entail definite relations between diagonal and off-diagonal triangular parameters (i.e. quark masses and LH rotation angles). Note that these simple relations between the hermitian and triangular forms are approximate, and follow from the hierarchical structure in Eq. (6).

We are thus led to the following procedure in analyzing hermitian texture zeros. First, we obtain all possible textures with certain number of zeros for $`Y^U`$ and $`Y^D`$ by referring to Eq. (10) and Eq. (14), as listed in Tables I and III. Then we can list all the 4-texture-zero pairs. Second, each pair is transformed into one of the ten triangular forms in the m.p.b., using common LH rotations when necessary. Possible relations between quark masses and mixing can then be read-off for each pair by referring to Table I of Ref. . Finally, we arrive at all the viable textures after confronting the prediction of each pair with data. In the next section, we apply this procedure to the analysis of hermitian matrices with four texture zeros.

## III Analyses of Matrices with Four Texture Zeros

The study of viable 4-texture-zero pairs is straightforward by following the procedure outlined above. For convenience in our analysis and presentation, we categorize the mass matrix pairs $`(M_U,M_D)`$ according to whether or not the $`(1,1)`$ matrix element is zero:
1) $`M_U^{11}=0,M_D^{11}=0`$,
2) $`M_U^{11}0,M_D^{11}=0`$,
3) $`M_U^{11}=0,M_D^{11}0`$,
4) $`M_U^{11}0,M_D^{11}0`$.
The results are presented in Tables I-IV, which contain all the viable 4-texture-zero pairs. It is seen from the Tables that almost all of the viable textures have $`M_D^{11}=0`$.

The first category involves eight types of hermitian matrices, which are listed in Table I along with their corresponding triangular forms. These hermitian matrices contain 1, 2 and 3 texture zeros for $`M_U`$ and/or $`M_D`$. Their corresponding triangular matrices can be simply obtained from Eqs. (10) and (14). They are then paired up to form possible patterns for the quark mass matrices. Note that $`M_2,M_3,M_6,M_7`$ and $`M_8`$ are the ones appearing in the analysis of five texture-zero matrices . One notable feature about the pairing is that $`M_U`$ and $`M_D`$ allow different texture zeros, as can be seen from Eqs. (10) and (14). For example, whereas $`M_D^{22}0`$ (unless we give up naturalness), $`M_U^{22}=0`$ is allowed due to the much larger mass hierarchy in the up-quark sector. As a second example, $`M_U^{11}=M_U^{22}=0`$ is allowed but the same relation is not valid for $`M_D`$ again because of the different mass hierarchies. As another tip for the pairing, a quick check on the CKM matrix element can be useful in the screening of possible candidates. For instance, the $`(M_4,M_4)`$ pair is out because both their (2,3) elements are zero, which in turns gives $`V_{cb}0`$. The same can not be said about hermitian pairs with vanishing (1,3) elements. Here, $`V_{ub}`$ and $`V_{td}`$ may or may not vanish since they are of higher order in $`\lambda `$, and may be induced from a combined $`R_{12}`$ and $`R_{23}`$ rotations. In this latter case, a further investigation is required.

The viability of each $`(M_U,M_D)`$ pair can be assessed by converting them, through common LH rotations if necessary, into one of the ten triangular pairs in the m.p.b. as listed in Table I of Ref. . Testable relations can then be obtained. To illustrate the method, we first study the $`(M_8,M_1)`$ pair in detail. The triangular forms for $`M_8`$ and $`M_1`$ are given by:

$$M_8M_U^t\left(\begin{array}{ccc}a_U\lambda ^8& b_U\lambda ^6& 0\\ 0& d_U\lambda ^4& e_U\lambda ^2\\ 0& 0& 1\end{array}\right),$$
(15)
$$M_1M_D^t\left(\begin{array}{ccc}a_D\lambda ^4& b_D\lambda ^3& c_D\lambda ^3\\ 0& d_D\lambda ^2& e_D\lambda ^2\\ 0& 0& 1\end{array}\right).$$
(16)
Here, $`d_U=|e_U|^2`$. This ($`M_U^t,M_D^t`$) pair has eleven nonzero elements, and two more zeros are needed to put it into one of the ten triangular pairs in the m.p.b.. There are several ways to achieve this. For example, a common LH 2-3 rotation with $`\theta _{23}e_U\lambda ^2`$ will set the $`(2,3)`$ element of $`M_U`$ to zero, and a common LH 1-2 rotation with $`\theta _{12}\frac{b_D}{d_D}\lambda `$ will generate a second zero at the $`(1,2)`$ position of $`M_D`$. In this way, we arrive at

$$M_U^t^{}=R_{12}(\frac{b_D}{d_D}\lambda )R_{23}(e_U\lambda ^2)M_U^t\left(\begin{array}{ccc}a_U\lambda ^8& b_U\lambda ^6\frac{b_D}{d_D}d_U\lambda ^5& 0\\ 0& d_U\lambda ^4& 0\\ 0& 0& 1\end{array}\right),$$
(17)
$$M_D^t^{}=R_{12}(\frac{b_D}{d_D}\lambda )R_{23}(e_U\lambda ^2)M_D^t\left(\begin{array}{ccc}a_D\lambda ^4& 0& c_D\lambda ^3\frac{b_D}{d_D}(e_De_U)\lambda ^3\\ 0& d_D\lambda ^2& (e_De_U)\lambda ^2\\ 0& 0& 1\end{array}\right).$$
(18)
Identifying $`(M_U^t^{},M_D^t^{})`$ with the 5th triangular pattern in Table I of Ref.  gives,

$$\lambda ^2b_U/d_U\lambda b_D/d_D=V_{us}/V_{cs}.$$
(19)
As $`M_U^{11}=M_D^{11}=0`$ in the hermitian form, we have two corresponding relations in terms of triangular parameters: $`a_U=\frac{|b_U|^2}{d_U}`$ and $`a_D=\frac{|b_D|^2}{d_D}`$. With these two relations, Eq. (19) can be written in the well known form,

$$\left|\frac{V_{us}}{V_{cs}}\right|\left|\sqrt{\frac{m_d}{m_s}}e^{i\delta }\sqrt{\frac{m_u}{m_c}}\right|,$$
(20)
where $`\delta \text{arg}[b_Ud_D/d_Ub_D]`$. From Eqs. (17-18) and Table I of Ref. , we note that the standard model $`CP`$-violation depends on additional phases besides $`\delta `$, thus leaving $`\delta `$ a free parameter. In this sense, Eq. (20) places a constraint on the mass matrices in fixing the phase $`\delta `$, but it is not a prediction in terms of physical quantities alone. We conclude that the $`(M_8,M_1)`$ pair is a viable texture since Eq. (20) can be made valid with a properly chosen phase $`\delta `$.

To assess the viability of any given texture, we will need to know about the quark masses and CKM elements. We use for the quark masses at $`m_Z`$ values taken from Ref. . The CKM matrix elements are taken from , except for $`V_{ub}`$ and $`V_{td}`$, for which we use a recent update :

$$\left|\frac{V_{ub}}{V_{cb}}\right|_{exp}=0.093\pm 0.014,$$
(21)
$$0.15<\left|\frac{V_{td}}{V_{ts}}\right|_{exp}<0.24.$$
(22)
Whereas Eq. (21) comes from an average of the LEP and SLD measurements, Eq. (22) comes from a standard model fit to the electroweak data. Note that Eq. (22) implies $`|V_{td}|<0.01`$.

After a straightforward analysis of all possible texture pairs, we arrive at the nine viable pairs of hermitian matrices for the first category. These are listed in Table II, together with their testable relations. Note that relations R1 and R2 allow two different solutions depending on the sign of $`d_U`$, and that some textures share the same quark mass-mixing relation.

Also listed in Table II is the particular pair (pattern 10) which is parallel in structure: $`(M_3,M_3)`$. This texture pair has been the center of focus in recent studies of hermitian 4-texture-zero quark mass matrices. To analyze the predictions of this pair, we first write down the triangular form,

$$M_3M_U^t\left(\begin{array}{ccc}a_U\lambda ^8& b_U\lambda ^6& 0\\ 0& d_U\lambda ^4& e_U\lambda ^2\\ 0& 0& 1\end{array}\right),$$
(23)
$$M_3M_D^t\left(\begin{array}{ccc}a_D\lambda ^4& b_D\lambda ^3& 0\\ 0& d_D\lambda ^2& e_D\lambda ^2\\ 0& 0& 1\end{array}\right).$$
(24)
We need one more zero to put the triangular pair $`(M_U^t,M_D^t)`$ in the m.p.b., and this can be attained by a common LH 2-3 rotation of $`\theta _{23}e_U\lambda ^2`$ ($`\theta _{23}e_D\lambda ^2`$) to set the $`(2,3)`$ element of $`M_U^t`$ ($`M_D^t`$) to zero. For example, we can have

$$M_U^t^{}=R_{23}(e_U\lambda ^2)M_U^t\left(\begin{array}{ccc}a_U\lambda ^8& b_U\lambda ^6& 0\\ 0& d_U\lambda ^4& 0\\ 0& 0& 1\end{array}\right),$$
(25)
$$M_D^t^{}=R_{23}(e_U\lambda ^2)M_D^t\left(\begin{array}{ccc}a_D\lambda ^4& b_D\lambda ^3& 0\\ 0& d_D\lambda ^2& (e_De_U)\lambda ^2\\ 0& 0& 1\end{array}\right).$$
(26)
Note that $`(M_U^t^{},M_D^t^{})`$ corresponds to the 2nd of the ten triangular pairs in the m.p.b. as listed in Table I of Ref . A simple comparison between them gives the following relations for the hermitian $`(M_3,M_3)`$ texture:

$$\frac{V_{ub}}{V_{cb}}\frac{b_U}{d_U}\lambda ^2,\frac{V_{td}^{}}{V_{cb}V_{cs}}\frac{b_D}{d_D}\lambda ,$$
(27)
and

$$\text{arg}\left[\frac{b_Ud_D}{d_Ub_D}\right]\alpha \text{arg}\left[\frac{V_{td}V_{tb}^{}}{V_{ud}V_{ub}^{}}\right],$$
(28)
where $`\alpha `$ is the $`\alpha `$-angle of unitarity triangle. Now the two zeros at the $`(1,1)`$ position of $`(M_3,M_3)`$ imply two relations among the triangular parameters: $`|b_U^2|=a_Ud_U`$ and $`|b_D^2|=a_Dd_D`$. This allows us to rewrite Eq. (27) in the form,

$$\left|\frac{V_{ub}}{V_{cb}}\right|\sqrt{\frac{m_u}{m_c}},$$
(29)
$$\left|\frac{V_{td}}{V_{ts}}\right|\sqrt{\frac{m_d}{m_s}}.$$
(30)
Using the unitarity of the CKM matrix, the two relations in Eq. (27) can be combined to give

$`{\displaystyle \frac{V_{us}}{V_{cs}}}`$ $`=`$ $`{\displaystyle \frac{V_{ub}}{V_{cb}}}+{\displaystyle \frac{V_{td}^{}}{V_{cb}V_{cs}}}{\displaystyle \frac{b_D}{d_D}}\lambda {\displaystyle \frac{b_U}{d_U}}\lambda ^2,`$ (31)
where we have assumed $`\mathrm{det}V^{\text{CKM}}=1`$ without loss of generality. Using $`|b_U^2|=a_Ud_U`$, $`|b_D^2|=a_Dd_D`$, and Eq. (28), Eq. (31) can be cast in the form,

$$\left|\frac{V_{us}}{V_{cs}}\right|\left|\sqrt{\frac{m_d}{m_s}}e^{i\alpha }\sqrt{\frac{m_u}{m_c}}\right|.$$
(32)
Note that Eq. (32) is similar but different from Eq. (20): while the phase in Eq. (32) is in principle fixed by the $`CP`$-violation parameter $`ϵ_K`$ to be $`\alpha \pi /2`$ , the phase $`\delta `$ in Eq. (20) can be varied freely as the standard model $`CP`$ violation depends on other phases as well. The relations in Eqs.(29), (30), and (32) were also obtained in Ref. .

Whereas Eqs. (30) and (32) are allowed by data, the prediction of Eq. (29), $`V_{ub}/V_{cb}=0.059\pm 0.006`$, is disfavored by the recent measurement (see Eq. (21) ). The same conclusion has been reached for three of the five hermitian pairs in our five-texture-zero analysis. It is interesting and surprising that going to four-texture-zero does not help with this problem of low $`V_{ub}/V_{cb}`$. Furthermore, even varying $`m_u`$ and $`m_c`$ in a reasonable range will not be able to accommodate a value of $`V_{ub}/V_{cb}=0.08`$ . In this regard, we are in disagreement with the analysis of Ref. , where much larger values ($`0.0830.099`$) for $`V_{ub}/V_{cb}`$ were obtained. We note also that texture zeros are not invariant under change of basis. Our result for the pair ($`M_3,M_3`$) is valid only in the hierarchical basis. Otherwise, larger values of $`V_{ub}/V_{cb}`$ are allowed, as in Eqs.(23-26) of Ref. . For these reasons, we exclude the $`(M_3,M_3)`$ pair from the “Yes” column in Table II.

It is worthwhile to compare in detail the results of our two examples. Both the pairs $`(M_8,M_1)`$ and $`(M_3,M_3)`$ have four texture zeros. But the former yields one relation that is not physically predictive in nature (Eq.(20)), while the later gives rise to two independent predictions of Eqs.(29-30) (note that Eq. (32) is not independent).

In general, for a pair of hermitian matrices with four texture zeros, there are eight real parameters plus one or more phases, so we expect at most one prediction. This is the case for the pair $`(M_8,M_1)`$, where we have more than one unremovable phases and no physical prediction. The situation is different with $`(M_3,M_3)`$. From Eqs.(25-26), using $`a_U=|b_U|^2/d_U`$ and $`a_D=|b_D|^2/d_D`$, the independent parameters are $`|b_U|`$, $`|b_D|`$, $`|d_U|`$, $`|d_D|`$, $`|e_Ue_D|`$, plus two overall mass scales and one physical phase. Thus, triangularization reveals that the pair $`(M_3,M_3)`$ actually contains one less parameter than expected, because only the combination $`(e_Ue_D)`$ enters. Put another way, we could set $`e_U=0`$ (or $`e_D=0`$) without any effect. This means that the physical contents of $`(M_3,M_3)`$ are the same as the five texture zero pair $`(M_7,M_3)`$ or $`(M_3,M_7)`$. In addition, we can see that the five-zero pair $`(M_8,M_3)`$ also reduces to the same triangular form by a LH 2-3 rotation. Thus, we conclude that all four pairs, $`(M_3,M_3)`$, $`(M_7,M_3)`$, $`(M_3,M_7)`$, and $`(M_8,M_3)`$, are physically equivalent. The analytical predictions of the last three pairs were studied in Ref. .

Similar analyses of the hermitian pairs with $`M_U^{11}0`$ and/or $`M_D^{11}0`$ follow directly. The results are presented in Tables III and IV. In Table III, we list the nine hierarchical hermitian matrices with $`(1,1)0`$ that are used in the construction of viable hermitian pairs. The corresponding triangular matrices are also given in the table. The viable hermitian 4-texture-zero pairs constructed from Tables I and III are listed in Table IV, together with their predictions for quark mixing. Whereas there is no viable pair with both $`M_U^{11}0`$ and $`M_D^{11}0`$, nine pairs with $`M_U^{11}0`$ and $`M_D^{11}=0`$ are found to be compatible with data. Of the hermitian matrices with $`M_U^{11}=0`$ and $`M_D^{11}0`$, one pair is viable, and the 2nd pair, $`(M_2,M_{12})`$, leads to the same prediction for $`V_{ub}`$ as that of the 5th five-texture-zero pair (see Eq.(15) in Ref. ). This pair allows two different relations depending on the sign of $`d_U`$. While the plus-sign relation is marginally acceptable, the minus-sign relation is disfavored by data.

Before we leave this section, we would like to point out that some of the disfavored hermitian pairs can give rise to nontrivial predictions. For example, the four pairs, $`(M_1,M_{17})`$ and $`(M_i,M_{13})`$ ($`i=2,4,5`$), all lead to the relation:

$`\left|V_{td}\right|`$ $``$ $`\sqrt{V_{us}^2m_c/m_t+m_u/m_t}=|V_{us}|\sqrt{m_c/m_t}\times (1+𝒪(\lambda ^2))>1\%.`$ (33)
This prediction contradicts the experimental limit of $`|V_{td}|<1\%`$, as can be derived from $`B\overline{B}`$ mixing (see also Eq. (22)). As another example, the hermitian pair $`(M_2,M_{11})`$ entails the relation:

$`|V_{cb}|`$ $``$ $`\sqrt{{\displaystyle \frac{m_c}{m_t}}}\left(1+{\displaystyle \frac{m_u}{m_cV_{us}^2}}\right)^{\frac{1}{2}}=\sqrt{{\displaystyle \frac{m_c}{m_t}}}\times (1+𝒪(\lambda ^2))0.06,`$ (34)
which is too large for this texture to be viable.

A notable feature about Tables I-IV is that some texture pairs lead to the same relation, like patterns 1, 2, 5, and 6. It can be shown that these texture pairs are related by weak basis transformation, and thus are physically equivalent.

## IV Conclusions

Because of the many redundancies in the quark mass matrices, a lot of work has been done in search for more restrictive patterns of matrices which are still compatible with experiments. In the literature, such considerations have focused on hermitian mass matrices with a certain number of texture zeros. It was found that the maximum allowed number of texture zeros is five. Amongst them one unique pattern has been identified which is most favorable with data. It is worthwhile to understand the situation for pairs of matrices with fewer texture zeros.

In this paper, we have investigated systematically hermitian hierarchical quark mass matrices with four texture zeros. By transforming the hermitian matrices into the triangular form in the minimal parameter basis, and using the fact that the latter contains only physical masses and CKM matrix elements, we can quickly rule out many of the possible pairs of mass matrices. For the remaining candidates, we showed that each pattern of texture zeros implies certain relations between quark masses and mixing. These relations can be used to determine the viability of each pair of mass matrices by comparing them with available data. In this way we identified nineteen pairs of mass matrices which are compatible with current data. One pair (pattern 21) is marginal, with the final verdict depending on more accurate experimental numbers. Of the nineteen viable textures, none has a parallel structure between the up and down mass matrices. In particular, the popular parallel-structure pair that has been a focus of much attention, pattern 10, leads to a low value of $`V_{ub}/V_{cb}`$, and is not favorable with present data. We also found that, by a proper LH rotation, pattern 10 can be shown to have the same physical contents as three (which are themselves equivalent) of the five pairs studied earlier . The asymmetry between up and down quark mass matrix textures could serve as a useful guideline in search for realistic models of quark-lepton masses.

###### Acknowledgements.

We would like to thank Sadek Mansour for helpful discussions. G.W. would like to thank the National Center for Theoretical Sciences of Taiwan for its hospitality during the final stage of this work. S.C. is supported by the Purdue Research Foundation. T.K. and G.W. are supported by DOE, Grant No. DE-FG02-91ER40681 and No. DE-FG03-96ER40969, respectively.