text:
# Axion Phantom Energy

## I Introduction

Cosmology has become the current commonplace where the greatest problems of physics are being concentrated or clearest manifest themselves. Along with the nature of dark matter, the origin of supermassive black holes and several unexplained huge amounts of released astrophysical energy are only some illustrative examples. However, perhaps the biggest problem of all of physics be the so-called dark energy problem  by which it is known that nearly the seventy percent of the total energy in the universe is in the form of an unobserved vacuum energy which is responsible for the present accelerating expansion of the universe  and whose nature remains being a mystery today . There has been a rather huge influx of papers in recent years trying to shed some light onto the existence and the kind of the possible stuff that may make up dark energy .

Four main candidates to represent dark energy have been hitherto suggested: A positive cosmological constant , the quintessence fields (which may  or may not  be tracked), some generalizations of the Chaplygin gas , and the so-called tachyon model of Padmanabhan et al. . Even though such models could be all accommodated to present observational constraints, these models are posing new or traditional problems in such a way that none of them become completely satisfactory. In particular, cosmological data from current observations do not exclude  but may be suggesting  values of the parameter $`\omega `$ in the perfect-fluid equation of state $`p=\omega \rho `$ of the most popular quintessence models which are smaller than -1. If this were the case, then dark energy would become what is now named phantom energy , an entity violating the dominant energy condition and thereby allowing the natural occurrence of wormholes and ringholes and even their corresponding time machines in the universe .

Caldwell, Kamionkowski and Weinberg have recently noticed  that in the framework of quintessence models phantom energy may lead to a doomsday for the universe - which would take place at a big rip singularity - once clusters, galaxies, stars, planets and, ultimately, nucleons and leptons in it are all ripped apart. Even though such a big rip does not take place in some phantom-energy models that use Chaplygin-like equations of state , it seems to be the simplest and most natural possibility stemming from phantom energy. On the other hand, the big rip cosmic scenario interestingly adds an extra qualitative feature to the known set of cosmological models as it introduces a curvature singularity, other than that at the big bang, at a finite, nonzero value of the cosmic time. However, in spite of the feature that the blowing up of the scale factor appears to be unavoidable in models with equation of state $`p=\omega \rho `$ and $`\omega <1`$, a complete account of the nature and properties of the field theory associated with such models has not been done yet. This paper aims at investigating the characteristics of the massless scalar field which allow emergence of a big rip singularity in the simplest of such phantom energy models, assuming a perfect-fluid equation of state. It will be shown that the stuff making up phantom energy can be interpreted to be a vacuum sea of cosmic axions which can be described in terms of the kind of rank-three tensor field strengths predicted in supergravity and strings theories.

This paper can be outlined as follows. In Sec. II we argue in favor of the idea that superlight axions are the source of cosmic phantom energy. A simple cosmological model accounting for a big rip singularity in the case that vacuum is filled with phantom energy is discussed in Sec. III. Sec. IV contains the solution of the phantom field theory within the simple cosmological model of Sec. III. We check that both the scalar field and its potential also have a singularity at the big rip. The results are summarized in Sec. V.

## II Axions as the source of phantom energy

The definition of the massless scalar field $`\varphi `$ which is assumed to make up dark energy is usually taken to be the conventional simplest one; that is, in terms of pressure $`p`$ and energy density $`\rho `$,

$$\rho =\frac{1}{2}\dot{\varphi }^2+V(\varphi ),p=\frac{1}{2}\dot{\varphi }^2V(\varphi ),$$
(1)
where $`V(\varphi )`$ is the field potential. From the equation of state $`p=\omega \rho `$ and Eq. (1) it immediately follows that

$$\rho =\frac{\dot{\varphi }^2}{1+\omega }.$$
(2)
Thus, if the weak energy condition $`\rho 0`$ is taken to be always satisfied , the requirement $`p+\rho <0`$, $`\omega <1`$ from phantom energy in this kind of models necessarily implies that the field $`\varphi `$ ought to be pure imaginary; that is to say, if the vacuum energy density for a phantom vacuum as referred to any timelike observer has to be positive, then the massless scalar field making up the phantom stuff should be pure imaginary. I will argue in what follows that, in the classical framework, such a massless, pure imaginary scalar field actually represents an axion describable as a rank-three antisymmetric tensor field, considering after the associated cosmic theory for such a field. In fact, the Lorentzian action that couples such an axion field to gravity can generally be written as

$$S=d^4x\sqrt{g}\left(\frac{R}{16\pi G}A^2+L_m\right),$$
(3)
where $`L_m`$ is the Lagrangian for observable matter and $`A^2A_{\mu \nu \alpha }A^{\mu \nu \alpha }`$ if we choose the axion to be given by a three-form $`A=dB`$ field strength so that $`dA=0`$, that is as a rank-three antisymmetric tensor field strength of a type arising in supergravity-theory motivated quantum-gravity solutions . The equations of motion derived from action (3) are

$$R_{\mu \nu }\frac{1}{2}g_{\mu \nu }R=16\pi G\left(3A_{\mu \nu }^2\frac{1}{2}g_{\mu \nu }A^2+T_{\mu \nu }^{(m)}\right)$$
(4)
$$dA=0,$$
(5)
in which $`A_{\mu \nu }^2=A_{\mu \alpha \beta }A_\nu ^{\alpha \beta }`$, the asterisk denotes Hodge dual, and $`T_{\mu \nu }^{(m)}`$ is the momentum-energy tensor for ordinary observable matter. One now can check that any explicit solution to these equations of motion subject to the usual Friedmann-Robertson-Walker (FRW) symmetry for the metric and the corresponding spherically symmetric ansatz for the axion antisymmetric tensor $`A=f(r)ϵ`$ (with $`ϵ`$ the volume form which, when integrated over a surface of constant radius, yields the area of the unit three-sphere $`A=2\pi ^2/\mathrm{\Gamma }(2)`$), if we set e.g. $`L_m=0`$ and $`f(r)=\mathrm{Const}.`$, is the same as the solution obtained from the equations of motion derived from the FRW action integral describing coupling of a real massless scalar field $`\varphi `$ to gravity where an extra boundary term accounting for the axion properties, $`a^3\dot{\varphi }\varphi |_0^T`$ ($`a`$ being the scale factor), is added; i.e. in the flat geometry case

$$S^{}=\frac{1}{16\pi G}_0^T𝑑ta^3\left[\frac{\dot{a}^2}{N}+8\pi G\frac{\dot{\varphi }^2}{N}\right]\frac{a^3\dot{\varphi }\varphi }{N}|_0^T,$$
(6)
where $`N`$ is the lapse function. The key point then is that the solution to the above equations of motion is even also preserved when we omit such an extra boundary term in the action for the scalar field $`\varphi `$, provided this scalar field is simultaneously rotated to the pure imaginary axis, $`\varphi i\mathrm{\Phi }`$. If the scalar field is equipped with a potential $`V(\varphi )`$, then an extra term $`_0^T𝑑tNA^3V(\varphi )`$ should be added to the action. In that case, the axionic action can be obtained by simply rotating $`\varphi i\mathrm{\Phi }`$ in an action containing the above potential term, without any extra boundary term.

Such an property, which was first noticed in Euclideanized solutions such as wormholes and other instantons describing nucleation of baby universes , is actually independent of the metric signature and, therefore, applies also to our Lorentzian cosmological context. Thus, one can interpret that the stuff that makes up phantom energy can be regarded to be a rank-three antisymmetric tensor axionic field. It is usually thought that axions can explain the absence of an electrical dipole moment for the neutron and thereby solve the so-called strong CP problem . The axions are chargeless and spinless particles with very tiny mass which interact with ordinary matter only very weakly. Such particles are believed to have been abundantly produced in the big bang. It is worth noticing that whereas relic axions are an excellent candidate for the dark matter in the universe , their vacuum quantum background could make cosmic phantom energy.

The coincidence and fine tuning problems could be thought to become exacerbated in the present scenario where one sets $`\omega `$ constant and $`<1`$. However, in dark energy models such as the generalized Chaplygin gas  and tachyon models , dark energy and dark matter are described as separate limiting cases from an existing unique field. Partly inspired by these models one could naively assume the existence of a unique axion field which, when excited, would make dark matter, and when at its vacuum ground state would be the source of phantom energy. Coincidence time could then be interpreted as the time when both the vacuum and the excited states are approximately equally populated. Of course, this would not solve the coincidence and fine tuning problems but provided some explanation to these problems and to the prediction of generating such strangely small amounts of homogeneously distributed axions. It remains nevertheless an intriguing possibility that such a small value of the cosmological constant be supplied by the potential energy density of an ”ultrainvisible” axion field which can be dubbed quintessence axion .

## III The cosmological model

We shall consider in what follows a simple model where the massless scalar field $`\varphi `$ is a quintessential field, equipped with a potential $`V(\varphi )`$, which is minimally coupled to Hilbert-Einstein gravity. The action for this model is

$$S_c=d^4x\sqrt{g}\left(\frac{R}{16\pi G}+\frac{1}{2}^i\varphi _i\varphi V(\varphi )+L_m\right),$$
(7)
where $`L_m`$ again is the Lagrangian for observable matter fields. If we again restrict ourselves to the case where (i) $`L_m=0`$, (ii) both the scalar field and Hilbert-Einstein gravity satisfy the Friedmann-Robertson-Walker symmetry, and (iii) we assume a perfect fluid equation of state, then by integrating the conservation law for cosmic energy, $`d\rho =3(p+\rho )da/a`$, it turns out that the energy density for the quintessence scalar field will be given by

$$\rho =Ra^{3(1+\omega )},$$
(8)
in which $`R`$ is an integration constant. From the Friedmann equation derived from action (7) subject to the FRW symmetry and $`L_m=0`$,

$$\left(\frac{\dot{a}}{a}\right)^2=Aa^{3(1+\omega )},A=\frac{8\pi GR}{3},$$
(9)
we can obtain the solution

$$a(t)=\left[a_0^{3(1+\omega )/2}+\frac{3(1+\omega )\sqrt{A}}{2}(tt_0)\right]^{\frac{2}{3(1+\omega )}},$$
(10)
where $`a_0`$ and $`t_0`$ are the initial radius and time, respectively. Note that for $`\omega >1`$ this solution describes an accelerating universe whose scale factor increases towards infinity as $`t\mathrm{}`$. We are here most interested in the observationally not excluded yet case where $`\omega <1`$ which corresponds to the so-called phantom dark energy for which the dominant energy condition is generally violated, i.e.

$$p+\rho <0,$$
(11)
even though the energy density is surprisingly ever increasing. Notice furthermore that in this case the scale factor blows up at a finite time,

$$t_{}=t_0+\frac{2}{3(|\omega |1)a_0^{3(|\omega |1)/2}},$$
(12)
giving rise to a ”big rip” . One would expect that in the vecinity of time $`t_{}`$ the size of all local geometrical objects increased in a rather dramatic way. We see that the larger $`a_0`$ and $`|\omega |`$ the nearer the doomsday, and hence the time expected to have macroscopically increased geometrical objects in our universe.

## IV The field theory

Now we shall consider the field theory which is associated with phantom energy, starting with the general formalism for any value of $`\omega `$. From the equation of motion for the scalar field $`\varphi `$,

$$\ddot{\varphi }+3H\dot{\varphi }=\frac{dV(\varphi )}{d\varphi },$$
(13)
where $`H\dot{a}/a`$ and $`V(\varphi )`$ again is the potential for the field $`\varphi `$, and the expression for kinetic energy of the field \[Eq. (2)\], we can obtain a relation between $`dV/d\varphi `$ and the scale factor which can be expressed in the form

$$\frac{3}{2}\sqrt{A(1+\omega )}(1\omega )a^{3(1+\omega )}=\frac{dV}{d\varphi }.$$
(14)
On the other hand, integrating the equation of motion for the field $`\varphi `$ we can derive the dependence of this field with time $`t`$.

$$\varphi (t)=\frac{2}{3\sqrt{1+\omega }\mathrm{}_p}\mathrm{ln}\left[\frac{a_0^{3(1+\omega )/2}+\frac{3(1+\omega )\sqrt{A}}{2}(tt_0)}{B_0}\right],$$
(15)
where $`\mathrm{}_p=\sqrt{R/A}`$ is the Planck length and $`B_0`$ is an integration constant.

From Eqs. (10) and (15) an expression for the scale factor in terms of the scalar field can now be derived. It is:

$$a(\varphi )=B_0^{\frac{2}{3(1+\omega )}}\mathrm{exp}\left(\frac{\mathrm{}_p\varphi }{\sqrt{1+\omega }}\right).$$
(16)
Using Eqs. (14) and (16) we finally obtain the expression of the field potential for any $`\omega <1/3`$ as a function of the field $`\varphi `$

$$V(\varphi )=V_0+\frac{1}{2}(1\omega )\sqrt{R}B_0^2\mathrm{exp}\left(3\sqrt{1+\omega }\mathrm{}_p\varphi \right),$$
(17)
in which $`V_0`$ is a constant. Potentials with this form naturally arise in supergravity models  and have been used in a variety of context, ranging from accelerating expansion models  to cosmological scaling solutions .

Together with the expression for the field $`\varphi `$ in terms of time $`t`$, this expression solves the field-theory problem for any value of $`\omega <1/3`$. According to our discussion above, for axionic phantom energy we have to perform the continuation $`\varphi i\mathrm{\Phi }`$ in addition to taking $`\omega <1`$. In such a case, the solution of the field theory problem would read

$$\mathrm{\Phi }(t)=$$

$$\frac{2}{3\sqrt{|\omega |1}\mathrm{}_p}\mathrm{ln}\left[\frac{a_0^{3(|\omega |1)/2}\frac{3(|\omega |1)\sqrt{A}}{2}(tt_0)}{B_0}\right]$$
(18)
$$V(\mathrm{\Phi })=V_0+\frac{1}{2}(1+|\omega |)\sqrt{R}B_0^2\mathrm{exp}\left(3\sqrt{|\omega |1}\mathrm{}_p\mathrm{\Phi }\right).$$
(19)
As to the physical motivation for the exponential axion potential (19), I will briefly comment on the influence that the matter-field sources may have on the stability of the present axion model of phantom energy. If, corresponding to the physically interesting case that $`V_m0`$, we would let the influence of matter sources on the expansion factor immediately before of phantom domination to have the form of a perturbation, and analyze the resulting model via phase space , it can be seen that if we choose $`V_0=0`$ then potential (19) and the solution given by Eq. (10) (with $`\omega <1`$) and Eq. (18), that is for

$$t=t_0+\frac{2\left[1\left(\frac{a_0}{a}\right)^{3(|\omega |1)/2}\right]}{3a_0^{3(|\omega |1)/2}\left(|\omega |1\right)\sqrt{A}},$$

correspond (see e.g. Ref. ) to a dynamical attractor for phantom energy domination and catastrophic big rip with $`\omega <1`$ and

$$\mathrm{\Gamma }=\frac{V(\mathrm{\Phi })^{}V(\mathrm{\Phi })^{\prime \prime }}{\left[V(\mathrm{\Phi })^{}\right]^2}=1,$$

at a critical point defined by

$$8\pi G\lambda _c=\frac{V(\mathrm{\Phi })^{}}{V(\mathrm{\Phi })}=3\sqrt{|\omega |1}\mathrm{}_p,$$

for which $`\omega =1\lambda _c^2/3`$ (See also phantom models with Born-Infeld type Lagrangians .)

Note that the potential $`V(\mathrm{\Phi })`$ takes on a value

$$V=V_0+\frac{1}{2}(1+|\omega |)\sqrt{R}a_0^{3(|\omega |1)}$$
(20)
at the initial time $`t=t_0`$ to steadily increase thereafter up to infinity when $`t=t_{}`$, at the big rip. From that moment on the potential would continuously decrease, tending to reach its constant minimum value $`V_0`$ as $`t\mathrm{}`$ (see solid curves in Fig. 1). That behaviour is in sharp contrast with that of potential $`V(\varphi )`$ for dark energy with $`\omega >1`$, which starts with a value

$$V=V_0+\frac{1}{2}(1\omega )\sqrt{R}a_0^{3(1+\omega )}$$

to monotonously decrease down to the value $`V_0`$ as $`t\mathrm{}`$, such as it is also shown by the dashed curve in Fig. 1. We have solved in this way the field theory which is associated with dark and phantom energy if the latter is interpreted as originating from the existence of an axionic field which is expressible as a massless, purely imaginary scalar field. The explicit emergence of axions in the cosmological context of phantom dark energy appears to be just another example where such rather elusive particles are invoked as candidates for inexorably needed cosmic constituents which are themselves defined to be unobservable. Actually, apart of being a key ingredient to solve the strong-CP problem  or to represent quantum-gravity topology changes , axions have already been widely used to solve a variety of astrophysical and cosmological shortcomings .

It appears of some interest to briefly comment next on a potentially fruitful consequence from the connection of axions with our cosmological model. If we are actually living in a universe with phantom energy where there will occur a big rip singularity which can never be circumvented by causally-violating connections to the future, then the contracting branch of solution (9) can be physically disregarded as the big-rip singularity is a true curvature singularity. In such a case, only the branch of potential $`V(\mathrm{\Phi })`$ before the singularity can have physical significance, and hence the global minimum of that potential would become placed at $`t=t_0`$ with the value given by Eq. (20). It corresponds to a value of the field

$$\mathrm{\Phi }_\mathrm{m}=\frac{2}{3\mathrm{}_p\sqrt{|\omega |1}}\mathrm{ln}\left(B_0a_0^{3(|\omega |1)/2}\right).$$
(21)
Interpreting potential (19) as that axion potential resulting from the QCD nonperturbative effects, we can take  $`\mathrm{\Phi }_\mathrm{m}=f_A\theta _{\mathrm{eff}}`$, with $`f_A`$ a constant and $`\theta _{\mathrm{eff}}`$ the effective $`\theta `$ parameter resulting after the diagonalization of the quark masses in the QCD Lagrangian. Setting then $`f_A=2/[3\mathrm{}_p\sqrt{|\omega |1}]`$, $`\theta _{\mathrm{eff}}=\mathrm{ln}\left(B_0a_0^{3(|\omega |1)/2}\right)`$ and $`B_0=1`$ it follows that, since $`\theta _{\mathrm{eff}}<10^9`$ , the initial radius of the phantom energy universe $`a_0`$ should be very close to unity, that is the onset of such a phantom energy regime, if it ever at all occurred, must always be placed at nearly just the present epoch. Actually, if $`\theta _{\mathrm{eff}}`$ is cancelled to completely solve the strong CP problem , then the choice $`B_0=1`$ implied that $`a_0`$ would exactly satisfy $`a_0=1`$. In such a case, the observational setting of the value of the state equation parameter $`\omega `$ would also set the value of $`a_0`$, and hence of time $`t_{}`$.

A potential problem with the proposal in this paper could be that whereas current axion theories have a cut off on the scale of inflationary energy, typically at around the GUT characteristic energy of $`10^{16}`$ GeV , in order to avoid phantom energy decaying into gravitons, a cut off of $`<100`$ MeV should be introduced in the effective phantom theory . Not with standing, while the former cut off appears relevant for dark matter axions which can initially be highly excited, the axion field proposed in this paper as the source of phantom energy can never be excited outside the vacuum. On the other hand, if we keep the current values for color anomaly of the Peccei-Quinn symmetry, the pion and quark masses and the pion decay constant, then the above Peccei-Quinn symmetry-breaking scale, $`f_A10^{19}(|\omega |1)^1`$ GeV, would imply a mass of the phantom axions

$$m_A0.62\mathrm{eV}\frac{10^7\mathrm{GeV}}{f_A}10^{12}\sqrt{|\omega |1}\mathrm{eV},$$

i.e. just at the extremal minimum of the allowed values for the axion mass, and clearly well below the cut off required by Carroll et al.  for phantom energy to be stable.

## V Summary

Imposing the weak energy condition, the negative kinetic energy for a phantom field is interpreted as being originated from super-light axions. We have built up a simple cosmological model encompassing constant equations of state for both $`\omega >1`$ and $`\omega <1`$. In the latter case the universe evolves toward a singularity at finite time. The field theory is then solved for all cases, deriving an increasing exponential potential for the phantom field which also diverges at the singularity and shows a dynamical attractor also when matter fields are present that inexorably leads to a catastrophic big rip.

###### Acknowledgements.

The author thanks Carmen L. Sigüenza for useful discussions. This work was supported by DGICYT under Research Project BMF2002-03758.


format: text

text:
# The Thermal Evolution following a Superburst on an Accreting Neutron Star

## 1. Introduction

Type I X-ray bursts from accreting neutron stars in low mass X-ray binaries involve unstable thermonuclear burning of accreted hydrogen (H) and helium (He) (Lewin, van Paradijs, & Taam 1995). In the last few years, long term monitoring of X-ray bursters by BeppoSAX and the Rossi X-Ray Timing Explorer (RXTE) has revealed a new class of very energetic Type I X-ray bursts, now known as ”superbursts” (see Strohmayer & Bildsten 2003; Kuulkers 2003 for reviews). The $`10^{42}\mathrm{erg}`$ energies and several hour durations of superbursts are $`100`$$`1000`$ times greater than usual Type I bursts. In addition, they are rare: so far 8 have been seen from 7 sources, with recurrence times not well-constrained, but estimated as $`1`$ year (Kuulkers 2002; in ’t Zand et al. 2003; Wijnands 2001), instead of hours to days for usual Type I bursts.

The current picture is that superbursts are due to unstable ignition of carbon at densities $`\rho 10^8`$$`10^9\mathrm{g}\mathrm{cm}^3`$. Hydrogen and helium burn at $`\rho 10^5`$$`10^6\mathrm{g}\mathrm{cm}^3`$ via the rp-process (Wallace & Woosley 1981), producing chiefly heavy elements beyond the iron group (including nuclei as massive as $`A=104`$; Schatz et al. 2001), but with some residual carbon (mass fraction $`X_C0.01`$$`0.1`$) (Schatz et al. 2003). Cumming & Bildsten (2001) (hereafter CB01) showed that this small amount of carbon can ignite unstably once the mass of the ash layer reaches $`10^{25}\mathrm{g}`$ (see also Strohmayer & Brown 2002). This fits well with observed superburst energies for $`X_C0.1`$ and an energy release from the nuclear burning of 1 MeV per nucleon. The low thermal conductivity of the rp-process ashes gives a large temperature gradient and ignition at the required mass (CB01). The heavy nuclei may also photodisintegrate to iron group during the flash, enhancing the nuclear energy release (Schatz, Bildsten, & Cumming 2003). Therefore superbursts offer an opportunity to study the rp-process ashes.

Previous authors used one-zone models to estimate the time-dependence of the flash (CB01; Strohmayer & Brown 2002). In this paper, we present the first multi-zone models of the cooling phase of superbursts. Unlike normal Type I bursts, the time to burn the fuel is much less than the convective turnover time. We therefore assume that the fuel burns locally and instantaneously in place, without significant vertical mixing. We do not calculate ignition conditions, but rather treat the amount of energy deposited and the thickness of the fuel layer as free parameters111This is a similar approach to Eichler & Cheng (1989) who studied the thermal response of a neutron star to energy deposition at different depths. However, the transient events they consider are less energetic than superbursts.. In §2, we describe our calculations of the subsequent thermal evolution of the layer, and present a simple analytic model which helps to understand the numerical results. At late times, the flux evolves as a power law in time rather than the exponential decay found by CB01 for a one-zone model. In §3, we use the long term flux evolution of the layer to predict the timescale of quenching of Type I bursts after the superburst, and compare to observations.

## 2. Time Evolution of the Superburst

After the fuel burns, the cooling of the layer is described by the entropy equation

$$c_P\frac{T}{t}=ϵ_\nu \frac{1}{\rho }\frac{F}{r}$$
(1)
where the heat flux is $`F=K(T/r)`$, and $`ϵ_\nu `$ is the neutrino energy loss rate. The layer remains in hydrostatic balance, in which case a useful independent coordinate is the column depth $`y`$ into the star (units: $`\mathrm{g}\mathrm{cm}^2`$), where $`dy=\rho dr`$, giving a pressure $`P=gy`$. The surface gravity is $`g=(GM/R^2)(1+z)`$, where $`1+z=(12GM/Rc^2)^{1/2}`$ is the gravitational redshift factor. In this paper, we assume $`M=1.4M_{}`$ and $`R=10\mathrm{k}\mathrm{m}`$, giving $`z=0.31`$ and $`g_{14}=g/10^{14}\mathrm{cm}\mathrm{s}^2=2.45`$.

To find the temperature profile just after the fuel burns, we deposit an energy $`E_{\mathrm{nuc}}=E_{17}10^{17}\mathrm{erg}\mathrm{g}^1`$ throughout the layer. Since carbon burning to iron gives $`10^{18}\mathrm{erg}\mathrm{g}^1`$, we expect $`E_{17}=1`$ to correspond to $`X_C0.05`$$`0.1`$, depending on how much energy is contributed by photodisintegration (Schatz et al. 2003). At each depth, we calculate the temperature of the layer, $`T_f`$, from $`_{T_i}^{T_f}c_P𝑑T=E_{\mathrm{nuc}}`$, where $`T_i`$ the initial temperature. Following CB01, a simple analytic estimate of $`T_f`$ is as follows. The electrons are degenerate and relativistic for $`\rho 10^7\mathrm{g}\mathrm{cm}^3`$, giving $`\rho Y_e=1.5\times 10^8\mathrm{g}\mathrm{cm}^3P_{26}^{3/4}`$, Fermi energy $`E_F=2.7\mathrm{MeV}P_{26}^{1/4}`$, and pressure scale height $`H=y/\rho =67\mathrm{m}P_{26}^{1/4}Y_e/g_{14}`$, where $`P=P_{26}10^{26}\mathrm{erg}\mathrm{cm}^3`$. The heat capacity $`c_P`$ is determined mainly by the electrons, $`c_P\pi ^2(Y_ek_B/m_p)(k_BT/E_F)=2.6\times 10^7\mathrm{erg}\mathrm{g}^1\mathrm{K}^1T_9Y_eP_{26}^{1/4}`$. Integrating and assuming $`T_iT_f`$, we find the convectively stable temperature profile $`T_f=3.6\times 10^9\mathrm{K}E_{17}^{1/2}P_{26}^{1/8}Y_e^{1/8}`$, insensitive to depth and depending mainly on $`E_{17}`$.

Our thermal evolution code uses the method of lines, in which the right hand side of equation (1) is differenced over a spatial grid, and the resulting set of ordinary differential equations integrated using a stiff integrator. We choose a uniform grid in $`\mathrm{sinh}^1(\mathrm{log}y/y_b)`$, which concentrates grid points around the base of the layer $`y_b`$, resolving the initial temperature discontinuity. We place the outer boundary at $`y=10^8\mathrm{g}\mathrm{cm}^2`$, and set flux $`T^4`$ there; at the inner boundary, typically $`y10^{14}\mathrm{g}\mathrm{cm}^2`$, we assume vanishing flux. We assume the layer is heated before the flash by a $`10^{21}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ flux from the crust, and that all the fuel burns to 56Fe. Our results are not sensitive to the details of the grid, or boundary conditions (for times longer than the thermal time at the top zone). We calculate the equation of state, opacity, neutrino emissivity, and heat capacity as described by Schatz et al. (2003).

Figure 1 shows temperature and flux profiles 10 minutes, 1 hour, 1 day, and 10 days after ignition for a model with $`E_{17}=1`$ and $`y_b=10^{12}\mathrm{g}\mathrm{cm}^2`$. Figure 2 shows a series of lightcurves for different $`y_b`$ and $`E_{17}`$. At early times, as the outer parts of the layer thermally adjust, the radiative flux depends mostly on $`E_{17}`$. At late times, after the cooling wave reaches the base of the layer, the flux depends mostly on $`y_b`$, and falls off as a power law $`Ft^{4/3}`$. Figure 3 shows the cumulative energy release for $`y_b=10^{12}`$ and $`10^{13}\mathrm{g}\mathrm{cm}^2`$ and $`E_{17}=1`$$`3`$. In the first few hours, the energy released from the surface is $`10^{42}\mathrm{ergs}`$, the exact value being mainly sensitive to $`E_{17}`$, rather than depth. A significant fraction of the heat is initially conducted inwards and released on a longer timescale, as pointed out by Strohmayer & Brown (2002).

The physical reason for the late-time power law flux decay is that as time evolves, the peak of the temperature profile moves to greater depths where the thermal timescale to the surface is longer222A similar problem is ohmic decay of crustal magnetic fields, where power law decay is also expected (Sang & Chanmugam 1987; Urpin, Chanmugam, & Sang 1994) (see Fig. 1). The simplest analytic model is a slab with constant thermal diffusivity $`D`$, whose temperature is perturbed close to the surface, for example at a depth $`x=a`$ (where $`x=0`$ is the surface). For a delta-function perturbation initially, the temperature evolution is given by the Green’s function333A simple way to obtain this result is to apply the method of images to the Green’s function for an unbounded domain $`T(x,t)t^{1/2}\mathrm{exp}\left(x^2/4Dt\right)`$. Eichler & Cheng (1989) derive a similar result for a power law dependence of conductivity on depth, which also shows self-similar behavior at late times (see Lyubarsky, Eichler, & Thompson 2002 for a recent application to cooling of SGR 1900+14 after an outburst).

$$T(x,t)=\frac{1}{\sqrt{\pi Dt}}\mathrm{sinh}\left(\frac{ax}{2Dt}\right)\mathrm{exp}\left(\frac{x^2+a^2}{4Dt}\right)$$
(2)
and the surface flux is $`F(T/x)_{x=0}t^{3/2}\mathrm{exp}(\tau /t)`$, where $`\tau `$ is the thermal time at the initial heating depth $`\tau =4a^2/D`$. For an initial ”top hat” temperature profile $`T(x<a)=1`$, $`T(x>a)=0`$, the surface flux is $`F(\tau /t)^{1/2}\left[1\mathrm{exp}(\tau /t)\right]`$. For $`t<\tau `$, before the cooling wave reaches the base of the layer, $`Ft^{1/2}`$; for $`t>\tau `$, the solution is independent of the initial temperature profile, and $`Ft^{3/2}`$.

The numerical results show a similar behavior, although with different power law indices. The relevant timescale in this case is $`t_{\mathrm{cool}}=H^2/D`$, where $`D=K/\rho c_P`$. The electron conductivity is $`K=\pi ^2n_ek_B^2T/3m_{}\nu _c`$, where $`m_{}=E_F/c^2`$, and $`\nu _c`$ is the electron collision frequency. When electrons dominate the heat capacity, the thermal diffusivity takes the particularly simple temperature-independent form $`D=c^2/3\nu _c`$. For electron-ion collisions, $`\nu _c=9.3\times 10^{16}\mathrm{s}^1P_{26}^{1/4}Z^2/A\mathrm{\Lambda }_{ei}/Y_e`$ (e.g. see Appendix of Schatz et al. 1999), giving

$$t_{\mathrm{cool}}=3.8\mathrm{hrs}y_{12}^{3/4}\left(\frac{Y_eZ^2/A\mathrm{\Lambda }_{ei}}{6}\right)\left(\frac{g_{14}}{2.45}\right)^{5/4}$$
(3)
(see also eq.  of CB01), where we insert the appropriate numbers for 56Fe composition. The simple “top hat” solution for constant conductivity motivates a fit to the numerical solutions,

$$F_{25}=0.2t_{\mathrm{hr}}^{0.2}E_{17}^{7/4}\left[1\mathrm{exp}\left(0.63t_{\mathrm{cool}}^{4/3}E_{17}^{5/4}t_{\mathrm{hr}}^{1.13}\right)\right],$$
(4)
where $`t_{\mathrm{hr}}=t/1\mathrm{hour}`$. For $`t>t_{\mathrm{cool}}`$, $`F_{25}=0.13E_{17}^{1/2}(t/t_{\mathrm{cool}})^{4/3}`$. The transition from $`Ft^{0.2}`$ to $`Ft^{4/3}`$ occurs when $`t/t_{\mathrm{cool}}E_{17}^{1.1}`$.

Equation (4) fits the numerical results to better than a factor of two for models without substantial neutrino emission. As emphasised by Strohmayer & Brown (2002), neutrino cooling is important for large carbon fractions: it depresses the flux at $`t5`$$`10`$ hours for the models with $`E_{17}=2`$ and $`3`$, $`y=10^{13}\mathrm{g}\mathrm{cm}^2`$ in Figure 2. Whenever neutrinos dominate the cooling, the peak temperature is large enough that emission is by pair annhilation. A good fit to the neutrino energy loss rate is $`ϵ_\nu 10^4\mathrm{erg}\mathrm{g}^1\mathrm{s}^1T_9^{12}y_{12}^{3/2}`$, giving a cooling time $`t_\nu =c_PT/ϵ_\nu =2.5\times 10^{12}\mathrm{s}y_{12}^{5/4}T_9^{10}`$. Inserting the peak temperature from equation (2) gives $`t_\nu 300\mathrm{hrs}E_{17}^5`$. Neutrinos dominate when $`t_\nu <t_{\mathrm{cool}}`$, or when $`E_{17}>2.3y_{12}^{3/20}`$.

## 3. Comparison to Observations

The cooling curves in Figure 2 compare well with observed lightcurves, including a rapid initial decay on hour timescales, followed by an extended tail of emission (as observed following some superbursts, e.g. KS 1731-260, Kuulkers et al. 2002; Ser X-1, Cornelisse et al. 2002). We will present a detailed comparison with the observed superburst lightcurves in a future paper. The initial decay from the peak depends mostly on $`E_{17}`$, and so it should be possible to constrain the amount of fuel consumed in the superburst. Our models do not resolve the peak itself, since this depends on the details of how the burning propagates out to the surface; however, for $`E_{17}2`$$`3`$, the flux exceeds the Eddington flux, $`F_{\mathrm{Edd}}=3\times 10^{25}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1/(1+X)`$, where $`X`$ is the H fraction, for timescales of minutes. Superburst peak luminosities are generally less than the Eddington luminosity (Kuulkers 2003), implying $`E_{17}2`$. The one exception is the superburst from 4U 1820-30, which showed dramatic photospheric radius expansion lasting for several minutes (Strohmayer & Brown 2002). This is consistent with the proposal that this source, which accretes and burns He rich material, produces large quantities of carbon (Strohmayer & Brown 2002; Cumming 2003a). The transition to the late-time power law occurs after $`t4\mathrm{h}E_{17}^{1.11}y_{12}^{3/4}`$ (eq. ), which corresponds to $`F_{25}0.13E_{17}^2`$. It may therefore be possible to measure the power law decay using superburst tails, although this depends upon being able to subtract out the underlying accretion luminosity, $`F_{\mathrm{accr},25}0.1(\dot{M}/0.1\dot{M}_{\mathrm{Edd}})`$, in a reliable way.

Another way to probe the late-time cooling is to use the remarkable observation that Type I bursts disappear (are “quenched”) for $`t_{\mathrm{quench}}`$ weeks following the superburst (e.g. Kuulkers 2003). CB01 proposed that the cooling flux from the superburst temporarily stabilizes the H/He burning. An estimate of the critical stabilizing flux, $`F_{\mathrm{crit}}`$, is as follows. The condition for temperature fluctuations to grow, and unstable He ignition to occur, is $`\nu ϵ_{3\alpha }=\eta ϵ_{\mathrm{cool}}`$ (Fushiki & Lamb 1987), where $`ϵ_{3\alpha }`$ is the triple alpha ($`3\alpha `$) energy production rate, $`ϵ_{\mathrm{cool}}`$ is a local approximation to the cooling rate, and $`\nu `$ and $`\eta `$ are the respective temperature sensitivities. For a large flux from below, the He burns stably before reaching this ignition condition, at a depth where the time to accumulate the layer equals the He burning time, $`y/\dot{m}=YQ_{3\alpha }/ϵ_{3\alpha }`$, where $`Y`$ is the He mass fraction, $`\dot{m}`$ is the local accretion rate per unit area, and $`Q_{3\alpha }=5.84\times 10^{17}\mathrm{erg}\mathrm{g}^1=0.606`$ MeV per nucleon is the $`3\alpha `$ energy release. At the transition from unstable to stable burning, both criteria are satisfied at the base of the H/He layer. Using the first condition to eliminate $`ϵ_{3\alpha }`$ from the second, and writing $`ϵ_{\mathrm{cool}}F/y`$, gives $`F=\nu \dot{m}Q_{3\alpha }Y/\eta =6.2\times 10^{22}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1(\dot{m}/\dot{m}_{\mathrm{Edd}})(Y/0.3)(\nu /4\eta )`$. Some of this flux is provided by hot CNO burning of accreted H, $`F_Hϵ_Hy=5.8\times 10^{21}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1y_8(Z/0.01)`$ (Cumming & Bildsten 2000; $`Z`$ is the metallicity); the remainder is $`F_{\mathrm{crit}}=FF_H`$. This estimate agrees well with a more detailed calculation using the ignition models of Cumming & Bildsten (2000), in which we find $`F_{\mathrm{crit}}\dot{m}Q_{3\alpha }0.7\mathrm{MeV}`$ per accreted nucleon, almost independent of $`\dot{M}`$. Therefore,

$$F_{\mathrm{crit},22}6(\dot{m}/\dot{m}_{\mathrm{Edd}})$$
(5)
(see also Paczynski 1983a; Bildsten 1995). Equation (5) in the limit $`tt_{\mathrm{cool}}`$ gives

$$t_{\mathrm{quench}}=38t_{\mathrm{cool}}F_{\mathrm{crit},22}^{3/4}E_{17}^{3/8}=6\mathrm{days}y_{12}^{3/4}F_{\mathrm{crit},22}^{3/4}E_{17}^{3/8},$$
(6)
which gives $`t_{\mathrm{quench}}`$ in terms of the thickness of the layer and energy release.

Figure 5 compares the predicted and observed quenching times. The observations of $`t_{\mathrm{quench}}`$ and accretion rates (used to find $`F_{\mathrm{crit}}`$ from eq. ) are taken from Kuulkers (2003) (except for 4U 1636-53, which has a revised upper limit of $`23`$ days, Kuulkers private communication). The observations are upper or lower limits only: nonetheless, the general agreement is very good and supports the quenching picture suggested by CB01. There is much to learn from a careful comparison of superburst lightcurves and the corresponding quenching times, separately constraining both $`E_{17}`$ and $`y_b`$.

## 4. Summary and Conclusions

We have presented the first multi-zone models of the cooling phase of superbursts. The flux decay is not exponential, but power-law (eq. ). For $`t<t_{\mathrm{cool}}`$, where $`t_{\mathrm{cool}}`$ is the cooling time at the base of the layer, the flux depends mostly on the energy release $`E_{17}`$, and is insensitive to depth: the inwards travelling cooling wave does not yet know that the layer has a finite thickness. For $`t>t_{\mathrm{cool}}`$, the flux decays as a power law $`Ft^{4/3}`$, independent of the initial temperature profile. The power law decay at late times gives predicted Type I burst quenching times of weeks (eq. ), consistent with observational limits. Future comparisons of both superburst lightcurves and quenching times with observations will constrain both the thickness of the fuel layer and the energy deposited, particularly when combined with models of normal Type I bursts from the same source (Cumming 2003a,b).

There is still much to be done in terms of theory. Perhaps the most important issues are the physics of the rise (which sets the initial condition for our simulations), and production of the fuel. Important clues to the first are the observed precursors to superbursts, which may be normal Type I bursts ignited by the superburst. Recent progress has been made on the second, with indications from both theory (Schatz et al. 2003; Woosley et al. 2003) and observations (in ’t Zand et al. 2003) that stable burning may be required to produce enough carbon to power superbursts. A self-consistent model of H/He burning, followed by accumulation and ignition of the ashes may require a better understanding of the transition from unstable to stable burning observed in normal Type I bursting (e.g. Cornelisse et al. 2003).

We thank P. Arras, E. Brown, R. Cornelisse, E. Kuulkers, G. Ushomirsky, S. Woosley, and J. in’t Zand for useful comments and discussions. AC is supported by NASA through Hubble Fellowship grant HF-01138 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. JM acknowledges support from DOE grant No. DE-FC02-01ER41176 to the Supernova Science Center/UCSC.


text:
# Bloch Structures in a Rotating Bose-Einstein Condensate

## Abstract

A rotating Bose-Einstein condensate is shown to exhibit a Bloch band structure without the need of periodic potential. Vortices enter the condensate by a mechanism similar to the Bragg reflection, if the frequency of a rotating drive or the strength of interaction is adiabatically changed. A localized state analogous to a gap soliton in a periodic system is predicted near the edge of the Brillouin zone.

The Bloch band is crucial to our understanding of the behaviors of periodic systems, such as electronic states in solids Ziman , atom dynamics in optical lattices Jessen , and Cooper-pair tunneling in small Josephson junctions Likharev . The underlying physics common to all of these is the Bragg reflection that occurs at the edges of the Brillouin zone. The resultant Bloch band leads to various interesting phenomena, such as Bloch oscillations Ziman  (a particle in a periodic potential driven by a weak constant force cannot be accelerated indefinitely but oscillates in real space) and formation of gap solitons Meystre  (localized wave packets arising from the balance between negative-mass dispersion and repulsive interaction); both phenomena have recently been observed with a Bose-Einstein condensate (BEC) in an optical lattice Morsch ; Eiermann .

In this Letter, we show that yet another system — a BEC confined in a rotating harmonic-plus-quartic potential — exhibits a Bloch band structure, and investigate the associated novel phenomena. Seemingly this system has no periodic structure, but may be considered to be a quasi-1D periodic system in the following sense. When the rotating frequency of the potential is high, the quartic potential $`r^4`$, together with the centrifugal one $`r^2`$, produces a Mexican-hat shaped potential Kasamatsu  whose minima form a quasi-1D toroidal geometry. Under such circumstances any perturbation $`V(\theta )`$ that is needed to drive the system into rotation by breaking the axisymmetry of the potential serves as a “periodic” potential since $`V(\theta )=V(\theta +2\pi )`$, where $`\theta `$ denotes the azimuthal angle.

In the present system, we will show that the Bragg reflection causes vortex nucleation. Since this process occurs adiabatically, in contrast to a method invoking dynamical instabilities Sinha ; Madison01 , we can obtain the vortex states without heating the atomic cloud. A stirring technique to produce vortices without heating has also been proposed in Refs. Damski ; Caradoc . Another scheme that adiabatically nucleates vortices is based on phase engineering techniques Dum ; Williams ; Nakahara . We will also show that a localized state is generated even with repulsive interactions, which is analogous to gap solitons in periodic systems Meystre .

We begin by discussing a BEC in a 1D ring to understand the essence of the phenomena. We assume that a rotating potential takes the form $`V(\theta ,t)=\epsilon \mathrm{cos}n(\theta \mathrm{\Omega }t)`$. In a frame rotating at frequency $`\mathrm{\Omega }`$, the potential $`V`$ becomes time-independent, and the Gross-Pitaevskii (GP) equation is given by

$$\left(\frac{^2}{\theta ^2}+i\mathrm{\Omega }\frac{}{\theta }+\epsilon \mathrm{cos}n\theta +\gamma |\psi _{1\mathrm{D}}|^2\right)\psi _{1\mathrm{D}}=\mu \psi _{1\mathrm{D}},$$
(1)
where energy and time are measured in units of $`\mathrm{}^2/(2mR^2)`$ and $`2mR^2/\mathrm{}`$ with $`m`$ and $`R`$ being the atomic mass and the radius of the ring. The wave function is normalized as $`_0^{2\pi }|\psi _{1\mathrm{D}}|^2𝑑\theta =1`$, and $`\gamma `$ characterizes the strength of interaction.

First, let us consider the noninteracting case ($`\gamma =0`$). Introducing a new variable $`\varphi e^{i\mathrm{\Omega }\theta /2}\psi _{1\mathrm{D}}`$, we transform Eq. (1) into the Mathieu equation $`(_\theta ^2+\epsilon \mathrm{cos}n\theta )\varphi =E\varphi `$, where $`E\mu +\mathrm{\Omega }^2/4`$ plays the role of the “total energy” in our Bloch band picture. Then the solution satisfies $`\varphi (\theta +\mathrm{\Theta }_0)=e^{i(p\mathrm{\Omega }/2)\mathrm{\Theta }_0}\varphi (\theta )`$, with $`\mathrm{\Theta }_02\pi /n`$ being the period of the potential $`V`$, and $`p`$ an integer. It follows that $`p\mathrm{\Omega }/2`$ may be regarded as the quasimomentum, indicating that we can move in the quasimomentum space by changing $`\mathrm{\Omega }`$. Since the unit “reciprocal lattice vector” is $`2\pi /\mathrm{\Theta }_0=n`$, the Bragg reflection occurs only between the states whose angular momenta differ by $`n`$, $`2n`$, … (Bragg’s law). Consequently there are $`n`$ independent branches corresponding to $`p=0,\mathrm{},n1`$ modulo $`n`$. Figure 1 shows $`E`$ and the angular momentum $`L=i_0^{2\pi }\psi _{1\mathrm{D}}^{}_\theta \psi _{1\mathrm{D}}d\theta `$ as functions of $`\mathrm{\Omega }`$ for the first and second Bloch bands with $`n=2`$ and $`\epsilon =0.1`$.

Two independent branches exist in this case, and are distinguished by the solid and dotted curves. The Bloch band structure in Fig. 1 indicates that if we adiabatically increase or decrease $`\mathrm{\Omega }`$ across $`\mathrm{\Omega }=2`$, the nonvortex state $`\psi _{1\mathrm{D}}1/\sqrt{2\pi }`$ transforms to the vortex state $`\psi _{1\mathrm{D}}e^{2i\theta }/\sqrt{2\pi }`$, and vice versa. The time scale of the adiabatic change must be much longer than the inverse energy gap $`\epsilon ^1`$.

When $`|\gamma |\epsilon `$, the system exhibits hysteretic behavior, and a loop structure emerges in the Bloch bands Wu  as shown in Fig. 1. We see that the first (second) band forms a loop for $`\gamma =1(1)`$; then the adiabatic nucleation of vortices is not possible with increasing (decreasing) $`\mathrm{\Omega }`$ because the nonlinear Landau-Zener transition Wu  occurs no matter how slow the change of $`\mathrm{\Omega }`$. It can be shown that the states on the loop have dynamical instability that breaks the symmetry of the system $`|\psi _{1\mathrm{D}}(\theta )|=|\psi _{1\mathrm{D}}(\theta +\mathrm{\Theta }_0)|`$. This instability is similar to the one that forms gap solitons in a BEC in an optical lattice Konotop .

We now consider a BEC in a harmonic-plus-quartic potential with a rotating stirrer. We assume that the axial trapping energy $`\mathrm{}\omega _z`$ is much larger than the other characteristic energies (a tight pancake-shaped trap) so that the system is effectively 2D. The external potential in the corotating frame of reference takes the form $`V(𝐫)=r^2/2+Kr^4/4+\epsilon r^2\mathrm{cos}2\theta `$, where $`K`$ is a constant and the third term $`\epsilon (x^2y^2)`$ is a stirring potential which can be produced by laser beams propagating in the $`z`$ direction. Such an anharmonic potential was theoretically considered in Refs. Kasamatsu ; Lundh ; Lundh04  and has recently been realized experimentally Bretin . We normalize the length, time, energy, and wave function by $`d_0(\mathrm{}/m\omega _{})`$, $`\omega _{}^1`$, $`\mathrm{}\omega _{}`$, and $`\sqrt{N}/d_0`$, respectively, with $`\omega _{}`$ and $`N`$ being the frequency of the radial harmonic trap and the number of atoms. The wave function is then normalized as $`|\psi _{2\mathrm{D}}|^2𝑑𝐫=1`$. The time-dependent GP equation in the frame corotating with the stirring potential at frequency $`\mathrm{\Omega }`$ is given by

$$i\frac{\psi _{2\mathrm{D}}}{t}=\left(\frac{1}{2}^2+i\mathrm{\Omega }\frac{}{\theta }+V(𝐫)+g|\psi _{2\mathrm{D}}|^2\right)\psi _{2\mathrm{D}},$$
(2)
where $`g=[\omega _z/(2\pi \omega _{})]^{1/2}4\pi Na/d_0`$ characterizes an effective strength of interaction in 2D Castin  with $`a`$ being the s-wave scattering length.

When $`\mathrm{\Omega }1`$, the system described by Eq. (2) can be approximated by a quasi-1D ring Kasamatsu . Assuming $`\psi _{2\mathrm{D}}f(r)e^{i\mathrm{}\theta }/\sqrt{2\pi }`$, we find that an effective potential for the radial wave function is given by $`\mathrm{}^2/(2r^2)+r^2/2+Kr^4/4`$. This potential has a minimum at $`r(\mathrm{}^2/K)^{1/6}`$ for $`\mathrm{}1`$, and the effective frequency around it is $`\omega _{\mathrm{eff}}\sqrt{6}(\mathrm{}K)^{1/3}`$. If $`\omega _{\mathrm{eff}}`$ is much larger than other characteristic frequencies, the dynamics of $`f(r)`$ can be ignored. After normalization of the time by $`\rho ^2r|f|^2r^2𝑑r`$, the equation of motion for $`\theta `$ reduces to Eq. (1), where $`\gamma `$ is given by $`2g\rho ^2r|f|^4𝑑r`$. Thus, the system described by Eq. (2) exhibits the quasi-1D circular flow for large angular momentum Kasamatsu , and we expect that a Bloch band structure emerges. We show below that this is true even for $`\mathrm{\Omega }1`$ and $`\mathrm{}1`$.

Figure 2 (a) shows the angular momentum $`L=i\psi _{2\mathrm{D}}^{}_\theta \psi _{2\mathrm{D}}d𝐫`$ and the density and phase profiles (insets) of the noninteracting stationary states of Eq. (2note .

We find that the behavior of $`L`$ is similar to that in Fig. 1 (b). When we start from the ground state with $`L=0`$ (1-1), two vortices enter at $`\mathrm{\Omega }1.57`$ (1-2), producing a doubly-quantized vortex (1-3) Lundh . Increasing the rotation frequency $`\mathrm{\Omega }`$ further, we can nucleate two more vortices in the condensate \[(1-4) and (1-5)\]. If we follow the branch starting from $`L=2`$ (2-1), the two vortices escape out of the condensate with increasing $`\mathrm{\Omega }`$ \[(2-2) and (2-3)\], and then four vortices enter at $`\mathrm{\Omega }1.7`$ (2-4). It should be noted that the Bloch band picture holds even for the nonrotating gas, while the ring-shaped profile manifests itself only if the gas is rotating \[cf. (1-1) and (2-3)\].

The frequencies at which the Bragg reflection occurs in Fig. 2 (a) are different from those in Fig. 1. The difference arises from the dispersion relation of the vortex states in 2D, i.e., the relation between the energy and angular momentum. In the 1D case, the energy is given by $`E_m^{1\mathrm{D}}=[\psi _m^{1\mathrm{D}}(_\theta ^2+i\mathrm{\Omega }_\theta )\psi _m^{1\mathrm{D}}+\gamma |\psi _m^{1\mathrm{D}}|^4/2]𝑑\theta =m^2\mathrm{\Omega }m+\gamma /(4\pi )`$ for $`\psi _m^{1\mathrm{D}}=e^{im\theta }/\sqrt{2\pi }`$. The Bragg reflection between $`\psi _m^{1\mathrm{D}}`$ and $`\psi _{m+2}^{1\mathrm{D}}`$ occurs when $`E_m^{1\mathrm{D}}=E_{m+2}^{1\mathrm{D}}`$, i.e., $`\mathrm{\Omega }=2m+2`$. In the 2D case, we employ a variational wave function $`\psi _m^{2\mathrm{D}}=(\pi |m|!)^{1/2}d_m^{|m|1}r^{|m|}\mathrm{exp}[r^2/(2d_m^2)+im\theta ]`$ to obtain $`E_m^{2\mathrm{D}}=(|m|+1)(d_m^2+d_m^2)/2+K(|m|+1)(|m|+2)d_m^4/4\mathrm{\Omega }m+g(2|m|)!/[2^{2|m|+2}\pi |m|!^2d_m^2]`$, where the variational parameter $`d_m`$ is determined from $`E_m^{2\mathrm{D}}/d_m=0`$. When $`K=1`$ and $`g=0`$, the condition for Bragg reflections to occur between $`\psi _m^{2\mathrm{D}}`$ and $`\psi _{m+2}^{2\mathrm{D}}`$ ($`E_m^{2\mathrm{D}}=E_{m+2}^{2\mathrm{D}}`$) for $`m=0`$, $`2`$, and $`4`$ is satisfied at $`\mathrm{\Omega }=1.58,1.83`$, and $`2.03`$, respectively, which are in good agreement with the numerical results in Fig. 2 (a) (corresponding to the points at which the red and blue curves cross with $`L`$ changing by $`\pm 2`$).

The energy gap can also be obtained by a variational method. The matrix element of the stirring potential, for example, between $`\psi _0^{2\mathrm{D}}`$ and $`\psi _2^{2\mathrm{D}}`$ is calculated to be $`|\psi _2^{2\mathrm{D}}\psi _0^{2\mathrm{D}}\epsilon r^2\mathrm{cos}2\theta d𝐫|0.46\epsilon `$ for $`K=1`$. The energy gap at $`\mathrm{\Omega }1.57`$ is then given by $`0.92\epsilon `$, which also agrees very well with the numerical result $`0.91\epsilon `$.

The Bloch band structure in Fig. 2 (a) indicates that if we prepare the nonvortex state at $`\mathrm{\Omega }1.55`$, and adiabatically increase $`\mathrm{\Omega }`$ to above $`\mathrm{\Omega }1.6`$, we obtain a doubly-quantized vortex state by following the red curve. Interestingly, we can nucleate vortices also by decreasing $`\mathrm{\Omega }`$, that is, by preparing a nonvortex state for $`1.6\mathrm{\Omega }1.7`$ and decreasing $`\mathrm{\Omega }`$ to below $`\mathrm{\Omega }1.55`$ by following the second Bloch band (blue curve). Figure 2 (b) illustrates these adiabatic processes. In order to examine the effect of dissipation on the adiabatic processes, we replace $`i`$ with $`i+\lambda `$ on the left-hand side of Eq. (2), where a constant $`\lambda `$ is taken to be $`0.03`$ Choi . We prepare the nonvortex ground state with $`g=0`$, and gradually ramp up $`\mathrm{\Omega }`$ from $`1.52`$ to $`1.6`$ \[red curve in Fig. 2 (b)\], or ramp it down from $`1.62`$ to $`1.54`$ (blue curve). As expected, two density holes with phase singularities come from infinity, and unite to form a doubly-quantized vortex. We note that the behaviors in Fig. 2 (b) are almost the same as for the dissipation-free case ($`\lambda =0`$), indicating that the process is robust against dissipation. The energy gap between the first and second Bloch bands at $`\mathrm{\Omega }1.57`$ is $`0.036`$ for $`\epsilon =0.04`$. According to the Landau-Zener formula, the transition probability between the energy gap, i.e., the probability that the nonvortex state remains so, is below 1% in the situation in Fig. 2 (b), in agreement with our numerical result.

In the case of $`g0`$, the interaction bends the Bloch bands as shown in Fig. 3 (a), yielding hysteresis as in the 1D case.

We note that the region of $`\mathrm{\Omega }`$ in which the Bragg reflection occurs is shifted due to the interaction. The shift is attributed to the difference in interaction energies between the nonvortex and vortex states. To show this we again employ the variational wave function $`\psi _m^{2\mathrm{D}}`$. We find that the condition for the Bragg reflection to occur between $`\psi _0^{2\mathrm{D}}`$ and $`\psi _2^{2\mathrm{D}}`$, i.e., $`E_0^{2\mathrm{D}}=E_2^{2\mathrm{D}}`$, is satisfied at $`\mathrm{\Omega }=1.53`$ for $`g=1.6`$ and $`\mathrm{\Omega }=1.64`$ for $`g=1.6`$, in qualitative agreement with the shifts shown in Fig. 3 (a). Since the interaction term of $`E_m^{1\mathrm{D}}`$ is independent of $`m`$, the shift does not occur in 1D (see Fig. 1).

The shift in the position of the Bragg reflection implies the possibility of the adiabatic nucleation of vortices by changing the strength of interaction using, e.g., the Feshbach resonance Inouye , at fixed $`\mathrm{\Omega }`$. For example, at $`\mathrm{\Omega }1.58`$, the first Bloch band has angular momentum $`L0`$ for $`g1`$ \[left panel in Fig. 3 (a)\], which continuously increases to $`L2`$ with an increase in the interaction to repulsive (red arrow). Similarly, there is a region in which $`L`$ changes from 0 to 2 as $`g`$ decreases from positive to negative \[blue arrow in the right panel in Fig. 3 (a)\]. No dynamical instability is present on the red and blue arrows. Figure 3 (b) illustrates these results, where the initial state is the nonvortex ground state, and the strength of interaction is gradually changed with $`\mathrm{\Omega }`$ held fixed (red and blue curves). We find that the two vortices are nucleated in a manner similar to that in Fig. 2 (b). The repulsive-to-attractive case is particularly interesting because the attractive interaction is usually considered to hinder vortex nucleation Lundh04 .

The adiabatic theorem breaks down when the path enters the region of dynamical instability. A notable example is shown as the green arrow in Fig. 3 (a). For $`g0.6`$ at $`\mathrm{\Omega }1.55`$, a twofold symmetric state \[like the inset (1-2) of Fig. 2 (a)\] becomes dynamically unstable against localization into the state shown in the inset of Fig. 3 (a). This localization is due to the interplay between repulsive interaction and negative-mass dispersion around the edge of the Brillouin zone. The localized state can therefore be regarded as an analog of the gap soliton in a periodic system Meystre . The green curve in Fig. 3 (b) demonstrates the time evolution corresponding to the path shown as the green arrow in Fig. 3 (a). The twofold symmetry as shown by the inset at $`t=300`$ is broken in the course of vortex nucleation, giving rise to two localized states as shown by the insets at $`t=400`$ and $`t=560`$.

Figure 3 (a) shows that a doubly-quantized vortex state with $`L2`$ is dynamically unstable (stable) in most of the region with $`g<0`$ ($`g>0`$). This instability also gives rise to a localized state that undergoes center-of-mass rotation. We find from the numerical analysis that a doubly-quantized vortex state with $`g<0`$ evolves into the localized state, followed by the oscillation between these two states in a manner similar to the split-merge oscillations reported in Ref. Saito .

In conclusion, we have shown that the Bloch band structure arises in a BEC confined in a harmonic-plus-quartic potential with a rotating drive, which enables us to nucleate vortices adiabatically. The physical mechanism of the vortex nucleation is very similar to the Bragg reflection at the edge of the Brillouin zone. Interestingly, we can nucleate vortices not only by increasing the stirring frequency but also by decreasing it, or by changing the strength of interaction at a fixed stirring frequency, which may be called a “Feshbach-induced Bragg reflection”. The adiabatic processes are robust against dissipation due, e.g., to a thermal cloud. Spontaneous localization of the rotating cloud is predicted, which signals the existence of a gap soliton.

This work was supported by the Special Coordination Funds for Promoting Science and Technology, a 21st Century COE program at Tokyo Tech “Nanometer-Scale Quantum Physics”, and a Grant-in-Aid for Scientific Research (Grant No. 15340129) from the Ministry of Education, Science, Sports, and Culture of Japan.


text:
# Non-reciprocal optical reflection from a bidimensional array of subwavelength holes in a metallic film

## Abstract

Using simulations and theoretical arguments we investigate the specular reflection of a perforated gold film deposited on a glass substrate. A square lattice of cylindrical holes is assumed to produce the periodic lateral corrugation needed to hybridize the surface plasmons with radiative states. It is shown that, contrasting transmission approaches, a knowledge of the reflection on either side of the film provides separate information on the gold-vacuum surface plasmons and on the gold-glass interface plasmons. Recent experimental data on a specific implementation of this system are reexamined; these show a good agreement between the measured reflections and the simulations in both directions of incident wave probes. This confirms the importance of taking into account the reflection asymmetry in the far-field assessment of surface plasmons properties.

In recent years, subwavelength hole arrays in metallic films have received considerable attention, mainly because of the experimental findings of Thomas EbbesenT.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998). Such systems show a particularly intriguing zeroth order optical transmission : first, the specularly transmitted beam shows an intensity higher than that of the incident beam filtered by the geometrical aperture of the holes; second, the transmitted beam develops a strong characteristic wavelength dependenceT.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998). Recent theoretical and experimental studies have tried to clarify the detailed mechanisms involved in these phenomenaT.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998); H.F. Ghaemi, T. Thio, D.E. Grupp, T.W. Ebbesen, H.J. Lezec (1998); E. Popov, M. Nevière, S. Enoch, R. Reinisch (2000); Martin-Moreno et al. (2001); A. Degiron, H.J. Lezec, W.L. Barnes, T.W. Ebbesen (2002); M. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003); E. Altewischer, M.P. van Exter, J.P. Woerdman (2003); W.L. Barnes, W.A. Murray, J. Dintinger, E. Devaux, and T.W. Ebbesen (2004). In this context, it was rapidly recognized that surface plasmons (SPs) could be a key to the understanding of the transmission measurementsT.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998); H.F. Ghaemi, T. Thio, D.E. Grupp, T.W. Ebbesen, H.J. Lezec (1998); E. Popov, M. Nevière, S. Enoch, R. Reinisch (2000); Martin-Moreno et al. (2001); A. Degiron, H.J. Lezec, W.L. Barnes, T.W. Ebbesen (2002); M. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003); E. Altewischer, M.P. van Exter, J.P. Woerdman (2003); W.L. Barnes, W.A. Murray, J. Dintinger, E. Devaux, and T.W. Ebbesen (2004).

To be more precise, when the incident light illuminates the array of apertures, many diffraction orders are generated. Some of them are evanescent and exhibit a resonant coupling with SPs. They contribute to the zeroth diffraction order as a result of multiscattering and the strong exaltation of the electromagnetic field associated with resonant processes results in an enhanced transmissionT.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998); H.F. Ghaemi, T. Thio, D.E. Grupp, T.W. Ebbesen, H.J. Lezec (1998); E. Popov, M. Nevière, S. Enoch, R. Reinisch (2000); Martin-Moreno et al. (2001); A. Degiron, H.J. Lezec, W.L. Barnes, T.W. Ebbesen (2002); M. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003); E. Altewischer, M.P. van Exter, J.P. Woerdman (2003); W.L. Barnes, W.A. Murray, J. Dintinger, E. Devaux, and T.W. Ebbesen (2004). In a recent workM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003) we have suggested redefining the role of SPs in Ebbesen experiments and placing it in the required context of resonant Wood anomaliesWood (1935). In so doing, we have shown that the transmission spectrum could be better depicted as a series of Fano profiles. These recognizable lineshapes result from the interference of non resonant transfers with resonant transfers which involve the film eigenmodes and evanescent diffraction orders. We were then able to draw attention to the fact that each transmission maximum-minimum pair can be viewed as a Fano profileFano (1938). On quantitative grounds, the interplay of such an assymetric Fano profile actually dissociates the presence of a maximum or minimum in the transmission spectrum from the expected location of an eigenmode resonanceM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003).

It is somewhat surprising that reflection has not been as widely studied as transmission. One of our early contributions to this field described the theoretical reflection of an array of holes in a chromium filmM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003) and, at about the same time, Altewischer et al reported experimental reflection studies which emphasized its non-reciprocal propertiesE. Altewischer, M.P. van Exter, J.P. Woerdman (2003). As is well-known, there is a significant difference between the reflection spectra for light beams coming from one side or the other side of the device. Such difference is not observed in the transmission spectra, in agreement with the reciprocity theoremPetit (1980). The dissymmetry observed in the reflection spectra usually contains very significant information and both reflection and transmission should be thoroughly understood before a full knowledge of the optical properties of the system can be claimedM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003); E. Altewischer, M.P. van Exter, J.P. Woerdman (2003). The present paper proposes a theoretical investigation of the device studied by Altewischer et alE. Altewischer, M.P. van Exter, J.P. Woerdman (2003). This will provide an opportunity to establish the role of SPs in these experiments and, at the same time, relate the far-field transfer coefficients to the very local SPs dynamics.

The simulations in this paper are based on a coupled modes method which combines a scattering matrix formalism with a plane wave representation of the fieldsM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003). The inherently slow convergence of this representation is here accelerated by a Li algorithmLi (1997). This technique provides a computation scheme for the amplitude and polarization ($`s`$ or $`p`$) of reflected and transmitted fields in any diffracted orderM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003).

Fig.1 describes the system under study, chosen to model the samples experimentally investigated by AltewischerE. Altewischer, M.P. van Exter, J.P. Woerdman (2003). These samples can be described as a square lattice grating of circular holes perforated in a thin gold film deposited on a glass substrate. In accordance with experiments, $`100`$ nm is set for the radius of each hole and the lattice parameter is taken to be $`a=700`$ nm. The film thickness is matched to the experimental value of $`200`$ nm. The gold permittivity is obtained from experimental tablesD.W. Lynch, W.R. Hunter (1991). The substrate slab is modelled by a semi-infinite glass medium with a refractive index of $`\sqrt{\epsilon _s}=1.51`$E. Altewischer, M.P. van Exter, J.P. Woerdman (2003). Note that the experimental device exhibited a very thin bonding layer of titanium which separates the gold film from the substrateE. Altewischer, M.P. van Exter, J.P. Woerdman (2003). The vanishing thickness of this film justifies that it will not be considered in our calculations. The linearly polarized incident light is normal to the interface. Its electric-field vector is oriented along the $`[1,0]`$ direction, which joins the centers of nearest-neighboring holes.

Fig.2(a) shows the specular reflection as a function of the incident wavelength in both directions of propagation, i.e. an incidence from vacuum or from glass. The transmission is also shown, with a 5-times expansion of the vertical scaleE. Altewischer, M.P. van Exter, J.P. Woerdman (2003). The reflection presents strong variations near wavelengths $`757`$ nm $`(1^{})`$, $`812`$ nm $`(2^{})`$ and $`1146`$ nm $`(3^{})`$. These variations are fully consistent with experimental spectral features found at $`747`$ nm, $`810`$ nm and, presumably, at just over $`1100`$ nmE. Altewischer, M.P. van Exter, J.P. Woerdman (2003). Fig.2(b), 2(c), and 2(d) reproduces the curves shown on Fig.2(a), but exposes them to a comparison with the corresponding experimental data. Despite the necessary idealization step in designing the modelled theoretical device, we note an excellent agreement in all these comparisonsE. Altewischer, M.P. van Exter, J.P. Woerdman (2003). In particular, the relative locations of transmission maxima and reflection minima are correctly predicted. The only discrepancies which may still raise questions seem to be the weak oscillations found in all calculated spectra between $`850`$ nm and $`1050`$ nm, not found in experiment. The spectral line shapes are reminiscent of those described for similar systems in recent computationM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003) and experimentalW.L. Barnes, W.A. Murray, J. Dintinger, E. Devaux, and T.W. Ebbesen (2004) reports. The difference between the reflections on both sides of the film, i.e. the reflection nonreciprocity, is evident from the comparison of Fig.2(c) and Fig.2(d). These non-reciprocal properties of the reflection were also mentioned by Barnes et alW.L. Barnes, W.A. Murray, J. Dintinger, E. Devaux, and T.W. Ebbesen (2004) when describing their optical properties of structured metallic films.

In the present investigation the coupling between the various evanescent p-polarized diffraction orders and the surface plasmons is of prime importance. Because of the periodic corrugation of the film, the diffracted surface plasmons are sampled by the diffracted probe field at many frequencies. The resonant coupling of surface plasmons to p-polarized diffraction orders at matching frequencies and surface wave vectors will lead to what will be called ”resonant diffraction orders”. For an isolated surface plasmon which develops on a single interface, the wavelength which characterizes each resonant diffraction order is given byRaether (1988)

$$\lambda _{sp}^{u/m\text{ }(i,j)}=\frac{a}{\sqrt{i^2+j^2}}\text{Re}\sqrt{\frac{\epsilon _u\epsilon _m}{\epsilon _u+\epsilon _m}}.$$
(1)
In this equation, $`\lambda _{sp}=2\pi c\mathrm{/}\omega _{sp}`$ where $`\omega _{sp}`$ expresses the resonant angular frequency, $`a`$ is the grating parameter and $`(i,j)`$ denotes the vector $`𝐠`$ of the reciprocal lattice such that $`𝐠=(2\pi /a)(i,j)`$. $`\epsilon _m`$ is the metal permittivity and $`\epsilon _u`$ is either the permittivity of glass ($`\epsilon _s`$) or vacuum ($`\epsilon _v=1`$).

For the resonant orders $`(\pm 1,0)`$ we get at substrate/metal and vacuum/metal interfaces $`\lambda _{sp}^{s/m\text{ }(1,0)}=1084.6`$ nm and $`\lambda _{sp}^{v/m\text{ }(1,0)}=720.5`$ nm respectively. Both wavelengths are indicated by vertical dashes $`(3)`$ and $`(1)`$ respectively in the Fig.2 and Fig.3. For the resonant orders $`(\pm 1,\pm 1)`$ at substrate/metal interface, we get $`\lambda _{sp}^{s/m\text{ }(1,1)}=787.7`$ nm indicated by dash $`(2)`$. Dash $`(0)`$ corresponds to the resonant orders $`(\pm 2,0)`$ at substrate/metal interface for $`\lambda _{sp}^{s/m\text{ }(2,0)}=608.1`$ nm. This resonance appears in transmission as the first peak at $`610`$ nm.

The contrast between the reflection spectra for both incidence directions allows additional information to be obtained on processes involving specifically the surface plasmons at both interfacesE. Altewischer, M.P. van Exter, J.P. Woerdman (2003). For light arriving from the vacuum, two dips at the wavelengths $`(1^{})`$ and $`(2^{})`$ are observed. By contrast, for light arriving from the glass side a significant dip shows up at the wavelength $`(3^{})`$ while $`(1^{})`$ almost disappears. Altewischer et alE. Altewischer, M.P. van Exter, J.P. Woerdman (2003) suggested that such an asymmetry in the reflection spectra could lead to the identification of the interface on which the SPs are dominantly excited. We think that, furthermore, these results indicate that the excited resonances are mainly localized on the interface which is illuminated first and thus their effects are more easily perceived with reflection than with transmission.

To confirm and detail this interpretation, we have computed the amplitudes of main resonant diffraction orders as a function of wavelength and illumination directions. The field amplitude is considered here as a measure of the degree of excitation of a surface plasmon located on a specific interface. Fig.3a shows the local amplitude of resonant orders $`(\pm 1,0)`$ at the substrate/metal interface. When the light crosses the substrate/metal interface at the first place, it shows a peak amplitude at the wavelength $`(3)`$, predicted by Eq.1. One notes a second peak at wavelength $`(3^{})`$ close to $`(3)`$, possibly related to interferences of a true resonance with a non-resonant process, as found in Fano resonancesM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003). These peaks can easily be associated with the dip at similar wavelengths in the reflection spectrum. By contrast, when the vacuum/metal interface is illuminated first, the amplitude of these resonances collapses with the result that there is an attenuation of the associated dip in the reflection spectrum.

Fig.3b shows resonant orders $`(\pm 1,0)`$ on the vacuum/metal interface. Here the resonance develops at the wavelength $`(1)`$ and we note the existence of two other peaks at wavelengths $`(1^{})`$ and $`(2^{})`$. It is likely that these peaks arise from a surface plasmon splitting due to the coupling of the $`(\pm 1,0)`$ plasmon orders at the vacuum/metal interface with the $`(\pm 1,\pm 1)`$ plasmon orders at the substrate/metal interface. If the vacuum/metal interface is illuminated first, amplitudes of peaks $`(1)`$ and $`(1^{})`$ are stronger than in the opposite illumination. This explains the nonreciprocity found in the reflection spectra for wavelengths close to $`(1^{})`$. The spectral region close to the wavelength $`(2^{})`$, of the resonant order amplitude, does not significantly change when changing the propagation direction, leading to a weak nonreciprocity effect in the reflection spectrum. In Fig.3c we show the amplitude of orders $`(\pm 1,\pm 1)`$ at the substrate/metal interface and focus on the resonance wavelength noted $`(2)`$. A major contrast is obtained, with a large field amplitude under the direct illumination of the substrate/metal interface and a much weaker amplitude under illumination from the vacuum side.

The amplitude of the resonances of orders $`(\pm 1,\pm 1)`$ at the substrate/metal interface and $`(\pm 1,0)`$ at the vacuum/metal interface are influencing each other because of the proximity of their resonance wavelengths. Indeed, the profile of the orders $`(\pm 1,0)`$ amplitude at vacuum/metal interface, just after the wavelength $`(1)`$, can be understood as a Fano profile. Such a Fano resonance can be interpreted as resulting from the interference of the orders $`(\pm 1,0)`$ resonance at the vacuum/metal interface, and the orders $`(\pm 1,\pm 1)`$ resonance at the substrate/metal interface which then appears as a result of multiscattering of $`(\pm 1,0)`$ into $`(\pm 1,\pm 1)`$. This explains, for instance, why the wavelengths $`(1^{})`$ and $`(2^{})`$ do not strictly correspond to the resonance wavelength of the different resonant orders. Fano profiles resulting from the interference of each order as a result of multiscattering lead to a visible resonance wavelength shiftM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003).

We note that this analysis of the comparison between surface plasmon features and reflection corroborate the interpretation of Altewischer et alE. Altewischer, M.P. van Exter, J.P. Woerdman (2003), from which useful information about surface plasmons can be drawn from the study of the reflection nonreciprocity.

In addition to transmission and reflection, it might be useful to investigate the loss response function for both directions of illumination. This is shown in Fig.3d. As light comes from the vacuum side, we observe two sharp peaks at wavelength $`(1^{})`$ and $`(2^{})`$. These peaks correspond to the surface plasmon resonance already described above. On the other hand, when the light enters through the glass side, the first peak collapses, whereas another large double peak appears at the wavelength $`(3^{})`$, in perfect correlation with the surface plasmon excitation. This, again, is in full agreement with the experimental results of Barnes et alW.L. Barnes, W.A. Murray, J. Dintinger, E. Devaux, and T.W. Ebbesen (2004) and previous computationsM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003), confirming that both loss peaks correspond to a surface plasmon resonance. The loss line shapes can easily be correlated with the dominant features of the reflection spectra. These features indicate that, when SPs are excited, a large part of their energy is dissipated via the Joule effect, leading to lower reflection intensities. In addition, a smaller part of SPs energy is transferred to the transmitted beams as a result of multiscattering. So, as SPs are essentially excited on the interface which is encountered first, the reflection spectrum is dominated by processes associated with the losses from the plasmons carried by this interface.

As already noted by BarnesW.L. Barnes, W.A. Murray, J. Dintinger, E. Devaux, and T.W. Ebbesen (2004), this is quite different for transmission, where SPs of both interfaces equally contribute to enhanced transmission whichever side is illuminated first. The interface which receives the incident beam provides strong SPs resonances. Such resonances only give a weak contribution to the transmission since the associated field is damped when propagating through the holes before providing transmitted intensity. On the opposite interface, SPs are excited by an incident light weakened by its propagation through the holes, so that SPs on the opposite interface are only weakly excited. Their contribution to the transmitted field is thus also weak. In other words, SPs of the illuminated interface are strongly excited and then damped after propagation through the holes, whereas SPs on the opposite interface are weakly excited as the incident light is first damped. The result is that SPs of both interfaces contribute to the transmission with similar magnitudes and that the choice of the direction of propagation does not influence the information content of the measured transmission.

In conclusion, the examination of several experimental studiesE. Altewischer, M.P. van Exter, J.P. Woerdman (2003); W.L. Barnes, W.A. Murray, J. Dintinger, E. Devaux, and T.W. Ebbesen (2004) in the light of the interpretation scheme presented elsewhereM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003) reinforces the point of view expressed in this scheme by exhibiting qualitative and quantitative agreements between measured and computed light filtering spectra of a perforated gold layer covering a glass substrate. A particular effort was made to study the nonreciprocity property of the reflection. It was found that the asymmetry relative to the light propagation direction is effectively caused by the localization of SPs resonances at gold film surfaces and that this can be exploited to investigate separately the gold/glass and gold/vacuum plasmon coupling to radiative waves. The transfer functions computed in detail in this paper agree well with the quoted experimental data, available to us at the time of writing. One of the outcomes of these experiments is the role played by the corrugated-surface plasmons in the build up of the light transport process. In order to shed some light on this issue, the present work also looked into the theoretical value of the excited field at the location of the metal/vacuum and metal/glass interfaces and determined the degree of excitation of the interface plasmons. The strong correlation found between the plasmon amplitude and the reflection variations shows unambiguously the role played by the metal surface collective excitations in the light transfer process. This conclusion suggests that the experiments could be taken some steps further by considering complementary measurements, at least on the outer vacuum/gold surface, of the far field in reflectance experiments and the near field with the scanning near-field optical microscope. A correlation between these separate investigations should be highly instrumental in elucidating the mechanisms of the light transfer processes.

###### Acknowledgements.

We thank M.P. van Exter and E. Altewischer for enlightening advice and for making extensive experimental data available to us during the preparation of this work and P. Kelly for his advice on the manuscript. We acknowledge the use of the Namur Interuniversity Scientific Computing Facility (Namur-ISCF), a joint project between the Belgian National Fund for Scientific Research (FNRS), and the Facultés Universitaires Notre-Dame de la Paix (FUNDP). This work was partially supported by the EU Belgian-French INTERREG III project “PREMIO”, the EU5 Centre of Excellence ICAI-CT-2000-70029 and the Inter-University Attraction Pole (IUAP P5/1) on “Quantum-size effects in nanostructured materials” of the Belgian Office for Scientific, Technical, and Cultural Affairs.


text:
# Alternative pathways of dewetting for a thin liquid two-layer film

## Abstract

We consider two stacked ultra-thin layers of different liquids on a solid substrate. Using long-wave theory, we derive coupled evolution equations for the free liquid-liquid and liquid-gas interfaces. Depending on the long-range van-der-Waals forces and the ratio of the layer thicknesses, the system follows different pathways of dewetting. The instability may be driven by varicose or zigzag modes and leads to film rupture either at the liquid-gas interface or at the substrate. We predict that the faster layer drives the evolution and may accelerate the rupture of the slower layer by orders of magnitude promoting thereby the rupture of rather thick films.

Instability phenomena in ultra-thin soft matter films with thicknesses below $`100`$ nm became relevant mainly because they obstruct the fabrication of homogeneous coatings Reiter (1992). The interest was further boosted by the possibility to control such processes and, to use them to manufacture functional layers on the nanometer scale Mertig et al. (1998); Thiele et al. (1998). The stability of ultra-thin films is dominated by the effective molecular interactions between the substrate and the film surface Ruckenstein and Jain (1974). They represent, for instance, long-range van-der-Waals forces which increase (decrease) the pressure in the film if they are attractive (repulsive) Israelachvili (1992). However, to determine the emerging length scale and pattern for unstable films a study of the film dynamics is required. Using a film thickness evolution equation obtained by long-wave approximation Oron et al. (1997), the dewetting of a single layer of liquid is now reasonably well understood (see e.g. Ref. Thiele (2003)).

However, little is known on the behaviour of two stacked ultra-thin layers of simple or polymeric liquid on a solid substrate (see Fig. 1). Such a two-layer film allows for richer dynamics than an one-layer system because, both, the free liquid-liquid and the free liquid-gas interface evolve in a coupled way. The evolution is driven by the effective molecular interactions between all the three interfaces separating the four material layers: substrate, liquid1, liquid2 and ambient gas. Although experimental studies investigated different aspects of dewetting for two-layer films like interface instabilities or the growth of holes Faldi et al. (1995); Lambooy et al. (1996); Pan et al. (1997); Sferrazza et al. (1997, 1998); David et al. (1998); Renger et al. (2000); Morariu et al. (2003) up to now no general theoretical description of the interface dynamics has been given not . The case of small interface deflections was investigated in Ref. Brochard-Wyart et al. (1993) for a thickness of the lower layer, $`d_1`$, much larger than that of the upper one $`(d_2d_1)`$.

The most intricate questions for the first stage of dewetting of a two-layer system are which interface will become unstable, where does the film rupture, and how long will it take. This will determine the observability of the instability and the final morphology of the film. Experiments found roughening of the liquid-liquid interface Sferrazza et al. (1998) or an instability of the liquid-gas interface Faldi et al. (1995); Morariu et al. (2003). Holes that evolve solely in the upper layer were also studied Lambooy et al. (1996); Pan et al. (1997).

In this Rapid Communication, we derive and analyse coupled long-wave evolution equations for the two interfaces that are valid for all interface deflections and thickness ratios. We show that solely by changing the thickness ratio of the layers one switches between different dominant instability modes. This leads to drastic changes of the pathway of dewetting from rupture at the substrate to rupture at the liquid-liquid interface (see Fig. 2 below). Remarkably, for systems composed of two layers of very different thickness, i.e. with very different time scales for the rupture of the individual layers, the faster layer drives the evolution and accelerates the growth of surface modulations of the slower layer by orders of magnitude. We illustrate our results for two-layer systems of polystyrene (PS) and polymethylmetacrylate (PMMA) with silicon (Si) or silicon oxide (SiO) as substrate like studied experimentally in Refs. Lambooy et al. (1996); Sferrazza et al. (1997); Morariu et al. (2003).

We believe, our model can be used not only for the description of two-layer experiments with simple or polymeric liquids, but using an appropriate free energy functional also for a liquid film on a substrate with a stable but soft coating like a polymer brush Reiter et al. (1996) and, including driving terms, for studies of the transport of liquid droplets in liquid-liquid microfluidic systems Velev et al. (2003).

Coupled film thickness equations: We obtain evolution equations for the film thicknesses $`h_1`$ and $`h_2`$ by simplifying the Navier-Stokes equations employing long-wave approximation Oron et al. (1997). Thereby a no-slip condition at the substrate, the continuity of the velocity field and the balance of the stress-tensors at the liquid-liquid and liquid-gas interfaces are used. Considering an isothermal two-layer system where both layer thicknesses are smaller than $`100`$ nm, we neglect gravity and solely focus on the effective molecular interaction. For simplicity we only regard non-retarded long-range van-der-Waals forces resulting from dipole-dipole interactions between apolar materials. However, inclusion of other forces, like e.g. short-range polar forces Sharma (1993), or of slip boundary conditions Oron et al. (1997) (that may be necessary for polymer films) is straightforward as for one-layer films. Details of the derivation will be presented elsewhere. We obtain not

$`{\displaystyle \frac{h_1}{t}}`$ $`=`$ $`\left(Q_{11}{\displaystyle \frac{\delta F}{\delta h_1}}+Q_{12}{\displaystyle \frac{\delta F}{\delta h_2}}\right)`$
$`{\displaystyle \frac{h_2}{t}}`$ $`=`$ $`\left(Q_{21}{\displaystyle \frac{\delta F}{\delta h_1}}+Q_{22}{\displaystyle \frac{\delta F}{\delta h_2}}\right),`$ (1)
where $`\delta F/\delta h_i`$ with $`i=1,2`$ denotes functional derivatives of the total energy of the system

$$F=[\rho _\text{s}+\rho _{\text{VW}}]𝑑𝒙.$$
(2)
It contains the densities of the surface energy $`\rho _\text{s}=\frac{1}{2}[\sigma _1(h_1)^2+\sigma _2(h_2)^2]`$, and of the energy for the van-der-Waals interaction $`\rho _{\text{VW}}=A_{g21s}/(12\pi h_2^2)A_{21s}/(12\pi h_1^2)A_{12g}/[12\pi (h_2h_1)^2]`$. The surface tensions $`\sigma _1`$ and $`\sigma _2`$ belong to the liquid-liquid and liquid-gas interface, respectively. $`A_{g21s}`$, $`A_{21s}`$ and $`A_{12g}`$ are four- and three-index Hamaker constants, with subscripts s, 1, 2 and g refering to the substrate, liquid1, liquid2 and gas, respectively not .

The symmetric matrix of the positive mobility factors $`Q_{ik}`$ reads

$$𝑸=\frac{1}{3\mu _1}\left(\begin{array}{cc}h_1^3& \frac{3h_1^2}{2}\left(h_2\frac{h_1}{3}\right)\\ \frac{3h_1^2}{2}\left(h_2\frac{h_1}{3}\right)& \frac{(h_2h_1)^3(\mu _1\mu _2)}{\mu _2}+h_2^3\end{array}\right)$$
(3)
where $`\mu _1`$ and $`\mu _2`$ are the viscosities of liquid1 and liquid2, respectively. Note, that for $`d_2d_1d_1`$ and for small surface deflections Eqs. (1) simplify to those of Ref. Brochard-Wyart et al. (1993). Assuming two identical liquids, Eqs. (1) reduce to the well known one-layer equation Oron et al. (1997).

To compare to the well understood one-layer systems we non-dimensionalize Eqs. (1) using scales derived from the upper layer as an effective one-layer system. We scale $`𝒙`$ with $`\lambda _{\text{up}}=4\pi (d_2d_1)^2\sqrt{\pi \sigma _2/|A_{12g}|}`$, $`h_i`$ with $`d_2d_1`$ and $`t`$ with $`\tau _{\text{up}}=48\pi ^2\mu _2\sigma _2(d_2d_1)^5/A_{12g}^2`$. The corresponding energy scale is $`|A_{12g}|/16\pi ^3(d_2d_1)^2`$. The ratios of the mean thicknesses, surface tensions and viscosities are $`d=d_2/d_1`$, $`\sigma =\sigma _2/\sigma _1`$ and $`\mu =\mu _2/\mu _1`$, respectively. To compare with the lower layer as effective one-layer system one introduces in an analogous way the length scale $`\lambda _{\text{low}}`$ and time scale $`\tau _{\text{low}}`$.

We simulate the coupled time evolution of $`h_1`$ and $`h_2`$, Eqs. (1), in an one-dimensional domain using a semi-implicit time integration scheme and periodic boundary conditions. Initial conditions consist of flat layers with an imposed noise of amplitude 0.001. Alternative pathways of dewetting that occur for different thickness ratios $`d`$ are presented in Fig. 2 using a Si/PMMA/PS/air system as an example. Fig. 2 (a) shows that for a relatively small $`d=1.4`$ the two interfaces start to evolve deflections that are in anti-phase indicating the dominance of a varicose mode. When the liquid-gas interface approaches the liquid-liquid interface the latter starts to move downwards due to dynamical effects. This pathway leads to rupture of the upper layer, i.e. at the liquid-gas interface. On the contrary, Fig. 2 (b) shows that for a larger $`d=2.4`$ the growing deflections of the two interfaces are in phase indicating the dominance of a zigzag mode. As a consequence, here the lower layer ruptures, i.e. rupture occurs at the substrate.

Note, that in both, Fig. 2 (a) and (b), at the moment of rupture the respective non-ruptured layer is also in an advanced stage of its evolution leading to subsequent rupture. This is remarkable because their time scales as effective one-layer systems are 15 times (Fig. 2 (a)) and 35 times (Fig. 2 (b)) slower than the time scales for the respective fast layer. The ratio of the time scales $`\tau _{\text{up}}/\tau _{\text{low}}`$ is proportional to $`(d1)^5`$, i.e. for a lower layer ten times thicker than the upper one the rupture time of the lower layer is about five orders of magnitude larger than the one of the upper layer. However, a simulation for a Si/PMMA/PS/air system with $`d_1=10`$ and $`d_2d_1=1`$ shows that at rupture of the upper layer at $`t=0.61\tau _{\text{up}}=3.99\times 10^5\tau _{\text{low}}`$ the lower layer already evolved a depression of one fourth of its thickness. If the lower layer is the fast one, the effect also exists but is less pronounced.

In both cases, the acceleration of the rupture of the slower layer is caused by the direct coupling of the layers via the liquid-liquid interface. The fast evolution of the thinner layer deforms the interface and brings the thicker layer beyond the slow linear stage of its evolution. If the upper layer is the driving layer the process is in addition dynamically enforced because the liquid-liquid interface is ’pushed away’ by the advancing liquid-gas interface.

Linear stability analysis: Deeper understanding of the different pathways can be reached by studying the linear stability of the initial flat layers. We linearize Eqs. (1) for small disturbances $`\chi _i\mathrm{exp}(\beta t)\mathrm{cos}(kx)`$ for $`i=1,2`$, where $`k`$, $`\beta `$ and $`𝝌=(\chi _1,\chi _2)`$ are the wave number, growth rate and amplitudes of the disturbance, respectively. The dispersion relation $`\beta (k)`$ is obtained solving the resulting eigenvalue problem $`(k^2𝑸𝑬(k)+\beta 𝑰)𝝌=0`$, where $`𝑸`$ is the scaled mobility matrix and $`𝑬`$ is derived from the free energy Eq. (2) as $`E_{ij}=_{h_ih_j}\rho _{\text{V}W}+\delta _{ij}\stackrel{~}{\sigma }_ik^2`$ ($`\stackrel{~}{\sigma }_1=1`$, $`\stackrel{~}{\sigma }_2=\sigma `$, and $`\delta _{ij}=1`$ for $`i=j`$ and zero otherwise). The stability threshold shown as solid line in Fig. 3 is determined by $`𝑬`$ for disturbances of infinite wavelength, i.e. $`k=0`$. The system is linearly stable for

$`det𝑬>0\text{and}E_{11}>0\text{at}k=0.`$ (4)
An instability sets in if at least one of the conditions (4) is violated. This implies, that the two-layer film can be unstable even if both, $`_{h_1h_1}\rho _{\text{V}W}`$ and $`_{h_2h_2}\rho _{\text{V}W}`$ are positive, i.e. if the effective one-layer systems related to these terms are both stable.

Fixing the Hamaker constants, i.e. the combination of materials, and changing $`d`$ one finds a line (trajectory) in the stability diagram Fig. 3 as shown for a variety of experimentally studied systems. Interestingly, for van-der-Waals interactions calculated as detailed in not  one can show that such a trajectory can not cross the stability threshold, i.e. it is not possible to stabilize such a two-layer system by solely changing $`d`$. For instance, for the Si/PMMA/PS/air system the second condition in (4) is violated for all $`d`$ and the system is always unstable. At $`d=1`$, i.e. for a vanishing upper layer, the system is on the boundary between the one- and the two-mode regions (dashed line). For $`1<d<2.3`$ the unstable mode is an asymmetric varicose mode. A corresponding dispersion relation $`\beta (k)`$ is shown for $`d=1.4`$ in Fig. 4 (cp. the time evolution in Fig. 2 (a)). For $`d>2.3`$, i.e. for smaller thicknesses of the lower layer, the unstable mode is an asymmetric zigzag mode. Fig. 4 gives $`\beta (k)`$ for $`d=2.4`$ corresponding to the time evolution shown in Fig. 2 (b). For the fastest mode the zigzag mode is strongly asymmetric, i.e. the deflection of the liquid-liquid interface dominates the linear stage of the evolution. Note, that $`\beta (k)`$ and the dominant mode type depend on $`\sigma `$ and $`\mu `$, whereas the stability does not.

Further on, the simultaneous action of the van-der-Waals forces between the three interfaces allows for dispersion relations with two maxima. An experimental system showing this unusual form of $`\beta (k)`$ can be realized with a substrate that is less polarisable than both layers. This is the case for the SiO/PMMA/PS/air system not . A dispersion relation showing maxima of equal height is given for $`d=2.16`$ and $`\sigma =10`$ in Fig. 4. The maxima at small and large $`k`$ correspond to strongly asymmetric zigzag and varicose modes, respectively. This implies that the larger (smaller) wavelength will predominantly be seen at the liquid-gas (liquid-liquid) interface (Fig. 5 (a)). Increasing (decreasing) the ratio of the surface tensions strengthens the smaller (larger) wavelength. This implies that solely changing $`\sigma `$ by adding an otherwise passive surfactant one can switch from an evolution entirely dominated by the liquid-liquid interface to one dominated by the liquid-gas interface. This illustrates Fig. 5 by single snapshots from the non-linear time evolutions for different $`\sigma `$.

To conclude, we have derived coupled evolution equations for a thin liquid two-layer film driven by long-range van-der-Waals forces. The system represents the most general form of coupled evolution equations for two conserved order parameter fields in a relaxational situation and is apt to describe a broad variety of experimentally studied two-layer systems Refs. Thiele, John and Bär (2004). Linear and non-linear analysis have shown that the mobilities have no influence on the stability threshold, but determine the length and time scales of the dynamics. We have shown that for a two-layer system both interface deflection modes - zigzag and varicose mode - may be unstable and lead to rupture at the substrate or the liquid-liquid interface. Remarkably, the faster layer accelerates the evolution of the slower layer even if the latter is rather thick implying that its rupture time may be shortened by orders of magnitude.

Fig. 1
Pototsky et al., PRE

Fig. 2
Pototsky et al., PRE

Fig. 3
Pototsky et al., PRE

Fig. 4
Pototsky et al., PRE

Fig. 5
Pototsky et al., PRE


text:
# Origin of asymmetric reversal modes in ferromagnetic/antiferromagnetic multilayers

## Acknowledgments

The authors thank B. Beschoten and G. Güntherodt for intense discussions. This work has been supported by the Deutsche Forschungsgemeinschaft through SFB 491.