text:
# 1 Introduction

## 1 Introduction

D-brane models provide an interesting supersymmetry breaking mechanism that I would like to study in an effective field theory approach . Parallel D-branes of the same dimensionality break half of the supersymmetry that is present in the bulk. In the effective field theory on the D-brane world-volume, this is seen from the fact that the fields only fill multiplets of the smaller supersymmetry algebra. Part of the supersymmetry is explicitly broken on the D-branes since the world-volume fields lack the corresponding superpartners. Supersymmetry can be completely broken by adding anti-D-branes (e.g., ), which preserve the other half of the bulk supersymmetry. Thus, supersymmetry is broken in a non-local way in such models. Each sector taken separately preserves part of the supersymmetry. A very similar situation arises in models containing stacks of D-branes at angles that intersect each other (e.g., ). There is an extended supersymmetry on each stack of D-branes, but only a fraction of this supersymmetry is preserved at each intersection. Supersymmetry is completely broken in models where different intersections break different fractions of supersymmetry.

To determine the couplings of the bulk fields to the boundary fields, it is important to note that the part of supersymmetry that is broken on the D-branes is still non-linearly realized. Although a rigorous proof for this statement is still missing, there is much evidence in favor of it, including the following facts: (i) The $`𝒩=4`$ vector multiplet of a single D3-brane in flat ten-dimensional space contains just the right number of fermions to provide the goldstinos required for the non-linear realization of the bulk $`𝒩=8`$ supersymmetry. The counting still works if the bulk supersymmetry is reduced to $`𝒩=4`$ or $`𝒩=2`$ by an orbifold projection. (ii) Consider an $`𝒩=1`$ supersymmetric $`U(1)`$ theory. The requirement that the gaugino be the goldstino for a second non-linearly realized supersymmetry uniquely determines the full non-linear action. It turns out  that this action agrees with the supersymmetric generalization of the Born-Infeld action, which is known to describe the world-volume theory of D-branes. (iii) Consistent gravitino couplings are very constrained and it is hard to imagine how the bulk gravitinos corresponding to the broken supersymmetries can satisfy these constraints without non-linear supersymmetry .

In this talk, I want to discuss a four-dimensional toy model, where an $`𝒩=2`$ $`U(1)`$ bulk vector $`(A_m,\lambda ^{(1)},\lambda ^{(2)},\varphi )`$ couples to two boundary sectors that preserve different halves of the $`𝒩=2`$ bulk supersymmetry . On the first boundary, there is an $`𝒩=1`$ chiral multiplet $`(\varphi ^{(1)},\psi ^{(1)})`$ which carries charge $`q_1`$ under the bulk vector and an $`𝒩=1`$ goldstino multiplet $`(\lambda _g,\stackrel{~}{A}_m)`$. On the second boundary, there is an $`𝒩=1^{}`$ chiral multiplet $`(\varphi ^{(2)},\psi ^{(2)})`$ which carries charge $`q_2`$ under the bulk vector and an $`𝒩=1^{}`$ goldstino multiplet $`(\lambda _g^{},\stackrel{~}{A}_m^{})`$. The boundary multiplets transform linearly under one supersymmetry and non-linearly under the other supersymmetry. The bulk vector multiplet transforms linearly under both supersymmetries.

$`\delta \varphi ^{(1)}`$ $`=`$ $`\sqrt{2}\xi ^{(1)}\psi ^{(1)}2i\kappa \lambda _g\sigma ^m\overline{\xi }^{(2)}_m\varphi ^{(1)}`$
$`\delta A_m`$ $`=`$ $`i\lambda ^{(1)}\sigma _m\overline{\xi }^{(1)}i\lambda ^{(2)}\sigma _m\overline{\xi }^{(2)}+\mathrm{h}.\mathrm{c}.`$ (1.1)
$`\delta \varphi ^{(2)}`$ $`=`$ $`\sqrt{2}\xi ^{(2)}\psi ^{(2)}2i\kappa \lambda _g^{}\sigma ^m\overline{\xi }^{(1)}_m\varphi ^{(2)}`$
To determine the Lagrangian describing this toy model, I will first develop a method to couple $`𝒩=1`$ multiplets to non-supersymmetric matter and then generalize it to couplings of $`𝒩=2`$ multiplets to $`𝒩=1`$ matter. Using the results of  on partially broken $`𝒩=2`$ supersymmetry, the pseudo-supersymmetry Lagrangian is then easily obtained. A supergravity approach to pseudo-supersymmetry is pursued in .

## 2 Non-linear supersymmetry

Let us start by briefly reviewing the formalism of non-linearly realized supersymmetry. A supersymmetry transformation acts as a shift on the superspace coordinates:

$`x^m`$ $``$ $`x^m=x^mi(\xi \sigma ^m\overline{\theta }\theta \sigma ^m\overline{\xi })`$
$`\theta `$ $``$ $`\theta ^{}=\theta +\xi `$ (2.1)
$`\overline{\theta }`$ $``$ $`\overline{\theta }^{}=\overline{\theta }+\overline{\xi }`$
An important observation is that the goldstino $`\lambda _g`$ of broken supersymmetry can be viewed as a hypersurface in superspace defined through

$$\theta =\kappa \lambda _g(x),$$
(2.2)
where $`\kappa `$ is a constant of mass dimension $`2`$ related to the supersymmetry breaking scale. This is very similar to the Goldstone bosons of a spontaneously broken global internal symmetry being hypersurfaces in the parameter space of the global symmetry group.

The requirement that the hypersurface (2.2) be invariant under supersymmetry transformations, $`\theta ^{}(x)=\theta (x^{})`$, implies the standard non-linear transformation law

$$\delta _\xi \lambda _g=\frac{1}{\kappa }\xi i\kappa (\lambda _g\sigma ^m\overline{\xi }\xi \sigma ^m\overline{\lambda }_g)_m\lambda _g$$
(2.3)
for the goldstino. This forms a non-linear realization of the supersymmetry algebra

$$[\delta _\eta ,\delta _\xi ]=2i(\eta \sigma ^m\overline{\xi }\xi \sigma ^m\overline{\eta })_m.$$
(2.4)
A matter field $`f(x)`$ is well-defined on the hypersurface (2.2) if, under a supersymmetry transformation, one has $`f^{}(x)=f(x^{})`$. This implies

$$\delta _\xi f=i\kappa (\lambda _g\sigma ^m\overline{\xi }\xi \sigma ^m\overline{\lambda }_g)_mf,$$
(2.5)
which is the standard non-linear transformation for a matter field in the goldstino background.

To construct an invariant action, one introduces covariant derivatives

$`D_m`$ $`=`$ $`\left(\omega ^1\right)_m^n_n,`$ (2.6)
$`\omega _m^n`$ $`=`$ $`\delta _m^ni\kappa ^2(_m\lambda _g\sigma ^n\overline{\lambda }_g\lambda _g\sigma ^n_m\overline{\lambda }_g),`$
and notes that $`det(\omega )`$ transforms as a density under the non-linear supersymmetry . It is then straightforward to verify that

$$S=d^4xdet(\omega )\left(\frac{1}{2\kappa ^2}D_mfD^mfV(f)\right)$$
(2.7)
is invariant under (2.3), (2.5).

## 3 Coupling $`𝒩=1`$ to $`𝒩=0`$

Consider an $`𝒩=1`$ superfield

$$\mathrm{\Phi }(x,\theta ,\overline{\theta })=\varphi (x)+\theta \psi (x)+\mathrm{}$$
(3.1)
on which supersymmetry is linearly realized. To be able to couple the lowest component $`\varphi `$ of this superfield to a sector where supersymmetry is non-linearly realized, we need to find a composite field $`\widehat{\varphi }`$ that transforms according to (2.5) and reduces to $`\varphi `$ in the limit $`\kappa 0`$. It is easy to see that the desired composite field is given by

$$\widehat{\varphi }(x)\mathrm{\Phi }(x,\kappa \lambda _g(x),\kappa \overline{\lambda }_g(x)).$$
(3.2)
Indeed, using the invariance of the hypersurface (2.2) under linear supersymmetry transformations, one finds

$$\delta _\xi \widehat{\varphi }=\widehat{\varphi }(x+\delta x)\widehat{\varphi }(x)=i\kappa (\lambda _g\sigma ^m\overline{\xi }\xi \sigma ^m\overline{\lambda }_g)_m\widehat{\varphi }.$$
(3.3)
As an illustrative example, consider a chiral multiplet $`\mathrm{\Phi }`$. One has

$$\widehat{\varphi }=\varphi \kappa \lambda _g\psi +𝒪(\kappa ^2).$$
(3.4)
The supersymmetric completion of the dilaton-like coupling $`(\varphi ^{}+\varphi )F_{mn}F^{mn}`$ is thus given by

$$d^4x\left((\varphi +\varphi ^{}\kappa \lambda _g\psi \kappa \overline{\lambda }_g\overline{\psi })F_{mn}F^{mn}+𝒪(\kappa ^2)\right).$$
(3.5)
## 4 Partially broken $`𝒩=2`$

Let us see how the above formalism of non-linearly realized supersymmetry generalizes to partially broken $`𝒩=2`$ supersymmetry. The goldstino is now the lowest component of a superfield $`\mathrm{\Lambda }_g(x,\theta ,\overline{\theta })`$ with respect to the unbroken supersymmetry. I concentrate on the case, where the superpartner of the goldstino is a $`U(1)`$ gauge boson:

$$\mathrm{\Lambda }_{g\alpha }=\frac{i}{2}\left(\lambda _{g\alpha }+(\sigma ^{mn}\theta )_\alpha F_{mn}+\mathrm{}\right)=\frac{1}{2}W_\alpha +𝒪(\kappa ^2).$$
(4.1)
Under the second supersymmetry, $`\mathrm{\Lambda }_g`$ transforms as

$$\delta ^{(2)}\mathrm{\Lambda }_g=\frac{1}{\kappa }\xi ^{(2)}i\kappa (\mathrm{\Lambda }_g\sigma ^m\overline{\xi }^{(2)}\xi ^{(2)}\sigma ^m\overline{\mathrm{\Lambda }}_g)_m\mathrm{\Lambda }_g.$$
(4.2)
The full non-linear invariant action for the goldstino superfield has been determined by the authors of . They find

$$S=d^4x\left[\frac{1}{4}d^2\theta W^2+\frac{1}{4}d^2\overline{\theta }\overline{W}^2+\frac{\kappa ^2}{8}d^2\theta d^2\overline{\theta }W^2\overline{W}^2+𝒪(\kappa ^4)\right].$$
(4.3)
The bosonic terms of (4.3) coincide with the Born-Infeld action,

$$S_{\mathrm{bos}}=\frac{1}{\kappa ^2}d^4x\left(1\sqrt{det(\eta _{mn}+\kappa F_{mn})}\right).$$
(4.4)
There is a chiral version of the non-linear supersymmetry transformation which acts on chiral superfields $`\stackrel{~}{\mathrm{\Phi }}`$ as

$$\delta ^{(2)}\stackrel{~}{\mathrm{\Phi }}=2i\kappa \stackrel{~}{\mathrm{\Lambda }}_g\sigma ^m\overline{\xi }^{(2)}_m\stackrel{~}{\mathrm{\Phi }},$$
(4.5)
where $`\stackrel{~}{\mathrm{\Lambda }}_g(x,\theta ,\overline{\theta })=\mathrm{\Lambda }(x^mi\kappa ^2\mathrm{\Lambda }_g\sigma ^m\overline{\mathrm{\Lambda }}_g,\theta ,\overline{\theta })`$.

Defining $`\mathrm{\Phi }(x,\theta ,\overline{\theta })=\stackrel{~}{\mathrm{\Phi }}(x^m+i\kappa ^2\mathrm{\Lambda }_g\sigma ^m\overline{\mathrm{\Lambda }}_g,\theta ,\overline{\theta })`$, one finds that an invariant action is of the form

$$S=d^4x\left[d^2\theta d^2\overline{\theta }\widehat{E}\mathrm{\Phi }^{}\mathrm{\Phi }+d^2\theta E_L𝒫(\stackrel{~}{\mathrm{\Phi }})+d^2\overline{\theta }E_R𝒫(\stackrel{~}{\mathrm{\Phi }}^{})\right],$$

where

$`\widehat{E}`$ $`=`$ $`1+{\displaystyle \frac{\kappa ^2}{8}}\overline{D}^2\overline{\mathrm{\Lambda }}_g^2+{\displaystyle \frac{\kappa ^2}{8}}D^2\mathrm{\Lambda }_g^2+𝒪(\kappa ^4),`$
$`E_L`$ $`=`$ $`1+{\displaystyle \frac{\kappa ^2}{4}}\overline{D}^2\overline{\mathrm{\Lambda }}_g^2+𝒪(\kappa ^4),E_R=E_L^{}.`$ (4.6)
## 5 Coupling $`𝒩=2`$ to $`𝒩=1`$

The method to couple an $`𝒩=2`$ multiplet to $`𝒩=1`$ matter is a straightforward generalization of the method explained in section 3. In this talk, I will just give the result for the coupling of an $`𝒩=2`$ vector $`(A_m,\lambda ^{(1)},\lambda ^{(2)},\varphi )`$ to $`𝒩=1`$ matter. The $`𝒩=2`$ vector can be split into an $`𝒩=1`$ vector $`V=(A_m,\lambda ^{(1)})`$ and an $`𝒩=1`$ chiral multiplet $`\mathrm{\Phi }=(\varphi ,\lambda ^{(2)})`$.

Using

$`\delta ^{(2)}V`$ $`=`$ $`{\displaystyle \frac{i}{\sqrt{2}}}\theta \sigma ^m\overline{\theta }(\xi ^{(2)}\sigma _m\overline{D}\overline{\mathrm{\Phi }}^{}+\overline{\xi }^{(2)}\overline{\sigma }^mD\mathrm{\Phi }),`$
$`\delta ^{(2)}\mathrm{\Phi }`$ $`=`$ $`i\sqrt{2}W\xi ^{(2)},`$ (5.1)
$`\delta ^{(2)}W_\alpha `$ $`=`$ $`{\displaystyle \frac{i}{\sqrt{2}}}\overline{\xi }^{(2)}\overline{D}D_\alpha \mathrm{\Phi }{\displaystyle \frac{i}{2\sqrt{2}}}\xi _\alpha ^{(2)}\overline{D}^2\mathrm{\Phi }^{},`$
one finds that

$$\widehat{\mathrm{\Phi }}\mathrm{\Phi }+i\sqrt{2}\kappa \stackrel{~}{\mathrm{\Lambda }}_gW\frac{1}{4}\kappa ^2\stackrel{~}{\mathrm{\Lambda }}_g\stackrel{~}{\mathrm{\Lambda }}_g\overline{D}^2\mathrm{\Phi }^{}$$
(5.2)
transforms as $`\delta ^{(2)}\widehat{\mathrm{\Phi }}=2i\kappa \stackrel{~}{\mathrm{\Lambda }}_g\sigma ^m\overline{\xi }^{(2)}_m\widehat{\mathrm{\Phi }}`$. This implies that the $`𝒩=2`$ supersymmetric generalization of the dilaton-like coupling $`\mathrm{\Phi }𝒲^\alpha 𝒲_\alpha `$, where $`𝒲_\alpha `$ is an $`𝒩=1`$ Maxwell superfield without $`𝒩=2`$ partner, is given by

$$S=d^4xd^2\theta E_L\widehat{\mathrm{\Phi }}𝒲^\alpha 𝒲_\alpha +h.c.$$
(5.3)
Similarly, one finds that

$$\widehat{V}=V+i\kappa \theta \sigma ^m\overline{\theta }\left(\overline{\mathrm{\Lambda }}_g\overline{\sigma }_mD\mathrm{\Phi }+\mathrm{\Lambda }_g\sigma _m\overline{D}\mathrm{\Phi }^{}\right)+𝒪(\kappa ^2).$$
(5.4)
transforms as $`\delta ^{(2)}\widehat{V}=i\kappa \left(\mathrm{\Lambda }_g\sigma ^m\overline{\xi }^{(2)}\xi ^{(2)}\sigma ^m\overline{\mathrm{\Lambda }}_g\right)_m\widehat{V}`$. This implies that the $`𝒩=2`$ supersymmetric generalization of the gauge coupling term $`\mathrm{\Phi }_\mathrm{b}^{}e^V\mathrm{\Phi }_\mathrm{b}`$, where $`\mathrm{\Phi }_\mathrm{b}`$ is a chiral superfield without $`𝒩=2`$ partner, is given by

$$S=d^4xd^2\theta d^2\overline{\theta }\widehat{E}\mathrm{\Phi }_\mathrm{b}^{}e^{\widehat{V}}\mathrm{\Phi }_\mathrm{b}.$$
(5.5)
## 6 Pseudo-Supersymmetry

We are now in a position to write down the explicit Lagrangian for the pseudo-supersymmetry toy model described in the introduction. The field content of this model is summarized in table 1.

Note that the $`𝒩=2`$ vector can either be split into an $`𝒩=1`$ vector $`V`$ and an $`𝒩=1`$ chiral multiplet $`\mathrm{\Phi }`$ or into an $`𝒩=1^{}`$ vector $`V^{}`$ and an $`𝒩=1^{}`$ chiral multiplet $`\mathrm{\Phi }^{}`$. The chiral fields on the first boundary couple to $`V`$ whereas the chiral fields on the second boundary couple to $`V^{}`$. The Lagrangian is

$$=_{\mathrm{bulk}}+^{(1)}+^{(2)},$$
(6.1)
with

$`^{(1)}`$ $`=`$ $`{\displaystyle d^2\theta d^2\overline{\theta }\widehat{E}^{(1)}\mathrm{\Phi }^{(1)}e^{q_1\widehat{V}}\mathrm{\Phi }^{(1)}}+{\displaystyle d^2\theta E_L^{(1)}\stackrel{~}{\mathrm{\Lambda }}_g\stackrel{~}{\mathrm{\Lambda }}_g}+h.c.`$
$`_{\mathrm{bulk}}`$ $`=`$ $`{\displaystyle d^2\theta d^2\overline{\theta }\mathrm{\Phi }^{}e^V\mathrm{\Phi }}+{\displaystyle \frac{1}{4}}{\displaystyle d^2\theta WW}+h.c.`$ (6.2)
$`^{(2)}`$ $`=`$ $`{\displaystyle d^2\stackrel{~}{\theta }d^2\overline{\stackrel{~}{\theta }}\widehat{E}^{(2)}\mathrm{\Phi }^{(2)}e^{q_2\widehat{V}^{}}\mathrm{\Phi }^{(2)}}+{\displaystyle d^2\stackrel{~}{\theta }E_L^{(2)}\stackrel{~}{\mathrm{\Lambda }}_g^{}\stackrel{~}{\mathrm{\Lambda }}_g^{}}+h.c.`$
This supersymmetry breaking mechanism has several interesting consequences:

* There are no quadratic divergences at one-loop. Only Feynman diagrams involving fields from both boundaries can contribute to supersymmetry breaking. But such diagrams only arise at two-loop.
* The scalar masses squared arising at two loops are expected to be $`(g/4\pi )^4M^2`$, where $`g`$ is the gauge coupling of the bulk vector and $`M`$ is the cut-off scale of the effective field theory.
* A non-vanishing vacuum energy only arises at three-loop. (See  for a calculation of scalar masses and vacuum energy in a 5D version of the model described in this talk.)
* Both goldstinos stay massless (at tree-level) even when coupled to supergravity because the $`U(1)`$ gauge symmetries of their superpartners are unbroken.
* Gravitino masses arise only at three-loop.

Acknowledgements

I would like to thank Michael Peskin for many stimulating discussions, Cliff Burgess, Elise Filotas and Fernando Quevedo for collaborating on a related work which influenced my interest in pseudo-supersymmetry and Ignatios Antoniadis, Shamit Kachru and Albion Lawrence for answering several questions I had. Special thanks go to my wife for her support and encouragement. This research is funded by the Deutsche Forschungsgemeinschaft.


format: text

text:
# A study of random laser modes in disordered photonic crystals

## Abstract

We studied lasing modes in a disordered photonic crystal. The scaling of the lasing threshold with the system size depends on the strength of disorder. For sufficiently large size, the minimum of the lasing threshold occurs at some finite value of disorder strength. The highest random cavity quality factor was comparable to that of an intentionally introduced single defect. At the minimum, the lasing threshold showed a super-exponential decrease with the size of the system. We explain it through a migration of the lasing mode frequencies toward the photonic bandgap center, where the localization length takes the minimum value. Random lasers with exponentially low thresholds are predicted.

A finite open system of scattering particles can be characterized by a set of quasi-stationary (leaky) optical modes. When gain, sufficient to compensate the leakage in at least one mode, is introduced in such system - random laser is formed Letokhov (1968); Cao et al. (1999); Wiersma (2000); Frolov et al. (1999); Misirpashaev and Beenakker (1998); Burin et al. (2002a); Balachandran et al. (1997); Siddique et al. (1996); Wiersma and Lagendijk (1996). Unconventional feedback mechanism in a random laser leads to properties interesting from both fundamental and practical point of view Letokhov (1968); Cao et al. (1999); Wiersma (2000); Frolov et al. (1999); Misirpashaev and Beenakker (1998); Burin et al. (2002a); Balachandran et al. (1997); Siddique et al. (1996); Wiersma and Lagendijk (1996); Jiang and Soukoulis (2000); Vanneste and Sebbah (2001); John and Pang (1996); Burin et al. (2002b). While previous studies concentrated mainly on disordered systems, the resent experiments by Shkunov et al Var  demonstrated random lasing peaks associated with the bandgap in partially ordered system. In this letter, we obtain the lasing threshold for a random laser with different degree of ordering. Important information about these states can be extracted from the scaling of the average lasing threshold with the system size Burin et al. (2002b), $`L`$. For increasing values of the disorder parameter we find 5 scaling regimes: (a) photonic band-edge, $`1/L^3`$, (b) transitional super-exponential, (c) bandgap-related exponential, (d) diffusive, $`1/L^2`$, and (e) disorder-induced exponential, due to Anderson localized modes, regimes. In this letter we show that by optimizing the disorderness of the sample one can dramatically reduce the threshold of random laser, to the values comparable to those of photonic bandgap defect lasers Hig . This finding can open up a road to practical applications of random lasers.

Finding the lasing threshold of a random laser theoretically is a difficult problem. In disordered photonic crystals one cannot take simplifying assumptions such as independent scattering approximation (low density limit), neglect the finite size of the scatterers, or the confined dimensions of the system. Random laser modes have the highest $`Q`$’s and therefore weakly contribute to the transport properties of the system. Among other methods McPhedran et al. (1999); Nojima (2001), finite difference time domain (FDTD) method has been shown to be a convenient tool in studying random laser modes in 1D Jiang and Soukoulis (2000) and 2D Vanneste and Sebbah (2001). It was recently shown Vanneste and Sebbah (2001); Jiang and Soukoulis (2000) that above the threshold, the shape of the lasing modes remain the same as in passive medium. In this work, we used FDTD method to find the modes with the smallest energy decay rate in an open passive 2D random medium system with different degree of order. These modes determine the value of the lasing threshold in an active medium. After an initial pulsed excitation with the frequency $`\omega _e`$ in vicinity of the photonic bandgap, electric energy of the system followed $`Re[\mathrm{exp}2i\omega _m(1+i/2Q_m)]`$ dependence. $`\omega _m`$ and $`Q_m`$ were identified as the frequency and quality factor of the random cavity. In what follows, we will refer to $`1/Q_m`$ as the lasing threshold of the system. Our main results can be summarized as follows: (i) for sufficiently wide bandgaps, the minimum of the lasing threshold develops at some finite value of disorder parameter; (ii) at this “optimal” disorderness, the cavity $`Q`$ is determined by the localization length similar to that of single defect in the ordered sample, leading to a similar random cavity $`Q`$; (iii) with an increase of the system size the optimal disorder decreases; (iv) near the minimum, lasing threshold scales super-exponentially on sample size, which we ascribe to the migration of mode frequency toward the bandgap center and, therefore the decrease of the “typical” localization length.

We consider a 2D photonic crystal made of $`NL^2`$ cylinders with diameter $`d=98nm`$ and refractive index $`n_0=2.2`$. The cylinders were arranged into hexagonal lattice with nearest-neighbor distance of $`a=140nm`$. In the absence of disorder, the infinite system with these parameters has full bandgap in $`[361nm,426nm]`$ range for TM modes (electric field along the cylinder axis). The disorder in the system was introduced in 2 ways: by uniformly randomizing the refractive index of different cylinders in the range $`[n_0w_n(n_01),n_0+w_n(n_01)]`$ and diameter - $`[d(1w_d),d(1+w_d)]`$. Care was taken in order to avoid the uncontrollable disorder due to discretization of the grid. Disorder in the system was characterized with parameter $`\delta \epsilon =<(\epsilon (𝐫)\epsilon _0(𝐫))^2𝑑𝐫>^{1/2}/(\epsilon _0^2(𝐫)𝑑𝐫)^{1/2}`$, where $`\epsilon _0(𝐫)`$ and $`\epsilon (𝐫)`$ is the dielectric constant distribution in ordered and disordered samples respectively. $`<\mathrm{}>`$ stands for the average over different disorder configurations. In the present work we studied the systems with $`11`$ different disorder strengths: $`1`$ to $`10`$ had $`w_n[0.1,1.0]`$ with the increment $`0.1`$ and $`w_d0`$, the $`11`$th had $`w_n=1.0`$ and $`w_d=0.43`$. This led to variations of dielectric constant, $`\delta \epsilon `$, from weak - $`0.08`$ to strong - $`0.95`$ disorder. Later we will discuss the effects of this particular choice of these types of the disorder. We used standard FDTD method to study the evolution of TM electromagnetic field in the described system system with $`\lambda /20`$ resolution grid. In the studied range of system sizes, $`N[75,600]`$, the required grid sizes varied from $`200\times 200`$ to $`500\times 500`$. To mimic an open system, a buffer layer of air of $`150nm`$ thickness was kept around the sample, followed by uniaxial perfectly matched layer Taflove and Hagness (2000). To excite the system we initially launched a short $`10fs`$ pulse at every grid point. The frequency of the pulse was at the center of the bandgap, $`391nm`$, of the ordered structure. In the frequency domain the full width at half maximum of the seed pulse was of the order of the bandgap width. Thus the seed pulse excited all states within the stop band and near the band edges. After initial excitation pulse, the system was left on its own. Field distribution in the system and total electric energy $`=1/2\epsilon (𝐫)𝐄^2(𝐫)𝑑𝐫`$ was monitored.

Right after the initial pulse, the competition between the modes Misirpashaev and Beenakker (1998); Burin et al. (2001) with the different lifetimes led to the complicated evolution of $``$. However, after a sufficient time only the mode with the longest lifetime (maximum quality factor) survived. $``$ followed a mono-exponential dependence $`Re[\mathrm{exp}2i\omega _m(1+i/2Q_m)]`$, from which we extracted the maximum quality factor $`Q_m`$ and frequency $`\omega _m`$ for this particular realization of the disorder. At the same time the spatial pattern, $`𝐄(𝐫)`$, stabilized and the mode profile could be seen. Generally, the time needed to reach the mono-exponential regime varied from about $`0.5ps`$ for the smallest system to $`10ps`$ for the largest. Finally, the $`1/Q_m`$, representing lasing threshold, was averaged over $`1000`$ ($`N=75`$), or $`100`$ ($`N=137,188,261,368,449,608`$) disorder realizations.

Fig. 1 shows the dependence of the average lasing threshold normalized by $`N`$ as a function of the disorderness, different curves correspond to different system sizes. This particular normalization makes it easy to see the deviation from diffusion predicted Letokhov (1968) dependence $`1/L^2N`$. One can see that significantly different scaling at different $`\delta \epsilon `$ leads to the appearance of a minimum at the finite disorder strength. Fig. 2 plots $`L`$ dependence of the threshold for the fixed disorder strengths.

The understanding of this behavior comes from the position of the mode frequencies - Fig. 3. For small $`\delta \epsilon `$ the frequencies are concentrated at lower (long-wavelength) band-edge, and they (as well as $`Q_m`$) are independent on the position $`\omega _e`$ of the excitation pulse. The reason for that is the way the disorder was introduced into the system. The long-wavelength modes are mostly concentrated in the high dielectric medium, which is being disordered by the index fluctuations. At $`w_n=0.1`$ the mode frequencies fell in the immediate vicinity of the band-edge, Fig. 3a. Lasing from the band-edge modes is well studied in the case of the ordered structures Dowling et al. (1994); Sakoda et al. (1999); Kopp et al. (1998); Susa (2001); Ryu et al. (2002); Burin et al. (2002b), with $`1/L^3`$ dependence of the threshold Burin et al. (2002b) giving a good fit with our results. At the increased disorder, $`w_n=0.2`$, we observe the decrease of the threshold that is noticeably faster than exponential, expected for localized modes Burin et al. (2002a). This unexpected behavior becomes clear from Fig. 3b. The value of the threshold can be estimated as $`1/Q\mathrm{exp}(L/2\xi (\omega _m))`$, where $`\xi (\omega _m)`$ is some “typical” value of the localization length at the frequency $`\omega _m`$. From Fig. 3b one can see that with increased size, the frequency of the modes advance toward the band center, where $`\xi `$ is the shortest, magnifying the exponential dependence of the threshold. The explanation for this frequency migration with the increase of the size comes from the fact that in the finite system the modes deep into the bandgap are not accessible due to low density of states there. The Urbach-like behavior can be expected McPhedran et al. (1999); Frigerio et al. (1993). The support for this, indirectly comes from Fig. 3b, where $`1/Q`$ follows rather well the exponential dependence (solid line) on the frequency shift from the band-edge of the ordered system. However, because of the special nature of the modes no direct conclusions can be made. Assuming Urbach-like dependence of the density of states, the advancement of the mode frequency can be estimated from the condition that total number of states, $`N`$, times the probability of having a state located $`\mathrm{\Delta }\omega (N)`$ away from the band-edge, $`\mathrm{exp}(\alpha (\delta \epsilon )\mathrm{\Delta }\omega (N))`$, is still of order of one. Here, $`\alpha (\delta \epsilon )`$ is exponential slope in the density of states, that should decrease with the increase of the disorder $`\delta \epsilon `$. For small disorders $`\alpha ^1(\delta \epsilon )\mathrm{\Delta }E_{PBG}`$, where $`\mathrm{\Delta }E_{PBG}`$ is the width of the photonic gap. Therefore, for weak disorder (or small system size) the bandedge-type modes have the highest $`Q`$. The crossover to the super-exponential dependence of the threshold occurs when the cavity $`Q`$ associated with the localized states with the smallest localization length $`\xi (\mathrm{\Delta }\omega (N))`$ available for this size $`N`$ exceeds that of the bandedge-type mode: $`\mathrm{exp}[N^{1/2}a/2\xi (\mathrm{\Delta }\omega (N))]1/N^{3/2}`$. Stronger size dependence in the super-exponential regime means that even $`N^{3/2}`$ band-edge type dependence observed at smaller disorder would eventually switch to the super-exponential one. However, the latter can be expected to saturate at larger sizes or disorders when the mode frequency reaches the bandgap center: $`N\mathrm{exp}(\alpha (\delta \epsilon )\mathrm{\Delta }E_{PBG}/2)1`$, where the localization length is the smallest. Therefore, we expect the limiting scaling of the lasing threshold to be exponential $`\mathrm{exp}[N^{1/2}a/2\xi (\mathrm{\Delta }E_{PBG}/2)]`$. Three important conclusions (i-iii) listed in the introduction immediately follow. (iii) is justified because $`\xi (\mathrm{\Delta }E_{PBG}/2)`$ is the smallest in the ordered sample. The saturated exponential dependence can be seen in our example for $`w_n=0.3`$ \- squares in Fig. 2.

In our system the sharp rise of the threshold at $`\delta \epsilon 0.4`$ on Fig. 1 is attributed to the removal of the bandgap. This could be seen as the loss of the hexagonal symmetry of the observed modes as well a sensitivity of the modes to the excitation pulse position $`\omega _e`$. In this regime, the frequency of the modes is not associated with the photonic gap, which does not exist anymore. However, in order to make a direct comparison with the ordered case, we kept the excitation pulse the same as before. The exact cause of the disappearance of the bandgap is being debated in the literature Ban ; McPhedran et al. (1999); Frigerio et al. (1993), and is not the subject of this study. In our particular case we found a simple explanation of the threshold behavior in the way the disorder was introduced. Indeed, the fluctuating index of refraction leads to the fluctuation of the frequency of the Mie resonances of the particle. For box distribution, there exists a value of $`w_n=0.6`$ when the defect Mie resonance falls into the gap, Fig. 4. This value matches the value of disorder parameter, when the sharp increase of the lasing threshold is observed on Fig. 1. Moreover, Fig. 3c shows that at this crossover disorder, the modes avoid the region of high scattering, which is the consequence of the sharp boundary in the distribution of $`n`$. It also indicates the presence of the residual bandgap, where the lasing frequencies are concentrated.

At $`w_n0.7`$ the photonic bandgap ceased to exist. Triangles on Fig. 2 correspond to $`w_n=0.7`$, and were successfully fitted with diffusion Letokhov (1968); Balachandran et al. (1997); Siddique et al. (1996); Wiersma and Lagendijk (1996); John and Pang (1996); Burin et al. (2002b) scaling dependence $`1/L^2`$. Deviations from this dependence can be seen from Figs. 1,2 at the largest sizes studied, where $`L>\xi _{Anderson}`$ the states become localized again by virtue of Anderson localization Burin et al. (2002a, b). The exponential dependence becomes especially pronounced at the largest disorder studied (see inverted triangles on Fig. 2), where the transition from $`1/L^2`$ to exponential dependence comes at small system sizes. We want to point out that even at such strong disorder, the obtained modes had a collective nature, rather than the single particle’s (high-order resonances), which were concentrated at higher frequencies. Comparing the localization length of these states to that of bandgap nature we see the factor of two difference, which makes the latter preferable - Fig. 1. This is in accordance with experiments of Ref. Var .

In conclusion, we studied the modes responsible for the lasing threshold of a random laser based on disordered photonic crystal with a complete bandgap. At different disorder strengths, we identified 5 scaling regimes of the lasing threshold. Difference in scaling rates led to the development of the minimum at “optimal” disorder, where the threshold takes minimum value. We, for the first time, showed the regime of super-exponential scaling, which was attributed to the migration of the lasing frequencies to the center of the bandgap. We also predict that the value of the “optimal” disorder should decrease with increase of the volume of the system. These results can be generalized to the case of the 3D disordered photonic crystals.

###### Acknowledgements.

A. L. Burin is acknowledged for fruitful discussions. This work is supported by the National Science Foundation under Grant No DMR-0093949. HC acknowledges the support from the David and Lucile Packard Foundation and Alfred P. Sloan Foundation.


text:
# Anomalous thickness dependence of the Hall effect in ultrathin Pb layers on Si(111)

## Abstract

The magnetoconductive properties of ultrathin Pb films deposited on Si(111) are measured and compared with density-functional electronic band-structure calculations on two-dimensional, free-standing, 1 to 8 monolayers thick Pb(111) slabs. A description with free-standing slabs is possible because it turned out that the Hall coefficient is independent of the substrate and of the crystalline order in the film. We show that the oscillations in sign of the Hall coefficient observed as a function of film thickness can be explained directly from the thickness dependent variations of the electronic bandstructure at the Fermi energy.

The electronic properties of crystalline bulk materials are well understood. However, as soon as the lateral dimension in any direction is reduced to a few atomic distances, each electron energy band splits into a set of discrete subbands and the spatial confinement leads to quantum-size effects (QSE). The quantum-size confinement affects in particular the electron transport properties. In an early theoretical approach Trivedi and Ashcroft TA88  considered transport of electrons in a metal film, confined in a two-dimensional potential well and scattered by impurities in the film and/or by surface roughness. They predicted a saw-tooth like variation of the conductivity with the film thickness with a period of half the Fermi wavelength. Such oscillations have not been observed in the experiment so far. Another source of QSE-induced conductivity oscillations is the layer-by-layer growth of, e.g., Pb films on Si(111)$`(7\times 7)`$ or Si(111)-Au$`(6\times 6)`$ when the in-plane conductivity shows oscillations with the period of the interlayer spacing, associated with the degree of roughness of the surface.JBKL ; PH01  An even more spectacular, but so far unexplained, observation was the reversal of sign in the Hall coefficient $`R_H`$ of epitaxial Pb films on Si(111), grown on a surface that was pre-covered with about one monolayer (ML) of Ag which shows a $`(\sqrt{3}\times \sqrt{3})`$ reconstruction.JHB96

In this report, we address the general issue of charge transport and Hall effect in ultrathin metallic layers both theoretically and experimentally. We focus on Pb films deposited on Si(111), although the conclusions should be generally valid for ultrathin metallic films in the region of QSE. In order to clarify the underlying mechanisms, we combined in our study theoretical investigations by density-functional calculations with measurements of conductivity and of magnetoconduction of ultrathin Pb films on a $`(7\times 7)`$ reconstructed Si(111) substrate.

Pb films with a thickness between 1 and 12 monolayers (ML) have been deposited in UHV on a well annealed Si(111)$`(7\times 7)`$ substrate at low temperatures (15 to 25 K). The thickness, crystalline order and defects have been monitored with high-resolution low-energy electron diffraction (SPA-LEED).Pet01  4 Mo contacts have been pre-deposited for the measurements. After the pseudomorphic first Pb layer, up to 4 Pb layers were amorphous, whereas thicker layers always turned out to be crystalline. By varying temperature and thus the growth conditions, the influence of order on the conductive and magnetoconductive properties was tested. For details on film preparation and ordering in the films, see Ref. Pfen02, .

Before starting the Hall effect measurement the films have been annealed close to temperature where the decrease of conductance indicated a break-up of the film into non-contiguous islands. Even after annealing, films with a thickness $`d<4`$ML were strongly disordered whereas thicker films grew epitaxially with the lattice constant of bulk Pb. While the conductance measurements, which provide information on the charge-carrier scattering, were to some extent sensitive to the degree of order, the oscillatory behavior of the Hall coefficient turned out to be remarkably insensitive. Therefore, it seems to be justified as a first approach to compare the experimental results of the Hall coefficient with those obtained theoretically for free-standing Pb slabs. Indeed, as we will show below, such calculations are able to identify the basic origin of the oscillations.

Conductance and magnetoconductance measurements were performed in a DC magnetic field up to 4 T at various temperatures down to about 7 K. In the van-der-Pauw arrangement for conductance vdP58  it has been checked that the magnetoconductance was completely symmetric. In the contact arrangement for Hall effect there was a sizable contribution from magnetoconductance due to deviations from a perfect quadratic arrangements of the contacts. Therefore the asymmetry was taken as the contribution of the Hall effect to the measured voltage. The Hall coefficient has been derived from this asymmetry at $`\pm 4T`$.

In the model calculations we treat the Pb layers as free-standing slabs and neglect, as mentioned, the effect of the substrate. Within the model we calculate, as a function of layer thickness, the densities of states at the Fermi energies, Fermi velocities, and effective-mass tensors. From these quantities, together with the experimental data for the conductivity, $`\sigma _0^{\mathrm{exp}}`$, the elastic scattering times are calculated and the Hall coefficient, $`R_H`$ is derived.

The Hall coefficient in the low-field limit is

$$R_H=\frac{\sigma _H}{\sigma _0^2}$$
(1)
where $`\sigma _0`$ and $`\sigma _H`$ are the electrical and Hall conductivities of a two-dimensional slab. Adapting the three-dimensional expressions SAT92  to two dimensions we have

$$\sigma _0=\frac{e^2}{V\mathrm{}^2}\underset{n,k}{}\tau _n(k)\left[_kϵ_n(k)\right]^2\left(\frac{f(ϵ)}{ϵ}\right)$$
(2)
and

$`\sigma _H={\displaystyle \frac{e^3}{V}}{\displaystyle \underset{n,k}{}}\tau _n(k)^2\stackrel{}{v}_n(k)[\mathrm{𝐓𝐫}(𝐌^1)𝐌^1]\stackrel{}{v}_n(k)`$
$`\times \left({\displaystyle \frac{f(ϵ)}{ϵ}}\right).`$ (3)
$`𝐌^1`$ is the inverse mass tensor with the elements

$$(𝐌^1)_{i,j}=\frac{1}{\mathrm{}^2}\frac{^2ϵ_n(k)}{k_ik_j},$$
(4)
$`V`$ the slab volume, $`\stackrel{}{v}_n(k)=\frac{1}{\mathrm{}}\stackrel{}{}ϵ_n(k)`$ the electron group velocity, $`ϵ_n(k)`$ the energy of an electron with momentum $`k`$ in the subband $`n`$, and $`f`$ the Fermi function. The sums in Eqs. (2) and (3) run over the wavevectors $`k`$ in the two-dimensional hexagonal Brillouin zone and over the electron subband index $`n`$. (Each atomic electron level splits into $`d`$ two-dimensional electron subbands where $`d`$ is the number of monolayers in the slab.) The spin degeneracy is included in the prefactors. As will be seen below, the relaxation time $`\tau _n(k)`$ of thin slabs at low temperatures is dominated by elastic scattering on the thickness fluctuations of the Pb film (roughness of the interface to the substrate). Therefore, we will take the limit of short-wavelength surface height fluctuations (with their correlation length $`\xi 1/k_F`$, $`k_F`$ being the Fermi wavevector). In this case the relaxation time is $`k`$independent.C90  Since no detailed information on the scattering mechanisms and on the interface structure is available, we assume that $`\tau `$ is also independent of the subband index $`n`$. Then, $`R_H`$ is independent of $`\tau `$ and is determined solely by the electron bandstructure at the Fermi energy $`ϵ_F`$.

The electronic bandstructure calculations were performed on relaxed 1 to 8 monolayers thick (111) oriented Pb slabs with the in-plane lattice constant $`a=3.50`$ Å, separated by $`10`$ Å of vacuum. The electron band energies and the total energy have been calculated using the full-potential linearized augmented plane-wave method within the WIEN2k code BSL97  in the local-density approximation.PW92  A Pb muffin-tin radius of 2.6 a.u. and a tetrahedral mesh of 623 $`k`$points in the irreducible part of the full Brillouin zone were used in the self-consistent electronic structure calculations. The kinetic-energy cutoff was set to $`E_{\mathrm{max}}^{\mathrm{wf}}=7.3`$ Ry and the plane-wave expansion cutoff to $`E_{\mathrm{max}}^{\mathrm{pw}}=100`$ Ry. Later, in the magnetoconduction calculations, we used about 4000 $`k`$ points in the irreducible part of the two-dimensional Brillouin zone.

Several electron subbands cross the Fermi energy and there are several Fermi lines in the two-dimensional Brillouin zone. As typical examples, Fig. 1 shows the Fermi lines of 3 and 6 ML thick slabs. The number of Fermi lines is an increasing function of the slab thickness. Only electrons close to the Fermi lines take part in the transport and $`\sigma _0/\tau `$ and $`\sigma _H/\tau ^2`$ are obtained after numerical integration of expressions (2) and (3) along the Fermi lines. Whereas $`\sigma _0`$ is always proportional to the density of states at the Fermi lines , $`\sigma _H`$ of a subband is negative if the Fermi line makes a closed loop around occupied states (electrons) and positive if the loop encircles empty states (holes) in a subband. It is the balance between the number of electron and hole states at $`ϵ_F`$ together with their effective masses which determines the sign of $`R_H`$. As a consequence, $`R_H`$ is sensitive to the details of the bandstructure.

Fig. 2 shows the experimental conductivities $`\sigma _0^{\mathrm{exp}}`$ of annealed Pb films on Si(111)$`(7\times 7)`$ together with the data of epitaxial Pb on Si(111) covered with Ag.JHB96  Three properties of these data indicate that interface and/or surface scattering is the dominant mechanism of (elastic) electron scattering in this system. Firstly, only a weak temperature dependence of $`\sigma _0^{\mathrm{exp}}`$ is found below 100 K. Secondly, there is only a weak dependence of conductivity on crystallinity and order in the films. Annealing of amorphous films grown at 15 K resulted in an increase of conductivity, but this increase never exceeded a factor of 2.5. Thirdly, the roughly linear thickness dependence of $`\sigma _0^{\mathrm{exp}}`$ is again consistent with the assumption just made that the dominant scattering mechanism in these ultra-thin films is diffuse scattering at interface fluctuations.C90 ; CL99  For thicker samples the scattering at imperfections within the film becomes noticeable, and the conductivity starts to deviate from the linear behavior and eventually saturates in the bulk value, where the scattering on lattice defects and at higher temperatures also on phonons limit the conductivity.

We measured $`\sigma _0^{\mathrm{exp}}`$ for arbitrary film thickness $`d`$. At non-integer thicknesses the film surface is covered with islands and their borders (i.e., steps) act as scattering centers for the electrons. Therefore one would expect a dip in $`\sigma _0^{\mathrm{exp}}`$ for non-integer thicknesses. However, a small effect of this kind has only been observed at $`d>12`$ ML. This suggests that the typical island sizes exceed the charge-carrier mean free path for very thin layers and electron scattering at the vacuum interface is not important for films with $`d<10`$.

Clearly seen is an offset in the conductivity towards $`d=1`$ ML. The onset of conductivity due to percolation of the evaporated Pb at a surface temperature of 15 K was found to be at 0.8 ML. The very low conductivity of this first monolayer indicates that the main transport mechanism in this layer is hopping between localized states. Indeed, the conductivity shows a temperature-activated behavior for $`d=1`$ ML.Pfen02  The first Pb monolayer on Si(111)$`(7\times 7)`$ has the periodicity of the substrate Pet01  and may incorporate the adatoms of the $`(7\times 7)`$ substrate. This might be the reason for the very low conductivity and for the localization of charge carriers.

The close agreement of the conductivity data obtained in our measurements with those of Ref. JHB96,  is not at all self-evident, since our measurements have been carried out on Pb films that are expected to be much more disordered than those grown epitaxially on Si(111)-Ag$`(\sqrt{3}\times \sqrt{3})`$. Since the thickness dependence is almost the same for Pb films prepared on both substrates, we must conclude that interface scattering plays the same dominant role on both substrates. This means on the other hand that the degree of order in these ultrathin films is of secondary importance for electronic transport. We also note that the Ag interface layer obviously does not contribute to the measured conductance.

In contrast to the experimentally measured $`\sigma _0^{\mathrm{exp}}`$, the calculations of $`\sigma _0/\tau `$ and $`\sigma _H/\tau ^2`$ were performed only for integer $`d`$. Surprisingly, $`\sigma _0/\tau `$ turned out to be almost independent of the slab thickness in the investigated range between 2 and 8 ML,

$$\frac{\sigma _0}{\tau }1.1\times 10^{21}(\mathrm{\Omega }\mathrm{ms})^1.$$
(5)
With increasing number of monolayers the number of subbands increases and thus the conductance but not the ratio $`\sigma _0/\tau `$. This means that the contributions of newly emerging subbands to $`\sigma _0/\tau `$ as a function of layer thickness are very similar. Eq. (5) implies that $`\tau `$ is proportional to $`\sigma _0`$. Using the experimental values for $`\sigma _0`$, we determined $`\tau `$, also shown in Fig. 2.

It is interesting to look at the charge-carrier mean free path $`\lambda =v_F\tau `$. From the bandstructure we get the Fermi velocity $`v_F1.2\times 10^6`$ m/s and, using the values for $`\tau `$ from Fig. 2, we see that $`\lambda `$ is of the order 1 to 2 film thicknesses. It is therefore tempting to say that the electrons scatter on the film surfaces. However, the charge carriers scatter only on imperfections (steps, misplaced atoms etc.) on the surfaces and not on a perfectly flat surface. For $`d=1`$ ML the mean free path is shorter than the lattice spacing and the wave-vector description of the electron states breaks down.

The relaxation time causes broadening of the electron energy levels in the vicinity of the Fermi energy. Therefore, $`R_H`$ was calculated as $`\sigma _0/\tau /\sigma _H/\tau ^2`$ where $``$ is an average over the Gaussian distribution with the half-width $`\delta \mathrm{}/\tau `$.

The experimental results (open symbols) of the Hall coefficient are compared in Fig. 3 with those calculated (solid circles). Whereas the calculations have been done for perfectly periodic slabs without substrate, the measurements refer to epitaxial films with different degrees of disorder. Clearly seen are the anomalous oscillations of $`R_H`$ as a function of film thickness $`d`$ for both sets of experimental data. Our experimental results (squares) agree very well within the experimental error with the results of Ref. JHB96,  (open circles in Fig. 3). After having found close agreement for the data of electrical conductance, the agreement of the results for $`R_H`$ are not surprising, since the disorder within the film plays a negligible role. On the other hand, this result justifies our assumption of scattering times independent of $`k`$ and $`n`$.

$`R_H`$, calculated in this paper without any adjustable parameter, reproduces the oscillating behavior seen in the experiment semi-quantitatively and demonstrates that the main reason for this behavior are the thickness dependent changes in the band structure. The deviations at $`d=2`$ and $`d=5`$ may be a consequence of the approximation of equal relaxation times for all subbands. Especially in the case $`d=2`$ also the $`7\times 7`$ periodicity of the substrate, which was neglected in the slab calculation, might play a role. Our calculations also show that $`\sigma _0/\tau `$ is a slowly varying function of energy. Therefore it is not very sensitive to the exact location of the Fermi energy within the band structure, and the conductance increases smoothly with film thickness, and with only small variations during deposition. In contrast, the calculated $`\sigma _H/\tau ^2`$ is a strongly energy-dependent function and most of the energy dependence of $`R_H`$ comes from $`\sigma _H/\tau ^2`$, therefore $`R_H`$ is sensitive to averaging.

To conclude, we have shown that only a quantitative bandstructure analysis, done here with DFT, allows an understanding of the anomalous magnetoconductive properties of ultrathin Pb films. Oscillations in $`R_H`$ as a function of film thickness are the consequence of competing – and to a large extent compensating – contributions of electrons from different subbands and cannot be explained with any model which does not take into account the details of the electronic bandstructure. The relaxation time and thus the interface fluctuations play a subordinate role, they affect $`R_H`$ through the width of the distribution function and through the subband-dependence of the relaxation time, which was neglected in this paper. Although our calculations with free-standing films already show all qualitative features observed in experiment, it is obvious that this method could be further refined by explicitly taking into account the underlying substrate. This might help to understand the specific role of the metal-semiconductor interface considered here. Generally, this interplay between calculations and experiments is of primary importance since inner interfaces cannot be varied easily in experiments. In this context, also the open problem of the very low conductivity and the high $`R_H`$ of the first monolayer may be solvable. This involves most likely hopping conduction between localized states in the Pb film on reconstructed Si(111) substrate.

The authors are grateful to Prof. Matthias Scheffler for interesting discussions and Prof. Peter Saalfrank for sending us his results prior to publication. I.V. would like to acknowledge the hospitality of the Institut für Festkörperphysik in Hannover. This work was supported in part by the Deutsche Forschungsgemeinschaft.


text:
# Thermodynamics and magnetic field profiles in low-𝜅 type-II superconductors

## Abstract

Two-dimensional low-$`\kappa `$ type-II superconductors are studied numerically within the Eilenberger equations of superconductivity. Depending on the Ginzburg-Landau parameter $`\kappa =\lambda /\xi `$ vortex-vortex interaction can be attractive or purely repulsive. The sign of interaction is manifested as a first (second) order phase transition from Meissner to the mixed state. Temperature and field dependence of the magnetic field distribution in low-$`\kappa `$ type-II superconductors with attractive intervortex interaction is calculated. Theoretical results are compared to the experiment.

Superconductors are divided into two groups, type-I and type-II, based on their behavior in external magnetic field. The underlying physics is especially transparent within the Ginzburg-Landau (GL) model. It can be shown that the difference in thermodynamics originates from the sign of vortex-vortex interaction, which can be attractive or repulsive. The sign of the interaction is controlled by GL parameter $`\kappa =\lambda /\xi `$, ratio of penetration depth and coherence length near $`T_c`$. In type-I superconductors ($`\kappa <1/\sqrt{2}`$) vortices attract each other. In type-II superconductors ($`\kappa >1/\sqrt{2}`$), vortex-vortex interaction is repulsive and the phase transition from Meissner to the mixed state at $`H_{c1}`$ is of the second order.

It has been verified experimentally since a long time that away from $`T_c`$ even in type-II superconductors with $`\kappa 1`$ interaction between vortices can be attractive. Vortex attraction is manifested as a first order phase transition from Meissner to the mixed state. This type of superconductivity we denote as type-IIa to distinguish it from the type-IIb superconductivity with purely repulsive vortex-vortex interaction. Obviously attractive interaction in type-IIa superconductors will turn into repulsive below some critical intervortex distance thus preventing vortices to collapse into the normal state as in type-I superconductors. Evidence for jump in magnetization and mixture of vortex and Meissner domains near $`H_{c1}`$ is found from many magnetic, calorimetric, decoration and neutron scattering measurements (experimental results are partially reviewed in Ref. lukyanchuk ). Here we refer just to a nicely performed experiment by Auer and Ullmaier auerandullmaier  in which all three types (type-I, IIa and IIb) of superconductivity were reproduced in TaN by tuning the value of $`\kappa `$.

Although low-$`\kappa `$ type-II superconductors are of less interest from the point of view of potential application, attractive vortex interaction is still a challenging problem not fully understood. Most of the theoretical efforts to understand interaction of vortices as a function of parameter $`\kappa `$ relied on the extended GL model. Utilization of Bogomolnyi method to treat GL equations for $`\kappa `$ close to $`1/\sqrt{2}`$ gave the further impetus for this research direction lukyanchuk ; mohamed . However, range of applicability of GL formalism is confined to the region near $`T_c`$. Therefore to calculate the full temperature dependence of $`\kappa _c(T)`$, boundary between type-IIa and type-IIb superconductivity, one should rather turn to the microscopic theory of superconductivity. An approximate curve $`\kappa _c(T)`$ has been calculated by Kramerkramerzph  based on the asymptotic solution of Eilenberger equations of superconductivity. It is argued that the calculated curve is just the lower limit for the true one. Brandt brandtjump  has obtained jumps at $`H_{c1}`$ by variational approximation of his free energy functional. The only attempt to calculate numerically $`\kappa _c(T)`$ without approximations besides those immanent to numerical procedure is due to Klein klein . In this thorough study clean and dirty isotropic 3D superconductors are analyzed.

The purpose of this work is two fold. First, we applied a new method to solve Eilenberger equations for a vortex lattice. It appears to be efficient and faster than the so-called “explosion method”klein ; ichioka . Second, we extended previous numerical study of 2D superconductorsichioka1  to the low-$`\kappa `$ case where the interesting attractive vortex-vortex interaction occurs. Though the quantitative results will be different from the 3D case, the qualitative $`(\kappa ,T)`$ phase diagram is quite similar. This is explicitly shown by calculating phase diagram $`(\kappa ,T)`$ of type-I, type-IIa and type-IIb superconductivity. We focused on the magnetic and temperature dependence of calculated magnetic field profiles and compare the result with experimental data.

Eilenberger version of BCS theory of superconductivity is proven to be a good starting point for both analytical and numerical calculations. For the isotropic gap Eilenberger equations are

$$\left[\omega +𝐮\left(+\mathrm{i}𝐀\right)\right]f=\mathrm{\Psi }g,$$
(1)
$$\left[\omega 𝐮\left(\mathrm{i}𝐀\right)\right]f^{}=\mathrm{\Psi }^{}g.$$
(2)
These are supplemented by the self-consistency equations for the gap function $`\mathrm{\Psi }`$ and vector-potential $`𝐀`$

$$\mathrm{\Psi }\mathrm{ln}t=2t\underset{\omega >0}{}\left[f\frac{\mathrm{\Psi }}{\omega }\right],$$
(3)
$$\times \times 𝐀=\frac{2t}{\stackrel{~}{\kappa }^2}\mathrm{Im}\underset{\omega >0}{}𝐮g.$$
(4)
Equations are written in following dimensionless units: order parameter is measured in units $`\pi T_c`$, length in units $`R=v_0/(2\pi T_c)`$, $`v_0`$ is Fermi velocity, magnetic field in units $`H_0=\mathrm{\Phi }_0/2\pi R^2`$, vector potential in units $`A_0=\mathrm{\Phi }_0/2\pi R`$, energy in units $`E_0=(\pi T_c)^2N(0)R^3`$. Eilenberger parameter $`\stackrel{~}{\kappa }`$ is

$$\stackrel{~}{\kappa }^2=2\pi N(0)\left(\frac{\pi }{\mathrm{\Phi }_0}\right)^2\frac{v_0^4}{(\pi T_c)^2}.$$
(5)
It is related to GL parameter $`\kappa `$ via $`\stackrel{~}{\kappa }^2=(7\zeta (3)/18)\kappa ^2`$ in 3D case and

$$\stackrel{~}{\kappa }^2=\frac{7\zeta (3)}{8}\kappa ^2.$$
(6)
in 2D case. Here $`\zeta `$ is Riemmann’s zeta function. Here, $`\omega =t(2n+1)`$ is Matsubara frequency with integer $`n`$, $`t=T/T_c`$ is reduced temperature, $`𝐮`$ is unit vector directed along Fermi velocity. Eilenberger Green’s functions $`f`$, $`f^{}`$ and $`g`$ are normalized so that $`g=\sqrt{1ff^{}}`$. Average over the isotropic cylindrical Fermi surface $`<\mathrm{}>`$ reduces to $`(1/2\pi )\mathrm{}𝑑\phi `$, average over polar angle $`\phi `$. Expression for the free energy density difference between superconducting and normal state is given by, simplified with the help of Eilenberger and self-consistency equations

$$F=\stackrel{~}{\kappa }^2\overline{(\times 𝐀)^2}t\underset{\omega >0}{}\overline{\frac{1g}{1+g}(\mathrm{\Psi }^{}f+\mathrm{\Psi }f^{})}.$$
(7)
We use the following notation for spatial average $`\overline{C}=(B/2\pi )_{\mathrm{cell}}C𝑑S`$, where $`B`$ is magnetic induction. Apart from the extreme cases $`HH_{c1}`$ and $`HH_{c2}`$, which allow for an analytical solution, the only tool at hand is numerical analysis. Standard way to solve Eilenberger equations for a vortex lattice is so-called “explosion” method (for details see Refs. klein ; ichioka ). Here, another approach which takes the advantage of periodicity of the vortex lattice has been adopted. The same approach is used by Brandtbrandt1  to solve numerically GL equations. For this purpose it will be useful to introduce auxiliary functions $`a`$ and $`b`$ through the following transformation schopohlandmaki

$$f=\frac{2a\mathrm{exp}(\mathrm{i}\varphi )}{1+ab},f^{}=\frac{2b\mathrm{exp}(\mathrm{i}\varphi )}{1+ab},g=\frac{1ab}{1+ab}.$$

where $`\varphi `$ is the phase of the order parameter, $`\mathrm{\Psi }=|\mathrm{\Psi }|\mathrm{exp}(\mathrm{i}\varphi )`$. Equations for auxiliary functions $`a`$ and $`b`$ are decoupled and have the form of Riccati’s differential equation

$$𝐮a=\left(\omega +\mathrm{i}\mathrm{𝐮𝐀}_s\right)a+\frac{|\mathrm{\Psi }|}{2}\left(1a^2\right),$$
(8)
$$𝐮b=\left(\omega +\mathrm{i}\mathrm{𝐮𝐀}_s\right)b\frac{|\mathrm{\Psi }|}{2}\left(1b^2\right).$$
(9)
In new gauge vector-potential $`𝐀_s=𝐀\varphi `$ is proportional to the superfluid velocity. It diverges as $`1/r`$ at the vortex center (we put index $`s`$ to denote its singular nature). Functions $`a`$ and $`b`$ are not independent but there is a relationship between them

$$a^{}(𝒓,𝒖,\omega )=b(𝒓,𝒖,\omega ).$$
(10)
Therefore it is sufficient to solve only the equation for auxiliary function $`a`$. Also, we don’t need to solve equation (8) for all Fermi velocity directions because

$$a(𝒓,𝒖,\omega )=a(𝒓,𝒖,\omega ).$$
(11)
Reflection at the $`y`$ axis will change the direction of the magnetic field. When combined with time reversal it leaves Hamiltonian of the superconducting system invariant. One can show, at least numerically, that the following equation holds

$$a(x,y,u_x,u_y)=a^{}(y,x,u_x,u_y).$$
(12)
Therefore it is enough to solve Eq. (8) for velocity directions $`0<\phi <\pi /2`$ to recover $`a`$ and $`b`$ for all $`𝒖`$ via symmetry relations (10), (11), and (12).

Singular vector-potential has the following property

$$\mathbf{}\times 𝑨_s=𝑩(𝒓)2\pi 𝒆_z\underset{𝒓_i}{}\delta (𝒓𝒓_i),$$
(13)
where $`𝒓_i`$ are positions of vortices and $`\delta (𝒓)`$ is 2D Dirac function. For an ideal vortex lattice $`z`$-component of the magnetic field $`𝐁(𝒓)=\mathbf{}\times 𝑨`$ and moduo of the gap function $`|\mathrm{\Psi }|`$ has the periodicity of the vortex lattice. The same is not true for $`\mathrm{\Psi }`$ and $`𝑨`$. However, superfluid velocity, i.e. $`𝑨_s`$, is invariant under the translation from one to another vortex lattice cell. Therefore one must assume that auxiliary functions $`a`$ and $`b`$ are also periodic functions. Then we expand $`a`$ in Fourier series

$$a(𝒓)=\underset{𝑸}{}\underset{¯}{a}(𝑸)\mathrm{exp}(i𝑸𝒓),$$
(14)
where summation goes over reciprocal vortex lattice vectors $`𝑸`$. Fourier transform of a periodic function $`f(𝒓)`$ is defined as

$$\left[f(𝒓)\right]_{\mathrm{FT}}=\underset{¯}{f}(𝑸)=\frac{B}{2\pi }\underset{\mathrm{cell}}{}f(𝒓)\mathrm{exp}(i𝑸𝒓)d^2𝐫.$$
(15)
We split the local magnetic field $`𝐁(𝐫)`$ into constant part, magnetic induction $`\overline{𝐁}`$, and periodic part $`𝐁^{}(𝐫)`$ with zero average over the vortex lattice cell. Similarly, vector-potential is splitted into two parts

$$𝐁(𝐫)=\overline{𝐁}+𝐁^{}(𝐫)=\times \overline{𝐀}+\times 𝐀^{}.$$
(16)
We will assume that $`𝐀^{}(𝐫)`$ is periodic with $`𝐀^{}=0`$, so that all uncertainty of the gauge is in $`\overline{𝐀}`$. Therefore, $`𝐀^{}`$ is periodic solution of the self-consistency equation (4)

$$^2𝐀^{}=\frac{2t}{\stackrel{~}{\kappa }^2}\mathrm{Im}\underset{\omega >0}{}𝐮g.$$
(17)
From the equation (13) one can construct singular vector-potential $`𝐀_s`$ through the periodic solution $`𝐀^{}`$, $`𝐀_s=𝐀^{}\varphi ^{}`$ where kleininternal

$$\varphi ^{}=\underset{𝐫_i}{}\frac{(y_iy)𝐞_x+(xx_i)𝐞_y}{(xx_i)^2+(yy_i)^2}.$$
(18)
Sum converges very slowly and the rate of convergence can be improved by Ewald construction

$`\varphi ^{}=B{\displaystyle \underset{𝐐0}{}}{\displaystyle \frac{Q_y𝐞_x+Q_x𝐞_y}{Q^2}}\mathrm{exp}(\gamma Q^2)\mathrm{sin}(\mathrm{𝐐𝐫})+`$
$`{\displaystyle \underset{𝐫_i}{}}{\displaystyle \frac{(y_iy)𝐞_x+(xx_i)𝐞_y}{(xx_i)^2+(yy_i)^2}}\mathrm{exp}\left[{\displaystyle \frac{(𝐫𝐫_i)^2}{4\gamma }}\right].`$ (19)
Parameter $`\gamma `$ is in principle arbitrary; we choose $`\gamma =\pi `$ to assure good convergence for both sums.

Self-consistency equations are solved by iteration. For a given vector-potential $`𝑨^{}`$ and gap function $`\mathrm{\Psi }`$, Eq. (8) is solved numerically. We start the iteration procedure with $`𝑨^{}=0`$ and $`\mathrm{\Psi }=\mathrm{const}`$. Functions $`a`$ and $`b`$ are then used to calculate new values of $`𝑨^{}`$ and $`\mathrm{\Psi }`$ with help of self-consistency equations, and the whole process is repeated until the self-consistency is achieved. Equation (8) is also solved by iteration. We solved it in Fourier space where unknown are coefficients $`\underset{¯}{a}`$. Iterative procedure for unknown coefficients $`\underset{¯}{a}`$ is

$$\underset{¯}{a}(𝐐)=\frac{\left[|\mathrm{\Psi }|(1a^2)/2\mathrm{i}\mathrm{𝐮𝐀}_sa\right]_{\mathrm{FT}}}{\omega +\mathrm{i}\mathrm{𝐐𝐮}}.$$
(20)
Maximum summation frequency was set as $`\omega _c=20t`$ ($`\mathrm{}\omega _c=20\pi k_BT_c`$ in physical units). Vortex lattice unit cell is divided into $`N\times N`$ mesh which turns into approximating auxiliary functions $`a`$ and $`b`$ by sum over $`N\times N`$ harmonics. We set $`N=128`$. Integration over polar angle $`\phi `$ in each quadrant is performed by 30-points Gauss-Legendre integration. Iteration procedure fails for extremely small $`\omega `$, i.e. for very low temperatures $`t0.05`$. In all other cases it is recommended to damp the iteration process via $`a_{new}^{}=(c_1a_{new}+c_2a_{old})/(c_1+c_2)`$ with $`c_1c_2`$. To manipulate with eqs. (14) and (15), fast Fourier transform algorithm should be used.ooura

Now we concentrate on type-II superconductors with $`\kappa 1`$. For a particular value of $`\kappa `$ and $`t`$ set of Eilenberger equations are solved numerically for various values of magnetic induction $`B`$, so we can construct free energy density $`F(B)`$. However, in the experiment we actually control applied field $`H`$ and not induction $`B`$. For a given applied filed $`H`$ proper thermodynamic potential is Gibbs energy $`G(\overline{𝐁},T)`$ which is at minimum when the thermodynamic equilibrium is achieved, $`G/B=0`$. This condition is rewritten via Helmholtz free energy $`F(B,T)`$

$$H(B)=\frac{1}{2\stackrel{~}{\kappa }^2}\frac{F}{B}.$$
(21)
If there are meta-stable states, as we expect in type-IIa superconductors, function $`B(H)`$ is multi-valued and to find a true magnetic induction one should look for global minimum of Gibbs energy. In Fig. 1 we plotted typical example of magnetization curve with meta-stable states, i.e. for some values of $`H`$ there are three different values of magnetic induction that correspond to the extrema of Gibbs energy. To find $`B`$ which corresponds to the minimum energy one has to construct Gibbs energy density:

$$G=F2\stackrel{~}{\kappa }^2BH.$$
(22)
For example, we found that the minimum Gibbs energy for $`H=0.34`$ correspond to the highest value of induction $`B=0.20`$. Therefore from the behavior of Gibbs energy we can reconstruct the $`B(H)`$ curve. From $`H=0`$ up to $`H=H_{c1}^{}`$ magnetic induction is zero and we have Meissner state. Note that $`H_{c1}^{}`$ is less than $`H_{c1}=\underset{B0}{lim}H(B)`$. At $`H=H_{c1}^{}`$ there is first order phase transition from the Meissner state $`B=0`$ to the mixed state $`B=B_0`$ (jump in magnetization). It means that energy of well separated vortices is bigger than vortex-lattice energy at $`B=B_0`$, due to attraction of vortices. For a fixed temperature the value of induction jump $`B_0`$ decreases with increasing $`\kappa `$ and at some critical value $`\kappa _c`$, $`B_0=0`$. Then for all $`\kappa >\kappa _c`$ we have a second order phase transition from the Meissner to mixed state, i.e. type-IIb superconductivity.

Calculations are performed for various temperatures $`t=0.20.9`$ to find the critical parameter $`\kappa _c`$ that separates type-IIa and type-IIb superconductivity. The results are shown on Fig. 2. Boundary line is qualitatively similar to the 3D case obtained by different numerical scheme. Boundary between type-I and type-II superconductivity can be obtained from the condition $`H_c=H_{c2}`$. Calculation itself is trivial and we did it merely to show the quantitative difference between 3D (shown by dashed line in Fig. 2) and 2D case. Compared to 3D case, boundary between type-I and type-II superconductivity in 2D superconductors is almost constant and changes only few percent down to $`T=0`$. In 3D case there is possibility for cross-over from type-I to type-II superconductivity as one goes from $`T_c`$ to $`T=0`$.

We obtained selfconsistent solution for the order parameter $`\mathrm{\Psi }(𝒓,t)`$ and field distribution $`𝑨(𝒓,t)`$ for various $`\kappa `$, $`B`$ and temperature $`T`$. Only one interesting aspect will be presented: magnetic field profile as viewed via nuclear magnetic resonance (NMR) and muon spin rotation ($`\mu `$SR). These are standard experimental techniques that probe the magnetic field profile in the mixed state. NMR absorption intensity in a small frequency interval $`\mathrm{\Delta }\omega `$ of the ac field is proportional to the fraction of the VL cell area where the field value is in the corresponding interval $`\mathrm{\Delta }\omega \mathrm{\Delta }H`$, i.e. signal is proportional to $`\delta (HH(𝒓))d𝒓`$. NMR line shape has the logarithmic singularity at the saddle point of magnetic field. In $`s`$-wave superconductor from the relative position of the NMR peak one can make conclusion about the VL structure.

Kung kung  measured NMR profile in V foil as a function of temperature. The applied field was perpendicular to the foil plane. In this geometry magnetic induction was fixed and almost temperature independent. It was found that the maximum magnetic field (field at the vortex center) is linear in $`T`$ down to the lowest attainable temperature $`0.1T_c`$ in that experiment. Delrieu delrieu  calculated magnetic field profile in the limit $`T=0`$ and small magnetization. He found singularities at the positions of field maximum and minimum which have conical shape compared to parabolic shape near $`T_c`$. It was shown that the slope $`\mathrm{d}H_{max}/\mathrm{d}T`$ is finite at $`T=0`$.

We calculated temperature dependence of $`H_{max}`$ in superconductor with $`\kappa =0.8`$ for fixed magnetic induction $`B=0.16`$. The result is shown in Fig. 3. Overall behavior $`H_{max}(T)`$ resemble the linear behavior reported in kung . The result is in agreement with the calculation of Delrieu if the curve is extrapolated to $`T=0`$. We note that in Ref. klein  $`H_{max}`$ has a more pronounced convex curvature. Therefore, in spite the value of magnetic penetration depth is saturated already below $`0.4T_c`$, magnetic field distribution is temperature dependent down to very low temperatures.

In summary, we presented the new method to solve selfconsistently Eilenberger equations in the mixed state. Its efficiency is demonstrated in the problem of low-$`\kappa `$ 2D superconductors. Being faster then previous techniques, we hope it will stimulate further numerical studies of type-II superconductors.


text:
# A phase transition due to thick vortices in SU(2) lattice gauge theory

## 1 Introduction

The role of thick vortices  as a mechanism for confinement in $`SU(2)`$ LGT has been extensively discussed in the last few years . Thick vortices are analogous to the domain walls in ferromagnets with a continuous symmetry–they are like thick Peierls contours. Thick vortices should be distinguished from thin vortices. Thin vortices have infinite action in the continuum limit. For this reason they are unlikely to play a major role in the continuum limit. On the other hand thick vortices can reduce their free energy by arbitrarily increasing their cross section . In addition to thin and thick vortices $`SU(2)`$ gauge theory also has thin and thick monopoles. A thin vortex can end in a thin $`Z(2)`$ monopole which (in three space-time dimensions) is defined on an elementary 3 dimensional cube. The thin $`Z(2)`$ monopole density on a cube $`c`$ is given by

$$\rho _1(c)=(1/2)(1sign(\underset{pc}{}trU(p))).$$
(1)
Analogously, a thick vortex can end in a thick $`Z(2)`$ monopole. Its density (in three space-time dimensions) is given by

$$\rho _d(c_d)=(1/2)(1sign(\underset{dc_d}{}trU(d))),$$
(2)
where the product is taken over all $`dXd`$ loops bordering a cube $`c_d`$ of side $`d`$. The subscript $`d`$ indicates that the density can be defined on any 3 dimensional cube of side $`d`$. In three space-time dimensions the monopoles are point like objects or they have a finite size. Both monopoles are $`SO(3)`$ invariant as they only depend on the $`SU(2)/Z(2)`$ coset of the link variables. A whole hierarchy of such monopoles can be defined with different thicknesses. Unlike the thin monopoles, which are suppressed by the action, the thick monopoles cost lesser energy because their energy is spread over a finite region. It was proposed in  that at lower and lower temperatures vortices and monopoles of thicker and thicker cross-section are important. These ideas were encapsulated in the effective $`Z(2)`$ theory of confinement  where it was proposed that the long distance properties of $`SU(2)`$ LGT can be described by a $`Z(2)`$ theory with an effective coupling $`\beta (d)`$ ( $`\beta (d)\mathrm{}asd\mathrm{}`$). Thick vortices are the vortices in the effective $`Z(2)`$ theory.

The presence of thick monopoles prevents the thick vortices from forming closed loops and makes it difficult to study their effects separately. If the thick monopoles are removed it might be possible to see if thick vortices are present at all and then study their effects. With this aim in mind we study the $`SU(2)`$ theory after eliminating the thin and thick monopoles. This model can be regarded as a generalisation of the Mack-Petkova model which eliminates only the thin monopoles. Just as in the Mack-Petkova model the elimination of the thin and thick monopoles changes only the short distance properties of the gauge theory.

Fig.1 The $`2X2`$ Wilson loop in the $`SU(2)`$ theory.

## 2 The model

Thin and thick monopoles can be removed by introducing two (large) chemical potential terms for the thin and thick monopole densities in the Wilson action. The action considered is

$$S=\frac{\beta }{2}\underset{p}{}trU(p)\lambda _1\underset{c_1}{}\rho _1(c_1)\lambda _2\underset{c_2}{}\rho _2(c_2).$$
(3)
$`\lambda _1,\lambda _2`$ are chosen to be large large so that the monopole densities are very small ( in practice $`\lambda _1=\lambda _220`$ resulted in identically zero densities for almost all thermalised configurations on $`12^3`$ lattices).

The simulation of this model presented its own share of difficulties. The environment of a single link is complicated because each link touches four thin monopoles and eight thick monopoles. Different strategies were attempted to simulate this model. Metropolis updating and a combination of heat bath and metropolis were both tried but metropolis updating was found to be more efficient provided the table of $`SU(2)`$ elements is tuned regularly to get a reasonable acceptance. One observes long metastabilities while thermalising the lattice and the simulation has to be kept running for unusually long times unless a good starting configuration is chosen. It was found that in the phase where thick vortices are in abundance an ordered start (which has zero density of monopoles and vortices) took a very long time to reach the equilibrium distribution. Similarly, in the phase where thick vortices are absent a random start ( which has a large number of thick vortices) took a long time to reach the equilibrium distribution.

Fig2. The plaquette in the monopole suppressed model.

These metastabilities donot arise in the pure $`SU(2)`$ theory or the Mack-Petkova model and they appear to be linked to the thick vortices present in this model. Several approaches were tried to deal with this metastability like starting from configurations ”closer” to the equilibrium configuration but the problem always remained. Nevertheless, with increased computational effort consistent results can be and were obtained. The simulation was performed in three dimensions mainly for reasons of computational speed. On the lattices used ($`12^3`$) it was found that simulating this model is more time consuming than simulating the four dimensional pure $`SU(2)`$ model.

## 3  Numerical results

It is well known that 3 dimensional $`SU(2)`$ LGT does not exhibit the crossover found in the four-dimensional theory but instead has a very smooth behaviour for the plaquette and other quantities. In Fig 1 we show this smooth behaviour in the $`2X2`$ Wilson loop. Eliminating thin monopoles leads to a slight jump in the plaquette due to the presence of thin vortices. Eliminating the thick monopoles leads to another jump at a lower temperature. Our study has focused on the properties of the system around this point. Fig 2 and Fig 3 show the sharp jump in the plaquette and the $`2X2`$ Wilson loop respectively. In the vicinity of this jump very long relaxation times are observed before a starting configuration reaches equilibrium. The sign of the $`2X2`$ Wilson loop changes abruptly across this transition. It should be emphasized that the thin vortices have all disappeared at this point (the sign of the plaquette is always a positive quantity here) and the fluctuations present are over longer distances.

Since $`Z(2)`$ fluctuations at longer length scales are observed it is natural to ask if these fluctuations are caused by thick vortices. Indeed they are! Thick vortex loops can be identified as closed loops which are composed of a co-closed set of $`2X2`$ plaquettes with a negative sign. A measurement of the thick vortices on either side of the transition shows that they form very long loops in one phase and smaller loops in the other phase. Hence, the thick vortices are producing the $`Z(2)`$ fluctuations at longer distances, in line with the ideas of the effective $`Z(2)`$ theory of confinement. The location of the phase transition was determined to be around $`\beta 2.5`$.

Fig3. The sign of the $`2X2`$ Wilson loop in the monopole suppressed model.

## 4 Conclusions

Thick vortices can be directly observed once the thick monopoles are eliminated. The condensate of thick vortices persists at higher values of $`\beta `$ after the thin vortices have all disappeared. They are responsible for the long range $`Z(2)`$ fluctuations which are a crucial ingredient in the effective $`Z(2)`$ theory of confinement. As $`\beta `$ increases the lattice spacing becomes smaller and vortices with a non-zero physical thickness can still be present. In the lattice model there will be a hierarchy of effective $`Z(2)`$ theories operating at larger and larger $`\beta `$ ; the model studied here is just the first member of the hierarchy. Eliminating monopoles of greater thickness will in principle unravel the entire hierarchy of $`Z(2)`$ like theories operating at larger lattice spacings . It is a challenge to show that in the continuum limit we are left with $`Z(2)`$ fluctuations at physical distances of the order of a fermi as predicted in .


text:
# Transverse acoustic nature of the excess of vibrational states in vitreous silica

## Abstract

We present a numerical simulation study of the density-dependence ($`\rho `$=2.2$`÷`$4.0 g/cm3) of the high-energy collective dynamics in vitreous silica at mesoscopic wavevectors ($`Q`$=1$`÷`$18 nm-1). The density-dependence of the longitudinal and transverse current spectra provides evidence that the excess modes observed in the density of states of this and many other glasses, i.e. the Boson Peak, arises from the high-$`Q`$ limit of the quasi-transverse acoustic branch. This conclusion emerges from the comparison of the numerical results with the experimentally observed energy-shift and intensity variation of the Boson Peak with increasing density.

Insulating disordered solids, when compared to their crystalline counterparts, exhibit some common peculiarities in their low-temperature thermal properties and low-energy spectroscopic features phill ; phimaggen , and in particular (i) a larger specific heat at temperatures up to $``$1 K, ascribed to tunneling processes anderson ; phillips , (ii) a much smaller thermal conductivity, which also shows a plateau in the temperature range $``$1$`÷`$10 K phill , (iii) a quasi-elastic light- and neutron-scattering intensity, and, most important, (iv) an excess of modes in the vibrational density of states, known as the Boson Peak (BP). An unambiguous understanding of the origin of these extra modes, and of their possible relation with other reported anomalies, is still missing, in spite of the extensive research effort primed by the pioneering work of Buchenau et al. buc86  and continued by many authors ben96 ; for96 ; dellanna ; duval ; schirm ; tarask97 ; fontana99 ; pilla ; tarask01 . No definite conclusion can yet be drawn as to the nature of the vibrational eigenvectors of the modes responsible for the BP: according to different authors, in fact, they are either spatially localized buc86 ; for96 , spatially delocalized and propagating ben96 ; dellanna ; pilla , spatially delocalized but diffusive in character feldman . Even more important, still we do not completely understand why should disorder accumulate vibrational eigenvalues in the same energy region, in such a broad variety of chemically and physically different materials. Some recent theoretical works have been devoted to this subject, and among them we recall (i) the work of Elliott and coworkers, who assigned the BP in glasses to the lowest-energy van Hove singularity of the corresponding crystal tarask01 ; elliotpm02 ; (ii) the work of Grigera et al., who interpreted the BP as the precursor of the dynamical instability expected in a disordered structure as function of density Grigera ; and (iii) the work of Götze and Mayr and that of Schilling et al., who obtained a spectral feature recalling the BP within a mode-coupling-like description of the high frequency dynamics of a model glass Goetze ; Schilling .

In this Letter, we present a Molecular Dynamics (MD) study of the high-frequency dynamics in vitreous silica (v-SiO2) at different densities in the range $`\rho `$=2.2$`÷`$4.0 g/cm3. We demonstrate that the BP in vitreous silica arise from the high $`Q`$ portion of the transverse acoustic branch. Our results quantitatively support and better specify i) the Elliott and coworkers’ proposal tarask01 ; elliotpm02  that the BP is connected to the van Hove singularity of crystalline quartz; and ii) the Buchenau finding buc86  that the modes at the BP are local rotation of SiO2 tetrahedra.

Experimentally, in SiO2 at room pressure the BP appears as a broad peak in the plot of the Debye-normalized density of states $`g(E)/E^2`$, centered at $``$5 meV and having a width of $``$ 5 meV (Fig.1a). By increasing the sample density up to $`\rho `$=4.0 g/cm3, what is observed sugai66 ; inamura ; inabuche ; jund  is a high-energy shift of the peak and its concomitant intensity decrease (Fig. 1b,c). In Fig. 1, the lines are experimental data, while symbols are the results of the present simulations (see below), and we would like to stress the excellent agreement between the two. It is important to realize that the excess of modes in the density of state itself $`g(E)`$, corresponding to the peak at $``$5 meV in $`g(E)/E^2`$ of Fig. 1, is actually centered (at room pressure) around $``$ 15 meV, the shift being due to the $`1/E^2`$ factor and to the very slow decay of the high frequency tail of the peak in Fig. 1 . With respect to this point, it is important to underline that -in the high $`Q`$ limit and in the one-excitation (one-phonon) approximation- the dynamic structure factor $`S(Q,E)`$ becomes proportional to $`g(E)/E^2`$, and hence the longitudinal current spectrum $`C_L(Q,E)`$=$`E^2/Q^2S(Q,E)`$ becomes $``$ $`g(E)`$. Therefore, the modes giving rise to the BP must be searched, in the current spectra, around 15 meV. This observation matches with the high $`Q`$ panels of Fig. 2, where the computed currents (vide infra) show a broad bump at about this energy value.

The investigated system consists of 680 SiO2 units ($`N`$=2040 ions), enclosed in cubic boxes of different lengths (from $`L`$=3.1359 nm, corresponding to $`\rho =2.2g/cm^3`$ for the glass at room pressure, down to $`L`$=2.5693 nm corresponding to a density of 4.0 $`g/cm^3`$), with periodic boundary conditions. The ions interact through the BKS vanbe  two-body interaction potential; the long-range interaction was treated by the Ewald sum technique. As it has already been demonstrated, this system reproduce quantitatively the high frequency dynamics of vitreous silica dellanna . The glass configuration at room pressure was obtained by standard MD methods for lowering the temperature down to $`300`$ K starting from a well equilibrated liquid configuration at $`T=6000`$ K, followed by a conjugate gradient geometrical minimization on the potential-energy hypersurface for an accurate location of the minimum. Such minimum configuration was taken as the starting point to generate a series of compressed systems. At each compression step, the box size was scaled by $``$1.5%, then the system was made to relax, and the new minimum configuration was searched by the conjugate gradient method. This procedure was repeated until the final density of 4.0 $`g/cm^3`$, corresponding to a sample under an hydrostatic pressure of about 35 GPa zha , was reached. A complete study of the structural and dynamical changes occurring during the compression will be reported in a forthcoming paper. We focus here on the changes in the high frequency dynamics that take place as the density is increased.

The vibrational dynamics in the minimum configurations was computed in the harmonic approximation by diagonalizing the dynamical matrix, to obtain the eigenvalues ($`E_p`$) and eigenvectors ($`e_p(i)`$) of the $`p`$-th normal mode ($`p=1÷3N`$). From these quantities all vibrational characteristics can be derived. In particular, we have computed the density of states $`g(E)`$ and the longitudinal and transverse currents spectra ($`C(Q,E)`$) which, in the one-excitation approximation, are given by:

$$C_{_{\alpha \beta }}^\eta (Q,E)=[K_BT/\sqrt{M_\alpha M_\beta }]\mathrm{\Sigma }_pW_p^\eta (Q)\delta (EE_p)$$

where $`\eta =L,T`$, $`\alpha ,\beta `$ indicate Si and O, and $`W_p^\eta (Q)`$ is the spatial power spectrum of the (longitudinal or transverse) component of the eigenvectors:

$`W_p^L(Q)`$ $`=`$ $`|\mathrm{\Sigma }_i(\widehat{Q}\overline{e}_p(i))\mathrm{exp}(i\overline{Q}\overline{R}_i)|^2`$ (1)
$`W_p^T(Q)`$ $`=`$ $`|\mathrm{\Sigma }_i(\widehat{Q}\times \overline{e}_p(i))\mathrm{exp}(i\overline{Q}\overline{R}_i)|^2.`$
Here $`\widehat{Q}`$=$`\overline{Q}/|Q|`$.

The calculated $`g(E)/E^2`$ are reported in Fig. 1 together with the corresponding experimental curves. The two sets of data compare favorably to each other, indicating the suitability of the employed potential model to follow the density dependence of the high frequency dynamics. In Fig. 2 we report longitudinal and transverse current spectra at selected $`Q`$ values in the uncompressed sample ($`\rho `$=2.2 $`g/cm^3`$). For $`Q`$ values larger than about 8 nm-1, both $`C^^L(Q,E)`$ and $`C^^T(Q,E)`$ show two distinct excitation maxima, a feature that becomes more and more evident at increasing $`Q`$. The excitation at higher energy disperses with $`Q`$ and is observed at all $`Q`$ values in the longitudinal current spectra, while it shows up as a weak shoulder in the transverse current spectra only at $`Q>10`$ nm-1. In agreement with previous findings dellanna ; horbach ; taraskin98 , we assign this feature to the longitudinal sound-like branch. The behaviour of the low-energy excitation is in some sense complementary: it is always present in the transverse current spectra, while it appears in the longitudinal current spectra only at $`Q>8`$ nm-1. At small $`Q`$, the low-energy peak disperses with a sound velocity of $`3800`$ $`m/s`$ (appropriate for the transverse sound mode), and becomes almost non-dispersing at $`Q>8`$ nm-1 (Fig. 3). We will call this low-energy feature -which is the main feature in the transverse current spectra- the transverse acoustic mode.

The presence of the signature of transverse dynamics in the longitudinal current spectra, and vice-versa, is only apparently surprising. Indeed, the polarization character of the modes (which is better and better defined at increasing wavelength, i. e. when the vibration sees the medium as an elastic continuum) becomes ill defined at short wavelengths (mixing phenomenon). This is at the basis of the growth of peaks associated to the opposite polarization modes in the current spectra, and of the increased visibility of these peaks at increasing $`Q`$ values. In Fig. 3 we report the computed dispersion curves of transverse excitations (open symbols) for samples of three densities ($`\rho =2.2,\mathrm{\hspace{0.17em}2.7},\mathrm{\hspace{0.17em}4.0}`$ g/cm3), and of longitudinal excitations (crossed circles) of the uncompressed system nota2 . In the same figure (full symbols), are also reported -at the available Q values, the experimental excitation energies (maxima of the current spectra) of the uncompressed SiO2 glass, as measured by IXS and INS fontana ; righetti . The agreement between MD and experimental data clearly indicates that the peaks observed in the INS experiment at $`Q`$ larger than 8 nm-1 for96 , must be associated with the transverse dynamics, which appears in the measured (longitudinal) spectra due to the mixing.

In Fig. 3, we also observe a strong positive dispersion of the velocity of the (longitudinal) sound waves, which continue to propagate at energy well above the BP energy. This dispersion, observable also in MD simulations of Lennard-Jones glasses relax-harm  and in previous v-SiO2 calculations horbach , can be due to the interaction of the sound waves with a relaxation process relax-harm  or with other modes. The investigation of the origin of the positive dispersion of $`v_L`$, and is possible interplay with the BP modes, is beyond the purpose of the present work.

For our present purposes, the first important result emerging from Fig. 3 is that the ”transverse” branches (open symbols) in low- and intermediate-density samples, at large $`Q`$’s, flatten at high $`Q`$’s to an energy value which increases with increasing density ($``$15 meV at room condition and $`Q`$=15 nm-1). The density of vibrational states ($`g(E)`$), associated with the branches which flatten, will have an excess of modes with respect to the Debye behaviour at these energies, reminiscent of the van Hove singularity of the corresponding crystal tarask01 . As mentioned, the corresponding peak in $`g(E)/E^2`$, will be red-shifted at lower energies ($``$5 meV at room condition). Therefore, we conclude that the Boson peak originates from the modes associated to the flat portion of the acoustic transverse dispersion curve.

This result is in agreement with the recent theoretical work of S. Taraskin et al., who associate the BP to the glassy counterpart of the lowest-energy van Hove singularity of the corresponding crystalline structure tarask01 . In this respect, it is worth noting that the transverse acoustic branch at $`Q`$ larger than $``$8 nm-1 is the glassy counterpart of a transverse optic phonon branch of $`\alpha `$-quartz (almost flat at $``$4 THz, i. e. $``$16 meV). This branch, in the extended Brillouin zone scheme which is more appropriate for disordered materials, is the prosecution of the transverse acoustic branch dorner . According to Boysen at al. boysen , this branch is (close to the M-point) strongly temperature-dependent, and its softening is responsible for the $`\alpha `$-$`\beta `$ transition in quartz. More importantly, the atomic displacements induced by the lattice modes of this branch in quartz, as determined in dorner , are very similar to the frustrated localized rotation of SiO4 tetrahedra, that in vitreous silica -according to Buchenau et al. buc86 \- are the modes contributing to the BP.

The assignment of the BP to the flattening of the transverse acoustic branch is strengthened by the behaviour of the high-energy transverse dynamics upon densification. Spectra similar to those of Fig. 2, but for the sample at the highest studied density ($`\rho `$=4.0 g/cm3), are reported in Fig. 4, while the corresponding dispersion curves are reported in Fig. 3. One can observe that (i) the two current spectra are now ”pure”, and no evidence of wrong-polarized modes is present; (ii) the positive dispersion of the L-branch is absent; and, more importantly, (iii) the T-branch no longer shows a flattening. This latter observation, together with the experimental and numerical finding that the BP intensity strongly decreases on increasing density inamura ; inabuche ; jund  (as evident in Fig. 1), gives a decisive support to the finding that the BP is produced by the flattening of the quasi-TA branch.

In conclusion, by comparing the spectra of the longitudinal and transverse current, and by studying their density dependence, we have shown that the BP in v-SiO2 is to be ascribed to the quasi-transverse acoustic modes, whose dispersion relation becomes $`Q`$-independent at high $`Q`$, and consequently, gives rise to an excess of modes with respect to the Debye behaviour. The whole picture presented here reconciles different previous studies on the origin of the BP: i) the high $`Q`$ part of the transverse acoustic branch in vitreous silica -which gives rise to the BP- is the counterpart of a low lying transverse optic branch in $`\alpha `$-quartz dorner ; ii) according to Boysen at al. boysen  the softening of this branch at the M-point give rise to the $`\alpha `$-to-$`\beta `$ transition; iii) in agreement with i) and ii), according to Taraskin et al. tarask01 ; elliotpm02 , the BP arises from the softening of the lowest-energy van Hove singularity of the corresponding crystals; iv) according to Dorner and coworkers the eigenvectors of this branch in quartz correspond to rotation of SiO4 tetrahedra; v) finally, in agreement with i) and iv), according to Buchenau et al. buc86  the modes contributing to the BP are localized rotation of SiO4 tetrahedra.

This work was supported by INFM Iniziativa di Calcolo Parallelo, and by MURST Progetto di Ricerca di Interesse Nazionale. One of us (GR) greatly acknowledge illuminating discussions with U. Buchenau and B. Dorner.

FIGURE CAPTIONS

FIG. 1: Normalized density of states $`g(E)/E^2`$ of vitreous silica at three different densities $`\rho `$: (a) 2.2 g/cm3 (full line and circles); (b) 2.7 g/cm3 (dashed line and squares); (c) 4.0 g/cm3 (triangles). Lines are experimental data from Ref. inabuche , symbols are the results of the present simulation.

FIG. 2: Selected examples of longitudinal (full lines) and transverse (dashed lines) current spectra at the indicated $`Q`$-values in nm-1, for the sample at $`\rho `$=2.2 g/cm3.

FIG. 3: Dispersion relation of the main peaks appearing in the current spectra. Crossed circles are the maxima of the longitudinal current at $`\rho `$=2.2 g/cm3 (the line is a guide for the eyes). The open symbols refer to the maxima of the transverse current at $`\rho `$=4.0 g/cm3 (down-triangles), $`\rho `$=2.7 g/cm3 (squares), $`\rho `$=2.2 g/cm3 (up-triangles). Full symbols are INS (squares) and IXS (circles and diamonds) experimental data for $`\rho `$=2.2 g/cm3 taken from Ref.s ben96 ; fontana ; righetti .

FIG. 4: Selected examples of longitudinal (full lines) and transverse (dashed lines) current spectra at the indicated $`Q`$-values (in nm-1) for the highest density sample ($`\rho `$=4.0 g/cm3).

Figure 1. Pilla et al.

Figure 2. Pilla et al.

Figure 3. Pilla et al.

Figure 4. Pilla et al.