text:
# Spin-orbit coupling and intrinsic spin mixing in quantum dots

## Abstract

Spin-orbit coupling effects are studied in quantum dots in InSb, a narrow-gap material. Competition between different Rashba and Dresselhaus terms is shown to produce wholesale changes in the spectrum. The large (and negative) $`g`$-factor and the Rashba field produce states where spin is no longer a good quantum number and intrinsic flips occur at moderate magnetic fields. For dots with two electrons, a singlet-triplet mixing occurs in the ground state, with observable signatures in intraband FIR absorption, and possible importance in quantum computation.

spin-orbit coupling, Rashba effect, quantum dots

The creation and manipulation of spin populations in semiconductors has received a great deal of attention in recent years. Conceptual developments that have motivated these efforts include prominently the Datta-Das proposal for a spin field-effect transistor,1  based on the Rashba spin-orbit coupling of electrons in a 2DEG,2  and the possibility of building quantum computation devices using quantum dots (QDs).3  It is then important for full control of spin-flip mechanisms in nanostructures that all spin-orbit (SO) effects be understood.

There are two main SO contributions in zincblende materials like $`A_3B_5`$: in addition to the structure inversion asymmetry (SIA) caused by the 2D confinement (the Rashba effect), there is also a SO term caused by the bulk inversion asymmetry (BIA) in those structures (the Dresselhaus term).5  Notice that additional lateral confinement defining a QD introduces another SIA term with important consequences, as we will see in detail. Although the relative importance of these two effects depends on the materials and structure design (via interfacial fields), only recently have authors begun to consider the behavior of spins under the influence of all effects. For example, a modification of the Datta-Das design was recently suggested to allow for a diffusive version of the spin FET, 4  and that proposal relies on the additional influence of the Dresselhaus SO coupling in the system.

Work in wide-gap materials (mainly GaAs), 6  uses a unitary transformation on the Hamiltonian of the system, 7  after which one gets an effective diagonal SO term which incorporates the Rashba effect in a perturbative fashion. That approach is valid since the SO coupling is small in GaAs. However, the approximation is not valid for all of the $`A_3B_5`$ structures, as it is the case for InSb, for example, where both SIA and BIA effects are anticipated to be much larger. Cardona  In this material, one needs to deal with the full Hamiltonian.

There are just a few works discussing SO effects in narrow-gap nanostructures. Among them, \[9, \] uses $`𝐤𝐩`$ theory in InSb QDs in order to include SIA SO terms from both the Rashba field and the lateral confinement which defines the QD. This last SIA term is considered in \[11, \], and since it is diagonal in the Fock-Darwin (FD) basis no level mixing is found nor expected. In contrast, level mixing events are clearly identified in \[9, \]. Experiments in InSb QDs have explored the FIR response in lithographically defined dots, 10a  and PL features of self-assembled dots.10b

The goal of this work is to show how important different types of SO couplings are in the spectra of parabolic QDs built in narrow-gap materials such as InSb. We consider the Rashba-SIA diagonal and SIA non-diagonal, as well as the Dresselhaus-BIA terms in the Hamiltonian, and proceed with its full diagonalization, in order to study features of the spectrum as function of magnetic field, dot size, $`g`$-factor, and electron-electron interaction. We draw attention to the appearance of strong level anticrossings (mixing) for moderate magnetic fields in typical QDs, and how this phenomenon (and ‘critical’ field where it occurs) is modified by the BIA terms not considered before.9  As the level mixing involves states with different spin, this induces strong intrinsic spin flips in the system, regardless of the strength of the SO coupling, providing an important channel for spin decoherence in these systems. Moreover, measurement of FIR absorption would yield direct access to the coupling constants; i.e., the dispersion of FIR absorption peaks and appearance of additional/split-off features are a direct consequence of the level mixing introduced by SO.

Model. Assuming a heterojunction or quantum well confinement $`V(z)`$ such that only the lowest $`z`$-subband is occupied, the Hamiltonian in the absence of SO interactions for a QD further defined by a lateral parabolic confinement is given by $`H_0=\frac{\mathrm{}^2}{2m}𝐤^2+V(\rho )+\frac{1}{2}g\mu _B𝐁\sigma `$, where $`𝐤=i+e𝐀/(\mathrm{}c)`$, and the in-plane vector potential $`𝐀=B\rho (\mathrm{sin}\theta ,\mathrm{cos}\theta ,0)/2`$ describes a perpendicular magnetic field $`𝐁=B𝐳`$; $`m`$ is the effective mass in the conduction band,remark\_non-parab  $`g`$ is the bulk $`g`$-factor, $`\mu _B`$ is Bohr’s magneton, $`V(\rho )=\frac{1}{2}m\omega _0^2\rho ^2`$ is the lateral confinement with frequency $`\omega _0`$, and in the Zeeman term $`\sigma _{X,Y,Z}`$ are the Pauli matrices. The analytical solution of $`H_0`$ yields the FD spectrum with energies $`E_{nl\sigma }=(2n+|l|+1)\mathrm{}\mathrm{\Omega }+l\mathrm{}\omega _c/2+g\mu _BB\sigma /2`$, with effective frequency $`\mathrm{\Omega }=\sqrt{\omega _0^2+\omega _C^2/4}`$, and cyclotron frequency $`\omega _C=eB/(mc)`$; states are given in terms of associated Laguerre polynomials. FD  The confinement, magnetic and effective lengths are, respectively, $`l_0=\sqrt{\mathrm{}/(m\omega _0)}`$, $`l_B=\sqrt{\mathrm{}/(m\omega _C)}`$ and $`\lambda =\sqrt{\mathrm{}/(m\mathrm{\Omega })}`$.

The SIA SO term with coupling parameter $`\alpha `$ is $`H_{SIA}=\alpha \sigma \left(V\times 𝐤\right)`$, where the total confinement potential is $`V(𝐫)=V(\rho )+V(z)`$. One can then write $`H_{SIA}=H_R+H_{SIA}^D`$, where the diagonal contribution coming from the lateral confinement in cylindrical coordinates is $`H_{SIA}^D=\alpha \frac{\mathrm{}\omega _0}{l_0^2}\sigma _Z\left(L_Z+\frac{\lambda ^2}{l_B^2}\frac{x^2}{2}\right)`$, with the adimensional radial coordinate $`x=\rho /\lambda `$, and $`L_Z=i/\theta `$.

The Rashba term coming from the perpendicular confinement field $`dV/dz`$ is

$$H_R=\alpha \frac{dV}{\lambda dz}\left[\sigma _+L_{}A_{}+\sigma _{}L_+A_+\right]\text{,}$$
(1)
where $`L_\pm =\mathrm{exp}(\pm i\theta )`$, $`\sigma _\pm =(\sigma _X\pm i\sigma _Y)/2`$, and operators $`A_\pm =/x+L_Z/x+x\lambda ^2/(2l_B^2)`$.

In zincblende structures one should also consider the BIA SO bulk Hamiltonian. 5  After averaging in the $`z`$-direction, due to quantization, one gets $`H_{BIA}=\gamma \left[\sigma _xk_xk_y^2\sigma _yk_yk_x^2\right]+\gamma k_z^2\left[\sigma _yk_y\sigma _xk_x\right]+\gamma \sigma _zk_z\left(k_x^2k_y^2\right)`$, where $`\gamma `$ is the coupling parameter, the resulting first (second) term is cubic (linear) in the in-plane momentum, and the last term is zero because $`k_z=0`$; also, $`k_z^2(\pi /z_0)^2`$, where $`z_0`$ is the $`z`$-direction confinement length. One may write the BIA SO term as $`H_{BIA}=H_D^L+H_D^C`$, where the linear Dresselhaus contribution is given by

$$H_D^L=i\frac{\gamma k_z^2}{\lambda }\left[\sigma _+L_+A_+\sigma _{}L_{}A_{}\right]\text{,}$$
(2)
while the cubic contribution $`H_D^C`$ can be expressed in terms of $`\sigma _{}L_\pm ^3`$ and $`\sigma _\pm L_\pm `$, and different powers in $`x`$, $`/x`$, and $`L_z`$. Destefani-thesis  Notice that under a finite magnetic field, the matrix elements with $`\sigma _\pm L_\pm `$ in $`H_D^C`$ are not hermitian, and one needs to symmetrize them; 7  if the field is zero, this problem does not occur. 13A

For the electron-electron interaction $`H_{ee}`$, an expansion in Bessel functions for $`|𝐫_1𝐫_2|^1`$ is employed. Destefani-thesis  The basis states are properly antisymmetrized, describing the unperturbed spin eigenstates.

The general form of the various SO terms in the Hamiltonian exhibit already interesting characteristics. For example, the magnetic field plays a role via its linear dependence in $`H_{SIA}^D`$, $`H_R`$, and $`H_D^L`$, or its $`B`$ to $`B^3`$ dependence in $`H_D^C`$. Destefani-thesis  Most interestingly, this form of the Hamiltonian yields selection rules explicitly, dictating which levels will be influenced by the SO effects. For example, at zero field the diagonal SIA term splits the levels according to the total angular momentum $`j`$. The Rashba term induces a set of anticrossings in the FD spectrum whenever $`\mathrm{\Delta }l=\pm 1=\mathrm{\Delta }\sigma `$ at finite field (due to the $`\sigma _\pm L_{}`$ terms; mostly negative $`l`$’s are affected since their magnetic dispersions allow for crossings); the lowest anticrossing is between $`\{n,l,\sigma \}=\{0,0,\}`$ and $`\{0,1,+\}`$. The cubic BIA terms (with $`\sigma _{}L_\pm ^3`$) induce a set of anticrossings, which obey $`\mathrm{\Delta }l=3`$ and $`\mathrm{\Delta }\sigma =\pm 1`$; the first one at low $`B`$-field involves the states $`\{0,1,\}`$ and $`\{0,2,+\}`$. Terms with $`\sigma _\pm L_\pm `$ in $`H_D^L`$ and $`H_D^C`$ do not induce anticrossings, but rather split and shift the spectrum due to matrix elements with $`\mathrm{\Delta }l=\pm 1=\mathrm{\Delta }\sigma `$. Notice that the matrix elements between states with different $`n`$’s are in general non-zero, so that the full diagonalization involves mixings with various $`n`$-values.

Results. The sequence of FD states of $`H_0`$ starts at zero $`B`$-field with $`\{n,l,\sigma \}=\{0,0,\pm \}`$, followed by the degenerate set of $`\{0,1,\pm \}`$ and $`\{0,1,\pm \}`$.FD  Spin and orbital degeneracies are broken by $`B`$ and the states with negative $`l`$ and positive $`\sigma `$ acquire lower energies because of the negative $`g`$-factor. The lowest energy level crossing is between states $`\{0,0,\}`$ and $`\{0,1,+\}`$, and the field where it occurs in the FD spectrum is

$$B_C^0=\frac{\stackrel{~}{m}}{\mu _B}\frac{\mathrm{}\omega _0}{\sqrt{\stackrel{~}{m}|g|(\stackrel{~}{m}|g|+2)}}\text{,}$$
(3)
where $`\stackrel{~}{m}=m/m_0`$. The moderate value of $`B_C^0`$ is a direct consequence of the large $`|g|`$ in InSb.parameters  For GaAs ($`|g|=0.44`$, $`\stackrel{~}{m}=0.067`$), for example, this level crossing appears only at $`B_C^{GaAs}9.4`$T for a much smaller confinement, $`\mathrm{}\omega _0=2`$meV, and in the region where Landau levels are well defined. Weaker confinement (smaller $`\omega _0`$) shifts this crossing to lower fields. Notice that for $`g<0`$, $`H_R`$ mixes these states (and $`H_D^L`$ shifts the crossing to higher fields). For $`g>0`$ it is $`H_D^L`$ that would produce relatively stronger level anticrossings (and $`H_R`$ would only shift the spectrum weakly), and it would then be absent in non-zincblende materials like silicon.

The energy spectrum for InSb QDs with typical characteristics, parameters  and for the full Hamiltonian is presented in Fig. 1A vs. $`B`$ field. The spectrum is obtained by direct diagonalization using a FD basis with $`n4`$ (or ten energy ‘shells’), i.e. 110 basis states. We have studied the progressive changes to the FD levels when including different SO terms in $`H`$.

The diagonal $`H_{SIA}^D`$ term shifts energies but does not change appreciably the position of the first crossing (3), shown by the dotted lines in Fig. 1A at $`2.6`$T. The energy shifts induce two new crossings at low fields (inset d), since the SO orders states according to their total angular momentum $`j=l+s`$; the highest (lowest) state at zero field has $`j=3/2`$ $`(1/2)`$ in the second shell. At about 0.2T one recovers the ‘normal’ sequence of states: $`\{0,1,+\}`$, $`\{0,1,\}`$, $`\{0,1,+\}`$, $`\{0,1,\}`$. This competition between SO and magnetic field is similar to the Zeeman and Paschen-Back regimes in atoms. Paschen-Back-ref  We should note that this level ordering is observed in \[11, \].

The non-diagonal Rashba contribution $`H_R`$ introduces strong state mixing for any value of the $`\alpha `$ parameter whenever FD levels with $`\mathrm{\Delta }l=\mathrm{\Delta }\sigma =\pm 1`$ cross. This mixing converts the crossings at $`B_C^0`$ to clear anticrossings. Higher levels which satisfy these selection rules also anticross at nearly the same field. The field-width and energy-amplitude (or level splitting) of the mixing is dictated by the value of $`\alpha `$, while the value of $`B_C`$ where the anticrossing occurs is nearly unaffected by $`\alpha `$.

The cubic Dresselhaus contribution $`H_D^C`$ induces anticrossings (via $`\sigma _{}L_\pm ^3`$) and zero-field splittings ($`\sigma _\pm L_\pm `$) in the FD spectrum. The splittings are much smaller than those induced by the Rashba term and practically unnoticeable in the spectrum, reflecting the smallness of the $`E_D^C`$ for these parameters. parameters  For the linear BIA contribution, however, the $`\sigma _\pm L_\pm `$ terms in $`H_D^L`$ have a much bigger impact on the zero-field splittings, which can in principle be ‘tuned’ by changing the effective $`z`$-size, $`z_0`$. $`H_D^L`$ alone induces such a strong mixing at low fields that one cannot identify the two Zeeman and Paschen-Back regimes.

Notice in the full spectrum of $`H`$ (Fig. 1A) that the first group of anticrossings (for $`n=0`$ levels) induced by $`H_R`$ is shifted to higher field due mostly to $`H_D^L`$, so that $`B_C^0B_C3.3`$T (box b and lower arrow). The set of anticrossings at $`5.5`$T is also due to $`H_R`$ and arises from the $`n=1`$ level manifold (box c and upper arrow). At low field, only a single crossing in the second shell at $`0.3`$T is present and dominated by $`H_D^L`$ (box a; compare with inset d and notice second crossing in box e). The sequence of the first excited levels at zero field is $`j=3/2`$ $`(1/2)`$ for higher (lower) energy, while at higher energies both SIA and BIA terms cooperate to produce anticrossings (not visible at the resolution in Fig. 1A).

Figure 2 illustrates the importance of the level anticrossings on the spin, as the expectation value of $`S_z`$ for each state is plotted vs. $`B`$. Figs. 2A and B include all states with $`E80`$meV (for full SO and only SIA terms, respectively), while figs. 2C and D focus only on the lowest seven levels. Although a large majority of states have $`S_z`$ close to $`\pm 1/2`$, as one expects for pure states, there are significant deviations. The various SO terms mix levels close to accidental degeneracy points in the FD spectrum and produce the large deviations seen in the figure. 2C shows how $`H_R`$ produces an intrinsic (i.e., no phonon-assisted) total collapse of the spin number for the low energy states in the QD. Although the ground state is nearly pure ($`S_z1/2`$, and more so at higher $`B`$), the first few excited states totally mix at $`B_C3.3`$T. 2D shows how a stronger Rashba field ($`dV/dz=2\times 10^3`$eV/Å) greatly widens the mixing region and lowers $`B_C2.8`$T.

One can further appreciate the intricate balance of SO terms under a magnetic field. We analyze how various quantities are affected by changes in the the lateral and vertical sizes, $`l_0`$ and $`z_0`$, or the Rashba field $`dV/dz`$, as shown on Fig. 1B. The zero-field splitting (dotted lines) is dominated by the linear BIA contribution for any value of $`l_0`$ here. Increasing $`z_0`$ strongly reduces the splittings because the Dresselhaus contribution weakens; the reduction is even more drastic if one increases $`dV/dz`$, which makes the $`H_R`$ contribution bigger and can then cancel or suppress better the splitting produced by $`H_{SIA}^D`$. Some authors have considered the possibility of tuning such SO terms to produce total cancellation of the zero-field splitting, although considering only $`H_R`$ and $`H_D^L`$.4  However, one also has to take into account $`H_{SIA}^D`$ and $`H_D^C`$ contributions, which may be important (the zero-field cancellation occurs at values of $`z_0`$ or Rashba field about ten percent smaller than with only the former terms). One should notice, in any event, that this change in parameters only eliminates the zero-field splitting but not the anticrossing at finite field, and measurement of both quantities on the same sample could yield information on the relative strength of the $`\alpha `$ and $`\gamma `$ parameters.

The anticrossing field $`B_C`$ (solid lines/symbols) decreases with QD size, roughly according to (3), $`B_C^0\omega _01/\sqrt{l_0}`$. A finite $`\alpha `$ slightly increases $`B_C`$, but the BIA contribution considerably upshifts it, as mentioned above. Increasing $`z_0`$ or $`dV/dz`$ decreases $`B_C`$. At $`l_0=320`$Å ($`\mathrm{}\omega _0=5`$meV), $`B_C=1.6`$T, while it shifts to 1.15T if $`dV/dz`$ is four times larger or to 0.85T if $`z_0`$ is doubled, both cases decreasing the BIA contribution. These values are comparable to those in \[9, \] without including BIA terms (adjusting for differences in system parameter values). Anticrossings at such low fields may be interesting for applications due to easier access.

The energy splitting at $`B_C`$ (inset in Fig. 1B), has main contribution from the Rashba term for any dot size considered, but the BIA reduces the splitting substantially. If $`z_0`$ is changed from $`40`$ to 80Å the splitting is enhanced slightly, but larger $`z_0`$ produces no significant changes. However, the splitting is drastically enhanced if one increases the Rashba field. Here, the splitting goes from 1 to 4.2meV if the interface field is increased fourfold.

Figure 3 illustrates the corresponding level structure for two electrons in the QD (full Hamiltonian $`H+H_{ee}`$; dashed lines show $`H_0+H_{ee}`$, the non-SO case). The repulsive interaction shifts the ground state upwards by $`5`$meV, and the exchange shifts the triplet down by 2meV. Most interestingly, the SO interaction introduces a strong mixing of the singlet and triplet transition at $`B2.7`$T. The fact that the mixing occurs at relatively low field makes that a possibly useful transition for the implementation of quantum computing devices. Moreover, the splitting will also be apparent in the FIR response of QDs, allowing the determination of the various SO coupling strengths.

We have shown that inclusion of all SO terms is essential in order to obtain a complete picture of the level structure in narrow-gap QDs. The combination of strong SO couplings and large (and negative) $`g`$ factor introduces strong intrinsic mixing of the low excitations for the single-particle spectrum. Consequently, the two-particle spectrum exhibits strong singlet-triplet coupling at moderate fields, with significant experimental consequences. Observation of FIR mode magnetic dispersion would allow the direct determination of coupling constants.

We acknowledge support from FAPESP-Brazil, US DOE grant no. DE-FG02-91ER45334, and the CMSS Program at OU.


format: text

text:
# A “poor man’s approach” to modelling of micro-structured optical fibres

## 1 Introduction

Since the first results on photonic crystal fibres  (PCF) many exciting phenomena have been reported (for a recent review see e.g. Ref.  and references therein). From the very start, a great emphasis has been on the modelling of the optical properties and frameworks of high complexity have been developed. In this work, we develop a “poor man’s approach” which allows for easy implementation and calculation of the propagation constant, $`\beta (\omega )`$, in arbitrarily microstructured fibres. Most approaches start from the fundamental vectorial wave equation for an isotropic dielectric medium

$$\mathbf{}\times \frac{1}{\epsilon (𝒓)}\mathbf{}\times 𝑯(𝒓)=k^2𝑯(𝒓)$$
(1)
where $`k=\omega /c`$ and $`\epsilon (𝒓)`$ is the dielectric function, which may be frequency dependent. For a fibre geometry with $`\epsilon (𝒓)=\epsilon (x,y)`$, i.e., translational invariance along the $`z`$-direction, we look for solutions of the form $`𝑯(𝒓)=𝒉(x,y)e^{i\beta z}`$. Substituting this ansatz into the wave equation results in an eigenvalue problem, which determines the dispersion $`\omega (\beta )`$. In the literature, it is often emphasized that in general a fully-vectorial approach is required for quantitative modelling of micro-structured fibres. Several fully-vectorial approaches have been reported including plane-wave methods , multi-pole methods , localised-function methods , and finite-element methods . The complicated implementation is common to all these methods. In this work, we present a “poor man’s approach” for calculating the propagation constant, $`\beta (\omega )`$, in arbitrary dielectric structures. Despite its obvious shortcomings, we believe it is far more useful for more qualitative studies as well as for teaching of the physics of micro-structured optical fibres at the introductory level.

## 2 Theory

### 2.1 The scalar Helmholtz wave-equation

We start from the scalar Helmholtz wave-equation

$$\left(\frac{^2}{x^2}+\frac{^2}{y^2}+\epsilon (x,y)k^2\right)\mathrm{\Psi }(x,y)=\beta ^2\mathrm{\Psi }(x,y)$$
(2)
with $`\mathrm{\Psi }`$ being the scalar field. In this approximation, we have neglected a logarithmic derivative of the dielectric function at the interfaces between e.g. SiO2 and air. We also note that polarisation effects and/or degeneracies may only be fully revealed by a fully-vectorial approach. The strength of the scalar approach is that it is straight forward to implement and thus serves as an excellent starting point for researchers and students entering the field of micro-structured optical fibres. Furthermore, the Helmholtz equation (2) allows for an easy incorporation of effects like material dispersion, $`\epsilon (\omega )`$.

### 2.2 The finite-difference approximation

Equation (2) constitutes an eigenvalue problem from which $`\beta (\omega )`$ can be calculated. We take the probably most simple approach based on a finite difference approximation of the Laplacian. For a quadratic grid with $`j`$ labeling the grid points with spacing $`a`$ we e.g. get the symmetrized representation of the Laplacian

$$\frac{^2}{x^2}f[x=ja]\frac{1}{a^2}\left(f[(j+1)a]+f[(j1)a]2f[ja]\right)$$
(3)
corresponding to nearest-neighbour coupling of the grid points. We can thus rewrite Eq. (2) as a matrix eigenvalue problem

$$𝚯𝚿=\beta ^2𝚿$$
(4)
with

$$𝚯_{ji}=\{\begin{array}{cc}4K^2+\epsilon _jk^2& j=i\\ K^2& j,i\mathrm{nearest}\mathrm{neighbours}\\ 0& \mathrm{otherwise}\end{array}$$
(5)
where $`K=1/a`$. For the dielectric function we have $`\epsilon _j=\epsilon (x_j,y_j)`$ where $`(x_j,y_j)`$ is the coordinates of the $`j`$th lattice point. The numerical task thus consists of finding eigenvalues of the matrix $`𝚯`$, which is easily done using standard numerical libraries. The matrix is highly sparse, symmetric and when the dielectric function is real, it is even Hermitian. The numerical accuracy of the approximation is of course increased by decreasing the lattice constant.

### 2.3 The homogeneous case

In order to estimate the required size of $`a`$, we first consider the homogeneous case where $`\epsilon _j=\epsilon `$. This problem is well-known from solid state theory; it can be diagonalized by a plane wave ansatz, which results in the usual cosine-band result

$$\omega ^2=\frac{c^2}{\epsilon }\left[\beta ^2+2K^2(2\mathrm{cos}\kappa _xa\mathrm{cos}\kappa _ya)\right].$$
(6)
This result also has the correct asymptotic behaviour of the homogeneous-space solution

$$\omega ^2\frac{c^2}{\epsilon }\left[\beta ^2+\kappa _x^2+\kappa _y^2\right]+𝒪(a^2)$$
(7)
and by choosing $`a`$ sufficiently small, the numerical discretisation is a good approximation to the solution of the exact problem. Because of the discretisation procedure, Eq. (6) has a finite band-width of $`\frac{c^2}{\epsilon }\times 8K^2`$, i.e.,

$`\mathrm{max}\{2K^2(2\mathrm{cos}\kappa _xa\mathrm{cos}\kappa _ya)\}`$ (8)
$`\mathrm{min}\{2K^2(2\mathrm{cos}\kappa _xa\mathrm{cos}\kappa _ya)\}=8K^2.`$
This means that only frequencies satisfying

$$\frac{c^2}{\epsilon }\beta ^2\omega ^2\frac{c^2}{\epsilon }\left[\beta ^2+8K^2\right]$$
(9)
can be accounted for by the discretisation procedure. In the low-frequency regime, the relative error of the finite-difference approximation becomes small, and typically we should be limiting ourselves to e.g.,

$$\omega ^2<\frac{c^2}{\epsilon }\left[\beta ^2+K^2\right]a<\left(\epsilon \omega ^2/c^2\beta ^2\right)^{1/2}$$
(10)
where the relative error of the finite-difference approximation is less than $`9\%`$, see Fig. 1. For higher frequencies, the finite-difference procedure becomes an inaccurate approximation to the exact problem, because of artificial band-structure effects.

### 2.4 The boundary problem

In principle $`𝚯`$ is infinite and for the implementation we obviously need to truncate the matrix. This truncation may affect the accuracy of the calculation. However, we are also faced with the problem of deciding how the finite-difference representation of the differential operators (in our case the Laplacian) are represented on the boundary of the calculation domain. This problem arises because calculation of finite differences on the boundary requires the use of points that fall outside the calculation domain. Thus, we have to determine a proper way of how these non-existing points should be treated.

In the definition of the $`𝚯`$-matrix, we have simply chosen to neglect the points that fall outside the calculation domain. This is done in order to limit the complexity of the $`𝚯`$-matrix, and thereby keep the formulation of the problem as simple as possible. The consequence of this simplification is that we impose the condition that the field has to be zero on the boundaries of the calculation domain. Naturally, this assumption will have an effect on the accuracy of the calculation, but the better the field is confined within the calculation domain, the better the truncated problem resembles the correct solution, since a zero amplitude on the boundary becomes a reasonable approximation in this case.

## 3 Modelling of Photonic Crystal Fibres

### 3.1 Numerical Implementation

Once the theory of the finite difference approximation has been established, the task of finding solutions to the scalar wave equation may be viewed as two subproblems. First, the $`𝚯`$-matrix has to be created from a given dielectric structure, and secondly the eigenvalues, $`\beta ^2`$, and the corresponding eigenvectors, $`𝚿`$, have to be found. We have chosen to make our implementation in *Matlab*, because it provides effective tools for solving both these problems.

To give a more precise description of our implementation, we consider an example where an arbitrary dielectric structure has been discretisized in a 100 $`\times `$ 100 grid. In this case $`𝚯`$ becomes a matrix with 10000$`\times `$10000 elements, which indicates that the finite difference approach is very demanding in terms of memory consumption. However, as most of the elements of $`𝚯`$ are zero, it is advantageous to store $`𝚯`$ as a sparse matrix, which is easily done with the *Sparse-type* in *Matlab*. For an n $`\times `$ n dielectric structure, we need to store $`n^4`$ elements in the full representation, while only $`5n^2`$ elements are required in the sparse representation. Obviously, this gives rise to a dramatic decrease in the memory consumption, and thereby a corresponding increase in the size of the dielectric structures that may be examined.

The second problem we are faced with, concerns the search for eigenvalues, $`\beta ^2`$, and their corresponding eigenvectors $`𝚿`$. In the sparse representation, a complete diagonalization of the $`𝚯`$-matrix may be performed, but this straightforward method has the disadvantage of being numerically heavy (and thus time consuming) since it calculates all $`n^2`$ eigenvalues. Typically, we are only interested in solving for a few eigenvalues, e.g the largest values of $`\beta ^2`$, and this is a common feature of several advanced eigensolver libraries. In our implementation we use the *EIGS* command in *Matlab* which is based on the ARPACK-library . As a further advantage, the *EIGS* command has the option of finding eigenvalues around a specified value, which may be particularly useful once a region with guided modes has been found.

### 3.2 Index-guiding Fibres

As a first example, the finite-difference method is applied to an index-guiding micro-structured fibre. Figure 2a shows the square dielectric structure used in the calculation, which we have chosen to discretisize in a 128 $`\times `$ 128 grid. The dielectric structure has a width of $`3\sqrt{3}\mathrm{\Lambda }`$, where $`\mathrm{\Lambda }`$ is the hole spacing, and it consists of air holes with a diameter of 0.4$`\mathrm{\Lambda }`$ placed in a silica background (n=1.45). A single air hole has been omitted in the center of the structure to create a waveguide core. In the case of index-guiding fibres, the fundamental mode corresponds to the eigenvector with the largest eigenvalue. Figure 2a shows the field distribution of the fundamental mode calculated at a normalized wavelength of $`\lambda /\mathrm{\Lambda }=0.15`$, and it is seen to be well confined to the core region.

In figure 3 we have mapped out the effective mode index of the fundamental mode for several values of the normalized wavelength. For comparison, we have included a finite-difference calculation, where the width of the calculation domain has been increased to $`6\sqrt{3}\mathrm{\Lambda }`$. Finally, we have included fully-vectorial results for a fully-periodic hole-structure, with a repeated defect, obtained in a plane-wave basis using periodic boundary conditions . The dielectric structure used in this calculation consist of 6 simple cells, and is thus similar to the structure in figure 2a.

By comparing the finite-difference calculations for the small and the large calculation domain, we are able to see the effect of the truncation. For the small calculation domain, we find that the calculated value of the effective mode-index decreases rapidly for long normalized wavelengths. This is because the field distribution extends to the edge of the calculation domain, and thus the assumption of a zero amplitude on the edge is no longer valid. Consequently, the field will penetrate into the air holes, and thereby cause a lowering of the effective mode index. By increasing the width of the calculation domain, we are able to shift the onset of this effect towards larger values of the normalized wavelength.

In the comparison between our scalar finite-difference approach and a full-vectorial method, we find that the scalar approach gives qualitative correct results and accounts well for the overall wavelength dependence. However, especially for the long wavelengths the simple approach becomes inaccurate because of the scalar approximation, and the scalar approach is seen to overestimate the value of the effective mode-index. For $`\lambda \mathrm{\Lambda }`$ the strong proximity of the air-holes require a correct treatment of the air-silica boundary conditions which can only be quantitatively accounted for by a fully-vectorial approach.

### 3.3 Photonic Bandgap-guiding Fibres

In the literature, it is often argued that accurate modelling of photonic bandgap materials requires the use of a full-vectorial method. This is true for all photonic crystals of practical interest, because it is required that they have a large index-contrast in order to obtain large bandgaps. However, this may lead to the incorrect conclusion that photonic bandgaps and defect modes are pure full-vectorial phenomena. From a theoretical point of view bandgaps do not arise as a consequence of coupling between field components at a dielectric interface, as it is the case in a full-vectorial method. Rather, they are a fundamental property of applying the wave equation to a periodic waveguide structure, and thus photonic bandgaps and defect modes can also exist in a scalar calculation.

The question is how well a scalar calculation actually approximates the results in photonic bandgap fibres, where the scalar wave equation is obviously not a correct representation of the actual physical problem. In order to examine this, we have chosen to apply our finite-difference method to the extreme case of airguiding photonic bandgap fibres. The dielectric structure used in our calculation is shown in figure 4a. It has a width of $`6\sqrt{3}\mathrm{\Lambda }`$ and consists of an air-silica cladding structure with a holediameter of $`D=0.88\mathrm{\Lambda }`$. The core defect is made by inserting a large air hole with a diameter of $`2.65\mathrm{\Lambda }`$. For the calculation, we discretized the structure in $`256\times 256`$ grid points. We have chosen this specific structure, because it is known to support guided modes in a full-vectorial calculation .

A disadvantage of this finite-difference implementation is that there is no simple way of calculating the position of the photonic bandgaps. Therefore, we have used a full-vectorial planewave-method to calculate the bandaps of an infinite triangular lattice with a hole diameter of $`0.88\mathrm{\Lambda }`$. Once the position of the photonic bandgaps have been located, it is possible to search for a defect mode. As already mentioned, the *EIGS*-command in *Matlab* has the useful ability to find eigenvectors around a specified eigenvalue. This is particularly useful for bandgap fibres, since the defect mode appears as an isolated eigenmode inside the boundaries of the photonic bandgap.

By choosing a normalized wavelength of $`\lambda /\mathrm{\Lambda }=0.7`$ and searching for an eigenvalue for which $`\beta k`$, the scalar finite-difference method finds a defect mode localized to the core region. The field distribution of this defect mode is shown in figure 4b. For simplicity we have chosen to depict the absolute value of the field distribution. Therefore, the 6 lobes surrounding the central defect do in fact have the opposite amplitude of the field inside the defect. These 6 resonances are a common feature of airguiding fibres, and they are also found in a full-vectorial calculation.

In figure 5 we have mapped out the effective mode-index for a range of the normalized wavelength. For comparison we have included the guided modes and the bandgap edges from a full-vectorial calculation. Both methods are found to predict the existence of a fundamental and a second-order mode, and a reasonable agreement is found between the results of the two methods. However, we generally find that the finite-difference method overestimates the values of the effective mode-index.

The major difference between the full-vectorial and the scalar calculation is that the latter is seen to predict significantly increased bandgaps. This naturally gives rise to a much wider range in which the structure supports a confined mode. The bandedges in the scalar calculation have been found as the modes that lie just above and just below the defect modes. As the bandedge-modes are infact cladding modes, and thus well distributed over the entire cross section of the dielectric structure, they are more affected by the truncation of the $`𝚯`$-matrix than the well confined defect modes. We have tried to increase the width of the calculation domain and also the number of grid points, but this is not found to have any influence on the overall result. It is therefore believed that the increased bandgap size is a consequence of the scalar approximation.

The final result indicates, that although a scalar approach provides great insight to the basic physics of photonic bandgap fibres, it cannot reveal the complete picture. This is not really surprising. However, we still believe that this method is of great interest, mainly due to it simple implementation. Also, the model can easily be expanded to include periodic boundary conditions and a full-vectorial implementation is feasible as well.

## 4 Conclusion

The field of photonic crystal fibres has by now existed for almost a decade and the results of the research efforts will probably soon move inside the classroom and also be found in text-book material on fibre optics and electro-magnetic theory of photonic crystals. This also calls for simple approaches to modelling of micro-structured optical fibres which are easy to implement and which without too much effort can produce results which are in qualitative agreement with those observed in real micro-structured fibres. We believe that the present work provides such a simple approach.

In order to limit the complexity of the mathematical formulation of the problem, we have considered a scalar wave-equation which is solved by means of the most simple numerical approach to differential equations; the finite-difference approximation. With appropriate boundary conditions this results in a matrix eigenvalue problem. The matrix, $`𝚯`$, of this problem is highly sparse and has very simple analytical matrix elements which only depend on the lattice spacing, the frequency, and the dielectric function. Thus, no implementation of complicated basis functions is required. For a given frequency $`\omega `$ the propagation constant $`\beta `$ can easily be found by finding the eigenvalues of $`𝚯`$. By the aid of standard numerical routines, we are able to solve for a specific number of eigenvalues closest to a specified value. This is particularly useful for bandgap guiding fibres, where the modes of interest corresponds to either the smallest or the largest eigenvalue, but are placed in an interval determined by the photonic bandgap edges.

In conclusion we find that the presented finite-difference method, apart from being simple to implement and despite the extremely simple model, is able to provide qualitative correct results for both index-guiding and photonic bandgap guiding fibres. The latter case is quite surprising, since modelling of photonic bandgap effects is normally associated with full-vectorial methods. Furthermore, we find that the method is robust and it is relatively simple to incorporate periodic boundary conditions, or to expand the model to a full-vectorial method. This holds interesting prospects for a future development of this method.

## Acknowledgements

We acknowledge useful discussions with M. D. Nielsen and B. T. Kuhlmey. J. L. is financially supported by the Danish Technical Research Council (STVF).

## References


text:
# Graphs with no 2𝛿+1 cycle

## 1 Introduction

Dirac  proved that any 2-connected graph with minimum vertex degree $`\delta `$ contains either a cycle of length at least $`2\delta `$ or a Hamilton cycle. In this paper, we shall characterize those graphs with minimum vertex degree $`\delta `$ which have no cycle of length greater than $`2\delta `$. This characterization was also motivated by Ali and Staton  in which similar results were given for graphs with no path exceeding $`2\delta +1`$. In particular, as a case in their work, they prove that if the number of vertices of a non-hamiltonian 2-connected graph is exactly $`2\delta +1`$, then the graph must be isomorphic to the join of a graph on $`\delta `$ vertices and a totally disconnected graph on $`\delta +1`$ vertices. Here, we will show that this same type of structure is present when the number of vertices exceeds $`2\delta +1`$.

## 2 Preliminaries

Unless specified otherwise, the terminology used here will follow Bondy and Murty . In particular, a graph may have loops and parallel edges. If $`v`$ is a vertex of a graph $`G`$, then the set of neighbors of $`v`$ will be denoted $`N(v)`$ or $`N_G(v)`$, and the minimum vertex degree of $`G`$ will be denoted $`\delta `$ or $`\delta (G)`$. If $`S`$ is a subset of the vertex set of $`G`$, then $`G[S]`$ denotes the induced subgraph of $`G`$ on $`S`$. The totally disconnected graph on $`m`$ vertices will be denoted $`N_m`$, and any subgraph of $`K_n`$ will be denoted $`H_n`$. The join of two graphs $`G`$ and $`H`$ is the graph $`GH`$ and is obtained by taking disjoint copies of $`G`$ and $`H`$ and adding edges joining every vertex of $`G`$ to every vertex of $`H`$. In particular, if $`|V(G)|=n`$ and $`|V(H)|=m`$, then the graph $`GH(E(G)E(H))`$ is isomorphic to $`K_{n,m}`$. Finally, if $`x_1x_2\mathrm{}x_n`$ is a path $`P`$ in a graph $`G`$, then $`P[x_i,x_j]`$, $`P[x_i,x_j)`$, $`P(x_i,x_j]`$, and $`P(x_i,x_j)`$, will denote the subpaths $`x_ix_{i+1}\mathrm{}x_j`$, $`x_ix_{i+1}\mathrm{}x_{j1}`$, $`x_{i+1}x_{i+2}\mathrm{}x_j`$, and $`x_{i+1}x_{i+2}\mathrm{}x_{j1}`$, respectively.

## 3 The Main Theorem

To establish the main result of this paper, we shall use the following theorem of Dirac .

###### Theorem 3.0.1

Let $`G`$ be a simple $`2`$-connected graph with minimum vertex degree $`\delta `$ and suppose $`|V(G)|3`$. Then $`G`$ contains either a cycle of length at least $`2\delta `$ or a Hamilton cycle.

Motivated by this result, we will now prove the following characterization for the graphs having no cycle longer than $`2\delta `$.

###### Theorem 3.0.2

Let $`G`$ be a simple $`2`$-connected graph with minimum vertex degree $`\delta `$. If $`G`$ has no cycle of length at least $`2\delta +1`$, then $`G`$ is Hamiltonian or $`GH_\delta N_m`$ where $`m>\delta `$.

Proof. Let $`G`$ be a graph with no Hamilton cycle and no cycle of length at least $`2\delta +1`$. By Theorem 3.0.1, G has a cycle $`C`$ of length $`2\delta `$. Now, since $`G`$ is not Hamiltonian, there is a vertex $`z`$ in $`V(G)V(C)`$, and by Menger’s Theorem , there are 2 paths $`Q_1`$ and $`Q_2`$ from $`z`$ to $`C`$ that have only the vertex $`x`$ in common and such that each $`Q_i`$ meets $`C`$ in exactly one vertex $`z_i`$.

###### Lemma 3.0.3

$`Q_1`$ and $`Q_2`$ must not meet consecutive vertices on $`C`$.

Proof. Suppose $`Q_1`$ and $`Q_2`$ meet $`C`$ at $`z_1`$ and $`z_2`$ where $`z_1z_2`$ is an edge of $`C`$. Then $`G`$ has a cycle $`C^{}`$ formed by the subpath $`C\{z_1z_2\}`$ together with the path $`z_1zz_2`$. Therefore, $`C^{}`$ contains $`2\delta +1`$ vertices of $`G`$, namely, all the vertices in $`V(C)\{z\}`$; a contradiction. $`\mathrm{}`$

Since there are paths from every vertex of $`V(G)V(C)`$ to $`C`$, choose a longest path $`P`$ subject to the condition that one of the endvertices of $`P`$ lies on $`C`$ and no other vertex of $`P`$ lies on $`C`$. We label the endvertex of $`P`$ that is not in $`V(C)`$ as $`x`$ and the vertex in $`V(P)V(C)`$ is labelled $`y`$. Now, since $`P`$ is a longest path of this type, it is clear that every neighbor of $`x`$ lies in $`V(P)V(C)`$. Moreover, we have the following lemma.

###### Lemma 3.0.4

$`x`$ has a neighbor in $`V(C)y`$.

Proof. Suppose not. Then every neighbor of $`x`$ lies on $`P`$. Since $`G`$ is $`2`$-connected, it is clear that there is a vertex of $`P`$ adjacent to some vertex on $`C`$. Viewing $`P`$ as a directed path from $`x`$ to $`y`$, let $`z`$ be the first vertex on $`P`$ that is adjacent to a vertex, $`z^{}`$, of $`C`$, and let $`u_x`$ be the first neighbor of $`x`$ after $`z`$ on $`P`$. Now, $`C`$ is partitioned into two paths from $`z^{}`$ to $`y`$, and we arbitrarily label them $`C_1`$ and $`C_2`$. Since $`C`$ has size exactly $`2\delta `$, one of $`C_1`$ and $`C_2`$ has at least $`\delta `$ edges. Without loss, assume that $`C_1`$ has at least $`\delta `$ edges.

Now, let $`P^{}`$ be the path containing $`z^{}z`$, $`P[z,x]`$, $`xu_x`$, and $`P[u_x,y]`$. Notice that $`P^{}`$ contains all the vertices of $`P`$ except those on $`P(z,u_x)`$. Since there are no neighbors of $`x`$ on $`P(z,u_x)`$, it is clear that in addition to $`x`$ itself, $`P^{}`$ contains all the neighbors of $`x`$; thus, $`P^{}`$ contains at least $`\delta +1`$ vertices from $`P`$ as well as the vertex $`z^{}`$ on $`C`$. Since $`C_1`$ contains $`\delta +1`$ vertices, we can combine $`C_1`$ and $`P^{}`$ to obtain a cycle $`C^{}`$ of size $`|V(C_1)|+|V(P^{})|2(\delta +1)+(\delta +2)2=2\delta +1`$; a contradiction. $`\mathrm{}`$

Having proved Lemma 3.0.4, let the vertices of $`Cy`$ that are adjacent to $`x`$ be labelled $`x_1,x_2,\mathrm{},x_k`$ in cyclic order on $`C`$, where $`y`$ lies between $`x_k`$ and $`x_1`$ on $`C`$, and relabel $`y`$ as $`x_{k+1}`$. It is clear that the cycle $`C`$ is partitioned into $`k+1`$ internally disjoint paths $`C(x_1),C(x_2),\mathrm{},C(x_{k+1})`$ where $`C(x_i)`$ is the subpath of $`C`$ from $`x_i`$ to $`x_{i+1}`$ that contains no other vertex in $`\{x_1,x_2,\mathrm{},x_{k+1}\}`$; here $`C(x_{k+1})`$ is the subpath of $`C`$ from $`x_{k+1}`$ to $`x_1`$. If $`k=1`$, then $`C`$ can be partitioned into two paths from $`x_1`$ to $`x_2`$. We arbitrarily label one of these $`C(x_1)`$ and the other $`C(x_2)`$. Since $`G`$ has no cycle containing $`2\delta +1`$ vertices, it is clear by Lemma 3.0.3 that Each $`C(x_i)`$ contains at least one vertex in its interior, and we let the vertex of $`C(x_i)`$ that is adjacent to $`x_i`$ be labelled $`a_i`$.

###### Lemma 3.0.5

There is no path from $`a_i`$ to $`a_j`$ avoiding $`xC`$.

Proof. Suppose $`Q`$ is a path joining $`a_i`$ to $`a_j`$ and avoiding $`xC`$. It is easy to see that $`G`$ has a cycle $`C^{}`$ formed by the edges of $`C\{a_jx_j,x_ia_i\}`$ together with the path $`x_jxx_i`$ and the edges of $`Q`$. Thus, $`C^{}`$ contains the vertices of $`C`$ and $`x`$; a contradiction. $`\mathrm{}`$

Now, consider $`C(x_k)`$ and $`C(x_{k+1})`$. These paths have the vertex $`x_{k+1}`$ as an endvertex. If one of these paths has less than $`|P|`$ internal vertices, then we can delete the path from $`C`$ and replace it with the path $`P`$ and the edge $`xx_k`$ or $`xx_1`$ creating a longer cycle. Thus, each of $`C(x_k)`$ and $`C(x_{k+1})`$ has length at least $`|P|+1`$. So, the length of $`C(x_j)`$ is greater than or equal to $`|P|+1`$ for $`j\{k,k+1\}`$, and by Lemma 3.0.3, the length of $`C(x_i)`$ is at least $`2`$ for $`i\{1,\mathrm{},k1\}`$. This implies that $`|C|2(|P|+1)+(k1)2`$ which means that $`2\delta 2|P|+2k`$. Therefore, $`\delta |P|+k`$.

Now, since $`d_G(x)|P|+k`$ and $`x`$ is adjacent to exactly $`k`$ vertices in $`V(C)x_{k+1}`$, the vertex $`x`$ must be adjacent to every vertex on $`Px`$. This implies that $`d_G(x)=|P|+k`$ which in turn means that $`\delta =|P|+k`$. Furthermore, this restriction shows that each $`C(x_i)`$ has length 2 when $`i\{1,\mathrm{}k1\}`$ and that both $`C(x_k)`$ and $`C(x_{k+1})`$ have length $`|P|+1`$.

###### Lemma 3.0.6

$`|P|=1`$. Moreover, $`\delta =k+1`$ and $`N_G(x)=\{x_1,\mathrm{},x_{k+1}\}`$.

Proof. Suppose $`|P|>1`$ and let the vertices of $`P`$ from $`x`$ to $`x_{k+1}`$ be labelled $`u_1,u_2,\mathrm{},u_{|P|}`$ in this order where $`u_{|P|}=x_{k+1}`$. Since $`xu_{i+1}`$ is an edge, we have a path $`P_i=P\{u_iu_{i+1}\}xu_{i+1}`$ for each $`u_i`$, and this path is a longest path having one endvertex on $`C`$ and the other end being $`u_i`$. Since $`d_G(u_i)\delta `$, the vertex $`u_i`$ must be adjacent to at least $`k`$ vertices on $`Cx_{k+1}`$, and using Lemmas 3.0.3 and 3.0.5, it is easy to see that the neighbors of $`u_i`$ on $`Cx_{k+1}`$ must be exactly the vertices in $`\{x_1,x_2,\mathrm{},x_k\}`$. So, $`u_i`$ must be adjacent to every member of $`P_i`$. But, $`(C\{a_1\})\{x_2xu_1x_1\}`$ is a $`2\delta +1`$ cycle, since $`u_1x_{k+1}`$. $`\mathrm{}`$

Now, by Lemma 3.0.6, every path from a vertex in $`V(G)V(C)`$ to $`V(C)`$ must have length $`1`$; so, since $`G`$ is connected, every vertex in $`V(G)V(C)`$ has a neighbor in $`V(C)`$. If $`x^{}`$ and $`x^{\prime \prime }`$ are distinct members of $`V(G)V(C)`$, then $`x^{}x^{\prime \prime }`$ is not an edge of $`G`$; otherwise if $`y`$ is a neighbor of $`x^{}`$ on $`C`$ then $`yx^{}x^{\prime \prime }`$ is a path of length 2 having only one common vertex with $`C`$. We conclude that $`G[V(G)V(C)]`$ is totally disconnected.

Upon combining Lemmas 3.0.3, 3.0.5, and 3.0.6, we see that $`N(x)=\{x_1,\mathrm{}x_{k+1}\}`$ for every $`x`$ in $`V(G)V(C)`$. Moreover, the set of neighbors for each $`a_i`$ is $`\{x_1,\mathrm{}x_{k+1}\}`$ since $`a_ia_j`$ is not an edge of $`G`$ and $`\delta =k+1`$. Therefore, $`G[\{a_1,\mathrm{},a_{k+1}\}(V(G)V(C))]`$ is totally disconnected with $`k+1+|V(G)V(C)|=m`$ vertices. Since $`V(G)V(C)`$ has at least one vertex, $`m`$ exceeds $`\delta `$, and the theorem is established. $`\mathrm{}`$

###### Corollary 3.0.7

Let $`G`$ be a $`2`$-connected graph with at least $`2\delta +1`$ vertices. Then, $`G`$ has no cycle of length at least $`2\delta +1`$ if and only if $`GH_\delta N_m`$ where $`m>\delta `$.


text:
# 1 Introduction

## 1 Introduction

Within general relativity it is well-known that there are three stages to the dynamical evolution of field perturbation on a black hole background: the initial outburst from the source of perturbation, the damping (quasi-normal) oscillations and asymptotic tails at very late time. The evolution significantly depends on the asymptotic behavior of the background: for different values of cosmological constant $`\mathrm{\Lambda }`$ there are three types of behavior depending on the sign of $`\mathrm{\Lambda }`$. Thus, the remarkable fact is that while for asymptotically flat backgrounds we have inverse power law tails, for asymptotically Anti-de Sitter background it is exponential low that governs the decay of a field at all times. The same exponential decay at very late times was observed by Brady et al  for asymptotically de Sitter background.

Apart from original interest coming from possibility to observe quasi-normal ringing with the help of gravitational wave detectors (see, e.g.,  for review), the quasi-normal modes bring now a lot of interest in different contexts: in Anti-de Sitter(de Sitter)/Conformal Field Theory correspondence , when considering thermodynamic properties of black holes in loop quantum gravity , in the context of possible connection with critical collapse .

The gravitational QNMs of Schwarzschild-de Sitter black hole were considered in . Recently, asymptotic values of QNMs were studied for nearly extremal Schwarzschild-de Sitter black hole and it was shown that in the near extremal limit the effective potential of the wave equation reduces to the Pöshl-Teller potential, and it is expected therefore the Pöshl-Teller method provides the best accuracy near the extremal value of $`\mathrm{\Lambda }`$ .

QNMs corresponding to decay of fields of different spin are extensively studied for Schwarzschild and Schwarzschild-anti de Sitter backgrounds. At the same time there is lack of such study of the Schwarzschild-de Sitter case. This motivated us to calculate QNMs of scalar, electromagnetic, gravitational and Dirac fields for Schwarzschild-de Sitter black hole. We also compare the results found here with those obtained recently in .

In section 2 we give the perturbation equations to describe fields of different spin in SdS background. Section 3 presents comparison of QNM frequencies obtained through sixth order WKB formula and approximation by the Pöshl-Teller potential.

There was shown that the frequencies of the massless Dirac field perturbations  are the same for both chiralities in SdS background as well. Near the extremal value of $`\mathrm{\Lambda }`$ the sixth order WKB formula gives the values which are very close to those obtained within the Pöshl-Teller potential approximation. This agrees with recent work of Cardoso and Lemos .

## 2 Perturbation equations and methods

It is well known that the Schwarzschild-de Sitter black hole of mass $`M`$ is described by the metric:

$$ds^2=f(r)dt^2\frac{dr^2}{f(r)}r^2d\sigma ^2;f(r)=1\frac{2M}{r}\mathrm{\Lambda }\frac{r^2}{3},d\sigma ^2=d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2.$$
(1)
A perturbation equation on this background can be reduced to the Schrödinger wave-like equations:

$$\left(\frac{d^2}{dr^2}+\omega ^2V(r^{})\right)\mathrm{\Psi }(r^{})=0,$$
(2)
where we assume exponential dependence of $`\mathrm{\Psi }`$ on time ($`e^{i\omega t}`$), $`\omega =\omega _{Re}i\omega _{Im}`$, $`r^{}`$ is the “tortoise” coordinate:

$$dr^{}=\frac{dr}{f(r)}.$$
(3)
Under the positive real part QNMs, by definition, satisfy the following boundary conditions:

$$\mathrm{\Psi }(r^{})C_\pm \mathrm{exp}(\pm i\omega r^{}),r\pm \mathrm{},$$
(4)
which correspond to purely out-going waves at infinity and in-going waves at the horizon.

The corresponding perturbation equations can be rewritten in form (2) with potentials

$`V_s`$ $`=`$ $`f(r)\left({\displaystyle \frac{l(l+1)}{r^2}}+{\displaystyle \frac{1}{r}}f^{}(r)\right);`$ (5)
$`V_{el}`$ $`=`$ $`f(r){\displaystyle \frac{l(l+1)}{r^2}};`$ (6)
$`V_{odd}`$ $`=`$ $`f(r)\left({\displaystyle \frac{l(l+1)}{r^2}}{\displaystyle \frac{6M}{r^3}}\right);`$ (7)
$`V_{even}`$ $`=`$ $`{\displaystyle \frac{2f(r)}{r^3}}{\displaystyle \frac{9M^3+3c^2Mr^2+c^2(1+c)r^3+3M^2(3cr\mathrm{\Lambda }r^3)}{(3M+cr)^2}},c={\displaystyle \frac{l(l+1)}{2}}1;`$ (8)
$`V_d`$ $`=`$ $`{\displaystyle \frac{\kappa f^{1/2}(r)}{r^2}}(\kappa f^{1/2}(r)\pm ({\displaystyle \frac{3M}{r}}1));`$ (9)
for massless scalar, electromagnetic, gravitational (odd end even parities) and Dirac fields respectively (see, e.g.,  and references therein), where $`l`$ is the multipole number, $`\kappa =l+1`$ for (+) sign and $`\kappa =l`$ for (-) sign in (9). All the potentials vanish at the event horizon and the cosmological horizon, having a peak near the event horizon. It allows to use the Pöshl-Teller potential approximation as well as WKB method.

The WKB approach is based on the analogy with the problem of waves scattering near the peak of the potential barrier $`V(r^{})`$ in quantum mechanics, where $`\omega ^2`$ plays a role of energy. The approach was used by Schutz and Will , developed to the third order by Iyer and Will , and recently developed by Konoplya to the sixth order beyond the eikonal approximation . The result has the form:

$$i\frac{\omega ^2V_0}{\sqrt{2V_0^{\prime \prime }}}\mathrm{\Lambda }_2\mathrm{\Lambda }_3\mathrm{\Lambda }_4\mathrm{\Lambda }_5\mathrm{\Lambda }_6=n+\frac{1}{2},$$
(10)
where $`V_0`$ is the height and $`V_0^{\prime \prime }`$ is the second derivative with respect to the tortoise coordinate of the potential at the maximum. $`\mathrm{\Lambda }_2`$ and $`\mathrm{\Lambda }_3`$ can be found in , $`\mathrm{\Lambda }_4`$, $`\mathrm{\Lambda }_5`$ and $`\mathrm{\Lambda }_6`$ are presented in ; the corrections depend on the value of the potential and higher derivatives of it at the maximum.

Generally, WKB method was effectively used for finding quasi-normal frequencies in a lot of papers . Comparing the values calculated by the different order WKB formula one can judge about the convergence to the some unknown accurate result as the WKB order increasing. We expect that the possible error is less then the difference between the sixth and the fifth order WKB values111The accuracy of values are to take as a half of the number last order (e. g. $`0.125`$ means $`0.125\pm 0.003`$), but it is much better for most of them.. Note, that standard machine precision in MATHEMATICA gives large errors due to numerical truncations. To exclude such “noises” we had to hold 40 digits in all intermediate values.

For comparison we will use the semi-analytical method for calculations proposed by B. Mashhoon  who used Pöshl-Teller approximate potential

$$V_{PT}=\frac{V_0}{\mathrm{cosh}^2(r^{}/b)}.$$
(11)
It contains two free parameters ($`V_0`$ and $`b`$) which are used to fit the height and the second derivative of the potential $`V(r^{})`$ at the maximum .

The quasi-normal modes of the Pöshl-Teller potential can be evaluated analytically:

$$\omega =\frac{1}{b}\left(\sqrt{V_0b^2\frac{1}{4}}(n+\frac{1}{2})i\right).$$
(12)
## 3 QNMs of SdS black holes

Both methods give best result for low overtone number but it is known that the higher modes give larger values of the imaginary part of frequencies . Therefore these modes decay faster then low-lying ones that are more interesting for us. So we calculated some first quasi-normal mode frequencies for massless scalar, electromagnetic and gravitational and Dirac fields using both of the methods. The results are presented below (all values are measured in the black hole mass (M) units).

QNMs of massless scalar field: $`l=1`$, $`n=0`$

QNMs of massless scalar field: $`l=2`$, $`n=0`$

QNMs of massless scalar field: $`l=2`$, $`n=1`$

QNMs of electromagnetic field: $`l=1`$, $`n=0`$

QNMs of electromagnetic field: $`l=2`$, $`n=0`$

QNMs of electromagnetic field: $`l=2`$, $`n=1`$

QNMs of gravitational field: $`l=2`$, $`n=0`$

QNMs of gravitational field: $`l=2`$, $`n=1`$

QNMs of gravitational field: $`l=3`$, $`n=0`$

QNMs of gravitational field: $`l=3`$, $`n=1`$

QNMs of massless Dirac field: $`\kappa =1`$, $`n=0`$

QNMs of massless Dirac field: $`\kappa =2`$, $`n=0`$

The relative error of the WKB method depends significantly on a field under consideration. Comparison of the sixth and fifth WKB order let us to conclude that the worst WKB convergence gives the Dirac field, where the expected error reaches 1% in the sixth WKB order. The results calculated in  by using the third WKB order differ from the sixth order values up to 4%. Thus, the possible relative error of the third order WKB results may exceed 4%.

We can also evaluate the relative error order for gravitational modes by comparison QNMs of odd and even perturbations. Because the potentials can be expressed in the form

$$V_{odd}=W_g^2+\frac{dW_g}{dr^{}}+\beta ;V_{even}=W_g^2\frac{dW_g}{dr^{}}+\beta ;$$
(13)
where

$$W_g=\frac{2M}{r^2}\frac{3+2c}{3r}+\frac{3c^2+2c^39\mathrm{\Lambda }M^2}{3c(3M+cr)}\frac{1}{3M}(c^2+c\frac{3\mathrm{\Lambda }M^2}{c});\beta =\frac{c^2(c+1)^2}{9M^2};$$

any solution for even parity can be found from the solution for odd parity :

$$\mathrm{\Psi }_{even}(r^{})=q\left(W_g(r^{})\frac{d}{dr^{}}\right)\mathrm{\Psi }_{odd}(r^{});q=const.$$

Obviously, if a solution $`\mathrm{\Psi }_{odd}`$ describing quasi-normal oscillations satisfies the boundary conditions (4), the $`\mathrm{\Psi }_{even}`$ solution satisfies the same ones describing quasi-normal oscillations too.

The same fact  takes place for massless Dirac QNMs of opposite chiralities:

$$V_d=W_d^2\pm \frac{dW_d}{dr^{}};W_d(r)=\frac{\kappa \sqrt{f(r)}}{r}.$$
(14)
Therefore, the QNM spectrum is the same for odd and even gravitational perturbations and for opposite chirality of Dirac perturbations in SdS background as well as in asymptotically flat universe. We see that the difference between the numerical results does not exceed our expectations.

We can see that the imaginary parts of QNMs, calculated using both sixth-order WKB and Pöshl-Teller approaches, depend on multipole and overtone number and black hole mass only near the extremal value of $`\mathrm{\Lambda }`$. This agrees with , where the QNMs of near extremal SdS space-time were found:

$$\omega b=\left(n+\frac{1}{2}\right)i+\sqrt{l(l+1)\frac{1}{4}}$$
(15)
for scalar and electromagnetic perturbations,

$$\omega b=\left(n+\frac{1}{2}\right)i+\sqrt{(l+2)(l1)\frac{1}{4}}$$
(16)
for gravitational perturbations. Here

$$b=\frac{54M^3}{(r_c3M)(r_c+6M)},r_c3M.$$

Using the same technique, one can find that

$$\omega b=\left(n+\frac{1}{2}\right)i+\sqrt{\kappa ^2\frac{1}{4}}$$
(17)
for massless Dirac field perturbations. It means in the extremal limit fields of different spin decay with the same rate. Note, that formulae (15), (16), (17) give right results for any overtone number $`n`$ and may be useful for studying the $`n\mathrm{}`$ limit.

The well-known approximate formula for large multipole number $`l`$ (or $`\kappa `$ in (9)) :

$`\omega `$ $`=`$ $`{\displaystyle \frac{1}{3\sqrt{3}M}}\left(l+{\displaystyle \frac{1}{2}}\left(n+{\displaystyle \frac{1}{2}}\right)i\right)+O\left({\displaystyle \frac{1}{l}}\right);`$ (18)
$`\omega `$ $`=`$ $`{\displaystyle \frac{1}{3\sqrt{3}M}}\left(\kappa \left(n+{\displaystyle \frac{1}{2}}\right)i\right)+O\left({\displaystyle \frac{1}{\kappa }}\right);`$
can be generalized for the case of non-zero $`\mathrm{\Lambda }`$. Making use of the first-order WKB method or the formula of the approximation by the Pöshl-Teller potential one can find for the non-extremal SdS black hole:

$`\omega `$ $`=`$ $`{\displaystyle \frac{\sqrt{19M^2\mathrm{\Lambda }}}{3\sqrt{3}M}}\left(l+{\displaystyle \frac{1}{2}}\left(n+{\displaystyle \frac{1}{2}}\right)i\right)+O\left({\displaystyle \frac{1}{l}}\right);`$ (19)
$`\omega `$ $`=`$ $`{\displaystyle \frac{\sqrt{19M^2\mathrm{\Lambda }}}{3\sqrt{3}M}}\left(\kappa \left(n+{\displaystyle \frac{1}{2}}\right)i\right)+O\left({\displaystyle \frac{1}{\kappa }}\right).`$
## 4 Conclusion

The QNMs of Schwarzschild-de Sitter black holes for fields of different spin have been calculated. The modes are determined by the black hole mass $`M`$ and the cosmological constant $`\mathrm{\Lambda }`$ only. The frequencies all have a negative imaginary part, which means that the black hole is stable against these perturbations. The presence of the cosmological constant leads to decrease of the real oscillation frequency and to a slower decay. The interesting problem that was outside our consideration is the search of asymptotic behavior of SdS QNMs at large imaginary part ($`n\mathrm{}`$) for which other approaches should be explored .

## Acknowledgements

I would like to thank R. Konoplya for proposing this problem to me and reading the manuscript.