text:
# DO–TH 00/12 hep-ph/0008044 Cancellations in Neutrinoless Double Beta Decay and the Neutrino Mass Matrix

## 1 Introduction

In the light of recent impressive experimental evidence on neutrino oscillations  the next fundamental question to be answered is the one regarding the neutrino character. From the theoretical side Majorana neutrinos are favored, since they pop out of almost every GUT and are for example the result of the very attractive see–saw–mechanism . In this case, heavy neutrinos are predicted. Experimental information on the existence of Majorana neutrinos might come from neutrinoless double beta decay (0$`\nu \beta \beta `$)  or from production of heavy Majorana neutrinos at colliders, see e.g.  and references therein. The most stringent limit however comes from the first process, investigating the effective Majorana mass of the electron neutrino, with a current bound of

$$m<0.2\mathrm{eV}.$$
(1)
Plans exist to build experiments exploring regions up to $`m0.002`$ eV . The variation within a factor of roughly 3 between different calculations of the required nuclear matrix elements has to be kept in mind. Results of oscillation experiments can be used to restrict the value of $`m`$ in different mass schemes and for the different solutions of the solar neutrino problem . Typical key scales for $`m`$ are 0.1 and 0.005 eV, thus lying in the range of current and forthcoming projects. However, $`m`$ has a form which includes possibilities for the contributions to cancel each other, namely via $`CP`$–violating phases or via the intrinsic $`CP`$–parities of the mass states, which exist in the case of $`CP`$ conservation. These effects were included in most of the above given references to get the maximal and minimal values of $`m`$. In  constraints on $`CP`$–violating phases were given by using a graphical representation for the complex mass matrix elements, although without applying numerical values. Numerical studies regarding the Majorana phases were given in . Some overlap with these works exists, however, we give many plots and statements, which are new to the literature and adopt an approach, which allows to investigate different situations regarding measurements of $`m`$ and/or $`m_0`$. Here $`m_0`$ denotes the common mass scale in a degenerate mass scheme. Among the topics discussed are a clarification of the kind of cancellation, i.e. a distinction between the 2 possibilities. When starting from a degenerate scheme and in a given situation the maximal allowed $`m_0`$ is comparable to the scale implied by the atmospheric anomaly, so the degenerate scheme is ruled out and the derived mass limit holds for the largest of the 3 mass eigenstates. We will find that only in very special cases definite statements about the phases or parities can be made and that alone from this fact there is a large uncertainty in the values of the other mass matrix entries. This uncertainty is ranging from 2 to factors around 20. The other entries are clearly needed to distinguish between different discussed models. The only information about the phases or parities can come from neutrinoless double beta decay: All other mass matrix entries are impossible to measure directly, since the respective branching ratios or cross sections they govern are way beyond experimental access . Apart from the obvious aspect that fundamental parameters of a model need to be measured, other interesting application of the phases exist: In  it was found that the Majorana phases play a crucial role in the stability of the mixing angles against quantum corrections. Ref.  finds that their values have influence on the magnitude of the lepton–asymmetry in the universe, which can be made responsible for the observed baryon–asymmetry. Thus, the precise knowledge of all mixing parameters in the lepton sector is certainly important.
The paper is organized as follows: In Section 2 we present the general framework and basic formulae for $`CP`$ violation and conservation, respectively. Section 3 sees a discussion of the connection between oscillation and $`m`$ in hierarchical schemes, whereas Sec. 4 concentrates on the degenerate scheme. Some special mixing matrices and the general case with experimentally favored values are discussed. For each case constraints on the phases, parities and $`m_0`$ are given. The paper is concluded in Sec. 5 with a summary of our results.

## 2 Formalism

Flavor eigenstates $`\nu _\alpha `$ ($`\alpha =e,\mu ,\tau `$) are connected to mass eigenstates $`\nu _i`$ ($`i=1,2,3`$) via a mixing matrix, i.e. $`\nu _\alpha =S_{\alpha i}\nu _i`$. A proper treatment of this issue can be found in , here we quote only the main points. This matrix is derived by diagonalising the Majorana mass term in the Langrangian:

$$_{\mathrm{mass}}=\frac{1}{2}M_{\alpha \beta }\overline{\nu _{\alpha L}}\text{ }\nu _{\beta R}^c\text{ }+\mathrm{h}.\mathrm{c}.$$
(2)
where $`U_{CP}\nu _{\alpha L}\text{ }U_{CP}^1=\nu _{\alpha R}^c`$ and $`M_{\alpha \beta }`$ is a symmetric matrix. For $`CP`$ invariance it is also real and can be diagonalized by an orthogonal matrix

$$O_{i\alpha }M_{\alpha \beta }O_{\beta j}^T=\eta _im_i\delta _{ij}$$
(3)
with $`\eta _i=\pm 1`$ and $`m_i0`$. Choosing $`\chi _i=_aO_{i\alpha }\nu _{\alpha L}\text{ }+\eta _iO_{i\alpha }(\nu _{\alpha L})^c\text{ }`$ we have $`2_{\mathrm{mass}}=m_i\overline{\chi }_i\chi _i`$ and $`U_{CP}\chi _iU_{CP}^1=\eta _i\chi _i`$. For $`CP`$ invariance it can be shown  that $`\eta _i`$ is connected to the intrinsic $`CP`$ parity of the Majorana, i.e. $`\eta _i=i\eta _i^{CP}`$. In addition it holds $`\chi ^c=\pm \eta _i\chi `$. Therefore the 0$`\nu \beta \beta `$ amplitude for $`CP`$ invariance is proportional to

$$_{CP}(0\nu \beta \beta )O_{ei}^2\gamma _{}\chi _i\overline{\chi _i^c}\gamma _{}\eta _im_iO_{ei}^2,$$
(4)
with $`2\gamma _\pm =1\pm \gamma _5`$. At the cost of a complex $`O`$ the $`CP`$ parities can also be absorbed in the mixing matrix via the identity $`\eta _i=e^{i\pi /2(\eta _i1)}`$. For complex $`M_{\alpha \beta }`$ and thus $`CP`$ violation222A geometrical description of $`CP`$ violation with Majorana neutrinos in terms of unitarity triangles can be found in . the diagonalization of the mass term is done by an unitary matrix,

$$U_{i\alpha }M_{\alpha \beta }U_{\beta j}^T=m_i\delta _{ij}.$$
(5)
Appropriate choice of $`U`$ can always make $`m_i0`$. Here, $`\eta _i=1`$ and the 0$`\nu \beta \beta `$ amplitude is proportional to

$$_{CP\text{/}}(0\nu \beta \beta )U_{ei}^2\gamma _{}\chi _i\overline{\chi _i^c}\gamma _{}m_iU_{ei}^2.$$
(6)
To conclude, the matrix in the relation $`\nu _\alpha =S_{\alpha i}\nu _i`$ is either orthogonal or unitarity and the quantity probed in 0$`\nu \beta \beta `$ is

$$m=\{\begin{array}{cc}|\underset{i}{}U_{ei}^2m_i|& CP\text{ violation}\\ |\underset{i}{}O_{ei}^2m_i\eta _i|& CP\text{ invariance}\end{array}.$$
(7)
For 2 flavors the two cases read:

$$m=\{\begin{array}{cc}|m_1\mathrm{cos}^2\theta +m_2\mathrm{sin}^2\theta e^{2i\varphi }|& CP\text{ violation}\\ |m_1\mathrm{cos}^2\theta +m_2\eta _1\eta _2\mathrm{sin}^2\theta |& CP\text{ invariance}\end{array}.$$
(8)
Hence, a single $`CP`$ parity has no physical meaning, only relative values are significant. Note though that the case $`\eta _1\eta _2=1`$ (opposite parities) and $`\varphi =\pi /2`$ (“maximal” violation) can not be distinguished. This can be also seen from the choice of the mixing matrix $`S_{\alpha i}`$: The case $`\eta _1=1`$ and $`S_{e1}=O_{e1}`$ is equivalent to $`\eta _1=+1`$ and $`S_{e1}=e^{i\pi /2}O_{e1}=iO_{e1}`$. This means that opposite parities with a real mixing matrix are equivalent to equal parities with a complex mixing matrix for a maximal phase. Hence, for “maximal” phases we can not tell from 0$`\nu \beta \beta `$ alone if there is $`CP`$ violation in the lepton sector and the answer if there is one at all has to come from long–baseline experiments . This is similar to the statement by Pal and Wolfenstein  for the decay $`\nu _2\nu _1\gamma `$: For equal (opposite) parities magnetic (electric) dipole radiation occurs. Opposite parities are equivalent to complex mixing matrix elements but do not imply $`CP`$ violation. This would be suggested only by the presence of both kinds of radiation. For a degenerate scheme with masses discussed in the following sections however, the life time of such a process is about $`\tau >10^{48}`$ s.

## 3 Oscillation Experiments, 0$`\nu \beta \beta `$ and $`m`$ in Hierarchical Schemes

The $`CP`$–violating mixing matrix $`U`$ can be parametrized as

$$\begin{array}{c}U=U_{\mathrm{CKM}}\mathrm{diag}(1,e^{i\alpha },e^{i(\beta +\delta )})\\ =\left(\begin{array}{ccc}c_1c_3& s_1c_3& s_3e^{i\delta }\\ s_1c_2c_1s_2s_3e^{i\delta }& c_1c_2s_1s_2s_3e^{i\delta }& s_2c_3\\ s_1s_2c_1c_2s_3e^{i\delta }& c_1s_2s_1c_2s_3e^{i\delta }& c_2c_3\end{array}\right)\mathrm{diag}(1,e^{i\alpha },e^{i(\beta +\delta )}),\end{array}$$
(9)
where $`c_i=\mathrm{cos}\theta _i`$ and $`s_i=\mathrm{sin}\theta _i`$. The orthogonal matrix $`O`$ is of course obtained by setting the phases to zero. From hereon we will always write $`U_{\alpha i}`$ for the mixing matrix with obvious changes for the $`CP`$ conserving case. Since the oscillation probability

$$P_{\alpha \beta }=\delta _{\alpha \beta }2\text{ Re }\underset{j>i}{}U_{\alpha i}U_{\alpha j}^{}U_{\beta i}^{}U_{\beta j}(1\mathrm{exp}i\mathrm{\Delta }_{ji}).$$
(10)
is invariant under such a multiplication of a diagonal matrix, the two additional Majorana induced phases are not observable in any oscillation experiment . Therefore, 0$`\nu \beta \beta `$ is the only probe to test these phases. In principle it would be possible to completely derive all phases by comparing the elements of the Majorana mass matrix

$$m_{\alpha \beta }=|(U\mathrm{diag}(m_1\eta _1,m_2\eta _2,m_3\eta _3)U^\mathrm{T})_{\alpha \beta }|.$$
(11)
However, the analogues of 0$`\nu \beta \beta `$ one needs to observe like $`K^+\pi ^{}\mu ^+\mu ^+`$ , $`\nu _\mu N\mu ^{}\mu ^+\mu ^+X`$  or $`e^+p\overline{\nu }_e\text{ }\mu ^+\tau ^+X`$  have branching ratios or cross section far too small to be measured .
All numerical analyses indicate a hierarchy in the mass squared differences of solar and atmospheric experiments:

$$\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_A^2.$$
(12)
Maximal values are $`(\mathrm{\Delta }m_{}^2)_{\mathrm{max}}10^4`$ eV2  and $`(\mathrm{\Delta }m_A^2)_{\mathrm{max}}10^2`$ eV2 . Without loss of generality we assume $`m_3m_2m_1`$. In a 3 flavor picture, three hierarchical schemes, where at least for one mass eigenstate $`m_i^2\mathrm{\Delta }m^2`$ holds, are capable of explaining the relation (12), they are called “completely hierarchical” ($`m_3\sqrt{\mathrm{\Delta }m_A^2}m_2\sqrt{\mathrm{\Delta }m_{}^2}m_1`$), “partially hierarchical” ($`m_3\sqrt{\mathrm{\Delta }m_A^2}m_2m_1`$) and “inverse hierarchical” ($`m_3m_2=\sqrt{\mathrm{\Delta }m_A^2}m_1`$). With this notation, the CHOOZ result demands a small $`|U_{e3}|`$ in the first two cases and a small $`|U_{e1}|`$ in the last one. The latter fact stems from the condition that for $`m_3m_2m_1`$ the completely and partially hierarchical schemes demand $`\mathrm{\Delta }m_{21}^2=\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{31}^2\mathrm{\Delta }m_{32}^2`$, whereas the inverse hierarchical scheme has to be given by $`\mathrm{\Delta }m_{32}^2=\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{31}^2`$. With a good approximation we can use 2 flavor fits, then the solar mixing angle gives

$$c_1^2,s_1^2\frac{1}{2}\left(1\pm \sqrt{1\mathrm{sin}^22\theta _{}}\right).$$
(13)
The data we use are taken from Refs.  and given in Table 1. As usual, SA (LA) denotes the small (large) angle MSW, LOW the MSW low mass and VO the vacuum solution. The MSW  resonance condition demands $`|U_{e1}||U_{e2}|`$. However, for the vacuum solution this might not be the case, we come back later on that point. In the inverse hierarchical scheme the resonance requires $`|U_{e3}||U_{e2}|`$. The CHOOZ experiment  gives unfortunately only a limit on $`|U_{e3}|`$ (in the inverse hierarchical scheme on $`|U_{e1}|`$), depending on the atmospheric mass scale, it reads at 90 $`\%`$ C. L.

$$|U_{e3}|^2<0.01\mathrm{}0.15\text{ for }\mathrm{\Delta }m_A^210^3\mathrm{}10^2\mathrm{eV}^2\text{[35]}.$$
(14)
For typical best–fit values of few$`\times 10^3`$ eV the bound is about $`|U_{e3}|^2<0.05`$. Even for the maximal allowed values of the mass squared differences and mixing angles $`m`$ is always below 0.2 eV:

$$\begin{array}{c}mm_1|U_{e1}|^2+m_2|U_{e2}|^2+m_3|U_{e3}|^2\\ \{\begin{array}{ccc}m_1+|U_{e2}|^2\sqrt{\mathrm{\Delta }m_{}^2}+|U_{e3}|^2\sqrt{\mathrm{\Delta }m_A^2}& m_1+0.022& \mathrm{completely}\mathrm{hierarchical}\\ (|U_{e1}|^2+|U_{e2}|^2)\sqrt{\mathrm{\Delta }m_{}^2}+|U_{e3}|^2\sqrt{\mathrm{\Delta }m_A^2}& 0.85m_1+0.015& \mathrm{partially}\mathrm{hierarchical}\\ |U_{e1}|^2m_1+(1|U_{e1}|^2)\sqrt{\mathrm{\Delta }m_A^2}& m_1+0.085& \mathrm{inverse}\mathrm{hierarchical}\end{array}.\end{array}$$
(15)
Note that there can be a additional contribution of $`m_1<\sqrt{\mathrm{\Delta }m_{}^2}`$, which can be traced back to the fact that only mass squared differences are measured. If an experimental sensitivity of $`10^2`$ eV on $`m`$ is achieved, this might be of importance. However, at the present time with $`m0.2`$ eV, all values of the phases and the parities are allowed. In order to make more definite statements about the effects of cancellations, we therefore apply a degenerate scheme in which $`m_3^2m_2^2m_1^2m_0^2`$. Per definition, it should hold $`m_0^2\mathrm{\Delta }m_{\mathrm{max}}^210^2`$. When in the following a maximal allowed mass scale of $`m_00.2\mathrm{}0.3`$ eV is derived so does that mean that the mass scheme is not “completely degenerate” but “slightly” hierarchical, i.e. an intermediate situation occurs. For example, if $`m_{0,\mathrm{max}}=0.2`$ eV, then another eigenstate has a mass of $`(0.2\epsilon )`$ eV, where $`\epsilon `$ is small. In order to get the atmospheric scale of $`10^2`$ eV2 one finds $`\epsilon 0.03`$ eV and the two eigenstates differ by about 15 $`\%`$. Then our $`m_{0,\mathrm{max}}`$ is the limit on the largest mass eigenstate $`m_3`$ with $`m_3\text{ or }>m_2`$. A limit of $`m_00.1`$ eV corresponds to a hierarchical scheme. However, due to the uncertainty in the value of $`m`$ and $`\mathrm{\Delta }m_A^2`$ a definite statement is somehow difficult to make. Using the best–fit point for the atmospheric scale means that the degenerate scheme fails for about 0.1 eV.

## 4 $`m`$ in the Degenerate Scheme

The upper bound on $`m_0`$ comes from the tritium spectrum, which limits

$$\underset{i}{}|U_{ei}^2|m_i^2=m_0^2(2.8\mathrm{eV})^2.$$
(16)
Forthcoming and ongoing projects intend to push the bound below 1 eV . In addition, cosmological observations might be interpreted in terms of a total neutrino mass of $`m_\nu `$ few eV . By measuring anisotropies in the cosmic microwave background, MAP and PLANCK may probe values down to $`m_\nu 0.5`$ eV . The interesting function to investigate is $`m/m_0\stackrel{~}{m}`$, which is depending on 4 parameters, either two angles and 2 phases or two angles and two relative parities. In the $`CP`$–violating case $`\stackrel{~}{m}`$ reads

$$\frac{m}{m_0}\stackrel{~}{m}=c_3^2\sqrt{c_1^4+s_1^4+t_3^4+2(s_1^2t_3^2c_{2(\alpha \beta )}+c_1^2(s_1^2c_{2\alpha }+t_3^2c_{2\beta }))}$$
(17)
where $`t_3=\mathrm{tan}\theta _3`$ and $`c_{2\alpha }=\mathrm{cos}2\alpha `$. For $`CP`$ invariance it can be written as

$$\stackrel{~}{m}=c_3^2(\eta _1c_1^2+\eta _2s_1^2+\eta _3t_3^2)=\{\begin{array}{ccc}1& (+++)& \alpha =\beta =\pi \\ |c_3^2(c_1^2s_1^2t_3^2)|& (+)& \alpha =\beta =\frac{\pi }{2}\\ |c_3^2(c_1^2s_1^2+t_3^2)|& (+)& \alpha =\frac{\beta }{2}=\frac{\pi }{2}\\ \mathrm{cos}2\theta _3& (+)& \alpha =2\beta =\pi \end{array},$$
(18)
where all 4 possible $`(\pm \pm \pm )`$ signatures with the corresponding $`CP`$–violating phases are given. In addition, there are ($`c_1`$–dependent) solutions for the phases, which give the same $`\stackrel{~}{m}`$ as special parity configurations, see below. Note that for $`(+)`$ the value is independent on the solar solution. The $`(+)`$ signature is in fact the minimal value for the general $`CP`$–violating case, provided that $`|U_{e1}||U_{e2}|`$. For the inverse situation, the $`(+)`$ signature gives the minimal $`\stackrel{~}{m}`$.
The different treatment of $`CP`$ parities and phases may seem somewhat artificial since the first option is a special case of the latter one. However, due to the maximal values certain configurations result in they deserve special attention.

### 4.1 Some special mixing matrices

Three intriguingly simple matrices have been widely discussed:

* Single maximal
Maximal atmospheric mixing and a vanishing angle in solar and the CHOOZ experiment, resulting in $`|U_{e1}|=1`$.
* Bimaximal
Both solar and atmospheric mixing is maximal, CHOOZ’s angle is zero. Then $`|U_{e1}|=|U_{e2}|=1/\sqrt{2}`$ and $`|U_{e3}|=0`$.
* Trimaximal
All elements have the same magnitude and $`|U_{e1}|=|U_{e2}|=|U_{e3}|=1/\sqrt{3}`$. The model gives a poor fit to the oscillation data.

Then one gets for the $`CP`$–conserving case

$$\stackrel{~}{m}=\{\begin{array}{cc}1,1,1& (+++)\\ 1,0,\frac{1}{3}& (+)\\ 1,0,\frac{1}{3}& (+)\\ 1,1,\frac{1}{3}& (+)\end{array},$$
(19)
for single, bi– and trimaximal mixing respectively. Single and bimaximal mixing are special cases of the mixing angles from Table 1 and might serve as a model to get a feeling for the bounds one obtains. With the current limit on $`m`$ we get for values of $`\stackrel{~}{m}`$ = $`1,\frac{1}{3}`$: $`m_00.2,0.6`$. Obviously, for $`\stackrel{~}{m}=0`$ it follows $`m=0`$. If there is $`CP`$ violation we have

$$\stackrel{~}{m}=\{\begin{array}{cc}1& \text{ single}\\ c_\alpha & \text{ bi}\\ \frac{1}{3}\sqrt{3+2(c_{2\alpha }+c_{2\beta }+c_{2(\alpha \beta )})}& \text{ tri}\end{array}.$$
(20)
Single maximal mixing means $`m_00.2`$ eV, unobtainable by currently planned experiments. For bimaximal mixing we get with the current limit and assuming $`m_0=1`$ eV: $`\alpha <78.5^0`$. For the trimaximal case Fig. 2 shows $`\stackrel{~}{m}`$ as a function of one phase for different values of the second phase. Note the constant value for $`\beta =\pi /2`$.

### 4.2 General Limits and Bounds

We go now from the special cases back to all experimentally allowed values of the mixing angels. The first thing to say is that for exactly vanishing $`U_{e3}`$ there is no way to find out about the second phase333Also, “normal” $`CP`$ violation in long–baseline experiments will be unobservable.. Without further effort we can say for the $`(+++)`$ signature, that $`m_00.2`$, i.e. is beyond experimental access in cosmology or spectrum measurements. More interesting is e.g. the $`(+)`$ case in which $`\stackrel{~}{m}`$ is depending on only one parameter, namely the angle bounded by CHOOZ. With Eqs. (1) and (14) we get the allowed range of $`m_0`$ depicted in Fig. 2. In this situation, the maximal neutrino mass is about 0.29 eV, using the best–fit point of SK gives $`m_00.22`$ eV. We also plot the range for a limit of $`m<0.1`$ eV, which further reduces the allowed values. The maximal $`m_0`$ is now 0.14 and 0.11 eV, respectively. As can be seen, if cosmology insists in $`m_01`$ eV, the $`(+)`$ configuration is ruled out.
There is still freedom in the ordering of $`|U_{e1}|`$ and $`|U_{e2}|`$ if the vacuum solution is correct. We denote the choice $`|U_{e1}|>|U_{e2}|`$ with VO1 and the other one with VO2. However, as can be seen from Eqs. (13) and (18) the case VO1 and ($`+`$) is equivalent to VO2 and $`(+)`$ as VO1 and ($`+`$) is to VO2 and ($`+`$). In Figures 4 and 4 we show the maximal values of $`m_0`$ as a function of $`|U_{e3}|^2`$ for the current $`m`$–limit of 0.2 eV. These maximal values scale with the $`m`$ limit. In  plots of the allowed $`|U_{e1}|^2`$ and $`|U_{e3}|^2`$ were given for $`m_0=1.7`$ eV and all 3 nontrivial parity configurations. Their conclusion that small $`|U_{e3}|`$ requires near–maximal mixing for the $`(+)`$ and $`(+)`$ case is consistent with Figs. 4 and 4. For $`(+)`$ they find that only small $`m_0`$ allows large mixing, which can be seen in Fig. 2 as well. The authors also give the allowed areas in $`|U_{e1}|^2`$$`|U_{e3}|^2`$ space for special values of the phases. However, to obtain these areas they included all $`(\pm \pm \pm )`$ possibilities, which is hard to compare with our approach. Reference  also plots $`|U_{e1}|^2`$ against $`|U_{e3}|^2`$ for different situations resulting in similar conclusions as ours. Our plots are thus different projections of the 5 dimensional parameter space giving complementary and additional information. Both mentioned works use cosmological arguments to set $`m_0`$ to 1.7 to 4 eV, whereas our approach (also used in ) allows to investigate different situations like positive results on $`m`$ and/or $`m_0`$.
In Figs. 6 and 6 we plot $`\stackrel{~}{m}`$ as a function of $`\alpha `$ for different values of $`\beta `$ whilst assuming $`|U_{e3}|^2=0.03`$. Maximal mixing (i.e. $`|U_{e1}|=|U_{e2}|`$) allows complete cancellation, whereas for other values a non–vanishing minimal $`\stackrel{~}{m}`$ is achieved. The dependence on the second phase is rather small which can be explained by the smallness of $`|U_{e3}|`$, confer with Fig. 2, where the dependence on the second phase is rather strong. In Fig. 8 we therefore show for $`m0.2`$ eV, $`\beta =\pi /2`$ and $`|U_{e3}|^2=0.03`$ the allowed area in the $`m_0`$$`\alpha `$ space for different $`\mathrm{sin}^22\theta _{}`$. The higher this angle is, the higher is the maximal allowed $`m_0`$. This value is approached for a solar mixing angle of about 0.7, which is also true for the $`(+)`$ and $`(+)`$ configurations. For the SA solution it follows $`c_11`$ in Eq. (18) so that $`\stackrel{~}{m}`$ is always about 1. This means if the SA solution turns out to be realized in nature direct mass searches will not find the neutrino mass scale because eutrinoless double beta decay then already demands $`m_00.2`$ eV. If $`m_0>0.5`$ eV is measured and SA turns out to be true, then neutrinos are Dirac particles. These are well known facts and are only given for the sake of completeness.
The best–fit value and the minimal respectively maximal allowed mixing angle gives lines too close together to be distinguishable in the plot. Regarding VO1 and VO2 it turns out that the difference for given phases is negligible (1 to 2 $`\%`$), so we plot in Fig. 6 only the VO1 option. The difference between the two options would only be sizable for large $`|U_{e3}|`$.
If the bounds on $`m`$ and $`m_0`$ are both further reduced, nothing special happens. Some more interesting possibilities however suggest themselves: First, a value for $`m_0`$ is found and $`m`$ continues to give an upper bound. Then we get a lower limit on $`\stackrel{~}{m}`$. Second, the reverse situation results in an upper limit on $`\stackrel{~}{m}`$. If both quantities are measured, then definite statements about the phases can be made. In connection with the then already obtained knowledge of the correct solar solution and a more stringent CHOOZ bound (or signal) some important conclusions could be drawn. In addition, we can assume that the precision in the mixing angles will be improved. To simplify things a bit we take now the best–fit values. Since these points of the VO and the LA solution are very close together we only plot the LA and the LOW case in Figs. 8 and 10, respectively. Recent data of the day–night spectrum seems to favor these solutions . Different values of the parameters lead to modifications of the results which are not difficult to do with the help of Eqs. (13,17,18). From the figures, we can easily analyze some specific configurations: For example, for $`\stackrel{~}{m}<0.6`$ the $`(+)`$ configuration is ruled out and for $`\stackrel{~}{m}>0.6`$ the $`(+)`$ and $`(+)`$ cases, making 0.6 a key scale, since it is independent on the solar solution and the value of $`|U_{e3}|`$. If we would know the value of $`\stackrel{~}{m}`$, even more could be said: The key value of 0.6 would mean that there is $`CP`$ violation in the leptonic sector. With the current limits on $`m`$and $`m_0`$ and the maximal sensitivities achievable by current experiments, only $`\stackrel{~}{m}<0.6`$ seems realistic. The highest value with currently planned experiments is in fact $`\stackrel{~}{m}=0.2/0.5=0.4`$, which for small $`|U_{e3}|`$ would rule out $`(+)`$ and $`(+)`$ for LOW. Also, again for small $`|U_{e3}|`$, the LA and VO solution would be out of the game. Note however that the mentioned uncertainty for the nuclear matrix elements might allow to use a higher value for $`m`$ , namely about $`\sqrt{3}`$$`m`$ which then might allow to reach a $`\stackrel{~}{m}`$ value of 0.6.
Conversely, the smallest number one can expect to measure is $`\stackrel{~}{m}=210^3/2.8=710^4`$, which forbids the degenerate scheme at all and demands hierarchical or intermediate schemes. We refer to Refs.  for the allowed values of $`m`$ in that case. In general, as can be seen from the Figures, a value of $`\stackrel{~}{m}`$ below 0.24 (0.06) rules out the degenerate scheme for the LA and VO (LOW) solution. Alternatively, if LA (LOW) turns out to be true and a $`\stackrel{~}{m}`$ bound of 0.24 (0.06) is achieved, neutrinos are Dirac particles. Note however that our discussion relies on precise knowledge of the mixing angles and $`\stackrel{~}{m}`$ and has thus to be taken with care. A realistic assumption about the minimal uncertainty in reach for these measurements might be given by 20 $`\%`$ for the former and a factor of 2 for the latter. We will show in the next section that the uncertainty for the other mass matrix elements is at least of this order.

### 4.3 Consequences for the Mass Matrix

We want to show that the uncertainty in the other mass matrix entries is at least as large as the one stemming from $`m_0`$, nuclear matrix elements and the solar mixing angle. To show this we shall assume that both quantities as well as $`|U_{e3}|`$ and the angle governing the atmospheric neutrino anomaly are known precisely. As an example we take $`\stackrel{~}{m}=0.3`$ and $`|U_{e3}|^2=0.1`$. This might be measured by $`m=0.18`$ eV and $`m_0=0.6`$ eV. However, different choice of the 2 phases can result in the same $`\stackrel{~}{m}`$, as can be seen in Fig. 10 for the LOW case, where “iso–$`\stackrel{~}{m}`$” lines are displayed. For example, $`(\alpha 1.177,\beta 1.864)`$ or $`(\alpha 1.7691,\beta 0.5225)`$ are possible solutions for $`\stackrel{~}{m}=0.3`$. The $`\mu \mu `$ entry of the mass matrix, $`m_{\mu \mu }=m_0|(U_{\mu 1}^2+U_{\mu 2}^2+U_{\mu 3}^2)|`$ is still a function of the third phase $`\delta `$. Assuming maximal atmospheric mixing, the two results for $`m_{\mu \mu }`$ differ by a factor of 4 (for $`\delta =0`$) to 8 for ($`\delta 2.8`$). It turns out that for a given $`\delta `$ some specific choices of the other two phases give for the resulting $`m_{\mu \mu }`$ relative differences of up to 15. Thus, to verify the real solution, one would need a measurement of $`m_{\mu \mu }`$. The charged kaon decay $`K^+\pi ^{}\mu ^+\mu ^+`$, which is depending on $`m_{\mu \mu }`$ as 0$`\nu \beta \beta `$ is on $`m`$, has a branching ratio given by

$$\frac{\mathrm{\Gamma }(K^+\pi ^{}\mu ^+\mu ^+)}{\mathrm{\Gamma }(K^+\mathrm{all})}10^{31}\left(\frac{m_{\mu \mu }}{\mathrm{eV}}\right)^2<810^{31},$$
(21)
which has to be compared with the experimental limit of $`3\times 10^9`$ . Even worse numbers hold for the other elements of the mass matrices and the respective processes they govern, see  for a compilation. The closer the measured $`\stackrel{~}{m}`$ value gets to the minimal value (0.1205 for our specific example) the smaller gets the allowed area or curve in $`\alpha `$$`\beta `$ space and therefore the resulting range for the other matrix elements. This minimal value of $`\stackrel{~}{m}`$ corresponds to the $`(+)`$ configuration. Of course, also for the maximal value of $`\stackrel{~}{m}=1`$ only one pair of phases is responsible, however, then holds $`m_00.2`$ eV. For a precision of 10 $`\%`$ in the measurement of the minimal $`\stackrel{~}{m}`$ we get a variation in $`m_{\mu \mu }`$ within a factor of 2. A slightly smaller number is obtained for a 5 $`\%`$ precision. In Fig. 10 the area for $`\stackrel{~}{m}=0.1205`$ is obtained for a precision of 5 $`\%`$. However, it seems to be questionable that these values of $`\stackrel{~}{m}`$ and the required precision are feasible, especially in the light of the different results of different nuclear matrix element calculations. Thus one needs to make compromises in terms of theoretical assumptions  in order to get the complete mass matrix. However, despite the large uncertainty in the other mass matrix entries, note that in most cases the mass spectrum of the eigenvalues can still be probed: from Fig. 26 of Ref.  it is seen that — provided the solar solution is known — some measured values of $`m`$ distinguish different schemes. The many possibilities for $`\alpha `$ and $`\beta `$ can be understood in the following way: In a three flavor scheme we have 9 parameters: 3 masses, 3 angles and 3 phases. Oscillation experiments can give two mass squared differences (equivalent to 2 masses), all 3 angles and 1 phase. The mass scale might be given by the tritium spectrum or from cosmology, so that we are left with 2 parameters and one observable, namely $`m`$. Only for minimal $`\stackrel{~}{m}`$ an unambiguous determination of the 2 phases is possible, which however requires an extremely precise measurement of this quantity. The smallness of $`|U_{e3}|`$ allows a broad range for $`\beta `$, but for the $`\mu \mu `$ entry $`\beta `$ contributes with $`c_2c_3e^{i\beta }1/\sqrt{2}e^{i\beta }`$, which can have a large effect on the magnitude of $`m_{\mu \mu }`$.
To conclude, for the $`(+)`$ case the uncertainty in $`m_{\mu \mu }`$ is at least of the same order as the one stemming from the solar angles and $`\stackrel{~}{m}`$. In general the difference can be up to factors around 20.
If we choose a smaller $`|U_{e3}|`$ then the uncertainty in $`\beta `$ is even larger and our arguments are strengthened. For $`m_0`$ smaller then 0.3 eV, Eq.(15) tells us that $`m`$ $`<m_1`$, which is the smallest eigenvalue, lying at the very end of the planned $`m`$ sensitivity. An analysis of the phases is then probably impossible and again our arguments are strengthened.

## 5 Conclusions

In a degenerate mass scheme with three neutrinos we did a full analysis of the allowed parameter range for the relevant function $`\stackrel{~}{m}=m/m_0`$. All $`CP`$–conserving and –violating possibilities for cancellation were considered and plots for observables in each case were given. When $`m_{0,\mathrm{max}}`$ is smaller then roughly $`0.2\mathrm{}0.1`$ eV, the limit applies no longer to a degenerate scheme but for the highest mass of a hierarchical or intermediate scheme. The other masses are then given by the measured mass squared differences. Then however a signal or improved bound from 0$`\nu \beta \beta `$ is needed to distinguish the different possibilities, i.e. to verify which two masses give which $`\mathrm{\Delta }m^2`$. For the SA solution the mass scale is below 0.2 eV, regardless of $`CP`$ violation or conservation. Effects of phases then lie in the $`\%`$ range. We may summarize the situation for $`CP`$ invariance such that the $`(+++)`$ case also means $`m_00.2`$ eV. If the $`(+)`$ configuration is realized (Fig. 2), the total mass scale is too small to be directly measured. Regarding the other two parity configurations (Figs. 4 and 4), maximal solar mixing allows a broad range for $`m_0`$, for small $`|U_{e3}|^2`$ even up to the tritium limit of 2.8 eV. The smaller the angle $`\mathrm{sin}^22\theta _{}`$ is, the smaller is the maximal allowed $`m_0`$ or the larger is the minimal value of $`\stackrel{~}{m}`$. For $`CP`$ violation two phases are present, in principle varying between 0 and $`\pi `$. Because of the smallness of $`|U_{e3}|^2`$ the dependence on one phase is small. For maximal solar mixing complete cancellation and therefore $`m_0`$ up to 2.8 eV is possible. From $`\mathrm{sin}^22\theta _{}0.7`$ on, $`m_{0,\mathrm{max}}`$ approaches 0.2 eV and the degenerate scheme is jeopardized. A key value for $`\stackrel{~}{m}`$ of 0.6 was identified, which would prove the existence of $`CP`$ violation or would rule out some parity configurations. This value can unfortunately only be reached if we loosen the $`m`$ bound with respect to the uncertainty in nuclear matrix element calculations. If we change the limit on $`m`$ by a factor of $`x`$ then all the maximal allowed values for $`m_0`$ obtained here have to be multiplied with this very factor $`x`$.
Provided neutrinos are Majorana particles, pushing $`\stackrel{~}{m}`$ below 0.24 would rule out the LA and VO solution, and further reduction below 0.06 would rule out also the LOW case. If luckily the minimal $`\stackrel{~}{m}`$ would be measured, an uncertainty of 10 $`\%`$ in $`\stackrel{~}{m}`$ translates into a variation of the other mass matrix entries of about a factor of 2. For non–minimal values many choices of the phases are possible. This will reflect on the result of the other entries in the mass matrices. The resulting uncertainty is as least as large as the one resulting from measurements of the solar angle and different calculations of nuclear matrix elements. Because the processes governed by these elements are way beyond experimental access, further input from the theoretical side is needed. Thus, with currently planned projects, experimental verification of a given mass matrix is very questionable.

Acknowledgments

I thank A. Joshipura, E. A. Paschos, W. G. Scott and L. Wolfenstein for helpful comments. This work has been supported in part by the “Bundesministerium für Bildung, Wissenschaft, Forschung und Technologie”, Bonn under contract No. 05HT9PEA5. Financial support from the Graduate College “Erzeugung und Zerf$`\ddot{\mathrm{a}}`$lle von Elementarteilchen” at Dortmund university is gratefully acknowledged.


format: text

text:
# Parametric Luminescence of Microcavity Polaritons

## Abstract

The spectral and dispersive emission properties are analytically determined for the two-dimensional system of exciton-polaritons in microcavities excited by a resonant and coherent optical pump. New collective excitations result from the anomalous coupling between one generic polariton state and its idler, created by the scattering of two pumped polaritons. The corresponding parametric correlation is stimulated by the emitter and idler populations and drives very efficiently the luminescence. The intrinsic properties of the collective excitations determine a peculiar emission pattern.

In their pioneering experiments, Weisbuch et al.  discovered a new kind of two-dimensional quasi-particles, resulting from the strong coupling between quantum well excitons and confined photons in a semiconductor microcavity. These peculiar particles, called microcavity polaritons, have a very sharp dispersion due to the very light mass of the cavity photon. Remarkably, this represents a condensed matter system of small mass quasi-particles which can be manipulated through laser beams both in frequency and momentum space. In principle, the low polariton density of states could allow large occupation numbers at relatively small densities of particles , well below the critical saturation value due to the fermionic nonlinearities. In other words, the polariton system could exhibit bosonic properties. Furthermore, unlike unbound electron-hole pairs in ordinary semiconductor lasers, microcavity polaritons have a relatively short time of recombination into external photons. All these ingredients are indeed very promising for applications in the domain of ultrafast all-optical switching and amplification.

Recently, Dang et al. measured photoluminescence spectra from a II-VI microcavity excited with a nonresonant pump. A threshold was observed in the dependence of the polariton luminescence intensity as a function of the input power. Similar results were reported by Senellart and Bloch in a III-V microcavity . These experiments have been interpreted in terms of enhanced scattering of reservoir excitons into the emitting polariton modes. The origin of the enhancement has been attributed to bosonic stimulation due to final state occupation.

More recently, great insight into the subject has been given by angle-resolved experiments under resonant excitation . In this kind of experiments, polaritons are optically excited at a desired energy and momentum, allowing a direct control of the polariton dynamics. In particular, Savvidis et al. have uncovered a new kind of polariton parametric amplifier through pump-probe experiments. The angular selection rules for the parametric conversion of pumped polaritons into the probe and idler modes are unambiguously given by the energy and momentum conservation for the polariton scattering. In the case of an applied probe, the nonlinearity can be explained in terms of phase-matched wave-mixing of polariton matter beams. Remarkably, giant nonlinearities occur also when the probe beam is absent and the spontaneous emission is detected . Spectrally, the luminescence shows a very peculiar distortion with respect to the bare polariton dispersion. Moreover, the emission is characterized by a finite angular width which appears to be an intrinsic property of the system. Indeed, the whole phenomenology suggests that the interacting polariton matter provides a new kind of collective excitations. The nature of these excitations and the connection with the stimulated scattering mechanism are not yet known. The explicit solution of this intriguing problem is the motivation and main result of our present investigation.

In this paper, we present a theoretical analysis of the polariton nonlinear emission in the case of resonant and coherent pumping. The coherence induced by the pump field allows an anomalous coupling between a generic polariton state and its idler, originated by the fission of two pumped polaritons. The anomalous coupling gives rise to new collective excitations. The luminescence is driven by the emitter-idler parametric correlation which is stimulated by the emitter and idler populations. We establish how the intrinsic properties of the collective excitations determine the rich emission features. In particular, (i) the peculiar dispersion, (ii) the lineshape features and (iii) the two-dimensional pattern are provided by our results.

In order to investigate the nonlinear emission of microcavities in the strong exciton-photon coupling regime, we work directly in the polariton basis. The polariton states have a two-dimensional character and are described by the in-plane wave-vector and their spin. In the following, each $`𝐤`$ will be a two-dimensional vector in the quantum well plane. Moreover, we will consider the case of a circularly polarized pump beam. Neglecting the spin relaxation, only the polariton states with a definite circular polarization will be retained. When showing numerical results, we will use the angular coordinates $`(\theta _x,\theta _y)`$ defined by $`(k_x,k_y)=\frac{\omega }{c}(\mathrm{sin}\theta _x,\mathrm{sin}\theta _y)`$, where $`\omega `$ is the emission frequency. Actually, these coordinates represent the angles of the far-field images in current experiments . Since we consider the situation of resonant pumping of the lower polariton branch, we can neglect the nonlinear contribution due to the upper branch. The destruction operator for a lower polariton with in-plane wave-vector $`𝐤`$ is $`p_𝐤=X_kb_𝐤+C_ka_𝐤`$, where $`b_𝐤`$ and $`a_𝐤`$ are the exciton and photon Bose operators respectively and $`X_k`$, $`C_k`$ the corresponding Hopfield coefficients ($`X_k>0`$ and $`C_k<0`$). The lower polariton Hamiltonian is $`H=H_{LP}+H_{PP}^{eff}+H_{qm}`$. The free term $`H_{LP}=_𝐤E_{LP}(k)p_𝐤^{}p_𝐤`$ contains the lower polariton energy dispersion $`E_{LP}(k)`$. The polariton-polariton interaction term reads

$$H_{PP}^{eff}=\frac{1}{2}\underset{𝐤,𝐤^{},𝐪}{}\frac{\lambda _X^2}{A}V_{𝐤,𝐤^{},𝐪}^{PP}p_{𝐤+𝐪}^{}p_{𝐤^{}𝐪}^{}p_𝐤p_𝐤^{},$$

where $`\lambda _X`$ is the two-dimensional exciton radius and A is the macroscopic quantization area. The interaction potential is $`V_{𝐤,𝐤^{},𝐪}^{PP}=\{2\frac{\mathrm{}\mathrm{\Omega }_R}{n_{sat}\lambda _X^2}(|C_{𝐤+𝐪}|X_𝐤^{}+|C_𝐤^{}|X_{𝐤+𝐪})+\frac{6e^2}{ϵ\lambda _X}X_{𝐤+𝐪}X_𝐤^{}\}X_{𝐤^{}𝐪}X_𝐤`$, with $`n_{sat}=7/(16\pi \lambda _X^2)`$ the exciton saturation density and $`ϵ`$ the quantum well dielectric constant. Notice that $`V_{𝐤,𝐤^{},𝐪}^{PP}`$ is positive and represents a repulsive interaction. The cavity system interacts with the external electromagnetic field through the standard quasi-mode coupling Hamiltonian $`H_{qm}=d\mathrm{\Omega }\{_𝐤gC_𝐤p_𝐤^{}\alpha _{𝐤,\mathrm{\Omega }}+H.c.\}.`$ The operator $`\alpha _{𝐤,\mathrm{\Omega }}`$ destructs an external photon with in-plane wave-vector $`𝐤`$ and frequency $`\mathrm{\Omega }`$. The spectrum of the polariton luminescence is

$$PL(𝐤,t,\omega )|C_k|^2e_0^+\mathrm{}𝑑\tau e^{i(\omega i0^+)\tau }p_𝐤^{}(t+\tau )p_𝐤(t),$$

where $`p_𝐤(t)`$ is the time-dependent polariton operator. Here and in the following we use the Heisenberg representation of time-dependent operators and take expectation values on the stationary state. In the steady-state regime, the two-time correlation $`p_𝐤^{}(t+\tau )p_𝐤(t)`$ depends only on the relative time $`\tau `$, and not on $`t`$. This implies that the stationary spectrum can be calculated at a fixed time $`t=t_0`$ when the steady-state regime is achieved. For simplicity, we choose $`t_0=0`$ and hence $`PL(𝐤,\omega )|C_k|^2e\{\stackrel{~}{p}_𝐤^{}(\omega )p_𝐤(0)\}`$, where $`\stackrel{~}{p}_𝐤(\omega )=_0^+\mathrm{}𝑑\tau e^{i(\omega +i0^+)\tau }p_𝐤(\tau )`$.

An applied cw-optical pump with in-plane wave-vector $`𝐤_𝐩`$ drives a polariton polarization $`p_{𝐤_𝐩}(t)=e^{i\omega _pt}|p_{𝐤_𝐩}|`$, where $`\omega _p`$ is the laser frequency resonant with the lower polariton mode. We are interested in the generic polariton operator $`p_𝐤(t)`$ with $`𝐤𝐤_𝐩`$. The polariton-polariton interaction couples $`p_𝐤(t)`$ to the pumped mode and to the idler polariton operator $`p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}(t)`$, where $`𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}=2𝐤_𝐩𝐤`$. This corresponds to the fission process $`\{𝐤_𝐩,𝐤_𝐩\}\{𝐤,𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}\}`$. Disregarding the pump noise, we can perform the replacement $`p_𝐤^{}(t)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)p_{𝐤_𝐩}(t)p_{𝐤_𝐩}(t)p_𝐤^{}(t)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)p_{𝐤_𝐩}(t)^2.`$ Furthermore, we neglect multiple scattering, that is interaction between modes other than the pumped one (i.e. $`\{𝐤,𝐤^{}\}\{𝐤+𝐪,𝐤^{}𝐪\}`$). The limit of validity of this approximation will be discussed later on. With these assumptions, the Heisenberg equations of motion read

$$i\mathrm{}\frac{d}{dt}\left(\begin{array}{c}p_𝐤(t)\\ p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)e^{i2\omega _pt}\end{array}\right)=\widehat{M}_𝐤^{par}\left(\begin{array}{c}p_𝐤(t)\\ p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)e^{i2\omega _pt}\end{array}\right)$$

$$+\left(\begin{array}{c}F_𝐤(t)\\ F_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)e^{i2\omega _pt}\end{array}\right),$$

where the coupling matrix is

$$\widehat{M}_𝐤^{par}=\left(\begin{array}{cc}\stackrel{~}{E}_{LP}(𝐤)i\gamma _k& E_𝐤^{int}𝒫_{𝐤_𝐩}^2\\ E_𝐤^{int}𝒫_{𝐤_𝐩}^2& 2\mathrm{}\omega _p\stackrel{~}{E}_{LP}(𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}})i\gamma _{k_{idler}}\end{array}\right).$$

The off-diagonal term contains the energy $`E_𝐤^{int}=(V_{𝐤_𝐩,𝐤_𝐩,𝐤𝐤_𝐩}^{PP}+V_{𝐤_𝐩,𝐤_𝐩,𝐤_𝐩𝐤}^{PP})/2`$ and the pump-induced polarization $`𝒫_{𝐤_𝐩}=\frac{\lambda _x}{\sqrt{A}}p_{𝐤_𝐩}`$. We point out that $`|𝒫_{𝐤_𝐩}|^2=n_p\lambda _X^2`$ is the coherent density of pumped polaritons in units of $`\lambda _X^2`$. The diagonal process $`\{𝐤,𝐤_𝐩\}\{𝐤,𝐤_𝐩\}`$ produces the renormalization of the energy dispersion (blueshift), namely $`\stackrel{~}{E}_{LP}(𝐤)=E_{LP}(k)+2V_{𝐤,𝐤_𝐩,\mathrm{𝟎}}^{PP}|𝒫_{𝐤_𝐩}|^2`$. Notice the equality of the diagonal elements of $`\widehat{M}_𝐤^{par}`$ is equivalent to the exact energy conservation for the process $`\{𝐤_𝐩,𝐤_𝐩\}\{𝐤,𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}\}`$. Finally, the coupling to the external photons is responsible for a radiative damping $`\gamma _k=2\pi g^2|C_k|^2/\mathrm{}`$ and a Langevin force $`F_𝐤(t)=𝑑\mathrm{\Omega }gC_𝐤^{}e^{i\mathrm{\Omega }t}\alpha _{𝐤,\mathrm{\Omega }}(0)`$.

From the coupled equations for the operators $`p_𝐤(t)`$ and $`p_𝐤^{}(t)`$, we can obtain the physical quantities of interest. Let us consider the polariton occupation number $`N_𝐤(t)=p_𝐤^{}(t)p_𝐤(t)`$. From the equations for $`p_𝐤(t)`$, we get:

$$\frac{d}{dt}N_𝐤(t)=\frac{2\gamma _k}{\mathrm{}}N_𝐤(t)+\frac{2}{\mathrm{}}m\{p_𝐤^{}(t)F_𝐤(t)$$

$$+E_𝐤^{int}𝒫_{𝐤_𝐩}^2e^{i2\omega _pt}p_𝐤^{}(t)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)\}.$$

This equation shows that the polariton population $`N_𝐤(t)`$ is driven by the anomalous quantum correlation $`p_𝐤^{}(t)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)`$. The anomalous correlation evolution reads:

$$i\mathrm{}\frac{d}{dt}p_𝐤^{}(t)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)=$$

$$[\stackrel{~}{E}_{LP}(𝐤)\stackrel{~}{E}_{LP}(𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}})i(\gamma _k+\gamma _{k_{idler}})]p_𝐤^{}(t)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)$$

$$+p_𝐤^{}(t)F_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)+F_𝐤^{}(t)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)$$

$$E_𝐤^{int}𝒫_{𝐤_𝐩}^2e^{i2\omega _pt}(1+N_𝐤(t)+N_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}(t))).$$

Notice that the emitter-idler correlation $`p_𝐤^{}(t)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t)`$ has a generation term proportional to the Boson enhancement factor $`(1+N_𝐤(t)+N_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}(t))`$. This means that there is a spontaneous inhomogeneous contribution and a stimulation term due to the emitter and idler populations.

In order to obtain the frequency spectrum of the parametric luminescence, we have to come back to the equations for $`p_𝐤(t)`$ and $`p_𝐤^{}(t)`$. In the frequency space (through the transformation $`_0^{\mathrm{}}e^{i(\omega +i0^+)t}𝑑t`$), the coupled equations read

$$\mathrm{}\omega \left(\begin{array}{c}\stackrel{~}{p}_𝐤(\omega )\\ \stackrel{~}{p}_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(2\omega _p\omega )\end{array}\right)=\widehat{M}_𝐤^{par}\left(\begin{array}{c}\stackrel{~}{p}_𝐤(\omega )\\ \stackrel{~}{p}_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(2\omega _p\omega )\end{array}\right)$$

$$+\left(\begin{array}{c}\stackrel{~}{F}_𝐤(\omega )+i\mathrm{}p_𝐤(t=0)\\ \stackrel{~}{F}_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(2\omega _p\omega )+i\mathrm{}p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t=0)\end{array}\right).$$

Neglecting the noise source $`\stackrel{~}{F}_𝐤(\omega )`$ due to the external photons, the solution of this linear inhomogeneous system is

$$\stackrel{~}{p}_𝐤(\omega )=\frac{\mathrm{\Delta }(\omega )i\mathrm{}p_𝐤(t=0)+E_𝐤^{int}𝒫_{𝐤_𝐩}^2i\mathrm{}p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(t=0)}{(\mathrm{}\omega E_{+,𝐤})(\mathrm{}\omega E_{,𝐤})},$$

with $`\mathrm{\Delta }(\omega )=\mathrm{}\omega +\stackrel{~}{E}_{LP}(𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}})+i\gamma _{k_{idler}}2\mathrm{}\omega _p`$. The complex pole energies $`E_{+,𝐤}`$ and $`E_{,𝐤}`$ are the eigenvalues of the parametric matrix $`\widehat{M}_𝐤^{par}`$ and describe the new collective excitations of the microcavity system. Through the expectation value $`e\{\stackrel{~}{p}_𝐤^{}(\omega )p_𝐤(0)\}=e\{p_𝐤^{}(0)\stackrel{~}{p}_𝐤(\omega )\}`$, we obtain the stationary luminescence as a function of the steady state population $`N_𝐤^s=p_𝐤^{}(0)p_𝐤(0)`$ and parametric correlation amplitude $`𝒜_𝐤^s=p_𝐤^{}(0)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^{}(0)`$. The result is

$$PL(𝐤,\omega )|C_k|^2e\left\{i\frac{\mathrm{\Delta }(\omega )N_𝐤^s+E_𝐤^{int}𝒫_{𝐤_𝐩}^2𝒜_𝐤^s}{(\mathrm{}\omega E_{+,𝐤})(\mathrm{}\omega E_{,𝐤})}\right\}.$$

Simple algebra shows that steady-state solutions of the coupled equations for population and parametric correlation are $`p_𝐤(t)p_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}(t)=𝒜_𝐤^se^{i2\omega _pt}`$ and $`N_𝐤(t)=N_𝐤^s`$. The anomalous correlation amplitude reads

$$𝒜_𝐤^s=\frac{E_𝐤^{int}𝒫_{𝐤_𝐩}^2\delta _𝐤}{|\delta _𝐤|^2\frac{(\gamma _k+\gamma _{k_{idler}})^2}{\gamma _k\gamma _{idler}}|E_𝐤^{int}𝒫_{𝐤_𝐩}^2|^2},$$

with $`\delta _𝐤=(2\mathrm{}\omega _p\stackrel{~}{E}_{LP}(𝐤)\stackrel{~}{E}_{LP}(𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}})i(\gamma _k+\gamma _{k_{idler}}))`$. Moreover, the steady-state population is

$$N_𝐤^s=\frac{1}{\gamma _k}m(E_𝐤^{int}𝒫_{𝐤_𝐩}^2𝒜_𝐤^s).$$

Of course, $`N_{𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}}^s=N_𝐤^s\gamma _k/\gamma _{k_{idler}}`$. As we will show below, the PL signal (as well as the parametric correlation amplitude $`𝒜_𝐤^s`$ and population $`N_𝐤^s`$) diverges when the pump density $`|𝒫_{𝐤_𝐩}|^2`$ reaches a threshold value. In fact, above threshold, the pump polarization can be no longer treated as a parameter. This means that the equation of motion for the pumped mode has to be included, accounting for the pump depletion. Moreover, when the transfer from the pump wave-vector to the other ones becomes very massive, multiple scattering can be no longer neglected. Hence, the results here presented are valid only below the threshold of the parametric luminescence.

For sake of illustration, we present our analytical results by using GaAs parameters for a realistic microcavity structure with Rabi splitting equal to 7 meV. We consider the situation of zero exciton-photon detuning, i.e. $`E_{exc}(0)=E_{cav}(0)`$. For simplicity, we consider a k-independent polariton linewidth $`\gamma _k0.4`$ meV (for large k, the decreasing of the radiative width is usually compensated by nonradiative broadening). The two-dimensional dispersion of the lower polariton branch is shown in Fig. 1. We consider a pump which resonantly excites the polariton branch at the critical wave-vector such as $`E_{LP}(0)+E_{LP}(2k_p)=2E_{LP}(k_p)`$, allowing the particular polariton fission $`\{𝐤_𝐩,𝐤_𝐩\}\{\mathrm{𝟎},\mathrm{𝟐}𝐤_𝐩\}`$ (with increasing pumping, one has to consider the renormalized dispersion $`\stackrel{~}{E}_{LP}`$ instead of $`E_{LP}`$). Fig. 2 contains the contour plot of the quantity $`|E_{LP}(𝐤)+E_{LP}(𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}})2E_{LP}(𝐤_𝐩)|`$ in k space ($`𝐤_𝐩`$ is along the x-direction). The white-dashed line represents the zero value contour, i.e. exact energy conservation for the fission $`\{𝐤_𝐩,𝐤_𝐩\}\{𝐤,𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}}\}`$. Remarkably, the conservation is only weakly forbidden in a large squeezed portion of $`𝐤`$-space (dark region) around $`𝐤_𝐩`$. This implies that the parametric coupling is efficient on a broad angular range. This is actually found in Fig.3(a), where a typical contour plot of the photoluminescence (log scale) is shown as a function of the angle $`\theta _x`$ ($`\theta _y=0`$) and of the emission energy. Moreover, Fig. 3(b) shows the two-dimensional pattern of the spectrally integrated emission (linear scale). The depression of the emission at large angles is due to the vanishing cavity fraction $`|C_𝐤|^2`$.

To determine the dependence on the pump intensity, we just need to analyze the evolution of the eigenvalues $`E_{\pm ,𝐤}`$. For a given $`𝐤`$ satisfying exactly the energy-momentum conservation for the parametric scattering, the expression for the pole energies $`E_{\pm ,𝐤}`$ is very simple, i.e. $`e(E_{\pm ,𝐤})=\stackrel{~}{E}_{LP}(𝐤)`$ and $`|m(E_{\pm ,𝐤})|=\left(\sqrt{\gamma _k\gamma _{k_{idler}}}\pm E_𝐤^{int}|𝒫_{𝐤_𝐩}|^2\right)`$. This means that the imaginary part of the pole energy $`E_{,𝐤}`$ gets smaller at the expense of the other pole $`E_{+,𝐤}`$. Consequently, the photoluminescence spectrum $`PL(𝐤,\omega )`$ gets spectrally narrower. The threshold density for the parametric luminescence is defined by $`Im(E_{,𝐤})=0`$, i.e. $`E_{int}|𝒫_{𝐤_𝐩}|^2=\sqrt{\gamma _k\gamma _{k_{idler}}}`$. For this value, the PL signal, the parametric correlation amplitude $`𝒜_𝐤^s`$ and steady-state population $`N_𝐤^s`$ diverge and therefore the pump depletion has to be included. All these features can be seen explicitly in Fig. 4. With increasing pump density the lineshape of the emission around $`k=0`$ blueshifts and narrows (Fig. 4(a)). The output intensity (see the log-log plot inset) shows a threshold behavior as a function of the coherent density of pumped polaritons $`n_p=|𝒫_{𝐤_𝐩}|^2/\lambda _X^2`$. Notice that when the polariton damping $`\gamma _k`$ is small enough, the threshold density $`n_{thr}`$ can be orders of magnitude smaller than the exciton saturation density $`n_{sat}`$. With the realistic material parameters employed here, we have $`n_{thr}0.05n_{sat}`$. Finally, Fig.4(b), shows explicitly the energy dispersions of the collective excitations. When the energy conservation for the parametric fission is well satisfied, $`e(E_\pm (𝐤))\stackrel{~}{E}_{LP}(𝐤)`$. On the other hand, when the energy conservation is strongly broken and consequently the parametric luminescence is negligible, $`e(E_{}(𝐤))`$ and $`e(E_+(𝐤))`$ bifurcate, tending to the diagonal elements of the matrix $`M_𝐤^{par}`$, i.e. $`\stackrel{~}{E}_{LP}(𝐤)`$ and $`2\mathrm{}\omega _p\stackrel{~}{E}_{LP}(𝐤_{\mathrm{𝐢𝐝𝐥𝐞𝐫}})`$. In the intermediate case, consistently with experiments, the dispersion $`e(E_\pm (𝐤))`$ shows peculiar features such a change of sign of the in-plane group velocity  around $`\mathrm{𝟐}𝐤_𝐩`$ and a flattening around $`𝐤=\mathrm{𝟎}`$ (see also the emission spectra in Fig. 3).

In conclusion, we have theoretically described the parametric luminescence of microcavity polaritons, showing the subtle interplay between parametric correlation and stimulated scattering. The new collective excitations resulting from the anomalous emitter-idler coupling produce a very peculiar luminescence pattern in frequency and momentum space, giving a key to recent experiments . Further investigations are in progress for the regime above threshold not treated in this paper, but investigated in new experiments. We wish to thank J.J. Baumberg, B. Deveaud, R. Houdré, M. Saba, P.G. Savvidis, R.P. Stanley, F. Tassone and C. Weisbuch for fruitful discussions.


text:
# Surface tension of the isotropic-nematic interface

## Abstract

We present the first calculations of the pressure tensor profile in the vicinity of the planar interface between isotropic liquid and nematic liquid crystal, using Onsager’s density functional theory and computer simulation. When the liquid crystal director is aligned parallel to the interface, the situation of lowest free energy, there is a large tension on the nematic side of the interface and a small compressive region on the isotropic side. By contrast, for perpendicular alignment, the tension is on the isotropic side. There is excellent agreement between theory and simulation both in the forms of the pressure tensor profiles, and the values of the surface tension.

The nematic liquid crystal phase is distinguished from the isotropic liquid by the existence of long-ranged orientational ordering of the constituent molecules . The structure of the planar interface between coexisting nematic (N) and isotropic (I) phases is of fundamental interest: along with the liquid-vapour surface, it is one of the simplest fluid-fluid interfaces. Molecular-scale interactions, capillary-wave fluctuations, and the coupling between translational and rotational degrees of freedom, may all influence the interface properties. The surface excess free energy, or interfacial tension $`\gamma `$, is a key parameter. The variation of this quantity with the direction of molecular ordering in the liquid crystal, the so-called director, dictates the molecular alignment which will be preferred by the surface: the phenomenon of anchoring, an essential feature of liquid crystal physics . Determining the surface tension of the N-I interface for a given molecular model by computer simulation, and predicting it from first principles by theoretical methods, is clearly desirable. For a planar interface, $`\gamma `$ may be calculated by integrating the difference between normal $`P_\mathrm{N}`$ and transverse $`P_\mathrm{T}`$ pressure tensor components with respect to position : $`\gamma =_{\mathrm{}}^{\mathrm{}}dz\left[P_\mathrm{N}(z)P_\mathrm{T}(z)\right]`$.

In this paper, we present the first discussion of the pressure tensor profile through the nematic-isotropic interface in the framework of a simple density functional theory in the Onsager approximation , modelling the molecules as freely-rotating and translating rod-like rigid bodies. We shall show that there is a significant difference between the two situations of most interest, when the director is oriented respectively parallel to, and normal to, the interface. Previous studies , calculating the surface tension from free energy differences, have indicated that, in the limit of infinite rod elongation, the minimum surface free energy is obtained when the rods lie parallel to the N-I interface; in other words, this interface favours planar anchoring. We provide additional insight into this result, through an examination of the normal and transverse pressure.

We back up these predictions by carrying out the first computer simulations to give the pressure tensor profiles for equilibrium N-I interfaces, confined between parallel walls under different anchoring conditions, using the same molecular models as in the theory, Our results are further supported by the first computer simulations of free standing nematic films in coexistence with the isotropic phase. From the simulation viewpoint, this study is much more challenging than determining the surface tension of the liquid-vapour interface  since the densities and other properties of the two phases are very similar. To our knowledge, the only other such studies have involved two phases maintained out of thermodynamic equilibrium .

We employ two very similar molecular models in this work, in both simulation and theory. Both models are axially symmetric, with molecular orientations represented by unit vectors $`𝐮_i`$ and centre-of-mass positions by vectors $`𝐫_i`$. In the first model, the molecules are taken to be hard ellipsoids of revolution of axial ratio $`\kappa =A/B=15`$ where $`A`$ is the length of the symmetry axis and $`B`$ that of the transverse axes; the bulk phase diagram of this model is well known , and the interfacial structure has been recently reported by one of us . For hard particle models, the temperature enters only trivially into the thermodynamics, and we define units of energy by setting $`k_\mathrm{B}T=1`$ in this work. Length units are fixed by taking $`B=1`$, mass units by setting particle mass $`M=1`$, and the moment of inertia was fixed at $`I/MB^2=50`$. In the second ‘soft ellipsoid’ model, a continuous pair potential $`v_{12}=4ϵ(s_{12}^{12}s_{12}^6)+1`$, for $`s_{12}<2^{1/6}`$, $`v_{12}=0`$ otherwise, was used. Here, $`s_{12}=(r_{12}\sigma _{12}+B)/B`$ is a scaled and shifted separation where the ‘diameter’ $`\sigma _{12}=\sigma (\widehat{𝐫}_{12},𝐮_1,𝐮_2;A,B)`$ is a standard approximation  to the contact distance of two ellipsoids (again with $`A/B=15`$ here) at given relative orientations; $`𝐫_{12}=𝐫_1𝐫_2`$, $`r_{12}=|𝐫_{12}|`$ and $`\widehat{𝐫}_{12}=𝐫_{12}/r_{12}`$. The potential is purely repulsive and varies quite strongly over a short range of interparticle separation. We choose an energy parameter $`ϵ=1`$, again take $`k_\mathrm{B}T=1`$, and set other parameters as for the hard ellipsoids.

In planar geometry, with translational invariance in the $`x`$ and $`y`$ directions, we adopt the Irving-Kirkwood convention  so as to write components of the pressure tensor in the form

$`P_{\alpha \alpha }(z)=k_\mathrm{B}T\rho (z)`$ (4)
$`{\displaystyle \frac{1}{2}}{\displaystyle d𝐫_{12}d𝐮_1d𝐮_2\left(r_{12}^\alpha \frac{v_{12}}{r_{12}^\alpha }\right)}`$
$`\times {\displaystyle _0^1}\mathrm{d}\zeta \varrho ^{(2)}(𝐫_{12},z\zeta z_{12},z+(1\zeta )z_{12};𝐮_1,𝐮_2)`$
$`\alpha =x,y,z.`$
Here, $`\rho (z)=d𝐮\varrho ^{(1)}(z,𝐮)`$ is the number density, $`\varrho ^{(1)}(z,𝐮)`$ the (position- and orientation-dependent) single-particle density, and $`\varrho ^{(2)}(𝐫_{12},z_1,z_2;𝐮_1,𝐮_2)`$ the pair density. The normal component $`P_\mathrm{N}(z)=P_{zz}(z)=P`$ is independent of $`z`$ throughout the system, while the transverse component $`P_\mathrm{T}(z)=P_{xx}(z)=P_{yy}(z)`$ is not. Far from the interface, the pressure tensor becomes isotropic, $`P_\mathrm{N}(z)=P_\mathrm{T}(z)=P`$.

For simulations of the soft-ellipsoid system, the pressure tensor was routinely computed from Eqn (4) via the intermolecular forces . For the hard-ellipsoid system, it was calculated by identifying nearly-contacting pairs of ellipsoids . Such pairs are defined as those which would overlap if their dimensions are uniformly scaled by a factor $`\lambda `$, very slightly larger than unity. Tests were carried out with several choices of $`\lambda `$: all of the results here are obtained with $`\lambda =1.01`$, the smallest value to give acceptable statistics, which was also found to be satisfactory in ref .

The construction and minimization of the Onsager free-energy functional followed closely the procedure described before . The excess part of the Helmholtz free energy $`=^{\mathrm{id}}+^{\mathrm{ex}}`$ is expressed in terms of $`\varrho ^{(1)}(𝐫,𝐮)`$:

$`^{\mathrm{ex}}[\varrho ^{(1)}]`$ $`=`$ $`{\displaystyle \frac{1}{2}}k_\mathrm{B}T{\displaystyle d𝐫_1d𝐮_1d𝐫_2d𝐮_2}`$ (6)
$`\times \varrho ^{(1)}(𝐫_1,𝐮_1)\varrho ^{(1)}(𝐫_2,𝐮_2)f(𝐫_{12},𝐮_1,𝐮_2).`$
where $`f(𝐫_{12},𝐮_1,𝐮_2)=\mathrm{exp}(v_{12}/k_\mathrm{B}T)1`$; this is the leading term in a low-density virial expansion. Fluctuation effects, which come in through higher-order terms, were neglected. The most important fluctuations are capillary waves, which only slightly broaden the profiles in small systems, and do not affect the interfacial tension. Since $`\varrho ^{(1)}`$ is independent of $`x`$ and $`y`$, $`f`$ may be integrated over these coordinates. The logarithm of the density is expanded in spherical harmonics, on a discrete grid in $`z`$ with an interval $`\delta z=B/5=A/75`$: $`\mathrm{ln}\varrho ^{(1)}(z_i,𝐮)=_{lm}C_{i;lm}Y_{lm}(𝐮).`$ Periodic boundary conditions apply in the $`z`$ direction. In the current application, this expansion was taken to fourth order in the $`Y_{lm}(𝐮)`$ (restricting to second order was found to change the coexistence pressure and densities by $``$1%). All relevant orientational integrals were expressed in terms of spherical harmonic coefficients of the density up to order 10 (reducing to order 8 was found to change the coexistence pressure and densities by $``$0.01%). The parameters $`C_{i;lm}`$ were varied to minimize the grand potential $`\beta \mathrm{\Omega }[\varrho ^{(1)}]\beta [\varrho ^{(1)}]\beta \mu N`$ of an I-N film system at the coexistence value of the chemical potential $`\mu `$, determined by a preliminary bulk calculation.

The normal pressure is constant through the system, and equal to $`\mathrm{\Omega }/V`$ in both bulk phases at coexistence. The transverse pressure tensor profile was calculated directly from the minimized density profiles. The definition consistent with Onsager’s approximation and the use of the Irving-Kirkwood convention in the simulations is easily derived from eqn (4), by inserting the low-density form $`\varrho ^{(2)}(𝐫_{12},z_1,z_2;𝐮_1,𝐮_2)\varrho ^{(1)}(z_1,𝐮_1)\varrho ^{(1)}(z_2,𝐮_2)\mathrm{exp}(v_{12}/k_\mathrm{B}T)`$ and carrying out an integration by parts. The result is

$`P_\mathrm{T}(z)=k_\mathrm{B}T\rho (z)`$ (9)
$`{\displaystyle \frac{1}{2}}k_\mathrm{B}T{\displaystyle d𝐫_{12}d𝐮_1d𝐮_2f(𝐫_{12},𝐮_1,𝐮_2)}`$
$`\times {\displaystyle _0^1}\mathrm{d}\zeta \varrho ^{(1)}(z\zeta z_{12},𝐮_1)\varrho ^{(1)}(z+(1\zeta )z_{12},𝐮_2).`$
All simulations were carried out at constant volume. A cuboidal simulation box was employed in all cases, with periodic boundaries in the $`x`$ and $`y`$ directions. Conventional, microcanonical ensemble molecular dynamics (MD) simulations of soft ellipsoids employed periodic boundaries also in the $`z`$ direction, and used a reduced timestep $`\delta t=0.002`$. Several systems of $`N=7200`$ particles in boxes of size $`L_x/A=L_y/A=2.5,L_z/A20`$ were studied. At average number densities $`0.24\rho AB^20.255`$, from an initially isotropic configuration, one or several nematic films were observed to assemble spontaneously, with planar interfaces normal to the $`z`$-direction and, in every case, the director lying parallel to the interface. In addition, from an initially nematic configuration aligned in the $`z`$ direction, a similar two-phase system nucleated with the nematic director turned through 90 degrees to lie parallel to the interface, within roughly 1 million time steps.

One of the systems at density $`\rho AB^2=0.255`$, containing a nematic film of width $`13B`$, was further equilibrated over 7 million time steps, monitoring the pressure tensor profiles; run averages were then accumulated over 4.7 million steps. In the absence of walls, the interfacial tension can be calculated simply from the difference between the system averages of $`P_\mathrm{N}`$ and $`P_\mathrm{T}`$. This has been done independently using the transverse component $`P_{yy}`$, which happens to be roughly parallel to the director, and the transverse component $`P_{xx}`$, which is perpendicular, giving respectively $`\gamma B^2/k_BT=0.022\pm 0.01`$, $`0.015\pm 0.01`$. Thus, the free film simulations yielded clear evidence that the anchoring at the interface is planar, provided a rough estimate of the interfacial tension, and indicated an unexpected difference between the transverse stress components $`P_{xx}`$ and $`P_{yy}`$. To provide more convenient control of the director orientation relative to the interface, a more systematic study was therefore conducted in a confined geometry.

Both hard and soft ellipsoids were simulated using conventional constant-$`NVT`$ Monte Carlo (MC)  using $`N8000`$ particles in a box of dimensions $`L_x/A=L_y/A=2.22`$, $`L_z/A=22.73`$. Flat walls with a short-range orienting field at the box ends allowed nematic films of width approximately $`8A`$ to be stabilized in a desired director orientation, separated by a central isotropic region of width approximately $`7A`$. Full details of the system, together with a description of the film structure, have been given elsewhere . In this paper we focus on the cases of normal (along $`z`$) and in-plane (along $`y`$) orientation of the director relative to the interface. For each model, and each director alignment, we report the results of averaging over eight independent runs. For the hard ellipsoids, more than $`10^6`$ Monte Carlo sweeps were allowed for equilibration, and the averages were accumulated over a further $`7.5\times 10^5`$ sweeps (one sweep is one attempted move per particle). For the soft ellipsoids, the corresponding figures are $`3.5\times 10^5`$ sweeps and $`6.25\times 10^5`$ sweeps. Profiles of number density $`\rho (z)`$, orientational order parameter $`S(z)`$ (the largest eigenvalue of the second rank order tensor ) and the pressure tensor were accumulated. To facilitate the precise comparison between simulation profiles and theoretical predictions, a correction was made for the fluctuations in the interface position, as described elsewhere . Error bars were estimated by assuming a normal distribution of the results of the eight separate runs. For the in-plane cases of director alignment along $`y`$ we detected (as for the free films) an extremely small discrepancy $`P_{zz}=P_{xx}P_{yy}`$ in the nematic phase far from the interface. We ascribe this to correlations of aligned molecules across the periodic box: to eliminate this finite-size effect from the very sensitive surface tension calculation we set $`P_\mathrm{T}(z)=P_{xx}(z)`$ in these cases. For normal alignment, we found the $`x`$ and $`y`$ directions to be equivalent and set $`P_\mathrm{T}(z)=(P_{xx}(z)+P_{yy}(z))/2`$.

Coexistence parameters from theory and simulation are compared in Table I. As is well known , for elongation $`A/B=15`$ the Onsager theory overestimates the coexistence pressures and densities by up to 30%. Profiles in the vicinity of the interface are compared in Figures 1 and 2. To help locate the interface, we show the number density curves for both simulation and theory, scaled by the bulk coexistence values $`\rho _\mathrm{I},\rho _\mathrm{N}`$ thus: $`\rho ^{}(z)=(\rho (z)\rho _\mathrm{I})/(\rho _\mathrm{N}\rho _\mathrm{I})`$, and similarly for the nematic order parameter $`S^{}(z)`$. The zero of $`z`$ is taken at the point of inflection of $`S^{}(z)`$; as discussed elsewhere , the density profile lies at $`zA/3`$ relative to this (on the nematic side). In the figures we show both the pressure tensor difference $`P_\mathrm{N}P_\mathrm{T}(z)`$ and the running integral $`_{z_0}^zdzP_\mathrm{N}P_\mathrm{T}(z)`$ taken from a point $`z_0`$ in the bulk nematic phase far from the interface. The results, for the two different director alignments, are dramatically different. For parallel, in-plane, alignment, there is a large tension (low transverse pressure) on the nematic side of the interface, followed by a small compressive region (high transverse pressure) on the isotropic side, as seen in simulations of the liquid-vapour interface of a simple fluid . The surface of tension defined by the first moment of $`P_\mathrm{N}P_\mathrm{T}(z)`$  lies on the denser side of the interface. For normal alignment, there is a compressive region on the nematic side, and a region under tension on the isotropic side: the surface of tension lies on the side of the less dense phase. Such behaviour is seen in the idealized penetrable sphere model of the liquid vapour interface  but to our knowledge this is the first time it has been observed for a molecular model with realistic interactions. The step change in the running integrals gives estimates of the surface tension, reported in Table I. We find that the planar orientation is favoured, as it has the lower excess free energy; the same result was obtained by Onsager theory in the infinite elongation limit , although the corresponding surface tensions in that case (appropriately scaled by $`A`$ and $`B`$) are about 40% larger than here.

Although the simulation results are subject to significant error bars, we have succeeded in resolving the variation of transverse pressure through the interface, and the profiles are very similar, both in form and magnitude, to those predicted by theory. Indeed, the level of agreement between theory and simulation results is remarkable, bearing in mind the absence of adjustable parameters in the Onsager theory. In principle, the location of the surface of tension depends on the convention adopted for the transverse stress, e.g. Irving-Kirkwood vs Harasima , but in the current study we find that the profiles are only slightly dependent on this choice; moreover, the results for such molecular shapes should be indicative of the behaviour for more elongated particles, so we expect them to be rather generally applicable.

This work was supported by the Engineering and Physical Sciences Research Council, the German Science Foundation (DFG), the Alexander von Humboldt Foundation, the Leverhulme Trust, and the British Council. AJM acknowledges support through grant ORS/99007016 of the Overseas Research Students Award. Helpful conversations with Bob Evans and René van Roij are gratefully acknowledged. Some of this work was performed while MPA and AJM visited the Max Planck Institute for Polymer Research, Mainz, and the Institute of Physics, University of Mainz; the authors are grateful to Kurt Kremer and Kurt Binder for their hospitality.


text:
# The one-dimensional 𝑋𝑋𝑍 model with long-range interactions

## Abstract

The one-dimensional $`XXZ`$ model ($`s=1/2`$, $`N`$ sites) with uniform long-range interactions among the transverse components of the spins is considered. The Hamiltonian of the model is explicitly given by $`H=J_{j=1}^N\left(s_j^xs_{j+1}^x+s_j^ys_{j+1}^y\right)(I/N)_{j,k=1}^Ns_j^zs_k^zh_{j=1}^Ns_j^z,`$ where the $`s^{x,y,z}`$ are half the Pauli spin matrices. The model is exactly solved by applying the Jordan-Wigner fermionization, followed by a Gaussian transformation. In the absence of the long-range interactions ($`I=0`$), the model, which reduces to the isotropic $`XY`$ model, is known to exhibit a second-order quantum phase transition driven by the field at zero temperature. It is shown that in the presence of the long-range interactions ($`I0`$) the nature of the transition is strongly affected. For $`I>0`$, which favours the ordering of the transverse components of the spins, the transition is changed from second- to first-order, due to the competition between transverse and $`xy`$ couplings. On the other hand, for $`I<0`$, which induces complete frustration of the spins, a second-order transition is still present, although the system is driven out of its usual universality class, and its critical exponents assume typical mean-field values.

1Departamento de Física da UFC, Campus do Pici, Cx. Postal 6030, 60451-970 Fortaleza (CE), Brazil

2Instituto de Física, Universidade de São Paulo, Cx. Postal 66318, 05315-970 São Paulo (SP), Brazil

The observed critical behaviour of magnetic materials, in the very low temperature limit, has renewed the interest in the study of magnetic quantum transitions (S.L.Sondhi,S.M.Girvin,J.P.Carini and D.Shahar, Rev. Mod. Phys. 69, 315(1997)) . Since these transitions, which are governed by quantum fluctuations, occur at T=0, the one-dimensional models play a important role in their study . Therefore, we will consider the exactly soluble one-dimensional $`XXZ`$ model ($`s=1/2`$), with a uniform long-range interaction among the spins along the $`z`$ direction. Due to the long-range interaction, the model also presents classical critical behaviour, with transitions of first and second order, and it has already been considered by Suzuki (J. Phys. Soc. Jpn. 21, 2140 (1966)). Since his study was restricted to the analysis of the classical second order transition of the model, and we are interested in its quantum transitions, the model will be considered again. In particular, we will be interested in the effect of the long-range interaction on its quantum critical behaviour.

The Hamiltonian of the model is given by

$$H=J\underset{j=1}{\overset{N}{}}(s_j^xs_{j+1}^x+s_j^ys_{j+1}^y)\frac{I}{N}\underset{j,k=1}{\overset{N}{}}s_j^zs_k^zh\underset{j=1}{\overset{N}{}}s_j^z,$$
(1)
where $`J>0,N`$ is the number of sites of the lattice and we assume periodic boundary conditions. By applying the Jordan-Wigner fermionization, followed by a Gaussian transformation, we can write the partition function of the model as

$$Z_N=\text{Tr}e^{\beta H}=C(\beta )e^{\frac{N}{2}z^2}\text{Tr}e^{\stackrel{~}{H}(z)}𝑑z,$$
(2)
with

$$\stackrel{~}{H}(z)=\frac{\beta J}{2}\underset{j=1}{\overset{N}{}}(c_j^{}c_{j+1}+c_{j+1}^{}c_j)\stackrel{~}{h}(z)\underset{j=1}{\overset{N}{}}c_j^{}c_j,$$
(3)
where $`\stackrel{~}{h}(z)=\beta (hI)+\sqrt{2\beta I}z`$, $`C(\beta )`$ depends only on the temperature, a boundary term has been neglected in $`\stackrel{~}{H}(z)`$, and the $`c_j`$ are fermion operators.

Introducing the Fourier transforms

$$c_j=\frac{1}{\sqrt{N}}e^{ikj}\widehat{c}_k,$$
(4)
we can rewrite $`\stackrel{~}{H}(z)`$ in the diagonal form

$$\stackrel{~}{H}(z)=\underset{k}{}\stackrel{~}{\epsilon }_k(z)\widehat{c}_k^{}\widehat{c}_k,$$
(5)
where $`\stackrel{~}{\epsilon }_k(z)=\beta J\mathrm{cos}k\stackrel{~}{h}(z)`$ and, due to the periodic boundary conditions, $`k=2\pi n/N`$ ($`n=1,\mathrm{},N)`$. The partition function is then given by

$$Z_N=C(\beta )e^{\frac{N}{2}z^2}\underset{k}{}\left[1+e^{\stackrel{~}{\epsilon }_k(z)}\right]dz,$$
(6)
which, in the thermodynamic limit ($`N\mathrm{}`$), can be evaluated by the saddle-point method. By explicit calculation, we conclude that

$$m\frac{1}{N}\underset{j=1}{\overset{N}{}}s_j^z=\frac{z_0}{\sqrt{2\beta I}}\frac{1}{2},$$
(7)
where $`z_0`$ is the value of $`z`$ which makes the integrand in Eq. (6) a maximum.

Noticing that $`z_0/\sqrt{2\beta I}`$ is just the average number of fermions per energy level, we can write the equation of state of the system,

$$m=\frac{1}{2\pi }_0^{2\pi }\frac{dk}{1+e^{\stackrel{~}{\epsilon }_k(m)}}\frac{1}{2},$$
(8)
where $`\stackrel{~}{\epsilon }_k(m)=\beta J\mathrm{cos}k+\beta (h+2Im)`$. In the limit $`T0`$ ($`\beta \mathrm{}`$), for $`(h+2Im)/J,`$ Eq. (8) takes the form

$$m+\frac{1}{2}=\frac{1}{\pi }\mathrm{arccos}\left(\frac{h+2Im}{J}\right),$$
(9)
which, for $`I=0`$, readily reduces to the well-known expression for the $`XX`$ chain . To analyze the behavior of the model near the quantum critical point, assuming $`h0`$, we define the order parameter  $`\sigma \frac{1}{2}m`$ and expand Eq. (9) to second order in $`\sigma 0^+`$, obtaining

$$\frac{\pi ^2}{2}\sigma ^2\frac{2I}{J}\sigma =\frac{h_ch}{J}0,$$
(10)
where $`h_c=JI`$. For $`I=0`$ we regain the usual $`XX`$ chain result

$$\sigma (h_ch)^{1/2},$$
(11)
while for $`I<0`$ we get the expected mean-field scaling form

$$\sigma (h_ch)^1.$$
(12)
Note that Eq. (10) cannot be satisfied for $`I>0`$, an indication that in this case the model undergoes a first-order transition at $`h=h_c`$ to a state where the transverse magnetization is saturated ($`m=1/2`$). In this case, there is a hysteresis cycle associated to the transition, which is due to the presence of the metastable states. These states can be identified by looking at the free energy which, for $`(h+2Im)/J`$, and in the limit $`T0`$, is explicitly given by

$$f=\frac{h}{2}\frac{J}{\pi }\left[\mathrm{sin}\phi +\left(\frac{h+2Im}{J}\right)\phi \right]+I(m^2+m),$$
(13)
where $`\phi `$ is defined as

$$\phi =\mathrm{arccos}\left(\frac{h+2Im}{J}\right).$$
(14)
Taking the limit $`h0`$ in Eq. (13), and by imposing that $`f(0)=f(1/2),`$ which are minima of the free energy, we can show that the systems presents spontaneous magnetization for $`I/J4/\pi `$ .

The previous analysis allows us to determine the phase diagram of the model at zero temperature, shown in Figure 1. Notice that there must be a finite temperature critical line ending at the point $`(h/J,I/J)=(1,0)`$, which is thus analogous to a bicritical point. The finite temperature behavior of the model will be considered in future work.

Acknowledments. This work was partially financed by the Brazilian agencies CNPq, FINEP and Fapesp. A. P. Vieira thanks T. A. S. Haddad and S. R. Salinas for useful discussions.

Figure caption.

Figure 1. Phase diagram of the model at $`T=0`$. The solid and dashed lines indicate second- and first-order phase transitions, respectively. The diagram has of course mirror symmetry with respect to the $`I/J`$ axis.