text:
# An SU(5) SUSY Model with 𝑅-Parity Violation and Radiatively Induced Neutrino Masses

## I Introduction

The origin of the neutrino mass generation is still a mysterious problem in the unified understanding of the quarks and leptons. The Zee model  is one of several promising models, because it has only 3 free parameters and it can naturally lead to a large neutrino mixing , especially, to a bimaximal mixing . However, the original Zee model is not on the framework of a grand unification theory (GUT). The most attractive idea  to embed the Zee model into GUTs is to identify the Zee scalar $`h^+`$ as the slepton $`\stackrel{~}{e}_R`$ in an $`R`$-parity-violating supersymmetric (SUSY) model. However, usually, it is accepted that SUSY models with $`R`$-parity violation are incompatible with a GUT scenario, because the $`R`$-parity-violating interactions induce proton decay .

In the present paper, in order to suppress this kind of proton decay, a discrete symmetry Z2 is introduced. The essential idea is as follows: At the unification scale $`\mu =M_X`$, we have two types of superfields $`\mathrm{\Psi }_{L(\pm )}=\overline{5}_{L(\pm )}+10_{L(\pm )}`$, which are transformed with $`\pm 1`$ under the discrete symmetry Z2 (we will call it “Z2-parity” hereafter). We consider that the particles $`\mathrm{\Psi }_{L(+)}`$ are decoupled after the SU(5) symmetry is broken (but Z2 is still unbroken), so that our quarks and leptons $`\overline{5}_L+10_L`$ belong to $`\mathrm{\Psi }_{L()}`$. The $`R`$-parity violating terms are given by the combinations $`\overline{5}_{(+)}\overline{5}_{(+)}10_{(+)}`$, $`\overline{5}_{()}\overline{5}_{()}10_{(+)}`$ and $`\overline{5}_{(+)}\overline{5}_{()}10_{()}`$, so that they basically do not contribute to the quarks and leptons with Z2-parity $`z=1`$, because of the Z2 symmetry. However, we assume that mixings of the members of $`\mathrm{\Psi }_{L(+)}`$ with those of $`\mathrm{\Psi }_{L()}`$ are caused in part after SU(5) is broken. As a result, the $`R`$-parity-violating interactions are in part allowed, so that the neutrino masses are radiatively generated, while the proton decay due to the $`R`$-parity-violating terms is still forbidden because the term $`d_R^cd_R^cu_R^c`$ is still exactly forbidden below $`\mu =M_X`$ in the present scheme. The details will be discussed in the next section.

The purpose of the present paper is to investigate the possible forms of the radiatively induced neutrino mass matrix under the Z2 symmetry. In Sec. III, we will give them, including a numerical study. In Sec. IV, we will give a comment on the Higgs scalars in the present scheme. Finally, Sec. V is devoted to our conclusion.

## II Z2 symmetry and the proton decay

We identify the Zee scalar $`h^+`$ as the slepton $`\stackrel{~}{e}_R^+`$, which is a member of SU(5) 10-plet sfermions $`\stackrel{~}{\psi }_{10}`$. Then, the Zee interactions correspond to the following $`R`$-parity-violating interactions

$`\lambda _{ijk}(\overline{\psi }_{\overline{5}}^c)_i^A(\psi _{\overline{5}})_j^B(\stackrel{~}{\psi }_{10})_{kAB}`$ (1)
$`={\displaystyle \frac{1}{\sqrt{2}}}\lambda _{ijk}\{\epsilon _{\alpha \beta \gamma }(\overline{d}_R)_i^\alpha (d_R^c)_j^\beta (\stackrel{~}{u}_R^{})_k^\gamma `$ (2)
$`[(\overline{e}_L^c)_i(\nu _L)_j(\overline{\nu }_L^c)_i(e_L)_j](\stackrel{~}{e}_R^{})_k`$ (3)
$`[(\overline{e}_L^c)_i(d_R^c)_j^\alpha (\overline{d}_R)_i^\alpha (e_L)_j](\stackrel{~}{u}_L)_{k\alpha }`$ (4)
$`+[(\overline{\nu }_L^c)_i(d_R^c)_j^\alpha (\overline{d}_R)_i^\alpha (\nu _L)_j](\stackrel{~}{d}_L)_{k\alpha }\},`$ (5)
where $`\psi ^cC\overline{\psi }^T`$ and the indices $`(i,j,\mathrm{})`$, $`(A,B,\mathrm{})`$ and $`(\alpha ,\beta ,\mathrm{})`$ are family-, SU(5)GUT\- and SU(3)colour-indices, respectively. The coefficients $`\lambda _{ijk}`$ are antisymmetric in $`i`$ and $`j`$. On the other hand, in SUSY GUT models, if the interactions (2.1) exist, the following $`R`$-parity-violating interactions will also exist:

$`\lambda _{ijk}(\overline{\psi }_{\overline{5}}^c)_i^A(\psi _{10})_{kAB}(\stackrel{~}{\psi }_{\overline{5}})_j^B,`$ (6)
$`={\displaystyle \frac{1}{\sqrt{2}}}\lambda _{ijk}\{\epsilon _{\alpha \beta \gamma }(\overline{d}_R)_i^\alpha (\stackrel{~}{d}_R^{})_j^\beta (u_R^c)_k^\gamma `$ (7)
$`[(\overline{e}_L^c)_i(\stackrel{~}{\nu }_L)_j(\overline{\nu }_L^c)_i(\stackrel{~}{e}_L)_j](e_R^c)_k`$ (8)
$`[(\overline{e}_L^c)_i(\stackrel{~}{d}_R^{})_j^\alpha (\overline{d}_R)_i^\alpha (\stackrel{~}{e}_L)_j](u_L)_{k\alpha }`$ (9)
$`+[(\overline{\nu }_L^c)_i(\stackrel{~}{d}_R^{})_j^\alpha (\overline{d}_R)_i^\alpha (\stackrel{~}{\nu }_L)_j](d_L)_{k\alpha }\},`$ (10)
which contribute to the proton decay through the intermediate state $`\stackrel{~}{d}_R`$. Also, the term $`(\overline{d}_R)_i^\alpha (d_R^c)_j^\beta (\stackrel{~}{u}_R^{})_k^\gamma `$ in the interactions (2.1) can contribute to the nucleon decay through the intermediate state $`\stackrel{~}{u}_R`$. The upper limits of the coupling constants $`\lambda _{ijk}`$ from proton decay experiments have been investigated by Smirnov and Vissani , and the values must be highly suppressed.

In order to forbid the contribution of the interactions (2.1) and (2.2) to the proton decay, we must consider that in the $`R`$-parity-violating interactions $`\overline{5}\times \overline{5}\times 10`$, the term $`d_R^cd_R^cu_R^c`$ is exactly forbidden, while the terms $`\nu _Le_Le_R^c`$ and/or $`\nu _Ld_R^cd_L`$ are in part allowed.

For such purpose, we introduce a discrete symmetry Z2, which exactly holds at every energy scale. At the unification scale $`\mu =M_X`$, we have two types of superfields $`\mathrm{\Psi }_{L(\pm )}=\overline{5}_{L(\pm )}+10_{L(\pm )}`$, which are transformed as $`\mathrm{\Psi }_{L(\pm )}\pm \mathrm{\Psi }_{L(\pm )}`$ under the discrete symmetry Z2. We consider that the particles $`\mathrm{\Psi }_{L(+)}`$ are basically decoupled after the SU(5) symmetry is broken, so that our quarks and leptons (and their SUSY partners) $`\overline{5}_L+10_L`$ are regarded as $`\mathrm{\Psi }_{L()}=\overline{5}_{L()}+10_{L()}`$. The $`R`$-parity-violating terms for quarks and leptons (and their SUSY partners) are basically forbidden under the symmetry Z2 below $`\mu =M_X`$, because the terms are composed of $`\overline{5}_{L()}\overline{5}_{L()}10_{L()}`$.

However, if we assume that mixings between the members of $`\mathrm{\Psi }_{L(+)}`$ and those of $`\mathrm{\Psi }_{L()}`$ in part take place after SU(5) is broken, $`R`$-parity-violating interactions $`\mathrm{\Psi }_{L(+)}\mathrm{\Psi }_{L()}\mathrm{\Psi }_{()}`$ become available at the low energy $`\mu =m_Z`$, too. For example, we assume a mixing

$$(2,1)_{Li}=(2,1)_{L()i}\mathrm{cos}\theta _i^A+(2,1)_{L(+)i}\mathrm{sin}\theta _i^A,$$
(11)
between the $`(2,1)`$ components of SU(2)$`\times `$SU(3) for the $`i`$-th family. (Hereafter, we will refer to the mixing (2.3) as a mixing of type Ai.) Then, the $`R`$-parity-violating interactions

$`\mathrm{sin}\theta _i^A\lambda _{ijk}\nu _{Li}d_{Rj}^cd_{Lk},\mathrm{sin}\theta _i^A\lambda _{ijk}e_{Li}d_{Rj}^cu_{Lk},`$ (12)
$`\mathrm{sin}\theta _i^A\mathrm{cos}\theta _j^A\lambda _{ijk}\nu _{Li}e_{Lj}e_{Rk}^c,`$ (13)
become available from the interactions

$$\lambda _{ijk}\overline{5}_{(+)i}\overline{5}_{()j}10_{()k},$$
(14)
above the unification scale $`\mu =M_X`$. Also, we can consider a mixing

$$(2,3)_{Lk}=(2,3)_{L()k}\mathrm{cos}\theta _k^B+(2,3)_{L(+)k}\mathrm{sin}\theta _k^B,$$
(15)
between the $`(2,3)`$ components of SU(2)$`\times `$SU(3) for the $`k`$-th family. (Hereafter, we will refer to the mixing (2.6) as a Bk-type mixing.) Then, the $`R`$-parity-violating interactions

$$\mathrm{sin}\theta _k^B\lambda _{ijk}^{}\nu _{Li}d_{Rj}^cd_{Lk}\mathrm{and}\mathrm{sin}\theta _k^B\lambda _{ijk}^{}e_{Li}d_{Rj}^cu_{Lk}$$
(16)
become available from the interactions

$$\lambda _{ijk}^{}\overline{5}_{()i}\overline{5}_{()j}10_{(+)k}.$$
(17)
On the other hand, note that the interaction

$$d_R^cd_R^cu_R^c$$
(18)
is exactly forbidden, independently of whether the mixings (2.3) and (2.6) occur or not, because those interactions have the Z2 parity $`z=1`$. Therefore, the proton decay due to the $`R`$-parity-violating terms is exactly forbidden because of the absence of the term $`d_R^cd_R^cu_R^c`$. On the other hand, the neutrino masses are radiatively generated through the interactions $`\nu _Ld_R^cd_L`$ and $`\nu _Le_Le_R^c`$ with $`z=+1`$. The possible forms of the radiative neutrino mass matrix will be discussed in the next section.

At present, we do not know a reasonable mechanism not only for such a mixing, but also for the decoupling of $`\mathrm{\Psi }_{L(+)}`$. In order to make $`\mathrm{\Psi }_{L(+)}=\overline{5}_{L(+)}+10_{L(+)}`$ heavy, the SU(2)L symmetry must be broken, but, of course, we cannot consider a scenario in which SU(2) is broken just after SU(5) is broken. In the present paper, we give only a phenomenological selection rule: if the superfield $`\mathrm{\Psi }`$ can make a five-body SU(5) singlet operator $`\mathrm{\Psi }\mathrm{\Psi }\mathrm{\Psi }\mathrm{\Psi }\mathrm{\Psi }`$ with the Z2 parity $`z=+1`$, then the superfield $`\mathrm{\Psi }`$ can be decoupled below $`\mu =M_X`$. Obviously, according to this selection rule, the superfield $`\overline{5}_{L(+)}`$ can be decoupled below $`\mu =M_X`$. Similarly, the superfield $`10_{L(+)}`$ is decoupled below $`\mu =M_X`$. However, note that those operators in the SU(5) singlets are symbolically expressed in terms of SU(2)$`\times `$SU(3) components as follows:

$`\left(\overline{5}_{L(+)}\right)^5`$ $`=`$ $`[(2,1)_{L(+)}]^2\times [(1,3)_{L(+)}]^3,`$ (19)
$`\left(10_{L(+)}\right)^5`$ $`=`$ $`(1,1)_{L(+)}\times [(1,\overline{3})_{L(+)}]^2\times [(2,3)_{L(+)}]^2`$ (21)
$`+(1,\overline{3})_{L(+)}\times [(2,3)_{L(+)}]^4,`$
and that even if the interchanges $`(2,1)_{L(+)i}(2,1)_{L()i}`$, and/or $`(2,3)_{L(+)k}(2,3)_{L()k}`$, are caused, the composite operators

$$\left(\overline{5}_L\right)^5=[(2,1)_{L()}]^2\times [(1,3)_{L(+)}]^3,$$
(22)
$`\left(10_L\right)^5`$ $`=`$ $`(1,1)_{L(+)}\times [(1,\overline{3})_{L(+)}]^2\times [(2,3)_{L()}]^2`$ (24)
$`+(1,\overline{3})_{L(+)}\times [(2,3)_{L()}]^4`$
still have $`z=+1`$. Such interchanges are possible only for the components $`(2,1)_L`$ and $`(2,3)_L`$. As a result, only the combination

$`\overline{5}_L+10_L=[(2,1)_{L(+)}+(1,\overline{3})_{L()}]`$ (25)
$`+[(1,1)_{L()}+(2,3)_{L()}+(1,\overline{3})_{L()}],`$ (26)
for the $`i`$-th family and/or

$`\overline{5}_L+10_L=[(2,1)_{L()}+(1,\overline{3})_{L()}]`$ (27)
$`+[(1,1)_{L()}+(2,3)_{L(+)}+(1,\overline{3})_{L()}],`$ (28)
for the $`k`$-th family survive below $`\mu =M_X`$ as the quarks and leptons (and their SUSY partners).

Of course, the above selection rule cannot be justified within the framework of the minimal SUSY standard model. At present, this is only an ansatz to select which components of SU(2)$`\times `$SU(3) can be interchanged.

## III Radiative neutrino masses

In a SUSY GUT scenario, there are many origins of the neutrino mass generations. For example, the sneutrinos $`\stackrel{~}{\nu }_{iL}`$ can have vacuum expectation values (VEVs), and the neutrinos $`\nu _{Li}`$ acquire their masses thereby (for example, see Ref.). Although we cannot rule out a possibility that the observed neutrino masses can be understood from such compound origins, we do not take such a point of view in the present paper, because the observed neutrino masses and mixings appear to be rather simple and characteristic. We simply assume that the radiative masses are only dominated even if there are other origins of the neutrino mass generations.

In the present scenario, the origins of the radiatively induced neutrino masses are two: one is induced by the $`R`$-parity-violating interactions $`\nu _Ld_R^c\stackrel{~}{d}_L`$ and $`\nu _L\stackrel{~}{d}_R^cd_L`$; the other one is induced by $`\nu _Le_L\stackrel{~}{e}_R^c`$ and $`\nu _L\stackrel{~}{e}_Le_R^c`$. Note that there is no Zee-type diagrams due to $`H_d^+`$$`\stackrel{~}{e}_R^+`$ mixing in this scheme.

First, we discuss the down-quark loop contributions. For simplicity, we assume that the masses $`\stackrel{~}{M}_{Li}`$ and $`\stackrel{~}{M}_{Ri}`$ of the squarks $`\stackrel{~}{d}_{Li}`$ and $`\stackrel{~}{d}_{Ri}`$ are approximately constant, independently of the flavours, although we consider the flavour-dependent structure for the mass terms $`\stackrel{~}{d}_L^{}\stackrel{~}{M}_d^2\stackrel{~}{d}_R`$. Then, the radiatively induced neutrino mass matrix due to the A-type mixing is given by

$$(M_\nu )_{ij}=m_0\lambda _{ikm}\lambda _{jln}(M_d^{})_{kn}(\stackrel{~}{M}_d^2)_{lm}+(ij),$$
(29)
where sine-factors have been drooped for simplicity, $`M_\nu `$ is defined by $`\overline{\nu }_LM_\nu \nu _L^c`$, and the coupling constants $`\lambda _{ijk}`$ are redefined by

$$\lambda _{ijk}[\overline{\nu }_{Li}d_{Rj}\stackrel{~}{d}_{Lk}^{}+\overline{\nu }_{Li}\stackrel{~}{d}_{Rj}d_{Lk}^c(\nu _Ld_R^c)].$$
(30)
Here, we have changed the definition of $`\lambda _{ijk}`$ from that in (2.1) as $`\lambda _{ijk}\lambda _{ijk}^{}`$ for the convenience of the expression of $`M_\nu `$ defined by $`\overline{\nu }_LM_\nu \nu _L^c`$. In the present paper, the unitary matrix $`U_\nu `$ used to diagonalize the Majorana mass matrix $`M_\nu `$ is defined as $`U_\nu ^{}M_\nu U_\nu ^{}=D_\nu `$. Then, the so-called Maki–Nakagawa–Sakata–Pontecorvo  matrix (we will simply call it the “lepton mixing matrix”) $`UU_{MNSP}`$ is given by $`U=U_L^eU_\nu `$. Usually, it is considered that the matrix form of $`\stackrel{~}{M}_d^2`$ is proportional to the form $`M_d`$. Then, the neutrino mass matrix (3.1) becomes, in a more concise form:

$$(M_\nu )_{ij}=m_0\lambda _{ikm}\lambda _{jln}(M_d^{})_{kn}(M_d^{})_{lm},$$
(31)
where we have redefined the common factor $`m_0`$ from that in (3.1). Of course, exactly speaking, for the Ai mixings with $`\mathrm{sin}\theta _i^A`$-factors, we should read the expression (3.3) as

$$(M_\nu )_{ij}=m_0s_i^As_j^A\lambda _{ikm}\lambda _{jln}(M_d^{})_{kn}(M_d^{})_{lm},$$
(32)
where $`s_i^A=\mathrm{sin}\theta _i^A`$ defined in (2.3). When we consider the Bk mixings, we read the expression (3.3) as

$$(M_\nu )_{ij}=m_0s_m^Bs_n^B\lambda _{ikm}^{}\lambda _{jln}^{}(M_d^{})_{kn}(M_d^{})_{lm},$$
(33)
where $`s_k^B=\mathrm{sin}\theta _k^B`$ defined in (2.6). For mixed-type mixings of Ai and Bk, we read (3.3) as

$`(M_\nu )_{ij}`$ $`=`$ $`m_0(s_i^Ac_m^B\lambda _{ikm}+c_i^As_m^B\lambda _{ikm}^{})`$ (35)
$`(s_j^Ac_n^B\lambda _{jln}+c_j^As_n^B\lambda _{jln}^{})(M_d^{})_{kn}(M_d^{})_{lm},`$
where $`c_i^A=\mathrm{cos}\theta _i^A`$ and $`c_k^B=\mathrm{cos}\theta _k^B`$.

The contributions from the charged lepton loops are essentially the same as (3.3), except for the absence of the B-type mixing and the replacement $`M_dM_e^T`$. For simplicity, we will continue the investigation for the case of the down-quark loop contributions.

For the phenomenological study of the mass matrix (3.3), it is convenient to take the basis on which the down-quark mass matrix $`M_d`$ is diagonal:

$$U_L^dM_dU_R^d=D_d\mathrm{diag}(m_1^d,m_2^d,m_3^d).$$
(36)
We consider that, on the basis with $`M_d=D_d`$, the charged lepton mass matrix $`M_e`$ is also approximately diagonal, $`U_L^eM_eU_RD_e=\mathrm{diag}(m_1^e,m_2^e,m_3^e)`$, so that the unitary matrix $`U_\nu `$ approximately gives the lepton mixing matrix $`U=U_L^eU_\nu `$. Then, we can express (3.3) as

$`(M_\nu )_{11}`$ $`=`$ $`(m_3^d)^2(\lambda _{133})^2+(m_2^d)^2(\lambda _{122})^2+2m_3^dm_2^d\lambda _{123}\lambda _{132},`$ (37)
$`(M_\nu )_{22}`$ $`=`$ $`(m_3^d)^2(\lambda _{233})^2+(m_1^d)^2(\lambda _{211})^2+2m_3^dm_1^d\lambda _{213}\lambda _{231},`$ (38)
$`(M_\nu )_{33}`$ $`=`$ $`(m_2^d)^2(\lambda _{322})^2+(m_1^d)^2(\lambda _{311})^2+2m_2^dm_1^d\lambda _{312}\lambda _{321},`$ (39)
$`(M_\nu )_{12}`$ $`=`$ $`(m_3^d)^2\lambda _{133}\lambda _{233}+m_3^dm_2^d\lambda _{123}\lambda _{232}`$ (41)
$`+m_3^dm_1^d\lambda _{131}\lambda _{213}+m_2^dm_1^d\lambda _{121}\lambda _{212},`$
$`(M_\nu )_{13}`$ $`=`$ $`(m_2^d)^2\lambda _{122}\lambda _{322}+m_3^dm_2^d\lambda _{132}\lambda _{323}`$ (43)
$`+m_3^dm_1^d\lambda _{131}\lambda _{313}+m_2^dm_1^d\lambda _{121}\lambda _{312},`$
$`(M_\nu )_{23}`$ $`=`$ $`(m_1^d)^2\lambda _{211}\lambda _{311}+m_3^dm_2^d\lambda _{232}\lambda _{323}`$ (45)
$`+m_3^dm_1^d\lambda _{231}\lambda _{313}+m_2^dm_1^d\lambda _{212}\lambda _{321},`$
where, for simplicity, we have dropped the common factor $`m_0`$. In order to give the best-fit values for the observed neutrino data

$$R\frac{\mathrm{\Delta }m_{21}^2}{\mathrm{\Delta }m_{32}^2}\frac{5.0\times 10^5\mathrm{eV}^2}{2.5\times 10^3\mathrm{eV}^2}=2.0\times 10^2,$$
(47)
$`\mathrm{sin}^22\theta _{solar}=\mathrm{sin}^22\theta _{12}=0.76,`$ (48)
$`[\mathrm{tan}^2\theta _{solar}=0.34],`$ (49)
$$\mathrm{sin}^22\theta _{atm}=\mathrm{sin}^22\theta _{23}=1.0,$$
(50)
we must seek a parameter set that gives $`(M_\nu )_{22}(M_\nu )_{33}`$ and $`(M_\nu )_{12}(M_\nu )_{13}`$ for the expression (3.8). When we consider the Ai-type mixings, we obtain

$$M_\nu \left(\begin{array}{ccc}\epsilon ^2& \epsilon & \epsilon \\ \epsilon & 1& 1\\ \epsilon & 1& 1\end{array}\right),$$
(51)
for $`\epsilon =s_1^A/s_2^A=s_1^A/s_3^A`$ and $`\lambda _{2j3}/\lambda _{3j2}m_2^d/m_3^d`$. Generally, when we consider only Ai-type mixings, the solutions are highly dependent on the fine-tuning among the coefficients $`\lambda _{ijk}`$, and, besides, since the mass matrix is too near to a rank-1 matrix, it is difficult to give the solution a small but sizeable value of $`R`$ (it leads to an extremely small value of $`R`$). Even if we take the charged lepton loop contributions into consideration, the situation is not improved unless we assume the special parametrization for $`\lambda _{ijk}`$ ($`\lambda _{2j3}/\lambda _{3j2}m_2^d/m_3^d`$ and so on).

Next, we consider the case of the Bk-type mixings. When we consider a B1 mixing, we obtain a simple mass matrix form

$$M_\nu =(s_1^B)^2(m_1^d)^2\left(\begin{array}{ccc}0& 0& 0\\ 0& (\lambda _{211}^{})^2& \lambda _{211}^{}\lambda _{311}^{}\\ 0& \lambda _{211}^{}\lambda _{311}^{}& (\lambda _{311}^{})^2\end{array}\right),$$
(52)
because the case gives the relations $`\lambda _{ij2}^{}s_2^B=\lambda _{ij3}^{}s_3^B=0`$. It is natural to consider that $`\lambda _{211}^{}\lambda _{311}^{}`$ since those come from the same interactions (2.8) \[not from (2.5)\]. Therefore, the mass matrix (3.13) can give a maximal mixing between $`\nu _\mu `$ and $`\nu _\tau `$. A small additional term will reasonably give a bimaximal mixing.

For example, we consider a mixed-type, A3 and B1, mixing. In this case, the $`\lambda _{ijk}`$ in (3.3) must be replaced according to Table I. It is natural to consider that $`\lambda _{ijk}\mathrm{const}\lambda `$ and $`\lambda _{ijk}^{}\mathrm{const}\lambda ^{}`$. Then, the mass matrix $`M_\nu `$ is reduced to the form

$$M_\nu (s_1^B)^2(m_1^d)\left(\begin{array}{ccc}0& 0& \epsilon (r_3+r_2)\\ 0& 1& 1+\epsilon r_3\\ \epsilon (r_3+r_2)& 1+\epsilon r_3& (1+\epsilon r_2)^2\end{array}\right),$$
(53)
where $`\epsilon =s_3^A\lambda /s_1^B\lambda ^{}1`$, $`r_2=m_2^d/m_1^d`$ and $`r_3=m_3^d/m_1^d`$, and we have assumed $`c_3^A1`$. This mass matrix can give a reasonable value of $`R`$ together with a nearly bimaximal mixing. For example, for the parameter value $`\epsilon =0.000374`$, we obtain the following numerical results:

$$m_1^\nu =0.0822m_0,m_2^\nu =0.3276m_0,m_3^\nu =2.2604m_0,$$
(54)
$$U=\left(\begin{array}{ccc}0.8726& 0.4822& 0.0776\\ 0.3925& 0.5978& 0.6990\\ 0.2907& 0.6404& 0.7109\end{array}\right),$$
(55)
i.e.

$$R=0.0201,$$
(56)
$$\mathrm{sin}^22\theta _{12}4U_{11}^2U_{12}^2=0.708,$$
(57)
$$\mathrm{sin}^22\theta _{23}4U_{23}^2U_{33}^2=0.988.$$
(58)
These values are in good agreement with the best fit values (3.9)–(3.11) for the observed neutrino data. Although it is difficult in the original Zee model to give a sizeable deviation of $`\mathrm{sin}^22\theta _{12}`$ from 1  (it must be $`\mathrm{sin}^22\theta _{12}=1.0`$), the present model can give a reasonable deviation from $`\mathrm{sin}^22\theta _{12}=1.0`$. The result

$$U_{13}^2=0.00602$$
(59)
is also consistent with the present experimental upper limit

$$|U_{13}|^2<0.03,$$
(60)
from the CHOOZ collaboration .

## IV Higgs sectors

In the present model, the quark and charged lepton mass matrices are generated by the VEVs of the Higgs scalars with $`z=+1`$. Therefore, even if SU(2)L is broken later, the Z2 symmetry still exactly holds.

Let us show the mass matrices $`M_f`$ for the case of the mixed-type mixing A3 and $`B_1`$ as an example of the explicit forms of $`M_f`$:

$$M_e=\left(\begin{array}{ccc}c_1& b_1& a_1\\ c_2& b_2& a_2\\ c_1^Ac_3& c_1^Ab_3& c_1^Aa_3\end{array}\right),$$
(61)
$$M_d=\left(\begin{array}{ccc}c_1^Bc_1& c_1^Bc_2& c_1^Bc_3\\ b_1& b_2& b_3\\ a_1& a_2& a_3\end{array}\right),$$
(62)
$$M_u=\left(\begin{array}{ccc}c_1^Bc_1^{}& c_1^Bc_2^{}& c_1^Bc_3^{}\\ b_1^{}& b_2^{}& b_3^{}\\ a_1^{}& a_2^{}& a_3^{}\end{array}\right).$$
(63)
Here, $`M_f`$ have been defined by $`\overline{f}_LM_ff_R`$ ($`f=u,d,e`$). Note that the mass matrix $`M_d`$ has a form different from $`M_e^T`$ because of the mixing factors. Usually, if we consider one type of Higgs scalar of SU(5) 5-plet ($`\overline{5}`$-plet), it is difficult to obtain realistic mass matrices $`M_f`$ ($`f=u,d,e`$). Therefore, the present model has a possibility to improve this problem. However, whether we can give reasonable mass matrix forms of $`M_f`$ or not is a future task to us.

Now, we would like to give a comment on the Higgs sectors. In Sec. II, we have assumed that although the superfield $`\overline{5}_{L(+)}`$ is decoupled below $`\mu =M_X`$, the components $`(2,1)_{(+)}`$ of $`\overline{5}_{L(+)}`$ can contribute to low energy phenomena through the mixing (2.3). If we consider the Higgs fields $`\overline{H}_{d(+)}`$ and $`\overline{H}_{d()}`$ with SU(5) $`\overline{5}`$-plet, and if we assume a situation similar to the matter fields $`\overline{5}_L`$, then we obtain

$$\overline{H}_d=(2,1)_{(+)}+(1,\overline{3})_{()},$$
(64)
where we have assumed a perfect interchange between $`(2,1)_{()}`$ and $`(2,1)_{(+)}`$, not a mixing. Note that in this scheme, the Higgs scalar component $`(3,1)_{()}`$ cannot couple to the fermions $`\overline{d}_Ru_R^c`$ with $`z=+1`$ independently of the mixings Ai and $`B_k`$, so that the scalar $`(3,1)_{()}`$ cannot contribute to the proton decay and it need not be super-heavy. (Although it couples to the fermions $`(\overline{\nu }_L^cd_L\overline{e}_L^cu_L)`$, these interactions cannot contribute to the proton decay.) We can consider a similar mechanism for the Higgs fields $`H_u`$ with SU(5) 5-plet. The interactions of $`(1,3)_{()}`$ with $`\overline{u}_L^cd_L`$ and $`\overline{e}_Ru_R^c`$ are absent, so that the scalar $`(1,3)_{()}`$ need not be super-heavy. The so-called $`\mu `$-terms are composed of $`H_{u(+)}\overline{H}_{d(+)}`$ and $`H_{u()}\overline{H}_{d()}`$.

## V Conclusion

In conclusion, we have investigated possible forms of radiatively induced neutrino mass matrix under an ansatz \[(2.3) and (2.6)\] within the framework of the SU(5) SUSY model. We have assumed two types of matter fields $`\mathrm{\Psi }_L=\overline{5}_L+10_L`$, $`\mathrm{\Psi }_{(\pm )}`$, which are transformed as $`\mathrm{\Psi }_{(\pm )}\pm \mathrm{\Psi }_{(\pm )}`$ under the Z2 symmetry. We have assumed that the Z2 symmetry exactly holds, even if SU(5) is broken. The essential ansatz is in the mixings (2.3) and (2.6). Although the origin of the mixings (2.3) and (2.6) is still an open question, if we admit this ansatz, we can obtain very interesting and simple neutrino mass matrix form, which can give satisfactory numerical results for the observed neutrino data. How we can justify the ansatz is a future task for us.

Acknowledgements

The work based on the Z2 symmetry was first begun by Dr. A. Ghosal and one of the authors (Y. K.), but the earlier version  (hep-ph/0203113) failed to suppress the proton decay sufficiently. This work is also based on the Z2 symmetry, but the suppression mechanism is completely different from the previous one. Y. K. would like to thank Dr. A. Ghosal for his collaboration at the earlier stage. Y. K. also wishes to acknowledge the hospitality of the Theory group at CERN, where this work was completed. The work of J. S. is supported in part by Grants-in-Aid for Scientific Research from the Ministry of Education, Science, Sports, and Culture of Japan, No.14740168, No.14039209 and No.14046217.


format: text

text:
# Modernizing the ESRF beamline software architecture with generic Python modules

## Introduction

This article describes the new application software architecture on the ESRF beamlines. It will show its basic elements and comment on the level of their advancement.

The work has been mainly done in the BLISS group (Beamline Instrument Software Support) and by CSS (Certified Scientific Software) for the Spec Server. The list of contributers to this project include G. Berruyer, A. Gobbo, L. Claustre,T. Jouve, J. Klora, E. Papillon, N. Pascal, V. Rey, A. Sole, D. Spruce, and G. Swislow.

## The global beamline control software structure

The beamline software architecture at the ESRF is built around the commercial data acquisition package SPEC. SPEC talks to our hardware mainly via the ESRF TACO device servers. The device servers can either be specific to one type of device or just represent the hardware channel to talk to a device (i.e. a serial line). For a small number of device types (i.e stepper motors) which represent however the majority of hardware items used on the beamlines (i.e. we have about 3000 motors on our beamline) the code to talk to these device has been added to SPECs internal codebase. For the rest of devices the access is done via SPECs macro language. The end user does normally not see a difference between these methods.

The data display routines in SPEC are not sufficient for our beamlines. We have therefore decided to add external data visualization tools (i.e. for CCD cameras). The communication between the different programs is done using shared memory with specific header informations.

GUI user interface are put on top of SPEC. In the past we were using a program which mixed the standard user input from the keyboard with the input from the GUIs. SPEC sent information back to this program via a dedicated pipe. Multiple GUIs can be used at the same time. The information sent has been distributed to the GUIs in some type of event.

Two types of such GUIs have been used at the ESRF. There is a general SPEC GUI written in Motif which had a built in editor for graphical panels. The relatively simple panels were used to gather information to start SPEC macros. No programming had to be done in the GUI. There exists also a GUI written in Tcl/Tk dedicated to Protein crystallography experiments. This GUI contains much more dedicated logic but will also in the end send commands to SPEC to carry out the experiments. It will however do many things in addition (i.e. add information to a central data base).

## The new building blocks

The modernization of the control system keeps the general structure described in the last paragraph. Individual components are replaced with modules and applications written in Python. We can therefore distinguish the following areas:

* The external program which mixes keyboard input and the GUI communication to SPEC has been replaced by a new server mode in SPEC. When started in server mode, SPEC listens to requests from clients and sends events to them. We have a Python client library which provides access to SPEC with an object oriented interface.
* A plotting framework (PyDVT) has been created. The framework defines classes in the following areas:

> + visualization (like 2D images or contour plots)
>     + data sources (like input from a file or from SPECs shared memory)
>     + user selections (like cuts)
>     + filters (like fitting routines or color selectors)

It makes it then possible to define the data flow between them.
* Applications for visualization have been created from this framework. They work either standalone with data from files or integrated in the online data visualization with data from SPEC shared memory. Examples of these applications are given later.
* A new configuration editor has been created. This editor can be freely configured. It integrates the functionality of SPECs config editor and our own in-house configuration tools.
* A general GUI framework application has been started in Python. The application allows to load Python classes and integrate them into the running system. There is an built in editor to create simpler panels interactively without any programming.

## Why Python

Our choice for Python was relatively simple to make. We wanted a scripting language which can be easily used by the scientists on the beamline. Python can be extended with modules in Python or other languages (i.e. C). The language is very compact, easy to read, and offers an enormous amount of builtin functionality. It is largely platform independent. You can read more on www.python.org.

Unluckily the choice of the toolkit for Python was less obvious. We started using the most widely used toolkit Tk (coming from Tcl/Tk). We did have to program many standard widgets for our projects and the final programs didn’t have the ”look and feel” we wanted. We decided therefore to try Qt with Python bindings (www.trolltech.com) and are very happy with the functionality. Python/Qt is however not very widely used. We tried therefore to be as independent of the graphical toolkit as possible in the basic classes.

## The Spec Server

SPEC when started in server mode accepts connections on a socket. Multiple clients can connect and send messages. Clients can execute commands, talk to motor or counter objects in the server, read variables (i.e. global variables or motor parameters) from the server, and be automatically informed about changes of these variables. Events can be sent from SPEC to the clients.

The SPEC server does execute one macro at a time. It will queue macro executions if the macro unit is currently occupied. This restriction is the easiest way to avoid uncontrolled mixing of experimental procedures. There are three points why this restriction is not very important for practical purposes.

1. Not many incoming requests will lead to the execution of a macro. For example, it is possible to ask the value of a global variable while a scan macro is executing.
2. The event nature of the communication leads to a reduction in the client server traffic. The client is automatically informed when something interesting for him changes (i.e. a motor stops).
3. The level of commands. The client tells the server for example to start a movement but will then be free to do other things. A graphical user interface can start a movement on a motor, start later another movement of another motors and all that while the server is counting on a detector.

The SPEC server and client can be started and stopped independently and will reconnect automatically. There currently exists a Python client library and other SPEC version an be used as clients too. It is for example possible to configure a motor in a client SPEC version to refer to a motor in a SPEC server. One client can talk to multiple servers.

The development of the SPEC server mode has been accompanied by a new method of writing pseudo motors and counters. Only these pseudo motors will be fully functional in server mode. The new method consists of calling predefined macro code from SPECs built in C code.

The server mode and pseudo motors are already pretty stable with only few functionalities missing.

## The plotting framework

PyDVT is a framework for visualization tools. It defines different objects and the data flow through them. The data flows from the ”source” to the graphical ”view” window. On its way it passes through ”filters”. There are special selection classes to select which data to take from the source and also selection classes which refer to the graphical selection done by the user with the mouse.

The views can be 2D display, contour plots, normal 1D plots, radial plots, 3D plot, data tables, or mesh plots. The underlying graphics kit is PLPLOT to which we added user interaction. Online display of fast 2D detectors requires a level of performance we could not get in this way. We therefore used a specificly written library in C and a pixmap widget for the 2D display.

For the moment the data can come from either files or SPEC shared memory. The source also includes a trigger mechanism when the data has been updated. The user of the tool kit does not have to be concerned with the update of the plot in this way.

It is relative easy to insert filters into the data-flow. There are filter fit-routines which take input from the data source and result in the fitted curve (which can then be displayed in a view possibly together with the raw data)

The toolkit is independent of the GUI toolkit. We currently use bindings for Qt and Tk.

## Visualization Applications

General graphical applications have been created with the Python toolkit. Examples of these applications are:

PyMCA is an application to work with spectra from multi channel analyzers. It can read the spectra from SPEC files or SPEC shared memory and display them. One can calibrate the spectra in user units (normally energy). Peak search and fitting are integrated.

PyDis takes two dimensional data from files or SPEC shared memory and displays it in different ways. The application can be extended with plug-ins. A poster presented at the workshop shows an example of PyDis extended for data analysis ()

Displays online or offline SPEC data from scans. Allows fitting of the data and all the common view operations.

The applications are finished but not still need some validating before we can use them on our beamlines on a large scale.

## A Configuration Editor

A common problem on our beamlines is the complexity of the configuration of its components. To add a device to the system we often have to add information in many different places. We decided therefore to create a configuration editor which could be easily adapted to input all the necessary information. It will replace the SPEC configuration editor and ESRF specific tools for configuration at our beamlines. It is also planned to use the editor in the NICOS control program of the Munich FRM II beamlines.

It consists of two parts:

1. A description editor which creates files which defines all the possible items to configure. It allows to create ”properties” (i.e. the motor parameter ”velocity”, or the unit). Every property has a type and a default value. One can then define the properties of the objects in the system (i.e. that all motor objects have a ”velocity” property).
2. A configuration editor which will use the information from the description editor to create the above objects and configure their properties. The editor allows to create groups (i.e. diffractometer or monochromator), enable and disable parts of the configuration, and generally assists the user during the configuration.

The program is in the final development phase. We still need to provide output filters to save the final configuration in the standard (”old”) file formats (i.e. SPEC config file).

## The general GUI framework

The general GUI framework is an attempt to integrate different standard graphical components with beamline specific parts. Starting from an empty application the individual elements are loaded interactively. This configuration is then saved. The general concept is born out of the experience from our generic SPEC GUI and specific GUI applications (i.e. ProDC for the protein crystallography beamlines). Both approaches had their advantages and disadvantages. While the specific interface pleased which its rich functionality we have to invest a lot of time to maintain it. We want to combine complicated graphical standard components with beamline specific panels. The framework consists therefore of two parts:

1. A wizard to load in Python classes into the application. These classes need to be Qt or Tk Widgets (both at the same time don’t work well at the moment). There is an optional mechanism to communicate with other classes with an ”eventhandler” (classes can send and receive events).

A special type of classes is the placer type class. These ”placers” once loaded into the applications accept other classes inside their window in a predefined way (a up/down placer might accept two windows it will place either in the upper or lower half of its window).

The wizard will therefore go through the following steps:

> * ask a filename and make the user select a Python class
>     * ask the user where to put it on the screen (either in its own window or inside another one)
>     * ask additional parameters and create and instance
2. An integrated editor for simple panels with simple graphical elements like input widgets, buttons, frames, labels, and so on. These panels will display information from SPEC and send commands to SPEC.

This project is in its final stage but needs some more work.

## FAQ

We need to provide fast (sometimes in hours) a 90% solution. Our highest priority is to make an experiment possible however inconvenient. We will later try to automate certain steps and only later add a graphical interface for convenience.

Apart from the obvious reason that SPEC had been more or less our only option when we started, there is a very import basic philosophy behind this decision - ”Configure - don’t program”. SPEC allows us to hide all the different server implementations behind a standard layer. It allows us to adapt the system to the different experimental stations. Writing a specific application for an experimental station is relatively simple - but it is extremely difficult to maintain these type of programs. This approach becomes also inefficient with the large number of experimental stations as at the ESRF.

## Conclusion

Python has been proven an excellent tool. It allows to quickly develop complex programs and helps to collaborate in the our group. On the down side we see that the installation of the finished software takes more effort and regret the absence of a powerful standard GUI toolkit. Our overall opinion stays however very positive.

The new concept is finished but some building blocks still need some more work. News about our progress will be available on http://www.esrf.fr/computing/bliss/bliss.html .


text:
# 1 𝑃_Λ^{𝑝𝑝} at √𝑆=62 GeV.

Hyperon Polarization from Unpolarized $`pp`$ and $`ep`$ Collisions

Yuji Koike

Department of Physics, Niigata University, Ikarashi, Niigata 950–2181, Japan

Abstract: Cross section formulas for the $`\mathrm{\Lambda }`$ polarization in $`pp\mathrm{\Lambda }^{}(\mathrm{}_T)X`$ and $`ep\mathrm{\Lambda }^{}(\mathrm{}_T)X`$ are derived and its characteristic features are discussed.

In this report we discuss the polarization of $`\mathrm{\Lambda }`$ hyperon produced in unpolarized $`pp`$ and $`ep`$ collisions relevant for the ongoing RHIC-SPIN, HERMES and COMPASS experiments. According to the QCD factorization theorem, the polarized cross section for $`pp\mathrm{\Lambda }^{}X`$ consists of two twist-3 contributions:

$`(A)`$ $`E_a(x_1,x_2)q_b(x^{})\delta \widehat{q}_c(z)\widehat{\sigma }_{abc},`$
$`(B)`$ $`q_a(x)q_b(x^{})\widehat{G}_c(z_1,z_2)\widehat{\sigma }_{abc}^{},`$
where the functions $`E_a(x_1,x_2)`$ and $`\widehat{G}_c(z_1,z_2)`$ are the twist-3 quantities representing, respectively, the unpolarized distribution in the nucleon and the polarized fragmentation function for $`\mathrm{\Lambda }^{}`$. $`\delta \widehat{q}_c(x)`$ is the transversity fragmentation function for $`\mathrm{\Lambda }^{}`$. $`a`$, $`b`$ and $`c`$ stand for the parton’s species, sum over which is implied. $`E_a`$ and $`\delta \widehat{q}_c`$ are chiral-odd. Corresponding to the above (A) and (B), the polarized cross section for $`ep\mathrm{\Lambda }^{}X`$ (final electron is not detected) receives two twist-3 contributions:

$`(A^{})`$ $`E_a(x_1,x_2)\delta \widehat{q}_a(z)\widehat{\sigma }_{eaa},`$
$`(B^{})`$ $`q_a(x)\widehat{G}_a(z_1,z_2)\widehat{\sigma }_{eaa}^{}.`$
The (A) contribution for $`pp\mathrm{\Lambda }^{}X`$ has been analyzed in , where it was shown that (A) gives rise to growing $`P_\mathrm{\Lambda }`$ at large $`x_F`$ as observed experimentally. Here we extend the study to the (B) term (see also ) at RHIC energy and also for the $`ep`$ collision.

The unpolarized twist-3 distribution $`E_{F,D}(x_1,x_2)`$ is defined in . Likewise the twist-3 fragmentation function for a polarized $`\mathrm{\Lambda }`$ (with momentum $`\mathrm{}`$) is defined as the lightcone correlation function as ($`w^2=0`$, $`\mathrm{}w=1`$)

$`{\displaystyle \frac{1}{N_c}}{\displaystyle \underset{X}{}}{\displaystyle \frac{d\lambda }{2\pi }\frac{d\mu }{2\pi }e^{i\frac{\lambda }{z_1}}e^{i\mu (\frac{1}{z_2}\frac{1}{z_1})}0|\psi _i(0)|\pi X\pi X|gF^{\alpha \beta }(\mu w)w_\beta \overline{\psi }_j(\lambda w)|0}`$
$`={\displaystyle \frac{M_N}{2z_2}}\left(\text{/}\mathrm{}\right)_{ij}ϵ^{\alpha \mathrm{}wS_{}}\widehat{G}_F(z_1,z_2)+i{\displaystyle \frac{M_N}{2z_2}}\left(\gamma _5\text{/}\mathrm{}\right)_{ij}S_{}^\alpha \widehat{G}_F^5(z_1,z_2)+\mathrm{}.`$ (1)
Note that we use the nucleon mass $`M_N`$ to normalize the twist-3 fragmentation function for $`\mathrm{\Lambda }`$. There is another twist-3 fragmentation functions which are obtained from (1) by shifting the gluon-field strength from the left to the right of the cut. The defined functions $`\widehat{G}_{FR}(z_1,z_2)`$ and $`\widehat{G}_{FR}^5(z_1,z_2)`$ are connected to $`\widehat{G}_F(z_1,z_2)`$ by the relation $`\widehat{G}_F(z_1,z_2)=\widehat{G}_{FR}(z_2,z_1)`$ and $`\widehat{G}_F^5(z_1,z_2)=\widehat{G}_{FR}^5(z_2,z_1)`$, which follows from hermiticity and time reversal invariance. Unlike the twist-3 distributions, the twist-3 fragmentation function does not have definite symmetry property. Another class of twist-3 fragmentation functions $`\widehat{G}_D^{(5)}(z_1,z_2)`$ is also defined from (1) by replacing $`gF^{\alpha \beta }(\mu w)w_\beta `$ by $`D^\alpha (\mu w)=^\alpha igA(\mu w)`$. Note, however, this is not independent from the above (1).

Following the method of  we present the analysis of the (C) term. The detailed analysis shows $`\widehat{G}_F(z,z)`$ appears as soft-gluon-pole contribution ($`z_1=z_2=z`$), while $`\widehat{G}_D(z_1,z_2)`$ appears as a soft fermion pole ($`z_1=0`$ or $`z_2=2`$). Physically, the latter is expected to be suppressed, and we include only the former contribution. This observation also applies to $`E_{F,D}(x_1,x_2)`$ relevant for the (A) term. In the large $`x_F`$ region, the main contribution comes from large-$`x`$ and large-$`z`$ (and small $`x^{}`$) region. Since $`E_F`$ and $`\widehat{G}_F`$ behaves as $`E_F(x,x)(1x)^\beta `$ and $`\widehat{G}_F(z,z)(1z)^\beta ^{}`$ with $`\beta `$, $`\beta ^{}>0`$, $`|(d/dx)E_F(x,x)||E_F(x,x)|`$, $`|(d/dz)\widehat{G}_F(z,z)||\widehat{G}_F(z,z)|`$ at large $`x`$ and $`z`$. In particular, the valence component of $`E_F`$ and $`\widehat{G}_F`$ dominates in this region. We thus keep only the valence quark contribution for the derivative of these soft-gluon pole function (“valence-quark soft-gluon approximation”) for the $`pp`$ collision. For the $`ep`$ case, we include all the soft-gluon pole contribution, since the calculation is relatively simple compared to the $`pp`$ case.

In general $`P_\mathrm{\Lambda }`$ is a function of $`S=(P+P^{})^22PP^{}`$, $`T=(P\mathrm{})^22P\mathrm{}`$ and $`U=(P^{}\mathrm{})^22P^{}\mathrm{}`$ where $`P`$ and $`P^{}`$ are the momenta of the two nucleons, and $`\mathrm{}`$ is the momentum of $`\mathrm{\Lambda }`$. In the following we use $`S`$, $`x_F=\frac{2\mathrm{}_{}}{\sqrt{S}}=\frac{TU}{S}`$ and $`x_T=\frac{2\mathrm{}_T}{\sqrt{S}}`$ as independent variables. The polarized cross section for the (B) term reads

$`E_\mathrm{\Lambda }{\displaystyle \frac{d^3\mathrm{\Delta }\sigma (S_{})}{d\mathrm{}^3}}={\displaystyle \frac{2\pi M_N\alpha _s^2}{S}}ϵ^{\alpha \mathrm{}wS_{}}{\displaystyle \underset{a}{}}{\displaystyle _{z_{min}}^1}{\displaystyle \frac{dz}{z^2}}{\displaystyle _{x_{min}}^1}{\displaystyle \frac{dx}{x}}{\displaystyle \frac{1}{xS+U/z}}{\displaystyle _0^1}{\displaystyle \frac{dx^{}}{x^{}}}`$
$`\times \delta \left(x^{}+{\displaystyle \frac{xT/z}{xS+U/z}}\right)`$
$`\times \{{\displaystyle \underset{b,c}{}}q^a(x)q^b(x^{})[z_1^2{\displaystyle \frac{}{z_1}}\widehat{G}_F^a(z_1,z)]_{z_1=z}({\displaystyle \frac{2p_\alpha }{T}}\widehat{\sigma }_{abc}^I+{\displaystyle \frac{2p_\alpha ^{}}{U}}\widehat{\sigma }_{abc}^{II})`$
$`+{\displaystyle \underset{b,c}{}}q^a(x)q^b(x^{})\left[z^2{\displaystyle \frac{d}{dz}}\widehat{G}_F^a(z,z)\right]{\displaystyle \frac{xp_\alpha +x^{}p_\alpha ^{}}{|xT+x^{}U|}}\left(\widehat{\sigma }_{abc}^I+\widehat{\sigma }_{abc}^{II}\right)`$
$`+q^a(x)G(x^{})\left[z_1^2{\displaystyle \frac{}{z_1}}\widehat{G}_F^a(z_1,z)\right]_{z_1=z}\left({\displaystyle \frac{2p_\alpha }{T}}\widehat{\sigma }_{aga}^I+{\displaystyle \frac{2p_\alpha ^{}}{U}}\widehat{\sigma }_{aga}^{II}\right)`$
$`+q^a(x)G(x^{})[z^2{\displaystyle \frac{d}{dz}}\widehat{G}_F^a(z,z)]{\displaystyle \frac{xp_\alpha +x^{}p_\alpha ^{}}{|xT+x^{}U|}}(\widehat{\sigma }_{aga}^I+\widehat{\sigma }_{aga}^{II})\},`$ (2)
where the lower limits for the integration variables are $`z_{min}=(T+U)/S=\sqrt{x_F^2+x_T^2}`$ and $`x_{min}=U/z(S+T/z)`$. The partonic hard cross sections are written in terms of the invariants in the parton level, $`\widehat{s}=(xp+x^{}p^{})^2=xx^{}S`$, $`\widehat{t}=(xp\mathrm{}/z)^2=xT/z`$ and $`\widehat{u}=(x^{}p^{}\mathrm{}/z)^2=x^{}U/z`$. They read

$`\widehat{\sigma }_{abc}^I={\displaystyle \frac{1}{36}}{\displaystyle \frac{\widehat{s}^2+\widehat{u}^2}{\widehat{t}^2}}\delta _{ac}+{\displaystyle \frac{7}{36}}{\displaystyle \frac{\widehat{s}^2+\widehat{t}^2}{\widehat{u}^2}}\delta _{bc}+{\displaystyle \frac{1}{54}}{\displaystyle \frac{\widehat{s}^2}{\widehat{t}\widehat{u}}}\delta _{ab}\delta _{ac},`$
$`\widehat{\sigma }_{abc}^{II}={\displaystyle \frac{7}{36}}{\displaystyle \frac{\widehat{s}^2+\widehat{u}^2}{\widehat{t}^2}}\delta _{ac}{\displaystyle \frac{1}{36}}{\displaystyle \frac{\widehat{s}^2+\widehat{t}^2}{\widehat{u}^2}}\delta _{bc}+{\displaystyle \frac{1}{54}}{\displaystyle \frac{\widehat{s}^2}{\widehat{t}\widehat{u}}}\delta _{ab}\delta _{ac},`$
$`\widehat{\sigma }_{a\overline{b}c}^I={\displaystyle \frac{1}{36}}{\displaystyle \frac{\widehat{s}^2+\widehat{u}^2}{\widehat{t}^2}}\delta _{ac}+{\displaystyle \frac{7}{36}}{\displaystyle \frac{\widehat{u}^2+\widehat{t}^2}{\widehat{s}^2}}\delta _{ab},\widehat{\sigma }_{a\overline{b}c}^{II}={\displaystyle \frac{1}{18}}{\displaystyle \frac{\widehat{s}^2+\widehat{u}^2}{\widehat{t}^2}}\delta _{ac}+{\displaystyle \frac{1}{18}}{\displaystyle \frac{\widehat{u}^2+\widehat{t}^2}{\widehat{s}^2}}\delta _{ab},`$
$`\widehat{\sigma }_{qgq}^I={\displaystyle \frac{1}{8}}\left(1{\displaystyle \frac{\widehat{s}\widehat{u}}{\widehat{t}^2}}\right)+{\displaystyle \frac{1}{288}}\left({\displaystyle \frac{\widehat{u}}{\widehat{s}}}+{\displaystyle \frac{\widehat{s}}{\widehat{u}}}\right){\displaystyle \frac{\widehat{s}}{16\widehat{t}}}{\displaystyle \frac{\widehat{u}}{16\widehat{t}}},`$
$`\widehat{\sigma }_{qgq}^{II}={\displaystyle \frac{9}{16}}\left(1{\displaystyle \frac{\widehat{s}\widehat{u}}{\widehat{t}^2}}\right)+{\displaystyle \frac{\widehat{u}}{32\widehat{s}}}{\displaystyle \frac{\widehat{s}}{4\widehat{u}}}+{\displaystyle \frac{9\widehat{u}}{16\widehat{t}}}.`$ (3)
Among these partonic cross sections, $`\widehat{\sigma }^I`$ becomes more important at large $`x_F`$ because of the $`1/T`$ factor in (2).

To estimate the above contribution, we introduce a model ansatz as $`\widehat{G}_F^a(z,z)=K_a\widehat{q}_a(z)`$ with twist-2 unpolarized fragmentation function $`\widehat{q}^a(z)`$, noting that the Dirac structure of $`\widehat{G}_F^a(z,z)`$ and $`\widehat{q}_a(z)`$ is the same . $`K_a`$’s are taken to be $`K_u=K_d=0.07`$ which are the same values used in the relation $`G_F(x,x)=K_aq^a(x)`$ to reproduce $`A_N`$ in $`p^{}p\pi X`$ observed at E704 . As noted before, $`\widehat{G}_F(z_1,z_2)`$ does not have definite symmetry property unlike the twist-3 distribution $`E_F(x_1,x_2)`$. Nevertheless we assume $`\left[(/z_1)\widehat{E}_F(z_1,z)\right]_{z_1=z}=(1/2)(d/dz)\widehat{E}_F(z,z)`$. The result for the $`\mathrm{\Lambda }`$ polarization $`P_\mathrm{\Lambda }^{pp}`$ at $`\sqrt{S}=62`$ GeV is shown in Fig. 1 together with the R608 data. There (A) (chiral-odd) contribution studied in  is also shown for comparison. (For the adopted distribution and fragmentation functions, see .) One sees that the tendency of $`P_\mathrm{\Lambda }^{pp}`$ from the (B)(chiral-even) contribution is quite similar to the R608 data. Rising behavior of $`P_\mathrm{\Lambda }^{pp}`$ at large $`x_F`$ comes from (i) the large partonic cross sections in (3) ($`1/\widehat{t}^2`$ term) and (ii) the derivative of the soft-gluon pole functions. With these parameters $`K_a`$, $`P_\mathrm{\Lambda }^{pp}`$ at RHIC energy ($`\sqrt{S}=200`$ GeV) is shown in Fig.2 at $`l_T=1.5`$ GeV. Fig. 3 shows the $`l_T`$ dependence of $`P_\mathrm{\Lambda }^{pp}`$ of the (B) term, indicating large $`\mathrm{}_T`$ dependence at $`1\mathrm{}_T3`$ GeV. Experimentally, $`P_\mathrm{\Lambda }^{pp}`$ grows up as $`\mathrm{}_T`$ increases up to $`\mathrm{}_T1`$ GeV and stays constant at $`1\mathrm{}_T3`$ GeV. So the $`P_\mathrm{\Lambda }^{pp}`$ observed at R608 can not be wholly ascribed to the twist-3 effect studied here which is designed to describe large $`\mathrm{}_T`$ polarization.

We next discuss the polarization $`P_\mathrm{\Lambda }^{ep}`$ in $`pe\mathrm{\Lambda }^{}(\mathrm{})X`$ where the final electron is not observed. In our $`O(\alpha _s^0)`$ calculation, the exchanged photon remains highly virtual as far as the observed $`\mathrm{\Lambda }`$ has a large transverse momentum $`\mathrm{}_T`$ with respect to the $`ep`$ axis. Therefore experimentally one needs to integrates only over those virtual photon events to compare with our formula.

Using the twist-3 distribution and fragmentation functions used to describe $`P_\mathrm{\Lambda }^{pp}`$, we show in Fig. 4 the obtained $`P_\mathrm{\Lambda }^{ep}`$ corresponding to (A’)(chiral-odd) and (B’)(chiral-even) contributions. Remarkable feature of Fig. 4 is that in both chiral-even and chiral-odd contributions (i) the sign of $`P_\mathrm{\Lambda }^{ep}`$ is opposite to the sign of $`P_\mathrm{\Lambda }^{pp}`$ and (ii) the magnitude of $`P_\mathrm{\Lambda }^{ep}`$ is much larger than that of $`P_\mathrm{\Lambda }^{pp}`$, in particular, at large $`x_F`$, and it even overshoots one. (In our convention, $`x_F>0`$ corresponds to the production of $`\mathrm{\Lambda }`$ in the forward hemisphere of the initial proton in the $`ep`$ case.) The origin of these features can be traced back to the color factor in the dominant diagrams for the twist-3 polarized cross sections in $`ep`$ and $`pp`$ collisions. Of course, the $`P_\mathrm{\Lambda }`$ can not exceeds one, and thus our model estimate needs to be modified. First, the applied kinematic range of our formula should be reconsidered: Application of the twist-3 cross section at such small $`\mathrm{}_T`$ may not be justified. Second, our simple model ansatz of $`E_F^a(x,x)\delta q^a(x)`$ (in (A) term) and $`\widehat{G}_F^a(z,z)\widehat{q}^a(z)`$ should be modified at $`x1`$ and $`z1`$, respectively. The derivative of these functions, which is important for the growing $`P_\mathrm{\Lambda }^{pp}`$ at large $`x_F`$, eventually leads to divergence of $`P_\mathrm{\Lambda }`$ at $`x_F1`$ as $`1/(1x_F)`$.

As a possible remedy for this pathology we tried the following: As an example for the (B) (chiral-even) contribution we have a model $`\widehat{G}_F^a(z,z)\widehat{q}_a(z)_{z1}(1z)^\beta `$ where $`\beta =1.83`$ in the fragmentation function we adopted. Tentatively we shifted $`\beta `$ as $`\beta \beta (z)=\beta +z^8`$, which suppresses the divergence of $`P_\mathrm{\Lambda }`$ at $`x_F1`$ but still keeps rising behavior of $`P_\mathrm{\Lambda }`$ at large $`x_F`$. This avoids overshooting of one in $`P_\mathrm{\Lambda }^{ep}`$ but reduces $`P_\mathrm{\Lambda }^{pp}`$ seriously. The result obtained by this modification is shown in Figs. 5 and 6.

To summarize we have studied the $`\mathrm{\Lambda }`$ polarization in $`pp`$ and $`ep`$ collisions in the framework of collinear factorization. Our approach includes all effects for the large $`\mathrm{}_T`$ production. One needs to be causious in interpreting the available $`pp`$ data at relatively low $`\mathrm{}_T`$ in terms of the derived formula. Determination of the participating twist-3 functions requires global analysis of future $`pp`$ and $`ep`$ data.

Acknowledgement: This work is supported in part by the Grant-in-Aid for Scientific Research of Monbusho.


text:
# I Introduction

## I Introduction

The spontaneous polarization in ferroelectrics can be switched by applying a strong electric field opposite to the polarization direction . It is crucial to understand the details of this switching process due to many applications in memory devices and electromechanical systems. It is now a well established fact that the switching behavior is strongly influenced by nucleation, growth and motion of domain walls. Due to the nucleation events, the experimentally measured coercive field is much smaller than the intrinsic value predicted by the simple 1-D Landau theory . In a recent work , we qualitatively showed that dipolar defects can act as nucleation centres for $`90^o`$ domain structures. However, for a constrained system, such as a grain in a ferroelectric ceramic, the surface itself is a heterogeniety that can cause localized nucleation. Although, there have been some studies investigating the role of a surface on ferroelectric phase transitions , the polarization switching in finite sized systems has not been well studied theoretically.

In this paper, we study the influence of a surface during polarization switching using the time dependent Ginzburg-Landau model described in . In particular, we are interested in understanding the role played by the surface in nucleating $`90^o`$ ferroelastic domains during the switching process. Recent experiments on $`BaTiO_3`$ single crystals  have demonstrated the importance of $`90^o`$ domains during polarization reversal. The theory incorporates elastic and electrostrictive effects in the form of an anisotropic long-range interaction obtained by eliminating strain fields, subject to the elastic compatibility requirement. Instead of using periodic boundary conditions (corresponding to a bulk system), we consider a finite sized constrained ferroelectric with the boundary condition of vanishing polarization at the surface. Physically, this corresponds to a system where the free charge at the surface is compensated. It is believed that this charge compensation occurs at the grain boundaries in ferroelectric ceramics. Due to the above discussed boundary conditions, we can neglect the effects of the depolarization field. The nonlocal interaction is represented by the appropriate polarization gradients terms in the free energy expansion.

The organization of this paper is as follows. In section 2, we give details of the time dependent Ginzburg-Landau model, starting from a free-energy functional for the polarization field. Section 3 gives details of our simulations of polarization switching for a finite square shaped system. In section 4, we describe simulations of switching behavior in a semi-infinite constrained system. We end the paper by a summary and discussion of our results in section 5.

## II Model

The 3-D Ginzburg-Landau free energy for ferroelectric systems has been given in Ref. . In this work, we restrict to a 2-D system undergoing a square to rectangle transition, for which the total free-energy is given by

$$F=𝑑\stackrel{}{r}\left[f_l+f_g+f_{el}+f_{es}+f_{ext}\right]$$
(1)
where $`f_l`$ is the local free-energy density given by

$`f_l={\displaystyle \frac{\alpha _1}{2}}(P_{x}^{}{}_{}{}^{2}+P_{y}^{}{}_{}{}^{2})+{\displaystyle \frac{\alpha _{11}}{4}}(P_{x}^{}{}_{}{}^{4}+P_{y}^{}{}_{}{}^{4})+{\displaystyle \frac{\alpha _{12}}{2}}P_{x}^{}{}_{}{}^{2}P_{y}^{}{}_{}{}^{2},`$ (2)
where $`P_x`$ and $`P_y`$ are the components of the polarization vector, $`\alpha _1`$, $`\alpha _{11}`$ and $`\alpha _{12}`$ are the free energy parameters. The term $`f_g`$ represents the gradient energy,

$`f_g={\displaystyle \frac{g_1}{2}}\left[\left({\displaystyle \frac{P_x}{x}}\right)^2+\left({\displaystyle \frac{P_y}{y}}\right)^2\right]+{\displaystyle \frac{g_2}{2}}\left[\left({\displaystyle \frac{P_x}{y}}\right)^2+\left({\displaystyle \frac{P_y}{x}}\right)^2\right]+{\displaystyle \frac{g_3}{2}}\left({\displaystyle \frac{P_x}{x}}\right)\left({\displaystyle \frac{P_y}{y}}\right),`$ (3)
where $`g_1`$, $`g_2`$ and $`g_3`$ are the gradient parameters. It has been shown that the gradient energy may be obtained as a local approximation to the nonlocal interactions . The term $`f_{el}`$ represents the usual elastic energy of the square system. We consider the bulk strain $`\varphi _1=(\eta _{xx}+\eta _{yy})/\sqrt{2}`$, deviatoric strain $`\varphi _2=(\eta _{xx}\eta _{yy})/\sqrt{2}`$ and shear strain $`\varphi _3=\eta _{xy}=\eta _{yx}`$. Here $`\eta _{ij}`$ is the linear elastic strain tensor given as $`\eta _{ij}=\frac{1}{2}(\frac{u_i}{x_j}+\frac{u_j}{x_i})`$, where $`u_i`$, $`(i=x,y)`$ represents the dispacement components. The elastic free-energy can then be written as

$`f_{el}={\displaystyle \frac{a_1}{2}}\varphi _{1}^{}{}_{}{}^{2}+{\displaystyle \frac{a_2}{2}}\varphi _{2}^{}{}_{}{}^{2}+{\displaystyle \frac{a_3}{2}}\varphi _{3}^{}{}_{}{}^{2},`$ (4)
where $`a_1`$, $`a_2`$ and $`a_3`$ are constants that can be expressed in terms of linear combinations of elastic constants of the system. The coupling energy $`f_{es}`$ may be written as,

$`f_{es}=q_1\varphi _1(P_{x}^{}{}_{}{}^{2}+P_{y}^{}{}_{}{}^{2})q_2\varphi _2(P_{x}^{}{}_{}{}^{2}P_{y}^{}{}_{}{}^{2})q_3\varphi _3P_xP_y.`$ (5)
Here $`q_1`$, $`q_2`$ and $`q_3`$ are constants that are related to the electrostrictive constants of the system. The term $`f_{ext}`$ is given by $`f_{ext}=\stackrel{}{E}_{ext}\stackrel{}{P}.`$ The amplitude of field $`\stackrel{}{E}_{ext}`$ is a tunable parameter in our simulations for studying polarization switching. We assume that the system reaches mechanical equilibrium very fast so that we may integrate out the elastic fields, subject to the elastic compatibility constraint . In terms of the strain components $`\varphi _1`$, $`\varphi _2`$ and $`\varphi _3`$, the elastic compatibility relation is given as

$$^2\varphi _1(\frac{^2}{x^2}\frac{^2}{y^2})\varphi _2\sqrt{8}\frac{^2}{xy}\varphi _3=0.$$
(6)
To incorporate the elastic compatibility constraint, we consider an effective elastic part of the free energy $`F_{eff}=F_{el}+F_{es}+F_{cons}`$, with Lagrangian multiplier $`\lambda `$, where

$$F_{el}=𝑑\stackrel{}{k}\left\{\frac{a_1}{2}|\varphi _1(\stackrel{}{k})|^2+\frac{a_2}{2}|\varphi _2(\stackrel{}{k})|^2+\frac{a_3}{2}|\varphi _3(\stackrel{}{k})|^2\right\},$$
(7)
$`F_{es}={\displaystyle 𝑑\stackrel{}{k}\left\{q_1\mathrm{\Gamma }_1(\stackrel{}{k})\varphi _1(\stackrel{}{k})+q_2\mathrm{\Gamma }_2(\stackrel{}{k})\varphi _2(\stackrel{}{k})+q_3\mathrm{\Gamma }_3(\stackrel{}{k})\varphi _3(\stackrel{}{k})\right\}}`$ (8)
$`F_{cons}={\displaystyle 𝑑\stackrel{}{k}\left\{\lambda (\stackrel{}{k})(k^2\varphi _1(\stackrel{}{k})+(k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2})\varphi _2(\stackrel{}{k})+\sqrt{8}k_xk_y\varphi _3(\stackrel{}{k}))\right\}}.`$ (9)
Here $`k^2=k_{x}^{}{}_{}{}^{2}+k_{y}^{}{}_{}{}^{2}`$ and $`\mathrm{\Gamma }_1(\stackrel{}{k})`$ , $`\mathrm{\Gamma }_2(\stackrel{}{k})`$ and $`\mathrm{\Gamma }_3(\stackrel{}{k})`$ are, respectively, the Fourier transforms of $`P_{x}^{}{}_{}{}^{2}+P_{y}^{}{}_{}{}^{2}`$, $`P_{x}^{}{}_{}{}^{2}P_{y}^{}{}_{}{}^{2}`$, and $`P_xP_y`$. The condition of mechanical equilibrium and the compatibility relation demands $`\frac{\delta F_{eff}}{\delta \varphi _i}=0`$$`(i=1,2,3)`$ and $`\frac{\delta F_{eff}}{\delta \lambda }=0`$.

$`a_1\varphi _1(\stackrel{}{k})`$ $`=`$ $`q_1\mathrm{\Gamma }_1(\stackrel{}{k})+k^2\lambda (\stackrel{}{k})`$ (10)
$`a_2\varphi _2(\stackrel{}{k})`$ $`=`$ $`q_2\mathrm{\Gamma }_2(\stackrel{}{k})(k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2})\lambda (\stackrel{}{k})`$ (11)
$`a_3\varphi _3(\stackrel{}{k})`$ $`=`$ $`q_3\mathrm{\Gamma }_3(\stackrel{}{k})\sqrt{8}k_xk_y\lambda (\stackrel{}{k}),`$ (12)
subject to the constraint

$$k^2\varphi _1(\stackrel{}{k})+(k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2})\varphi _2(\stackrel{}{k})+\sqrt{8}k_xk_y\varphi _3(\stackrel{}{k})=0.$$
(13)
Using the above equations, the multiplier $`\lambda (\stackrel{}{k})`$ is given as

$$\lambda (\stackrel{}{k})=\frac{k^2q_1\mathrm{\Gamma }_1(\stackrel{}{k})/a_1+(k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2})q_2\mathrm{\Gamma }_2(\stackrel{}{k})/a_2+\sqrt{8}k_xk_yq_3\mathrm{\Gamma }_3(\stackrel{}{k})/a_3}{k^4/a_1+(k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2})^2/a_2+8k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2}/a_3}.$$
(14)
At this stage we make a simplification by assuming that the polarization electrostrictively couples only to the deviatoric strain. This assumption will not effect the qualitative behavior of the model in context surface effects, which are the main focus of this paper. We should also keep in mind that for the present case of a square to rectangle ferroelectric transition, the deviatoric strains are the dominant mode of deformation. The bulk and shear strains are usually localized at the domain walls for this case. Including the coupling with bulk and shear deformations may introduce intermediate metastable states with more complex patterns. Exclusion of these couplings makes the computation faster due to faster equilibriation. Finally, we should also mention that the neglect of this coupling in the free energy does not meen that bulk and shear strains will not be generated in this model at all. In fact, in this model, the strains are related to each other through the elastic compatibility relations which leads to the creation of bulk and shear strains, even under the above described approximation. In the limit of no coupling of polarization with shear and bulk strains, i.e. $`q_10`$ and $`q_30`$, using the equilibrium conditions and the lagrange multiplier $`\lambda (\stackrel{}{k})`$, the free energy $`F_{eff}`$ can be expressed in Fourier space as

$$F_{eff}=\frac{q_{2}^{}{}_{}{}^{2}}{2a_2}𝑑\stackrel{}{k}H(\stackrel{}{k})|\mathrm{\Gamma }_2(\stackrel{}{k})|^2,$$
(15)
where $`H(\stackrel{}{k})=(\frac{h_{1}^{}{}_{}{}^{2}(\stackrel{}{k})}{\alpha }+h_{2}^{}{}_{}{}^{2}(\stackrel{}{k})2h_2(\stackrel{}{k})+\frac{h_{3}^{}{}_{}{}^{2}(\stackrel{}{k})}{\beta })`$. The quantities $`h_1(\stackrel{}{k})=k^2Q(\stackrel{}{k})`$, $`h_2(\stackrel{}{k})=\{1(k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2})Q(\stackrel{}{k})\}`$ and $`h_3(\stackrel{}{k})=\sqrt{8}k_xk_yQ(\stackrel{}{k})`$, with $`Q(\stackrel{}{k})`$ defined as

$$Q(\stackrel{}{k})=\frac{(k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2})}{(k^4/\alpha +(k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2})^2+8k_{x}^{}{}_{}{}^{2}k_{y}^{}{}_{}{}^{2}/\beta )}.$$
(16)
We have introduced dimensionless constants $`\alpha =a_1/a_2`$ and $`\beta =a_3/a_2`$.

The effective interaction derived above is strongly direction dependent and is crucial to describe the domain wall orientations. Similar anisotropic interactions have been considered in context of martensitic transformations .

As we wish to study heterogeneous nucleation, a dynamical formalism is needed to describe the non-equilibrium effects associated with domain switching. The time dependent Ginzburg-Landau model for the polarization fields is used. We introduce rescaled time variables $`t^{}=(t|\alpha _1|L)`$$`(\alpha _1<0)`$( $`L`$ is the kinetic coefficient ), space variable $`\stackrel{}{r^{}}=\stackrel{}{r}/\theta `$, the wave vector $`\stackrel{}{k^{}}=\stackrel{}{k}\theta `$ $`(\theta =\sqrt{g_1/a|\alpha _1|})`$, where $`a`$ is a dimensionless constant. The polarization field is transformed as $`P_x=P_Ru`$ and $`P_y=P_Rv`$, where $`P_R=\sqrt{|\alpha _1|/\alpha _{11}}`$ is the remnant polarization of the homogeneous state without elastic effects. With this set of parameters, the dimensionless time dependent Ginzburg-Landau equations are given as

$`u,_t^{}=uu^3duv^2+au,_{x^{}x^{}}+bu,_{y^{}y^{}}+cv,_{x^{}y^{}}+e_x\gamma u{\displaystyle }d\stackrel{}{k^{}}H(\stackrel{}{k^{}})\mathrm{\Gamma }(\stackrel{}{k^{}})exp(i\stackrel{}{k^{}}\stackrel{}{r^{}})`$ (17)
$`v,_t^{}=vv^3dvu^2+av,_{y^{}y^{}}+bv,_{x^{}x^{}}+cu,_{x^{}y^{}}+e_y+\gamma v{\displaystyle }d\stackrel{}{k^{}}H(\stackrel{}{k^{}})\mathrm{\Gamma }(\stackrel{}{k^{}})exp(i\stackrel{}{k^{}}\stackrel{}{r^{}}).`$ (18)
Here, $`\stackrel{}{e}=\stackrel{}{E}_{ext}/(P_R|\alpha _1|)`$ is the rescaled external electric field. The notation $`u,_t`$ represents the time derivative. The defined constants are: $`b=(g_2/|\alpha _1|\theta ^2)`$, $`c=(g_3/|\alpha _1|\theta ^2)`$, $`d=(\alpha _{12}/\alpha _{11})`$ and $`\gamma =(q_{2}^{}{}_{}{}^{2}/\alpha _{11}a_2)`$. The quantity $`\mathrm{\Gamma }(\stackrel{}{k^{}})`$ is the Fourier transform of $`(u^2v^2)`$.

## III Simulation of polarization switching for a constrained square shaped ferroelectric

We now describe the details of our simulations on the switching behavior for a finite constrained ferroelectric. The long-range interaction is conveniently defined in Fourier space, which allows us to use the Fourier Transform method. This enforces the need to have periodic boundary conditions at the surfaces. To define a finite system using this scheme, we consider an $`88\times 88`$ square region within a full simulation grid of size $`128\times 128`$. The constraining is implemented by assuming that the polarization variables $`u(\stackrel{}{r^{}})`$ and $`v(\stackrel{}{r^{}})`$ vanish at the boundary and outside of the $`88\times 88`$ region. This configuration represents an electromechanically clamped ferroelectric in which all displacements vanish at the surface. Equation (15) is discretized using the Euler scheme on the above described grid. The space discretization step $`\mathrm{\Delta }x^{}=\mathrm{\Delta }y^{}=1`$ and the time interval $`\mathrm{\Delta }t^{}=0.02`$. The free energy parameters chosed by us are $`d=2`$ and $`\gamma =0.05`$. The gradient coefficients are chosen as $`a=b=c=1`$. Similarly, we choose the elastic parameter ratios as $`\alpha =\beta =1`$. We first start from a random small amplitude initial conditions for the fields $`u(\stackrel{}{r^{}})`$ and $`v(\stackrel{}{r^{}})`$ and an external field $`e_x^{}=10`$. We solve equation (15) to obtain a single domain configuration with the polarization oriented towards the $`+x^{}`$ gradually decaying at the surface of the square system. Figure 1(a) shows this configuration corresponding to an electric field value $`e_x^{}=10`$. In this figure, the arrows represent the local polarization vector ( we have not shown all the polarization vectors for the sake of legibility). Notice the gradual decay of the amplitude of the vectors in the surface region. This configuration is used as the initial state for the simulations of polarization switching. Equations (15) are now solved at each time step, starting from the initial configuration depicted in figure 1(a) with a time dependent external field $`e_x^{}(t^{})`$. The time dependence is given as $`e_x^{}=e_0[12t^{}/T]`$, where $`e_0=10`$ is the initial electric field value and $`T=5000`$ is the total number of simulated switching time steps. As the field decreases, the magnitude of the polarization vectors decreases uniformly every where in the simulated system at the beginning (keep in mind that the polarization is fixed to zero at the surface for all times). As the external electric field approaches zero, we observe that near two of the corners, the polarization vectors are distorted towards the  directions. Figure 1(b) shows the pattern corresponding to $`e_x^{}=0.22`$ where we can see the diagonal orientation of the polarization vectors near the corners. This indicates that shear strains are generated near the corners. The result is in agreement with recent simulations by Jacobs for a constrained ferroelastic system undergoing a square to rectangle transition where it was demonstrated that nondeviatoric strains are generated at the corners . These diagonal shear strains at the corners are responsible for nucleation of $`90^o`$ domain walls, as seen figure 1(c)$`(e_x^{}=0.199)`$ and figure 1(d)$`(e_x^{}=0.0)`$. In figure 1(e) $`(e_x^{}=0.199)`$, we can clearly see a well developed twinned structure with $`90^o`$ domain walls. As we further decrease the field to $`e_x^{}=0.38`$, domains polarized along $`[01]`$ direction grow at the expense of domains polarized along $`[10]`$ directions, as can be seen by comparing figure 1(e) and figure 1(f). In figure 1(f) $`(e_x^{}=0.38)`$, we can also observe that the dipoles in the vicinity of the corners and domain walls have started reorienting towards the $`x^{}`$ direction. Thereafter, the reversed domains grow at the expense of domains polarized along $`[01]`$, and eventually all the dipoles are oriented along the $`x^{}`$ direction. This growth can be clearly seen from figure 1(g) $`(e_x^{}=0.46)`$ and figure 1(h) $`(e_x^{}=10)`$. Figure 1 illustrates the special role played by the corners during switching. The corners act as nucleation centres for $`90^o`$ domain structures. Unlike the classically held belief of $`180^o`$ reorientation of polarization, the switching takes place by sucessive $`90^o`$ reorientations of the dipoles. Recently, evidence for such non $`180^o`$ switching has been observed in rhombohedral single crystals of PZN-PT . The reverse cycle can be simulated by using the configuration in figure 1(h) as the initial condition and increasing the field as $`e_x^{}=e_0[12t^{}/T]`$, where $`e_0=10`$ and $`T=5000`$ is the total number of switching time steps. Figure 2 shows the hysteresis loops obtained by plotting the average polarization $`<u>`$ with the field $`e_x^{}`$. The solid line corresponds to the situation depicted in figure 1 and the dashed line corresponds to this loop for the case without any surface (periodic boundary conditions). We can observe that for the case of constrained square, the simulated hysteresis loop shows a long ”tail” region where the forward and reverse cycles curves do not close up. This corresponds to the appearence and disappearence of the $`90^o`$ domain pattern. In contrast, for the periodic boundary condition case, hysteresis loop has no ”tail” region as there are no nucleation events. We can also observe that the saturation polarization and remnant polarization field are lower for the constrained system than that for the case with periodic boundary conditions. However, there is not much difference between the coercive fields of the two cases. We have simulated the hysteresis loops for different system sizes to study the size dependence of the switching phenomena. In figure 3, we show the simulated hysteresis loops for system sizes indicated at the top of each frame. A comparison of the loops for sizes $`L=88`$ and $`L=68`$ reveals that the coercive fields and the saturation polarization are almost the same for both systems. Also, the $`90^o`$ domain structure is shorter lived for the smaller system since the ”tail” region of the hysteresis loop for the $`L=68`$ system is smaller than that for the $`L=88`$ system. For smaller sizes, the domain pattern becomes increasingly shortlived during the switching process, as can be concluded from the loop for $`L=48`$ where the tail has moved further inwards. However, the saturation polarization and the coercive fields still do not show much variation. The tail almost disappears for the $`L=20`$ hysteresis loop and we see a small bump near the coercive field region. For this size, we can see a reduction in the coercive field as well as the saturation polarization due to the surface effects. Below this size, nucleation of $`90^o`$ domains does not take place. The switching occurs by $`180^o`$ reorientations that preferentially starts at the surface. The hysteresis loops in this regime have shape similar to the periodic boundary condition case. It is clear from figure 3 that on further decreasing the size, the saturation polarization, remnant polarization and coercive field decrease rapidly. We can see that the hysteresis is very small for $`L=6`$ which is the smallest size simulated. For the present set of parameter values, we do not observe a ferroelectric to paraelectric transition for $`L=6`$. Such a transition was observed for a simulation with much higher value of the gradient coefficients. The behavior observed in figure 3 can be understood in terms of size dependence of ferroelectric behavior. It is experimentally and theoretically established that the spontaneous polarization progressively decreases on decreasing size for a system in which the polarization is constrained to decay at the surface . Thus, the remnant polarization, saturation polarization and coercive field decrease on decreasing the system size. The anisotropic nonlocal interaction responsible for the nucleation of $`90^o`$ domain structures also decreases as the size is decreased. Thus for small sizes, the hysteresis loops are free of $`90^o`$ nucleation. A similar transition from multidomain to single domain states on decreasing the size has been discussed in earlier works . The main effect of size reduction comes from the surface layer under the vanishing polarization boundary condition. The transition temperature is effectively reduced and the energy barrier between the low temperatures states is also lowered.

## IV Simulation of polarization switching for a semi infinite constrained ferroelectric

In this section, we describe switching behavior for a semi-infinite constrained ferroelectric system. Here we use periodic boundary conditions along the $`y`$ direction and vanishing polarization conditions at the surface for $`x`$ direction. Figure 4 shows the dynamics of polarization switching for the semi-infinite case. Except the boundary conditions, all other parameter values, including the time dependence of the switching field are the same as that for the case of the constrained square described in section 3. Figure 4(a) shows the fully polarized single domain initial configuration for $`e_x^{}=10`$. In figure 4(b), we can see that a single domain persists even when the applied field is lowered to $`e_x=0.319`$, although the polarization magnitude has decreased by a very small amount. As the field is further decreased to $`e_x^{}0.36`$, nucleation of the reversed polarization domains from the two surfaces takes place, as depicted in figure 4(c). The reversed domains grow larger and eventually into a single domain crystal with reversed polarization. This reversed domain growth is clear from figure 4(d) ($`e_x^{}=0.399`$) and figure 4(e) ($`e_x^{}=0.44`$). Figure 4(f) shows the final state corresponding to $`e_x^{}=10`$. Figure 4 illustrates the fact that for the corner free semi-infinite geometry, the reversal occurs only by $`180^o`$ reorientations. The hysteresis loop corresponding to the situation shown in figure 4 can be computed. In figure 5, we plot the hysteresis loop for the case shown in figure 4 along with the loop for the surface free case with the periodic boundary conditions. Unlike the case shown in figure 2, the hysteresis loops for the surface free case and the semi-infinite case have the same shape. However, the saturation polarization, the remnant polarization and the coercive field are slightly smaller for the semi-infinite case due to the suppressed transition at the surface. We have also simulated the size dependence for this case and found that the hysteresis loops become progressively narrower with decreasing size without any shape change. This is due to the fact that there are no $`90^o`$ domains and the long-range elastic interaction does not play any role in this case.

## V Summary and discussion

We have simulated the switching behavior in constrained ferroelectrics with the boundary condition of vanishing polarization at the surface. A square shaped system is considered along with a system with semi-infinite geometry. We find that the corners of the square shaped system nucleate $`90^o`$ domains and the switching is dominated by the motion of $`90^o`$ domain walls. In contrast, for the semi-infinite geometry, there is no nucleation of $`90^o`$ domains and the switching proceeds by $`180^o`$ reversal. Clearly, the corners in the constrained square play a crucial role during the switching process. We believe that the anisotropy of the long-range elastic interaction interacts with the corners to produce inhomogenieties that result in $`90^o`$ nucleation. Our results suggest that the switching behavior and the hysteresis loops in real systems may depend on the shape of the grain and its boundary conditions.

We have also simulated the size dependence of hysteresis loops. The decrease of the saturation polarization, remnant polarization and the coercive field with decreasing size is common to both configurations. This can be understood in terms of the supression of the transition in an electromechanically constrained system as the system size approaches the correlation length of order parameter. The free energy for a constrained system is a function of the system size. For the constrained square case, the polarization reversal process is dominated by intermediate $`90^o`$ domain pattern changes for larger systems and becomes totally via $`180^o`$ reorientations for size less than a critical dimension. This transition is due to the weakening of the long-range elastic interactions with decreasing system size.

Finally, one should note that the behavior of decreasing coercive field with decreasing size indicated by the simulation results does not apply to thin films , but only to clamped grains in ferroelectric polycrystals. For thin films, not all the dimensions are reduced simultaneously, which is different than our condition here. In addition, the misfit strain with the substrate is important, which has to be incorporated in the theory .

## Acknowledgement

This research is sponsored by the Hong Kong Research Grant Council under grant number CityU 1056/02P.

Figure captions:

Figure. 1: Pattern evolution for the simulated switching in a constrained square. The snapshots correspond to (a) $`e_x^{}=10`$; (b) $`e_x^{}=0.22`$; (c) $`e_x^{}=0.199`$; (d) $`e_x^{}=0`$; (e) $`e_x^{}=0.199`$; (f) $`e_x^{}=0.38`$; (g) $`e_x^{}=0.46`$; (h) $`e_x^{}=10`$.

Figure. 2: Plots of $`<u>`$ vs $`e_x`$ showing the simulated hysteresis loops for the constrained square (solid line) and the periodic boundary condition case (dotted line).

Figure. 3: Simulated hysteresis loops for different sizes of the constrained square.

Figure. 4: Pattern evolution for the simulated switching in a semi-infinite constrained system. The snapshots correspond to (a) $`e_x^{}=10`$; (b) $`e_x^{}=0.319`$; (c) $`e_x^{}=0.36`$; (d) $`e_x^{}=0.399`$; (e) $`e_x^{}=0.44`$; (f) $`e_x^{}=10`$.

Figure. 5: Plots of $`<u>`$ vs $`e_x`$ showing the simulated hysteresis loops for the semi-infinite constrained system (solid line) and the periodic boundary condition case (dotted line).