text:
# Halting planet migration by photoevaporation from the central source

## 1 Introduction

It is believed that gas giants do not generally form at small orbital distances from the central star (Boss, 1995). Thus, a natural explanation for extra solar planets orbiting close to the central star is that these planets formed further away in the protoplanetary disk and migrated inward to where they are now observed. A variety of mechanisms have been proposed to explain planet migration: the interaction between a planet and a planetesimal disk (Murray et al., 1998), the gravitational interaction between two or more Jupiter mass planets (Rasio & Ford, 1996), and the tidal gravitational interaction between the planet and the surrounding disk gas (Goldreich & Tremaine, 1979, 1980). The last mechanism, focused on in this paper, is expected to be dominant at early times since the surrounding gaseous disk is required for the formation of planets.

If the perturbation exerted on the disk by the planet is small, the disk structure is not greatly altered, and the planet moves inward relative to the surrounding gas (Ward, 1997). This type of migration is referred to as ‘type-I’. However, if the planet is large it may open a gap in the disk (Goldreich & Tremaine, 1980). The planet is locked to the disk and moves either inward or outward in lock-step with the gaseous disk. This slower migration is referred to as ‘type-II’.

We propose a mechanism for halting type-II migration: photoevaporation driven by radiation from the central star. The planet’s final location is consistent with the Solar System and the growing class of extra solar planets with nearly circular orbits outside of a few AU (see Tinney et al., 2002, figure 4). Photoevaporation by the central star was proposed by Shu et al. (1993) and Hollenbach et al. (1994) as a way to remove a gas disk. Hollenbach et al. (2000) generalized the discussion, describing the variety of possible disk removal mechanisms: accretion, planet formation, stellar encounters, stellar winds or disk winds, and photoevaporation by ultraviolet photons from the central source or massive external stars. Hollenbach et al. (2000) concluded that the dominant mechanisms for a wide range of disk sizes are viscous accretion and photoevaporation, operating in concert within the disk. In this paper, we consider photoevaporation by the central source and viscous accretion.

## 2 Model

The model for disk removal used here is developed in a paper by Matsuyama, Johnstone, & Hartmann (2003) and is similar to that used by Clarke, Gendrin, & Sotomayor (2001). Additionally, we assume a planet with a large mass, which opens a narrow gap in the disk, and assume that planet migration proceeds in lock-step with the disk evolution (type-II migration). The gas disk orbiting the central star with the local Keplerian circular velocity, $`v_k`$, is axisymmetric and geometrically thin. Considering angular momentum and mass conservation of a disk annulus at a radius, $`R`$, with kinematic viscosity, $`\nu `$, the disk surface density evolution can be described by (Pringle, 1981)

$$\frac{\mathrm{\Sigma }}{t}=\frac{3}{R}\frac{}{R}\left[R^{1/2}\frac{}{R}\left(\nu \mathrm{\Sigma }R^{1/2}\right)\right],$$
(1)
where $`\mathrm{\Sigma }`$ is the disk surface density and $`t`$ is the disk evolutionary time. We adopt the standard $`\alpha `$-parameterization of Shakura & Sunyaev (1973) and write $`\nu =\alpha c_sH,`$ where $`c_s`$ is the sound speed at the disk mid-plane and $`H`$ is the disk thickness. Hartmann et al. (1998) estimate $`\alpha 10^210^3`$ and a disk mass $`0.010.2M_{\mathrm{}}`$ for T Tauri stars (TTSs). For the modeled disk temperature distribution, $`T_dR^{1/2}`$ (D’Alessio et al., 1998), the viscosity takes the simple form, $`\nu (R)R`$. Given a solution for equation 1, we find the drift velocity,

$$v_R=\frac{3}{\mathrm{\Sigma }R^{1/2}}\frac{}{R}(\nu \mathrm{\Sigma }R^{1/2}),$$
(2)
and describe the evolution of the disk stream lines.

The EUV ($`\lambda <912`$Å) photons from the central star and the accretion shock (we will refer to these photons as the EUV photons from the central source) are capable of ionizing hydrogen and evaporating material from the disk surface. Photoevaporation forms an ionized atmosphere above the thin viscous disk. Hollenbach et al. (1994) found analytic solutions for the photoevaporation mass loss rate by EUV photons from the central source. These photons are attenuated by recombined hydrogen atoms and scattered in the ionized atmosphere, providing a source of diffuse EUV photons that dominates the flux onto the disk at radii much larger than the size of the central star. Disk material is gravitationally bound to the central star inside the gravitational radius, $`R_g=GM_{}/3c_s^2`$, and it flows out of the disk surface at the sound speed of the heated material, $`c_s^{}`$, beyond this radius. Given the number density of ionized hydrogen at the base of the ionized layer, $`n(R)`$, we can calculate the evaporation rate:

$$\dot{\mathrm{\Sigma }}_{ph}(R)=\{\begin{array}{cc}2\mu m_pn(R)c_s^{},\hfill & \text{if }R>R_g\text{;}\hfill \\ 0,\hfill & \text{otherwise,}\hfill \end{array}$$
(3)
where $`\mu =0.68`$ is the mean molecular weight for the ionized material.

Assuming ionization equilibrium Hollenbach et al. (1994) found the number density at the base of the ionized layer, for $`R>R_g`$:

$$n(R)=3.1\times 10^5\left(\frac{\varphi }{10^{40}\text{s}^1}\right)^{1/2}\left(\frac{R_g}{1AU}\right)\left(\frac{R}{1AU}\right)^{5/2}\text{cm}^3;$$
(4)
where $`\varphi `$ is the total ionizing flux from the central source. We assume that half of the accretion luminosity is radiated as hot continuum from accretion shocks and that the EUV flux can be characterized as blackbody emission at $`T_{as}=10^4`$K (Gullbring et al., 2000; Johns-Krull et al., 2000; Kenyon et al., 1989). There are currently no strong observational constraints on the ionizing flux from the central star, $`\varphi _{}`$. The EUV flux for the sun is $`2.6\times 10^{37}\text{s}^1`$ (see Wilhelm et al., 1998, figure 10b). Since TTSs are chromospherically active (D’Alessio et al., 1998), their EUV fluxes must be higher by few orders of magnitude. We estimate the ionizing flux from the FUV spectra of TTSs (Valenti et al., 2000, figure 3) and the corresponding distances to the star (Johns-Krull et al., 2000, table 1). The typical values are $`10^{39}\text{s}^1`$ for classical T Tauri stars. Extrapolating the FUV flux underestimates the EUV flux; thus the EUV flux might be as high as $`10^{40}\text{s}^1`$.

Combining photoevaporation with viscous accretion is done numerically. At each time step, photoevaporation induced mass loss and viscous diffusion induced disk evolution are solved. For details on the disk dispersal calculation refer to Matsuyama et al. (2003). The simulation is stopped when the disk mass is $`1M_J`$.

## 3 Discussion

At a given disk radius, a gap forms when the mass transported by viscous accretion is equal to the mass removed by photoevaporation. This occurs when the photoevaporation time scale, $`\mathrm{\Sigma }/\dot{\mathrm{\Sigma }}_{ph}`$, is equal to the viscous diffusion time scale, $`R^2/3\nu `$. Therefore, we can estimate the surface density for gap formation at the gravitational radius,

$$\mathrm{\Sigma }_{gap}(2\text{ g\hspace{0.17em}cm}^2)\left(\frac{\alpha }{10^3}\right)^1\left(\frac{\varphi _{}}{10^{40}\text{s}^1}\right)^{1/2}\left(\frac{R_g}{2.4\text{AU}}\right)^1,$$
(5)
Figure 1 shows snapshots of the disk surface density distribution for a disk with initial mass, $`M(0)=0.03M_{}`$ and $`\alpha =10^3`$. A gap starts forming at the gravitational radius when the surface density reaches $`\mathrm{\Sigma }_{gap}2\text{ g\hspace{0.17em}cm}^2`$ during the last stages ($`t_{gap}4\times 10^7`$yr) of the disk evolution and the disk is divided into two annuli. The subsequent evolution is dominated by two competing effects; viscous diffusion attempts to spread both annuli and remove the gap structure while photoevaporation removes material predominantly at the gravitational radius, reopening the gap. The material in the inner annulus is quickly removed both at the inner edge, by accretion onto the central star, and at the outer edge, by the combination of viscous spreading and photoevaporation. In contrast, the outer annulus loses material primarily from its inner edge. The inner edge of the outer annulus recedes from the star as photoevaporation removes disk material, thereby reducing the mass removal rate in accordance with equations (3) and (4).

We model a planet which has opened a narrow gap in the disk. The migrating planet will open a gap if its mass is greater than the critical mass, $`M_cM_{}\left(81\pi \alpha c_s^2/8v_k^2\right)`$, (see, e.g., Lin & Papaloizou, 1993). The planet is initially tidally locked to the disk and moves with the gas, migrating outward or inward depending on its initial location. This migration may stop for two reasons. First, the combination of photoevaporation from the central source and viscous spreading removes disk mass, reducing the gravitational interaction between the planet and the disk and allowing the planet to decouple its orbit from that of the gas. Second, the planet will stop if it reaches the gap created by photoevaporation, as the tidal torques are unable to influence the planet across such a large opening. In order to describe planet migration, we show in figure 2 the disk stream lines which start between $`4`$ and $`18`$AU, at $`t=10^4`$yr, where there is enough mass to form a Jupiter mass planet in an annulus of width $`4`$AU. After the planet opens a gap it migrates along the disk stream lines; therefore, we can follow planet migration in figure 2 for different planet gap opening locations and times. The disk mass outside the gravitational radius is reduced to Jupiter’s mass in $`10^8`$yr (see figure 3); therefore, a migrating planet with Jupiter’s mass stops migrating at any disk location soon after this time. However, a migrating planet may reach the inner disk boundary ($`R=10^2`$AU) or the photoevaporation gap before the disk mass is reduced to Jupiter’s mass.

Figure 2 provides information on which initial locations, and times, allow a planet tidally locked to the disk to survive. It is difficult to time the gaseous disk removal in such a manner as to maintain planets at large distances from the central star. A planet which opens a narrow gap at $`t=10^4`$yr only survives if its initial disk location is greater than the critical radius $`14`$ AU, because the planet reaches the outer edge of the photoevaporation gap after the gap has formed. Since the disk expands rapidly at early times, this critical radius increases with later planet gap opening timescales; i.e. the later the planet forms, the further away from the central star it has to be in order to survive. However, this behavior changes after the photoevaporation induced gap forms. If the planet forms at a time similar to when the photoevaporation gap forms, $`t_{gap}4\times 10^7`$, it survives as long as it originates outside the gravitational radius.

## 4 Conclusions

In this letter, we have studied different scenarios in which a planet undergoing type-II migration survives due to the formation of a gap by photoevaporation from the central source. If the planet gap opening occurs close to when the photoevaporation gap forms (at $`t4\times 10^7`$yr for our fiducial model), the final semi-major axis of the planet is near the gravitational radius ($`2.4`$AU) for a wide range of initial semi-major axes. This scenario requires long time scales for planet formation ($`10^610^8`$yr), similar to those predicted by the core accretion model (Pollack et al., 1996; Bodenheimer et al., 2000).

We can also predict the final semi-major axis when considering rapid planet formation by disk instability (Boss, 2000). Our disk initially spreads outward under the influence of viscosity, allowing planets that form early at small radii to move out beyond $`R_g`$ before the photoevaporation gap opens, and enhancing their survival rate. The rate and amount of viscous spreading in real disks depends on the initial surface density profile and any continuing accretion onto the disk. If the planet induced gap opening occurs at early stages of the disk evolution (for example, at $`t10^4`$yr), all the planets initially located inside $`14`$AU will reach the disk inner boundary at the central star. The only planets which survive are those which form at a narrow range of semi-major axes ($`14R18`$AU). However, viscous diffusion spreads this narrow range during planet migration, and the final semi-major axis distribution is $`3R22`$AU. This final distribution of planets has similarities to our own Solar System and it suggests that the gas giants could have formed at $`1418`$ AU and migrated to their current locations. Furthermore, the initial planet locations are in agreement with recent simulations of planet formation by disk fragmentation (Mayer et al., 2002).

The minimum surface density for gap opening by photoevaporation, $`\mathrm{\Sigma }_{gap}`$, and the disk size, $`R_d`$, determine the disk mass outside the gravitational radius. In our fiducial model, the mass outside the gravitational radius at the time of gap opening is $`2\pi \mathrm{\Sigma }_{gap}R_g(R_dR_g)3M_J`$; therefore it is possible to consider migration of a Jupiter mass planet. However, for lower ionizing fluxes (i.e. $`\varphi _{}10^{40}\text{s}^1`$), the minimum surface density for gap opening by photoevaporation becomes smaller (see equation 5) and the halting of planet migration becomes marginal. For example, for $`\varphi _{}=10^{38}\text{s}^1`$, the minimum surface density for gap opening and the disk mass outside the gravitational radius are smaller by one order of magnitude ($`\mathrm{\Sigma }_{gap}0.2\text{gcm}^2`$ and $`M_{out}=0.3M_J`$).

We note that the gravitational radius, $`R_g`$, roughly coincides with the location of the asteroid belt in the solar system. It is well known that the surface density of heavy elements reaches a severe local minimum in the asteroid belt; this is usually attributed to long term gravitational effects from Jupiter, although it is not clear that the amount of depletion that results is sufficient (Liou & Malhotra, 1997; Holman & Murray, 1996). Simple estimates show that the planetesimal formation time scale is much shorter than the gap formation time, so it is not clear that the formation of a photoevaporation gap would significantly reduce the surface density of planetesimals, but the coincidence is intriguing.

The photoevaporation gap opening timescale ($`4\times 10^7`$yr) is long compared to the typically observed gaseous disk lifetimes ($`10^610^7`$yr). However, the disk lifetime for stars outside stellar clusters is considerably longer than the disk lifetime for stars in the hostile stellar clusters (Matsuyama et al., 2003). In addition, there is evidence that some stars retain disks between $`10^710^8`$yr (see Hillenbrand, 2002, figure 2).

It is also possible to halt both type-I and type-II migration in stellar clusters by removing the disk. This is a likely situation for the disks surrounded by ionization fronts in the Trapezium Cluster (Johnstone et al., 1998; Störzer & Hollenbach, 1999; Bally et al., 2000). Störzer & Hollenbach (1999), Hollenbach et al. (2000) and Matsuyama et al. (2003) show that the disk lifetime due to photoevaporation by external stars is in the range $`10^510^6\text{yr}`$ in the neighborhood of massive O stars, and between $`10^610^7\text{yr}`$ at reasonable distances ($`0.030.3`$pc) from external O stars. These timescales are in agreement with observational estimates for the lifetime of protoplanetary disks in stellar clusters (Haisch et al., 2001). Since the outer disk is removed first, there is no photoevaporation gap formation. This scenario is possible when considering the rapid planet formation by disk instability; on the other hand, it is uncertain whether gas giant planets form by core accretion in these short lived disks.

The research of I.M. is supported by the University of Toronto fellowship and the international recruitment award. The research of D.J. is partially supported through an NSERC grant held at the University of Victoria. The research of N.M. is supported by the Canada Research Chair Program and by NSERC.


format: text

text:
# Physical Review Letters Version of The Geometry of Developing Flame Fronts: Analysis with Pole Decomposition

## Abstract

The roughening of expanding flame fronts by the accretion of cusp-like singularities is a fascinating example of the interplay between instability, noise and nonlinear dynamics that is reminiscent of self-fractalization in Laplacian growth patterns. The nonlinear integro-differential equation that describes the dynamics of expanding flame fronts is amenable to analytic investigations using pole decomposition. This powerful technique allows the development of a satisfactory understanding of the qualitative and some quantitative aspects of the complex geometry that develops in expanding flame fronts.

The study of growing fronts in nonlinear physics  offers fascinating examples of spontaneous generation of fractal geometry. Advancing fronts rarely remain flat; usually they form either fractal objects with contorted and ramified appearance, like Laplacian growth patterns and diffusion limited agregates (DLA), or they remain graphs, but they “roughen” in the sense of producing self-affine fractals whose “width” diverges with the linear scale of the system with some characteristic exponent. The study of interface growth where the roughening is caused by the noisy environment, with either annealed or quenched noise, was a subject of active research in recent years. These studies met considerable success and there is significant analytic understanding of the nature of the universality classes that can be expected. The study of interface roughening in system in which the flat surface is inherently unstable is less developed. One interesting example that attracted attention is the Kuramoto-Sivashinsky equation  which is known to roughen in 1+1 dimensions but is claimed not to roughen in higher dimensions. Another outstanding example is Laplacian growth patterns. This Letter is motivated by a new example of the dynamics of outward propagating flames whose front wrinkles and fractalizes . We will see that this problem has many features that closely resemble Laplacian growth, including the existence of a single finger in channel growth versus tip splitting in cylindrical outward growth, extreme sensitivity to noise, etc. In the case of flame fronts the equation of motion is amenable to analytic solutions and as a result we can understand some of these issues.

The physical problem that motivates this analysis is that of pre-mixed flames which exist as self-sustaining fronts of exothermic chemical reactions in gaseous combustion. It had been known for some time that such flames are intrinsically unstable . It was reported that such flames develop characteristic structures which includes cusps, and that under usual experimental conditions the flame front accelerates as time goes on . In recent work Filyand et al.  proposed an equation of motion that is motivated by the physics and seems to capture a number of the essential features of the observations. The equation is written in cylindrical geometry and is for $`R(\theta ,t)`$ which is the modulus of the radius vector on the flame front:

$`{\displaystyle \frac{R}{t}}`$ $`=`$ $`{\displaystyle \frac{U_b}{2R_{0}^{}{}_{}{}^{2}(t)}}\left({\displaystyle \frac{R}{\theta }}\right)^2+{\displaystyle \frac{D_M}{R_{0}^{}{}_{}{}^{2}(t)}}{\displaystyle \frac{^2R}{\theta ^2}}`$ (1)
$`+`$ $`{\displaystyle \frac{\gamma U_b}{2R_0(t)}}I(R)+U_b.`$ (2)
Here 0$`<\theta <2\pi `$ is an angle and the constants $`U_b,D_M`$ and $`\gamma `$ are the front velocity for an ideal cylindrical front, the Markstein diffusivity and the thermal expansion coefficient respectively. $`R_0(t)`$ is the mean radius of the propagating flame:

$$R_0(t)=\frac{1}{2\pi }_0^{2\pi }R(\theta ,t)𝑑\theta .$$
(3)
The functional $`I(R)`$ is best represented in terms of its Fourier decomposition. Its Fourier component is $`|k|R_k`$ where $`R_k`$ is the Fourier component of $`R`$.

Numerical simulations of the type reported in ref. are presented in Fig.1. The flame front $`R(\theta ,t)`$ is shown at four equal time intervals. The two most prominent features of these simulations are the wrinkled multi-cusp appearance of the fronts and its acceleration as time progresses. One observes the phenomenon of tip splitting in which new cusps are added to the growing fronts between existing cusps. Both experiments and simulations indicate that for large times $`R_0`$ grows as a power in time

$$R_0(t)=(const+t)^\beta ,$$
(4)
with $`\beta >1`$, (of the order of $`1.5`$) and that the width of the interface $`W`$ grows with $`R_0`$ as

$$W(t)R_0(t)^\chi ,$$
(5)
with $`\chi <1`$ (of the order of 2/3). The understanding of these two features and the derivation of the scaling relation between $`\beta `$ and $`\chi `$ are the main aims of this Letter.

Equation (2) can be written as a one-parameter equation by rescaling $`R`$ and $`t`$ according to $`rRU_b/D_M`$, $`\tau tU_b^2/D_M`$. Computing the derivative of Eq.(2) with respect to $`\theta `$ and substituting the dimensionless variables one obtains:

$$\frac{u}{\tau }=\frac{u}{r_0^2}\frac{u}{\theta }+\frac{1}{r_0^2}\frac{^2u}{\theta ^2}+\frac{\gamma }{2r_0}I\{u\}.$$
(6)
where $`u\frac{r}{\theta }`$. To complete this equation we need a second one for $`r_0(t)`$, which is obtained by averaging (2) over the angles and rescaling as above. The result is

$$\frac{dr_0}{d\tau }=\frac{1}{2r_0^2}\frac{1}{2\pi }_0^{2\pi }u^2𝑑\theta +1.$$
(7)
These two equations are the basis for further analysis

Following  we expand now the solutions $`u(\theta ,\tau )`$ in poles whose position $`z_j(\tau )x_j(\tau )+iy_j(\tau )`$ in the complex plane is time dependent:

$`u(\theta ,\tau )`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}\mathrm{cot}\left[{\displaystyle \frac{\theta z_j(\tau )}{2}}\right]+c.c.`$ (8)
$`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{2\mathrm{sin}[\theta x_j(\tau )]}{\mathrm{cosh}[y_j(\tau )]\mathrm{cos}[\theta x_j(\tau )]}},`$ (9)
$$r(\theta ,\tau )=2\underset{j=1}{\overset{N}{}}\mathrm{ln}\left[\mathrm{cosh}(y_j(\tau ))\mathrm{cos}(\theta x_j(\tau ))\right]+C(\tau ).$$
(10)
In (10) $`C(\tau )`$ is a function of time. The function (10) is a superposition of quasi-cusps (i.e. cusps that are rounded at the tip). The real part of the pole position (i.e. $`x_j`$) describes the angle coordinate of the maximum of the quasi-cusp, and the imaginary part of the pole position (i.e $`y_j`$) is related the height of the quasi-cusp. As $`y_j`$ decreases (increases) the height of the cusp increases (decreases). The physical motivation for this representation of the solutions should be evident from Fig.1.

The main advantage of this representation is that the propagation and wrinkling of the front can be described now via the dynamics of the poles and of $`r_0(t)`$. Substituting (8) in (6) we derive the following ordinary differential equations for the positions of the poles:

$$r_0^2\frac{dz_j}{d\tau }=\underset{k=1,kj}{\overset{2N}{}}\mathrm{cot}\left(\frac{z_jz_k}{2}\right)+i\frac{\gamma r_0}{2}sign[Im(z_j)].$$
(11)
After substitution of (8) in (7) we get, using (11) the ordinary differential equation for $`r_0`$,

$$\frac{dr_0}{d\tau }=2\underset{k=1}{\overset{N}{}}\frac{dy_k}{d\tau }+2\left(\frac{\gamma }{2}\frac{N}{r_0}\frac{N^2}{r_0^2}\right)+1.$$
(12)
In the case of flame fronts propagating in channels of width $`L`$ ref. presented a rather complete analysis of the available stationary solutions. Some aspects of this analysis are important also for our case of cylindrical geometry, and we therefore briefly summarize the main results of . These are: (i) In noiseless conditions the total number of poles $`N_T`$ is conserved by the dynamics. This is also the case in the present problem. (ii) There is only one stable stationary solution which is geometrically represented by a giant cusp (or equivalently one finger) and analytically by $`N(L)`$ poles which are aligned on one line parallel to the imaginary axis. (iii) The reason for this behaviour is the existence of an attraction between the poles along the real line, and the resulting dynamics merges all the $`x`$ positions. The $`y`$ positions are distinct, and the poles are sitting above each others in positions $`y_{j1}<y_j<y_{j+1}`$ with the maximal $`y_{N(L)}`$. (iv) If one adds an additional pole to such a solution, this pole (or another) will be pushed to infinity along the imaginary axis. If the system has less than $`N(L)`$ poles it is unstable to the addition of poles, and any noise will drive the system towards this unique state. The number $`N(L)`$ is

$$N(L)=\left[\frac{1}{2}\left(\frac{L}{2\pi \nu }+1\right)\right],$$
(13)
where $`\left[\mathrm{}\right]`$ is the integer part, and in the (different) parametrization of ref. $`\nu `$ is the coefficient of the viscous term. (v) The height of the cusp is proportional to $`L`$. we will refer to the solution with these properties as the Thual-Frisch-Henon (TFH)-cusp solution.

In our problem the outward growth introduces important modifications to the channel results. The number of poles in a stable configuration is proportional here to the radius $`r_0`$ instead of $`L`$, but the former grows in time. The system becomes therefore unstable to the addition of new poles. If there is noise in the system that can generate new poles, they will not be pushed toward infinite $`y`$. It is important to stress that any infinitesimal noise (either numerical or experimental) is sufficient to generate new poles. These new poles do not necessarily merge their $`x`$-positions with existing cusps. Even though there is attraction along the real axis as in the channel case, there is a stretching of the distance between the poles due to the radial growth. This may counterbalance the attraction. Our first new idea is that these two opposing tendencies define a typical scale denoted as $``$. if we have a cusp that is made from the $`x`$-merging of $`N_c`$ poles on the line $`x=x_c`$ and we want to know whether a $`x`$-nearby pole with real coordinate $`x_1`$ will merge with this large cusp, the answer depends on the distance $`D=r_0|x_cx_1|`$. There is a length $`(N_c,r_0)`$ such that if $`D>(N_c,r_0)`$ then the single cusp will never merge with the larger cusp. In the opposite limit the single cusp will move towards the large cusp until their $`x`$\- position merges and the large cusp will have $`N_c+1`$ poles.

This finding stems directly from the equations of motion of the $`N_c`$ $`x`$-merged poles and the single pole at $`x_1`$. First note that from Eq.4 (which is not explained yet) it follows that asymptotically $`r_0(\tau )=(a+\tau )^\beta `$ where $`r_0(0)=a^\beta `$. Next start from 11 and write equations for the angular distance $`x=x_1x_c`$. It follows that for any configuration $`y_j`$ along the imaginary axis

$$\frac{dx}{d\tau }\frac{2N_c\mathrm{sin}x[1\mathrm{cos}x]^1}{(a+\tau )^{2\beta }}=\frac{2N_c\mathrm{cot}(\frac{x}{2})}{(a+\tau )^{2\beta }}.$$
(14)
For small x we get

$$\frac{dx}{d\tau }\frac{4N_c}{x(a+\tau )^{2\beta }}.$$
(15)
The solution of this equation is

$$x(0)^2x(\tau )^2\frac{8N_c}{2\beta 1}(a^{12\beta }(a+\tau )^{12\beta }).$$
(16)
To find $``$ we set $`x(\mathrm{})>0`$ from which we find that the angular distance will remain finite as long as

$$x(0)^2>\frac{8N_c}{2\beta 1}a^{12\beta }.$$
(17)
Since $`r_0a^\beta `$ we find the threshold angle $`x^{}`$

$$x^{}\sqrt{N}_cr_0^{\frac{(12\beta )}{2\beta }},$$
(18)
above which there is no merging between the giant cusp and the isolated pole. To find the actual distance $`(N_c,r_0)`$ we multiply the angular distance by $`r_0`$ and find

$$(N_c,r_0)r_0x^{}\sqrt{N_c}r_0^{\frac{1}{2\beta }}.$$
(19)
To understand the geometric meaning of this result we recall the features of the TFH cusp solution. Having a typical length $`L`$ the number of poles in the cusp is linear in $`L`$. Similarly, if we have in this problem two cusps a distance $`2`$ apart, the number $`N_c`$ in each of them will be of the order of $``$. From (19) it follows that

$$r_0^{\frac{1}{\beta }}.$$
(20)
For $`\beta >1`$ the circumference grows faster than $``$, and therefore at some points in time poles that appear between two large cusps would not be attracted toward either, and new cusps will appear. We will show later that the most unstable positions to the appearance of new cusps are precisely the midpoints between existing cusps. This is the mechanism for the addition of cusps in analogy with tip splitting in Laplacian growth.

We can now estimate the width of the flame front as the height of the largest cusps. Since this height is proportional to $``$ (cf. property (v) of the TFH solution), Eq.(20) and Eq.(5) lead to the scaling relation

$$\chi =1/\beta .$$
(21)
This scaling law is expected to hold all the way to $`\beta =1`$ for which the flame front does not accelerate and the size of the cusps becomes proportional to $`r_0`$.

Next we shed light on phenomenon of tip splitting that here is seen as the addition of new cusps roughly in between existing ones. We mentioned the instability toward the addition of new poles. We argue now that the tip between the cusps is most sensitive to pole creation. This can be shown in both channel and radial geometry. For example consider a TFH-giant cusp solution in which all the poles are aligned (without loss of generality) on the $`x=0`$ line. Add a new pole in the complex position $`(x_a,y_a)`$ to the existing $`N(L)`$ poles, and study its fate. It can be shown that in the limit $`y_a\mathrm{}`$ (which is the limit of a vanishing perturbation of the solution) the equation of motion is

$$\frac{dy_a}{d\tau }=\frac{2\pi \nu }{L}(2N(L)+1)1y_a\mathrm{}.$$
(22)
Since $`N(L)`$ satisfies (13) this equation can be rewritten as

$$\frac{dy_a}{d\tau }=\frac{4\pi \nu }{L}(1\alpha )y_a\mathrm{},$$
(23)
where $`\alpha =(L/(2\pi \nu )+1)/2N(L)`$. Obviously $`\alpha 1`$ and it is precisely $`1`$ only when $`L`$ is $`L=(2n+1)2\pi \nu `$. Next it can be shown that for $`y_a`$ much larger than $`y_{N(L)}`$ but not infinite the following is true:

$`{\displaystyle \frac{dy_a}{d\tau }}`$ $`>`$ $`\underset{y_a\mathrm{}}{lim}{\displaystyle \frac{dy_a}{d\tau }}x_a=0`$ (24)
$`{\displaystyle \frac{dy_a}{d\tau }}`$ $`<`$ $`\underset{y_a\mathrm{}}{lim}{\displaystyle \frac{dy_a}{d\tau }}x_a=\pi `$ (25)
We learn from these results that there exist values of $`L`$ for which a pole that is added at infinity will have marginal attraction ($`dy_a/dt=0`$). Similar understanding can be obtained from a standard stability analysis without using pole decomposition. Perturbing a TFH-cusp solution we find linear equations whose eigenvalues $`\lambda _i`$ can be obtained by standard numerical techniques. Fig.2 presents $`Re(\lambda _i)`$ as a function of systems size, and shows very clearly that (i) all $`Re(\lambda _i)`$ are non-positive. (ii) at the isolated values of $`L`$ for which $`L=(2n+1)2\pi \nu `$ $`Re(\lambda _1)`$ and $`Re(\lambda _2)`$ become zero (note that due to the logarithmic scale the zero is not evident) (iii) There exists a general tendency of all $`Re(\lambda _i)`$ to approach zero from below as $`L`$ increases. This indicates a growing sensitivity to noise when the system size increases. (iv) There exists a Goldsone mode $`\lambda _0=0`$ due to translational invariance.

The upshot of this discussion is that finite perturbations (i.e. poles at finite $`y_a`$) will grow if the $`x`$ position of the pole is sufficiently near the tip. The position $`x=\pi `$ (the tip of the finger) is the most unstable one. In the channel geometry this means that noise results in the appearance of new cusps at the tip of the fingers, but due to the attraction to the giant cusp they move toward $`x=0`$ and disappear in the giant cusp. In fact, one sees in numerical simulations a train of small cusps that move toward the giant cusp. Analysis shows that at the same time the furthest pole at $`y_{N(L)}`$ is pushed towards infinity. Also in cylindrical geometry the most sensitive position to the appearance of new cusps is right between two existing cusps independently if the system is marginal (the total number of poles fits the radius) or unstable (total number of poles is too small at a given radius). Whether or not the addition of a new pole results in tip splitting depends on their $`x`$ position. When the distance from existing cusps is larger than $``$ the new poles that are generated by noise will remain near the tip between the two cusps and will cause tip splitting.

Lastly, we note that without the noisy generation of new poles acceleration is impossible. This is seen directly from Eq.(12) in which all the terms on the RHS but unity go to zero with $`r_0\mathrm{}`$, and the velocity saturates. We need the noisy appearance of new poles to achieve acceleration. The precise connection between the noise amplitude, system size and acceleration that leads to the computation of the exponent $`\beta `$ is beyond the scope of this letter and will be discussed in a forthcoming publication.


text:
# CERN-TH/2003-027 Natural relations among physical observables in the neutrino mass matrix

## 1 Introduction

To date, 3 different kinds of experiments are being used or can be conceived to get informations on the physical parameters in the mass matrix of 3 Majorana neutrinos: i) oscillation experiments, which can make access to the 2 independent squared mass differences $`\mathrm{\Delta }m_{32}^2`$, $`\mathrm{\Delta }m_{21}^2`$ ($`|\mathrm{\Delta }m_{32}^2|>\mathrm{\Delta }m_{21}^2>0`$), to the 3 mixing angles $`\theta _{23}`$, $`\theta _{12}`$, $`\theta _{13}`$, and to a CP violating phase $`\delta `$; ii) $`\beta `$ and $`0\nu 2\beta `$ decay experiments, probing 2 appropriate combinations of masses and oscillation parameters, $`m_\beta `$ and $`m_{ee}`$; iii) cosmological and/or astrophysical experiments, which might detect, e.g., the sum of the neutrino masses $`_im_i`$. An increasing level of difficulty is involved in the different experiments, however, to reach the needed sensitivity, crucially depending on the actual neutrino mass spectrum. For a hierarchical spectrum, either “normal” ($`m_3m_2>m_1`$) or “inverted” ($`m_1m_2m_3`$), as we shall consider in the following111The degenerate spectrum is not considered since it does not satisfy the naturalness criterion as defined below., one can tentatively assume, at least for illustration purposes, that all the six oscillation observables will be measured, some of them with significant precision, perhaps together with the $`0\nu 2\beta `$ mass $`m_{ee}`$. This makes a total of seven observables, which a theory of neutrino masses should be able to correlate among each other.

A simple minded conjecture is that some of these correlations could arise from the economy in the number of independent basic parameters. As the simplest example, one may consider the possibility that the neutrino mass matrix $`M_\nu `$, in the flavour basis, has a maximum number of negligibly small entries . Inspection shows that, out of the 6 independent elements of the $`3\times 3`$ symmetric matrix $`M_\nu `$, only 2 of them could at most vanish consistently with current observations

$$\begin{array}{cc}\hfill 0.02<R\mathrm{\Delta }m_{21}^2/\mathrm{\Delta }m_{32}^2& <0.04,0.35<\mathrm{tan}^2\theta _{12}<0.55,\hfill \\ \hfill 0.6<\mathrm{tan}^2\theta _{23}<1.4,& \mathrm{sin}^2\theta _{13}0.03.\hfill \end{array}$$
(1)
This leaves a total of 15 different mass matrices to be examined, all of which dependent on 4 real parameters and 1 phase that can affect the 6 oscillation observables and $`m_{ee}`$. One expects therefore 2 relations between these observables, a priori different in every case.

Although interesting, the limit of this approach is that the elements of the neutrino mass matrix in the flavor basis are, in general, only combinations of the basic parameters in the relevant piece of the lagrangian, $`_m^\nu `$. For sure $`M_\nu `$ is generally influenced also by the charged lepton mass matrix or, if the seesaw mechanism is operative, by the mass matrix of the right handed neutrinos. A natural question to ask therefore is what happens if the economy in the number of parameters is not in $`M_\nu `$ but rather in the basic lagrangian. This is the question we address in this paper.

To answer this question in full generality would require examining a huge number of different possibilities, of which the 15 cases for $`M_\nu `$ mentioned above are only a small subset. In the following we shall only consider the possibilities that describe the currently observed pattern of the data (1) in a “natural” way, i.e. only by the relative smallness (vanishing) of some parameters in $`_m^\nu `$, barring special relations or accidental cancellations among them. This should in particular be the case when accounting simultaneously for the smallness of $`R\mathrm{\Delta }m_{21}^2/\mathrm{\Delta }m_{32}^2`$ and for the largeness of $`\theta _{23}`$, the most peculiar feature of the data so far, even though an accidental cancellation might also produce the same feature. As illustrated below, it turns out that none of the 15 cases for $`M_\nu `$ mentioned above satisfy this criterion in a strict sense. Note that naturality is not defined here in terms of any symmetry. In particular we do not require that the relative smallness (vanishing) of some parameters in the basic lagrangian be understood in terms of an explicit symmetry, exact or approximate. The coming into play of symmetries could add new cases to our list. Viceversa, it could prove hard to understand some of the cases we consider as due to a flavour symmetry of any sort.

As anticipated, we shall independently consider:
i) the direct non see-saw case (NSSC) for the 3 Majorana neutrinos $`N_L^T=(n_1,n_2,n_3)`$, described by the mass lagrangian (after electroweak symmetry breaking)

$$_m^\nu (\text{NSSC})=\overline{E_L}M_EE_R+N_L^TM_NN_L+\text{h.c.},$$
(2)
where $`E_L^T=(e,\mu ,\tau )_L`$, $`E_R^T=(e,\mu ,\tau )_R`$;
ii) the see-saw case (SSC) with the lagrangian mediated by 2 or 3 right handed neutrinos222In a 3 neutrino context, the 2 neutrino case corresponds to the limit in which one or more entries of the heavy neutrino mass matrix become much larger than the others., collectively denoted by $`N`$

$$_m^\nu (\text{SSC})=\overline{E_L}M_EE_R+N^TM_{RL}N_L+N^TM_RN+\text{h.c.}.$$
(3)
The neutrino mass matrix in the flavour basis, $`M_\nu `$, and the leptonic analog of the Cabibbo Kobayashi Maskawa matrix $`V_\nu `$ are given respectively by

$$M_\nu =U_l^TM_NU_l,V_\nu =U_l^{}U_\nu $$
(4)
in terms of the diagonal matrices

$$D_E=U_l^{}M_EV_l,D_N=U_\nu ^TM_NU_\nu ,$$
(5)
where, in the SSC, $`M_N=M_{RL}^TM_R^1M_{RL}`$. To count easily the number of effective parameters in $`M_N`$, in the SSC it is convenient to define

$$M_{RL}=mA$$
(6)
where $`m_{ij}=m_i\delta _{ij}`$ and $`(AA^{})_{ii}=1`$ with $`i=1,2`$ ($`1,2,3`$) for 2 (3) right handed neutrinos. Note that this decomposition is unique, with all $`m_i`$’s positive. It is then

$$M_N=A^T\mu ^1A\text{with}\mu ^1mM_R^1m$$
(7)
($`A`$ is adimensional, while $`\mu `$ has the dimension of an inverse mass). Of course the naturalness requirement applies equivalently to $`M_R`$, $`M_{RL}`$ or to $`\mu `$, $`A`$. The number of effective parameters in $`M_N`$ is the sum of the parameters in $`\mu `$ and $`A`$.

Within this framework, in the next section we describe all cases that i) are consistent with current data, as summarized in (1); ii) lead to some definite testable correlation between the different observables.333A case which leads to $`\mathrm{sin}\theta _{13}=0`$, $`m_{ee}=0`$ is disregarded, since we do not see how it can receive experimental support. More precisely, to realize ii) we stick to the cases that give a neutrino mass matrix in the flavor basis involving at most 4 real parameters and one phase.

## 2 Summary of results

Our results for the NSSC and for the SSC with 2 and 3 right-handed neutrinos are summarized in Table 1 and illustrated in Fig. 2 and Fig. 2. Table 1 is meant to be self explanatory. There we show the definite testable correlations between $`\theta _{13}`$, $`m_{ee}`$ and $`\theta _{23}`$, $`\theta _{12}`$, $`\delta `$ and $`R`$, that can arise from all possible cases, both in the NSSC and in the SSC, and are consistent with present data. There are only four cases (A, B, C, D) that allow to connect $`\mathrm{sin}\theta _{13}`$ and $`m_{ee}`$ with $`\theta _{23}`$, $`\theta _{12}`$ and $`R`$, whereas in case E the correlation also involves the CP-violating phase $`\delta `$. Quite a few “natural” cases have been left out from Table 1 because they are not compatible with current data at 90% confidence level: some give $`\theta _{12}45^{}`$, while two predict a too large value of $`\mathrm{sin}\theta _{13}`$, $`\mathrm{sin}\theta _{13}=\mathrm{sin}\theta _{12}R^{1/4}`$ or $`\mathrm{sin}\theta _{13}=\mathrm{tan}\theta _{12}/2(R\mathrm{cos}\theta _{12})^{1/4}`$. One case has been omitted because it is consistent, within the present uncertainty, with $`\theta _{13}=0`$ (F in the following).

The relations quoted in Table 1 are obtained through an expansion in $`R`$. The higher orders in the expansion are suppressed by $`R^{1/2}`$ in all cases except $`E`$, where the leading corrections to the relations in Table 1 are of order $`R`$. In all cases, anyhow, the exact relations can be obtained from Tables 2, 3, and 4, where the sets of parameters that originate these correlations are shown. They will be motivated and discussed in Sect. 3. All cases can be obtained in the see-saw context. Table 1 specifies the number of right-handed neutrinos involved and whether the individual cases can also originate from the NSSC. It also gives the form of the rotation on the charged lepton sector. The form of the light neutrino mass matrix before the charged lepton rotation is given in the Appendix. E is the only case that leads to an “inverted” spectrum. For cases A1, B1, E1, the independent possibility exists where $`\mathrm{tan}\theta _{23}\mathrm{cot}\theta _{23}`$, $`\mathrm{cos}\delta \mathrm{cos}\delta `$, denoted in the following by A2, B2, E2, respectively. Case A is discussed in Ref. .444For recent studies of textures in neutrino masses, also in connection with leptogenesis, see also .

Given the present knowledge of $`\theta _{23}`$, $`\theta _{12}`$ and $`\mathrm{\Delta }m_{21}^2`$, including the recent Kamland result , the ranges of values for $`\mathrm{sin}\theta _{13}`$ are shown in Fig. 2 at 90% confidence level for the different cases. It is interesting that all the ranges for $`\mathrm{sin}\theta _{13}`$, except in case D, are above $`0.02`$ and some can saturate the present limit. Long-baseline experiments of first or second generation should explore a significant portion of this range while reducing at the same time the uncertainties of the different predictions at about 10% level . Note that, in cases E, although the determination of $`\mathrm{sin}\theta _{13}`$ requires the knowledge of the CP violating phase as well, the allowed range is still limited, being $`\mathrm{sin}\theta _{13}0.10`$. Furthermore, the requirement of not exceeding the present experimental bound on $`\mathrm{sin}\theta _{13}`$ gives a lower bound on $`|\mathrm{cos}\delta |`$ (and therefore an upper limit on CP-violation) that we can quantify as

$$|\mathrm{cos}\delta |>0.8\text{at 90 \% CL }$$
(8)
given the present uncertainties. Notice that $`\mathrm{cos}\delta <0`$ ($`>0`$) in case E1 (E2). Verifying the prediction for $`\mathrm{sin}\theta _{13}`$ in case D would require the measurement of $`\theta _{23}45^{}`$; a bound on $`|1\mathrm{sin}^22\theta _{23}|`$ only sets an upper bound on $`\mathrm{sin}\theta _{13}`$, as shown in Fig. 1.555The inclusion of higher corrections in $`R`$ does not change this conclusion in any significant way.

While the prediction for $`\mathrm{sin}\theta _{13}`$ are in an experimentally interesting range, the expectations for the $`0\nu 2\beta `$-decay effective mass are mostly on the low side, except, as expected , in the only inverted hierarchical case E. The ranges for each individual cases with non-vanishing $`m_{ee}`$ are shown in Fig. 2. The challenge of detecting a non zero $`m_{ee}`$, when applicable, is therefore harder than for $`\mathrm{sin}\theta _{13}`$, with a better chance for the only inverted hierarchical case E.

## 3 Table 1 justified

In this Section we describe how we arrive to select the relatively small number of cases enumerated in Sect. 2. We build up the overall picture by commenting on the individual cases, which are summarized in Tables 2, 3 and 4. It is non trivial to make sure that we are not missing any possibility. This is obtained partly by general considerations and partly by direct inspection of the individual cases.

### 3.1 Non See-Saw Case, $`M_E`$ diagonal

As an illustration, before describing the positive cases, let us show why the NSSC with diagonal $`M_E`$ does not lead to any acceptable example, according to our rules. The neutrino mass matrix in the flavour basis, $`M_\nu `$, coincides with $`M_N`$ in eq. (2). This justifies considering the matrices $`M_\nu `$ with 2 vanishing entries also in a natural sense. The only form of $`M_\nu `$ which gives the zeroth order pattern characterized by $`R=0`$ and $`\theta _{23}`$ large without correlations between the different entries, is (PD = Pseudo-Dirac)

$$M_\nu ^{\text{PD}}=m_0\left(\begin{array}{ccc}0& c& s\\ c& 0& 0\\ s& 0& 0\end{array}\right),$$
(9)
where $`m_0`$ is the overall mass scale and $`c`$, $`s`$ stand for the cosine and the sine of a mixing angle, without loss of generality. Perturbations to (9) can lead to a $`M_\nu `$ consistent with the data with an inverted hierarchical spectrum, except for one difficulty: the angle $`\theta _{12}`$ is too close to $`45^{}`$ degrees. To get inside the range (1) requires a fine-tuning. For this reason we have not considered these cases here. To get solutions in the NSSC without fine-tuning requires a non-trivial rotation from the charged lepton sector (see Table 4).666To account for the atmospheric neutrino observations, $`s/c`$ in eq. (9) should be one within about 15%. This *is* a relation among parameters that we assume here and in the following without offering any explanation of it. We also assume that it is worth accounting for any other correlation among parameters in a natural way, as defined above.

### 3.2 See-Saw Case (2 $`N`$’s), $`M_E`$ diagonal

Again because $`\theta _{12}`$ is too close to $`45^{}`$, one does not get a natural solution with an inverted spectrum even in the SSC with a diagonal $`M_E`$. On the other hand, it is easy to obtain a normal hierarchical pattern. A $`2\times 2`$ matrix $`\mu `$ defined in eq. (7) with one dominant diagonal element and a non zero determinant may lead, in fact, to a fully natural description of the pattern (1. Taking $`\mu _{22}`$ as the dominant element, without loss of generality, one is readily convinced that the matrix $`A`$ must have the form

$$A=\left(\begin{array}{ccc}0& s& c\\ c^{}& s^{}e^{i\varphi }& 0\end{array}\right)\text{or}A=\left(\begin{array}{ccc}0& s& c\\ c^{}& 0& s^{}e^{i\varphi }\end{array}\right),$$
(10)
where we explicitly included the only phases ineliminable by a redefinition of the neutrino fields. By counting the number of real free parameters, which should not exceed 4, the only ambiguity is where one places the single small but not vanishing entry in $`\mu `$ other than $`\mu _{22}`$. There are 2 possibilities: $`\mu _{11}0`$, which leads to cases A in Table 1, and $`\mu _{12}=\mu _{21}0`$, which leads to $`\mathrm{sin}\theta _{13}=\mathrm{sin}\theta _{12}R^{1/4}`$, outside the currently allowed region. Notice that case A splits in two phenomenologically similar cases A1 and A2 according to which of the two possibilities in eq. (10) is chosen. Clearly, case A2 can be obtained from case A1 by a $`\nu _\mu \nu _\tau `$ exchange in the mixing matrix, which corresponds to $`\theta _{23}\pi /2\theta _{23}`$ in the predictions.

### 3.3 See-Saw Case (3 $`N`$’s), $`M_E`$ diagonal

In the SSC with 3 $`N`$’s and a diagonal $`M_E`$, the atmospheric angle can originate either from the heavy neutrinos Majorana mass matrix $`M_R`$ or from the Dirac mass matrix $`M_{RL}`$. In order to give rise to a large $`\theta _{23}`$, the right-handed neutrino mass matrix must correspond to a $`\mu `$ in the form

$$\mu /\mu _0=\left(\begin{array}{ccc}a& c& s\\ c& 0& 0\\ s& 0& 0\end{array}\right)+\text{1 small entry,}$$
(11)
where $`a=𝒪\left(1\right)`$ or $`a=0`$. The large $`\theta _{12}`$ angle and a non-vanishing $`\mathrm{\Delta }m_{21}^2`$ can also be obtained from (11) by adding a small entry in the “23” submatrix, provided that $`a0`$. No special structure is demanded to the Dirac matrix, which has to be diagonal in this case ($`A=\mathrm{𝟏}`$) since the four available parameters have been already used in $`\mu `$.777The texture in which the small entry in eq. (11) is in the 22 or 33 position appears also in . On the other hand, if $`a=0`$, the solar mixing angle must be provided by a non-diagonal Dirac matrix:

$$A=\left(\begin{array}{ccc}c^{}& s^{}e^{i\varphi }& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\text{or}\left(\begin{array}{ccc}c^{}& 0& s^{}e^{i\varphi }\\ 0& 1& 0\\ 0& 0& 1\end{array}\right).$$
(12)
A non-vanishing entry in the “23” submatrix of $`\mu `$ is still necessary in order to have a non zero determinant888The small determinant limit is in principle as meaningful as the large determinant limit leading to the SSC with 2 $`N`$’s. However, it does not give rise to any interesting case..

As for the cases in which the $`\theta _{23}`$ rotation arises from the Dirac mass matrix, one needs

$$\mu /\mu _0=\left(\begin{array}{ccc}a& 1& 0\\ 1& 0& 0\\ 0& 0& \sigma \end{array}\right)$$
(13)
Again, the solar angle can originate from (13) if $`a0`$999Here, as above, we have conventionally chosen a labeling for the three right-handed neutrinos that corresponds to a given order of the rows in the matrix $`A`$. and

$$A=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& s& c\end{array}\right)\text{or}\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& s& c\end{array}\right).$$
(14)
Otherwise, if $`a=0`$, the solar angle must also be provided by a matrix $`A`$ in the form

$$\left(\begin{array}{ccc}c^{}& s^{}e^{i\varphi }& 0\\ 0& 1& 0\\ 0& s& c\end{array}\right)\text{or}\left(\begin{array}{ccc}c^{}& 0& s^{}e^{i\varphi }\\ 0& 1& 0\\ 0& s& c\end{array}\right)\text{or}\left(\begin{array}{ccc}c^{}& s^{}e^{i\varphi }& 0\\ 0& 0& 1\\ 0& s& c\end{array}\right)\text{or}\left(\begin{array}{ccc}c^{}& 0& s^{}e^{i\varphi }\\ 0& 0& 1\\ 0& s& c\end{array}\right).$$
(15)
All the SSCs with 3 $`N`$’s and a diagonal $`M_E`$ are summarized in Table 3, when they lead to a case consistent with current observations. Cases B2 can be obtained from cases B1 by exchanging the last two columns of $`A`$ (and possibly a relabeling of $`N_1`$, $`N_2`$, $`N_3`$). The corresponding predictions are reported in Table 1.

### 3.4 $`M_E`$ non diagonal

With $`M_E`$ non-diagonal, the mixing matrix $`V_\nu =U_l^{}U_\nu `$ receives a non trivial contribution also from the rotation matrix in the charged sector $`U_l`$. Within our hypotheses, none of the 4 general parameters in $`U_l`$ can be correlated with the charged lepton masses. It is on the contrary possible to obtain an arbitrary $`U_l`$ with the only limitation that a large rotation in the 23 sector, $`R_{23}^l`$, if present at all, should precede all the rotations in the other sectors $`R_{12}^l`$ and $`R_{13}^l`$. Explicitly $`U_l=R_{23}^lR_{13}^lR_{12}^l`$, up to phases, is a natural form, with any of the $`R_{ij}^l`$ possibly reduced to $`\mathrm{𝟏}`$. This follows from the hierarchy of the charged lepton masses, ordered in the usual way, $`(1,2,3)=(e,\mu ,\tau )`$. Note that $`R_{12}^l`$ must be close to the identity since otherwise, when commuted with $`R_{23}^l`$ to obtain the conventional order in $`V_\nu `$, a large $`R_{13}`$ rotation would also be generated.

#### 3.4.1 $`U_l=R_{23}^l`$

$`U_l=R_{23}^l`$ (up to phases) must be considered in association with the NSSC and SSC cases for $`M_N`$. By inspection one shows that in both cases it is only in association with a hierarchical pattern of neutrino masses that interesting cases can be produced, if we disregard the further generation of cases with inverted spectrum and the $`\theta _{12}45^{}`$ problem already encountered. The only predictive and viable form of $`M_N`$ is

$$M_N=\left(\begin{array}{ccc}0& \sigma & 0\\ \sigma & 0& ϵ\\ 0& ϵ& 1\end{array}\right)$$
(16)
($`\sigma `$, $`ϵ`$ are small corrections), which leads to case C in Table 1 and 4.

In the SSC with 2 $`N`$’s, a $`R_{23}^l`$ contribution from the charged lepton sector opens the possibility

$$\mu =\mu _0\left(\begin{array}{cc}\sigma & 0\\ 0& 1\end{array}\right),A=\left(\begin{array}{ccc}0& 0& 1\\ c^{}& s^{}& 0\end{array}\right).$$
(17)
This gives a variation of case A which predicts $`\mathrm{sin}\theta _{13}=0`$ and $`m_{ee}`$ as in the other cases A, probably too small to be detected and, as such, not included in the summary of Table 1.

With 3 $`N`$’s, the only viable possibility is

$$\mu /\mu _0=\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& \sigma \end{array}\right),A=\left(\begin{array}{ccc}c^{}& 0& s^{}e^{i\varphi }\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)$$
(18)
leading to case C again. The case

$$\mu /\mu _0=\left(\begin{array}{ccc}0& \sigma & 0\\ \sigma & 0& 0\\ 0& 0& 1\end{array}\right),A=\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 1\\ s^{}& c^{}& 0\end{array}\right)$$
(19)
gives the same predictions as (17).

Since a $`\nu _\mu \nu _\tau `$ exchange can be reabsorbed in a redefinition of the $`\theta _{23}`$ mixing angle, the cases discussed in this subsection do not split in subcases. The experimental observation that $`\mathrm{tan}^2\theta _{12}`$ not only deviates from 1 but is lower than 1 serves to cut out a few otherwise acceptable cases.

#### 3.4.2 $`R_{12}^l`$ or $`R_{13}^l0`$

As already mentioned, the rotations in $`U_l`$ other than $`R_{23}^l`$, if present at all, must be small. On the other hand, they introduce new parameters. Therefore, the only way in which they can give rise to a new realistic and predictive case is when they act in association with a $`M_N`$ that depends at most on 3 effective parameters and is close enough to the fully realistic situation. On this basis, one is readily convinced that an interesting case is represented by any $`M_N`$, frequently encountered, which gives an inverted spectrum and $`\theta _{12}^\nu 45^{}`$. The role of the charged lepton contribution $`R_{12}^l`$ or $`R_{13}^l`$ is essential in providing the required correction to $`\theta _{12}=45^{}`$ in a natural way.

In this way one is led to consider mixing matrices of the form

$$V_\nu =R_{12}^l(\alpha )R_{23}\left(\begin{array}{ccc}e^{i\varphi }& & \\ & 1& \\ & & 1\end{array}\right)R_{12}(\theta _{12}^\nu )$$
(20)
where $`\alpha `$ is a free small mixing parameter. Note that the case of a $`R_{13}^l`$ contribution from the lepton sector can be obtained from eq. (20) with a $`\nu _\mu \nu _\tau `$ exchange (and a redefinition of $`R_{23}`$), which corresponds to $`\theta _{23}\pi /2\theta _{23}`$ and $`\delta \delta +\pi `$ in the predictions. To obtain $`\theta _{13}`$, it is necessary to bring $`V_\nu `$ to the standard form, i.e. to commute the $`R^l`$ rotations, either $`R_{12}^l`$ or $`R_{13}^l`$, with $`R_{23}`$. To first order in $`\alpha `$, it is

$$R_{12}^l(\alpha )R_{23}\left(\begin{array}{ccc}e^{i\varphi }& & \\ & 1& \\ & & 1\end{array}\right)R_{12}(\theta _{12}^\nu )=R_{23}R_{13}(\alpha s_{23})\left(\begin{array}{ccc}e^{i\varphi }& & \\ & 1& \\ & & 1\end{array}\right)R_{12}(\theta _{12}^\nu +\alpha c_{23}\mathrm{cos}\varphi )$$
(21)
which leads to case E1 in Table 1 and 3 (where the relations are exact in $`\alpha `$). The corresponding case E2 is obtained for $`U_l=R_{13}^l`$.

The explicit realizations of the mechanism are obtained in the NSSC case by adding a single correction to $`M_N`$ in eq. (9). In the SSC case with 2 $`N`$’s, they correspond to adding a single correction to $`\mu /\mu _0`$ or $`A`$ in

$$\mu =\mu _0\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),A=\left(\begin{array}{ccc}1& 0& 0\\ 0& c& s\end{array}\right),$$
(22)
while in the SCC with 3 $`N`$’s they correspond to

$$\mu /\mu _0=\left(\begin{array}{ccc}0& \sigma & 0\\ \sigma & 0& 0\\ 0& 0& 1\end{array}\right)\text{and}A=\left(\begin{array}{ccc}0& c& s\\ 1& 0& 0\\ 1& 0& 0\end{array}\right),\left(\begin{array}{ccc}0& c& s\\ 1& 0& 0\\ 0& 1& 0\end{array}\right)\text{, or}\left(\begin{array}{ccc}0& c& s\\ 1& 0& 0\\ 0& 0& 1\end{array}\right).$$
(23)
Starting from similar configurations, $`U_l=R_{12}^lR_{23}^l`$ and $`U_l=R_{13}^lR_{23}^l`$ lead to the same cases as $`U_l=R_{12}^l`$ and $`U_l=R_{13}^l`$.

The leptonic rotations $`R_{12}^l`$ and $`R_{13}^l`$ can also play a role in hierarchical cases in which the neutrino contribution to $`V_\nu `$ gives $`\theta _{12}>45^{}`$. The charged lepton contribution is then used to bring $`\theta _{12}`$ back to the allowed range below $`45^{}`$. This possibility is only realized in the SSC with 3 $`N`$’s for

$$\mu /\mu _0=\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& \sigma \end{array}\right)\text{and}A=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& s& c\end{array}\right),\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& s& c\end{array}\right),$$
(24)
giving respectively cases F1a, F1b for $`U_l=R_{12}`$ and cases F2a, F2b for $`U_l=R_{13}`$; or adding a single correction in the 2-3 submatrix of $`\mu /\mu _0`$ in

$$\mu /\mu _0=\left(\begin{array}{ccc}0& c& s\\ c& 0& 0\\ s& 0& 0\end{array}\right),A=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),$$
(25)
giving cases F1a, F1b, F1c for $`U_l=R_{12}`$ and cases F2a, F2b or F2c for $`U_l=R_{13}`$. As anticipated in Section 2, cases F do not appear in Table 1 because their predictions for $`\mathrm{sin}\theta _{13}`$ are either inconsistent with present data (Fa, Fc) or not conclusive, since it is at present compatible with zero (Fb). To get to this conclusion it is necessary to take into account not only the leading order prediction for $`\mathrm{sin}\theta _{13}`$, which is common to the three cases F1a, F1b, F1c,

$$\mathrm{cos}\delta \mathrm{sin}\theta _{13}^{(0)}=\mathrm{tan}\theta _{23}\frac{1\mathrm{tan}\theta _{12}}{1+\mathrm{tan}\theta _{12}}$$
(26)
but to include as well higher order corrections in $`R`$:

$$\begin{array}{cc}\hfill \mathrm{cos}\delta \mathrm{sin}\theta _{13}& =\mathrm{tan}\theta _{23}\frac{t\mathrm{tan}\theta _{12}}{1+t\mathrm{tan}\theta _{12}}+q\left(\frac{R}{2}\right)^{1/3}\hfill \\ \hfill t& =1+\frac{p}{2}\left(\frac{R}{2}\right)^{1/3},\text{with}\hfill \\ \hfill p& =(\mathrm{tan}\theta _{23})^{4/3}\text{}(\mathrm{tan}\theta _{23})^{4/3}\text{}1\hfill \\ \hfill q& =(\mathrm{tan}\theta _{23})^{1/3}\text{}(\mathrm{tan}\theta _{23})^{1/3}\text{}(\mathrm{tan}2\theta _{23})^1\hfill \end{array}$$
(27)
in cases F1a, F1b, F1c respectively.

## 4 Conclusions

The economy in the number of basic parameters could be at the origin of some correlations between the physical observables in the neutrino mass matrix. At the present state of knowledge, the variety of the possibilities for the basic parameters themselves is large: $`M_N`$, $`M_E`$, $`M_R`$, $`M_{RL}`$ are the matrices that might be involved. Finding the minimal cases that describe the present pattern of the data in a natural way could be a first step in the direction of discriminating the relevant $`_m^\nu `$. This we have done with the results summarized in Tables 14 and illustrated in Fig. 2 and Fig. 2. It is remarkable that the number of possible correlations between the physical observables is limited (Table 1), with a relatively larger number of possibilities for the basic parameters (Tables 24). The relatively best chance for selecting experimentally one out of the few relevant cases is offered by $`\mathrm{sin}\theta _{13}`$. Its predictions, in Table 1, should have a 10% uncertainty with the improved determination of the other parameters foreseen in long-baseline neutrino experiments . Combining this with independent studies of leptogenesis or of lepton flavour violating effects could lead to the emergence of an overall coherent picture.

We have insisted on “naturalness” both in solving the “large $`\theta _{23}`$-small $`R`$” problem for the normal hierarchy case and in obtaining a significant deviation of $`\theta _{12}`$ from $`45^{}`$, with small $`R`$, in the inverted hierarchy case. Both features might be due to an accidental tuning of parameters. Nevertheless, explaining these features in a natural way offers a possible interesting guidance for model building. Note, in this respect, the use in cases E1, E2 of a relatively small charged lepton rotation in the 12 or 13 sectors to solve the second fine tuning problem mentioned above. This is analogous to the use of a large charged lepton rotation in the 23 sector to account for the “large $`\theta _{23}`$-small $`R`$” problem in a natural way .

## Acknowledgments

This work has been partially supported by MIUR and by the EU under TMR contract HPRN–CT–2000–00148. Part of the work of A.R. was done while at the Scuola Normale Superiore. We thank A. Strumia for useful comments.

## Appendix A General forms of the $`M_N`$ matrices for the various cases


text:
# Submitted to Physical Review Letters, February 2003. Structure and bonding of dense liquid oxygen from first principles simulations

## Abstract

Using first principles simulations we have investigated the structural and bonding properties of dense fluid oxygen up to 180 GPa. We have found that band gap closure occurs in the molecular liquid, with a “slow” transition from a semi-conducting to a poor metallic state occurring over a wide pressure range. At approximately 80 GPa, molecular dissociation is observed in the metallic fluid. Spin fluctuations play a key role in determining the electronic structure of the low pressure fluid, while they are suppressed at high pressure.

The study of pressure-induced transformations in solids and liquids has become an active field of research in the last decade GeneralReview , due to key progress both in experimental methods (diamond anvil cell (DAC) and shock wave experiments) and in accurate computational tools. In particular, molecular systems such as hydrogen Review\_H , oxygen Review\_O  and nitrogen Review\_N  have been investigated with a variety of techniques. Nevertheless, fundamental questions on the nature of these systems in e.g. the dense liquid state are yet unanswered.

In this Letter, we focus on the liquid state of oxygen. In particular, we examine the structural and bonding changes induced by pressure, as well as the mechanism leading to the insulator-to-metal transition (ITMT) observed experimentally in shock wave experiments Bastea01 .

Metallization of oxygen was first observed in the solid state by Desgreniers et al. Desgreniers  who reported a Drude-type metallic behavior of its reflectivity for pressures (P) above 96 GPa. These results were later confirmed by Shimizu and coworkers shimizu95  who observed a change in the sign of the resistance-temperature curve at 96 GPa, and also discovered superconductivity in solid, dense oxygen shimizu98 . Following measurements by Akahama et al. Akahama , Weck et al. Weck02  showed that in solid oxygen the ITMT is accompanied by a structural change between two different crystallographic arrangements of molecules. The insulating phase below 96 GPa was also studied with infrared absorption spectroscopy by Gorelli et al. Gorelli  who suggested the existence of a lattice of O4 molecules at high pressure. However, this hypothesis is not confirmed by recent ab initio calculations NA02 , showing that a geometrical configuration built from extended “herringbone”-type chains is energetically favored over a lattice of O4 molecules. The structure predicted theoretically is consistent with earlier suggestions based on DAC experiments Agnew87 .

In fluid oxygen, metallization was recently reported by Bastea et al. Bastea01  who used reverberation shock wave techniques to generate compression states on a quasi-isentrope. For a given shock speed, the thermodynamic states reached in these experiments are at higher pressure and lower temperatures ($`T`$) than Hugoniot states obtained in earlier shock wave measurements Ha88 . When increasing P from $``$ 30 to 120 GPa (and correspondingly $`T`$ by several thousands K), Bastea et al. Bastea01  observed an increase of the electrical conductivity over six orders of magnitude. At about 120 GPa and 3.4-fold compression, the measured conductivity attains values characteristic of a poor metal. As $`P`$ is further increased from 120 to 190 GPa, the conductivity levels off and is nearly constant as a function of $`P`$. These experimental results, which indicate that the ITMT occurs at higher P in the liquid than in the solid, are rather surprising. In general, one would expect disorder or dissociation effects to lead to a lower metallization pressure in the fluid than in the ordered solid.

In order to investigate the mechanism by which the ITMT occurs in the fluid, we have carried out a series of ab initio molecular dynamics (MD) simulations CP85  in the liquid state for different pressures. Our findings shed light on the microscopic structure and electronic properties of the dense fluid state and on the key role of spin fluctuations present in fluid before metallization occurs. Our predicted metallization pressure in the fluid is lower than that reported experimentally Bastea01 , as well as lower than the metallization pressure measured in the solid. Our findings further indicate that the system becomes metallic in a predominantly molecular state, as opposed to other fluids like hydrogen HG02  and nitrogen mazevet03 .

We have carried out ab initio MD simulations CP85  using a gradient corrected energy functional PBE  in the density ($`\rho `$) and temperature range explored in recent shock wave experiments: $`2.8\rho 4.5\mathrm{g}\mathrm{cm}^3`$ and $`1000T6400\mathrm{K}`$. In some of our simulations we have treated the spin variable explicitly within the gradient corrected local spin density approximation (GCSDA). We performed constant volume microcanonical simulations Jeep  with cells containing 54 and 108 oxygen atoms. For each $`\rho `$ and $`T`$ the total simulation time was 4–5 ps, after equilibration.

Our calculated pair correlations functions ($`g(r)`$) are shown in Fig. 1. Finite size effects were found to be negligible compared to temperature and pressure effects. At $`\rho =2.8\mathrm{g}\mathrm{cm}^3`$ and $`T=1000\mathrm{K}`$ weck\_melting , the first sharp peak of $`g(r)`$ at 1.20 Å, followed by a minimum very close to zero, is indicative of a molecular fluid. The computed pressure is 24 GPa. The first maximum is at a distance very close to the bond length of an oxygen molecule in the gas phase in its triplet ground state (1.207 Å) PBE\_molecules . As the temperature is increased to 3000$`\mathrm{K}`$ (and correspondingly the pressure to 34 GPa) the first molecular peak is broadened while its position remains unchanged. This trend continues as the fluid is compressed to $`\rho `$ = 3.6$`\mathrm{g}\mathrm{cm}^3`$ at T=3000$`\mathrm{K}`$ ($`P`$=84 GPa). Upon further compression to $`\rho `$ = 4.5$`\mathrm{g}\mathrm{cm}^3`$ at T=6400$`\mathrm{K}`$ (P=179 GPa), the first molecular peak becomes less intense and the second broad peak moves to shorter distances. It is interesting to note that the position of the intramolecular peak remains constant over a pressure range of 150 GPa and does not shift to larger separations as e.g. in dense liquid hydrogen Ga99 .

The ground state of the O2 molecule is a triplet, and this has been shown to give rise to interesting structural effects in the low density fluid Pasquarello02 . It is therefore interesting to explore the effect of explicit spin treatment on the structural properties of the high pressure fluid as well. We found that spin effects are negligible on average structural properties both at low and high pressure, consistent with recent ab initio calculations for solid oxygen NA02  in the same pressure range. For example in Fig. 1, $`g(r)`$ calculated with and without explicit treatment of spin are found to be indistinguishable. However, spin effects are key in determining the optical excitation spectrum of the low-P fluid, as we will discuss below when we present our results for the electronic conductivity.

The microscopic structure of the liquid can be further examined by a cluster analysis. We define a set of atoms as belonging to a given cluster if these atoms are separated by less than a given cut-off distance $`r_c`$. (We chose $`r_c=1.50`$ Å, which corresponds to the first $`g(r)`$ minimum in the molecular regime.) We then define cluster lifetimes by simply computing the duration for which a given cluster remains intact. This analysis yields a very simple picture for our simulations at low density (2.8$`\mathrm{g}\mathrm{cm}^3`$) and 1000 and 3000$`\mathrm{K}`$: we find that the liquid consists entirely of O2 molecules. A density increase to 3.6$`\mathrm{g}\mathrm{cm}^3`$ at 3000$`\mathrm{K}`$ leads to a small fraction of dissociated molecules. Although the liquid is still predominantly molecular (94% O2), we find 1% single atoms and 5% O3-like clusters. At this density, the average lifetime of O2 molecules is about 15 vibrational periods of the isolated O2 molecule (21.1 fs), while O1 and O3 exist only for about one molecular vibration. When the temperature is raised to 5000$`\mathrm{K}`$ at 3.6$`\mathrm{g}\mathrm{cm}^3`$, the number of molecules decreases to 79%, leading to 7% atoms, 10% O3, and 4% O4-like complexes. The lifetime of molecules decreases to about three O2 vibrations, indicating the presence of a dynamical dissociation equilibrium with molecules still being the predominant and most stable species. At $`\rho `$=4.5$`\mathrm{g}\mathrm{cm}^3`$ and T=6400$`\mathrm{K}`$, we find 9% atoms and 58% molecules and also 16% O3 clusters. The remaining 17% of the nuclei belong to short-lived complexes composed of four or more atoms. Further inspection reveals that O4-like clusters arrange predominantly into short zig-zag chains, consistent with the finding of Ref. NA02  for the compressed solid. At this high density, the O2 lifetime is only about 1.5 vibrational periods. The short lifetime of O3 at all densities clearly indicates that triatomic configurations are closer to scattering states than to stable chemical species. Nevertheless, it is interesting to note that the “bond” angles of O3 clusters are not randomly distributed but vary between 80 and 150 with an average of 120. Therefore, the transient structure of O3 clusters resembles that of the ozone molecule.

As mentioned above, while spin effects are negligible for the structural properties of the fluid, it is crucial to include them explicitly in order to understand the electronic properties of liquid oxygen. In the low density regime, the liquid is composed of a disordered arrangement of O2 molecules carrying a net spin of one as in the gas phase. The left panel of Fig. 2 shows the difference ($`\mathrm{\Delta }\rho _{}`$) between spin-up ($`\rho _{}`$) and spin-down densities ($`\rho _{}`$), for a snapshot at $`\rho =2.8\mathrm{g}\mathrm{cm}^3`$ and 1000$`\mathrm{K}`$, showing a random distribution of up and down spins. At higher P (right panel of Fig. 2), $`\mathrm{\Delta }\rho _{}`$ has almost vanished throughout the whole system, indicating that spin fluctuations are suppressed. The effect shown qualitatively in Fig. 2 can be characterized quantitatively by the correlation function of the nuclei $`𝐫_i`$ and the electronic magnetic moment, $`m_z(𝐫)=\mu _B[\rho _{}(𝐫)\rho _{}(𝐫)]`$,

$$\rho _eg_{\mathrm{O},m_z}(|𝐫|)=\frac{1}{N}\underset{i=1}{\overset{N}{}}m_z(𝐫_i𝐫),$$
(1)
where $`\rho _e`$ is the electron density, $`N`$ the number of atoms and $`\mathrm{}`$ denotes an average over ionic configurations and the solid angle of $`𝐫`$. This correlation function is similar to the nuclei-electron pair correlation function. The first, intense peak in the low density correlation function shown in Fig. 3 indicates that the majority of molecules are in the triplet state. At higher density, the intensity of the first peak is decreased by almost an order of magnitude, showing that most of the pairs in the fluid do not carry a net spin.

In order to qualitatively understand the effect of pressure on the spin state of the oxygen molecule, we carried out calculations of the total energies of O2 in the singlet ($`E_{S=0}`$) and triplet ($`E_{S=1}`$) state in the presence of a confining potential mimicking an average pressure effect. Interestingly, we found that the difference $`E_{S=0}E_{S=1}`$ can be substantially reduced by applying an external confining potential. This difference decreases as a function of the confinement strength.

In order to compare our results with recent shock wave data, we computed the conductivity ($`\sigma `$) of the liquid for selected snapshots taken from our GCSDA simulations at different pressures, using the Kubo-Greenwood formula Si99 ; collins01 . We note that in the low pressure regime, where spin fluctuations are present, GCSDA and GGA results are very different. In particular, the GGA yields much smaller optical gaps, consistent with the incorrect treatment of excitation in the isolated molecule if spin is not explicitly taken into account. On the other hand, at higher densities, the computed GGA and GCSDA conductivities are very similar.

Our results are reported in Fig. 4 along with the experimental results of Ref. Bastea01  and Ha88 . The comparison between theory and experiment is not straightforward: the temperature in the measurements by of Hamilton et al. Ha88  is in general different from that of Bastea et al.’s data Bastea01 , and in both experiments, T was not measured but estimated on the basis of empirical equation of state (EOS) models. However, the two sets of experimental data show the same qualitative rise of $`\sigma `$ as a function of P and T. In particular, Bastea et al.’s data show a six order of magnitude rise in conductivity with increasing shock compression, corresponding to a P and T increase from about 30 to 120 GPa and $`1000`$ to $`5000\mathrm{K}`$, respectively. Our calculations within GCSDA show only two orders of magnitude increase in $`\sigma `$ in the same P and T range. This increase of $`\sigma `$ as a function of $`T`$ occurs primarily in a state where the fluid is entirely composed of stable O2 molecules ($`\rho =2.8\mathrm{g}\mathrm{cm}^3,\mathrm{T}=10003000\mathrm{K}`$). This indicates that in oxygen, metallization is not accompanied by molecular dissociation as found in other first row molecular fluids, in particular hydrogen collins01  and nitrogen mazevet03 . Above 100 GPa both theory and experiment show an almost constant $`\sigma `$ as a function of compression, but the theoretical values are about a factor of 7 larger than the experimental ones.

In our calculations, overestimates of conductivity values could be due to the use of the GCSDA, which is known to underestimate electronic gaps in many systems. In order to qualitatively investigate this effect, we used a scissors operator and arbitrarily increased the value of the computed gaps by a factor of two. This lowered the computed conductivity by one order of magnitude at $`\rho `$ = 2.8$`\mathrm{g}\mathrm{cm}^3`$, both at 1000 and 3000 K. However, even after this “gap correction”, our theoretical data show values of $`\sigma `$ approaching those of a poor metal at a pressure much lower (30–40 GPa) than found experimentally. The computed conductivity is then slowly changing as a function of P, and it reaches metallic-like values at approximately 80 GPa. This “slow” metallization of the fluid as a function of pressure may be related to spin fluctuations and their “slow” suppression as a function of pressure. At low P, the molecular fluid, with each oxygen molecule in a triplet state, can be viewed as a disordered spin system. As the pressure is increased, this system shows a glass-like behavior due to the presence of a multitude of local minima in its potential energy surface, corresponding to a multitude of different spin arrangements. Eventually, the fluid reaches metallization when spin fluctuations are almost completely suppressed. In this respect, compressed liquid oxygen is unique among first row molecular fluids and rather different from, e.g. liquid hydrogen and liquid nitrogen, composed of molecules with spin zero.

In summary, our first principles simulations show that above 30 GPa, over a wide pressure range of about 50 GPa, liquid oxygen transforms from a semi-conducting to a metallic fluid by closure of the band gap in a disordered molecular state. For $`P`$ 80 GPa, molecular dissociation is observed, giving rise to short-lived atomic- and ozone-like species. We did not observe any stable O4 molecule for the thermodynamic conditions considered here, which is consistent with recent theoretical findings for solid oxygen NA02 . Our results also show that spin fluctuations are suppressed at high pressure, with no significant magnetic moment present in the high pressure liquid, as opposed to the semi-conducting low pressure regime.

Our findings for the electrical conductivity point at the need for further experimental and theoretical work. Experimentally, in shock wave experiments only P and $`\sigma `$ are measured. Data are then analyzed using a model EOS as input for hydrodynamic codes which compute densities. Finally, the temperature is inferred from additional empirical chemical model calculations. This type of analysis may introduce uncertainties in the estimated $`\rho `$ and T at a given P and $`\sigma `$. Measurements of the oxygen EOS obtained with shock wave experiments and possibly with other techniques, e.g. DAC measurements, may help clarify the existing discrepancies between theory and experiment. Finally, we note that theoretical inaccuracies on calculated conductivity values may arise from the use of an approximate exchange and correlation functional. Although computationally still out of reach, quantum Monte Carlo calculations for liquid oxygen under pressure might be used in the future to address this issue HG02 .

###### Acknowledgements.

We thank M. Bastea for her support of this project and for many useful discussions, D. Young for providing us with EOS tables from chemical models, and L. Benedict, V. Recoules, E. Schwegler, and T. Ogitsu for useful suggestions. This work was performed under the auspices of the U.S. Dept. of Energy at the University of California/LLNL under contract no. W-7405-Eng-48.


text:
# Amorphous vortex glass phase in strongly disordered superconductors

(December 18, 2002)

## Abstract

We introduce a model describing vortices in strongly disordered three-dimensional superconductors. The model focuses on the topological defects, i.e., dislocation lines, in an elastic description of the vortex lattice. The model is studied using Monte Carlo simulations, revealing a glass phase at low temperatures, separated by a continuous phase transition to the high temperature resistive vortex liquid phase. The critical exponents $`\nu 1.3`$ and $`\eta 0.4`$ characterizing the transition are obtained from finite size scaling.

Disorder has a profound impact on the phase diagram of type-II superconductors. The long range order of the Abrikosov vortex lattice in very clean samples is destabilized upon the introduction of disorder Larkin  and a glass state can appear instead M.P.A.Fisher ; FFH . The effects of weak disorder are relatively well understood and may be treated within an elastic theory of the vortex lattice Feigelman-89-Nattermann ; Giamarchi-Doussal ; D.Fisher . The resulting phase, known as the Bragg glass or the elastic vortex glass, is characterized by an algebraic decay of translational correlations and perfect orientational order, and the absence of large dislocation loops D.Fisher .

For stronger disorder dislocations start to become energetically favorable and begin to proliferate. What happens then is still a matter of controversy. One possibility is that any kind of order in the system is lost and the system becomes thermodynamically equivalent to a vortex line liquid (although perhaps a very viscous one) with a non-zero resistivity FFH . The other possibility is that the dislocations themselves remain pinned and the system freezes into a different genuine thermodynamic amorphous glass phase, full of dislocation loops, and with zero linear resistivity. Such a disorder driven transition from the Bragg glass to the vortex glass have been proposed to occur as the magnetic field is increased Giamarchi-Doussal ; Gingras-Huse .

Unfortunately, the theoretical methods based on pinning of elastic manifolds in random media, that were so successful for the Bragg glass, does not easily extend to incorporate the effects of dislocations. In this Letter we introduce and study a model that goes beyond a purely elastic description by including dislocations, and is well suited for numerical simulations.

Experimentally, the evidence for a vortex glass comes mainly from the scaling relations obeyed by the current-voltage characteristics and resistivity at a continuous phase transition FFH . The first order vortex melting line, observed in very clean samples at low magnetic fields, is replaced by a smooth continuous transition upon increasing the field or disorder. The universal critical exponents extracted from measurements in this regime provide an important characterization of the glass transition that can be compared with theoretical predictions. A large number of studies have thus been able to obtain good scaling behavior in different materials, and this fact has been taken to support the existence of an amorphous glass phase vortexglass-experiment . This interpretation has, however, been questioned recently Strachan . In practice it may be extremely difficult to discriminate between a true thermodynamic glass and a very slowly moving non-equilibrium viscous vortex liquid (as window glass), with a small, but non-zero, resistivity due to thermally activated flux creep.

On the theoretical side a lot of interest has focused on the so called gauge glass model gauge-glass . This approach, which is inspired by spin glasses in random magnets, combines the $`U(1)`$ symmetry of a superconducting order parameter with quenched disorder and frustration, which are believed necessary to create a glass phase, but has no other connection with the microscopics. The model, which consists of an $`XY`$-model in a quenched random gauge field, has a glass phase at low temperatures in three dimensions gauge-glass  (but probably not in two). Once screening of the vortex interaction, due to gauge field fluctuations, is accounted for, the glass phase disappears gauge-glass-screening . However, the way disorder enters the model, via quenched random fluxes instead of a random core energy coupling to the vortices, is not very realistic. More importantly, perhaps, the model lacks a length scale, the inter-vortex distance, set by the external magnetic field.

Attempts to study more realistic models have come to contradictory conclusions. Simulations of frustrated 3D $`XY`$-models with quenched disorder show evidence for a direct disorder driven transition from the Bragg glass to a vortex liquid at low disorder and fillings, and find no clear signs of any amorphous glass phase Olsson-Teitel  (See, however, Ref. Kawamura, ). On the other hand recent simulations Vestergren  of a random pinning vortex model for much stronger disorder and very high field, corresponding to a half-filled system, do find a continuous transition, but with critical exponents quite far from experiments. None of these models included screening. Including screening in a similar model appears to destroy the glass phase Kawamura , just as in the gauge glass. Langevin simulations of vortex lines also found no glass phase molasses . Recently, Kierfeld et al. Kierfeld  invoked a Landau theory in terms of the dislocation density, assuming the dislocation lines were directed, and proposed a phase diagram containing an amorphous glass phase. Some support for a glass phase has also been found using the boson analogy Kihong .

The approach taken in this Letter starts from an elastic theory of the vortex lattice, but modified to include dislocations. The basic degrees of freedom in our model is the two-dimensional displacement vector $`𝐮_i`$ and a tensor field $`𝐁`$ that accounts for the dislocations. Since the clean Abrikosov vortex lattice usually has hexagonal symmetry we define the model on a stacked triangular lattice with $`N=L^3`$ sites and unit lattice constant, and with the displacements $`𝐮_i`$ living on the vertices and the field $`𝐁_{ij}`$ living on the links connecting nearest neighbors $`i`$ and $`j`$. The Hamiltonian reads

$$=\frac{1}{2}\underset{ij}{}\left(𝐰_{ij}𝐅_{ij}\right)^2,$$
(1)
where $`ij`$ denote all nearest neighbor pairs on the stacked triangular lattice, $`𝐰_{ij}=𝐮_i𝐮_j+𝐁_{ij}`$, and $`𝐮_i^2`$. The $`𝐁_{ij}`$ are integer multiples of the elementary triangular lattice vectors $`𝐚_1=(1,0)`$, $`𝐚_2=(\frac{1}{2},\frac{\sqrt{3}}{2})`$ in the plane, i.e., $`𝐁_{ij}=m𝐚_1+n𝐚_2`$, $`m,n`$, and contain all the information about the dislocations. More precisely the dislocation density is given by the lattice curl of $`𝐁`$, implying that dislocations form lines with conserved Burgers vector as they should. A dual representation in terms of interacting dislocation lines can be obtained by integrating out $`𝐮`$, but Eq. (1) turns out to be more convenient for our purposes. The coupling to disorder occurs via the random variable $`𝐅_{ij}`$ 111This corresponds to a random stress. Other ways of incorporating disorder would also be interesting to consider.. For simplicity we take the maximum possible disorder strength, with $`𝐅_{ij}`$ independently and uniformly distributed in the planar unit cell of the triangular lattice.

In the absence of dislocations all $`𝐁_{ij}`$ are zero and the model reduces to an ordinary elastic model of the vortex lattice with compression modulus $`c_{11}=3/2`$, tilt modulus $`c_{44}=1`$, and shear modulus $`c_{66}=3/2`$ 222 To avoid factors of $`\sqrt{3}/2`$ extensive quantities are normalized using $`N`$ instead of the volume $`\mathrm{\Omega }=N\sqrt{3}/2`$. . The generalization to other values of these (bare) elastic constants is straightforward via the addition of a term proportional to $`_{ij}\left(𝐰_{ij}𝐞_{ij}\right)^2`$ \[$`𝐞_{ij}`$ being the unit vector connecting nearest neighbor sites $`i`$ and $`j`$\] to Eq. (1) and/or by changing the coupling in the $`z`$-direction to make it anisotropic. In this work we will, however, restrict to the isotropic model defined above.

The presence of thermally excited dislocations will renormalize the elastic constants and eventually cause a transition at which they drop to zero. The renormalized elastic constants can be obtained by taking the second derivative of the free energy with respect to a deformation, $`𝐰_{ij}𝐰_{ij}+\epsilon _{\mu \nu }(𝐞_{ij}𝐞_\mu )𝐞_\nu `$. This gives

$`C_{\mu \nu |\rho \sigma }`$ $`=`$ $`{\displaystyle \frac{^2f}{\epsilon _{\mu \nu }\epsilon _{\rho \sigma }}}={\displaystyle \frac{1}{2}}\delta _{\mu \rho }\delta _{\nu \sigma }\left(3\delta _{\mu ,\widehat{z}}\right)`$ (2)
$`{\displaystyle \frac{N}{T}}\left(\sigma _{\mu \nu }\sigma _{\rho \sigma }\sigma _{\mu \nu }\sigma _{\rho \sigma }\right),`$
where $`\mathrm{}`$ denotes a thermal average, and the $`3\times 2`$ stress tensor $`\sigma `$ is given by

$$\sigma _{\mu \nu }=\frac{1}{N}\underset{ij}{}\left(𝐞_{ij}𝐞_\mu \right)\left(𝐞_\nu \left[𝐁_{ij}𝐅_{ij}\right]\right).$$
(3)
Unfortunately, both of these vanish upon averaging over the disorder realizations, $`\left[\sigma \right]_{\mathrm{av}}=\left[C\right]_{\mathrm{av}}=0`$, since $`\left[f\right]_{\mathrm{av}}`$ becomes independent of $`\epsilon _{\mu \nu }`$ for the fully disordered model, where the averaging includes all possible deformations. A non-zero $`\sigma `$ is still indicative of a phase with rigidity, and hence we use

$$\sigma _{\mathrm{rms}}^2=\left[\underset{\mu \nu }{}\sigma _{\mu \nu }^2\right]_{\mathrm{av}}=\left[\underset{\mu \nu }{}\sigma _{\mu \nu }_\alpha \sigma _{\mu \nu }_\beta \right]_{\mathrm{av}}$$
(4)
instead, where the second form, which involves two different replicas $`\alpha `$ and $`\beta `$ to avoid any bias in the expectation value, is the one used in the actual simulations. Alternatively, $`\sigma _{\mathrm{rms}}^2`$ can be calculated from Eq. (2) using that $`\left[C\right]_{\mathrm{av}}=0`$,

$$\sigma _{\mathrm{rms}}^2=\left[\underset{\mu \nu }{}\sigma _{\mu \nu }^2\right]_{\mathrm{av}}\frac{8T}{N}.$$
(5)
In the solid low temperature phase (if it exists) $`\sigma _{\mathrm{rms}}^2`$ will acquire a non-zero expectation value, while in the liquid it approaches zero. Close to a continuous phase transition it obeys the finite size scaling form,

$$\sigma _{\mathrm{rms}}=L^{1d}\widehat{\sigma }(tL^{1/\nu }),$$
(6)
where $`t=(TT_c)/T_c`$ is the reduced temperature, $`\nu `$ the correlation length exponent, $`d=3`$ the dimensionality, and $`\widehat{\sigma }`$ is a scaling function 333 The absence of any nontrivial anomalous dimension in Eq. (6) is a result of the underlying periodicity in the model. At a fixed point with such periodicity the scaling dimension of $`𝐮`$ must vanish, and hence $`𝐰`$ and $`𝐁`$ scales as $`L^1`$ and $`\sigma `$ as $`L^{1d}`$..

Lattice order may be probed by the structure function $`S_𝐊=\frac{1}{N}\left|_ie^{i𝐊𝐮_i}\right|^2`$ at a reciprocal lattice vector $`𝐊`$ of the triangular lattice, but for strong disorder, when the system might freeze into a highly disordered glassy state, this is not useful. To quantify order in a glass one may instead consider an overlap order parameter,

$$q_𝐊=\frac{1}{N}\underset{i}{}e^{i𝐊\left(𝐮_i^\alpha 𝐮_i^\beta \right)},$$
(7)
between two replicas $`\alpha `$ and $`\beta `$. The order parameter susceptibility, given by

$$\chi _𝐊=N\left[\left|q_𝐊\right|^2\right]_{\mathrm{av}},$$
(8)
obeys the finite size scaling relation

$$\chi _𝐊=L^{2\eta }\widehat{\chi }(tL^{1/\nu }),$$
(9)
which defines the critical exponent $`\eta `$.

We study this model via Monte Carlo (MC) simulations. One MC sweep consists of $`N`$ randomly chosen heat bath updates of the $`𝐮_i`$ and $`12N`$ Metropolis updates of $`𝐁_{ij}`$, where $`N=L^3`$ is the number of sites of the lattice. In order to overcome the slow glassy dynamics in the fully disordered system an exchange (or parallel tempering) MC algorithm exchange  is used, where several systems with identical disorder configurations but different temperatures are simulated in parallel and sometimes exchanged. The temperatures were chosen in a range around $`T_c`$, such that the acceptance rate of temperature exchanges were no less than $`15\%`$.

We begin with a brief discussion of the clean model (all $`𝐅_{ij}=0`$). Figure 1 shows a sharp drop in the structure function $`S_𝐊`$ and in the renormalized shear modulus $`c_{66}`$, and a discontinuous jump in the average energy, with a corresponding approximate delta-function peak in the specific heat $`C_V`$. The very sharp variation in these quantities at $`T_c0.076`$ signifies a discontinuous first order transition in accordance with experiments first-order-melting-expr  as well as numerical simulations of frustrated $`XY`$-models Hu-1997  and vortex line models helicity ; Nordborg-Blatter . This is also confirmed by looking at histograms for the energy or the structure function, which show an increasingly growing double-peak structure developing for the largest system sizes.

We now turn to the disordered model. In the maximally disordered model the dynamics becomes very slow at low temperatures, and it is important to check the convergence of the MC simulation. To do so we monitor several quantities averaged over increasingly longer simulation times, after discarding an equally long time for equilibration. In Fig. 2 we plot, as an example, $`\sigma _{\mathrm{rms}}^2`$ calculated from Eq. (4) and 5, and $`\chi _𝐊`$ as a function of MC time for the size $`L=5`$ and the lowest $`T`$ used in the exchange MC for this size. As seen the curves saturate and become time independent which indicates that the simulation has converged at around 100000 sweeps. Similar checks are done for all sizes. Typically, 4000–7000 disorder realizations are used to form disorder averages.

In the inset of Fig. 3 we show MC data for $`L^4\sigma _{\mathrm{rms}}^2`$ vs. $`T`$ for several different system sizes, which according to Eq. (6) should become independent of $`L`$ exactly at $`T=T_c`$. Indeed, the curves cross, which indicates that the model has a glass phase below this point, with a non-zero $`T_c0.018`$. The main part shows a finite size scaling collapse of the same data, which gives an estimate of the critical exponent $`\nu `$. The structure function $`S_𝐊`$ stays constant $`1`$ independent of temperature and system size, implying a total lack of crystalline order throughout the glass phase. Also, no signature of the transition is found in the density of dislocations, which remains finite down to the lowest temperatures. In Fig. 4 we plot a finite size scaling collapse using Eq. (9) of $`\chi _𝐊`$ (averaged over the smallest reciprocal lattice vectors). As a measure of the quality of the data collapses we use the error in a $`\chi ^2`$–fit of the scaled data to a polynomial (indicated by the dotted lines in the figures). The critical parameters are estimated by adjusting them to minimize this error, and the errorbars are estimated by varying them until the data no longer scales. The final results for the critical parameters are summarized in Table 1.

In summary, we have introduced a model of vortex lattices subjected to disorder. The model goes beyond an elastic description by including dislocations. Such a description should be valid as soon as the screening length is larger than the inter-vortex distance, $`\lambda a_v`$. The model is studied via MC simulations and finite size scaling. Without disorder the model shows a first order melting transition. For strong disorder, the model has a low-temperature glass phase where the dislocation loops are pinned by the disorder. The glass melts via a continuous transition characterized by the critical exponents in Table 1. We propose that the model gives an effective description of the vortex glass transition in type-II superconductors. Experiments usually find a correlation length exponent in the range $`1<\nu <2`$ vortexglass-experiment , which compares well with the present value $`\nu =1.3\pm 0.2`$. Very recent simulations on disordered 3D $`XY`$-models find evidence for a glass transition with $`\nu =1.5\pm 0.12`$ or $`\nu =1.1\pm 0.2`$, $`\eta =0.5\pm 0.1`$ Olsson-Kawamura , consistent with the present study. It would be interesting to extend the study to include dynamic quantities, and to map out the phase diagram for intermediate disorder strengths and other values of the bare elastic constants.